RD-E：4300 多項式EOSを用いた理想気体のモデル化

この例題の目的は、数値的圧力、内部エネルギー、理想気体材料則の音速をプロットすることにあります。

理論界との比較がなされます。AbsoluteとRelative定式化のコントロールカードが使用されます。

多項式EOSはRadioss流体力学的圧力を計算や理想気体のモデル化に用いられます。これは圧縮で3次、膨張で線形です。(1)

P = C_{0} + C_{1} μ + C_{2} μ^{2} + C_{3} μ^{3} + (C_{4} + C_{5} μ) E

ここで、(2)

E = \frac{E_{int}}{V_{0}}

および (3)

μ = \frac{ρ}{ρ_{0}} - 1

これらの定式化の出力を理論解と比較するため圧縮 / 膨張の単純なテストが行われます。

使用されるオプションとキーワード

理想気体
多項式 EOS（/EOS/POLYNOMIAL）
絶対 / 相対定式化
圧力シフト
流体材料（/MAT/LAW6 (HYDROまたはHYD_VISC))
強制変位（/IMPDISP）
Lagrangian境界条件（/ALE/BCS）
それぞれの面の節点は強制変位で移動されます（/IMPDISP）
要素圧力、密度と内部エネルギー密度は、時刻歴応答ファイルに保存されます。

入力ファイル

必要なモデルファイルのダウンロードについては、モデルファイルへのアクセスを参照してください。

モデル概要

このテストは球状の膨張と圧縮を受ける理想気体の要素体積からなります。

初期条件を以下に記します：

$P_{0} = 1 e 5 [Pa]$
$V_{0} = 1000 [m^{3}]$
$ρ_{0} = 1.204 [\frac{k g}{m^{3}}]$
$μ_{0} = 0$

流体は理想気体と仮定します。体積は純圧縮 $- 1 < μ < 0$ とそれに続く膨張 $μ > 0$ （図 3）を考慮するために3方向に変化させます。

このテストは1つのALE要素（8節点ソリッド）と多項式EOSでモデル化されます。

圧力、内部エネルギーと音速の変化が数値計算出力と理論解の間で比較されます。

長さは/IMPDISPで修正され、その

V

と図 3への影響がプロットされます（図 3）。

シミュレーションの反復

1つのALEソリッド要素が用いられます。材料は、ソリッドの節点がLagrangeであると定義することにより要素内に閉じ込められます。それぞれの面に対して、4節点に強制変位が法線方向に与えられます。

多項式 EOS

多項式EOSが/EOS/POLYNOMIALの中で流体力学的圧力の計算に用いられます。これは圧縮で3次、膨張で線形です。(4)

P = C_{0} + C_{1} μ + C_{2} μ^{2} + C_{3} μ^{3} + (C_{4} + C_{5} μ) E

Cパラメータは流体力係数と呼ばれる入力パラメータです。材料の挙動に対する仮説により、これらの係数を決めることが可能となります：

非圧縮性気体
線形弾性材料
理想気体

この例題は理想気体のモデル化のみに焦点を当てます。

結果

理論解

このセクションの目的は、圧力、内部エネルギーと音速を1つのパラメーター

V

または

μ

の関数としてプロットすることです。

圧力：
理想気体の圧力は以下で与えられます：(5) $P V = (γ - 1) E_{int}$

したがって：(6) $d P (V, E_{int}) = {\frac{\partial P}{\partial V} |}_{E_{int}} d V + {\frac{\partial P}{\partial E_{int}} |}_{V} d E_{int}$

Radiossは内部エネルギー変化の計算に等エントロピー過程の仮説を仮定します：(7) $d E_{int} = - P d V$

この理論は以下の微分方程式を与えます：(8) $\frac{d P}{d V} = - \frac{γ P}{V}$

これは $y' + \frac{γ}{x} = 0$ という形式を持ち、一般解は次のとおりです：(9) $y = C s t . x^{- γ}$

圧力はポリトロープな状態にもあります：(10) $P V^{γ} = P_{0} V_{0}^{γ}$ (11) $P (V) = P_{0} {(\frac{V_{0}}{V})}^{γ}$

ここで、 $γ$ は材料定数です（熱容量の比）。原子の気体では、 $γ$ =1.4。空気は主に2原子の気体からなるので、通常空気のガンマは1.4に設定します。
内部エネルギー：
式 5 と式 11から直接の結果を導けます：(12) $E_{int} (V) = \frac{P_{0} V_{0}^{γ}}{(γ - 1) V^{γ - 1}}$
音速：
理想気体の音速は：(13) $c = \sqrt{\frac{γ P}{ρ}}$

式 11 は体積に関する表現を与えます：(14) $c = \frac{γ P {}_{0}}{ρ_{0}} {(\frac{V_{0}}{V})}^{γ - 1}$

圧力、内部エネルギーと音速は

V

と

μ

の両方で表現されます。

表 1. 理論解
圧力（Pa）		内部エネルギー密度（J）		音速（m/s）
$P^{R E F} (V)$	$P^{R E F} (μ)$	$ρ e^{R E F} (V)$	$ρ e^{R E F} (μ)$	$c^{R E F} (V)$	$c^{R E F} (μ)$
$P_{0} {(\frac{V_{0}}{V})}^{γ}$	$P_{0} {(1 + μ)}^{γ}$	$\frac{P_{0}}{(γ - 1)} {(\frac{V_{0}}{V})}^{γ - 1}$	$\frac{P_{0}}{(γ - 1)} {(1 + μ)}^{γ - 1}$	$\sqrt{\frac{γ P_{0}}{ρ_{0}} {(\frac{V_{0}}{V})}^{γ - 1}}$	$\sqrt{\frac{γ P_{0}}{ρ_{0}} {(1 + μ)}^{γ - 1}}$

対応するプロットは：

材料コントロールカード

/MAT/LAW6 (HYDROまたはHYD_VISC) と/EOS/POLYNOMIALは静水圧の計算にこの式を使用します。絶対値または相対変化を考慮することが可能です。材料は理想気体と想定します。以下のケースが調査されます。

表 2. 理想気体のモデリング定式化
ケース	数学的モデリング	圧力	エネルギー
1	$P (μ, E)$	絶対値	絶対値
2	$Δ P (μ, E)$	相対値	絶対値
3	$Δ P (μ, Δ E)$	相対値	相対値
4	$P (μ, Δ E)$	絶対値	相対値

音速と時間ステップ

LAW6は流体の通常表現で音速を計算します：(15)

c^{2} = \frac{d P}{d ρ}

これは

μ

の関数として書くことができます：(16)

μ = \frac{ρ}{ρ_{0}} - 1 \Rightarrow \frac{1}{d ρ} = \frac{1}{ρ_{0}} \frac{1}{d μ}

したがって：(17)

c^{2} = \frac{1}{ρ_{0}} \frac{d P}{d μ}

P

の内部エネルギー

E

と

μ

に関する全微分は：(18)

d P (μ, E) = {\frac{\partial P}{\partial μ} |}_{E} d μ + {\frac{\partial P}{\partial E} |}_{μ} d E

等エントロピー変換の場合（即ち可逆で断念の場合、体積

V

での内部エネルギー

E_{int}

と圧力

P

は以下で与えられます：(19)

d E_{int} = - P d V

E_{int}

と

E

の関係を用いて以下が導かれます：(20)

d E = - \frac{P}{V_{0}} d V

μ

は体積比に関して以下のように表せます：(21)

μ = \frac{v_{0}}{v} - 1

対壁変化に対する関数の変化も以下のように表せます：(22)

d μ = - \frac{V_{0}}{V^{2}} d V = - \frac{{(1 + μ)}^{2}}{V_{0}} d V

単位体積当たりの内部エネルギー変化

E

はその時：(23)

d E = - \frac{P}{{(1 + μ)}^{2}} d μ

(24)

\frac{d P (μ, E)}{d μ} = {\frac{\partial P}{\partial μ} |}_{E} + \frac{P}{{(1 + μ)}^{2}} {\frac{\partial P}{\partial E} |}_{μ}

最後に、音速は以下のように与えられます：(25)

c^{2} = \frac{1}{ρ_{0}} {\frac{\partial P}{\partial μ} |}_{E} + \frac{P}{ρ_{0} {(1 + μ)}^{2}} {\frac{\partial P}{\partial E} |}_{μ}

この式を用いて、状態

P (μ, E)

の特定の方程式の音速を計算できます。理想気体の場合は、各タイプの定式化（絶対または相対）について、EOSを次のように記述できます：(26)

P (μ, E) = C_{0} + C_{1} μ + (C_{4} + C_{5} μ) E

式 25 は音速の計算に用いられます：(27)

{\frac{\partial P}{\partial μ} |}_{E} = C_{1} + C_{5} E

(28)

{\frac{\partial P}{\partial E} |}_{μ} = C_{4} + C_{5} μ

(29)

c^{2} = {\frac{C_{1} + C_{5} E}{ρ 0}}_{} + \frac{C_{0} + C_{1} μ + (C_{4} + C_{5} μ) E}{ρ_{0} {(1 + μ)}^{2}} (C_{4} + C_{5} μ)

この計算は次にそれぞれの4ケースに適用されます。

表 3. 数値解析結果の音速対理論表現
ケース	C₀	`C`₁	`C`₄	`C`₅	c² 以下より：式 25	理論値との比較
1	0	0	$γ - 1$	$γ - 1$	$\frac{γ (γ - 1) E}{ρ_{0}}$	$c = c^{R E F}$
2	0	0	$γ - 1$	$γ - 1$	$\frac{γ (γ - 1) E}{ρ_{0}}$	$c = c^{R E F}$
3	$E_{0} (γ - 1)$	$E_{0} (γ - 1)$	$γ - 1$	$γ - 1$	$\frac{γ (γ - 1) (E + E_{0})}{ρ_{0}}$	$c = c^{R E F}$
4	$E_{0} (γ - 1)$	$E_{0} (γ - 1)$	$γ - 1$	$γ - 1$	$\frac{γ (γ - 1) (E + E_{0})}{ρ_{0}}$	$c = c^{R E F}$

それぞれの4つの定式化に対して、Radiossにより計算された音速は理論値と一致しています。時間ステップとサイクル数は影響されません。

ケース 1：圧力とエネルギーの両方が絶対値

圧力：
状態方程式(30) $P = (γ - 1) \frac{E_{int}}{V} = (γ - 1) (1 + μ) \frac{E_{int}}{V_{0}}$

ここで、 ${E_{int} |}_{t = 0} = E_{0} V_{0} = \frac{P_{0} V_{0}}{γ - 1}$

多項式の係数を特定すると以下が導かれます：(31) $P = (C_{4} + C_{5} μ) E$

ここで、 $C_{4} = C {}_{5} = (γ - 1); E_{0} = \frac{P_{0}}{γ - 1}; P_{s h} = 0$

対応する入力：

#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/MAT/HYD_VISC/1
Polynomial EOS-Absolute Pressure-Absolute Energy
#              RHO_I               RHO_0
               1.204                   0
#                Knu                Pmin
                   0                   0
/EOS/POLYNOMIAL/1
Polynomial EOS-Absolute Pressure-Absolute Energy
#                 C0                  C1                  C2                  C3
                   0                   0                   0                   0
#                 C4                  C5                  E0                 Psh               RHO_0
                  .4                  .4              250000                   0               1.204
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|

結果の出力：

表 4.
時刻歴応答	寸法	初期値	Unit
/TH/BRICK ( $P$ )	$P$	$P_{0}$	圧力
/TH (IE)	$E_{int} (= E \cdot V_{0})$	$E_{0} V_{0}$	エネルギー
/TH/BRICK (IE)	$E_{int} / V$	$E_{0}$	圧力

理論界との比較
理想気体圧力の数値解析結果は、時刻歴で与えられます。要素時刻歴（/TH/BRICK）でこれを表示できます。この結果を理論解と比較します。曲線は重ねて示されます。

図 7. 数値的圧力、モデル1

内部エネルギーは2つの異なる方法で得ることができます。1つ目は、要素時刻歴（/TH/BRICK）に記録された内部エネルギー密度（ $E_{int} / V$ ）です。2つ目は、全体の時刻歴からの内部エネルギー $\sum_{e l e m e n t} E_{int}$ で、これは、モデルが単一要素であるためです。

図 8. 数値的内部エネルギー、モデル1

ケース 2：圧力が相対値でエネルギーが絶対値

圧力：
状態方程式(32) $P (μ, E) = (γ - 1) (1 + μ) \frac{E_{int}}{V_{0}}$

相対圧力：(33) $Δ P = P (μ, E) - P_{0} = (γ - 1) (1 + μ) \frac{E_{int}}{V_{0}} - P_{0}$

多項式の係数を特定すると以下が導かれます：(34) $Δ P = P = P_{s h} = - P_{s h} + (C_{4} + C_{5} μ) E$

ここで、 $C_{4} = C {}_{5} = (γ - 1); E_{0} = \frac{P_{0}}{γ - 1}; P_{s h} = 0$

最小圧力： (35) $P_{\min} = - P_{0}$

$P > 0 \Rightarrow Δ P \geq - P_{0}$ のため、最小圧力は0以外の値に設定される必要があります。

対応する入力：

#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/MAT/HYD_VISC/1
Polynomial EOS-Relative Pressure-Absolute Energy
#              RHO_I               RHO_0
               1.204                   0
#                Knu                Pmin
           1.5256E-5             -100000
/EOS/POLYNOMIAL/1
Polynomial EOS-Relative Pressure-Absolute Energy
#                 C0                  C1                  C2                  C3
                   0                   0                   0                   0
#                 C4                  C5                  E0                 Psh               RHO_0
                  .4                  .4              250000              100000               1.204
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|

結果の出力：

時刻歴応答	寸法	初期値	Unit
/TH/BRICK ( $P$ )	$Δ P$	0	圧力
/TH (IE)	$E_{int} (= E \cdot V_{0})$	$E_{0} V_{0}$	エネルギー
/TH/BRICK (IE)	$E_{int} / V$	$E_{0}$	圧力

理論界との比較
理想気体圧力の数値解析結果は、時刻歴で与えられます。要素時刻歴（/TH/BRICK）でこれを表示できます。この結果を理論解と比較します。曲線は重ねて示されます。

要素時刻歴（/TH/BRICK）は、P_shに対する相対圧力です。結果の曲線は、P_sh値でシフトされ、0から開始します。

図 9. 数値的圧力、モデル2

内部エネルギーは2つの異なる方法で得ることができます。1つ目は、要素時刻歴（/TH/BRICK）に記録された内部エネルギー密度（ $E_{int} / V$ ）です。2つ目は、全体の時刻歴からの内部エネルギー $\sum_{e l e m e n t} E_{int}$ で、これは、モデルが単一要素であるためです。

図 10. 数値的内部エネルギー、モデル2

ケース 3：圧力とエネルギーの両方が相対値

圧力：
状態方程式(36) $P = (γ - 1) (1 + μ) \frac{E_{int}}{V_{0}}$

初期内部エネルギーを導入することができます：(37) $E_{int} = E_{int} + ({E_{int} |}_{t = 0} - {E_{int} |}_{t = 0})$

参照値からの圧力は以下を与えます：(38) $P - P_{0} = Δ P = (γ - 1) (1 + μ) (Δ E + E_{0}) - P_{0}$

ここで、 $Δ E = \frac{E_{int} - {E_{int} |}_{t = 0}}{V_{0}}; E_{0} = \frac{{E_{int} |}_{t = 0}}{V_{0}}$ 。

多項式の係数を特定すると以下が導かれます：(39) $Δ P = P - P_{s h} = C_{0} + C_{1} μ + (C_{4} + C_{5} μ) Δ E - P_{s h}$

ここで、 $C_{0} = C_{1} = E_{0} (γ - 1)$ 、 $C_{4} = C_{5} = γ - 1$ 、 $Δ E_{0} = 0$ および $P_{s h} = P_{0}$ 。

最小圧力：(40) $P_{\min} = - P_{0}$

$P \geq 0 \Rightarrow Δ P \geq - P_{0}$ のため、最小圧力は0以外の値に設定される必要があります。

対応する入力：

#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/MAT/HYD_VISC/1
Polynomial EOS-Relative Pressure-Absolute Energy
#              RHO_I               RHO_0
               1.204                   0
#                Knu                Pmin
           1.5256E-5             -100000
/EOS/POLYNOMIAL/1
Polynomial EOS-Relative Pressure-Absolute Energy
#                 C0                  C1                  C2                  C3
              100000              100000                   0                   0
#                 C4                  C5                  E0                 Psh               RHO_0
                  .4                  .4                   0              100000               1.204
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|

結果の出力：

時刻歴応答	寸法	Unit
/TH/BRICK ( $P$ )	$Δ P$	圧力
/TH (IE)	$E_{int} (= Δ E \cdot V_{0})$	エネルギー
/TH/BRICK (IE)	$Δ E_{int} / V$	圧力

理論界との比較
理想気体圧力の数値解析結果は、時刻歴で与えられます。要素時刻歴（/TH/BRICK）でこれを表示できます。この結果を理論解と比較します。曲線は重ねて示されます。

要素時刻歴（/TH/BRICK）は、P_shに対する相対圧力です。結果の曲線は、P_sh値でシフトされ、0から開始します。

図 11. 数値的圧力、モデル3

内部エネルギーは2つの異なる方法で得ることができます。1つ目は、要素時刻歴（/TH/BRICK）に記録された内部エネルギー密度（ $E_{int} / V$ ）です。2つ目は、全体の時刻歴からの内部エネルギー $\sum_{e l e m e n t} E_{int}$ で、これは、モデルが単一要素であるためです。数値解析結果の内部エネルギーは初期エネルギーからの相対値で、絶対値から $E_{0} V_{0}$ 値だけシフトされ、0から開始します。

図 12. 数値的内部エネルギー、モデル3

ケース 4：圧力が絶対値でエネルギーが相対値

圧力：
状態方程式(41) $P = (γ - 1) (1 + μ) \frac{E_{int}}{V_{0}}$

初期内部エネルギーを導入することができます：(42) $E_{int} = E_{int} + ({E_{int} |}_{t = 0} - {E_{int} |}_{t = 0})$

与えられた参照値からの圧力： (43) $P = (γ - 1) (1 + μ) (Δ E + E_{0})$

多項式の係数を特定すると以下が導かれます：(44) $P = C_{0} + C_{1} μ + (C_{4} + C_{5} μ) Δ E$

ここで $C_{0} = C_{1} = E_{0} (γ - 1)$ ； $C_{4} = C_{5} = γ - 1$

対応する入力：

#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/MAT/HYD_VISC/1
Polynomial EOS-Relative Pressure-Absolute Energy
#              RHO_I               RHO_0
               1.204                   0
#                Knu                Pmin
           1.5256E-5                   0
/EOS/POLYNOMIAL/1
Polynomial EOS-Relative Pressure-Absolute Energy
#                 C0                  C1                  C2                  C3
              100000              100000                   0                   0
#                 C4                  C5                  E0                 Psh               RHO_0
                  .4                  .4                   0                   0               1.204
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|

結果の出力：

時刻歴応答	寸法	初期値	Unit
/TH/BRICK ( $P$ )	$P$	$P_{0}$	圧力
/TH (IE)	$E_{int} (= Δ E \cdot V_{0})$	0	エネルギー
/TH/BRICK (IE)	$Δ E_{int} / V$	0	圧力

理論界との比較
要素時刻歴（/TH/BRICK）は、絶対圧力を与えます。この結果を理論解と比較します。曲線は重ねて示されます。

図 13. 数値的圧力、モデル4

内部エネルギーは2つの異なる方法で得ることができます。1つ目は、要素時刻歴（/TH/BRICK）に記録された内部エネルギー密度（ $Δ E_{int} / V$ ）です。2つ目は、全体の時刻歴からの内部エネルギー $\sum_{e l e m e n t} E_{int}$ で、これは、モデルが単一要素であるためです。数値解析結果の内部エネルギーは初期エネルギーからの相対値で、絶対値から $E_{0} V_{0}$ 値だけシフトされ、0から開始します。

図 14. 数値的内部エネルギー、モデル4