高速度侵徹の侵徹効率

前回の記事では、高速度侵徹におけるエネルギー効率を考えました。この例では装甲の抵抗が衝突速度と同じオーダーで増加するので、その極限はエネルギー効率は密度から定まる非0な値となりました。

よく知られているように、高速度侵徹の高速度側では侵徹深さは侵徹体の速度に依存しなくなります。ということで、運動エネルギーに対する侵徹深さの効率は高速度側で0になります。ここで、侵徹効率$\eta_p$をどのように考えるかは議論の余地がありますが、

\[ \eta_p(V_0) = P(V_0)\times R_t(V_0)/KE(V_0) \]

として定義してみます。ここで、$V_0$は衝突速度であり、$R_t$は標的強度項です。つまり、投入された運動エネルギーのうち、どれだけ標的強度項が消耗したエネルギーを占めるかという定義になります。

以下に、Anderson-Walkerモデルを用いて計算した$\eta$、$\eta_p$と$V_0$の関係を示します。ここで、標的は鉄($\rho$:7800g/cm3, Y:1GPa)、侵徹体はWHA($\rho$:17800, Y:1GPa)です。

結論としては見てのとおりですね。

高速度侵徹のエネルギー効率

侵徹体の運動エネルギーが標的の侵徹にすべて利用される低速度侵徹と異なり、高速度侵徹では侵徹体が消耗するために初期の運動エネルギー$KE0$は標的が受け止めるエネルギー$W$とは一致しません。

そこで、高速度側の極限では初期の運動エネルギーの何割くらいが侵徹に利用されているのかふと気になりました。

高速度侵徹の極限では侵徹深さ$P$は、侵徹体長さ$L$、侵徹体密度$\rho_p$、標的密度$\rho_t$を用いて

$$P = \sqrt{\frac{\rho_p}{\rho_t}}L $$

と表されます。一方、衝突時の侵徹体の運動エネルギー$KE$は、衝突速度$V_0$を用いて

$$KE = \frac 1 2 L\rho_p V_0^2$$

と表されます。標的の侵徹に使用されるエネルギーは、少し考えどころです。エネルギーは力×距離ですから、侵徹深さがわかってる現在、標的が作る力がわかれば直ちに求められます。標的が作る力はベルヌーイの式を参考にすれば、標的の侵徹に使用されるエネルギー$W$は

$$ W= \frac{1}{2}\rho_t u^2\times P = \frac{1}{2}\rho_t u^2\times \sqrt{\frac{\rho_p}{\rho_t}}L $$

として簡単に得られます。$W$と$KE$の比がエネルギー効率ですから、比を取って適当に解くと、エネルギー効率$\eta$は

$$ \eta = \frac{W}{KE} = \frac{\rho_{t} \sqrt{\frac{\rho_{p}}{\rho_{t}}}}{\rho_{p} \left(\sqrt{\frac{\rho_{t}}{\rho_{p}}} + 1\right)^{2}} $$

として簡単に得られます。まあ適当に、侵徹体としてタングステンの密度(19.3g/cm3)、

標的の密度として鉄の密度(7.8g/cm3)を代入してみれば、$\eta=0.24$程度となり、

高速度侵徹のエネルギー効率は24%程度になることがわかります。

以上、ふと思って計算してみましたが、たぶん昔の偉い人が似たような計算をしている気がします。

AWモデルにおけるDynamic cavitation analysisの解析解

完全にメモ書きです。

ちょっと頑張ってL/D効果も含めた侵徹深さの評価をしたい話で

\[ \left(1+\frac{\rho_tu^2}{Y_t}\right)\sqrt{K_t-\rho_tu^2\alpha^2}= \left(1+\frac{\rho_tu^2\alpha^2}{2G_t}\right)\sqrt{K_t-\rho_tu^2} \]

に対してニュートンラフソン法、あるいは二分法を用いれば$\alpha$を求められると書いていました。これは、この式が平方根を含むために直接解析解を求めることができないためで、AWモデルを用いた計算の大部分を$\alpha$を求めるために消費していました。

一方で、$\alpha$の値域は$\alpha \gt 1$の実数というのは明らかですので、

\[ \left(1 + \frac{\rho u^{2}}{Y}\right)^{2} \left(K_{t} - \alpha^{2} \rho u^{2}\right) = \left(1 + \frac{\alpha^{2} \rho u^{2}}{2 G}\right)^{2} \left(K_{t} - \rho u^{2}\right) \]

とした解を求めても基本的には等価です。ここで$\alpha^2$については明らかに解析解が存在していて、sympyで適当に求めると、

…ここに示しておきたかったんですが長すぎて数式に表示されませんでした。

import sympy 
rho,u2, Y, K_t, alpha, G = sympy.symbols("rho, u^2, Y, K_t, alpha^2, G")
lh = (1 + rho*u2 /Y) * (1 + rho*u2 / Y)* (K_t - rho*u2 * alpha)
rh =  (1 + rho*u2 * alpha * 0.5 / G) * (1 + rho*u2 * alpha /2 / G) * (K_t - rho*u2)
eq = sympy.Eq(lh,rh)
s1,s2 = sympy.solve(eq,alpha)
s2

と書くことが出来ます、とだけ示しておきます。これの平方根を取れば、$\alpha$になります。数値解との対応は以下のとおりです。

個人的には計算時間が4msから1msまで減らせたので満足です。

ノイマン効果ってなんだ？な話

ノイマン効果ってなんだ

いつの間にか、ブログ記事の多くが終末弾道の話になっていて、自分の出自と無関係な方向にどんどんと進んでいるな、という感じがあります。
まぁ、それはともかくとして、終末弾道を考えるとき簡単の場合には侵徹体はバルクと考えてよく、侵徹体は砲から発射されたままの形状を初期条件として侵徹を考えることができます。このような例として、低速度のAPからAPFSDSまでが該当します。
一方で、侵徹体をバルクとして考えられない例として、侵徹体をその場で形成するHEATが考えられます。

HEATの侵徹体の形成について説明を行う際には、モンロー効果、ノイマン効果といった効果が引き合いに出されるかと思います。
こういう効果の説明を簡潔に行っているのは、やはりWikipediaですので、ノイマン効果に関する記述を引っ張ってきますと、

ノイマン効果（英: Neumann effect）とは、ドイツの科学者、エゴン・ノイマン（Egon Neumann）が1910年に発見した、モンローの円錐形のくぼみに金属板で内張り（くぼみと同じ形の金属の円錐をはめ込むこと）をすると穿孔力がさらに強くなる現象。

爆轟が進行して金属の内張りに達すると、爆轟波によりライナーは動的超高圧に晒されユゴニオ弾性限界を超える圧力に達すると固体の金属でも可塑流動性を持つようになり、液体に近似した挙動を示すようになる。これにより、融着体と呼ばれる金属塊となって前方へ超音速で飛び出していく。

モンロー/ノイマン効果-Wikipedia

のように説明されるようです。結局の所、ノイマン効果というのは「円錐形のくぼみに金属板で内張り（くぼみと同じ形の金属の円錐をはめ込むこと）をすると穿孔力がさらに強くなる現象。」として説明される現象と考えていいようです。
さて、では実際にどのような現象が起こっているかという話ですがWikipediaの記事はなかなか言っていることが難しくてどういうことが起きているのかなかなかわかりません。
そこで、今回の記事では1948年にBirkhoff, MacDougall, Pugh, Taylorにより提案された最も単純なモデルに則ってノイマン効果がどのようにして生じているかを概観したいと思います。
今回の記事はBallisticsのShaped chargeの項を元にしたものになります。

ライナーの中で成り立つ幾何学的関係

今回の記事では以下の図のような状況を考えます。

 Fig. 1 BMPTモデルで考える前提

爆轟速度$D$[m/s]の炸薬中に角度$2\alpha$[rad]の開き角で挿入された円錐状のライナーが挿入されており、爆轟は炸薬中を平面波的に右へ進んでいます。
爆轟が通過し終わったライナーは円錐中央部に叩き出されます。爆轟はライナー形状に関係なくまっすぐ右に進んでいますので、ライナーもまた円錐中央部から爆轟面まで一定の角度をなして分布していると考えられます。そこで、爆轟によって叩き出されたライナーの角度を$2\beta$[rad]とします。

Fig. 2 BMPTモデルで成り立つ各部分の速度と角度の関係

BMPTモデルでは、Fig. 1に示した初期ライナー角と、叩き出されたライナー角の関係から、Fig. 2のような速度の関係を幾何学的に示しました。

ここで、点Aは叩き出されたライナーが円錐中央部で集まっている点(ジャンクション)であり、$V_1$[m/s]はジャンクションが前方に進む速度です。また、点Pは今まさに爆轟が到達した点であり、点Bに位置するライナーは速度$V_0$[m/s]で点Bに向かって叩き出されます。この時、点Pに立った読者は点Aに向かって$V_2$[m/s]で接近します。
また、$D$[m/s]は炸薬の爆轟速度です。

$V_0$は炸薬によって叩き出されたライナーの速度であり、この図形中で唯一不明なパラメータです。$V_0$はGurney methodという炸薬のエネルギーと、ライナー質量/成形炸薬質量比の関係と成形炸薬の特性で決まるとされていますが、今回の記事ではこれは魔法の力によって所与のものとしましょう。

さて、そうすると未だ明らかでないパラメータは$\theta$[rad]と$V_1$、$V_2$、$\beta$[rad]となりますが、$\theta$についてはFig. 2から明らかに

$$\theta = \frac{\pi}{2} -\frac{\beta -\alpha}{2} $$

となります。そうすると、$V_1$は正弦定理から

$$\frac{V_0}{\sin \beta} = \frac{V_1}{\sin\left(\frac{\pi}{2} -\frac{\beta -\alpha}{2}\right) } = \frac{V_1}{\cos \left(\frac{\beta -\alpha}{2}\right)} $$

が直ちに得られます。
さらに$V_2$はオレンジ色の補助線に従うと、

$$V_2 = V_1\cos \beta + V_0 \cos \left(\frac{\pi}{2} -\frac{\beta -\alpha}{2}\right) =V_1\cos \beta + V_0 \sin \left(\frac{\beta -\alpha}{2}\right)  $$

が同様に得られます。
最後に、$\beta$についてですが、これは右側の三角形について考えると、

$$\frac{V_0}{{\sin \left(\pi-2\theta\right)}}=\frac{\frac{D}{\cos \alpha}}{\sin \theta}$$
$$\frac{V_0\cos\left(\left(\beta-\alpha\right)/2\right)}{\sin\left(\beta-\alpha\right)} = \frac{D}{\cos \alpha}$$

という$\beta$のみが未知の式が得られるので、これを数値的に解くことで得られます。
以上でFig. 2の項目はすべて明らかになりました。

ジャンクションでの速度関係

次に、ジャンクションで成り立つ関係に移りましょう。以下にジャンクションで成り立つ速度の関係について示します。

Fig. 3 ジャンクションで成り立つ速度の関係

点Pから速度$V_2$でライナーが点A(ジャンクション)に流れこむ時、ある部分は真後ろ(AC方向)に、ある部分は真正面(AB方向)に打ち出されるでしょう。この時、理想的な円錐であればその他の速度成分は持ちません。これは円錐系のライナーが軸対称に流れ込むために、軸方向以外の速度成分がすべてキャンセルされるためです。この、ジャンクションから放出されるライナーは軸方向以外の速度成分を持たず、真正面に打ち出されるという現象こそがノイマン効果をもたらす現象と言えるでしょう。

さて、一般に用いられるライナー材料は柔らかく、材料の強度項は無視していいでしょう。一方で、ジャンクションに流れ込む速度$V_2$はそれほど高くないので*、ここでは非圧縮性の流体として扱ってみましょう。

そうすると、ジャンクション(点A)に入る圧力と出る圧力はベルヌーイの式から釣り合っており、AC方向にでる速度$V_{AC}$と、AB方向に出る速度$V_{AB}$と$V_2$の間に成り立つ関係は、

$$\frac 1 2 \rho V_2^2 = \frac 1 2 \rho V_{AC}^2 = \frac 1 2 \rho V_{AB}^2$$

となり、Fig. 3の関係が得られます。

メタルジェット速度、スラグ速度

さて、Fig. 3の関係から、ジャンクションから出ていく物質の速度はAC方向、AB方向、いずれの方向にも$V_2$であることがわかりました。

一方で、ジャンクションは$V_1$で動いているので、AB方向に出る物質と、AC方向に出る物質を外から見ると、

$$V_{AB} = V_1 + V_2$$

$$V_{AC}= V_1 - V_2 $$

だけの速度で動いていることがわかります。後ろ向きの方向(AC方向)の速度をスラグ速度$V_s$と呼び、前向きの方向(AB方向)の速度をメタルジェット速度$V_j$と呼んでみると

$$V_{AB} =V_j = V_1 + V_2$$

$$V_{AC}=V_s = V_1 - V_2 $$

と書けます。

以上がBMPTモデルに基づいたメタルジェット速度とスラグ速度の導出です。

BMPTモデルは円錐中央部に流れ込む速度関係を考えるという単純な仮定で見通しよく成形炸薬で生じるスラグとメタルジェットの速度を求められるのでいいですね。

メタルジェット速度と侵徹深さ

口径6 in.のTNTにライナー径5 in.の鉄ライナーを$\alpha$=45度の角度で埋め込んだときのメタルジェット速度を以下に示します。計算に必要なデータはBallisticsより抜粋しました。

Fig. 4 口径6 in.のTNTにライナー径 5 in.の鉄ライナー($\alpha$=45度)のメタルジェット速度とスラグ速度

メタルジェットの速度は期待したとおり高速であることがわかります。このような速度域の侵徹深さを計算してみると、ベルヌーイの式で期待される、侵徹深さが密度のみで決定される侵徹様式に到達していることが以下の図からわかります**。

Fig. 5 軟鋼に鉄侵徹体を衝突させたときの侵徹深さ(修正ベルヌーイの式により計算)

結局ノイマン効果ってなんだ

と言われれば、現象だけを見ればやはり、
炸薬に叩き出された円錐型のライナーが円錐中心部で真正面に高速で打ち出されることで高速なメタルジェットを形成し、そのメタルジェットが深い侵徹孔を形成するという現象がノイマン効果と言えるでしょう。
そしてなぜそのような現象が生じるか、というのは、まさに本記事で議論した内容が最も基礎にあると考えていいと思います***。

本記事の内容は以上です。

*炸薬をTNT、ライナー材料を鉄とした場合、$V_2$はたかだか2km/s程度です。
**ここでは形成されるメタルジェット長さや実際の侵徹深さの議論は一度置いておいて、メタルジェットくらいの速度域での侵徹様式にのみ注目しています。
***本モデルや個々から発展したモデルはスラグとメタルジェットの質量比やメタルジェットの長さについても議論できるんですが筆者がそこまで興味がないのでこのへんで。

自作コンテナクラスでboost odeintのcontrolled stepperを使うメモ書き

先日

えすだぶさんの高射砲の本が面白くて、どれくらい当てるの難しいんだろうかと思ったのでざっくり書いて見たけど、全然あたんないですね。高射砲の初速度1000 m/sで、グラフの左から右にフラフラしながら飛んでる飛行機を撃ち落としたい pic.twitter.com/lvEJQDHBxW

— 鹿部 等 (@c_curve1870) 2019年3月16日

というような計算をしたときに、軌道を計算するためにscipyのodeintを使って計算したところ、思ったより簡単に、ちゃんとした公式の元で積分が行えるのでちょっと感動しました。

で、C95で公開したプログラムもodeintで実装できると助かるなぁと思ったんですが、scipy.integrate.odeintは引数がリストのみなので、Anderson-Walkerモデルみたいに13個もある式を真面目にやろうとすれば絶対どこかでindexを間違える未来が見えます。

そういうわけでscipyのodeintを調べているとboost::odeintというものがあり、こちらは自作のコンテナクラスも使えるということを知りました。自作コンテナクラスが使えるのであれば、メンバーにstd::vectorを持ち、それらのindexに対してsetter, getterを作れば、計算部分を書く時は構造体っぽく使えるので良さそうだな、と思いました(自作のPoint classも使えるんですが、operatorを何個も手打ちしたくなかったので今回はstd::vectorをメンバに持つクラスを作りました)。
ここまでの段階ではC++で常微分方程式：boost::odeint入門およびboost odeintのドキュメントを多く参考にしました。

さて、せっかくodeintのようなパッケージが使えるんですから、Controlled stepperを使って可変時間刻みの元で計算をしたいわけですが、

void lorenz(const state_type &x, state_type &dxdt, const double t)
{
  const double sigma(10.0);
  const double R(28.0);
  const double b(8.0 / 3.0);

  dxdt[0] = sigma * (x[1] - x[0]);
  dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
  dxdt[2] = -b * x[2] + x[0] * x[1];
};
using namespace boost::numeric::odeint;

int main()
{
  state_type x;//git hubにあるサンプルコードに従って作った自作のコンテナクラス
  x[0] = 5.0;
  x[1] = 10.0;
  x[2] = 10.0;
  state_type xold = x;

  BOOST_STATIC_ASSERT(is_resizeable::value == true);

  using base_stepper_type = runge_kutta_dopri5;
  auto Stepper = make_controlled(1e-6, 1e-6, base_stepper_type());

  double t = 0.;
  double dt_init = 1e-2;
  double dt = dt_init;
  double dt_max = 1.;
  double tmax = 10.;
  boost::numeric::odeint::controlled_step_result res;
  while (t < tmax)
  {

    xold = x;
    res = Stepper.try_step(lorenz, x, t, dt);
    if (res)
    {
      res = Stepper.try_step(lorenz, x, t, dt);
    }
    else
    {
    }
    if (dt > dt_max)
    {
      dt = dt_max;
    }
    std::cout << x.x() << x.y() << x.z() << std::endl;
  };
}

みたいなことをやると、
/usr/include/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp:89:40: error: cannot convert 'boost::numeric::odeint::norm_result_type<AwState, void>::type {aka AwState}' to 'boost::numeric::odeint::default_error_checker<double, boost::numeric::odeint::range_algebra, boost::numeric::odeint::default_operations>::value_type {aka double}' in return
         return algebra.norm_inf( x_err );
以下、大量のエラーを吐いて爆死します。で、これを回避したい、というのが今回の話題。
C++全然わからないマンなりにいろいろ試した結果、git hubにあるサンプルコード

template< size_t MAX_N >
class my_vector
{
    typedef std::vector< double > vector;

public:
    typedef vector::iterator iterator;
    typedef vector::const_iterator const_iterator;
 
public:
    my_vector( const size_t N )
...

を

template< size_t MAX_N >
class my_vector
{
    typedef std::vector< double > vector;
public:
    typedef vector::iterator iterator;
    typedef vector::const_iterator const_iterator;
    typedef vector::value_type value_type; //ここを追加

public:
    my_vector( const size_t N )
...

とすればなんの問題もなく通るようになりました。
cannot convert 'boost::numeric::odeint::norm_result_type<AwState, void>::type {aka AwState}' to 'boost::numeric::odeint::default_error_checker<double, boost::numeric::odeint::range_algebra, boost::numeric::odeint::default_operations>::value_type {aka double}' in return
         return algebra.norm_inf( x_err );
よく見たら、typeをvalue_typeに変換できないって書いてましたね。
そんなメモ書きでした。

Indesignで隠しノンブルを一気に入れる時のスクリプト

Indesignで、目次までのページ番号をローマ数字で振って、それ以降は普通の数字でページ番号を振りたいということは多々あると思います。
一方で、印刷所に出すときは一貫したノンブルを隠しノンブルとして入れないといけないんですが、ページ数が多いと結構面倒です。
そういうわけで一貫した隠しノンブルを自動で挿入したいんですが、意外にそういうスクリプトを見つけられなくて困りました。
というわけで備忘録がてらおいておきます。

pdfに出力するぞ、というタイミングで使うのがオススメです

PageObj = app.activeDocument.pages;
pStyle = app.activeDocument.paragraphStyles.item("Normal");
a = [Justification.rightAlign,Justification.leftAlign] //右綴じの場合。左綴じは中身を入れ替える
for (i=0; i<PageObj.length; i++)
{
        total = PageObj[i].textFrames.add();
        w = app.activeDocument.documentPreferences.pageWidth;
        h = app.activeDocument.documentPreferences.pageHeight;
        top = h-8
        if (i==0) {
            x=3
            }else{
                x=w
            }

     
       total.visibleBounds = [top.toString() + " mm",(x+(Math.pow(-1, i+1))).toString()+" mm",h.toString()+" mm",(x+(Math.pow(-1,i+1)*4*4)).toString ()+"  mm"];
       total.contentType = ContentType.textType;
       total.contents = (i+1).toString();
     

       total.paragraphs[0].applyParagraphStyle(pStyle);
       total.paragraphs[0].justification = a[i%2];
}

スペース順砲雷撃戦サークル一覧

舞鶴砲雷撃戦が近づいてきましたが、サークル一覧がスペース順じゃなくて見にくいですね。僕は見にくいです。というわけでpythonで適当に作ったので、サークル一覧のページのhtmlを落としてからご活用ください(いいのか？

import re
with open("落としてきた.html", "r", encoding="euc_jp") as f:
    data = f.read()
patt = r'\s*<table .*?SP-No.([A-Z]-[0-9/]*).*?</table>'
table_patt = r'\s*<table.*</table>'
spaces = list(re.finditer(patt, data, re.S))
table = re.search(table_patt, data, re.S)
spaces_sorted = sorted(spaces, key=lambda x: x.group(1))
sspace_concat_str = "\n".join([_.group() for _ in spaces_sorted])
data = data.replace(table.group(), sspace_concat_str)

data = data.replace('src="./',
                    'src="https://sdf-event.sakura.ne.jp/mform/cutfile/')

with open("sorted_蓬莱.html", "w", encoding="euc_jp") as f:
    print(data, file=f)