お手軽に高速度侵徹を解析したい話

夏コミも終わって小休止と言った感じであんまり本を読めていません。
鹿部です。


今回はベルヌーイの式とエロージョンを伴う侵徹の覚え書きで得た結果から実際に計算するプログラムを作る話です。夏コミの本で行った解析は今回のコードを基本にしたコードで行われていて、今回のコードを使えば高L/D比な弾芯による高速度侵徹のモデルを大まかに扱えます。
基礎的なモデルは当該記事を参照のこと。

ちなみに、Solution of the Long Rod Penetration EquationsにはFortran版のソースコードが載っていますので、今回の話はあんまり目新しいものはありません。

何で書くか

今回はPythonで書いてます。理由としては、自分がPythonしか使えないから最後に出てくる数値積分を簡単に扱うためです。

前提条件

Python3、Numpy、Scipyが使えることを前提としています。

作図とか色々考えるとAnacondaを入れると快適だと思います。

計算をする上で、今回はベルヌーイの式とエロージョンを伴う侵徹の覚え書きのY,Rともに強度がある時の部分を実装しています。それ以外は解析的に簡単に求められるので省略します。

解きたいもの

今回は、
$$P =  \frac{\rho_p}{Y_p} \int_0^{V_0}U\ l\ dV$$
という積分の式を、
侵徹速度$U$と侵徹体後端速度$V$の関係
$$U = \frac{1}{1-\mu^2} \Biggr(V - \mu \sqrt{ V^2 +A}\Biggl)$$
と、侵徹体が初速$V_0$から$V$まで減速したとき、どれくらい消耗しているかを示す式
$$\frac{l}{L} = \Biggr( \frac{ V+\sqrt{V^2+A}}{V_0+\sqrt{V_0^2+A}} \Biggl)^{\bigr(\frac{\rho_p \mu A}{2(1-\mu^2\ )}\ \ \ \bigl)} \exp \Biggr( \frac{\rho_p \mu}{2Y_p(1-\mu^2)} \biggr( V\sqrt{V^2+A} - \mu V^2 - \bigr[ V_0\sqrt{V_0^2+A} - \mu V_0^2 \bigl] \biggl) \Biggl)$$
から解きます。$A,\mu$はそれぞれ$A=\frac{2(R_t-Y_p)(1-\mu^2)}{\rho_t},\mu = \sqrt{\frac{\rho_t}{\rho_p}}$です。

$R,Y,\rho$はそれぞれ標的強度、侵徹体強度、密度で、添字のtとpは標的と侵徹体を指しています。

まずは書いてみる

必要なモジュールとしてとりあえずnumpyとscipyだけimportしておきます。

import numpy as np
import scipy

今の侵徹体後端速度$V$から侵徹速度$U$を求める関数は簡単に

def calc_u(V,  mu, A):
    u = 1./(1.-mu**2.)*(V-mu*np.sqrt(V**2.+A))
    if u != u:#nan返したら0にしとくと無難です。
        u = 0
    elif u < 0.:#u<では侵徹は停止しているので0にしときます。
        return 0.
    return u
 
def calc_A(mu, Rt, Yp, rhot):
    return 2.*(Rt - Yp)*(1.-mu**2.)/rhot
 
def calc_mu(rhot,rhop):
    return np.sqrt(rhot/rhop)

と書けます。$A,\mu$を計算する関数も一緒に書いておきました。
次に侵徹体後端速度と侵徹体残存長さ$l/L$を求める式を書き下してみます。今、関数の形に注目してみると、侵徹体残存長さ$l/L$というのはある関数$f(v)$を用いると

$$\frac{l}{L} = \frac{f(V)}{f(V_0)}$$

として書かれることがわかります。式を見比べてみると、関数$f(v)$というのは

$$f(v) = \left(v+\sqrt{v^2+A}\right)^{\frac{\rho_p \mu A}{2(1-\mu^2\ )}} \ \ \ \exp\left(\frac{\rho_p \mu}{2Y_p(1-\mu^2)}\left( v\sqrt{v^2+A} - \mu v^2\right) \right)$$

なので、これを素直に書き下せば、侵徹体残存長さ$l/L$というのは関数calc_lのように得られます。

def f(V,  mu, A,  Yp, rhop):
    return (V+np.sqrt(V**2.+A))**(rhop*A/2./(1.-mu**2.)*\
            np.exp(rhop*mu/2./Yp/(1.-mu**2.)*(V*np.sqrt(V**2.+A)-mu*V**2.))
 
def calc_l(V, V0, mu, A, Yp, rhop):
    #f(V)/f(V0)を計算
    return f(V, mu, A, Yp, rhop)/f(V0, mu, A, Yp, rhop)

ここまでくればあと少しで、これを
$$P =  \frac{\rho_p}{Y_p} \int_0^{V_0} u\ l\ dV$$
に突っ込めば侵徹長さが出ます。後は$U$と$l$の積をscipyのquad関数に投げればよくて、

def ul(V, V0, mu, A, rhop, Yp):
    return calc_l(V, V0,mu,A,rhop)*calc_u(V, mu, A)*rhop/Yp
def calc_P(V0, mu, A, rhop, Yp):
    a = scipy.integrate.quad(ul, 0, V0, args=(V0, mu, A, rhop, Yp))
    return a[0]

として、侵徹深さ$P/L$を求める事ができます。
侵徹体としてタングステン、標的として高強度鋼を想定して、$\rho_p = 17000 kg/m^3,\ \rho_t = 7900 kg/m^3,\ Y_p = 2\times10^9 Pa,\ R_t = 6\times10^9 Pa$として初速$V_0$を0m/sから2500m/sの間で変化させながら計算してみると、以下のような図を得ます。

これはタングステン弾芯と高強度鋼標的の衝突ではありがちな曲線です。

この形の問題点

この形の問題点は簡単で、今は侵徹体が標的より軟らかい時のみを考えており、$Y_p \gt R_t$では正しい値を返しません。というわけで$Y_p\gt R_t$な時の状況を考えてみましょう。

$Y_p \gt R_t$の時の解 

この時の解は前回は与えていませんが、C92で頒布したWeekend ballisticsでは説明しているので、今回はその説明を使って進めます。

$Y_p \gt R_t$の時、侵徹は以下の図のように進行します。

ここで、$V^{\prime}_{lim}$はエロージョンによる侵徹が止まる速度で、$V^{\prime}_{lim}$以下では一般的な徹甲弾と同じように弾体形状を保ったまま侵徹が進行します。そこで、計算を二段階に分けます。

1. 侵徹体後端速度が$V\gt V^{\prime}_{lim}$である時、侵徹体は消耗しながら侵徹し、$V^{\prime}_{lim}$で侵徹体長さは$l_r/L$だけ残して消耗が止まる。

2. $V^{\prime}_{lim}$以下の速度では侵徹体は消耗せず侵徹を行う。

2.のステップについては運動方程式を解くことにより、

$$P/l_r = \frac{1}{\mu^2}\log \left(1+\frac{\rho_t V^{\prime 2}_{lim}}{2R_t}\right)$$

として簡単に与えられます。

$V^{\prime}_{lim}$を計算するのは簡単で、ベルヌーイの式とエロージョンを伴う侵徹の覚え書きの圧力の釣り合いの式で$V=U$になるような$V$を求めればよく、

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

から、

$$V^{\prime}_{lim} = \sqrt{\frac{2(Y-R)}{\rho_t}}$$

で与えられます。また、$V^{\prime}_{lim}$での侵徹体残存長さ$l_r/L$はcalc_lの引数Vを$V^{\prime}_{lim}$とするだけで大丈夫です。そうすると関数calc_Pは

def calc_vlim(rhot, Yp, Rt):
    return np.sqrt(2.*(Yp-Rt)/rhot)
 
def calc_P(V0, mu, A, rhop, rhot, Yp, Rt):
     if Yp =< Rt:
         a = scipy.integrate.quad(ul, 0, V0, args=(V0, mu, A, rho, Yp))
         return a[0]
    elif  Yp > Rt:
        vc = calc_vlim(rhot, Yp, Rt)
        lr = 1
        b = 0.
        if V0 > vc:
            lr = calc_l(vc, V, mu, A, Yp, rhop)
            a = scipy.integrate.quad(ul, vc, V0, args=(V0, mu, A, rhop, Yp))
            b += a[0]
        elif V0 <= P.vc:
            vc = V0
        b += lr*np.log(1 +0.5*rhot*vc**2./Rt)*1/mu**2.
        return b

のような形に書き換えれば良いことになります。これを用いて侵徹深さを求めたのが、以下の図になります。計算条件は図中の通りで、図中丸印は実験結果です。このようにベルヌーイの式を用いることで実験結果を精度良く再現することが可能です。

計算する時に気をつけること

そんなわけで侵徹体と標的の密度と強度さえ設定すれば色々遊ぶことができるこのコードですが、気をつけることが色々あるように思います。

1. 強度はあくまでも計算パラメーター
タングステン弾芯と高強度鋼標的について計算した際、Rtを6GPaとして与えました。これはJ. ZukasのHigh velocity impact dynamicsよりパラメーターを拝借したものですが、高強度鋼は普通このような強度をユゴニオ弾性限界であっても示しません。なので、特に標的強度については材料的な意味を求めるのは中々難しいものがあります。

2. Primary phaseのみを考慮している
高速度侵徹は大きく分けて3つに侵徹過程を分けることができますが、今回のコードではその中で最も主要な役割を果たすPrimary phaseのみを考慮しています。極端な低速度/高速度はあまり見ないほうが良いような気がします。

3. L/D比を考慮することができない
高速度侵徹はL/D比によって侵徹深さが変わるのはよく知られた事実ですが、この現象をうまく整理するのはなかなか難しいものがあります(試みられてはいます)。


こういう事情がありますが、dticなどにある実験結果を用いてフィッティングして、全体の曲線を見積もるなどすると楽しめる気がします。