Siegelモジュラー形式環と計算機代数

Rungeという人の1993年の論文

On Siegel modular forms. Part I.
https://doi.org/10.1515/crll.1993.436.57

に、Siegelモジュラー形式環の構造が、有限群の不変式の計算と、環の商体内での正規化(整閉包)によって計算できるという結果が書いてある(Theorem3.12)。有限群が作用する環は、テータ関数(テータ零値)の多項式たちで生成され、単純な多項式環というわけではない。偶数ウェイトのSiegelモジュラー形式環は、もう少し、簡単な計算でいける(Corollary3.17)。


Siegelモジュラー形式は1930年代終わり頃、Siegelが導入した。Siegelモジュラー形式に関する概略が、以下の論説の1〜2節などに書かれてる。
コンパクト化の今昔
https://www.jstage.jst.go.jp/article/sugaku1947/51/2/51_2_129/_article/-char/ja/

歴史的に見ると、(19世紀頃の)初期のモジュラー形式論は、個々のモジュラー形式を理解しようという特殊関数論的(というか古典解析的)な動機が強かったように見える。一方で、Siegelモジュラー形式環を調べる動機の一つは、離散群による商空間という扱い難い対象を理解しようというところにあって、解析的側面は薄れてて、代数幾何学的な話になっている。

【モジュライ空間史】On the early history of moduli and Teichmüller spaces
https://arxiv.org/abs/1602.07208
を読んだところ、Siegelは、1935年の論文で、Torelli空間とSiegel上半平面、Siegel基本領域を定義し、Torelli空間から、Siegel上半平面への複素構造を保つ写像を定義したということが書いてあった。Torelli空間は、おそらく、Weilによる命名で、今では、タイヒミュラー空間のTorelli群による商と定義されることが多い。タイヒミュラー空間の論文は、この少し後に書かれたらしいけど、TeichmüllerがSiegelの仕事を知ってたのかは謎。論文は読んでたとしても、当時、Teichmüllerは、ナチの反ユダヤ主義を支持し、Siegelは、ユダヤ人だったのかは知らないけど、1940年にアメリカに亡命してるので、仲良く交流したとは思えない。それより前の時代だと、Schottkyなんかは、純粋に解析的なアプローチを取っていたようだし、1913年にTorelliの定理を示したTorelliも、周期行列の集合に良い幾何構造が入るという認識は持ってなかったっぽい。1920年代に、対称空間の分類とかが行われてて、Siegelに、その影響があったかは不明。モジュライ空間を厳密に定式化するという点で、Siegelは大きな貢献があったようである



私は、コンパクト化の議論とか読んでも、そんなんええから、Siegelモジュラー形式環の生成元と関係式を書いてくれへんかなという気持ちになる。現在、種数3までのSiegelモジュラー形式環の構造が決定できているっぽい。種数1のSiegelモジュラー形式環は、19世紀には、既に知られていた通り、4次と6次のEisenstein級数で生成される。この環と、有限群の不変式環との間に自然な同型が存在するというのは、誰が最初に気付いたのか分からないけど、有限群の不変式環とある種のモジュラー形式環に自然な同型があるという現象は、1970年頃には知られてたっぽい(もっと古い起源があるのかは不明)。

例えば、1972年のフランス語の論文

Polynômes des poids de certains codes et fonctions thêta de certains réseaux
http://www.numdam.org/item/ASENS_1972_4_5_1_157_0/

には、位数192の群の不変式環が、\mathbf{C}[E_4,\Delta]と自然に対応するということが書かれてるように見える(フランス語は読めない)。

そういう話とは関係なく、Siegelモジュラー形式環を頑張って決めようとした人もおり、種数2のSiegelモジュラー形式環は、井草(1962)
On Siegel modular forms of genus two
https://www.jstor.org/stable/2372812

種数3のSiegelモジュラー形式環は、露峰(1986)の結果がある。
On Siegel Modular Forms of Degree Three
https://www.jstor.org/stable/2374517?seq=1

種数3の場合は、計算が大変なようで、生成元がminimal setかどうか分かってなかったらしいけど、2019年に、minimal setを決めたと主張している論文が出ているよう(内容は未精査)。19個の生成元と55個の関係式があるらしいので、結果をトレースするだけでも、手計算は辛そう

Siegel modular forms of degree three and invariants of ternary quartics
https://arxiv.org/abs/1907.07431


Rungeの結果を使って、モジュラー形式環を決定するには、テータ零値の間に存在する代数的関係式を決める必要があるけど、一般的な方法は知られてないっぽい。m=(m',m'')\in \{0,1\}^g \times \{0,1\}^gに対して
\theta_m(z,\tau) = \displaystyle \sum_{x \in \mathbf{Z}^g} \exp 2\pi i\left(\dfrac{1}{2} {}^{t}[x+\dfrac{1}{2}m'] \tau [x+\dfrac{1}{2}m'] + \langle x+ \dfrac{1}{2}m',z+\dfrac{1}{2}m'' \rangle \right)
とすると、zについて偶関数であるものが、2^{g-1}(2^g+1)個ある。z=0の時しか考えないので、oddな方は、どうでもいい。Rungeの論文では、テータ関数の二乗が、
f_a(\tau) = \theta_{(a,0)}(0,2\tau) = \displaystyle \sum_{x \in \mathbf{z}^g} \exp 2\pi i \left( {}^{t}[x+\dfrac{a}{2}] \tau[x+\dfrac{a}{2}] \right)
を使って書けることから、この形のテータ零値の作る環を中心に考える。この形のテータ零値は、どれも0でなく、2^g個の線形独立な関数となっている(Rungeの論文p59)

テータ関数やテータ零値には、沢山の恒等式があって、関係式の間の関係式が必要なほど、ややこしい。必要なテータ零値の間の関係式も、知られている関係式から得られるものもある。Rungeの論文§6では、種数3の場合に、既知のテータ関係式から、f_aの関係式を得て、Siegelモジュラー形式環の構造を決定している(Theorem6.2)

cf)A remark on a theorem of Runge
https://link.springer.com/article/10.1007/s000130050220
[PDF]https://www.mathi.uni-heidelberg.de/~freitag/preprints/runge.pdf
の式(2)も参照


もっと高い種数では、既知のテータ恒等式に帰着できない関係式が出るのか、よく分からない。

A theta relation in genus 4
https://projecteuclid.org/euclid.nmj/1114631553
では、種数4の場合に、一つの関係式を見つけて、"It is a natural question, whether our relationis a consequenceof Riemann's theta relations. We could not decide this."と書いてる(2001年付けだから、少し前だけど、凄く古いわけではない)。全部が既知のテータ恒等式に帰着できるなら、原理的には、変数消去の計算で、全ての関係式を見つけることができる。実際上は、evenなテータ関数が2^{g-1}(2^g+1)個あって、(τではなく)2τの関数であるf_aが、2^g個あるので、単純に考えると、g=4で152変数、g=5で560変数の多項式の変数消去計算をやらないといけないことになって、現在のコンピュータとアルゴリズムで計算できるかは怪しいけど

【参考文献】arXiv:1401.5368やarXiv:1510.02699には、色々なテータ恒等式の発見に、誰が寄与したかという話が少し書いてある

【テータ恒等式とPlucker関係式】2x4行列の2x2小行列式D_{ij} (1 \leq i \lt j \leq 4)で生成される環は、多項式環を、一個のPlucker関係式で生成されるイデアルで割ったものになる。
\mathbf{C}[D_{12},D_{13},D_{14},D_{23},D_{24},D_{34}] \simeq \mathbf{C}[X_{12},X_{13},X_{14},X_{23},X_{24},X_{34}]/(X_{12}X_{34}-X_{13}X_{24}+X_{14}X_{23})
行列式D_{ij}に対して、4x4対角行列h \in H_4 \subset GL(4)の作用を
\rho(h)(D_{ij})(x) =D_{ij}(xh)
で定める。\rho(h)(D_{12}D_{34}) = det(h) D_{12}D_{34}etc.などが成立するけど、h \in H_4の作用で、det(h)多項式倍だけ変化する元の全体は環をなし、
\mathbf{C}[D_{12}D_{34} , D_{13}D_{24} , D_{14}D_{23}] \simeq \mathbf{C}[y_0,y_1,y_2]/(y_0-y_1+y_2)
のようになる。一般の群と指標に対して、この手の環は、半不変式環と呼ばれてることがある。一方、種数1のテータ零値には
\theta_{00}^4 = \theta_{10}^4 + \theta_{01}^4
という関係式があるので、環同型\mathbf{C}[\theta_{00}^4,\theta_{10}^4,\theta_{01}^4] \simeq \mathbf{C}[D_{12}D_{34} , D_{13}D_{24} , D_{14}D_{23}]が成立する。

代数的に書くと、大した意味はなさそうだけど、これは偶然でなく、19世紀から知られる超幾何関数と楕円関数の不思議な繋がりが背景にある。\mathbf{C}[\theta_{00}^4,\theta_{01}^4,\theta_{10}^4]は、モジュラー群\Gamma(2)のモジュラー形式環M(\Gamma(2))になっている。このモジュラー形式環は、楕円曲線単位元以外の等分点e_1,e_2,e_3を使うと、
M(\Gamma(2)) = \mathbf{C}[e_1 , e_2 , e_3]
のようにも書ける。テータ零値で書くと
e_1=-\dfrac{1}{3}(\theta_{00}^4 + \theta_{10}^4)
e_2=\dfrac{1}{3}(\theta_{00}^4 + \theta_{01}^4)
e_3=-\dfrac{1}{3}(\theta_{01}^4 - \theta_{10}^4)
で、e_1+e_2+e_3=0が成立する。SL(2,\mathbf{Z})/\Gamma(2) \simeq S_3は、e_1,e_2,e_3の置換として作用し、M(SL(2,\mathbf{Z})) \simeq M(\Gamma(2))^{S_3}で、M(SL(2,\mathbf{Z}))の生成元として、e_1e_2+e_2e_3+e_3e_1e_1e_2e_3が取れて、それぞれ、定数倍を除いてEisenstein級数と一致する(E_4=-3(e_1e_2+e_2e_3+e_3e_1) , E_6 = \dfrac{27}{2}e_1e_2e_3)(逆に、M(\Gamma(2))M(SL(2,\mathbf{Z}))上の自由加群になってて、この構造を具体的に書いてみるのも面白い)。一方、\mathbf{C}[D_{12},D_{13},D_{14},D_{23},D_{24},D_{34}]のProjは、Grassmann多様体Gr(2,4)で、これをH_4で"簡約"した空間が、Gauss超幾何関数の住む空間と思える。

"私説超幾何関数"(ISBN:978-4-320-01576-0)という本を読むと、種数2でも類似の同型が作れると書いてある。今度は、Grassmann多様体Gr(3,6)を考える。3x6行列の3x3行列式D_{ijk}(1 \leq i < j < k \leq 6)は、20個あって、代数的トーラスのdetに関する半不変式環Z(3,6)は、10個の生成元と5個の独立な関係式からなる。種数2のテータ零値は10個あるけど、テータ零値の2次式を10個うまく選ぶと、それらの二乗から生成される環(適当なモジュラー群のモジュラー形式環となるらしい)と半不変式環Z(3,6)とが同型になるようだ。種数3以上に一般化する話は知られてないっぽいし、根源的な理由も分からないので、モジュラー形式というより、テータ零値の恒等式(の一部)に、何か良い構造があるという可能性もあるんじゃないかと勝手に思ってる



それ以外のステップ。有限群の不変式環を計算するアルゴリズムは以前に実装した

有限群の不変式環の生成元を計算するアルゴリズム
https://m-a-o.hatenablog.com/entry/20130316/p2

asir_misc/finvar.rr
https://github.com/vertexoperator/asir_misc/blob/master/finvar.rr

にも、Risa/Asirで実装したものが置いてある。しかしまぁ、過去の私の記述によれば、"有限群の不変式の計算は、SingularやMagamaにも実装されてるらしいので、計算するだけなら、それ使うのが一番確実かもしれない"らしいので、Singular使えばいい気がする。

正規化の方は、Risa/Asirには、ないので、使いたければ自分で実装するしかない(と思う)けど、Singularには、既に実装されている

D.4.18 normal_lib
https://www.singular.uni-kl.de/Manual/4-0-3/sing_1269.htm#SEC1344

つまり、総合的に見て、Singular使えばいいのでは。一応、Singularで使ってる正規化アルゴリズムは、以下の論文に書かれているようだ。

Normalization of Rings
https://arxiv.org/abs/0904.3561

以前測定した限りでは、グレブナー基底の計算は、Singularの方が、Risa/Asirより、速い(数倍くらい違うこともあった気がする)。ただ、Risa/Asirは、何故か、代数的数のサポートが割と手厚い気がするので捨てがたい。あと、Singularは、(Mac以外)インストールめんどくさい



種数1の場合は、容易に計算できるし、Rungeの論文では、Remark3.6の直後に、結果が書いてある。種数1の場合は、群がShephard-Toddの(既約擬鏡映)群No.8というものになる

cf)List of irreducible complex reflection groups
https://en.wikipedia.org/wiki/Complex_reflection_group#List_of_irreducible_complex_reflection_groups

【用語について】正方行列Aが擬鏡映とは、Aが有限位数で、rank(I-A)が1であることを指す。英語だと、pseudo-reflection。位数が2の時は、ある超平面に関する鏡映になる。同じ概念に対して、complex reflection,unitary reflectionとかいう呼び方をされることもある。擬鏡映は、条件から、1でない固有値を一つだけ持ち、それは、1のべき根になってるから、こういう名前が付いてるのだと思う。+1でない固有値が実であれば、-1しかありえず、Aは鏡映になり、普通の鏡映は、実鏡映ということになる。個人的には、擬鏡映という名前が一番わかり易い。擬鏡映で生成される群が、擬鏡映群


私の実装でも

Z2 = newalg(x^2+1);
T=newmat(2,2,[[1,1],[1,-1]])*(1+Z2)/2;
D=newmat(2,2,[[1,0],[0,Z2]]);

finvar(grp([T,D]),[f0,f1]);

とかやって、一瞬で正しい答えが出た。
f_0 = \theta_{00}(2\tau) , f_1 = \theta_{10}(2\tau)
で、
f_0^2 = \dfrac{\theta_{00}^2 + \theta_{01}^2}{2}
f_1^2 = \dfrac{\theta_{00}^2 - \theta_{01}^2}{2}
が成立し、 2つの独立な不変式は
g_4 = f_0^8+14f_1^4 f_0^4+f_1^8
g_6 = f_0^{12} - 33f_0^4f_1^4(f_0^4+f_1^4)+ f_1^{12}
となる。Eisenstein級数
E_4 = \dfrac{1}{2} ( \theta_{00}^8 + \theta_{10}^8 + \theta_{01}^8)
E_6 = \dfrac{1}{2} (-3\theta_{10}^8(\theta_{00}^4+\theta_{01}^4) + \theta_{00}^{12} + \theta_{01}^{12})
なので、両者が一致することは、単なる計算問題


種数2の場合。群H_2の生成元は、整係数Symplectic群の生成元から来ていて、Jから来る方はいいけど、整係数対称行列は、無限個ある。しかし
S = \left( \begin{matrix} a & b \\ b & c \end{matrix} \right)
と置くと
D_S = diag(1 , i^c , i^a , i^{a+2b+c})
になって、a,b,cは、0~3の範囲を舐めれば十分。有限個の生成元が決まったので、群の要素を列挙しようと思ったけど、pythonのsetに相当するものがないRisa/Asirの実装では、まる一日かけても、46080個の元を列挙するのは無理だった。仕方ないので、pythonで列挙して、ファイルから群の要素を読み込むという方法で対応(計算時間は、夜寝ている間に終わる程度)。

Rungeの論文には、これは擬鏡映群だと書いてあるけど、群に含まれる擬鏡映をコード書いて列挙してみると60個ある。位数は全て2で、どうやら、Shephard-Toddの群No.31っぽい(ちゃんとチェックはしてないけど、擬鏡映を列挙できてるので、比較すれば、力ずくで確認できるはず)。Molienの公式でHilbert級数を計算する(これはRisa/Asirで実装してある。昼寝してたら終わってた)と、
\dfrac{1}{(1-t^8)(1-t^{12})(1-t^{20})(1-t^{24})}
になっている。なので、この群を、4変数多項式環に作用させた不変式環は、8,12,20,24次の代数的に独立な不変式で生成される(8*12*20*24=46080で、群の位数が得られる)。具体的な式は、Rungeの論文p77冒頭に記載されているけど、Reynolds作用素の計算を、コンピュータにやってもらえばできる。論文にも書いてあるけど、20次と24次の不変式の空間は、一次元ではないので、定数倍を除いても、論文と同一の式にならない可能性はある。というか、私はならなかった。

モジュラー形式としては、ウェイト4,6,10,12に対応し、(Corollary3.17によって)種数2、偶数ウェイトのSiegelモジュラー形式環は、この4つの元で生成される。奇数ウェイトのモジュラー形式として、ウェイト35のものがあり、この5つの元が、種数2のSiegelモジュラー形式環を生成する


種数3の場合は、Rungeの群の位数が約3.7億なので、全部要素を列挙して計算するというのは望みがない。多分、私のPCでは、メモリに載らない。詳細は、
On Siegel modular forms. II
https://projecteuclid.org/euclid.nmj/1118775400
の方にあって、群の構造を調べて何とかしてるようだ。

現在のアルゴリズムは、Molienの公式やReynolds作用素の直接計算に依存してるので、位数の大きな群に対して使うには、どうしても限界がある。けど、種数が大きくなるごとに試行錯誤してたらキリがないので、超でかい有限群の不変式と戦えるアルゴリズムが欲しい。次元が上がると、今度はグレブナー基底の計算が辛くなるという二重苦が待っている

解析力学から見た古典/量子Kepler系の解法

古典Kepler系の解析解
https://vertexoperator.github.io/2016/12/01/kepler2d_classical.html

2次元量子Kepler系の束縛状態
https://vertexoperator.github.io/2016/12/01/kepler2d_quantum.html

で、2次元のKepler系を、古典力学量子力学で解いている。多少見かけは違ってるけど、同じ方法で空間変数の取り直しを行って、調和振動子に帰着(全く同じなので、書いてないけど、正エネルギー/散乱状態の場合は、inverted harmonic oscillatorというものに帰着する)、周波数の形も古典/量子で同じ形になる

ところで、古典Kepler系の方は、時間変数も取り直しているのに対して、量子Kepler系では、時間変数の変換は出てきてない。量子系の方は、時間非依存のシュレディンガー方程式を考えているのだから当然だけど、量子Kepler系のHamiltonianに対して、時間依存シュレディンガー方程式を考えて、時間変数の取り換えを行ってから、(新しい時間変数に対する)定常シュレディンガー方程式を考えると、結局、Hamiltonianは、左からrを掛けるという変換を受けることになる。

"量子Kepler系の束縛状態"の解法2の方で、「量子Kepler系のシュレディンガー方程式
r(\hat{H} - E)\Psi = 0
と置くと・・・」と書いてあるのは、自明な式変形を書いてるのではなく、シュレディンガー方程式をこのように見ることが、古典系に於ける時間変数の取り換えと同等の操作を行うことに相当している。また、spectrum generating algebraに含まれるのは、Hamiltonianそのものではなく、r(\hat{H} - E)であることからも、この変換の必要性が理解できる。

ただ、これだけだと、ややad hocな印象もある



Duru-Kleinertによる、量子Kepler系の経路積分に関する以下の論文には、"pseudo Hamiltonian"というものが書いてある(論文中の式(10))。彼らは、これ以前にも同じような論文を書いてるけど、そっちには、このような説明はないので全然参考にならない

Quantum Mechanics of H-Atom from Path Integrals
https://doi.org/10.1002/prop.19820300802
[PDF]http://users.physik.fu-berlin.de/~kleinert/83/83.pdf

r(\hat{H} -E)は、Duru-Kleinertが"pseudo Hamiltonian"と呼んでいるものに相当する。


で、Duru-Kleinertの論文の式(15)を見ると、このpseudo Hamiltonianが住んでるのは、拡大相空間であることが読み取れ、pseudo Hamiltonianに対して正準方程式を考えると、普通の運動方程式+時間変数の取り換えの式が出ている。Duru-Kleinertは、拡大相空間とは言ってないけど、昔の私が調べた所によると、1949年にLanczosが導入した比較的新しい概念らしい。昔の私が書いてることを信じるなら、H-Eはextended Hamiltonianと呼ばれることもあるそう

cf)extended phase space
https://formalgroup.tumblr.com/post/105356647155/extended-phase-space

古典力学では、H-E=0は定義みたいなものであり、量子力学では、H-Eを波動関数に作用させると0になるというのは、シュレディンガー方程式そのもの。


拡大相空間を考えるなら、(明らかに、拡大相空間も、自然なsymplectic構造を持ってるので)正準変換も"拡大"されるべきと思うのだけど、

Stackel系の全ての保存量を保つ離散化
https://repository.kulib.kyoto-u.ac.jp/dspace/handle/2433/42747

とかいう関係なさそうな文献を見ると、そのような変換を考えた人はいるらしい。Duru-Kleinertのpseudo Hamiltonianは、extended Hamiltonianに、一般化した正準変換を施したものであると理解できる。

DuruとKleinertは、経路積分のために、pseudo Hamiltonianを考えたけど、経路積分自体は、どうでもよく、背景に、こういう構造があることが重要なのだと思う



古典力学では、ハミルトニアンは関数なので、右から関数を掛けるか左から関数を掛けるか気にする必要はないけど、量子力学では、事情が違う。Kepler系では、r(\hat{H}-E)としたのだけど、(\hat{H}-E)rでない理由があるだろうかという疑問が残っている。一般には
f_{l}(\hat{H}-E)f_{r}
という形の変換が考えられる。この形の変換を取り扱った文献で、一番古いのは、調べた限り

On the generalization of the Duru-Kleinert-propagator transformations
https://link.springer.com/article/10.1007/BF01318171

じゃないかと思う。

Kleinertの本が、Kleinert自身によって公開されているけど、

New Path Integral Formula forSingular Potentials
http://users.physik.fu-berlin.de/~kleinert/b5/psfiles/pthic12.pdf

この形の変換が扱われている。式(12.26)とか。本を読めばわかるけど、新しい時間変数と、時間変数の関係は、f_lf_rの積のみで決まる。式(12.37)など。


この形の変換で一番良く見るのは、相似変換で、その場合は、時間変数を変える必要はない。相似変換については、当たり前すぎて、注意を払われることもないけど、それより、ずっと一般的な変換を行っても、ツケを時間変数に押し付けてしまえるというのは、言われれば当たり前だけど、何かに使えることもあるかもしれない。Kepler系の場合は、rを掛けることによって、Coulombポテンシャルの特異性が見かけ上消えるという見方もできる


【おまけ】Kleinertは、1967年に、A.O. Barutと共に
Transition Probabilities of the Hydrogen Atom from Noncompact Dynamical Groups
https://doi.org/10.1103/PhysRev.156.1541

という論文を書いていて、3次元量子Kepler系の束縛状態の空間が、so(4,2)の既約表現空間になってることを、(私が知る限り)初めて明確に書いた人でもある。この論文には、まだ"pseudo Hamiltonian"は書かれていない。もう少し後のBarutの論文では、Hamiltonianそのものより、"pseudo Hamiltonian"(とは呼んでないけど)を重視する書き方をしているものがあったと記憶している

Group Theory and Orbital Fluctuations of the Hydrogen Atom
http://users.physik.fu-berlin.de/~kleinert/kleiner_re222/barut92.html

のAbstractで、Kleinertは、以下のように書いている。

The group theoretic development was triggered by A.O. Barut who suggested to me the search for a dynamical group larger than SO(4). In this way he became partly responsible also for the later and most recent path integral development.

Chapman-Jouguet理論

物理で定義とか真面目に考え出すと、辛いことになるけど、"法線方向の速度が不連続に変化する波面の伝播"を衝撃波と総称する立場からすると、発熱反応を伴う圧縮衝撃波がデトネーション波ということになる。発熱反応を伴う衝撃波では、膨張波も存在することができ、こっちは、デフラグレーション波と呼ぶ。日本語では、デトネーションは爆轟と訳され、デフラグレーションは爆燃と訳されている。

※)接線方向の速度が不連続に変化する不連続面は、渦層と呼ばれ、衝撃波面とは異なる。同じ不連続面とはいえ、使われる手法とかも全然別のよう。

普通の衝撃波は、化学反応とかは伴わないので、そのへんがデトネーション/デフラグレーションとの違い。デフラグレーションは、不連続波面を持つけど、伝播速度は音速より遅いので、衝撃波のイメージに、そぐわない。まぁ、あまり深く考えず、気相爆轟=超音速で伝わるガス爆発くらいの認識でいいっぽい。

気相爆轟は、19世紀の終わり頃から研究されるようになり、理論的な研究は、Mikelson(1890),Chapman(1899),Jouguet(1905)などが独立に行ったらしい。これらは、同じ結果を与えるらしく、現在では、Chapman-Jouguet理論(以下、CJ理論)として知られる。現在、CJ理論は、衝撃波のRankine-Hugoniot理論を拡張した形で展開される。

どうでもいいけど、CJ理論のChapmanとChapman-Enskog展開のChapmanは別人。Chapman-Enskog展開の人は、Sydney Chapman(1888-1970)という名前で、CJ理論の人は、David Chapman(1869-1958)らしいけど、David Chapmanには、Sydney John Chapman(1871-1951)という弟がいて経済学者だったらしい。ややこしいわ

CJ理論では、波面の前後が、それぞれ未燃ガス、既燃ガスになっていて、反応は一瞬で終了するものとされている。未燃ガス、既燃ガスは共に、理想気体として記述するのが一般的。基本的に、代数的な計算しか出てこないので、難しいわけでないけど、既燃ガスの比熱比は常温のものと違うことに気付かず、計算してみると、値が合わないという罠にはまった

後になって気付いたけど、デトネーションのことは、ランダウ・リフシッツの流体力学の教科書に載っていて、CJ理論が説明されてる。この本も、測定値と理論値の比較とかはないので、助けにはならなかったろうけど。

Fluid Mechanics §121
https://archive.org/details/FluidMechanics/page/n491



罠にはまった時に、Chapmanの原論文を読んだ。

On the rate of explosion in gases
https://www.tandfonline.com/doi/abs/10.1080/14786449908621243

Chapmanの議論は、現在の定式化とは違っていて分かりにくいけど、TABLE IIIに、"Velocity of Explosion"の式として

V^2 = \dfrac{2RJ}{\mu C_{v}^2} \left( {(m-n)C_p + m C_{v} }C_{p} t_0 + (C_p+C_v)\mathrm{h} \right)

が書かれていて、様々な混合気体の燃焼反応に対して、Vを計算し、測定値と比較している。単位とかが書かれてないけど、私が解読した限りでは、多分、以下のようになっている。

Chapmanの論文には、多くの混合気体に対する測定値と計算値が載ってるけど、全部やるのは大変なので、水素と酸素がモル比で1:1に混合した気体の燃焼反応(2モルの水素と2モルの酸素を燃焼させて、1モルの水蒸気が発生し、1モルの酸素が反応せず残るという状況)
2\mathrm{H}_2 + 2\mathrm{O}_2 = 2 \mathrm{H}_{2}\mathrm{O} + \mathrm{O}_2
を例に取ると、

RJ = 8.314462618(J/mol/K) : 気体定数
μ = 68(g) : 単位量あたりの気体の質量(反応前後で保存されるはずなので、水素2モル+酸素2モルの質量であり、また、水2モルと酸素1モルの質量でもある)
m = 3(mol) : 反応前のモル量(水2モルと酸素1モル)
n = 4(mol) : 反応後のモル量(水素2モルと酸素2モル)
h = 117000(cal) : 反応熱=1molのH2O(g)につき、245kJ弱の計算
t0 = 273.15(K) : 基準温度
Cv = 11.63(cal/mol/K) : 生成ガスの定積モル比熱(表に記載されてる値)
Cp = 13.62(cal/mol/K) : 生成ガスの定圧モル比熱(Cvから、マイヤーの関係式を使って計算)
V = 2341.8(m/s) : Chapman-Jouguet速度(論文では2340と計算されている)

だと思う。


生成ガスの比熱について、論文には、以下のようなexcuseが長々と書いてあって、時代を感じる

A few words are necessary regarding explosions in which water is formed. If the specific heat of steam is taken as 3/2 × specific heat of the diatomic gases, the found rates of explosion fall below the calculated rates when the dilution with inert gas is great, and vice versd when the dilution is small.

It is possible to account for this by two theories. The first theory is that at high temperatures the water is dissociated, whereas at low temperatures the combination of hydrogen and oxygen is complete. The second theory is that the specific heat of steam rises more rapidly with the temperature than the specific heat of the diatomic gas.
(中略)
We are therefore encouraged to test the first theory, i.e. that the specific heat of steam rises more rapidly with the temperature than that of the diatomic gases.

水蒸気の比熱を、2原子分子の1.5倍とすると、うまくいかないので、高温では、水蒸気の比熱は、2原子分子より急速に上昇するのだろうとか。Chapmanが書いている1.5倍というのは、多分、エネルギー等分配則の誤った適用による(比熱が一分子あたり原子数に比例すると考えたんだろう)ものだろう。Chapmanは、もう一つの可能性として、現在でいうところの弱電離プラズマのようなものを考えて却下している。



反応直後のガスは、高圧で、温度は数千度に達するらしいので、常温の場合と違って、水や酸素でも、分子振動の比熱への寄与が無視できなくなる。分子内振動を調和振動子で近似して、独立なモードに分離すると基準振動が得られる。基準振動波数は、赤外分光やラマン分光で測定でき、水素分子の基準振動波数は、ν=4160(/cm)程度。分光業界(?)では、波数の記号は、νの上にチルダを付けたものを使うことが多い気がするけど、以下では、普通のνを使う。回転スペクトルを測るマイクロ波分光はGHzを使ってるっぽいので、赤外/ラマンも、THz単位に揃えろやと思うけど。水の場合、適当に検索した数値によると、1595(/cm),3756(/cm),3657(/cm)。

世の中には合成法が確立してなかったり合成難度が高い化合物もあるけど、その場合でも(低分子なら)割と高速かつ、それなりの精度で、基準振動波数を数値計算することもできる。例えば、現時点で、合成法が確立してない低分子の一つに、テトラヘドランC4H4があり、分子量などからすると、常温常圧で気体の可能性がある(構造異性体の1,2−シクロブタジエンや、1,3-シクロブタジエンの沸点が40度前後という情報がある)

cf)Tetrahedrane
https://en.wikipedia.org/wiki/Tetrahedrane

振動モードが一個の場合、分子振動の分配関数は
q_{\mathit{vib}} = \left( \dfrac{e^{-\beta c \mathit{h} \nu/2}}{1 - e^{-\beta c \mathit{h} \nu}} \right)^N
なので、一個の振動自由度による定積比熱への寄与は
C_{v,\mathit{vib}} = N k_{B} \dfrac{x^2 e^{x}}{(e^{x}-1)^2}
と計算できる。但し、x = \dfrac{c h \nu}{k_{B}T}

ν=4160(/cm)、T=300(K)だと、x=19.95で、y=\dfrac{ C_{v,\mathit{vib}} }{N k_{B}}は、8.62e-7とほぼ無視できる。T=3000(K)だと、x=1.995で、y=0.725になって、無視できない寄与があることが分かる。

十分高い温度では、振動の自由度1つにつき、モル比熱は気体定数Rだけ増加する。振動自由度を高温近似で扱うと、高温の二原子分子理想気体では、定積比熱が3.5R,定圧比熱が4.5Rになり、高温の非直線型三原子分子理想気体では、定積比熱は6R,定圧比熱は7Rになる(はず)。

6R=50(J/mol/K)=12(cal/mol/K)程度で、上のChapmanの計算に出ている11.63(cal/mol/K)に近い値になっている。実際に上の計算で扱ってるのは、水蒸気2モルと酸素1モルの混合気体の定積比熱なので、それぞれ理想気体として、高温近似が成立するならば、定積比熱は、10.27(cal/mol/K)程度という計算。完全に一致とは言えないけど、それなりに合ってそうではある


100年以上前の論文の数値と理論値だけを頼りにするのもどうかと思うので、もうちょっと新しいデータを見ておく(数千度とかで比熱を測定できるのか知らないけど、以下のデータは、測定値というわけではない)

Water@NIST Chemistry Webbook
https://webbook.nist.gov/cgi/cbook.cgi?ID=C7732185&Units=SI&Mask=1&Type=JANAFG&Table=on#JANAFG

を見ると、水の6000(K)に於けるCpの値は、60.59(J/mol/K)となっていて、これは、7R=58.20(J/mol/K)と近い。4100(K)では、丁度、58.20(J/mol/K)となっている。二原子分子気体の例として、H2(g)を見ると

Hydrogen@NIST Chemistry Webbook
https://webbook.nist.gov/cgi/cbook.cgi?ID=C1333740&Units=SI&Mask=1&Type=JANAFG&Table=on#JANAFG

に、3200(K)で、Cp=37.53(J/mol/K)とある。これは、4.5R=37.415(J/mol/K)と、ほぼ一致している。上の基準振動の数値に基づく計算を考慮すると、3200(K)で、この数値に達するのちょっと早くね?という気もするけど。6000(K)に於けるCpの値は、41.97(J/mol/K)となっている。

この手の小さな多原子分子では、3000〜4000(K)で、比熱への振動自由度の寄与を高温近似で扱うことが許されそうにも思える。幸い、デトネーションに於ける生成ガスの温度は、この水準にある

これらのデータが、どれくらい信用できるものなのか分からないけど、並進、回転、分子内振動を単純に考慮しただけのモデルとは少しズレがあるけど、実験値と比較しないことには、何も言えそうにない。


【おまけ】大学一年の時、エネルギー等分配則による気体の比熱の説明を習った際、結構ずれてる気体もあって胡散臭いと思った記憶がある(条件とか分子は忘れた)。ファインマンの教科書の

The Feynmann Lectures on Physics Vol.I 40-6 The failure of classical physics
http://www.feynmanlectures.caltech.edu/I_40.html#Ch40-S5

を見ると、古典統計力学で気体の比熱を説明できないことは、丁度Chapmanの論文が書かれた時代の物理学者も不可解に思ってたと書いてある。そんな説明を21世紀にもなって教えないで欲しかった。ファインマンは、分子振動を調和近似して量子統計力学で扱えばいいと言ってるけど、計算例は一つもないし、他で見たこともないので、やってみたのが、以下の表。

分子 Cp(J/mol/K) Fr Cp/R-2.5-Fr/2 基準振動(/cm) X R*(2.5+Fr/2+X)
N2(g) 29.124 2 0.0028 2358.57 0.001477 29.113
O2(g) 29.376 2 0.0331 1580.19 0.028395 29.337
HCl(g) 29.136 2 0.0043 2990.95 0.000112 29.102
HI(g) 29.156 2 0.0067 2309.01 0.001798 29.116
Br2(g) 36.048 2 0.8356 325.32 0.8177 35.899
CO2(g) 37.129 2 0.9656 [667,667,1388,2349] 0.9565 37.053
N2O(g) 38.617 2 1.1446 [589,589,1285,2224] 1.1428 38.603
NO2(g) 36.974 3 0.4470 [750,1320,1617] 0.4653 37.127
O3(g) 39.238 3 0.7192 [719.6 ,1062.1 ,1157.8] 0.6745 38.866
H2S(g) 34.192 3 0.1124 [1183,2615,2626] 0.1098 34.171
NH3(g) 35.652 3 0.2880 [950,1627,1627,3337,3444,3444] 0.2671 35.479
ND3(g) 38.225 3 0.5974 [793,1225,1225,2495,2652,2652] 0.5254 37.626
PH3(g) 37.102 3 0.4623 [992.1,1118.3,1118.3,2322.9,2327.7,2327.7] 0.4656 37.130
BF3(g) 50.446 3 2.0673 [480.7,480.7,696.7,888,1463.3,1463.3] 2.0673 50.446
C2H2(g) 44.095 3 1.3034 [612,729,729,1974,3289,3374] 1.2929 44.077
CH4(g) 35.639 3 0.2864 [1306,1306,1306,1534,1534,2917,3019,3019] 0.2864 35.639
CH3Cl(g) 40.731 3 0.8988 [732,1017,1017,1355,1452,1452,2937,3039,3039] 0.8998 40.739
SF6(g) 96.960 3 7.6616 [346,346,346,524,524,524,...](略) 7.6648 96.987
C3H8(g) 73.60 3 4.8520 [216, 268, 369, 748,...](略) 4.524 70.874

定圧モル比熱Cpの出典は、JANAF熱化学データ表かNIST Chemistry WebBookで、298.15(K) , 0.1(MPa)の値。Frは回転の自由度で、Rは気体定数。ND3は、アンモニアの水素を全部重水素に置換したもの。

表の値Xは、全基準振動波数について、以下の和を計算したもの(温度Tは、298.15K)
X = \displaystyle \sum_{i} f(\dfrac{c h \nu_i}{k_{B}T})
f(x) = \dfrac{x^2 e^{x}}{(e^{x}-1)^2}
並進と回転運動からのみ比熱への寄与がある理想気体を高温近似で考えると、Cp=(5+Fr)R/2になる。分子内振動からの比熱への寄与を量子統計力学に基づいて計算し、加えると、Cp=(2.5+Fr/2+X)*Rになる。これを測定値と比較すればいい

基準振動の出典は、以下のどれか
(1)JANAF熱化学データ表
(2)Vibrational Modes in HITRAN2012 Database
https://www.cfa.harvard.edu/hitran/vibrational.html
(3)Tables of molecular vibrational frequencies, consolidated volume I (Takehiko Shimanouchi)
https://nvlpubs.nist.gov/nistpubs/Legacy/NSRDS/nbsnsrds39.pdf
https://archive.org/details/tablesofmolecula39shim

(2)のアセチレンのデータは、縮退度にミスがあるので修正してある。また、六フッ化硫黄の基準振動は(2)から、
[346,346,346,524,524,524,615,615,615,643,643,775,948,948,948]
を使い、27個あるプロパンの基準振動は、(3)に従って、以下の数値を使った
[216,268,369,748,869,922,940,1054,1158,1192,1278,1338,1378,1392,1451,1462,1464,1472,1476,2887,2887,2962,2967,2968,2968,2973,2977]

計算値は、これなら、過去の自分も納得しただろうと思う程度には合ってる。どっちみち、調和近似が完全に正確な結果を与えるはずはないので、多少のずれは、あって然るべきもの



Chapmanの論文にあるCJ速度の式は、実際の数値を見ると、大抵の場合、第二項が卓越しており、更にマイヤーの関係式を使うと
V^2 \approx \dfrac{2(C_p - C_v)}{\mu C_v^2} (C_p + C_v) h = \dfrac{2(C_p^2-C_v^2)}{C_v^2} \dfrac{h}{\mu}
と書ける。この近似式は、現代の標準的なCJ理論でも導くことができる。冒頭のランダウ・リフシッツの教科書だと式(121.9)

水素2モルと酸素1モルの爆発の場合、生成ガスは、水蒸気2モルのみで、
h = 242(kJ/mol)
μ = 18(g/mol)
を使うことにする。更に、高温の非直線型三原子分子理想気体に対するCp=7R,Cv=6Rを使うと、CJ速度は
V = 3116(m/s)
となる。

デトネーション伝播の基礎
https://www.jstage.jst.go.jp/article/jcombsj/55/174/55_317/_pdf/-char/ja

の表1を見ると、CJ速度は2841(m/s)で、"既燃ガスの実効的な比熱比"は、1.129となっている。観測されるデトネーション速度も、こんなものらしい


上で使用した反応熱は298K,1atmのもので、0K換算の値を使うべきと思う。0Kに於ける生成エンタルピーは、定圧比熱を、0(K)から298(K)まで積分した値を引けば得られる。面倒なので、(振動の影響を無視した常温での)定圧比熱と基準温度の積で見積もることにすると、H2(g)の場合、定圧比熱が3.5Rとすれば、8.672(kJ/mol)くらい引けばいい。O2(g)とH2O(g)に対して、同様に熱補正を考慮すると、0(K)に於ける反応エンタルピーは、上の値から3.1(kJ/mol)ほど引いたものになる。まぁ影響は小さい


もうひとつ、上の近似式は、反応ガス、生成ガスが共に理想気体と仮定されていて、定圧比熱が温度に依らず一定であることは重要。もうちょっとよく考えると、CJ理論のRankine-Hugoniotの式を使う定式化では、波面前後でのエネルギー保存則が、
h_1(T_1) + \dfrac{u_1^2}{2} = h_2(T_2) + \dfrac{u_2^2}{2}
と書かれる。添字の1と2は、反応ガスと生成ガスの物理量を表し、h_i,u_i,T_iは、比エンタルピーと流速、温度。q = h_1(0)-h_2(0)が、0Kに於ける(単位質量当たりの)反応熱(Chapmanの論文の記号に従うと、h/μに相当する量)。で、生成ガスが理想気体と仮定すれば
h_2(T_2) = h_2(0) + \dfrac{\gamma_2}{\gamma_2-1} \dfrac{N k_B T}{\mu}
となる。\gamma_2は生成ガスの比熱比。\mu,Nは、反応前後の質量と、反応後の分子数。

T_2は高温であるが、単原子分子気体の場合を除けば、分子内振動の影響で、比熱は全温度領域で一定でありえない。単原子分子気体は反応性が低く、デトネーションの文脈では多原子分気体が、いつでも重要だろう。多原子分気体では、上のエンタルピーの式を、そのままで使うのは、どう考えても良くない

考えられる処方として
C = \lim_{T_2 \to \infty} \left( h_2(T_2) - h_2(0) - \dfrac{\gamma_2}{\gamma_2-1} \dfrac{N k_B T}{\mu} \right)
を定義して、高温領域で
h_2(T_2) \approx (h_2(0) + C) + \dfrac{\gamma_2}{\gamma_2-1} \dfrac{N k_B T}{\mu}
と近似すれば、q = h_1(0)-h_2(0)-Cと取りなおすだけで、Chapman-Jouguet速度に関する解析的な式を使うことも許されるはず

少し手を抜いて、並進と回転成分については、全領域で高温近似が妥当だとみなし、振動成分からくる比熱のみ、量子統計力学と古典統計力学で、差があると近似する。基準振動モードが一つで、波数νの場合は
C = \dfrac{N k_B T}{\mu} \displaystyle \int_{0}^{\infty}(f(\dfrac{c h \nu}{k_B T}) - 1)dT
f(x)=\dfrac{x^2 e^{x}}{(e^{x}-1)^2}
で、積分を実行すると
 C = -\dfrac{N}{\mu} \dfrac{c h \nu}{2}
となる。基準振動モードが複数ある場合は足すだけ。H2O(g)の場合、1595(/cm),3756(/cm),3657(/cm)を使うと

Cμ = -6.02214076e23 * 2.99792458e8 * 6.62607015e-34 *(1595e2+3756e2+3657e2)*0.5 / 1000.0 = -53.88(kJ/mol)

となる。符号が負なので、CJ速度は、より大きく評価されるようになり、ますます測定値と合わなくなる(ダメ)けど、少なくとも、Chapmanの原論文の意図からすると、このような数値を使うのが正しい気がする。Cの値は熱力学データから決まらない(比熱の温度依存性からフィッティングを行えば、基準振動を知ることができるかもしれないけど、普通、そんなことはしない)せいか、普通、教科書を読んでも、こうは書いてない

まだ、生成ガスの比熱比の値が、考慮に入れてない効果によって、7/6より小さいかもしれないという希望は残ってるけれど


Lasing action and the relative populations of vibrationally excited carbon monoxide produced in pulse-discharged carbon disulfide-oxygen-helium mixtures
https://pubs.acs.org/doi/abs/10.1021/j100639a019
では、CS2-O2-He混合ガス火炎中に生成するCOの振動準位の分布が、点火後の経過時間で、どう変化するか計測している。論文Figure9によると、点火後、数十μsまでは、ボルツマン分布と異なり、特定の振動準位にピークがあり、200μs以上かけて、ボルツマン分布に近付くらしい。デトネーションでも、平衡に達するのに、同じような時間スケールが必要だとすると、例えば、速度3000(m/s)で、200μsの間に、60cm進む。これは、平均自由行程より遥かに長いから、化学反応が一瞬で終わって、短時間で平衡に達するという近似の方が問題かもしれない。


【補足】Chapman-Jouguet理論で必要な熱力学パラメータは、未知化合物でも、量子化学の力によって原理的には計算できる。0(K)に於ける生成エンタルピーは、電子エネルギー+分子内振動の零点エネルギーで計算すればいい。電子エネルギーは、量子化学計算で、ほぼ必ず計算する量で、電子系のハミルトニアンの最低固有値に相当する。分子内振動の零点エネルギーは、分子内振動を調和振動子で近似して最低エネルギーを見る。例えば、H2(g)の場合、基準振動波数ν=4160(/cm)とすれば、
0.5*6.02214086e23*2.99792458e8*6.62607004e-34*4160e2/1000 = 24.88(kJ/mol)
くらいの影響がある。実際上の問題としては、反応熱には、複数の分子の電子エネルギーの差が寄与することになるけど、通常、その大きさは電子エネルギーと比較すると非常に小さい。従って、電子エネルギー自体が1%程度の誤差で計算できても、反応熱の誤差は、割と大きくなり得る 

higher genus AGMは超幾何級数表示を持つか

a,bの算術幾何平均を、M(a,b)として、FがGauss超幾何関数の時

\dfrac{1}{M(1,x)} = F(\dfrac{1}{2} , \dfrac{1}{2},1 ; 1-x^2)

であることは、よく知られてる。


19世紀後半に、Borchardtという人が、種数2のテータ関数(テータ定数)の倍角公式に基づいて、算術幾何平均の一般化を作ったようだ。

[PDF] 算術平均の周辺:Borchardtによる4項算術幾何平均とThomaeの公式
https://www2.tsuda.ac.jp/suukeiken/math/suugakushi/sympo18/18_2shiga.pdf

Borchardtの論文は、オンラインでは見つけられなかった

Über das arithmetisch-geometrische Mittel aus vier Elementen
https://books.google.co.jp/books/about/%C3%9Cber_das_arithmetisch_geometrische_Mitt.html?id=iA0inQAACAAJ&redir_esc=y

通常の算術幾何平均では、楕円曲線y^2=4(x-a_1)(x-a_2)(x-a_3)(a_1 \lt a_2 \lt a_3)に対して、2つの周期
\omega_1 = \dfrac{\pi}{M(\sqrt{a_3-a_1} , \sqrt{a_3-a_2})}
\omega_2 = \dfrac{i \pi}{M(\sqrt{a_3-a_1} , \sqrt{a_2-a_1})}
を計算できるのと同様、Borchardt AGMの場合は、2つの2x2周期行列の行列式が計算できる。BorchardtのAGMは、テータ関数側から作ってるので、超幾何関数に相当するようなものは、調べた限り、明示的に書かれているものはない。時期的には、Borchardtの論文は1876年で、Appelの超幾何級数の論文が1880年、Lauricellaの超幾何級数の論文が1893年(Wikipedia情報)なので、多変数超幾何関数の方が新しい


最近は(?)、他にも、色んな"平均"が構成されていて、これらは、超幾何関数の変換公式を利用している

[PDF] 超幾何関数の関数等式と算術調和平均およびその拡張
https://www.riam.kyushu-u.ac.jp/fluid/meeting/18ME-S5/papers/Article_No_38.pdf

超幾何関数の変換公式と平均反復
http://fe.math.kobe-u.ac.jp/Movies/cm/2009-01-08-matsumoto.pdf


このへんを、何気なく、見てると、とある四項平均が、三変数のLauricella超幾何関数の二乗で書けるという式が書いてあったので、Borchardt AGMも、意外と簡単な形で書けたりするのだろうかと疑問に思った。


頑張れば、超楕円積分をLauricella超幾何関数で書いて解析的に計算できるかもしれないけど、別に重要な問題とも思えないので、そこまでモチベーションがない。Mathematicaなら、超楕円積分を適当な超幾何関数で書くくらい、やってくれないかと期待したけど、無理だったので、数値計算でやってみることにした。

何もなしには出来ないので、いくつか仮定を置く。W(a,b,c,d)をBorchardt AGMとして、
\dfrac{1}{W(1,x_1,x_2,x_3)} = F(1-x_1^2 , 1-x_2^2 , 1-x_3^2)
となる冪級数Fがあり、係数は、有理数で書けるものとする。(x1,x2,x3)=(1,1,1)の時、左辺は1なので、F(0,0,0)=1である。収束半径とかは分からないけど、どうせ原点近くでの振る舞いしか見ないので、気にしなくていいだろう

Borchardt AGMの計算は、平方根を使うので、有理数だけで済ますことはできず、倍精度程度の浮動小数点数では全然精度が足りない。ということで、任意精度演算をやるために、pythonのmpmathを使う。やることは簡単なので、別になんでもいいけど

最初は、3変数の内、2つは0とする。
F(z,0,0) = 1 + c_{100}z + c_{200} z^2 + \cdots
で、zが十分小さければ、高次の項の寄与は、殆ど無視できるはず

例えば、一次の項を決めるのに、以下のようなコードを実行する

from mpmath import *
mp.dps = 500


def M(x,y,z):
   a,b,c,d = 1,x,y,z
   max_iter = 1000
   for n in range(max_iter):
      a2 = (a+b+c+d)/4
      b2 = (sqrt(a*b) + sqrt(c*d))/2
      c2 = (sqrt(a*c) + sqrt(b*d))/2
      d2 = (sqrt(a*d) + sqrt(b*c))/2
      a,b,c,d=a2,b2,c2,d2
   return a


if __name__=="__main__":
  for eps in [1.0e-5 , 1.0e-10, 1.0e-15, 1.0e-20]:
    x1,x2,x3=1-mpf(eps),1,1
    z1,z2,z3=1-x1*x1,1-x2*x2,1-x3*x3
    c = (1/M(x1,x2,x3)-1)/z1
    print(str(c)[:50])

結果は、

0.125001171883984440643884854175753723800930555241
0.125000000011718750000898437927005455667004923082
0.125000000000000117187500000000098949851548046299
0.125000000000000000001171874999999999935735474360

とかなるので、まぁ、0.125=1/8に収束しそう、と分かる。どこに収束するかは勘で判断してるだけので、合ってるかどうか証明とかはない

同様に、z_1^2の係数は、以下のようなコードで決める

if __name__=="__main__":
  for eps in [1.0e-5 , 1.0e-10, 1.0e-15, 1.0e-20,1.0e-30]:
    x1,x2,x3=1-mpf(eps),1,1
    z1,z2,z3=1-x1*x1,1-x2*x2,1-x3*x3
    c=(1/M(x1,x2,x3)- (1+z1/8))/(z1*z1)
    print(str(c)[:50])

結果は

0.058594492194493161915357892508282237042562173272
0.058593750007421875000699310573188476983608793132
0.058593750000000074218750000000075698227920533509
0.058593750000000000000742187499999999959300436760
0.058593750000000000000000000000074218750000000006

で、0.0589375 = 15/256に収束しそう

続ける。

if __name__=="__main__":
  for eps in [1.0e-5 , 1.0e-10, 1.0e-15, 1.0e-20,1.0e-30]:
    x1,x2,x3=1-mpf(eps),1,1
    z1,z2,z3=1-x1*x1,1-x2*x2,1-x3*x3
    c=(1/M(x1,x2,x3)- (1+z1/8+15*z1*z1/256))/(z1*z1*z1)
    print(str(c)[:50])

で、

0.037109910207646803092967250526967109948749607120
0.037109375005352020264233913616668406369705086927
0.037109375000000053520202636718810362650910452159
0.037109375000000000000535202026367187470651540126
0.037109375000000000000000000000053520202636718754

とかになるので、0.037109375=19/512だろう

だんだん辛くなってきたけど

if __name__=="__main__":
  for eps in [1.0e-5 , 1.0e-10, 1.0e-15, 1.0e-20,1.0e-30]:
    x1,x2,x3=1-mpf(eps),1,1
    z1,z2,z3=1-x1*x1,1-x2*x2,1-x3*x3
    c=(1/M(x1,x2,x3)- (1+z1/8+15*z1*z1/256+19*z1*z1*z1/512))/(z1*z1*z1*z1)
    print(str(c)[:50])

0.026760516142735359850255545761188337409031621442
0.026760101322507572174538118990437334607904260568
0.026760101318359416481971740722706058694628176862
0.026760101318359375000414819717407226539753154096
0.026760101318359375000000000000041481971740722659

は、0.026760101318359375=7015/262144らしい。などとやっていくと、2つの変数が0の場合は決まる。

0でない変数が複数ある場合は、面倒くさくなる。とりあえず、二次、三次の係数は、普通にやればいい。z_1^3z_2なんかは、例えば

def F3(z1,z2,z3):
   v1 = (z1+z2+z3)/8
   v2 = 15*(z1*z1+z2*z2+z3*z3)/256 + 3*(z1*z2+z2*z3+z3*z1)/128
   v3 = 19*(z1*z1*z1+z2*z2*z2+z3*z3*z3)/512 + 3*(z1*z1*z2+z2*z2*z1+z1*z1*z3+z2*z2*z3+z3*z3*z1+z3*z3*z2)/256 + 3*z1*z2*z3/512
   return 1+v1+v2+v3

if __name__=="__main__":
  mp.dps = 1000
  for eps in [1.0e-20,1.0e-40,1.0e-80,1.0e-120]:
    for ratio in [0.9,0.0011,1.1e-30,1.1e-60]:
       x1,x2,x3=1-mpf(eps),1-mpf(eps)*mpf(ratio),1
       z1,z2,z3=1-x1*x1,1-x2*x2,1-x3*x3
       diff = (1/M(x1,x2,x3)- F3(z1,z2,z3) - 7015*(z1**4+z2**4+z3**4)/262144)
       if diff < z1*z1*z1*z2 and diff > z1*z1*z2*z2:
          c = diff/(z1*z1*z1*z2)
          print(str(c)[:50])

とやると

0.007620766506805420299495398378483141958519102149
0.007620766506805419922275006967238197476573541568
0.007614135745958588340065696042324398959715200301
0.007620766506805419922275006967238197472801337654
0.007614135742187500000000000000006621551513671874
0.007614135742187500000377108834006569608349931553
0.007620766506805419922275006967238197472801337654
0.007614135742187500000000000000006621551513671874
0.007614135742187500000000000000000000000000000000

なので、0.0076141357421875=499/65536と求まる


そんな感じで計算してみた表が、以下。合ってるのかどうか分からないけど、分母が、2のべきになったので、多分合ってる。以下が計算値のまとめ

i j k c(i,j,k)
1 0 0 1/8
2 0 0 15/256
1 1 0 3/128
3 0 0 19/512
2 1 0 3/256
1 1 1 3/512
4 0 0 7015/262144
3 1 0 499/65536
2 2 0 789/131072
2 1 1 201/65536
5 0 0 43487/2097152
4 1 0 11685/(2^21)
3 2 0 4161/(2^20)
3 1 1 1059/(2^19)
2 2 1 15123/(2^21)
6 0 0 282387/(2^24)
5 1 0 36591/(2^23)
4 2 0 49185/(2^24)
7 0 0 946293/(2^26)
8 0 0 3323953575/(2^38)

意外と単純な係数が出たけど、超幾何関数そのものではないし、超幾何関数の二乗とかいう可能性も、なさそう。3変数のままやると辛いので、2つの変数は0にすると
F(z,0,0) = 1 + z/8 +\dfrac{15}{256} z^2 + \dfrac{19}{512} z^3 + \dfrac{7015}{262144} z^4 + \cdots
なので、形式的冪級数として、平方根を計算ことができるが、
 \sqrt{F(z,0,0)} = 1 + \dfrac{1}{16} z + \dfrac{7}{16^2} z^2 + \dfrac{69}{16^3} z^3 + \dfrac{6267}{2^19} z^4 + \cdots
で、だめっぽい。

なんか、それらしいLauricella超幾何級数のいくつかの積になったりしないだろうかと思って、脳死状態で、コンピュータに計算してもらったけど、候補を見つけることはできなかった。まぁ、あんまり単純な式はないのかもしれない

とりあえず、数値計算で、冪級数の係数を予想するというのは初めてやったけど、意外と、うまくいくらしい


種数3以上の超楕円積分に対して、Borchardt AGMの一般化が、以下で与えられている。
Hyperelliptic integrals and generalizedarithmetic–geometric mean
https://link.springer.com/article/10.1007/s11139-011-9353-7

防御力の物理〜cavity expansion model

弓矢の貫通力
https://m-a-o.hatenablog.com/entry/2019/09/22/185608

に書いたcavity expansion modelは、由来を、ちゃんと調べなかったけど、起源は、1945年の論文

[1] The theory of indentation and hardness tests
https://doi.org/10.1088/0959-5309/57/3/301

にあるらしい。オリジナルのcavity expansion modelは、延性破壊を起こす材料を対象とした、比較的単純なモデルだけど、その後、色々な魔改造の末、土壌や脆性材料への適用や、時間依存性を考慮した計算(dynamic cavity expansion model)なども現れたという経緯っぽい

著者は、R.F.Bishop , R.Hill , N.F.Mottとなっていて、Mottは、Mott散乱公式やMott転移、Mott絶縁体に名前が残っている物理学者Mottと同一人物っぽい。Mottは、自伝を書いていて(絶版になってるっぽいけど、日本語訳も出ている)、戦時中は、軍事研究に携わったことが書かれているけど、この論文について、明示的な言及は別にない

cf)科学に生きる―ネビル・モット自伝
https://www.amazon.co.jp/gp/product/4532062748

R.Hillは、おそらくRodney Hillで、塑性力学の分野では、そこそこ知られた名前のようである。Rodney Hillについては、以下に詳しく書かれている

Rodney Hill. 11 June 1921 — 2 February 2011
https://royalsocietypublishing.org/doi/abs/10.1098/rsbm.2014.0024

"In 1943 he moved to the Armament Research Department at Fort Halstead in Kent, for three years. Here he was involved in, for example, the modelling of armour penetration by projectiles."だそうな。cavity expansion modelも、装甲の弾丸による貫通を調べる過程で生まれたのかもしれない。

R.Hillは、1950年に以下の本を出版している。

The Mathematical Theory Of Plasticity
https://www.amazon.co.jp/dp/0198503679

日本語訳もあったが、絶版らしい

塑性学 (1954年)
https://www.amazon.co.jp/dp/B000JB7J1I

R.F.Bishopという人のことは、よく分からない。


論文[1]のタイトルにある"indentation and hardness tests"は、押し込み硬さ試験(Indentation hardness tests)のことではあるんだろうけど、通常の押し込み硬さ試験より荷重を大きくしていって、押し込み深さが非常に深くなるような、穿孔に近い状況を調べているっぽい。

いわゆる押し込み硬さは、Brinell硬度とかVickers硬度とかで、くぼみができるまで圧子を押し付けて測定する。圧子の大きさと比べて、そんなに深くは押し込まないと思う。押し込み硬さの歴史はそんなに古くなくて、1900年に、スウェーデンのBrinellが導入したBrinell硬さが、現在使われてる中で、一番古い基準らしい。何故か有名なモース硬度とかは、もっと古くからあるけど、あんまり定量的ではない。


論文[1]で扱われてる対象は、穿孔機(先端形状は、円錐で、先端の角度は複数ある)を、非常に深く押し込んだ時の(単位面積当たりの)抗力で、論文[1]のFigure3とFigure4で与えられた実測結果から、深さが十分であれば、それはある一定値に収束すると見なせそうではある。この抗力の大きさは、論文で使っている材質(hardened copperとannealed copperと書いてある)では、降伏応力の4〜5倍に達する。

Figure3とFigure4の結果では、極限抗力は、多少幅があるけど、おおまかには、100tons/sq.inch.程度っぽい。単位がわかりづらいけど、100[tons/sq.inch]=1.0e3*100*9.8/(25.4*25.4)=1520[N/mm^2] = 1520[MJ/m^3]程度。これは、銅なので、前回の銃弾によるスチールの貫通抗力よりは、小さくなっている。

前回と同じく、比較として、単位体積の銅を融解するのに必要な熱エネルギーを見ておくと、比熱(※)は、常温以上の金属では、大体、3R≒25(J/mol/K)で、融点が1000度強らしいので、融点まで熱するのに、ざっくり25(kJ/mol)で、融解熱が13.26(kJ/mol)とある。モル質量が63.55(g/mol)で、密度が8960(kg/m^3)から、8960*(25+13.26)*1.0e-3/63.55e-3≒5400(MJ/m^3)程度だろう。1520/5400=0.28…で、割合としては、球状弾でスチールを貫いた時より、多少小さいけど、大きな差はない

※)面倒なので、Dulong-Petit則が使える定積比熱を用いてるけど、定圧比熱を使うほうが適切な気はする。常温の固体では、定積比熱と定圧比熱は大差なくても、高温では、それなりの差が出る可能性がある。けど、オーダーを見るだけなので気にしない

穿孔の手段として、現代だとレーザーなんかもあると思うけど、エネルギー的には、レーザーより、ずっと効率がよさそう。


この抗力の大きさを見積もるための論文[1]の方針は、以下のように書かれている。

Denoting ths maximum load by p_0 A, where A is the area of the cross-section of the punch, it is found that p_0, for a lubricated punch is about twice the hardness, or five times the yield stress, of the work-hardened material. A theoretical method is given for calculating p_0, as follows : the pressures p_c, and p_s, required to enlarge a cylindrical and a spherical hole in a material showing any kind of strain hardening can be calculated. It is plausible to assume that p_0, should be between p_c, and p_s, and since p_s, is only slightly greater than p_c, an approximate theoretical estimate of p_0, is obtained. This is in good agreement with experiment.

抗力の大きさは、hardnessの二倍とか書いてある。このhardnessが何を指してるのか明言されてないけど、何らかの押し込み硬さだろうと思う

cavity expansion modelで計算する力は、物体内に作られた円筒状あるいは球状の穴の径を拡大する時に受ける力であって、穿孔機がプレートを貫く時の抗力とかとは(少なくとも、見かけ上は)違うものである。plausibleとか言ってるので、Bishop-Hill-Mottの主張も、それほど強い根拠があってのものでもなさそう

なので、安易に同一視しないため、論文では使われてない用語であるけど、p_cをcylindrical cavity expansion pressure、p_sをspherical cavity expansion pressureと呼んでおく。

また、摩擦については本文内に言及はあるものの見積もれないと書いてある。


論文[1]の計算では、線形硬化弾塑性体(降伏応力以下では、線形弾性体で、降伏応力以上でも、ひずみと応力に一次の関係式が成立する。降伏応力以上で、ひずみの応力に対する比例係数を塑性係数と呼ぶ)を仮定した計算を実行しているけど、他のモデルを利用することもできる。ただ、塑性係数を0にしても比較的影響は小さいようなので、真面目に扱ってみても、大した意味はないかもしれない。

計算の主要な部分、特に、cylindrical cavityに関する部分は、1931年出版の

Plasticity: A Mechanics Of The Plastic State Of Matter
https://archive.org/details/in.ernet.dli.2015.166225/page/n213

という本に既に書かれている。この本は、塑性力学に関する、(少なくとも、英語で書かれたものとしては)世界最初の教科書らしい。The thick-walled tubed under internal pressureというタイトルの章がある。論文[1]独自の部分には、妥当性の疑わしい議論が含まれると思う。

論文は、古い時代のもあって、よく分からない計算が含まれているので、現代の適当な教科書を見たほうがいい。塑性力学の教科書は沢山はないけれど、日本語で読めるものから探すと、必要な計算の(全てではないが)殆どは、例えば、

基礎塑性力学 (実用理工学入門講座)
https://www.amazon.co.jp/dp/481730149X
http://www.nissinpb.co.jp/jitsuyorikogaku-kikai.htm#%E5%9F%BA%E7%A4%8E%E5%A1%91%E6%80%A7%E5%8A%9B%E5%AD%A6

の6.4節、6.5節とかに書いてある。


spherical cavityの場合、球の中心からの距離をsとして、cavity半径がaであるとし、 a \leq s \leq cを塑性領域、c \leq sを弾性領域とする。

Poisson比を\mu < 0.5として、弾性領域c \leq sでは、問題なく、consistentな解
u(s) = \dfrac{(1+\mu)Y}{3E}\left( \dfrac{c^3}{s^2} \right)
\sigma_r(s) = -\dfrac{2Y}{3}\left( \dfrac{c}{s} \right)^3
\sigma_{\theta}(s) = \dfrac{Y}{3}\left( \dfrac{c}{s} \right)^3
\epsilon_{r}(s) = -\dfrac{2(1+\mu)Y}{3E}\left( \dfrac{c}{s} \right)^3
\epsilon_{\theta}(s) = \dfrac{(1+\mu)Y}{3E}\left( \dfrac{c}{s} \right)^3
が見つかる。uはradial displacementで、Eはヤング率、Yが降伏応力etc.。\sigma_{r},\sigma_{\theta},\epsilon_r,\epsilon_{\theta}は、radial stress, tangentical stress, radial strain, tangentical strain。

現在の大抵の教科書は、球全体の半径bを有限としているが、式が複雑になるので、単に、b→∞の極限を取った。弾性領域の計算は、古典的な弾性論通りなので、特に問題はない

これらの関数が満たすべき条件は、いくつかある(Hookeの法則とか平衡方程式とか)けど、全部書くと面倒なので省略。降伏条件は、Misesの降伏条件でもTrescaの降伏条件でも同じで
\sigma_{\theta}(c) - \sigma_{r}(c) = Y
を採用している。

弾性領域と塑性領域の境界に於いて
u(c) = \dfrac{(1+\mu)Y}{3E}c
であるけど、金属などでは、ヤング率は降伏応力よりずっと大きいので、u(c)/cは小さい。

塑性領域の体積が空洞形成前後で不変だと仮定すると
\dfrac{4\pi}{3}(c^3 - a^3) = \dfrac{4\pi}{3} (c-u(c))^3
で、近似的に
u(c) \approx \dfrac{a^3}{3c^2}
となって、
\left( \dfrac{c}{a} \right)^3 = \dfrac{E}{(1+\mu)Y}
が導ける。

塑性領域a \leq s \leq cで平衡方程式
 \dfrac{d \sigma_r}{ds} = \dfrac{2(\sigma_{\theta} - \sigma_r)}{s}
が成立して、塑性係数が0なら、塑性領域に於いて
\sigma_{\theta}(s) - \sigma_{r}(s) = Y
となって、平衡方程式を積分して、応力が弾塑性領域境界で連続だと仮定すると
\sigma_r(s) = \sigma_{r}(c) - 2Y \log(\dfrac{c}{s}) = -\dfrac{2Y}{3} - 2Y \log(\dfrac{c}{s})
と計算できる。

知りたいcavity expansion pressureは、cavityと塑性領域の境界位置s=aに於ける-\sigma_{r}(a)の値なので、塑性領域に於ける歪みの詳細が分からなくても、欲しい答え
P = \dfrac{2Y}{3} + 2Y \log(\dfrac{c}{a})
が得られる(\dfrac{c}{a}は既に決定した)。

論文[1]では、塑性係数が0とは仮定していないので、かなり乱暴に、塑性領域での歪みを決め打ちしているけど、弾塑性領域境界での連続性も成立してないし、空洞・塑性領域境界では発散しているので、あまり鵜呑みにできない。しかしまぁ、計算された結果を信じるなら、
 P = \dfrac{2Y}{3} + 2Y \log(\dfrac{c}{a}) +\dfrac{2\pi^2}{27}A
となる。但し、Aは塑性係数。



前回計算した"破断エネルギー(密度)"と比較されるべき量は、(cavity expansion pressureではなく)弾性変形領域に蓄えられている歪みエネルギーと塑性仕事の和を、cavityの体積で割った量になるはず。cavity形成に投入されたエネルギー(銃弾による貫通の場合は、着弾時運動エネルギー)は、弾性域のひずみエネルギー+塑性仕事+その他散逸に等しい。その他散逸というのは、表面摩擦熱などだけど、今は無視する

塑性仕事の計算には、塑性歪み増分を知る必要があるけれど、塑性係数が0の時は、適当な仮定を置くことで、塑性域の歪み(とradial displacement)をモデル化する方法がある(論文[1]には書いてないけど、教科書"基礎塑性力学"には書いてある)。

塑性領域では
\epsilon_r = \dfrac{1}{E} \left(\sigma_r - 2\nu \sigma_{\theta}\right) - \dfrac{2}{3}(\sigma_{\theta} - \sigma_{r})\phi
\epsilon_{\theta} = \dfrac{1}{E} \left( (1-\nu)\sigma_{\theta} -\nu\sigma_r \right) + \dfrac{1}{3}(\sigma_{\theta} - \sigma_{r})\phi
の形の歪みを仮定しても矛盾は、ない。\phiは関数で、それぞれ、歪みの第一項は、Hookeの法則と同じ形なので、弾性歪み成分と呼ばれ、第二項は塑性歪み成分と呼ばれる。

特に、この形の歪みは、
\epsilon_r + \epsilon_{\theta} = \dfrac{1-2\nu}{E} (\sigma_r + 2\sigma_{\theta})
を満たす。radial displacementは
\epsilon_r(s) = \dfrac{du}{ds} , \epsilon_{\theta}(s) = \dfrac{u}{s}
を満たすように決める。既に、応力は計算したので、これらの条件から、歪みを決められる。具体的には
\phi(s) = \dfrac{3(1-\nu)}{E} \left( \dfrac{c^3}{s^3} - 1 \right)
となり、歪みの塑性成分は、
\epsilon_r^{(p)}(s) = \dfrac{-2Y}{3} \phi(s)
\epsilon_{\theta}^{(p)}(s) = \dfrac{Y}{3} \phi(s)
のように決まる。

計算したい量は、半径aのcavityを準静的に形成する塑性仕事W(a)で、塑性仕事増分は、
\dfrac{dW}{da} = \displaystyle \int_{a}^{c} 4\pi s^2 \left( \sigma_r(s) \dfrac{d\epsilon_{r}^{(p)}}{da} + 2\sigma_{\theta}(s)\dfrac{d\epsilon_{\theta}^{(p)}}{da} \right) ds
と求まるだろう。
c^3 = \dfrac{E}{(1+\nu)Y} a^3
に注意すると、(当然、W(0)=0なので)私の計算では
W(a) = \dfrac{4 \pi a^3}{3} \dfrac{6(1-\nu)Y}{1+\nu} \log(\dfrac{c}{a})
となった。 ポアソン比が、0.5なら、
W(a) = \dfrac{4 \pi a^3}{3} (2 Y \log(\dfrac{c}{a}))
となって、弾性変形による歪みエネルギー(この計算は、弾性論の範疇なので省略)をcavity体積で割った量が、丁度、2Y/3なので、足すとspherical cavity expansion pressureに等しくなる。

金属とかだと、ポアソン比は、0.3〜0.4くらいだったりするので、塑性仕事+弾性歪みエネルギー/空洞体積の値は、spherical cavity expansion pressureより、もう少し大きい値になる。文献では見当たらない結果なので、見なかったことにしよう(まぁオーダーは合ってるし...

実際のとこ、c/aは定数としているけど、弾丸がプレートを貫通する前とかは、a=c=0で、これらの仮定の妥当性が疑われる。また、c/aが上で決めた値と違っても、定数であれば、W(a)は空洞体積に比例するので、c/aが素早く定数値に近付くのであれば、結果は少し変わるけど、計算自体に変更は要らない。例えば
\left( \dfrac{c}{a} \right)^3 = \dfrac{E}{3(1-\nu)Y}
であれば、
W(a) = \dfrac{4 \pi a^3}{3} (2 Y \log(\dfrac{c}{a}))
が(ポアソン比に依らず)常に成り立つ。



加えられた塑性仕事の大部分は、熱として散逸する(らしい)。もし、空洞形成が(熱伝導より)素早く起きるなら、塑性仕事を塑性域の体積で割って、更に比熱で割ると、塑性域の温度が何度くらい上昇するか見積もれそうに思える。金属の場合、塑性域の体積は空洞体積よりずっと大きいので、簡単のため、塑性域体積+空洞体積で割る。どうせ概算なので、ポアソン比は0.5としておく。すると
W(a)/(\dfrac{4 \pi c^3}{3}) = \dfrac{3Y^2}{E} \log(\dfrac{c}{a})
で、E=200(GPa)とY=500(MPa)を採用すると、19(MJ/m^3)くらいの塑性発熱があるだろう。鉄の場合は、比熱が3.5[MJ/m^3/K]程度なので、ざっくり5度くらいの温度上昇がある計算になる。

ヤング率も降伏応力も、温度によって変化する(なので、鉄を溶かせなかった時代でも、熱した鉄を叩いて加工できた)が、5度くらいなら影響を無視できるだろう

現在では、高速な塑性変形で材料加工する際の温度変化を計測して、塑性変形量を計測する技術というものが存在するっぽい

高速変形試験(4)~高速変形時の加工発熱分布の計測技術~
https://www.jfe-tec.co.jp/jfetec-news/27/5p.html

cylindrical cavityでもspherical cavityでも、cavity周りに塑性領域が存在すると考える点では同じだが、c/a(塑性領域半径/空洞半径)の値は、ヤング率、降伏応力、ポアソン比が同じ場合でも結構違う。例えば、上のE=200(GPa)、Y=500(MPa)、ν=0.5だと、cylindrical cavityでは、c/a≒15.2で、spherical cavityでは、c/a≒6.4なので、弾丸が装甲を貫いた時に、穴の周囲に、どれくらいの大きさの塑性領域が形成されるか分かれば、どちらのcavity modelに近いか判定する基準になるかもしれない



論文[1]では、(cold-workedとあるので、冷間加工した)銅に対して、
降伏応力Y=285(MPa)
塑性係数A=100(MPa)
Young率E=130(GPa)
Poisson比ν=0.34
という値を使って(応力などの単位は、tons/sq.inchなので、換算した)いて、spherical cavity expansion pressureはP=1375(N/mm^2)程度となる。穿孔機での実測値と比較するとオーダーは合ってる

論文[1]のFigure7で、軟鋼に対する応力歪み曲線が掲載されていて、鋼のヤング率は、大体200GPaくらい、ポアソン比は0,3くらいで、それ以外の値は、私が目視で読み取った限り、
降伏応力Y=275(MPa)
塑性係数A=350(MPa)
Young率E=200(GPa)
Poisson比ν=0.3
で、spherical cavity expansion pressure=1600(MPa)=1600(N/mm^2)となる。論文著者たちの計算だと、120[tons/sq.inch] ≒1960[N/mm^2]程度。計算に使った値が書かれてないけど、まぁ、こんなものだろう



前回のマスケット銃に対するGrazテストでは、ブリネル硬さ290のcold-worked mild steelを使ったとか書いてあった

出典:Material Culture and Military History: Test-Firing Early Modern Small Arms
https://journals.lib.unb.ca/index.php/MCR/article/view/17669/22312

引張強度(近似値)換算で、965(MPa)程度と思われるが、降伏応力の値は分からない。引張強度が400(N/mm^2)のSS400という鋼材だと、降伏応力が235~245(N/mm^2)と規格で決まってるので、ざっくり、引張強度の60%程度である。965(MPa)の60%で、Y=580(MPa)として、A,E,νは上記の値を流用すると、spherical cavity expansion pressure P=2800(MPa)と計算され、Grazテストの破断エネルギーとオーダーは合ってる(勿論、偶然かもしれない)

同じパラメータを使って、cylindrical cavity expansion pressureだと、P=2350(MPa)程度。

表を再掲。前回書き忘れた(元論文には書くまでもないのか何も言及されてない)けど、銃弾の素材は、(密度を計算すると)鉛っぽい。

銃器 試験 運動エネルギー "破断エネルギー" 着弾速度(m/s)
"Doppelhaken" G 284 PS100 1779.6J 3138(MJ/m^3) 305
"Doppelhaken" G 358 PS100 2992.7J 2335(MJ/m^3) 349
Matchlock LG 1514 PS30 1241.7J 3467(MJ/m^3) 378
Matchlock LG 1514 PS100 605.7J 3382(MJ/m^3) 264
Wheellock RG 33 PS30 2333.2J 3347(MJ/m^3) 394
Wheellock RG 33 PS100 1238.0J 2664(MJ/m^3) 287
Wheellock RG 117 PS30 660.2J 2779(MJ/m^3) 349
Wheellock RG 117 PS100 307.0J 2539(MJ/m^3) 238
Wheellock RG 272 PS30 1874.9J 3898(MJ/m^3) 342
Wheellock pistol RP 2895 PS30 602.4J 2754(MJ/m^3) 355
Flintlock STG 1287 PS30 2269.8J 2560(MJ/m^3) 406
Flintlock STG 1287 PS100 1166.1J 2630(MJ/m^3) 291
Flintlock STG 1288 PS30 2032.8J 3131(MJ/m^3) 390
Flintlock STG 1288 PS100 1055.3J 2438(MJ/m^3) 281
Flintlock STG 1316 PS30 2458.3J 3368(MJ/m^3) 391
Flintlock STG 1316 PS100 1324.5J 2722(MJ/m^3) 287
Flintlock STG 1318 PS30 2806.5J 2917(MJ/m^3) 426
Flintlock STG 1318 PS100 1438.6J 2991(MJ/m^3) 305
Flintlock pistol STP 1128 PS30 753.8J 2633(MJ/m^3) 323
アサルト・ライフル58 PS100 2801.5J 4824(MJ/m^3) 770
アサルト・ライフル77 PS100 1375.0J 5987(MJ/m^3) 874
Glock80ピストル PS30 467.9J 3444(MJ/m^3) 342

Grazテストでは、スプルース(木材の一つ)をターゲットにしたデータもあるので、同様の計算をしてみる

銃器 試験 運動エネルギー "破断エネルギー" 着弾速度(m/s)
"Doppelhaken" G 284 PW100 1779.6J 41.0(MJ/m^3) 305
"Doppelhaken" G 358 PW100 2992.7J 49.4(MJ/m^3) 349
Matchlock LG 1514 PW30 1241.7J 53.0(MJ/m^3) 378
Matchlock LG 1514 PW100 605.7J 40.6(MJ/m^3) 264
Wheellock RG 33 PW30 2333.2J 52.9(MJ/m^3) 394
Wheellock RG 33 PW100 1238.0J 66.6(MJ/m^3) 287
Wheellock RG 117 PW30 660.2J 42.1(MJ/m^3) 349
Wheellock RG 117 PW100 307.0J 30.8(MJ/m^3) 238
Wheellock RG 272 PW30 1874.9J 46.4(MJ/m^3) 342
Wheellock RG 272 PW100 1083.6J 43.7(MJ/m^3) 260
Wheellock pistol RP 2895 PW30 602.4J 45.5(MJ/m^3) 355
Flintlock STG 1287 PW30 2269.8J 89.0(MJ/m^3) 406
Flintlock STG 1287 PW100 1166.1J 63.4(MJ/m^3) 291
Flintlock STG 1288 PW100 1055.3J 61.0(MJ/m^3) 281
Flintlock STG 1316 PW30 2458.3J 51.8(MJ/m^3) 391
Flintlock STG 1316 PW100 1324.5J 37.0(MJ/m^3) 287
Flintlock STG 1318 PW30 2806.5J 63.8(MJ/m^3) 426
Flintlock STG 1318 PW100 1438.6J 52.5(MJ/m^3) 305
Flintlock pistol STP 1128 PW30 753.8J 46.2(MJ/m^3) 323
アサルト・ライフル58 PW100 2801.5J 120.0(MJ/m^3) 770
アサルト・ライフル77 PW100 1375.0J 187.8(MJ/m^3) 874
Glock80ピストル PW30 467.9J 54.7(MJ/m^3) 342

材料力学的パラメータは
Wood, Panel and Structural Timber Products - Mechanical Properties
https://www.engineeringtoolbox.com/timber-mechanical-properties-d_1789.html
に、
ヤング率E〜10(GPa)
降伏応力Y〜70(MPa)
とあるけど、降伏応力が破断エネルギーより大きい。そんなわけ無いと思うのだけど、木材の強度は、方向依存性があって、繊維に直交する方の降伏応力は、もっとずっと弱いようだ。cavity expansion modelは等方性を仮定したモデルなので、そのままでは適用できない


論文[1]の実験だと、抗力がサチるには、(穿孔機の直径と比較して)ある程度の深さまで、押しこむ必要があるっぽい。一方、Grazテストでは、直径1cm程度の弾丸で、数mmのプレートを貫通するという状況がメインなので、最大抗力が発揮されるのか謎である

論文[1]より少し前に、BetheやG.I.Taylor、プレートの厚みが比較的薄い状況を念頭に置いてモデル化を試みている。

まず、Betheは1941年、

Attempt of a Theory of Armor Penetration
https://books.google.co.jp/books/about/Attempt_of_a_Theory_of_Armor_Penetration.html?id=AYzkwgEACAAJ&redir_esc=y

という報告を書いてるけど、オンラインでは読めないっぽい。これに対して、G.I.Taylorは、Notes on H.A. Bethe's "Theory of armor penetration"という文書を書いていて、ネット上で読むことができる。G.I.Taylorの結果は、Betheのものとは少し違うようだけど、極端な差はないらしい

Notes on H.A. Bethe's "Theory of armor penetration", I. Static penetration,
https://apps.dtic.mil/dtic/tr/fulltext/u2/1009963.pdf

Notes on H.A. Bethe's "Theory of armor penetration", II. Enlargement of a hole in a flat sheet at high speed
https://apps.dtic.mil/dtic/tr/fulltext/u2/1009965.pdf

この2つの報告より、戦後に書かれた以下の論文を読む方がいいかもしれない。

The formation and enlargement of a circular hole in a thin plastic shee
https://doi.org/10.1093/qjmam/1.1.103

cavity expansion modelでは、空洞の周囲に、塑性変形領域と弾性変形領域の二層ができると考えていた。G.I.Taylorは、プレートの厚みが薄いことにより、空洞の周囲では、厚みが増大するとし、厚みの増大が無視できない塑性変形領域、厚みの増大が無視できる塑性変形領域、弾性変形領域の三層で考えている(Betheも、同様の仮定を採用したようだ)。

一番内側の厚みの増大が無視できない塑性変形領域では、軸中心からの距離sに於ける厚さh(s)を、新しい変数として導入している。プレートは薄く、完全に貫通している状況を考えているので、計算するべきは、抗力ではなく、穴を形成するのに必要なエネルギーということになる。

G.I.Taylorの議論の詳細は見てないけど、結論は、"the work done in expanding to a given radius"が式(81)に書いてあるけど、単位体積あたりの破断エネルギー/降伏応力は、1.33になるらしい。これだと、Grazテストの結果とは合わなそうに思える。モデル化自体に問題はない気がするので、塑性仕事の計算を間違えてるんじゃないかと思うけど、面倒くさいので、詳細は確認してない



ちなみに、G.I.Taylorの1946年の講義録(?)

THE TESTING OF MATERIALS AT HIGH RATES OF LOADING
https://www.icevirtuallibrary.com/doi/pdf/10.1680/ijoti.1946.13699

には、以下のように書いてある。

Lead bullets weighing from 1.8 gram to 2.7 grams were fired from B 0.22-inch rifle at the mid-point of the anvil. We found that the lead behaved like a fluid, spreading out in a thin film over the plane surface of the anvil. me measured the force exerted by the bullet and found it to be, in fact, just what a fluid jet of diameter 0.22 inches and the density of lead would exert.

I have already mentioned that a lead bullet fired at a speed of about 1,OOO feet per second produces a thin spray of lead on striking a steel target, just as a jet of water would. To understand in a quantitative way what happens when impact occurs at such speeds would, in general, be beyond the power of our present analytical methods, even if we knew beforehand the stress-strain characteristics of the material.



もうひとつ、上の話とは直接関係ないけど、800m/s〜2500m/sの高速度領域ではあるが、タングステンカーバイド製の非常に硬い球状の弾丸による貫通深度の実測値が以下の論文に記載されているのを見つけた。

[2]Penetration of HSLA-100 steel with tungsten carbide spheres at striking velocities between 0.8 and 2.5 km/s
https://www.sciencedirect.com/science/article/pii/S0734743X03000800

ターゲットのHSLAプレートは、Table3によれば、
ヤング率E=197(GPa)
降伏応力Y=730(MPa)
Poisson比ν=0.29
らしい。流石に塑性係数は記載されてないので、0としておくと、spherical cavity expansion pressureはP〜3100(MPa)である。

この論文で使用している弾丸は、直径6.4mmの球状弾丸で、タングステンカーバイドの密度は、14.9(g/cm^3)とTABLE3に記載されているので、弾丸質量は、2.045(g)程度

実際の実測値は、Table2を見ると、まず、速度が1000m/s程度までは、形成されたクレーター直径は、弾丸直径と殆ど変わらないようであるが、それ以上では、クレーター直径の増大が起きるらしい。高速になるにつれて、貫通深度は伸び悩むのに、クレーター直径は、ますます大きくなる。

データに含まれるクレーター体積のデータを使って、衝突時の運動エネルギー/クレーター体積によって"破断エネルギー密度"を計算してみる。データは、12個の異なる速度のものがあるけど、そのうち、8個を気分で選んで計算した結果が以下

Shot 速度(m/s)  体積(ml) 破断エネルギー 深さ(mm)
1 830 0.11 6404(MJ/m^3) 4.57
3 980 0.16 6138(MJ/m^3) 5.46
4 1270 0.27 6872(MJ/m^3) 7.09
6 1500 0.35 6574(MJ/m^3) 8.51
8 1910 0.5 7461(MJ/m^3) 8.53
9 2150 0.73 6475(MJ/m^3) 8.41
11 2460 0.93 6654(MJ/m^2) 9.50
12 2550 1.0 6649(MJ/m^3) 9.60

破断エネルギー密度のばらつきは、Grazテストの計算結果よりも小さく、ほぼ一定っぽい。spherical cavity expansion pressureと比較すると、まぁ、二倍くらい?

これくらいずれるとなると、Grazテストの時とは、大きく異なることが起きてる可能性があるが、弾丸材質が違うせいということもありえる(プレートの材質も微妙に違うけど、こっちの差は小さそうに思える)

[3]Impact Phenomena at High Speeds
https://aip.scitation.org/doi/10.1063/1.1722215

[4]Crater Formation in Metallic Targets
https://aip.scitation.org/doi/abs/10.1063/1.1723437

などの古い論文を読むと、様々な材質の球を1000m/sから5000m/s(もしくは750m/sから2250m/s)まで加速して、同じターゲットに衝突させると、形成されたクレーター体積と弾丸の運動エネルギーの比は、弾丸の材質が同じなら一定値を取るということが書いてある。どちらの論文も、単位が、m^3/jouleになっているが、逆数を取ると

弾丸材質 [4]TABLE I [3]TABLE II
亜鉛 1700MJ/m^3 5000MJ/m^3
960MJ/m^3
1700MJ/m^3
3600MJ/m^3
Steel 9000MJ/m^3
200MJ/m^3 610MJ/m^3
Brass 1000MJ/m^3
マグネシウム 1390MJ/m^3
1250MJ/m^3
アルミニウム 1300MJ/m^3
Wax

という感じ。かぶってる材質が殆どないが、亜鉛、鉛では、比率が大体、1:3になっていて、この違いは、プレート側の材質に依存するものだろう。鉄とSteelは一応、区別した

鉛と鉄、あるいは、鉛とsteelでは、10倍以上違いが出てるので、論文[2]の結果とGrazテストの結果の差は、弾丸速度の差によるものかなぁという気がする


高速衝突で見られる破壊モードは、いくつかあるとか書いてあったりする(プラグ破壊とかスポール破壊とか...)けど、意味のある分類なのか、よく分からない

cf)スポール破壊条件の研究
https://www.jstage.jst.go.jp/article/jsms1963/41/468/41_468_1396/_article/-char/ja/

The interrelation of failure modes observed in the penetration of metallic targets
https://www.sciencedirect.com/science/article/pii/0734743X84900010

以下には、アルミニウム合金プレートに、900m/sから1500m/sの弾丸を打ち込んだ時、衝突速度の増大に伴って、破壊モードが、プラグ破壊からdiscingに変化したとか書いてある

Thick AA7020-T651 plates under ballistic impact of fragment-simulating projectiles
https://doi.org/10.1016/j.ijimpeng.2015.08.001

(論文[2][3]を信じるなら)ある程度の高速度領域に於いて、装甲と弾丸材質が同じなら、着弾速度に依存しない破断エネルギー密度が出るのは、不思議なことではあるが、破壊モードの違いにも関わらず、破断エネルギー密度が一定なのだろうかとか、色々疑問は沸くけど、細かいことは、よく分からない。

これらの値を(解析的に)計算できるモデルは、ないようである。また、運動エネルギーが同じなら、柔らかい弾丸の方が、大きなクレーターを形成できるようだけど、これもなんだか不思議な気がする

低次視覚野の計算モデル

今まで、(低次)視覚神経系のモデルについて、見聞きした雑多な話題として、

(1)カブトガニの眼で側方抑制がどうたら...(1950年代後半〜1960年代初頭、Hartline & Ratliff)
(2)ネコの一次視覚野(V1野)に、特定方向の線分に反応する細胞があるとか何とか(1959,Hubel & Wiesel)→1981年のノーベル医学・生理学賞
(3)V1野の方位選択性の機能は、Gaborフィルターでモデル化できる(Jones & Palmer 1987)
(4)infomax原理がどうとか

くらいがあったけど、それぞれの話の繋がりがどうなってるのか、よく分かってなくて、理解しようと思った。


(1)は、生理学とか生物学の結果で、画像処理的には、エッジ強調/画像鮮鋭化に対応する処理が行われているとされる。論文としては
Inhibitory interaction of receptor units in the eye of Limulus.
https://www.ncbi.nlm.nih.gov/pubmed/13398569

など。元々の話は、カブトガニを使ったものだったけど、マッハバンドとかハーマン格子のような古くから知られる錯視を説明でき、人間の網膜にも存在するよう。いずれにせよ、基本的には、網膜で処理されるレベルの話のよう。

マッハバンドは、Wikipediaには「微妙に濃淡の異なるグレーの領域が接触している場合に、暗い方の領域の境界付近はより暗く、明るい方の領域の境界付近はより明るく強調されて見える、錯視の一種である」と説明されている。実際の処理は、一種のhigh-passフィルタでモデル化できると考えられ、このようなhigh-passフィルタは、比較的単純な神経回路で作ることができるということが、基本的には(1)で解明されたらしい。



網膜にある視細胞で受けた情報は、次に、LGN(外側膝状体)という領域に送られる。視細胞には、rod cellとcone cellの二種類があるけど、基本的には、視細胞→双極細胞→網膜神経節細胞(RGC)→(視神経)→LGN→一次視覚野(V1)と伝わっていくらしい。rod cellは、明暗にのみ反応し、cone cellは色情報を担当するとされ、前者は、片目で一億個くらいあるのに対して、後者は500〜700万個くらいらしい。

cone cellも一様に分布してるわけじゃなく、外縁部は低密度らしいけど、平均的には10000x10000くらいの解像度があるということで、8Kモニターより上。cone cellは、大多数の人類が、三色認識するので、一色あたり200万くらいと考えると、色情報についてはフルHD程度の解像度(?)。

視細胞〜網膜神経節細胞には、水平細胞とアマクリン細胞というのもいて、沢山種類があるけど、機能の全容が解明されているわけではないらしい

視覚系の神経細胞は、それぞれが、視野の中の一部の領域からの入力にのみ反応するようになっていて、この領域を受容野と呼んでいる。画像処理で、Gaussian blurとかのフィルター処理をかける際には、サイズが3x3とか5x5の領域の情報から、一ピクセルの色情報を計算するのと同じような話。一般に、受容野は、低次視覚野の神経細胞では狭く、高次視覚野では、広くなっていくらしい。

RGC細胞数が約100万個とか書かれてたりするので、rod cellの数との比は1:100くらい。単純に考えると、10x10くらいのrod cellからの入力が一つにまとめられている計算になる(?)。この比率は、網膜位置によって、かなり異なるらしいので、目安程度にしかならない。

低次視覚野は、局所的な幾何学的特徴に反応するのに対して、高次視覚野では、顔や、特定の物体に反応する細胞とかがあるとされている。


で、神経細胞は、受容野に入ってきた入力をまとめて、ある大きさの応答を起こす。RGCやLGNの応答は、DoGフィルタとしてモデル化でき、一次視覚野の単純細胞の応答は、Gaborフィルタでモデル化できるとされている。というのが、(2)や(3)の話。DoGフィルタによるモデル化は、誰が言い出したのか分からないけど、よく書いてある

と言われても、フ〜ンという感想しかなかったので、実際に、DoGフィルタとかGaborフィルタを通した結果が、どんな感じになるか見てみた。

DoGフィルタやGaborフィルタには、いくつかのパラメータが存在し、動物や人間の神経系のモデルとして、どのようなパラメータを選ぶべきか、通常、何の説明もない。
Analysis of multidimensional difference-of-Gaussians filters in terms of directly observable parameters
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3789628/

Modeling lateral geniculate nucleus response with contrast gain control. Part 1: Formulation
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3918962/

などでは、DoGカーネルを3つのパラメータで規定される
f(x,y) = \dfrac{1}{2 \pi \sigma_C^2} \exp(-\dfrac{x^2+y^2}{2\sigma_C^2}) - \dfrac{\beta}{2 \pi \sigma_S^2} \exp(-\dfrac{x^2+y^2}{2\sigma_S^2})
の形に取っている。βをbalance parameterと呼び、無次元量である。他の2つのパラメータは、長さの次元を持つので、面倒だけど、例えば、2つの比は無次元量である。積分すると、1-βになる。

CV業界では、β=1として、残り2つのパラメータ比は、1.6にすることが多いっぽい。上の論文は、観測可能な量から各パラメータを見積もろうというノリらしいけど、実際の値は与えられていない。後者の論文の§2に於ける記述によれば、LGNにおけるβの値は、生物種や細胞種ごとに、0.95から0.52まで、色んな見積もりがあるらしい。

DoGフィルタは、フィルタを通しただけでは、あまり意味はなくて、ゼロ交差の検出をしなければならない。このへんは、工夫の余地があるのかもしれないけど、何も考えない実装になっている。

あと、適当に検索すると、gray画像のピクセル値をuint8のまま、2つのGaussianフィルタに通して(cv2.filter2Dを使うと、適当に近い整数値を選んで、出力もuint8のままにしてくれる)、そのまま引くという実装をしている人がいた。大体の場合は、引き算した値は、0に近い値になるけど、uint8同士の引き算なので、オーバーフローを起こして、たまに大きい値になる。これが、ゼロ交差判定の代替として、うまく働くのか、こっちの方が、真面目に計算したものより、結果が綺麗に見えるっぽい。あまり意味はないけど、参考までに実装した

また、(1)の話からすると、網膜では、エッジ強調が行われるらしいので、簡単なエッジ強調フィルタも通してみることにした。で、以下がコード

import cv2
import numpy as np
import math


def sharpen(img):
    kernel = np.array([[-1,-1,-1], [-1,9,-1], [-1,-1,-1]], np.float32)
    return cv2.filter2D(img, -1, kernel)


def DoG(gray , W, sigma1, sigma2 ,beta = 1.0):
    assert(W > 0)
    assert(gray.dtype==np.uint8)
    k1 = cv2.getGaussianKernel(W , sigma1)
    k2 = cv2.getGaussianKernel(W , sigma2)
    dog_kernel = np.dot(k1 , k1.T) - beta*np.dot(k2 , k2.T)
    res = cv2.filter2D(np.float64(img), -1, dog_kernel)
    minNN = cv2.morphologyEx(res, cv2.MORPH_ERODE, np.ones((3,3)))
    maxNN = cv2.morphologyEx(res, cv2.MORPH_DILATE, np.ones((3,3)))
    zeroCross = np.logical_or(np.logical_and(minNN >= 0 , maxNN <= 0), np.logical_and(minNN<=0 , maxNN>=0))
    return np.uint8(zeroCross)*255


def DoG_int(gray , W, sigma1, sigma2 ,beta = 1.0):
    assert(W > 0)
    assert(gray.dtype==np.uint8)
    k1 = cv2.getGaussianKernel(W , sigma1)
    k2 = cv2.getGaussianKernel(W , sigma2)
    G1,G2 = np.dot(k1,k1.T) , np.dot(k2,k2.T)
    res = cv2.filter2D(img, -1, G1) - cv2.filter2D(img , -1, G2)
    return res


def gabor(gray,lam=10.0):
   assert(gray.dtype==np.uint8)
   sigma = lam*0.4
   gamma = 1.0
   imgs = [cv2.filter2D(gray , -1 ,cv2.getGaborKernel((int(lam), int(lam)), sigma, np.radians(5*n), lam, gamma, 0)) for n in range(36)]
   ret = 0*gray
   for img in imgs:
      ret += img
   return np.clip(ret , 0, 255)


if __name__=="__main__":
   img = cv2.cvtColor(cv2.imread("Lenna.png"),cv2.COLOR_BGR2GRAY)
   contours = []
   contours.append( DoG(img, 5, 2.2 , 2.2*1.6) )
   contours.append( DoG(sharpen(img), 5, 2.2 , 2.2*1.6) )
   contours.append( DoG_int(img, 5, 2.2 , 2.2*1.6) )
   contours.append( DoG_int(sharpen(img), 5, 2.2 , 2.2*1.6) )
   contours.append( DoG(img, 7, 7*math.sqrt(0.2459)/math.pi , 7*math.sqrt(0.9837)/math.pi , 0.891) )
   contours.append( DoG(sharpen(img), 7, 7*math.sqrt(0.2459)/math.pi , 7*math.sqrt(0.9837)/math.pi , 0.891) )
   contours.append( DoG_int(img, 7, 7*math.sqrt(0.2459)/math.pi , 7*math.sqrt(0.9837)/math.pi , 0.891) )
   contours.append( DoG_int(sharpen(img), 7, 7*math.sqrt(0.2459)/math.pi , 7*math.sqrt(0.9837)/math.pi , 0.891) )
   contours.append( cv2.Canny(img,100,200) )
   contours.append( cv2.Canny(sharpen(img),100,200) )
   contours.append( gabor(img) )
   contours.append( gabor(sharpen(img)) )
   contour = np.concatenate([np.concatenate(contours[0:4], axis=1),np.concatenate(contours[4:8], axis=1), np.concatenate(contours[8:12], axis=1)] , axis=0)
   cv2.imwrite("result.png",contour)


結果発表。合計12パターンの計算をやってみた
f:id:m-a-o:20190922185224p:plain

1行目:左2つが、真面目に計算したDoGで、エッジ強調なしとあり。右2つは、オーバーフローするDoGで、エッジ強調なし、あり
2行目:パラメータを変えて、一行目と同じ組み合わせの計算をしたもの
3行目:1列目と2列目は、Cannyエッジフィルタ(エッジ強調なし、あり)。3列目と4列目は、Gaborフィルタを5度ずつ回転して足したもの(エッジ強調なし、あり)


まぁ、エッジ強調なしで、Canny検出器が一番いいように思える。Cannyは、エッジ強調すると、(当然だが)全然ダメになる。Gaborも割といい結果になってる気がする。こっちは、エッジ強調なしでもありでも、それなりの成績を出す。

肝心のDoGフィルタは、まぁ、どう解釈したものか、よく分からない。とりあえず、エッジ強調の影響が全然見られない。2行目の1,2列目とか見ると、全然取れてないってこともないので、パラメータ次第では、よい結果が出ることもあるのかもなーという感じではある。

以上は、生理学側からのアプローチで、既に形成された神経回路の機能を解析した話。



一方で、そもそも、こういう都合の良い神経回路が、どうやって自然に形成されるのかという問題もある。

脳内情報表現への情報理論的アプローチ
http://www.jnns.org/previous/niss/1999/text/sakakaba.pdf

によれば、「1970年代には,視覚野に見られる特徴抽出細胞が生後の感覚体験を通じて形成されることを示唆する実験報告が相次いだ」らしい。

これに対して、色々なモデルを作って、シミュレーションで再現するという研究が色々出たっぽい。上の総説は2000年頃のもののようだけど、現在でも、この問題が解決してるのかは、よく分からない。

大体、基本になってるのはHebb則で、加えて、いくつかの補助的な機構が仮定されるという感じで、網膜〜一次視覚野でHebb則による学習が存在するのかは知らないけど、Hebb則自体は、海馬なんかでのLTPの発見によって、脳神経の学習の適切なモデル化だと考えられている。


1980年代後半に、Linskerという人も、この手のシミュレーションをして、Linskerは、シミュレーション結果から一歩進んで、Hebb則による学習が何を最適化しているのかという視点から、infomax原理を提唱した。今では、"入出力間の相互情報量を最大化する"と説明されてることが多い

From basic network princples to neural architecture:Emergence of spatial-opponent cells
https://www.pnas.org/content/83/19/7508

From basic network princples to neural architecture:Emergence of orientation-selective cells
https://www.pnas.org/content/83/21/8390

From basic network princples to neural architecture:Emergence of orientation columns
https://www.pnas.org/content/83/22/8779

Self-organization in a perceptual network
https://ieeexplore.ieee.org/document/36
http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.183.1750


Linskerの論文は、結局、主成分分析が出るみたいな話に見えるけど、主成分分析とEdge detectorやGabor filterが、どう繋がってるのかは定かでない

Linskerが、infomax原理からPCAを導くに当たって、二次の統計量までを考慮するという近似と、noiseはgaussianであるという性質を仮定したが、1996年の以下の論文では、その点を批判している。

The "independent components" of natural scenes are edge filters
https://www.ncbi.nlm.nih.gov/pubmed/9425547

論文の主な主張は、PCAではなく、ICAをやると、エッジフィルターが出るという話。PCAとICAの比較がされていて、PCAだと、エッジフィルターとは思えないものが得られている

論文でも、Gabor-likeという表現が使われているけど、本当にGaborなのか?とか、ICAだと何故Gaborが出るのか?とかは説明がない。この論文の結果を追試しようかと思ったけど、論文を読んでから長時間放置してたら、やる気がなくなったので、まぁ、そのうち

【補足】使用した画像データなんかは、既に入手できなくなってるっぽいけど、2013年の以下の論文でも同じようなことをやっていて、こっちはデータが現時点でまだ入手できるので、こっちのデータを使えば、まぁいいかという感じ

Double-Gabor Filters Are Independent Components of Small Translation-Invariant Image Patches
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3693455/


1990年代とかに、どうだったのか知らないけど、PCAを適用してみた事例は、今では沢山あって、まぁ概して、特にうまくいってるという感じもない。ICAだったら、どうかっていうのは、よく分からない。

同じ著者による論文で、自然音声に対して、同様のテクニックを適用したという論文が以下。

Learning the higher-order structure of a natural sound.
https://www.ncbi.nlm.nih.gov/pubmed/16754385

以下の論文も同じようなことをやってるっぽいけど、Gammatoneフィルター(?)っぽいものが得られるとか書いてある
Efficient coding of natural sounds.
https://www.ncbi.nlm.nih.gov/pubmed/11896400

Efficient auditory coding.
https://www.ncbi.nlm.nih.gov/pubmed/16495999

Gammatone filter
https://en.wikipedia.org/wiki/Gammatone_filter

以下は、二次視覚野もICAでイケソウというような話っぽい
Statistical model of natural stimuli predicts edge-like pooling of spatial frequency channels in V2.
https://www.ncbi.nlm.nih.gov/pubmed/15715907

弓矢の貫通力

12〜16世紀くらいのヨーロッパでは、プレートアーマーが流行って、騎士たちは、全身金属板の鎧を付けて戦ったと言われる。プレートアーマーの重量は、表面積と厚みと鉄の密度の積となる。人間の表面積は、2(m^2)前後で、鉄の密度が7800(kg/m^3)くらいだから、厚さ2mmとしても、重量は31.2kgという計算。実際のプレートアーマーも、これくらいの重さはあったようで、また、これ以上厚くするのは困難だったようである(Wikipediaによれば「鎧は種類にも拠るが、重量は数十キログラムにも及び、鎧だけでも20~30kg、兜や武器を含めると35kgを超えた」らしい)

一方、武器は色々あるけど、クロスボウやロングボウは、プレートアーマーを貫通することもできたらしい。百年戦争(1337〜1453)では、ロングボウを主体とするイングランドと、クロスボウを主体とするフランスの戦いで、連射性能に劣るクロスボウが負けたみたいな話があって、当時のヨーロッパではポピュラーな飛び道具だったらしい(特定の兵器の性能が勝敗を分けたみたいな逸話は、あまり信用できないとは思うけど、ロングボウは、一分間に12発撃てないと、一人前とは見なされなかったらしい。一方のクロスボウは、一分間に2〜3発しか撃てなかったと言われる)。携帯可能な火薬兵器であるマスケット銃が登場したのは、フス戦争(1419~1439)あたりからのよう。

ロングボウの射手は、180m先の標的に当てることが出来、90mなら、鋼鉄を貫いたとも言われる。現代の日本でも、堅物射貫という鉄板を射貫く競技があるそう(距離は、2〜30mで、厚さは1〜2mmとかっぽい?)だから、そういうこともあったのだろう。とはいえ、マスケット銃の銃弾と弓矢の矢弾では、弾丸速度や運動エネルギーに、かなりの差があるのも事実なので、矢弾の貫通力について、もう少し調べてみる



とりあえず、飛び道具は、飛んでいる間に、空気抵抗で減速するので、どれくらい減速するか知らないといけない。抗力は、空気密度と速度の二乗と断面積に比例するとしていい。銃弾の場合は、音速を超えるので、造波抵抗が効いてきて面倒くさいから、音速よりは、大分遅い矢が主な対象。で、時刻tにおける速度v=v(t)は、
 m v'(t) = -\dfrac{C_{D} \cdot \rho S}{2} v^2
の形の微分方程式を満たす。C_Dは抵抗係数で、ρは空気の密度、Sは投射物の断面積。抵抗係数は、Reynolds数に依存するけど、今は定数としておく。抵抗係数の計算は、流体力学が必要だけど、多分、簡単な形状でも手計算は難しい。

この微分方程式の厳密解は簡単にわかる。古典力学の演習問題とかでありそうなのに見たことないから、お前積分できたのかという感じだ。右辺のv^2にかかっている定数係数を、符号を除いて、k>0とすると
 v(t) = \dfrac{1}{kt/m + 1/v_0}
となる。もう一回積分すると、時刻tに於ける位置座標もわかる。一応、単位を確認しておくと、kについては[kg/m]となっているので、問題ない

断面積を見積もるために、矢の半径を知りたい。

Medieval Arrowheads from Oxfordshire
http://oxoniensia.org/volumes/2008/wadge.pdf

を見ると、シャフト部分(矢の胴体部分)の直径1cmくらいは、ありそうに見える

Applications of Physics to Archery
https://arxiv.org/abs/1511.02250
のTABLE1に、以下のショップで購入したらしい矢のデータが載っている。

Easton
https://eastonarchery.com/targets/

シャフト半径3.63mmで重量20.62gとある。現代の軽い素材っぽいので、重量は、さほど参考にならない。シャフト半径は、5mmくらいのものもあるっぽいので、直径1cmは、特に異常な数値というわけではなさそう。また、論文には、抵抗係数1.94±0.14とある。球の抵抗係数は、Reynolds数が大きい時、0.44くらいなのを考えると、かなり、大きいけど、他の論文を見ても、矢の抵抗係数は、2.0前後が多いので、2.0とする。

矢の重量や初速は、検索して見る限り、50~100gで、50m/s程度という数値が多いようなので、このあたりを採用する

空気の密度は、1.2(kg/m^3)を使うと、k=1.885e-5(kg/m)となる。矢の質量m=50(g)とすれば、k/m=0.000377(/m)。矢は落下するので、水平方向と鉛直方向は分けるべきだが、とりあえず水平に打ち出し、初速は、50m/sとした場合を考える。重力は無視すると、一秒後の速度は、49m/sくらいで、2秒後でも48m/sくらい。殆ど減速してない。大体、一秒後に50m飛び、二秒後に100m飛ぶ。m=20(g)とすると、2秒後の速度は、45.7m/sまで減る。ある程度の重量があれば、抵抗による減速は、あまり問題にならないっぽい。



次に、貫通する深さを単純に見積もる。基本的には、貫通する深さは、運動エネルギーによって決まる。微視的には、貫通する時、貫通される物質の分子結合を切断していってるので、運動エネルギーが、結合を切断するエネルギーに変換される。そういうミクロな詳細を考えなくても、ニュートン力学に基づく単純化したモデルで、貫通する深さが、運動エネルギーに比例することを示せる。

具体的には、「質量Mの物体に、質量mの物体を速度vで水平に打ち込んで、一定の抗力Fを受けながら、深さxまで侵入して停止した時の深さxを計算する(侵入される物体は十分厚いとして、反対側に貫通する可能性は考慮しない)」という問題を考えればいい。これは、単純化したモデルなので、現実の現象では、適切でない可能性があるけど、それほど悪い問題設定にも思えない。これを解くのは、高校物理の演習問題でしかない。

質量mの物体を、弾丸と呼称することにすると、弾丸が内部で停止した後、質量(M+m)の物体が、一定速度Vで動くことになる(侵入された方の物体が、床などから受ける摩擦などは考慮しない)。すると、単に、運動量保存則から
m v = (M + m)V
が成立する。弾丸が行う仕事は、侵入直前の弾丸の運動エネルギーと停止後の全体の運動エネルギーの差に等しい。つまり
 F x = \dfrac{1}{2} m v^2 - \dfrac{1}{2}(M+m) V^2
となる。結局
 x = \dfrac{M}{(M+m)F} \cdot \dfrac{1}{2} m v^2
を得る。Mが、mよりも十分大きいならば、実質的に
 F x \approx  \dfrac{1}{2} m v^2
として差し支えない。どちらにせよ、貫通する深さは、弾丸の運動エネルギーに比例する。

実際には、弾頭の変形や、熱などで散逸する分もあるし、抗力一定の仮定が成立しているかも全く明らかではないけど、矢と銃の(着弾時の)運動エネルギーを見ておくと、矢の場合は、重量100g、着弾速度50m/sとした場合、125(J)となる。銃は、ピンキリだけど、以下の拳銃を対象にすると、初速は345m/sで、9mm弾は直径9mm、重さが7.5〜9gの物があるようで、発射直後の運動エネルギーは、8gと350m/sを採用して、490(J)となる。
9mm拳銃
https://ja.wikipedia.org/wiki/9mm%E6%8B%B3%E9%8A%83

運動エネルギーのみの比較では、矢弾は、9mm拳銃弾の1/4以下しかない。直径も、矢のシャフト断面積と9mm弾断面積で大きな差はなさそう。

銃弾に使用される火薬の種類や量の記載は、見当たらなかったけど、現在、使用されてる火薬や爆薬の爆発エネルギーは、概ね、TNTの半分〜2倍程度の範囲に収まっており、拳銃の火薬重量は1g前後が多いようである。1TNT換算グラムは、4184(J)と定義されているので、ざっくり、10%程度のエネルギーが銃弾の運動エネルギーに変換されているということだろうか



少し、実測データを確認してみる。

[1] Material Culture and Military History: Test-Firing Early Modern Small Arms
https://journals.lib.unb.ca/index.php/MCR/article/view/17669/22312

に、16〜19世紀のマスケット銃(及び、比較として、20世紀の銃器)の性能が掲載されている。この中のGlock 80 pistolは、弾丸直径9.3mm、弾丸重量8.0g、初速360m/sで、30m地点に於ける速度は342m/sとなっている(運動エネルギーにして、468J)。この性能で、30m地点においたsteelを2mm貫通するらしい(Table2)。このsteelは、"2.8-3.0-mm thick cold-worked mild steel (hardness 290 HB)"と書かれている。

硬さ換算表(SAE J 417) 1983年改訂
https://www.nbk1560.com/~/media/PDF/ja-JP/products/pdf/30_Conversion_Table_of_Hardness.ashx?la=ja-JP
によると、ブリネル硬さ290は、引張強度(近似値)にすると、965(MPa)とされている

中世ヨーロッパの鎧は、この実験で使用したsteelほどは硬くなかったかもしれない。鎧の硬さは、上論文[1]で引用されていた以下の論文などで確認することができる。

The metallography and relative effectiveness of arrowheads and armor during the middle ages
https://www.sciencedirect.com/science/article/pii/104458039290109U

この論文のTable2のThe Hardness of armorでは、1500年前後のBreast plateのビッカース硬度が、120〜220と書いてある。上の硬さ換算表の下限を超えているけど、ビッカース硬度120は、引張強度400MPa程度に相当すると思われる。ビッカース硬度220は、730MPa程度になっている。Table1には、鎧の厚さのデータもあるが、2〜3mmというのは、概ね正しいらしい。

一応、中世鎧の材質は、論文[1]で使われた鋼鉄プレートの50〜75%くらいの強度しかなかったと言って良さそう(?)。既に見たように、矢弾の運動エネルギーは、Glock 80 pistolの銃弾の1/4程度以下なので、例え、鎧の強度が半分だったとしても、矢弾と銃弾で抗力に差がなければ、矢は、厚さ1mmの鎧を貫通するのが、精々だっただろう。実際には、矢弾と銃弾は形状が違うので抗力が同じではないだろう(直感的に、矢の方が尖ってるので刺さりやすい)し、抗力が速度に依存する場合は、矢弾は銃弾より随分遅いので、抗力も小さいだろう。


矢弾のことは少し置いておいて、多くの試験データがある銃弾について、抗力一定の仮定の妥当性を検討する。論文[1]には、他にも多くの試験データがあり、このデータから、プレート着弾時の運動エネルギー/破断した体積を、算出した結果、以下のようになった。

"破断エネルギー"は、体積当たりの破断に要したエネルギーで、平均的な抗力を表していると考えられる(単位面積当たりの力と、単位体積あたりのエネルギーは、同じ次元を持つことに注意。3000MJ/m^3は、3000N/mm^2に相当する。また、破断エネルギーは、勝手に作った造語で、一般に通用している名称ではない)。マスケット銃の弾丸は、総じて、重く、でかい傾向にあるので、破断エネルギーを見たほうがいい

銃器 試験 運動エネルギー "破断エネルギー" 着弾速度(m/s)
"Doppelhaken" G 284 PS100 1779.6J 3138(MJ/m^3) 305
"Doppelhaken" G 358 PS100 2992.7J 2335(MJ/m^3) 349
Matchlock LG 1514 PS30 1241.7J 3467(MJ/m^3) 378
Matchlock LG 1514 PS100 605.7J 3382(MJ/m^3) 264
Wheellock RG 33 PS30 2333.2J 3347(MJ/m^3) 394
Wheellock RG 33 PS100 1238.0J 2664(MJ/m^3) 287
Wheellock RG 117 PS30 660.2J 2779(MJ/m^3) 349
Wheellock RG 117 PS100 307.0J 2539(MJ/m^3) 238
Wheellock RG 272 PS30 1874.9J 3898(MJ/m^3) 342
Wheellock pistol RP 2895 PS30 602.4J 2754(MJ/m^3) 355
Flintlock STG 1287 PS30 2269.8J 2560(MJ/m^3) 406
Flintlock STG 1287 PS100 1166.1J 2630(MJ/m^3) 291
Flintlock STG 1288 PS30 2032.8J 3131(MJ/m^3) 390
Flintlock STG 1288 PS100 1055.3J 2438(MJ/m^3) 281
Flintlock STG 1316 PS30 2458.3J 3368(MJ/m^3) 391
Flintlock STG 1316 PS100 1324.5J 2722(MJ/m^3) 287
Flintlock STG 1318 PS30 2806.5J 2917(MJ/m^3) 426
Flintlock STG 1318 PS100 1438.6J 2991(MJ/m^3) 305
Flintlock pistol STP 1128 PS30 753.8J 2633(MJ/m^3) 323
アサルト・ライフル58 PS100 2801.5J 4824(MJ/m^3) 770
アサルト・ライフル77 PS100 1375.0J 5987(MJ/m^3) 874
Glock80ピストル PS30 467.9J 3444(MJ/m^3) 342

最後の3つの銃器は、20世紀のもので、他とは、銃弾形状が異なっている。他のマスケット銃は、当時に合わせて、球状の弾丸を使用したらしい。PS30とPS100は、距離30mから撃ったか、距離100mから撃ったかの差で、衝突時の速度も論文に記載されている。破断した体積は、弾丸の半径の2乗×円周率×貫通した長さで計算した

貫通した厚さは、mm以下の単位がなく、1〜4mmのものが殆どで、破断体積あたりのエネルギー(=単位面積あたりの抗力)は、若干ばらついてはいる。元々、マスケット銃は、弾丸がまっすぐ飛ばないらしいので、正面から当たっているものも少ない可能性が高い。現代でも、プレート厚さなどの条件を揃えて、50%の確率で貫通できる速度V50という指標を、よく使うようで、古い銃だと、なおさら、ばらつきが大きくても仕方ないのかもしれない。それでも、アサルト・ライフルを除けば、2000〜4000(MJ/m^3)に集中している。マスケット銃や拳銃のような低速弾に於いては、抗力一定の仮定に基づく貫通深度計算は、そこそこ妥当らしいことが、窺える(抗力の大きさそのものは、弾丸形状に依存するだろうけど)

ちなみに、参考までに、鉄を溶かすのに、必要なエネルギー量を計算すると、融点が1538度、25度に於ける熱容量が25.10J/mol/Kなので、融点まで熱するのに、ざっくりと、37.7kJ/mol程度だろうと推測し、それ以外に、融解熱は13.81kJ/molらしいので、合計では、51.51kJ/molとなる。上の破断エネルギーと単位を合わせると、原子量は、55.85なので、モル質量を55.85g/molとして、密度は、7874kg/m^3だから、7874*51.51e-3/55.85e-3=7262(MJ/m^3)となる。



抗力について、もう少し定量的な議論が存在する

[2] Analysis on the resistive force in penetration of a rigid projectile
https://www.sciencedirect.com/science/article/pii/S2214914714000658

の式(1)に、

 F = \dfrac{\pi d^2}{4}(A \sigma_y N_1 + B \rho V^2 N_2)

という式が書いてある。S=\dfrac{\pi d^2}{4}は、断面積なので、特に問題はない。

Vは着弾速度で、A,Bは、材質に依存する無次元定数。

\sigma_y,\rhoは、降伏応力(yield stress)、密度で、やはり材質に依存する。

N_1,N_2は、弾丸形状及び摩擦係数から決まる無次元定数としている。

第一項をquasi-static resistive forceで、第二項をdynamic resistive forceと呼んでいる。第二項が0であれば、抗力一定で計算するのと、同じことになる。


この式の導出は、以下の論文で行われたようである。空洞膨張理論(cavity exansion theory)という、土木系の業界で使われてるらしい解析法に基づいているっぽいけど、そんなものは知らないので、理論的観点から見た、この式の妥当性については、私には判断できない
Penetration into soil targets
https://www.sciencedirect.com/science/article/pii/0734743X9290167R


元々、土壌の破壊を想定して作られた理論らしいので、金属の破壊に適用するのは、正しくないようにも思える(ベースになっているMohr-Coulombの破壊基準というのは、脆性破壊を起こすような物質に適用されるものらしいから、コンクリートに銃弾を打ち込むとかなら、正しい扱いかもしれない)けど、この式は妥当だと仮定する。プレート内部に侵入した弾丸は急激に減速していくので、上式のパラメータをまとめて
F = \alpha + \beta V^2
と書いた時、貫通する深さxは
x = \dfrac{m}{2 \beta} \log \left(1+\dfrac{\beta V^2}{\alpha} \right)
と計算できる。βー>0の時は、\log \left(1+\dfrac{\beta V^2}{\alpha} \right) \approx \dfrac{\beta V^2}{\alpha}となり、
x = \dfrac{1}{\alpha}\frac{m V^2}{2}
という既に計算したのと同じ結果を得る。


論文[2]の式(1)のパラメータを第一原理的に決めるのは、しんどいけど、論文[2]のTABLE2のデータを利用して、概算値を得ることができる。TABLE2の値以外に、論文[2]の式(2a)(2b)(2c)に従って、nose shape factorを計算することが必要になる。なお、式(2c)に何故か負符号がついてるけど、これは間違い

これらの値は、弾丸形状に依存するので、簡単のため、球状弾を仮定すると、前方半分で(2a)の積分を計算して
N_1 = 1 + \dfrac{8 \mu}{d^2} \displaystyle \int_{0}^{d/2} \sqrt{\dfrac{d^2}{4} - x^2} dx = 1 + \dfrac{\pi \mu}{2}
となる。鉄と鉄の摩擦係数は、0.5くらいらしいので、大体、このnose shape factor N1は、鉄の球状弾では、約1.8くらいになるだろうと見積もれる。

この値と、TABLE2の値から、降伏応力400MPaの鋼鉄に於いて、断面積をSとすると、α/S=3130(MPa)=3130(MJ/m^3)という値を得る。これは、論文[1]の試験結果で見た、マスケット銃や拳銃の破断エネルギーと近い値である(材質が同一というわけではないのだけれど)。同様に計算して、球状弾では
N_2 = \dfrac{1}{2} + \dfrac{\pi \mu}{8}
を得る。TABLE2の値と合わせて、降伏応力400MPaの鋼鉄に於いて、β/S=6200(kg/m^3)程度と分かる。

従って、特に、β/αは、大体、2.0e-6程度である(単位は、s^2/m^2で、速度の-2乗の次元を持つ)。この比を採用すれば、アサルト・ライフルと同程度の着弾速度V=800m/sに於いて、\beta V^2/\alphaは、1.28となり、これは、dynamic resistive forceとquasi-static resistive forceの比を表している。弾丸の速度は減速していくので、dynamic resistive forceは、どんどん小さくなるが、ある程度、減速するまでは、dynamic resistive forceが無視できない寄与をすることが分かる。

一方、矢弾の速度V=50m/sの場合は、\beta V^2/\alphaは0.005であり、着弾速度が50m/sであれば、dynamic resistive forceは、最大でもquasi-static resistive forceの0.5%くらいしか寄与せず、抗力一定の仮定は十分妥当だと結論付けられる。これは、球状弾の場合の試算であるけど、矢弾でも、それほど大きく事情は変らないだろう。


論文[1]の試験で用いたプレートの降伏応力は、400MPaではないので、念の為、論文[1]の試験結果からも、β/αを見積もってみる。この試験に於いては、着弾速度が、マスケット銃や拳銃と、アサルト・ライフルで大きく違い、それに応じて、破断エネルギーも異なっているように見える。マスケット銃や現代の拳銃は、着弾速度が300m/sであり、破断エネルギーは3000MJ/m^3前後であるのに対して、アサルト・ライフルでは、着弾速度は800m/sであり、破断エネルギーは5〜6000MJ/m^3程度である。仮に、着弾速度が300m/sと800m/sで、他の条件が同一なら、(破断エネルギーの計算は、着弾時運動エネルギーを、貫通深さと弾丸断面積の積で割ったものだったので)貫通深度が3.5倍になるとした場合、β/αは、6.88e-6程度であることが計算できる。オーダーは、さっきの計算と合っている。本当は、ライフル弾は、球状ではないので、α、βの比だけでなく、個別の大きさの違いも、考慮しなければならない。


マスケット銃や矢弾程度の速度では、抗力一定の仮定は妥当そうであるけど、cavity expansion modelでは、先端形状の違いが、抗力に与える影響を見積もることが出来る。具体的には、論文[2]の式(2.a)によって、先端形状の違いによるquasi-static resistive forceの大きさの比を見積もる。矢の先端は、尖っており、円錐形状であるとして計算すればいい。積分区間を、弾丸の場合と揃えないとアンフェアな結果になってしまう(長い領域で積分するほど、抗力は大きくなる)ので、シャフト根本から、矢の先端の長さをsとして、sがd/2以下の時と、d/2以上の時に分けて計算する(dはシャフトの直径で、シャフト部分は、単なる円柱なので、y=y(x)は定数関数)
N_1 = \begin{cases} 1+ 2\mu(1- \frac{s}{d}) & (s < d/2) \\ 1+\mu\dfrac{d}{2s} & (s \geq d/2) \end{cases}

全く尖っていない状況(s=0の時)では、N_1 = 1 + 2\muとなって、球状弾丸より大きな値を取る。境界部分にあたるs=d/2の時は、N_1=1+\muくらいで、更にsが大きくなるにつれて、1に近付いていく。ざっくりと、球状弾丸と比較すると、75〜50%くらいの範囲になると推定される。何にせよ、quasi-static resistive forceについて、矢弾は銃弾より小さい(定性的には意外でも何でもないけど)。

上の方で、「鎧の強度が半分だったとしても、矢弾と銃弾で抗力に差がなければ、矢は、厚さ1mmの鎧を貫通するのが、精々だっただろう」というように書いたけど、実際には、矢弾の方が、抗力が小さくなり、運動エネルギーの小ささを補ってくれる。暫定的な結論として、矢弾の貫通力は、その先端形状による抗力の減少に起因する部分が大きいのだろうとなる。

あと、そもそも厚さ2〜3mmだと、穴はあくけど、矢が通り抜けられないという状況も発生するはずで、これも「貫通した」と言えなくもないけど、人体へのダメージという観点からすれば、矢を防いだと言えなくもない。穴があくだけというのと、矢が完全に通り抜けるのとでは、穴の断面積が異なるので、抗力の大きさも異なってくる。「貫く」というのが、どちらの意味で使われてるのか、はっきりしないけど、プレートアーマーが流行ってたということは、穴が開くだけで済むという事態は、結構あったのでないかと思われる。銃弾でも同じ問題はあるけど、こちらは、貫通する=通り抜けてしまうというイメージがある

もっとデータが欲しいとこだけど、下手に実験すると、銃刀法違反になってしまいそうで、うざい



近代戦では、徹甲弾で、戦艦や戦車の装甲を貫くということが行われるらしい。Wikipediaに、一部の徹甲弾のデータが載っていた
APFSDS#貫徹力
https://ja.wikipedia.org/wiki/APFSDS#%E8%B2%AB%E5%BE%B9%E5%8A%9B

実戦で使われる装甲は、様々なものがあり、簡単に評価はできないけど、「装甲を貫く力は、均質圧延鋼装甲(RHA:Rolled Homogeneous Armor)を貫ける厚さで表現される」とある。RHAは"全体が均質な圧延鋼板で作られた装甲"らしいので、鋼鉄性のプレートみたいなもんと思っていいのだろう。引張強度は、1GPaとか書いてあったりする(何種類かあるっぽい?)ので、論文[1]で使われたプレートと、ほぼ同程度の強度だろう

1961年、ソ連が開発したBM-3は、飛翔体重量4kg、口径115mm、初速1615m/secで、距離2000mに於いて、RHAを236mmを貫通する、とある。着弾速度は不明だけど、1500m/secとして、論文[1]のデータに対して、やったのと同じ要領で、(弾丸直径は口径に等しいとして)、破断エネルギーを計算すると、1835.8(MJ/m^3)という値を得る。論文[1]のデータと比較すると、アサルト・ライフル以上に高速であるのに、破断エネルギーは、むしろ小さくなっている。この情報が本当かどうかも容易に確認することはできないけど、別のモデルを必要とすると考えるべきかもしれない


衝突速度が上昇すると、侵徹メカニズムが変わるらしきことは、以下のような論文でも書かれている。どこらへんで変わるかは、弾丸材料と装甲材料に依存するっぽいけど、大体、1000(m/s)くらいっぽい(?)。ただ、高速になっても、破断エネルギーは小さくならないように思えるのだけど、まぁ、データがWikipediaの記述のみという信頼性低い状態なので、何とも言えない

Transition from Nondeformable Projectile Penetration to Semihydrodynamic Penetration
https://ascelibrary.org/doi/abs/10.1061/(ASCE)0733-9399(2004)130:1(123)

Long-rod penetration: the transition zone between rigid and hydrodynamic penetration modes
https://www.sciencedirect.com/science/article/pii/S2214914714000397

hydrodynamic modelとか、hydrodynamic penetrationというのは、Garrett Birkhoffらが提案した、成形爆薬の侵徹理論(1948年)に遡り、現象論的に改良されたモデルが、AlekseevskiとTateによって独立に提案され(1960年代)て、Alekseevski-Tateモデルと呼ばれることもある。

[Birkhoffの論文] Explosives with Lined Cavities
https://aip.scitation.org/doi/10.1063/1.1698173

[Alekseevskiの論文] Penetration of a rod into a target at high velocity
https://link.springer.com/article/10.1007%2FBF00749237

成形爆薬のメタルジェットは、金属の音速(〜5000m/s)を超える速度で装甲を侵徹するとされる。これくらいの速度だと、金属を流体のように扱うことが適切になるとして、モデル化されているので、hydrodynamic modelと呼ばれる。

これらは大分単純化されたモデルだけど、cavity expansion modelでは、弾丸は無限に硬い剛体として扱ったのに対して、Alekseevski-Tateモデルでは、装甲側も弾丸側も、有限の硬さを持つ物質としてモデル化されている(Birkhoffらのモデルは、完全に流体として扱って、硬さは影響しないというモデルらしい)。

しかし、徹甲弾の場合は、速度が2000m/sを超えることは、ほぼないと思われ、金属中の音速に比べると大分遅く、Alekseevski-Tateモデルを徹甲弾の侵徹に適用していいのかは、疑問ではある

On the hydrodynamic approximation for long-rod penetration
https://www.sciencedirect.com/science/article/pii/S0734743X9800044X


現在の所、人類は、主に大気圏内で兵器を打ち合っているが、未来には、宇宙で兵器を打ち合っているかもしれない。

Fred Whippleという人は、1947年に、Meteorites and space travelという短い記事の中で、ミリグラム単位の微小なmeteoroidから、宇宙船を防御するための提言を行い、現在、Whippleシールドと呼ばれている。宇宙空間では空気抵抗がないので、この場合の衝突速度は、3000〜20000(m/s)程度までの領域が想定されるらしい。

SFとかで、運動エネルギー弾を撃ちあっている場合、衝突速度は、これ以上のオーダーに達しているかもしれない。