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で計算してやったんですが、積分もきれいにしてくれて滅茶苦茶便利ですね。

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

侵徹中に標的が作る抵抗 2 -動的Cavity expansion analysis-

 前回力の釣り合いの式は以下に示すように、右辺が0となっていました。

$$\frac{d\sigma_{rr}}{dr}+\frac{2}{r}\left(\sigma_{rr}-\sigma_{\theta\theta}\right) = 0$$

上式は力の釣り合いの式なわけですが、標的の各部分についての運動方程式とも言えます。

力が釣り合っていて、右辺が0ということは、侵徹中にも関わらず標的の各部分は静止をしていることを意味しており、

その意味で、前回の結果は静的cavity expansion analysisとも呼ばれます。

では、力が釣り合わず、標的の各部分が運動する状況を考えると、その$\sigma_{rr}$はどうなるでしょうか？というのが今回の話題です。

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

まず、力の釣り合いの式を、見慣れた運動方程式の形にすると、

$$\frac{d\sigma_{rr}}{dr}+\frac{2}{r}\left(\sigma_{rr}-\sigma_{\theta\theta}\right)  =- \rho_t \left( \frac{\partial v}{\partial t} +v\frac{\partial v}{\partial r} \right)$$

となります。ここで、右辺が妙な形になっているのは、オイラー座標系での時間微分になっているためです。$v$は各部分の粒子速度であり、$r, u$と同じ向きで外側が正です。また、$v$は$u$の時間微分ですので、同じような手続きを取ると、

$$v = \left( \frac{\partial u}{\partial t} +v\frac{\partial u}{\partial r} \right)$$

$$\frac{\partial u}{\partial t} = v \left(1 - \frac{\partial u}{\partial r} \right)$$

として$u$と$v$の関係が得られます。

今、前回の記事から$u$の形は以下の通り

$$u=r\left[1-\left(1-\frac{a^3}{r^3}\right)^{1/3}\right]$$

として与えられているので、$a$を

$$a = a(t)$$

なる関数だと思って微分して

$$\frac{\partial u}{\partial t} = v \left(1 - \frac{\partial u}{\partial r} \right)$$

に代入すれば、

$$v = \frac{\dot{a} a^{2}{\left(t \right)}}{r^{2}}$$

を得ます。ここで、$\dot{a}$は$a$の時間微分です。

この結果を用いて、前回と同じくまずは降伏領域($a\geq r \geq c$)の$\sigma_{rr}$を求めましょう。

本項第2式に$\sigma_{rr}-\sigma_{\theta\theta} = Y$と$v$を代入して、

$$\frac{d}{d r} \sigma_{rr}{\left(r \right)} = - \frac{2 Y}{r} - \rho_{t} \left(\frac{\ddot{a} a^{2}{\left(t \right)}}{r^{2}} + \frac{2 \dot{a}^{2} a{\left(t \right)}}{r^{2}} - \frac{2 \dot{a}^{2} a^{4}{\left(t \right)}}{r^{5}}\right)$$

という形で得られます。ここで$\ddot{a}$は$a$の時間に対する2階微分です。

前回の記事と同じように積分してやれば、

$$\sigma_{rr}(a) = \sigma_{rr}(c)-2Y\ln{\frac{a}{c}} +\rho_t\left[ -(\ddot{a}a^2+2a\dot{a}^2)\left(\frac{1}{c} - \frac{1}{a}\right) + \frac{a^4\dot{a}^2}{2}\left(\frac{1}{c^4}-\frac{1}{a^4}\right) \right]$$

として$\sigma_{rr}(c)$を得ます。それでは、$\sigma_{rr}(c)$を求めるために、前記事と同じく弾性変形領域($c\leq r$)を解きましょう。

やることは大して変わらなくて、結局、

$$\frac{d\sigma_{rr}}{dr}+\frac{2}{r}\left(\sigma_{rr}-\sigma_{\theta\theta}\right)  =- \rho_t \left( \frac{\partial v}{\partial t} +v\frac{\partial v}{\partial r} \right)$$

に

$$\sigma_{rr}-\sigma_{\theta\theta} = \frac{4Ea^3}{3r^4}$$

を代入して、

$$\sigma_{rr}(a) = \frac{4Ea^3}{9c^3} -2Y\ln{a/c}+ \rho_t\left(\ddot{a}a+\frac{3\dot{a}^2}{2}\right)$$

となります。前回示したように、

$$\left(\frac{a}{c}\right)^3 =  \frac{3Y}{2E}$$

ですから、

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

となります。また、$\ddot{a}$は極めて小さいこと、$\dot{a}$、すなわちCavityの拡大速度は侵徹速度$V_z$のCavityの$r$方向成分なので$\dot{a}=V_z\cos{\theta}$となることを用いて

$$\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での$\sigma_{rr}$を得ます。

侵徹中に標的が作る抵抗 1 -静的Cavity expansion analysis-

 侵徹過程における侵徹体の運動は、侵徹体が標的から受ける抵抗によって決まることに疑いの余地はありません。

とはいえ、その抵抗をどうやって求めるか、はなかなか難しい問題です。

本記事と次回、次々回の記事で、標的が作る抵抗を求め、侵徹深さを直接求めてみようと思います。

まず、今考えている状況を図に示します。下の図は前回示した侵徹体の模式図に、侵徹体に加わる応力を加えています。

さて、侵徹体は標的を押し広げるように変形させながら侵徹をするので、侵徹体周辺で標的は以下の図に示すように、侵徹体の周囲に塑性変形領域があり、その外側には弾性変形領域があります。ここで$a$は侵徹体半径で、$c$は塑性変形領域の半径ですがまだわかっていません。

今、この$\sigma_{rr}$を求めてみましょうというのが今回の話題です。

図に示したように侵徹体先端付近の標的は塑性/弾性変形しています。$a \leq r \leq c$では標的は塑性変形しており、$c\gt r$では弾性変形のみが生じています。侵徹体の運動を考えるためには、変形する標的の侵徹体/標的界面で作る侵徹体を押し返す力($F_z$)を求める必要があります。

本項と次回では、単位面積あたりに動径方向に働く力、$\sigma_{rr}$を求めてみます。次々回では、$\sigma_{rr}$から、$z$方向に働く力$F_z$を求めます。

さて、本項では、球対称な極座標系(f(r))での標的の変形を考える必要がありますが、出発点となるのは以下の2式です。

$$\frac{d}{dr}\left[\left(r-u\right)^3\right] = 3r^2$$

$$\frac{d\sigma_{rr}}{dr}+\frac{2}{r}\left(\sigma_{rr}-\sigma_{\theta\theta}\right) = 0$$

ここで、第1式は標的の各部分の変形前後での質量(体積)保存則、第2式は標的の各部分に成り立つ力の釣り合いの式です。

第1式では、圧力に対して標的の体積は不変であるという、標的を非圧縮性の材料と仮定しています。

また、$r,u,\sigma_{rr},\sigma_{\theta\theta}$はそれぞれ、侵徹体中心からの距離、$r$における変位、動径方向の応力、円周方向の応力です。

また、動径方向のひずみ、円周のひずみをそれぞれ、

$$\varepsilon_{rr} = \ln{\left(1-\frac{du}{dr}\right)} \approx -\frac{du}{dr}$$

$$\varepsilon_{\theta\theta} = \ln{\left(1-\frac{u}{r}\right)} \approx -\frac{u}{r}$$

と定義します。

$\sigma_{rr}$と$\sigma_{\theta\theta}$、標的の降伏強度の間にはトレスカの降伏条件から、

$$\sigma_{rr}-\sigma_{\theta\theta} = Y$$

の関係があります。(加工硬化を仮定した形も解くことが出来るが、簡単のため、完全弾塑性体を仮定します。)

さて、侵徹体/標的界面($r=a$)では、標的は$a$まで押しのけられているので$u(r=a)=a$です。

また、$r<a$には標的はないので考える必要はありません。

この初期条件を用いて、第1式を$r:a \rightarrow r (r\lt c)$まで解くと、

$$(r-u)^3 = r^3-a^3$$

となるから、$u$は$r$について

$$u=r\left[1-\left(1-\frac{a^3}{r^3}\right)^{1/3}\right]$$

となります。(今回は使わないのですが、一応導出しておきます。完全弾塑性でない仮定を置くと使います。)

次に、$\sigma_{rr}$を求めたいわけですが、最初に述べたように、標的は塑性変形領域と弾性変形領域に分かれています。塑性変形領域については、トレスカの条件から第2式中の$\sigma_{rr}-\sigma_{\theta\theta}$は$Y$であることがわかっているので、直ちに解くことが出来ます。

第2式に$\sigma_{rr}-\sigma_{\theta\theta}=Y$を代入して$r:a\rightarrow c$まで積分すれば、

$$\frac{d\sigma_{rr}}{dr} = \frac 2 r Y$$

$$\sigma_{rr}(a) = \sigma_{rr}(c) -2Y\ln{\frac{a}{c}}$$

を得ます。

$\sigma_{rr}(c)$については、弾性変形領域について考えることで明らかになります。

ということで、$r\gt c$の弾性変形領域について考えてみます。

弾性変形領域ではひずみが小さい(<0.2%以下)なので、常に$r \gg u$です。ということで、第6式中の$u$の2次以上の項を無視すれば、

 

  $$r^3-3r^2u+3ru^2-u^3 = r^3-a^3$$

  $$-3r^2u = -a^3$$

  $$u = \frac{a^3}{3r^2}$$

が得られます。ここで得られた$u$を第3,4式のひずみの定義(第3項)に代入すれば、

$$\varepsilon_{rr} = -\frac{du}{dr} = \frac{2a^3}{3r^3}$$

$$\varepsilon_{\theta\theta} = -\frac u r = - \frac{a^3}{3r^3}$$

となります。さて、$\varepsilon_{rr},\varepsilon_{\theta\theta}$間のモールの応力円を考えると、最大せん断ひずみは$\varepsilon_{rr}-\varepsilon_{\theta\theta}$ですから、

$$\sigma_{rr} - \sigma_{\theta\theta} = 2k = 2G\left(\varepsilon_{rr}-\varepsilon_{\theta\theta}\right)$$

となり、上式に$\varepsilon_{rr},\varepsilon_{\theta\theta}$を代入すると、

$$\sigma_{rr}-\sigma_{\theta\theta} = \frac{4Ga^3}{9r^4} = \frac{4Ea^3}{3r^4}$$

となります。ここで、今、標的は非圧縮性であることを仮定しているのでポアソン比$\nu=0.5$であること、$G=\frac{E}{2(1+\nu)}$であることを用いています。

これで、$\sigma_{rr}-\sigma_{\theta\theta}$を第2式を代入すると、

$$\frac{d\sigma_{rr}}{dr} = -\frac{4Ea^3}{3r^4}$$

となり、これを$r:c\rightarrow \infty$について積分すると、

$$\sigma_{rr}(c)=\frac{4Ea^3}{9c^3}$$

を得ます。

今、明らかに$c$において$\sigma_{rr} - \sigma_{\theta \theta}= Y$であることから、

$$\sigma_{rr} - \sigma_{\theta \theta} =\frac{3}{2} \sigma_{rr} = Y$$

$$\sigma_{rr}(c) = \frac{2}{3} Y = \frac{4Ea^3}{9c^3}$$

$$\left(\frac{a}{c}\right)^3 =  \frac{3Y}{2E}$$

となります。これらを、第9式に代入すると

$$\sigma_{rr}(a) = \frac{4Ea^3}{9c^3}+2Y\ln{\frac{a}{c}} = \frac{2Y}{3}+\frac{2Y}{3}\ln{\frac{3Y}{2E}}$$

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

が得られます。

これが、準静的な仮定から得られる動径方向の応力です。

この一連の導出をcavity expansion analysisと呼びます。

この結果と、侵徹深さの間の関係は、また次回(次々回)

侵徹体の形状から侵徹体の重量を求める

 本記事では典型的な侵徹体の質量を、その形状因子のみから計算することを目的とします。質量は体積が明らかであれば密度を掛ければ算出できるので、ここでは特に体積の導出方法について検討します。(逆に形状と質量があれば平均的な密度を求めることが出来ます。

今回はsympyで展開したので、興味がある方はこちらをご覧ください。

侵徹体の形状は様々なものがありますが、ここでは以下のように、半径$a$、長さ$L$の円柱部に半径$s$の円を侵徹体軸上で回転させた先端部を貼り付けた侵徹体を考えましょう。この侵徹体形状は低速度侵徹を解析的に解くことが可能な形状であり、侵徹の試験において多く用いられています。

このような侵徹体の先端形状を定義する際に、CRH(Calibre-Head radius)と呼ばれる以下の指標を用います。

$$\psi = \frac{s}{2a}$$

この指標を用いることで、侵徹を簡潔に表すことができます。(その辺の詳細はまた後日

本項の目標は、上述の特徴を有する侵徹体の体積$V$を

$$V = \pi a^2 (L+ka)$$

のように、円柱部の断面積を用いて簡潔に表現することです。この$k$を求めてみましょう。

今、上の図には、$\theta_0, l$が定義されていますが、その具体的な形はわかっていません。まず、図から明らかに、

$$\sin{\theta_0} = \frac{s-a}{s} = \frac{2\psi-1}{2\psi}$$

なので、$\theta_0$は$\sin^{-1}{\frac{2\psi-1}{2\psi}}$であることがわかります。また、$l$は明らかに$s\cos{\theta_0}$ですが、

$$\cos{\theta_0} = \sqrt{1-\sin^2{\theta_0}}=\sqrt{\frac{4\psi^2-(4\psi^2-4\psi+1)}{4\psi^2}}=\frac{\sqrt{4\psi-1}}{2\psi}$$

であること、また、$s=2\psi a$であることから、

$$l=a\sqrt{4\psi-1}$$

であることがわかります。

では、侵徹体先端部の重さ(=体積)を求めましょう。

侵徹体軸上に、$\theta_0 \rightarrow \pi/2$の範囲で存在する半径$s$の円の体積を求めるためには、

$$y = s \sin{(\theta)} - (s-a)$$

$$x = s \cos{(\theta)}$$

$$V_o = \int_{0}^{l} \pi y^2 dx$$

を解けばいいわけですが、$x$は$\theta$の関数なので、$\theta$についての導関数を求めて置換すると、

$$dx=-s\sin(\theta) d\theta$$

$$l:\ 0 \rightarrow l$$

$$\theta:\ \pi/2 \rightarrow \theta_0$$

から、

$$V_o = \int_{\pi/2}^{\theta_0} -\pi s\sin(\theta)  y^2 d\theta$$

となり、これを積分して気合いで整理すると、

$$V_o = \pi a^{3} \left(- 8 \psi^{3} \operatorname{asin}{\left(\frac{\sqrt{4 \psi - 1}}{2 \psi} \right)} + 4 \psi^{2} \operatorname{asin}{\left(\frac{\sqrt{4 \psi - 1}}{2 \psi} \right)} + \sqrt{4 \psi - 1} \left(- 4 \psi^{2} + \frac{4 \psi}{3} - \frac{1}{3}\right)\right)$$

となります。ここまでくれば、

$$V = \pi a^2 (L+ka)$$

と式の形を合わせるように、$k$を

$$k = \left(- 4 \psi^{2} \left(2 \psi - 1\right) \operatorname{asin}{\left(\frac{\sqrt{4 \psi - 1}}{2 \psi} \right)} + \sqrt{4 \psi - 1} \left(- 4 \psi^{2} + \frac{4 \psi}{3} - \frac{1}{3}\right)\right)$$

と定めれば、当初欲しかった$k$が得られます。

後は$V$に侵徹体の材料の密度$\rho$を掛ければ、侵徹体の重量がわかりますし、逆に重量がわかっていれば、侵徹体の平均的な密度がわかります。