Ponceletの式のエネルギー保存則

よく知られているように、低速度侵徹ではエネルギー保存則が成り立ちます。(例えば、低速度侵徹のエネルギー効率)

Robinsの式ではその成立性は明らかに自明ですが、Ponceletの式では必ずしも自明ではなく、前回の記事では数値計算でお茶を濁していました。

最近遠出していた際に、時間があったので手で式展開したらきれいに導出できたので、ここに記事にしておきます。

侵徹に要するエネルギーは、結局装甲が消費するエネルギーなので、Ponceletの式での、侵徹体の加速度から、侵徹終了時までに消耗するエネルギー$E$は

$$a =C+BV^2$$

$$F  = ma  = m(C+BV^2)$$

$$E = \int_0^P Fdx$$

となります。ここで、$P$は侵徹深さで、$a$は抵抗側を正に取っています。

さて、$E$は$F$の$dx$に対する積分ですが、$V$は$a$に対して、

$$V = V_0 + \int_0^t a dt$$

なので、$t$と$x$の変数変換が必要です。なお、$V_0$は初速度です。

一見、私には非常に難しいのですが、実は低速度侵徹の侵徹深さの求め方 - Forrestalの式 -にその解き方は書かれていて、侵徹深さを求める定法から、

$$a = \frac{dV}{dt}=\frac{dx}{dt}\frac{dV}{dx}= V\frac{dV}{dx} = C+BV^2$$

です。これを適当に積分すれば、

$$x: 0\rightarrow  x$$

$$V: V_0 \rightarrow V$$

ですから、

$$-x = \int_{V_0}^V \frac{V}{C+BV^2}dV$$

から、

$$V^2= \frac{1}{B}\left(\left(C+BV_0^2\right)\exp{\left(-2Bx\right)}-C\right)$$

となります。この関係を$F$と$E$の関係に代入すれば、

$$E=\int_0^PFdx = m\int_0^Padx =m\int_0^P\left(C+B\times \frac{1}{B}\left(\left(C+BV_0^2\right)\exp{\left(-2Bx\right)}-C\right) \right)   dx$$

となり、$P= \frac{1}{2B}\ln{\left(1+\frac{BV^2}{C}\right)}$を代入して、具体的に積分すれば、装甲が消費するエネルギーは

$$E=m\frac{\left(BV_0^2+C\right)\exp{\left(-2BP\right)}}{2B} =\frac{mV_0^2}{2}$$

となります。すっきりですね。

酒飲みながら書いているので式展開は途中適当なところがあるかもしれませんが、結論は正しいはずです。

APFSDSとユゴニオ弾性限界と標的強度項$\sigma_s$

ユゴニオ弾性限界と標的強度項$\sigma_s$

低速度侵徹の侵徹深さの求め方 - Forrestalの式 -やAPFSDSとAPの差異とその由来では、APFSDSとAPのいずれの場合でもCavity expansion analysisから求められる$\sigma_s(=\sigma_{rr}(a))$が重要であることを見ました。
ところで、$\sigma_{rr}(a)$の導出の仮定で全範囲の関数を求めているので、$\sigma_{rr}(r)$を求めることが出来ます。降伏強度1GPaの鉄について計算した結果を示します。

 装甲は侵徹体の近傍で非常に大きな応力を生じて受け止めていることがわかりますね。図には、2本の縦線と、1本の横線が書かれていますが、縦線はそれぞれ、侵徹体径$a$、塑性変形領域径$c$に対応しており、横線はユゴニオ弾性限界に対応しています($\nu=1/3$)。

侵徹に寄与する応力は$r/a=1$の値なので、ユゴニオ弾性限界よりも遥かに大きな応力が生じていることがわかります。これは、ユゴニオ弾性限界を超えた後の応力ひずみ線図の覚え書きで示したように、周囲から拘束のある材料は塑性変形後も変形に必要な応力が急峻に増加するためです。

と、いうようなことを思うと、そろそろユゴニオ弾性限界はいいのではないかという気もしますが、どうでしょうか。

APFSDSとAPの差異とその由来

気が付いてみれば、2017年2月にユゴニオ弾性限界の覚え書きを投稿してから早4年、
ベルヌーイの式とエロージョンを伴う侵徹の覚え書き、ユゴニオ弾性限界を超えた後の応力ひずみ線図の覚え書き、衝撃インピーダンスってなんだ？な話、お手軽に高速度侵徹を解析したい話、ちょっと頑張ってL/D効果も含めた侵徹深さの評価をしたい話、APFSDSのRHA侵徹深さを簡単に見積もりたい話、AWモデルにおけるDynamic cavitation analysisの解析解、高速度侵徹のエネルギー効率、高速度侵徹の侵徹効率と高速度侵徹について多くの記事を書いてきました。(実は高速度侵徹はZukasのHigh velocity impact dynamicsに由来し、低速度侵徹はTerminal ballisticsの4章のEroding penetratorを高速度侵徹と訳した時に、どうRigid penetratorを訳すかと考えたときに生まれたオレオレ訳に過ぎないものです。)

最初のモチベーションは、ユゴニオ弾性限界に代表される特殊な変形挙動についての興味でしたが、多くの勉強をすることが出来ました。また、AWモデルという現代高速度侵徹モデルに触れたのちには、Cavity expansion analysisに感銘を受け、侵徹体の形状依存性を考慮してCavity expansion analysisを解くForrestalの式を追うことで、低速度侵徹の代表的なモデルにも触れることが出来ました。その過程は侵徹中に標的が作る抵抗 1 -静的Cavity expansion analysis-、侵徹中に標的が作る抵抗 2 -動的Cavity expansion analysis-、侵徹体の形状から侵徹体の重量を求める、低速度侵徹の侵徹深さの求め方 - Forrestalの式 -、Jacob de Marreの式と(低速度)徹甲弾の一般化した侵徹の解釈、低速度侵徹のエネルギー効率に示した通りです。

エネルギー保存の観点から見て、APFSDSのような高速度侵徹と、APのような低速度侵徹はまったく異なる振る舞いをします。低速度侵徹では常にエネルギーは保存されますが、高速度侵徹ではそうではありません。それは高速度侵徹のエネルギー効率で触れたように、「高速度侵徹では侵徹体が消耗するため」です。

さて、今回はこの「高速度侵徹では侵徹体が消耗する」ことと、「標的の抵抗はCavity expansion analysisで求められる」ことを出発点として高速度侵徹について考え、低速度侵徹と高速度侵徹の間の差異、ひいては、多くのミリタリー趣味者が心惹かれたに違いない、APFSDSのメカニズムがAPと決定的に異なる点について考えてみます。

結論

APFSDSのような高速度侵徹とAPのような低速度侵徹の間の本質的な差異とは、ユゴニオ弾性限界、衝撃インピーダンスや「流体としてふるまい（塑性流動）、相互侵食を起こして機械的強度を無視し、装甲を貫徹する。」といったことではなく、侵徹体が消耗するかどうかに起因します。

準備

とはいえ、抜き身では大変なので、高速度侵徹の特徴をここで一つ示してしまいましょう。

以下の図に、AWモデルで計算された、L=700mm, D=20mmの降伏強度0.5GPaのタングステン(WHA)製侵徹体を2000 m/sで降伏強度1GPaの鋼製(C>0.02)の標的に衝突させた際の侵徹体の先端速度と後端速度を示します。見て取れるように衝突直後から侵徹先端はある速度まで急速に減速し、その後は侵徹の大部分をその速度で行うことがわかります。

図1 侵徹体の速度プロファイル

ここで、侵徹体先端速度は以下の図の通り侵徹体が標的中を突き進む速度であり、侵徹深さを定める上でカギとなる値です。

 図2 単位時間あたりに標的を侵徹する距離

仮定

さて、図1に示した通り、高速度侵徹では侵徹体先端と侵徹体後端は異なる速度で運動しています。これはすなわち、侵徹体は塑性変形し、消耗しているということを示しているのですが、侵徹過程において、侵徹体は全体が変形しているわけではなく、標的と侵徹体の界面近傍のわずかな領域のみが塑性変形しています。ここでは今、この侵徹によって以下の図に示すように

1. 侵徹先端のわずかな領域が塑性変形し、(半球状になり)

2. その時、標的が作る$\sigma_z$はN=1/3のForrestalの式に等しい

と仮定します(N=1/3の妥当性はさておきましょう。)。図中$R_t$は$\sigma_s$のことです。

図3 侵徹体、標的の諸特性と侵徹体先端/標的界面の模式図

さて、今この侵徹体先端の塑性変形部と標的の間にどのような関係が成り立っているかを考えてみましょう。

衝突直後

衝突直後の侵徹体先端速度$u$の速度については議論の余地がありますが、ここでは$u_0$として、所与のものとしましょう。

Forrestalの式で見た通り、材料を動かすと、それに比例して$\frac{1}{2}\rho v^2$の圧力が生じます。

侵徹体と標的の界面(移動速度$u$)に着目すれば、侵徹体が$v-u$の速度で次々に界面に押し寄せているわけですから、

$$\sigma_{\rightarrow} = \frac{1}{2}\rho_p(v-u)^2$$

となります。また、侵徹体は自身が$Y$の強度を持っているとして、$\sigma_{\rightarrow}$は$\sigma_z$と同じように

$$\sigma_{\rightarrow} = Y+\frac{1}{2}\rho_p(v-u)^2$$
と置いておきます.

侵徹中期

さて、ここから$u$の加速度$\dot{u}$を求めてみましょう。$u$の加速度は、低速度侵徹でいう$\dot{V_z}$に相当します(低速度侵徹では$v=u$なので)。右向きを正とすれば

$$\begin{eqnarray} \dot{u} &=&\sigma_{\rightarrow} -\sigma_z &=&Y+\frac{1}{2}\rho_p(v-u)^2-\left(\frac{2Y}{3}\left(1+\ln{\frac{2E}{3Y}}\right)+\frac{1}{2}\rho_t u^2\right)\end{eqnarray}$$

となります。$\dot{u}$を直接解いてもいい(簡易的なAWモデル)のですが、ここでは、この式が時間発展した時の$\dot{u}$の振る舞いを考えてみましょう。

$\sigma_s, Y$が速度に依存しないなら、$\dot{u}$は$u, v$の関数となります。今、$u$が減速すると、$(v-u)^2$が大きくなり、$u$の減速は緩やかになります。一方で$u$が小さくなりすぎて$\dot{u}$が正になってしまって(なりませんが)加速するようになると逆に$\dot{u}$が負になるような変化となるので、結局十分長い時間ののちには、

$$\dot{u} = 0 =Y+\frac{1}{2}\rho_p(v-u)^2-\left(\frac{2Y}{3}\left(1+\ln{\frac{2E}{3Y}}\right)+\frac{1}{2}\rho_t u^2\right)$$

となります。$\dot{u}=0$は侵徹体先端速度$u$が変化しないということになりますが、これはまさに、図1の速度一定の領域のことを示しています。APFSDSのようなL/D比が大きな侵徹体では、まさにここで述べたような$\dot{u}=0$の領域が侵徹過程の大半を占めていることが知られています。

さて、$\dot{u}=0$となる侵徹体先端速度$u$は、上式が2次式であることから、解析的に解くことが出来て、

$$u = \frac{1}{1-\mu^2} \Biggr(v- \mu \sqrt{v^2 +\frac{2(\sigma_s-Y)(1-\mu^2)}{\rho_t}}\Biggl)$$ $$\mu = \sqrt{\frac{\rho_t}{\rho_p}}$$ $$A=\frac{2(\sigma_s-Y)(1-\mu^2)}{\rho_t}$$

を得ます。

これはベルヌーイの式とエロージョンを伴う侵徹の覚え書きで得られた形とまさに同じものです。

この式に、基づいて、侵徹を計算してみた結果を以下に示します(MBEが上式の結果を示しています)。

図4 各モデルで計算された侵徹体先端、後端速度

見てのとおり、両者のモデルは極めてよく一致しており、APFSDSのような高速度侵徹でも、侵徹の挙動にはCavity expansion analysisが極めて重要な役割を果たしていることがわかります。
ここで非常に注目されるべきことは、現代の高速度侵徹モデルであるAWモデルと、Forrestalの式に簡潔な仮定を加えて修正した高速度侵徹モデルが非常な一致を示すことです。

このことは何を意味しているのでしょうか？

これは、APFSDSのような高速度侵徹とAPのような低速度侵徹の間の本質的な差異とは、ユゴニオ弾性限界、衝撃インピーダンスや「流体としてふるまい（塑性流動）、相互侵食を起こして機械的強度を無視し、装甲を貫徹する。」といったことではなく、侵徹体が消耗するかどうかに起因することを意味しています。
侵徹体が消耗する挙動にこれらの要素が関わってきたり、斜め衝突での過渡的な挙動ではいろいろ難しい話があると思うのでいろいろあると思うんですが、まあ初手で出しても、よくわからないまま進むような気がしますね。
ちなみに、高速度侵徹のエネルギー効率、高速度侵徹の侵徹効率、低速度侵徹のエネルギー効率と同じ議論から、侵徹中に標的の強度が影響する割合$f_{R_t}$を、

$$f_{R_t} = P(V_0)\times R_t(V_0)/W$$

と定義して求めると、以下のような挙動になります。APFSDS程度の高速度侵徹では装甲の強度はまだまだ重要なことがわかります(そうでないとセラミックスを使わないという話もありますが)。面白いですね。

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

高速度侵徹のエネルギー効率については、高速度侵徹のエネルギー効率、高速度侵徹の侵徹効率で触れました。

これらの記事で興味深いのは、APFSDSの侵徹効率が最大となる速度は侵徹深さの衝突速度依存性とは全く異なることでした。せっかく低速度侵徹のモデルができましたので、低速度侵徹のエネルギー効率について考えてみましょう。

装甲が消費するエネルギー$E$は標的が作る$z$方向応力$\sigma_z$の侵徹深さ$P$についての積算なわけですが、

$$\sigma_z = \sigma_s + \frac{3}{2}N\rho_t V_z^2$$

$$W= \int_0^P\sigma_z dz$$

を解きましょう、ということになります。$\frac{3}{2}N\rho_t V_z^2$がなければ、簡単に

$$W= \int_0^P\sigma_z dz = \sigma_s P = \frac{2Y}{3}\left(1+\ln{\frac{2E}{3Y} } \right)$$

となます。初期(断面積当たり)運動エネルギー$KE0$を$\frac{1}{2}\rho_p(L+ka)V_0^2$として等式を立ててみれば、

$$\frac{1}{2}\rho_p(L+ka)V_0^2 =  \sigma_s P$$

$$P = \frac{\rho_p(L+ka)V_0^2}{2\sigma_s}$$

となって、これは低速度侵徹の侵徹深さの求め方 - Forrestalの式 -の最後の式に一致します。

では、動的な場合について求めてみましょう。となるわけですが、$\frac{3}{2}N\rho_t V_z^2$の$z$依存性を出すのは結構大変なので、数値的に計算した結果から$V-z$の関係を出し、積分することにします。高速度侵徹の場合と同じように、

$$\eta = \frac{W}{KE0}$$

$$\eta_p = \frac{\sigma_sP}{KE0}$$

と定義してプロットをしてみると

となって、低速度侵徹ではエネルギーは保存されていることがわかります。一方、侵徹効率は衝突速度が高くなるにしたがって低下している、すなわち、$\frac{3}{2}\rho_tNV_z^2$の寄与が大きくなることがわかります。

とはいえ、高速度侵徹に比べればその効率はずっといいですね。

Jacob de Marreの式と(低速度)徹甲弾の一般化した侵徹の解釈

 GWがまさかブログ連投ウィークになるとは誰が予想したでしょうか。

いつかまとめなきゃなぁと思っていたんですが、sympyを使うと割と簡単にlatexの形に出力できることに気が付いたために今回の量産が出来ています。ありがとう、sympyとMathJax。

Cavity expansion analysisを使った解析は現代の侵徹理論では滅茶苦茶重要なので、名前だけでも憶えてもらえると嬉しいです。

今回の結論

Jacob de Marreの式では侵徹深さが衝突速度$V$に対して$V^{1.4}$とかになるのはなぜ？→そもそもそういう形をしていないから

本文

さて、このブログを読んでいらっしゃる方はJacob de Marreの式はご存じかと思います。

Jacob de Marreの式は侵徹深さの実験的な評価式で、材質によっていろいろな定式があるかと思いますが、滅茶苦茶おおざっぱに言えば侵徹深さ$P$と、衝突速度$V_0$の間に、

$$P \propto V_0^n$$

の関係があるという式で、$n$はだいたい1.4~1.5の値を取るという理解です。侵徹体寸法や質量は今回は見ないことにします。

$n$が1.4~1.5に比例するというのは奇妙な話です。例えば、我々がプリミティブに侵徹深さを求めようとして、侵徹体の運動エネルギーと、装甲が消費するエネルギーの間に成り立つ関係を立ててみれば、

$$\frac{1}{2}mV_0^2 = YSP$$

$$P \propto V^2$$

となって、指数は2になるからです。

この差をどうとらえるかは人によって異なるかと思いますが、個人的には大きな差だと考えています。なぜなら、衝突速度が500m/sと1000m/sの場合を比較して考えようとすると、$n$が2と1.4では、1000m/sの時の侵徹深さは500m/sの時の侵徹深さの4倍、2.6倍となり、大きな乖離が生じるためです。

$n$が2でないことの理由として、

1. 侵徹はエネルギー保存則を守っていない
2. 装甲の強度以外の別の抵抗が存在する

が考えられますが、1はありえないので、2.について考えるのが妥当でしょう。

今、我々はすでに、標的が作る抵抗と標的の強度とある時刻の侵徹速度について、

$$\sigma_{z} = \frac{3NV_{z}^{2} \rho_{t} }{2} + \frac{2Y}{3}  \left(\log{\left(\frac{2 E}{3 Y} \right)} + 1\right)$$

という関係があることがわかっています。そして、この抵抗と侵徹深さ$P$の間には、

$$P = \rho_p \frac{L+ka}{3N\rho_t}\ln{\left(1+ \frac{3N\rho_t}{2\sigma_s}V_0^2 \right)}$$

というForrestalの式という関係があることがすでに分かっています。

今一度、Jacob de Marreの式とForrestalの式を比較してみましょう。

$$P \propto V_0^n$$

$$P = \rho_p \frac{L+ka}{3N\rho_t}\ln{\left(1+ \frac{3N\rho_t}{2\sigma_s}V_0^2 \right)}$$

この二つの式が似ていると感じる人は少ないのではないかと思いますが、実際この式の二つをプロットすると、その形は大きく異なっています。以下に、Forrestalの式、Jacob de Marre、エネルギー保存則の結果を示します。ここで、Jacob de Marreの式、エネルギー保存則の式はForrestalの式をフィッティングすることでプロットしています。

速度が低い領域では、Jacob de Marreの式、エネルギー保存則の式はいずれもForrestalの式に追従していますが、速度が上昇するにつれてForrestalの式の傾きは小さくなり、両者から乖離していることが見て取れます。しかし、Jacob de Marreの式はその指数が小さいためにより高速度域まで追従していることがわかります。

Forrestalの式がこのように高速度域で侵徹深さが増加しなくなる理由は、

$$\sigma_{z} = \frac{3NV_{z}^{2} \rho_{t} }{2} + \frac{2Y}{3}  \left(\log{\left(\frac{2 E}{3 Y} \right)} + 1\right)$$

の第1項に起因しています。衝突速度が低いとき、第1項は無視できるほど小さく、その侵徹深さは

$$P = \rho_p \frac{L+ka}{2\sigma_s}V_z^2$$

に漸近します。一方、衝突速度が極めて高いとき、

$$\ln{\left(1+ \frac{3N\rho_t}{2\sigma_s}V_0^2 \right)} \approx \ln{\left( \frac{3N\rho_t}{2\sigma_s}V_0^2 \right)}$$

となって、侵徹深さは衝突速度$V_0$の$\ln{V_0^2}$に比例します。そこで、Forrestalの式の高速度域について$\ln{V_0^2}$にフィッティングした結果を以下に示します。

このことから、低速度侵徹では標的強度が支配的な領域と、衝突速度の2乗に比例する慣性項ともいえる項が支配的な領域の二つに分かれており、高速度になるとともに慣性項の影響が大きくなり、侵徹深さの衝突速度依存性は鈍くなることがわかります。

よって、Jacob de Marreの式のように、$V$の指数は2よりも小さくなることがわかります。

ということで、Jacob de Marreの式の速度依存性が$V^2$でない理由は、そもそも低速度侵徹の速度依存性は$V^n$の形をしていないから、でした。

また、CRH、標的強度$\sigma_s$、標的密度$\rho_t$が変わると、$n$はどんどんと変わります。まあJacob de Marreの式の指数はざっくりそのくらい、と考えるのがいいのではないでしょうか。

そうするとJacob de Marreの式でのエネルギー保存則ってどうなっているの？と思うのが人情ですが、それは機会があればまた。(ちゃんと保存されています。)

Forrestalの式の有効性の確認

 前回の記事で低速度侵徹の解析モデルを導出しました。

まあそれって実際どれくらい会うの？というのが気になるところだと思います。

試験結果はここにあるらしいので、これを使ってみましょう。gsysで手拾いしたプロットを、Forrestalの式の結果をプロットすると、ざっくり以下の通りになります。(プロットの詳細は議論の余地がありますが（そういう論文

合うといっていいのではないでしょうか。

低速度侵徹の侵徹深さの求め方 - Forrestalの式 -

 前記事までで、$\sigma_{rr}$の具体的な形を求めることが出来ましたので、最後に侵徹体軸方向に生じる力を求めましょう。

計算ノートはgithubに置いています。

以下の図に今回用いる関係を示します。

まず、応力$\sigma$と力$F$の間には面積$S$を用いて$\sigma =FS$の関係があるので、角度$\theta$での面積を求める必要があります。

とはいえ、似たようなことは$k$の計算でもやっています。

角度$\theta$における球の表面積$dS$は$dS = 2\pi rdL$ですが、$r,dL$はそれぞれ

$$r = s\sin{\theta}-\left(s-a\right)$$

$$dL = sd\theta$$

ですから、 

$$dS = 2\pi s\left(s\sin{\theta}-\left(s-a\right)\right) d\theta$$  

となります。よって、角度$\theta$での$dFz$は

$$dF_z = \sigma_{rr}^z dS = 2\pi \sigma_{rr}s^2\left(\sin{\theta}\cos{\theta} -\left(\frac{s-a}{s}\right)\cos{\theta} \right) d\theta$$

となります。ここで、$\sigma_{rr}$の$z$方向成分$\sigma_{rr}^z$は$\sigma_{rr}\cos{\theta}$であることを用いています。

また、$\sigma_{rr}$は$V_z,\theta$の関数$\sigma_{rr}\left(V_z, \theta\right)$です。

あとは上式を$\theta: \theta_0 \rightarrow \pi/2$まで積分して終わりです。ここで、$\theta_0$は$\sin^{-1}\left(\frac{2\psi-1}{2\psi}\right)$です。

ここでは、$\sigma_{rr}(a)$を

$$\sigma_{rr}(a) =\frac{2Y}{3}\left[ 1+\ln{\frac{2E}{3Y}}\right] + \rho_t\left(\frac{3(V_z\cos{\theta})^2}{2}\right)$$

として、動的Cavity expansion analysisの解について解いてみましょう。

定義に代入して、適当に積分すれば

$$F_{z} = \pi a^{2} \left(\frac{V_{z}^{2} \rho_{t} \left(8 \psi - 1\right)}{16 \psi^{2}} + \frac{2 Y \left(\log{\left(\frac{2 E}{3 Y} \right)} + 1\right)}{3}\right)$$

として、$F_z$を得ます。ここで、

$$\frac{1}{16} = \frac{3}{2}\times \frac{1}{24}$$

$$N = \frac{8\psi-1}{24\psi^2}$$

として、上式を置換することで、

$$F_z = \pi a^2\left[ \frac{2Y}{3}\left(1+\ln{\left(\frac{2E}{3Y}\right)}  \right)+ \frac{3}{2}N\rho_t V_z^2 \right]$$

として簡単な形で$F_z$が得られます。上式は動的Cavity expansion analysisについての結果ですが、静的Cavity expansion analysisの結果は第一項のみを持ってきて、

$$F_z = \pi a^2\left[ \frac{2Y}{3}\left(1+\ln{\left(\frac{2E}{3Y}\right)}  \right) \right]$$

として得られます。

これで、侵徹体に掛かる力$F_z$と質量$m$を導くことができましたので、侵徹体に関する運動方程式を

$$m\dot{V}_z = -F_z$$

としておくことが出来るようになりました。ここで、$\dot{V}_z$は$V_z$の時間微分で、侵徹体の加速度です。

今、$m$は

$$m= \pi \rho_p a^2\left(L+ka\right)$$

としてわかっているので、

$$-\rho_p \left(L+ka\right)\dot{V_z} =  \frac{2Y}{3}\left(1+\ln{\left(\frac{2E}{3Y}\right)}  \right)+ \frac{3}{2}N\rho_t V_z^2 = \sigma_s + \frac{3}{2}N\rho_t V_z^2$$

と具体的な形がわかります。

このような、

$$a \dot{V}_z = b +c V_z^2$$

の積分は、

$$a\frac{dV_z}{dt} = a \frac{dz}{dt}\frac{dV_z}{dz} = aV_z\frac{dV_z}{dz} = b+cV_z^2$$

$$dz = \frac{aV_z}{b+cV_z^2} dV_z$$

$$z: 0 \rightarrow P$$

$$V_z : V_0 \rightarrow 0$$

となり、また、$f=b+cV_z^2$とおけば、分子の$aV_z$は$\frac{a}{2c}f^{\prime}$であるので、

$$\int_{V_0}^0 \frac{f^{\prime}}{f} dV= \ln{f(V_0)/f(0)}$$

から、

$$P = \frac{a}{2c}\ln{\left(1+\frac{cV_0^2}{b}\right)}$$

として、求められます。$a,b,c$に代入すれば、

$$P = \rho_p \frac{L+ka}{3N\rho_t}\ln{\left(1+ \frac{3N\rho_t}{2\sigma_s}V_0^2 \right)}$$

となり、これが動的Cavity expansion analysisから求められる侵徹深さであり、Forrestalの式と呼ばれるものです。

また、静的Cavity expansion analysisの結果からは、

$$P = \rho_p \frac{L+ka}{2\sigma_s}V_z^2$$

が得られます。

以上が低速度で最も基本的な解析モデルです。

この辺の計算は今回全部sympyで計算してやったんですが、積分もきれいにしてくれて滅茶苦茶便利ですね。

去年くらいに必死に手で導出した苦労は一体...(大切ですが