エムデン方程式 - 空気の場合

レーン・エムデン方程式(Lane-Emden equation)を 空気 (2原子分子)について解いてみる。中心が27°, 1気圧 という条件で具体的な数字を出した。
2原子分子なので ɤ = 7/5, N = 5/2. とする。水素の場合と同じようにオイラー法でスプレッドシートを使って数値的に解く。
宇宙空間に水星の質量の半分程度の空気をあつめると自己重力でまとまって、中心は人間が生存できるとなった。だれか験算しておくれ。

エムデン方程式の具体的な形

水素の場合をexcelで計算してグラフに描いた。

球対称で平衡状態にある星(自己重力ガス球)の構造は、密度と圧力の間にポリトロピック関係が成り立つとした場合に次のレーン-エムデン方程式 (Lane-Emden equation) であらわされる。ガス球が水素からできている(単原子理想気体)とするとɤ=5/3 (N=3/2) となる。と、ここまでは、教科書[1]に書いてある話。

では、具体的にはどうなるかを計算してみた。
N=3/2の場合のLane-Emden equationをMaximaを使って変形して、Dの二階微分をDとDの微分であらわすとこうなる。

これをスプレッドシートでオイラー法を使って数値的解くとこんな感じ。より分かりやすく、Dを圧力Pと密度ρに読み替えてみた。


[1] 福江 純・和田桂一・梅村雅之、宇宙流体力学の基礎 (日本評論社、2014)

自己重力等温ガス球の構造

構造が中心に対して球対象だとし、時間的に変化しないとして計算。
星内部に半径方向に微小な円柱を考え、上面と下面の全圧力と円柱内のガスへの重力のつり合いを考えると次の式になる。
(3.29)
また、半径rより内側に含まれる質量Mrは、
(3.31)
さらに、ガスの温度がどこでも一定として、理想気体の状態方程式を当てはめると次の式になる。
(3.46)
ここで、次のような変数変換をする。
(3.47)
(3.48)
すると式はこうなる。[1]
(3.49)
以上 ・・・
いや、これではどんな形状だかさっぱりわからん。
この2階微分方程式は解析的には解けないということなので数値計算する。と言っても、オイラー法ならスプレッドシートでも計算できるので、まずはそれであたりをつけてみる。
まずは、式3.49 を変形して、Dの2階微分 (D")をDの1階微分(D')とDの式 f(D', D)表現する。高校数学のレベルだがチョンボはしたくないので数式処理システム(maxima)に手助けしてもらう。なお、maximaでDと書くと定数扱いされてしまうので、ξの関数 D(ξ) として表現している。

あとは、次の式でξを少しずつ増やしながら計算していけばよい。なお、教科書では、境界条件としてξ=0 の時、D=1で D' = 0と書いてあるが、D"の計算でξ=0は特異点になっているので使えない。まあ、中心付近は重力の効果はほぼ無くて一様だろうから、わずかに離れたところでもD=1でD'=0だとして計算できる。
D (ξ + d) = D(ξ) + d * D'(ξ)
D' (ξ + d) = D'(ξ) + d * D"(ξ)
D" (ξ) = f (D'(ξ + d), D(ξ + d))

ということで計算してみたのだが、どうもDの減少の仕方がが思ったより小さい。グラフにするとこうなる。

青線がDで赤線が1/ξ^2 。D∝1/ξ^2 だとξが増えるに比例して質量Mが増え続けるということなので、どこかが等温ガス球の境界とは言えないということだよね。
定常状態にある理想気体からなる自己重力等温ガス球というのは星とは言えないのかぁ。


[1] 福江 純・和田桂一・梅村雅之、宇宙流体力学の基礎 (日本評論社、2014)

気体の平均分子量の計算

Q. 空気(分子量28のN分子78%, 分子量32のO分子21%, 原子量40のAr 1%)の平均分子量はいくらか.

誤答: 0.78*28 + 0.21*32+0.01*40 ≒ 28.96
正解: 1.0/ ( 0.78/28 + 0.21/32 + 0.01/40) ≒28.84

なぜか？ それは、出題者が78%,21%,1%を重量パーセントとして出題したから。粒子数あるいは体積のパーセントなら最初の式になる。

音波を対象としたPMLの実験的考察 その3. 計算実験結果

3. 実験手法

PMLの性能は、反射波の大きさを剛な境界の場合と比較することで評価する。単純に受音点で音圧を観測した場合、反射が小さい場合に直接波が尾を引く残余と反射波との区別が難しくなることが予想される。そこで、測定時間内には反射波が戻ってこないほど境界の位置を遠のけた場合の測定値を基準として、それとの差を反射波として抽出する。
斜めに音波が入射する場合の性能を評価するため、点音源を用い受音点を3か所設けた。剛な境界の場合、境界に対して音波が垂直に入射し反射して到達する点 P1、およそ４５度で入射して反射が到達する点 P2、およそ３０度で入射して反射が到達する点 P3 の３点である。音波はいったんPML層の内部に入って反射してくる。PML層が厚くなると厳密には反射角度が異なってくると想定されるが、ここでは無視できるものとして扱う。
具体的には、X 999 cm, Y 599 cm, Z 799 cm の天井以外は剛な境界で囲まれた空間に、下図のように音源と受音点を配置した。Spacial stepは1cm、 time step は 0.01ms。音速cは330m/s、空気密度は1.2047291 kg/m3を使用。CFLは0.33 である。

Fig 3-1 配置図
反射波の大きさには、音圧の絶対値の最大値を使用する。より具体的には、P1では計算開始から11.50msから14.50msまでの3.0ms、P2では16.60msから19.60ms、P3では23.90msから26.90ms の音圧の絶対値の最大値を用いた。
音源は、音源となるセルの音圧を1kHzのsin波で1周期分だけ強制的に与えるという簡便な方法を採用した。
計算は単精度で行った。


4. PMLの層数と反射量

PMLの層数のみを変更し、反射量の関係を調べる。
PMLの層数を8から256まで変え、その他は次の条件で一定とした。Absorption coefficients σには式5を採用し、乗数Mを2, σMAX は 1.0 。音源は、1kHzのsin波1周期。
計算した各受音点での反射波の絶対値のピークを剛な境界の場合を1として現したものが次の図である。

Fig 4-1 PML層数と減衰

Fig.4-1 で示すようにPML層数を指数関数的に増加させることで反射波は指数関数的に減少する。100層あたりで反射波の減少が底を打つが、これはPML層に入射しつつある音圧に対して、単精度浮動小数点の有効桁数(十進で6桁強)では不足するほど反射波の音圧が小さくなりすぎたことにより発生するノイズが卓越するためと考える。
Fig.4-1で興味深い点は、PMLが10層以下であるとPMLへ浅い角度で入射する音波の減衰が特に悪いが、層数を増やすことで顕著に改善される点である。
具体例として、受音点P2における16.6msから19.6msまでの反射波の音圧波形を示す。波形を見やすくするため、音圧のスケールはすべて異なる。青線でしめす剛な境界からの反射波形と比べると、PMLが8層の場合の反射波形はおおむね元の形状をとどめているが、PMLが64層になると反射波形はノイズまみれとなっている。

Fig. 4-2 受音点 P2 での反射波の波形

5. Absorption coefficientsと反射量の関係

Absorption coefficients σに式５を採用した場合、乗数Mとσ max の二つのパラメータとPMLでの反射量の関係を検討する。
fig.5-1は、σmaxを1.0に固定し乗数を1から4まで変えた時の反射量を示す。PMLの層数は16で固定である。

Fig. 5-1
入射角が直角(P1)の場合は乗数が増えるほど反射量が減るが、角度が浅くなるほど反射量が最低となる乗数があることが見える。σ maxが1.0の場合、角度によらない結果を求めた場合では乗数として2が最適で、その場合の反射は1.0E-4以下 (吸収率99.99% 以上) と読み取れる。さらに、σ maxを変えて変えてみる。実は、Absorption coefficents は1.0を超えても吸収層として働く。

Fig. 5-2

Fig. 5-3
fig. 5-2, fig. 5-3は、σ maxを1.2にした場合とσ maxを2.0にした場合の図である。
σ maxを大きくするほど入射角が浅いケースでの反射量が最低となる乗数が増えるとともに、反射量の最低値も減少することが見て取れる。
Fig 5-4は、σ max によって減衰量がどのように変化するかと異なる乗数ごとに追った図である。

Fig. 5-4

角度に因らない結果を求めP1, P2, P3 の3点での最悪値を最良とすることを考えよう。
乗数２の場合、σ max が0.8 付近がほぼ最良で、それよりもσ maxを大きくした場合 P3での値は改善してもP2, P1での値が悪化している。乗数が3の場合、グラフからσ max が 1.4～1.6 あたりが最良で、その時の反射は 1.0E-5 と推測される。
乗数4の場合、σ max 2.0 より大きいところに最良点がありそうに見える。では、σ maxは どこまで大きくすれば良いか？　式２に従えば、absorption coefficients が2.0を超えると発振してしまう。が、式４を使っている場合、音圧および粒子速度計算点のσが 2.0を超えなければ問題ないはずである。そこで、計算に使用するσの最大値を横軸にとって反射量を図にしたものが、fig 5-5である。図から、乗数が４の場合 σ max に 2.1 にすることが最良と読み取れる。では、乗数をさらに増やすとどうか。乗数が５の場合、P3地点 (浅い入射角)での値を十分下げきれない。Fig 5-3 と合わせてみるに、乗数が4が最良と読み取れる。

Fig. 5-5

結果
PML (Perfectly matched layer) を具体的に構成するために、absorption coefficients のσ max と 乗数M の最適値を計算実験で探った。
M とσ max は、相互に関係しておりセットで最適値を探索する必要があることを見出した。
単精度の計算でPMLが16層である場合、乗数M=４、σ max = 2.1 が最適で、剛な境界の場合に対して反射を 1.7 E-6 程度に減衰させることができることを見出した。


[1] Ludovic Métivier, Laurence Halpern. Perfectly Matched Layers equations for 3D acoustic wave propagation in heterogeneous media. WAVES 2011, Jul 2011, Vancouver, Canada. ⟨hal-00826617⟩
[2] 豊田正弘 編著. FDTD法で視る音の世界, コロナ社 (2015), p12-14
[3] 豊田正弘 編著. FDTD法で視る音の世界, コロナ社 (2015), p9-10
[4] 宇野 亨編著　何 一偉, 有馬 卓司 共著、数値電磁界解析のためのFDTD法 -基礎と実践-、コロナ社(2016)