ベルヌーイの式とエロージョンを伴う侵徹の覚え書き

ユゴニオ弾性限界と来たら次はベルヌーイの式だろうと思って最近ちまちま調べ物していました。

今回はいつも以上に覚え書きです。
また、流体で扱うベルヌーイの式を探して来られた方には合わない内容となっています。

高速度の侵徹(特にHEAT)を議論する時、大前提として侵徹長さは密度比で決定されるとよく言われるように思います。

これを議論する時、しばしばベルヌーイの式が出てきますが、そこからどうやったら出てくるのかがよくわからなかったので適当に式をいじってみるか、というのが今回の記事の目的です。誰かに示すための数式とかもう何年も書いたことないので読みにくいと思います。

ベルヌーイの式

自分は流体の基礎が全く無いのでまず式ありきになってしまうんですが、いわゆる一般的な材料強度を考慮したベルヌーイの式は
$$Y_p +\frac{1}{2} \rho_p (V-U)^2 = \frac{1}{2}\rho_t U^2 + R_t$$
で与えられます。左辺は弾芯が作る圧力、右辺は標的が作る圧力で、$Y_p$は弾芯強度、$\rho_p, \rho_t$は弾芯密度、標的強度、$V,U$は弾芯後端速度(衝突直後は衝突速度)と侵徹速度、$R_t$は標的強度です。$Y_p$,$R_t$が与える影響の議論は今回は省略します。先端が速度$U$で侵徹し、後端は速度$V$で進行することから、弾芯が消耗する速度は以下の式で与えられます。
$$\frac{dl}{dt} = -(V-U)$$
ここで、$l$は弾体の長さ、$t$は時間です。
さて、侵徹の際、先端部のみが塑性変形をし、残部は剛体のように振る舞えば反力によって減速しますが、この際に弾芯後端が受ける応力は(完全弾塑性体であれば)、$Y_p$に等しいことがわかります。そこで、弾芯後端の運動方程式を
$$-Y_p = \rho_p l \frac{dv}{dt}$$
と立てます*。

これが今回使う式です。

弾芯と標的の両方に強度がないケース

一般的に解く前に、特殊な場合について考えます。$Y_p=R_t=0$の時、ベルヌーイの式は

$$\frac{1}{2} \rho_p (V-U)^2 = \frac{1}{2}\rho_t U^2$$
となります。
この時、侵徹速度$U$は単純に解いて
$$U = \frac{V}{1+\sqrt{\frac{\rho_t}{\rho_p}}}$$
が得られます。$Y_p$が0ですから、
$$-Y_p = \rho_p l \frac{dv}{dt}$$
から$V$は一定であり、
$$\frac{dl}{dt} = -(V-U)$$
弾芯の消耗速度は一定であることがわかります。そこで、侵徹の持続時間$\tau$は弾芯が消耗しきるまでの時間ですから、弾芯長さを$L$として、
$$\tau = \frac{L}{V-U}$$
とするのが妥当です。 すると侵徹深さ$P$というのは、時間$\tau$の間に侵徹速度$U$だけ進むことと同じなので、
$$P =U\tau = \frac{UL}{V-U} = \frac{ \frac{VL} {1+\sqrt{\frac{\rho_t}{\rho_p} } }} {\frac{V\sqrt{\frac{\rho_t}{\rho_p}}} {1+\sqrt{ \frac{\rho_t} {\rho_p} } }} = L\sqrt{\frac{\rho_p}{\rho_t}}$$
が得られます。この結論は重要で、理想的には弾芯長さで規格化された侵徹深さ$\frac{P}{L}$は標的と侵徹体の密度の比の平方根$\sqrt{\frac{\rho_p}{\rho_t}}$で一定となることがわかります。

弾芯に強度がないケース

弾芯に強度がない時、これは先程の例と同じように$Y_p$が0ですから、
$$-Y_p = \rho_p l \frac{dv}{dt}$$
から$V$は一定であり、
$$\frac{dl}{dt} = -(V-U)$$
弾芯の消耗速度は一定であることがわかります。ただし、ベルヌーイの式は少し変わって
$$\frac{1}{2} \rho_p (V-U)^2 = \frac{1}{2}\rho_t U^2  +R_t$$
となります。この場合、侵徹体が強度がない一方で、標的には強度があるので、速度が一定の速度まではエロージョンを伴った侵徹を起こさない(侵徹速度$U=0$)ことがわかります。そこで、$U=0$について上式を解けば、
$$V = \sqrt{\frac{2R_t}{\rho_p}}$$
が得られ、衝突速度がこれ以上であれば侵徹が生じることがわかります。この速度を特別に$V_c$と置いておきます。

そこで、$V_c$以上の速度について、侵徹速度$U$を平方完成して求めると侵徹速度$U$は
$$(\rho_p - \rho_t)U-2\rho_p UV + \rho_p V^2 = 2R_t$$
$$(\rho_p -\rho_t) (U-\frac{\rho_p}{\rho_p- \rho_t}V)^2 - \frac{\rho_p^2 - \rho_p (\rho_p-\rho_t)}{\rho_p-\rho_t} = 2R_t$$
から、
$$U =  \frac{ \rho_p V-\sqrt{ \rho_p\rho_t V^2 +2R_t(\rho_p-\rho_t)}}{ \rho_p-\rho_t}$$
が得られます。あるいは、$\mu = \sqrt{\frac{\rho_t}{\rho_p}}$とおけば、
$$U = \frac{1}{1-\mu^2} \Biggr(V - \mu \sqrt{ V^2 +\frac{2R_t(1-\mu^2)}{\rho_t}}\Biggl)$$
と整理することが出来ます。

後は先程と同様に、
$$\tau = \frac{L}{V-U}$$
から、
$$P =U\tau = \frac{UL}{V-U} = \frac{ L \frac{1}{1-\mu^2} \Biggr(V - \mu \sqrt{ V^2 +\frac{2R_t(1-\mu^2)}{\rho_t}}\Biggl) } { \frac{1}{1-\mu^2} \Biggr(-\mu^2 V + \mu \sqrt{ V^2 +\frac{2R_t(1-\mu^2)}{\rho_t}}\Biggl) } = L\frac{\Biggr(V - \mu \sqrt{ V^2 +\frac{2R_t(1-\mu^2)}{\rho_t}}\Biggl)}{ \Biggr(-\mu^2 V + \mu \sqrt{ V^2 +\frac{2R_t(1-\mu^2)}{\rho_t}}\Biggl)}$$
が得られます。

ここで重要なのは、弾芯に強度がない場合でも、一定の速度以上であれば強度のある標的を侵徹することが出来るという点です。この時、最も単純には標的にはユゴニオ弾性限界を超えた圧力がかかっていて、互いにエロージョンを起こしながら侵徹することがわかります**。

弾芯と標的共に強度があるケース

この場合、$Y_p > R_t$と$Y_p < R_t$の場合が考えられますが、エロージョンによる侵徹を考える立場から、$Y_p < R_t$の場合のみを考えます。この場合、ベルヌーイの式は

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

であたえられます。今回は先の2つの例と異なり$Y_p \neq 0$であるため、

$$-Y_p = \rho l \frac{dv}{dt}$$
から、弾芯後端速度$V$は連続的に変化することがわかります。そこで、これについて考慮しつつ考えていきます。
何はともあれ、ベルヌーイの式は$V$と$U$の関係を与えるので、ここから手を付けていくことにします。この場合、$Y_p$を右辺に移してしまえば後は標的にのみ強度がある場合と同じであることは明らかで、
$$U = \frac{1}{1-\mu^2} \Biggr(V - \mu \sqrt{ V^2 +\frac{2(R_t-Y_p)(1-\mu^2)}{\rho_t}}\Biggl)$$
となります。$A=\frac{2(R_t-Y_p)(1-\mu^2)}{\rho_t}$と置いてしまえばシンプルに、
$$U = \frac{1}{1-\mu^2} \Biggr(V - \mu \sqrt{ V^2 +A}\Biggl)$$
が得られます。これで$V$と$U$の関係がわかったので、
$$-Y_p = \rho_p l \frac{dv}{dt}$$
に
$$\frac{dl}{dt} = -(V-U)$$

を突っ込んで解いてやれば解けそうな気がします。そこで素直に突っ込んでやると、

$$Y_p = \rho_p l (V-U) \frac{dV}{dl}$$

が得られるので、片々整理して

$$\frac{dl}{l} = \frac{\rho_p}{Y_p} (V-U) dV.$$

$U$も顕に書けば、

$$\frac{dl}{l} = \frac{\rho_p}{Y_p} \Bigr( \frac{1}{1-\mu^2}\bigr(-\mu^2 V +\mu\sqrt{V^2+A} \bigl) \Bigl) dV$$

が得られるので、これを積分すればいいだろう、ということがわかります。左辺は残った弾体長さ$l$についての積分で、不定積分は

$$I_l = \log l +C_l$$

$V=V_0$のとき、侵徹は始まっていないので$l=L$でしょう。また、任意の速度$V$の時は$l$としておきます。なので左辺は

$$左辺=\log(l/L)$$

になります。右辺の積分は少し手強いです。積分の範囲は衝突速度$V_0$から任意の速度$V$までと取ることにして、まずは不定積分をしてみます。$\frac{1}{1-\mu^2}(-\mu^2 V)$は簡単で、

$$I_0 = \frac{1}{2(1-\mu^2)}(-\mu^2 V) +C_0$$

です。残った$\frac{1}{1-\mu^2}(\mu\sqrt{V^2+A})$の積分は、このページを参考にすると

$$I_1 = \frac{1}{2(1-\mu^2)}\Bigr(\mu\bigr(V\sqrt{V^2+A} + A \log(V+\sqrt{V^2+A})\bigl)\Bigl)  +C_1$$

が得られるので、これらを$V_0$から$V$まで積分すると、

$$右辺=\frac{\rho_p \mu}{2Y_p(1-\mu^2)}\Bigr( V\sqrt{V^2+A} + A \log(V+\sqrt{V^2+A}) - \mu V^2 -\big[ V_0\sqrt{V_0^2+A} + A \log(V_0+\sqrt{V_0^2+A}) - \mu V_0^2 \bigl] \Bigl)$$

が得られます。右辺と左辺をまとめてやると、

$$\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)$$

ここまでくれば侵徹深さを求めるまであと一息です。侵徹深さは先程と同じように侵徹速度$U$の時間積分ですが、このままでは扱いにくいので、$-Y_p = \rho_p l \frac{dv}{dt}$から、$dt = -\frac{\rho_p l}{Y_p}$となるので、

$$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$$

簡潔にまとめて、

$$P =  \frac{\rho_p}{Y_p} \int_0^{V_0} u\ l\ dV$$

を適当に積分することで得られます。ただ、個人的な雑感ですが、最後の積分の下端は0ではなく、エロージョンが開始する速度$V_c= \sqrt{\frac{2(R_t-Y_p)}{\rho_p}}$としたほうが見通しがいいかなと思います。

以上の式はAlekseevskiiとTateによって196年代後半にかけて展開されたモデルであり、非常に見通しの良い解を簡便に得ることが出来ます***。近いうちにこのモデルからわかるいくつかの簡単なことについてネタに出来たらなと思っています。

ただ、この式は侵徹の過程初期と侵徹後期では厳密なシミュレーションから外れることが知られているので、あくまでも初期のモデルの一つと理解するのがいいのかな？という感触です。侵徹深さに対する侵徹体長さの影響は明快ですが侵徹体径がどう影響するかは不明ですし、セルフシャープニングの影響を取り込むことも出来ません。わかりやすいだけに注意して取り扱うのがいいのかな、と考えています。

今回のネタを書くに当たり

Z. Rosenberg, E. Dekel, Terminal Ballistics, (2016), Springer

PJ Hazell: Armour: Materials, Theory, and Design, (2015), CRC Press

を底本にしました。

*この式に至る前にはもう一本式があるらしいですが細かいことはおいておきます。
**しかし$V_c$より僅かに速い程度では、材料強度の$R_t$項の寄与が大きく、侵徹はほとんど起こらないこともこの式から明らかです。
***密度が一緒だったら発散しないの？とかいうことは考えないようにしましょう。

ユゴニオ弾性限界の覚え書き

　お久しぶりです。冬コミは委託での参加でしたが、盛況だったようで何よりです。

　最近は鉄からちょっと離れてずっと気になってた、終末弾道というか高速度衝撃について調べていて、ユゴニオ弾性限界についての覚書を今回は。前回のセルフシャープニングの覚書の続きみたいなものですね。

本記事は、徹甲弾の侵徹の議論で出てくるユゴニオ弾性限界について定義を紹介するものであり、まじめに衝撃波について考えるものではありません。

ユゴニオ弾性限界

ユゴニオ弾性限界(Hugoniot Elastic Limit, HEL)というのはHEAT、あるいはAPFSDSの侵徹過程を特徴づけるものとして頻繁に話される用語だと思います。HEATやAPFSDSではユゴニオ弾性限界より高い圧力が生じるので、流体のように振る舞うことで侵徹が進行する、みたいな感じでしょうか？

Wikipediaを見ると「固体が塑性変形を開始し流体のように振舞う領域に入る境界線となる圧力である。」とありますが*、これは少し不思議な表現です。普通、弾が装甲に侵入する段階では、どちらかが変形するためには塑性変形する必要がありますから、侵徹の過程では塑性変形により流体のように振る舞っているはずです(低速度でも塑性変形に必要な応力を流動応力、Flow stressと呼ぶことは多々あります。)。

　そうすると、ユゴニオ弾性限界を特徴づける由来はもう少し別のところにある気がしてきます。

固体中の衝撃波

　ユゴニオ弾性限界の出自は固体の変形よりはむしろ衝撃波の研究にあります**。歴史的には、ユゴニオ弾性限界の出自は気体や液体について盛んに行われていた衝撃波の研究を固体についても同様に扱おうとする試みの中にあるかと思います。衝撃波の特徴として強調したいのは、以下の図に示すように圧力の変化が衝撃波の前面で不連続に起こるという点です***。

図1 衝撃波の模式図

　この点を元に、固体の変形挙動から大雑把に衝撃波が生じる条件について考えてみます。以下に、通常の応力ひずみ線図を示します。以下の式に示すように固体中の音速というのは応力ひずみ線図の各部分での傾きの平方根に比例します。
$$c \propto \sqrt{\frac{d\sigma}{d\varepsilon}}$$
　降伏応力以下ではフックの法則が成り立つため、傾きはひずみに依存せずヤング率で与えることが出来ますが、降伏応力より高い応力がかかったところでは図のように傾きが連続的に変化するため、それぞれのひずみ量がそれぞれの速度で動くことになります。そして、この図から普通に変形させた時には変形すればするほど変形した部分の移動速度は遅くなることがわかります。

図2 通常の応力ひずみ線図

というわけで、固体を普通に変形させたときには、小さなひずみ(低圧力)が先行し、その後から大きなひずみ(高圧力)が追従するため、衝撃波が発生しないことがわかります。

しかし、普通に変形させなければ話は別です。高速度衝撃で生じる様相を以下の図に示します**。

図3 高速衝撃で生じる一軸ひずみ領域

図にあるように、このような衝突では衝突目標と衝突体の界面付近に1D srain(1軸ひずみ領域)が生じますが、この一軸ひずみ領域の変形は、先程の応力ひずみ線図の場合の一軸引張/圧縮変形とは挙動が大きく異なります。圧縮することが出来る方向が一軸に決められているため、一軸ひずみ変形では軸方向の応力の他に、非常に大きな別の応力成分が生じることが雰囲気わかります。雰囲気雰囲気。なので、その応力ひずみ線図も一軸応力の応力ひずみ線図とは異なることが雰囲気察せられます。そこで、一軸ひずみの応力ひずみ線図(のようなもの)を以下の図に示します。

図4 一軸ひずみの応力ひずみ線図

この図から、大きくひずむほど傾きが大きく、つまり塑性変形の速度は塑性変形するほど大きくなることがわかります。これは期待通りの結果であり、このような変形過程では塑性変形について衝撃波(塑性変形の遷移)が生じることがわかります。この応力ひずみ線図は基本的には熱力学の要請によって決定され、この曲線はHugoniot curve、ユゴニオ曲線と呼ばれます(ユゴニオ曲線はユゴニオ曲線であって、厳密には応力ひずみ線図とは異なりますが…)****。

ユゴニオ弾性限界の見積もり

そうするとユゴニオ弾性限界は衝撃波が発生する一軸ひずみ変形が生じる場合の降伏応力と読み替えることが出来るでしょう。そうすると、HEATやAPFSDSの侵徹は図3のような一軸ひずみの存在が侵徹過程を通して影響を与える速度域で進行するために、ユゴニオ弾性限界が重要になってくるのかもしれません。

最後に、フックの法則を用いて一軸ひずみの降伏応力(HEL)と一軸応力の降伏応力(Yd)の関係を示して終わります。一般化フックの式で引張成分は以下で与えられますが、

$$\sigma_{11} = \frac{ \nu E }{( 1 + \nu )( 1- 2 \nu ) }(\varepsilon_{11}+\varepsilon_{22}+\varepsilon_{33})+ \frac{E}{1+ \nu } \varepsilon_{11}$$

$$\sigma_{22} = \frac{ \nu E }{(1+\nu) 1- 2\nu}(\varepsilon_{11}+\varepsilon_{22}+\varepsilon_{33})+ \frac{E}{1+\nu}\varepsilon_{22}$$

今回は一軸ひずみを考えているので、ε11以外はゼロとなり、書き下せば以下のようになります。(一軸応力であればσ11だけが非ゼロでその他の応力はゼロになります。vはポアソン比です。
$$\sigma_{11} = \frac{ \nu E }{( 1 + \nu )( 1- 2 \nu ) }(\varepsilon_{11})+ \frac{E}{1+ \nu } \varepsilon_{11}$$
$$\sigma_{22} = \frac{ \nu E }{(1+\nu) (1- 2\nu)}(\varepsilon_{11})$$

この二つの応力と、塑性変形が開始する条件を主応力の差が一軸引っ張りの降伏応力と一致する時(トレスカの条件)とすると、一軸ひずみ変形の降伏条件を以下のように得ることが出来ます。
$$\sigma_{11} - \sigma_{22} = Y$$

$$\sigma_{11}-\sigma_{22} = \frac{E}{1+\nu } \varepsilon_{11} = \frac{1-2\nu}{1-\nu}\sigma_{11}=Y$$

一軸ひずみ変形の降伏応力をHELとおくと、HELは一軸応力の降伏応力の(1-v)/(1-2v)倍となることがわかります。

$$\mathrm{HEL} = \frac{1-\nu}{1-2\nu}Y$$

ここから、一軸応力での降伏応力がわかればHELは自然とわかり、その逆もそうなります。

まぁユゴニオ弾性限界という語句自体が意味するところはこれくらいで、ここからダイレクトに侵徹にどうこう、と持っていくのはすこし間にギャップがあるような気がします。HEATやAPFSDSでよく言われる式として材料の強度を考慮に入れた修正ベルヌーイの式?がありますが、これはHELが条件というよりはむしろ、弾芯がエロージョンを起こしながら侵徹することが条件になってるような気がします。HELは超えているがエロージョンが起きないような状況は普通に起こり得て、以下の図にアルミ合金にタングステンを500 m/sで衝突させたときの応力の時間変化を計算した結果を示しますが**、

図4 WプレートをAlプレートに500m/sで衝突させたときの長軸方向の圧力変化

このようなエロージョンが問題にならない場合でも応力はHELを大きく上回ることがわかります。なので、HELだけがHEATやAPFSDSを特徴づけると言われるとなんだか不思議な気分になるなぁと思った話でした。

まぁ自分はこの辺全然詳しくないので間違ってたらスミマセン。

*ユゴニオ弾性限界 – Wikipedia

** PJ Hazell: Armour: Materials, Theory, and Design, (2015), CRC Press

***衝撃波 – Wikipedia

****:ユゴニオ曲線はP-V線図なのにどうなったら応力ひずみ線図になるんだ、という点は、横軸を(V0-V)/V0かln(V0/V)とすると一応応力ひずみ線図として示せるかな、と思ってます。