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年3月3日日曜日
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)
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侵徹深さを簡単に見積もりたい話
お久しぶりです。
これを線形近似して、傾きを求めると、冒頭の
$$DoP[\mathrm{mm}] = 43 KED[\mathrm{kJ/mm^2}] +10$$
が得られます。(切片を0にしたほうが物理的な意味はありますが
この間コード公開した記事アップロードばっかりだから～と高をくくっていたら前回の記事は去年の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$は
というのも、運動エネルギーは$\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の低下から生じる問題を解決しようと思えば、それなりの運動エネルギーが必要になります。そういう意味で大口径は七難隠す、と言えるかもしれませんね（１月頭くらいにツイッターでそんな話をしてました)。
その辺の問題を与えられた口径の範囲で調整した結果が今ある多くの侵徹体の形、と考えるのが自然な形かもしれません。
そういうある種の最適化問題に興味はあるんですが、侵徹体径の効果を考える問題はいろいろ難しくて手が回りません。
しかし、例えば装甲材料の強度を変化させたとき、侵徹深さにどのような影響を与えるか、といったことは妥当性を持って言えるでしょう。
というわけで適当に標的強度を変化させたときの侵徹深さの依存性を求めると、以下のようになります。強い装甲は強い、ということですね。
とはいえ、適当な侵徹体が与えられたときに、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
もちろん、砲に要求される能力は上の式に断面積を掛けたものになるので、高L/D比に伴うDの低下から生じる問題を解決しようと思えば、それなりの運動エネルギーが必要になります。そういう意味で大口径は七難隠す、と言えるかもしれませんね（１月頭くらいにツイッターでそんな話をしてました)。
その辺の問題を与えられた口径の範囲で調整した結果が今ある多くの侵徹体の形、と考えるのが自然な形かもしれません。
そういうある種の最適化問題に興味はあるんですが、侵徹体径の効果を考える問題はいろいろ難しくて手が回りません。
逆に何が言えるか
今度はここから逆に何が言えるかを考えてみると、実のところほとんど言えません。それは、複合装甲ではあれこれがどうこうで…とかなんだか難しそうな話ではなく、シンプルにするために様々な単純化を行っているためです。しかし、例えば装甲材料の強度を変化させたとき、侵徹深さにどのような影響を与えるか、といったことは妥当性を持って言えるでしょう。
というわけで適当に標的強度を変化させたときの侵徹深さの依存性を求めると、以下のようになります。強い装甲は強い、ということですね。
とはいえ、適当な侵徹体が与えられたときに、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$を計算するルーチンで参考にした論文にミスがあったことがわかったので修正しています。
ダウンロードはこちらからどうぞ。
環境としては、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モデルは以下の式で表されるように侵徹体と標的界面での圧力の釣り合いを考え、
\[ 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モデルのこういう困難から、個人的にこれ以上弄ることの難しさを感じていました。
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
あたりを参考にしていますが、実際の侵徹体についてはあんまり詳しくないので明らかにおかしい点があれば指摘いただければ幸いです。
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に強く影響します*。
例えば、タングステン侵徹体と高強度鋼標的の侵徹深さを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とし、それ以外は同じにしています。
また、これらのことから、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
あたりを参考にしていますが、実際の侵徹体についてはあんまり詳しくないので明らかにおかしい点があれば指摘いただければ幸いです。
2017年10月29日日曜日
お手軽に高速度侵徹を解析したい話
夏コミも終わって小休止と言った感じであんまり本を読めていません。
鹿部です。
今回はベルヌーイの式とエロージョンを伴う侵徹の覚え書きで得た結果から実際に計算するプログラムを作る話です。夏コミの本で行った解析は今回のコードを基本にしたコードで行われていて、今回のコードを使えば高L/D比な弾芯による高速度侵徹のモデルを大まかに扱えます。
基礎的なモデルは当該記事を参照のこと。
ちなみに、Solution of the Long Rod Penetration EquationsにはFortran版のソースコードが載っていますので、今回の話はあんまり目新しいものはありません。
自分がPythonしか使えないから最後に出てくる数値積分を簡単に扱うためです。
今の侵徹体後端速度$V$から侵徹速度$U$を求める関数は簡単に
次に侵徹体後端速度と侵徹体残存長さ$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のように得られます。
ここまでくればあと少しで、これを
\[P = \frac{\rho_p}{Y_p} \int_0^{V_0} u\ l\ dV \]
に突っ込めば侵徹長さが出ます。後は$U$と$l$の積をscipyのquad関数に投げればよくて、
侵徹体としてタングステン、標的として高強度鋼を想定して、$\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の間で変化させながら計算してみると、以下のような図を得ます。
これはタングステン弾芯と高強度鋼標的の衝突ではありがちな曲線です。
鹿部です。
今回はベルヌーイの式とエロージョンを伴う侵徹の覚え書きで得た結果から実際に計算するプログラムを作る話です。夏コミの本で行った解析は今回のコードを基本にしたコードで行われていて、今回のコードを使えば高L/D比な弾芯による高速度侵徹のモデルを大まかに扱えます。
基礎的なモデルは当該記事を参照のこと。
ちなみに、Solution of the Long Rod Penetration EquationsにはFortran版のソースコードが載っていますので、今回の話はあんまり目新しいものはありません。
何で書くか
今回は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}}$です。
\[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などにある実験結果を用いてフィッティングして、全体の曲線を見積もるなどすると楽しめる気がします。
1. 強度はあくまでも計算パラメーター
タングステン弾芯と高強度鋼標的について計算した際、Rtを6GPaとして与えました。これはJ. ZukasのHigh velocity impact dynamicsよりパラメーターを拝借したものですが、高強度鋼は普通このような強度をユゴニオ弾性限界であっても示しません。なので、特に標的強度については材料的な意味を求めるのは中々難しいものがあります。
2. Primary phaseのみを考慮している
高速度侵徹は大きく分けて3つに侵徹過程を分けることができますが、今回のコードではその中で最も主要な役割を果たすPrimary phaseのみを考慮しています。極端な低速度/高速度はあまり見ないほうが良いような気がします。
3. L/D比を考慮することができない
高速度侵徹はL/D比によって侵徹深さが変わるのはよく知られた事実ですが、この現象をうまく整理するのはなかなか難しいものがあります(試みられてはいます)。
こういう事情がありますが、dticなどにある実験結果を用いてフィッティングして、全体の曲線を見積もるなどすると楽しめる気がします。
2017年8月28日月曜日
衝撃インピーダンスってなんだ？な話
お久しぶりです。鹿部です。
前回の記事からほとんど5ヶ月ぶりですね。終末弾道ネタであれやこれやしてたらまるで面白いネタが見つからなくて全然書くことがありませんでした。
前回とか前々回何を書いてたのかな、と思うとユゴニオ弾性限界の定義とか、修正ベルヌーイの式の具体的な形とかを与えていて、これを使った話を…と書いてるんですが、これはまぁ今回の話ではありません。
前前前回のユゴニオ弾性限界の覚え書きは読んでおくといいかもしれません
今回の話は衝撃インピーダンスってなんだ？というところをササッと。
今回の結論
衝撃インピーダンスと言うのは、物体内を衝撃波が伝播する際、その衝撃波とともに物体の各領域がある速度(粒子速度)で動くが、その時各領域にどれだけの応力が発生するか、を与える係数。
衝撃インピーダンスは物質の密度$\rho$と衝撃波速度$c$の積$\rho c$で表され、衝撃波が弱い時の衝撃インピーダンスは代表的な物質について以下の表の$Z_{s0}$のように表される。
衝撃インピーダンスが出て来るまで
簡潔に言って、音響インピーダンスあるいは衝撃インピーダンスと言うのは、物体内を音波が伝播する時、その音波とともに物体の各領域がある速度(粒子速度)で動く時、その各領域にどれだけの応力が発生するか、を与える係数です。
簡単な例として音叉やガラスのコップのような叩くと音が出るものを考えてみましょう。これらを破壊しないような強さで叩くと、当然ながら音が出ますが、これは叩かれた時に物体内に生じる振動が物体を伝播し、それを受けて空気が一定の振動数で振動させられていることにほかなりません。波動方程式の一般解だとか、そういう面倒な話を取り除き、物体中の各領域が振動する速度(粒子速度)$u(x,t)$が時間$t$とともに$x$軸方向にどう変化するかを単純に次の式で与えます。
\[ u(x,t) =u_0 \sin\left(\omega\left( t - \frac{1} {c} x\right)\right) \]
ここで、$\omega$は角振動数,$c$は物体の音速です。$u_0$は最初に与えられる最大の粒子速度です。
まぁここでは、物体の各領域での粒子速度っていうのは時間と場所で周期的に変化していて音速$c$で伝播していくんだ程度の認識でいます。
物体中の各領域で成り立つ運動方程式は
\[\frac{\partial p}{\partial x} =- \rho \frac{\partial u}{\partial t} \]
で表されます。$\rho$は密度です。なんとなく両辺を$x$で積分すると、右辺の$\frac{\partial u}{\partial t}$というのは$\omega u_0 \cos (\omega ( t- \frac{1}{c} x) )$なので
\[ p = -\rho \omega u_0 \int\left( \cos \left(\omega ( t- \frac{1}{c} x \right) \right) dx \]
を計算すれば良くて、これは$\cos (ax)$の積分が$\frac{1}{a} \sin (ax)$であることから、
\[ p = -\rho \omega u_0 \left( \frac{-c}{\omega}\right) \sin \left(\omega \left( t- \frac{1}{c} x \right) \right) = \rho c u\]
と簡単に表されます。ところで、この$p= \rho c u$に注目してみると、密度$\rho$は定数、音速$c$は弾性領域では定数となり、振動により生じる応力というのは粒子速度の定数倍になります。そういうわけで$\rho c$をひとまとめにしてインピーダンスと呼び、このインピーダンスが高い物質は同じ粒子速度でも高い応力が生じますよ、というのがこの記事で書きたいことです。衝撃波についても、音速$c$を衝撃波速度に置き換えれば同じように衝撃インピーダンスを定義できます。
ちょっと気をつけること
ところで、音速というのはその音のモードによっていくつか種類がありますが、衝撃インピーダンスを音響インピーダンスで近似したいときは静水圧波の音速(ちょっと日本語でなんていうのがわからないんですがBulk sound speedのことです)を使うのが普通です。その理由はユゴニオ弾性限界の覚え書きとユゴニオ弾性限界を超えた後の応力ひずみ線図の覚え書きの通り、弱い衝撃波が生じた時の衝撃波速度は体積弾性率の平方根に比例するからです。
これは少し注意する点で、いわゆる音速で調べると縦波の速度が書かれていることがあります。金属同士の比較であればポアソン比が近いので(比を取れば)そんなに大きな問題にはなりませんが、セラミックスと金属を比較するとポアソン比が大きく異なるので意味がわからない比較になります。
結論
というわけで静水圧波で整理したいくつかの物質の衝撃インピーダンスを以下の表に示します。
この表の出典はC92新刊のWeekend ballisticsで、[2],[3],[8]はそれぞれ
A. Primer, Shock wave compression of condensed matter, (2012), springer
P.J. Hazell, Armor, (2016), CRC press
SP. Narsh, LASL Shock Hugoniot Data, (1980), University of california press.
です。たぶん日本語でも整理されてる本はいっぱいあると思いますが、手元にたまたまあったのがこれらでした。[8]はインターネットで自由に見れるので参考になるかと思います。ここで、$c_0$が体積弾性率から決まる音速で$Z_{s0}$が弱い衝撃波での衝撃波インピーダンスです。これは衝撃波速度の速度が粒子速度に対して$c_0+Su$なる関係を持つとして整理された表です。これらを用いると衝撃インピーダンスの粒子速度が明らかな2つの任意の物質について衝突時の応力を求めることができます。
表を見ると案外セラミックスの衝撃インピーダンスが低いのが興味深いですが、これはセラミックスのポアソン比が低いためにほかなりません(ヤング率とポアソン比の関係はユゴニオ弾性限界の覚え書き参照)。
侵徹との関係
ところで、均一な板を侵徹する場合、衝撃インピーダンスの寄与は小さくL/D=10で1割程度と言われています。これは衝突後純粋に衝撃による界面の移動が生じる時間は極短い時間に限られるためです(この過程、特に除荷側の過程は自分はよく知らないので曖昧ですみません）。その後は前々回のベルヌーイの式とエロージョンを伴う侵徹の覚え書きで示した修正ベルヌーイの式が主であるとよく言われます。
ある種の構造を持った板を侵徹する時、境界ごとに侵徹体は圧力のマッチングをし直す必要があるので(ベルヌーイの式の釣り合いが崩れるので)、衝撃インピーダンスに関連して起こる現象の寄与が大きくなることが想像されます。
まぁなので一次元の侵徹も難しいなぁ、と最近思っています。いい本募集中です。
登録:
投稿 (Atom)