2019年9月1日日曜日

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

ノイマン効果ってなんだ


いつの間にか、ブログ記事の多くが終末弾道の話になっていて、自分の出自と無関係な方向にどんどんと進んでいるな、という感じがあります。
まぁ、それはともかくとして、終末弾道を考えるとき簡単の場合には侵徹体はバルクと考えてよく、侵徹体は砲から発射されたままの形状を初期条件として侵徹を考えることができます。このような例として、低速度の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程度です。
**ここでは形成されるメタルジェット長さや実際の侵徹深さの議論は一度置いておいて、メタルジェットくらいの速度域での侵徹様式にのみ注目しています。
***本モデルや個々から発展したモデルはスラグとメタルジェットの質量比やメタルジェットの長さについても議論できるんですが筆者がそこまで興味がないのでこのへんで。

2019年4月7日日曜日

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

先日

というような計算をしたときに、軌道を計算するために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に変換できないって書いてましたね。
そんなメモ書きでした。

2019年3月3日日曜日

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];
}

2019年2月8日金曜日

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

舞鶴砲雷撃戦が近づいてきましたが、サークル一覧がスペース順じゃなくて見にくいですね。僕は見にくいです。というわけで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)

2019年2月1日金曜日

APFSDSのRHA侵徹深さを簡単に見積もりたい話

お久しぶりです。
この間コード公開した記事アップロードばっかりだから~と高をくくっていたら前回の記事は去年の3月でした。
いろいろ恐怖です。

さて、今回は表題の通りAPFSDSのRHA侵徹深さを簡単に見積もるにはどうするか、という話です。
それぞれの侵徹体の寸法と速度がわかれば前回の記事に基づいて割合いいところまで求まるんですが、それを毎回はしたくない、というのが正直なところです。

というわけでスペックが与えられたときにぱっと概算できるような指標があるといいですね、というのが今回の記事です。結論から言えば、侵徹深さ$DoP$を求めるには断面積あたりの運動エネルギー$KED$を求めて、

$$DoP[\mathrm{mm}] = 43 KED[\mathrm{kJ/mm^2}] +10$$

くらいの関係を使うのよさそうです。

簡便な指標

侵徹深さが何に基づくか、というのは悩ましい問題です。例えば、以下に横軸に運動エネルギー、縦軸に侵徹深さを取った図を示してみます。




図中のプロットは実弾薬のデータを適当に拾ってきたもので*、図中の線は、侵徹体のLD比をいくつか設定して、与えられた運動エネルギーの中で侵徹深さを最大化させる条件での侵徹深さを計算で求めたものです。
運動エネルギーが大きくなるに連れて侵徹深さは大きくなる傾向がありますが、一方で実弾薬のプロット点とはまるで相関が見られません。
この理由は大雑把に、以下の模式図から説明されます。

このことから、運動エネルギーは、侵徹孔全体を侵徹するのに必要なエネルギーを説明しているのであり、侵徹深さ以外にも断面積が寄与するということです(ということは侵徹孔体積と運動エネルギーはいい相関がある、ということですね、と書いていて気が付きました)。

しかし、今は侵徹深さだけに興味があるので、運動エネルギーを侵徹体断面積で割った値を横軸に取って見ると、非常にいい相関が見られます。運動エネルギーと直径がわかればいいので、大変便利ですね。



これを線形近似して、傾きを求めると、冒頭の

$$DoP[\mathrm{mm}] = 43 KED[\mathrm{kJ/mm^2}] +10$$

が得られます。(切片を0にしたほうが物理的な意味はありますが
ところで、運動エネルギーを断面積で割るってことは直径0の侵徹体があると無限に侵徹できるってこと?ともしかしたら思われるかもしれません。が、残念ながらそうはなりません。別に太くてもいいですし、細くてもいいです。
というのも、運動エネルギーは$\frac{1}{2}mV^2$ですが、$m$には侵徹体寸法と密度が関与する余地があります。そうすると、運動エネルギー$KE$と断面積あたりの運動エネルギー$KED$は

$$KE = \frac{1}{2}mV^2 = \frac{1}{2}\frac{\rho \pi LD^2}{4}V^2$$
$$KED = \frac{4KE}{\pi D^2} = \frac{2}{\pi}\rho LV^2 $$

となり、RHA侵徹深さを求めるには、侵徹体の材質、侵徹体長さ、速度があればいい、ということになります。簡単ですね。

もちろん、砲に要求される能力は上の式に断面積を掛けたものになるので、高L/D比に伴うDの低下から生じる問題を解決しようと思えば、それなりの運動エネルギーが必要になります。そういう意味で大口径は七難隠す、と言えるかもしれませんね(1月頭くらいにツイッターでそんな話をしてました)。
その辺の問題を与えられた口径の範囲で調整した結果が今ある多くの侵徹体の形、と考えるのが自然な形かもしれません。
そういうある種の最適化問題に興味はあるんですが、侵徹体径の効果を考える問題はいろいろ難しくて手が回りません。

逆に何が言えるか

今度はここから逆に何が言えるかを考えてみると、実のところほとんど言えません。それは、複合装甲ではあれこれがどうこうで…とかなんだか難しそうな話ではなく、シンプルにするために様々な単純化を行っているためです。
しかし、例えば装甲材料の強度を変化させたとき、侵徹深さにどのような影響を与えるか、といったことは妥当性を持って言えるでしょう。
というわけで適当に標的強度を変化させたときの侵徹深さの依存性を求めると、以下のようになります。強い装甲は強い、ということですね。



とはいえ、適当な侵徹体が与えられたときに、RHA換算でどれくらいかを見積もりたいときに便利な式だと思います。

それだけ?と言われるとこれだけなんですが、まぁでもインターネット見ててもなかなか見つからなくて困った(どこかの軍事評論家がこの評価をしてる海外の記事**から侵徹体の諸元を書いた表だけ引っ張ってましたね。肝心の部分は無くて残念でした。)のと、一番基礎のこういうところの話しないままに先端的なところを見てもよくわからないので整理する意味で、といったところです。

まぁむにゃむにゃ難しそうな横文字とかいろんな用語を並べたら、読んで楽しい侵徹小説を書いてみたりいろいろできるとは思うんですが、自分はそういうやりかたで何かの役に立つような定量的な評価を導き出せないので、本ブログでは簡単なところからゆっくりやっていけたらなと思います。

* 侵徹深さと侵徹体のデータはWikipedia及び、 
http://www.inetres.com/gp/military/cv/weapon/M256.html,https://en.wikipedia.org/wiki/M829 
https://www.orbitalatk.com/defense-systems/armament-systems/120mm/ 
http://sturgeonshouse.ipbhost.com/topic/1086-tanks-guns-and-ammunition/
https://www.researchgate.net/profile/Mariusz_Magier/publication/317409282_ANALYSIS_OF_ENERGETIC_PARAMETERS_FOR_ANTITANK_KINETIC_AMMUNITION_OF_CONTEMPORARY_BATTLEFIELD/links/5a1bf02ea6fdcc50adeca812/ANALYSIS-OF-ENERGETIC-PARAMETERS-FOR-ANTITANK-KINETIC-AMMUNITION-OF-CONTEMPORARY-BATTLEFIELD.pdf
あたりを参考にしていますが、実際の侵徹体についてはあんまり詳しくないので明らかにおかしい点があれば指摘いただければ幸いです。

** https://www.researchgate.net/profile/Mariusz_Magier/publication/317409282_ANALYSIS_OF_ENERGETIC_PARAMETERS_FOR_ANTITANK_KINETIC_AMMUNITION_OF_CONTEMPORARY_BATTLEFIELD/links/5a1bf02ea6fdcc50adeca812/ANALYSIS-OF-ENERGETIC-PARAMETERS-FOR-ANTITANK-KINETIC-AMMUNITION-OF-CONTEMPORARY-BATTLEFIELD.pdf

2018年3月29日木曜日

ATモデル及びAWモデルで侵徹挙動を見るコードの公開

前回のちょっと頑張ってL/D効果も含めた侵徹深さの評価をしたい話と、ベルヌーイの式とエロージョンを伴う侵徹の覚え書きで書いていた内容を割と簡単に扱えるような形に書き下すことができた(と思った)ので、大したものではないですが公開しておきます。

ダウンロードはこちらからどうぞ。

環境としては、Windowsの方はAnaconda入れたら動くと思います。
(実際に使っているのはPandas, numpy, matplotlibの3つです)

jupyter notebookで使い方.ipynbを開くと一連の操作を実際に動かすことができます。
他の環境の方は、src内でpython setup.py build_ext -iを打てばその環境でのpydができると思うんですが、すいません、こういう配布には詳しくないのでわかりません。適当にぐぐってください。

あとdocstringは適当です。

計算ルーチンは間違ってはないと思うんですが、計算ルーチン部分に不備があれば、その箇所をコメントに書いていただければ対応します(適当に直して使ってもらっていただいても大丈夫です)。

しっかりしたプログラムに関する教育を受けたことがないのでプログラムが下手なのは多めに見ていただければ幸いです。

再配布したい人はいないと思いますが、一応再配布はしないでいただけると幸いです。

3/30追記
$\alpha$を計算するルーチンで参考にした論文にミスがあったことがわかったので修正しています。

2018年3月25日日曜日

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


鹿部です。
久しぶりにTerminal ballisticsを読んでいて、ふと、趣味でAPFSDSとかを調べている人にとってAnderson-Walkerモデルがすごくいいモデルなんじゃないかと思い先週からちまちま調べ物をしていました。


今回はその辺のネタです。

この記事でできたこと

最初に、今回の記事の要点を以下に示します。
以下にAWモデルを用いてLD比を3から100まで変えて侵徹深さと速度の関係を計算した図を以下に示します。
パラメーターとして、侵徹体はタングステンを想定した物性値を用い、$\sigma_p$として2GPaを用い、標的は鋼を想定した物性値を用い、$Y_t$として1GPaを用いました。
ATモデルについて計算する際には$R_t$を5.5GPaとし、それ以外は同じにしています。


このようにLD比によって侵徹深さは大きく異なります。このことについて考えるために、衝突速度2000 m/sで衝突した際の侵徹中の侵徹速度(破線)と侵徹体後端速度(実線)の推移を以下に示します。




図から明らかなように、LD比が大きい時、侵徹速度は速やかに低下し、ATモデルとほとんど同じ侵徹速度と後端速度の推移を示します。
一方で、LD比が3と小さいとき、侵徹速度はATモデルの侵徹速度に漸近するものの、全過程の20%程度まではATモデルの侵徹速度よりも高い値を示しています。
このため、低LD比な侵徹体は高LD比な侵徹体に比べて高P/Lを示したと言えます。
また、遷移過程が一瞬で終了する高L/D比な侵徹体は、侵徹の殆どの領域に渡ってATモデルから計算される侵徹速度と同じ値を示しており、ATモデルの有効性が改めて確認できます。また、これらのことから、ATモデルはLD比が∞の侵徹体の侵徹過程を示しているとざっくりみなすことができます。

せっかくなので実用の弾薬とRHA換算侵徹深さを比較した時の結果も見てみます。破線上にプロットが乗っていればいるほど計算と実際がよく一致していると言えます*****。



割といい線行ってるのでは、と思いました。


趣味で調べてる人にとってのATモデルの困難さ

本ブログでも何度も取り上げた(1,2)Alekseevskii-Tateモデル(修正ベルヌーイの式: ATモデル)は非常に優れたモデルであり、ATモデルに基づいて多くの研究がなされています。

以前記事にしたように、ATモデルは以下の式で表されるように侵徹体と標的界面での圧力の釣り合いを考え、

\[ Y_p +\frac{1}{2} \rho_p (V-U)^2 = \frac{1}{2}\rho_t U^2 + R_t \]

そこから侵徹速度を弾芯後端速度の関数として求め、

\[ U = \frac{1}{1-\mu^2} \left(V - \mu \sqrt{ V^2 +\frac{2(R_t-Y_p)(1-\mu^2)}{\rho_t}}\right) \]

同時に侵徹体後端での運動方程式

\[ -Y_p = \rho l \frac{dv}{dt} \]

と合わせて以下の積分を行うことで侵徹深さを求めるモデルでした。

\[ P = \int^t_0 u dt = -\frac{\rho_p}{Y_p} \int^0_{V_0} uldV = \frac{\rho_p}{Y_p} \int_0^{V_0} u\ l\ dV \]

しかしその一方で、普通に解くと運動方程式を時間ではなく侵徹体長さで置換して解くため、

- 侵徹体寸法を考慮できない。
- 侵徹初期及び終末期の挙動の表現が十分でない。
- 標的強度項$R_t$の取り扱いが難しい

といった問題があります。侵徹体寸法が考慮できないというのは、例えば侵徹体のL/D比などは考えられないということですが、実際には以下の図のようにL/D比は侵徹深さP/Lに強く影響します*。

侵徹初期及び終末期の表現が十分でないというのは以下の図のように侵徹初期の減速過程を考慮できず、また侵徹終末期の減速を過小評価してしまう、といった点です*。



私達のような素人にとって、標的強度項$R_t$が標的強度と大きく異なることは侵徹について考える上で大きな障害です。
例えば、タングステン侵徹体と高強度鋼標的の侵徹深さをATモデルで解析するとき、侵徹体強度$Y_p$には侵徹体の動的強度である2GPaをそのまま適用できるのに対して、標的強度$R_t$はおよそ5.5GPa程度の値が用いられます**。

以前の記事で書いたように、ユゴニオ弾性限界はポアソン比と降伏強度で得られますが、降伏強度1GPaの鋼でのユゴニオ弾性限界は2GPaであり、$R_t$の5.5GPaとはまだまだ隔たりがあります。$R_t$がこのような大きな値を取る理由は、侵徹体/標的界面で侵徹が進行する際、周囲の塑性変形していない部分からの拘束が加わるためです***。
このように、標的強度$R_t$というのは材料強度と相関がありつつも適切な値を定めるのが難しいパラメーターです****。

で、ATモデルのこういう困難から、個人的にこれ以上弄ることの難しさを感じていました。

AWモデル


Anderson-Walkerモデル(AWモデル)はATモデルと同じく一次元の解析的な侵徹モデルです。
前述のような理由から、ATモデルを使って踏み込んで考えるのは難しいなあと考えていたわけですが、一方でAWモデルは侵徹体と標的界面での運動量の釣り合いを考えることで、侵徹速度についても素直に時間発展を考えることができるようにしたモデルで、これは楽しめるんでないか、と思った次第です***。
とはいえ、自分も***の記事とTerminal ballisticsを読み込んでいるわけではないので、モデルの十分な説明はここでは行なえません。(そのうち読んでどこかでまとめると思います)
しかしながら、計算自体は比較的簡単に行えるので、ここで簡単に紹介し、AWモデルに基づいて手元で行った計算結果を示します。

計算に必要な式はChocron et al., Int. J. Impact Eng. 28(2003), pp. 391と***の2つの記事を合わせて読むとわかりやすいと思いますが、大本となる式は

\[ \rho_p \dot{v}(L-s)+\dot{u}\left\{\rho_ps+\rho_tR\frac{\alpha-1}{\alpha+1}\right\}+\frac{\rho_ps^2}{2}\frac{d}{dt}\left(\frac{v-u}{s}\right)+\rho_t\dot{\alpha}\frac{2Ru}{\left(\alpha+1\right)^2}=\frac{1}{2}\rho_p(v-u)^2-\left\{\frac{1}{2}\rho_tu^2+\frac{7}{3}Y\ln{\alpha} \right\} \]

です。この式のうち、現時点ではわからない項目をどうにかして決定して、$\dot{u}=$の形に整理すれば時間に対して顕に侵徹速度を求められるということになります。シンボルは元記事を適当に参照していただければ。

一番簡単に決めれるのは$L$で、これは$\dot{L} = -(v-u)$から簡単に求めることができます。$\dot{v}$は***中でいくつかの過程を置くことで

\[ \dot{v} = -\frac{\sigma_p}{\rho_p(L-s)}\left\{1+\frac{v-u}{c}+\frac{\dot{s}}{c}\right\} \]

として得られます。$R$は侵徹体半径$R_p$を用いて実験結果にフィッティングする形で、衝突速度$V_0$ km/sと以下の式で与えられています(ここではW/Feの組み合わせだったはずです)。

\[ R = R_p\left(1+0.287V_0+0.148V_0^2\right) \]

$\alpha$はCavity expansion analysisを使うことで定式化できるらしいのですが、求めるのはちょっと厄介で、

\[ \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$をニュートンラフソン法なりで適当に求めることで得られます。左辺平方根中を見れば明らかなように、$\alpha$が実数解を取りうる範囲は、$(-\sqrt{\frac{K_t}{\rho_t u^2}},\sqrt{\frac{K_t}{\rho_t u^2}})$で、更に***の記事Fig. 6から$\alpha>1$であることは明らかなので、$(1,\frac{K_t}{\rho_t u^2})$の範囲で解を求めればいいと思います(ただ、関数としては単調減少関数なのでそこまで厄介ではないと思います)。

ここでの$K_t$はある応力で圧縮されている標的の体積弾性率で、

\[ K_t = K_0\left(1+k\frac{u}{c_0}\right)^2 \]

として表せます。$k$は衝撃波速度を$U_s = c_0+k u_p$として表した時の$k$です。少し補足すると、標的界面の粒子速度というのは結局侵徹速度なので、カッコの中身$(1+k\frac{u}{c_0})$というのは$(c_0+ku)/c0=U_s/c0$で、$c_0=\sqrt{K_0/\rho_t}$です。すると上の式は

\[ K_t = \rho_t U_s^2 \]

あるいは、

\[ U_s = \sqrt{K_t/\rho_t} \]

となりますから、界面ではこういう衝撃波が標的内部に生じているということになります(あるいはその衝撃波で圧縮された部分の体積弾性率を求めているということ?)。

ここまで求まると、$s$は以下の式に

\[ s = \frac{R}{2}\left(\frac{v}{u}-1\right)\left(1-\frac{1}{\alpha^2}\right) \]

それぞれ代入して得られます。
残るは$\frac{d}{dt}\left(\frac{v-u}{s}\right)$ですが、これは***中で

\[ \frac{d}{dt}\left(\frac{v-u}{s}\right)=\frac{4}{(\alpha^2-1)R}\left\{\frac{\alpha^2\dot{u}}{2}-\frac{u\alpha\dot{\alpha}}{\alpha^2-1}\right\} \]

という形で与えられていますので、これを使います。$\dot{\alpha}, \dot{s}$だけは前の計算からの差分で与えるしか無いと思います。
後はこれらを適当な時間刻みで回せばいいですが、C++で適当に書いただけでも時間刻み1e-7 sで1msまで計算しても実行した瞬間には終わってる(5msくらい)ので、そんなにストレスなく回せると思います。また、時間について回す一方で、u<=0, L<=0で侵徹は終了しているので、その辺に気をつけるといいと思います。
ここで重要なのが、AWモデルでは、初期条件がいずれも材料強度と物性値で、ATモデルで私達を苦しめた$R_t$のようなパラメーターが消え、その他のパラメーターはモデルから導かれているという点です。
モデルの簡略化の限界を除けば、AWモデルは私達に材料強度、物性値と侵徹過程についていい見通しを与えてくれるように思います。

まず、以下にAWモデルを用いてLD比を3から100まで変えて侵徹深さと速度の関係を計算した図を以下に示します。パラメーターとして、侵徹体はタングステンを想定した物性値を用い、$\sigma_p$として2GPaを用い、標的は鋼を想定した物性値を用い、$Y_t$として1GPaを用いました。ATモデルについて計算する際には$R_t$を5.5GPaとし、それ以外は同じにしています。

このようにLD比によって侵徹深さは大きく異なります。このことについて考えるために、衝突速度2000 m/sで衝突した際の侵徹中の侵徹速度(破線)と侵徹体後端速度(実線)の推移を以下に示します。


図から明らかなように、LD比が大きい時、侵徹速度は速やかに低下し、ATモデルとほとんど同じ侵徹速度と後端速度の推移を示します。一方で、LD比が3と小さなとき、侵徹速度はATモデルの侵徹速度に漸近するものの、全過程の20%程度まではATモデルの侵徹速度よりも高い値を示しています。このため、低LD比な侵徹体は高LD比な侵徹体に比べて高P/Lを示したと言えます。また、遷移過程が一瞬で終了する高L/D比な侵徹体は、侵徹の殆どの領域に渡ってATモデルから計算される侵徹速度と同じ値を示しており、ATモデルの有効性が改めて確認できます。

また、これらのことから、ATモデルはLD比が∞の侵徹体の侵徹過程を示していると考えることができます。

AWモデルの際立った特徴は以上ですが、せっかくなので実用の弾薬とRHA換算侵徹深さを比較した時の結果も見てみましょう。破線上にプロットが乗っていればいるほど計算と実際がよく一致していると言えます。

割といい線行ってるのでは、と思うのですがどうでしょうか。


白状すると、この計算では標的の$Y_t$を0.7GPaにしていい感じに合うような値を用いています(侵徹体強度は2GPaで揃えています)。

今回は以上です。
ありがとうございました。

(そういえばまだ他の記事の計算結果と摺り合せたValidation?みたいなのをしていませんでした。今度して変だったら訂正します。)


ところで当然ですが、このモデルももっと厳密な数値計算モデルの結果と比較すると異なるところがある点、傾斜した板に対する評価はできない点など、色々難しいところもあるので3割引くらいの目で見たいところです。


*Z. Rosenberg, E. Dekel, Terminal ballistics, (2016), Springer
**例えばJA. Zukas, High velocity impact dynamics (1990)
***JD. Walker, CE. Anderson, Int. J. Impact Eng. 16(1995), pp. 19-48
**** $R_t$と材料強度の間の関係はHillのCavity expansion analysisを用いることで表せますが*、ちょっと微妙さが残るのと、他の困難は解決できません。
***** 侵徹深さと侵徹体のデータはWikipedia及び、
http://www.inetres.com/gp/military/cv/weapon/M256.html,https://en.wikipedia.org/wiki/M829
https://www.orbitalatk.com/defense-systems/armament-systems/120mm/
http://sturgeonshouse.ipbhost.com/topic/1086-tanks-guns-and-ammunition/
https://www.researchgate.net/profile/Mariusz_Magier/publication/317409282_ANALYSIS_OF_ENERGETIC_PARAMETERS_FOR_ANTITANK_KINETIC_AMMUNITION_OF_CONTEMPORARY_BATTLEFIELD/links/5a1bf02ea6fdcc50adeca812/ANALYSIS-OF-ENERGETIC-PARAMETERS-FOR-ANTITANK-KINETIC-AMMUNITION-OF-CONTEMPORARY-BATTLEFIELD.pdf
あたりを参考にしていますが、実際の侵徹体についてはあんまり詳しくないので明らかにおかしい点があれば指摘いただければ幸いです。