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

理想気体のモデル化に多項式状態方程式が用いられます。圧力またはエネルギーは、絶対値または相対値であることが可能です。材料則LAW6（/MAT/HYDRO）がこれらのケースの材料カードの構築に用いられます。

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

多項式EOSはRadioss流体力学的圧力を計算するためにしばしば用いられます。これは圧縮で3次、膨張で線形です。(1) ex43_polynomial_eq

ここで、(2)

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

および (3)

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

材料LAW6（/MAT/HYDRO）は静水圧の計算にこの式を使用します。絶対値または相対変化を考慮することが可能です（表 1）。この例題では以下のケースのそれぞれののために材料コントロールカードをどのように構築するかを紹介します：

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

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

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

理想気体
多項式 EOS（/EOS/POLYNOMIAL）
絶対 / 相対定式化
圧力シフト
流体材料（/MAT/LAW6 (HYDROまたはHYD_VISC))
強制変位（/IMPDISP）
それぞれの面の節点は強制変位で移動されます
境界条件（/ALE/BCS）
境界節点はLagrangianとして定義されます

要素圧力、密度と内部エネルギー密度は、時刻歴ファイルに保存されます。

入力ファイル

本例題で使用される入力ファイルは下記のとおり：

Model 1: <install_directory>/hwsolvers/demos/radioss/example/43_perfect_gas_polynomial_eos/01-Pabsolute_Eabsolute/*
Model 2: <install_directory>/hwsolvers/demos/radioss/example/43_perfect_gas_polynomial_eos/02-Prelative_Eabsolute/*
Model 3: <install_directory>/hwsolvers/demos/radioss/example/43_perfect_gas_polynomial_eos/03-Prelative_Erelative/*
Model 4: <install_directory>/hwsolvers/demos/radioss/example/43_perfect_gas_polynomial_eos/04-Pabsolute_Erelative/*

モデル概要

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

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

$P_{0}$: 1e5 Pa
$V_{0}$: 1000 m³
$ρ_{0}$: 1.204 $[\frac{kg}{m^{3}}]$
$μ_{0}$: 0

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

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

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

多項式 EOS

多項式EOSが材料則6（/MAT/HYDRO）の中で流体力学的圧力の計算に用いられます。これは圧縮で3次、膨張で線形です。(4)
ex43_polynomial_eq

ここで、

$P$: 流体力学的圧力

(5)

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

および(6)

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

{C_{i}} i = 0.5

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

一般的なケースはMie-Guneisen EOSに一致（Theory ManualのAppendix B: Basic Relations of Elasticityを参照のこと）
非圧縮性気体
線形弾性材料
理想気体

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

モデリングの方法

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

材料則LAW6（/MAT/HYDRO）が用いられ、流体力学的粘性流体材料を記述します。

(1)	(3)	(5)	(7)
/MAT/LAW6/`mat_ID`/`unit_ID`または/MAT/HYDRO/`mat_ID`/`unit_ID`
`mat_title`
$ρ_{i}$
$υ$
`C`₀	`C`₁	`C`₂	`C`₃
`P`_min	`P`_sh
`C`₄	`C`₅	`E`₀

圧力シフト

材料則LAW6ではフラグP_shを導入し、多項式状態方程式の圧力計算シフトを可能にしています：(7)
ex43_pressure_shift

Radioss Engineは、C₀フラグと計算された圧力 $P (μ, E)$ を-P_shのオフセットでシフトします。

最小圧力

(8)

P_{\min} = \lim_{μ \to - 1} P (μ, E) - P_{s h}

理論解は $P_{m i n} = 0 P a$ （絶対圧力）で、デフォルト値を-10³⁰とすることで、相対圧力定式化の負の値を可能にしています。

フラグは-P_shで手動にてオフセットする必要があります。

結果

理論解

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

V

または

μ

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

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

したがって：(10) $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は内部エネルギー変化の計算に等エントロピー過程の仮説を仮定します：(11) $d E_{int} = - P d V$

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

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

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

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

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

理論解を下の表に記します。圧力、内部エネルギーと音速は

V

と

μ

の両方で表現されます。

圧力（Pa）		内部エネルギー密度（J）		音速（m/s）
P^REF(V)	P^REF( $μ$ )	$ρ$ e^REF(V)	$ρ$ e^REF( $μ$ )	c^REF(V)	c^REF( $μ$ )
$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}}$

対応するプロットは：

材料コントロールカード

材料は理想気体と想定します。以下のケースが調査されます：

ケース1：圧力とエネルギーの両方が絶対値： $P (μ, E)$
ケース2：圧力が相対値でエネルギーが絶対値： $Δ P (μ, E)$
ケース3：圧力とエネルギーの両方が相対値 $Δ P (μ, Δ E)$
ケース4：圧力が絶対値でエネルギーが相対値： $P (μ, Δ E)$

音速と時間ステップ

材料則6は流体の通常表現で音速を計算します：(19)

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

これは

μ

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

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

したがって：(21)

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

P

の内部エネルギー

E

と

μ

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

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

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

V

での内部エネルギー

E_{int}

と圧力

P

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

d E_{int} = - P d V

E_{int}

と

E

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

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

μ

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

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

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

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

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

E

はその時：(27)

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

(28)

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

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

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

この式を用いて、状態

P (μ, E)

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

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

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

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

(32)

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

(33)

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ケースに適用されます。

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

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

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

状態方程式
状態方程式は以下のように書くことができます：(34) $P = (γ - 1) \frac{E_{int}}{V} = (γ - 1) (1 + μ) \frac{E_{int}}{V_{0}}$

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

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

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

対応する入力

(1)	(3)	(5)	(7)
/MAT/LAW6/`mat_ID`/`unit_ID`または/MAT/HYDRO/`mat_ID`/`unit_ID`
`AbsolutePRESSURE_AbsoluteENERGY`
$ρ_{i}$
$υ$
0	0	0	0
0	0
`C`₄ = $γ - 1$ 1	`C`₅ = $γ - 1$	$E_{0} = \frac{P_{0}}{γ - 1}$

結果の出力

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

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

図 7. 数値的圧力、モデル1： $P (μ, E)$

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

図 8. 数値的内部エネルギー、モデル1： $P (μ, E)$

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

状態方程式
理想気体の状態方程式は：(37) $P (μ, E) = (γ - 1) (1 + μ) \frac{E_{int}}{V_{0}}$

参照値からの圧力計算は相対圧力を与えます：(38) $Δ P = P (μ, E) - P_{0} = (γ - 1) (1 + μ) \frac{E_{int}}{V_{0}} - P_{0}$

この式を展開して多項式の係数を特定すると以下が導かれます：(39) $Δ P (μ, E) = P (μ, E) = 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$
最小圧力(40) $P (μ, E) \geq 0 \Rightarrow Δ P \geq - P_{0}$
最小圧力は0でない値 $P_{\min} = - P_{0}$ に設定される必要があります。

対応する入力

(1)	(3)	(5)	(7)
/MAT/LAW6/`mat_ID`/`unit_ID`または/MAT/HYDRO/`mat_ID`/`unit_ID`
`RelativePRESSURE_AbsoluteENERGY`
$ρ_{i}$
$υ$
0	0	0	0
`-P`₀	`P`₀
`C`₄ = $γ - 1$	`C`₅ = $γ - 1$	$E_{0} = \frac{P_{0}}{γ - 1}$

結果の出力

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

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

図 9. 数値的圧力、モデル2： $Δ P (μ, E)$

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

図 10. 数値的内部エネルギー、モデル2： $Δ P (μ, E)$

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

状態方程式
理想気体の状態方程式は：(41) $P (μ, E) = (γ - 1) (1 + μ) \frac{E_{int}}{V_{0}}$

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

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

ここで、(44) $Δ E = \frac{E_{int} - E_{int} {}_{| t = 0}}{V_{0}} E_{0} = {\frac{E_{int} {}_{| t = 0}}{V_{0}}}_{}$

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

ここで、 $C_{0} = C_{1} = E_{0} (γ - 1)$ $C_{4} = C_{5} = γ - 1$ $Δ E_{0} = 0$ $P_{s h} = P_{0}$
最小圧力(46) $P (μ, E) \geq 0 \Rightarrow Δ P \geq - P_{0}$
最小圧力は0でない値 $P_{\min} = - P_{0}$ に設定される必要があります。

対応する入力

(1)	(3)	(5)	(7)
/MAT/LAW6/`mat_ID`/`unit_ID`または/MAT/HYDRO/`mat_ID`/`unit_ID`
`RelativePRESSURE_RelativeENERGY`
$ρ_{i}$
$υ$
$E_{0} (γ - 1)$	$E_{0} (γ - 1)$	0	0
-P₀	P₀
`C`₄ = $γ - 1$	`C`₅ = $γ - 1$	0

結果の出力

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

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

図 11. 数値的圧力、モデル3： $Δ P (μ, Δ E)$

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

図 12. 数値的内部エネルギー、モデル3： $Δ P (μ, Δ E)$

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

状態方程式
理想気体の状態方程式は：(47) $P (μ, E) = (γ - 1) (1 + μ) \frac{E_{int}}{V_{0}}$

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

ここに、以下が導かれます：(49) $P (μ, E) = (γ - 1) (1 + μ) (E_{0} + Δ E)$

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

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

対応する入力

(1)	(3)	(5)	(7)
/MAT/LAW6/`mat_ID`/`unit_ID`または/MAT/HYDRO/`mat_ID`/`unit_ID`
`AbsolutePRESSURE_RelativeENERGY`
$ρ_{i}$
$υ$
$E_{0} (γ - 1)$	$E_{0} (γ - 1)$	0	0
0	0
`C`₄ = $γ - 1$	`C`₅ = $γ - 1$	0

結果の出力

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

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

図 13. 数値的圧力、モデル4： $P (μ, Δ E)$

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

図 14. 数値的内部エネルギー、モデル4： $P (μ, Δ E)$

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

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

入力ファイル

モデル概要

多項式 EOS

モデリングの方法

圧力シフト

最小圧力

結果

理論解

材料コントロールカード

音速と時間ステップ

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

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

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

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

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

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

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

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