NASA-CEA 利用方法
三好 明
三好 明
CEA は NASA Glenn Research Center の S. Gordon と B. J. McBride らによって開発された化学平衡計算プログラムです. 平衡計算プログラム本体 (CEA2) 以外に, 熱力学データライブラリや, 熱力学データの回帰や推定を行う PAC99 プログラムなどが含まれています.
problem tp p(atm) = 1 3 t(k) = 300 350 900 2000 reactant name = CH4 moles = 0.11 name = O2 moles = 0.21 name = N2 moles = 0.79 end
tp
)
平衡状態を計算します.
圧力 1, 3 atm, 温度 300, 350, 900, 2000 K の8つの組合せについて
計算を行います.
prob[lem]
, reac[tant]
はそれぞれ, 反応物を除く計算条件, 反応物 に関するブロックの開始を示します.
計算は理想気体を仮定して行われるので, 条件には一定とする2つの状態量を
記号で与えます. この例では温度 (t
) と圧力 (p
)
を一定としていますが, 他に hp
(エンタルピー・圧力一定 = 定圧断熱),
uv
(内部エネルギー・体積一定 = 定容断熱) などを指定することが
できます.
moles
) もしくは, 質量比 (wt
)
で与えることができます. モル分率や重量 % のように合計が 1 あるいは 100 で
ある必要はありません.
"THERMODYNAMIC EQUILIBRIUM ... PROPERTIES AT
ASSIGNED"
と書かれた行以降に出力されています.
主要部分の抜粋を以下に示します.THERMODYNAMIC EQUILIBRIUM PROPERTIES AT ASSIGNED TEMPERATURE AND PRESSURE CASE = REACTANT MOLES ENERGY TEMP KJ/KG-MOL K NAME CH4 0.1100000 0.000 0.000 NAME O2 0.2100000 0.000 0.000 NAME N2 0.7900000 0.000 0.000 O/F= 0.00000 %FUEL= 0.000000 R,EQ.RATIO= 1.047619 PHI,EQ.RATIO= 0.000000 THERMODYNAMIC PROPERTIES P, BAR 1.0132 1.0132 1.0132 1.0132 3.0397 3.0397 3.0397 3.0397 T, K 300.00 350.00 900.00 2000.00 300.00 350.00 900.00 2000.00 RHO, KG/CU M 1.3484 0 9.6032-1 3.7013-1 1.6648-1 4.1432 0 3.0809 0 1.1104 0 4.9956-1 H, KJ/KG -3290.30 -2962.43 -2263.54 -683.63 -3321.66 -3061.66 -2263.55 -688.02 U, KJ/KG -3365.44 -3067.94 -2537.30 -1292.25 -3395.03 -3160.33 -2537.31 -1296.51 G, KJ/KG -5210.59 -5559.99 -10010.3 -20168.8 -5129.20 -5445.26 -9709.59 -19500.2 S, KJ/(KG)(K) 6.4010 7.4216 8.6076 9.7426 6.0251 6.8103 8.2734 9.4061 M, (1/n) 33.194 27.581 27.335 27.323 33.998 29.495 27.335 27.328 MW, MOL WT 27.339 27.581 27.335 27.323 27.335 27.503 27.335 27.328 (dLV/dLP)t -1.03649 -1.00001 -1.00000 -1.00026 -1.01180 -1.17474 -1.00001 -1.00013 (dLV/dLT)p 1.6437 1.0002 1.0000 1.0090 1.2079 3.5514 1.0001 1.0046 Cp, KJ/(KG)(K) 4.2185 1.1279 1.3161 1.6012 2.2991 11.7487 1.3162 1.5554 GAMMAs 1.1415 1.3649 1.3006 1.2395 1.1674 1.1466 1.3006 1.2458 SON VEL,M/SEC 292.9 379.5 596.7 868.6 292.7 336.4 596.7 870.7 MOLE FRACTIONS CH4 0.00007 0.00450 0.00000 0.00000 0.00001 0.00307 0.00000 0.00000 *CO 0.00000 0.00000 0.00330 0.01235 0.00000 0.00000 0.00330 0.01220 *CO2 0.08937 0.09460 0.09491 0.08582 0.08929 0.09292 0.09491 0.08599 *H 0.00000 0.00000 0.00000 0.00013 0.00000 0.00000 0.00000 0.00007 *H2 0.00000 0.00003 0.01455 0.00597 0.00000 0.00001 0.01454 0.00589 H2O 0.02872 0.18917 0.18188 0.19008 0.00935 0.12677 0.18188 0.19032 NH3 0.00000 0.00001 0.00000 0.00000 0.00000 0.00000 0.00001 0.00000 *NO 0.00000 0.00000 0.00000 0.00014 0.00000 0.00000 0.00000 0.00008 *N2 0.70545 0.71170 0.70536 0.70497 0.70537 0.70969 0.70536 0.70515 *O 0.00000 0.00000 0.00000 0.00001 0.00000 0.00000 0.00000 0.00000 *OH 0.00000 0.00000 0.00000 0.00045 0.00000 0.00000 0.00000 0.00026 *O2 0.00000 0.00000 0.00000 0.00008 0.00000 0.00000 0.00000 0.00003 C(gr) 0.00879 0.00000 0.00000 0.00000 0.00891 0.00283 0.00000 0.00000 H2O(L) 0.16760 0.00000 0.00000 0.00000 0.18707 0.06470 0.00000 0.00000 * THERMODYNAMIC PROPERTIES FITTED TO 20000.K PRODUCTS WHICH WERE CONSIDERED BUT WHOSE MOLE FRACTIONS WERE LESS THAN 5.000000E-06 FOR ALL ASSIGNED CONDITIONS *C *CH CH2 CH3 CH2OH CH3O CH3OH CH3OOH *CN CNN ...
RHO, KG/CU M
)
にある "9.6032-1"
は限られた桁数に多くの有効数字を出力するための
CEA 特有の書式で "9.6032E-01"
(9.6032 × 10−1) と同じです.
( )
(括弧)
のついていないものは気相の化学種で H2O(L) は水 C(gr) はグラファイトを表しています.
正確な情報は CEA をインストールしたフォルダにある thermo.inp
ファイルで確認することができます.
reac
データセットは
以下のように指定しても (この場合, 推奨される入力ではありませんが...)
計算結果は同じです.name = C(gr) moles = 0.11 name = H2 moles = 0.22 name = NO moles = 0.42 name = N2 moles = 0.58
"Air"
で指定することもできます.
problem tp p(atm) = 1 3 t(k) = 300 350 900 2000 reactant fuel = CH4 moles = 0.11 oxid = Air moles = 1.00 end
name
で指定しましたが,
この例のように fuel
と oxid[ant]
で指定することも可能です. name
による指定と
計算結果は変わりませんが, 以下の燃料と酸化剤の比率に関する出力のみに影響します.
O/F= 16.41389 %FUEL= 5.742541 R,EQ.RATIO= 1.050163 PHI,EQ.RATIO= 1.050240(例 1 の結果では
O/F
, %FUEL
, PHI,EQ.RATIO
が 0 になっているのと比較して下さい)
prob tp phi(eq.ratio) = 1.05 p(atm) = 1 3 t(k) = 300 350 900 2000 reac fuel = CH4 oxid = Air end
fuel
と oxid
で与えておくと phi(eq.ratio)
を使って, この例のように, 指定した当量比になるように燃料と酸化剤の比を
与えることができます. 燃料と酸化剤の比率に関する出力は以下のようになります.
O/F= 16.41764 %FUEL= 5.741306 R,EQ.RATIO= 1.049924 PHI,EQ.RATIO= 1.050000
prob tp phi(eq.ratio) = 1.05 p(atm) = 1 3 t(k) = 300 350 900 2000 reac fuel = CH4 oxid = Air output trace = 1e-10 end
outp[ut]
キーワードで始まる
output dataset の指定の例を示しています. trace
は平衡計算の結果でモル分率が, ここで指定した値のものまで出力することを指示します.
(デフォルトは 5E-06
= 5 × 10−6)
ここでは 1E-10
を指定していますので, 指定した温度・圧力のいずれかで
1 × 10−10 以上のものはすべて出力されます.
以下のこの場合の出力の一部を示します. 出力の形式も上の 例1 の出力で説明した
CEA 独特の形式に変わります.
MOLE FRACTIONS *Ar 8.3594-3 8.4370-3 8.3583-3 8.3546-3 8.3584-3 8.4097-3 8.3583-3 8.3563-3 CH4 6.857 -5 4.709 -3 5.288 -8 8.665-16 7.266 -6 3.077 -3 4.752 -7 7.427-15 *CO 1.022-11 1.496 -9 3.473 -3 1.289 -2 5.828-12 9.006-10 3.472 -3 1.275 -2 *CO2 8.9170-2 9.4656-2 9.4964-2 8.5502-2 8.9098-2 9.2724-2 9.4965-2 8.5668-2 COOH 0.000 0 2.070-32 1.089-14 3.007 -9 0.000 0 2.337-32 1.886-14 5.188 -9 *H 0.000 0 5.730-33 1.434-11 1.279 -4 0.000 0 2.052-33 8.275-12 7.338 -5
prob tp phi(eq.ratio) = 1.05 p(atm) = 1 3 t(k) = 300 350 900 2000 reac fuel = CH4 oxid = Air outp trace = 1e-10 plot p t C(gr) H2O(L) H2O end
plot
以降に並べた変数の値が, 拡張子が plt
のファイル
(この例の場合は CH4air_tp_d.plt) に出力されます.
変数名 p
, t
はそれぞれ, 圧力, 温度で,
化学種の名前を指定すると, そのモル分率が出力されます.
以下に CH4air_tp_d.plt の内容を示します.
この形式のテキストファイルは Excel などで読み込むことができます.
コマンドプロンプト用の改変版 (FCEA2m.exe) を利用している場合は
拡張子が csv
のファイルも生成します.
このファイルは直接 Excel で開くことができます.
# p t C(gr) H2O(L) H2O 1.0132E+00 3.0000E+02 9.2118E-03 1.6748E-01 2.8712E-02 1.0132E+00 3.5000E+02 0.0000E+00 0.0000E+00 1.8870E-01 1.0132E+00 9.0000E+02 0.0000E+00 0.0000E+00 1.8108E-01 1.0132E+00 2.0000E+03 0.0000E+00 0.0000E+00 1.8969E-01 3.0397E+00 3.0000E+02 9.3334E-03 1.8695E-01 9.3430E-03 3.0397E+00 3.5000E+02 3.2418E-03 6.4611E-02 1.2673E-01 3.0397E+00 9.0000E+02 0.0000E+00 0.0000E+00 1.8108E-01 3.0397E+00 2.0000E+03 0.0000E+00 0.0000E+00 1.8993E-01 # p t C(gr) H2O(L) H2O
prob
dataset に hp
(エンタルピー・圧力一定 =
定圧断熱) を指定した例を示しています. 燃料と空気(酸素)を反応物に指定した場合,
これは断熱火炎の計算に相当します.
prob hp phi(eq.ratio) = 1.05 p(atm) = 1 3 reac fuel = CH4 t(k)= 298.15 oxid = Air t(k)= 298.15 end
prob
dataset
に hp
を指定したときは, reac
dataset
に燃焼前の気体の温度を与えなければなりません. 出力の一部を以下に示します.
断熱火炎温度は 1 atm で 2231 K, 3 atm で 2249 K であることを示しています.
THERMODYNAMIC PROPERTIES P, BAR 1.0132 3.0397 T, K 2230.98 2249.09 RHO, KG/CU M 1.4929-1 4.4487-1 H, KJ/KG -271.06 -271.06
problem hp p(atm) = 1 phi(eq.ratio) = 1 reac fuel= C3H6,propylene t(k)= 298.15 oxid= Air t(k)= 298.15 end
C3H6,propylene
としなければなりません.
シクロプロパンは C3H6,cyclo-
という名称で登録されています.
THERMODYNAMIC PROPERTIES P, BAR 1.0132 T, K 2332.68 RHO, KG/CU M 1.4913-1 H, KJ/KG 26.047 U, KJ/KG -653.41 G, KJ/KG -22336.7 S, KJ/(KG)(K) 9.5867 M, (1/n) 28.545 (dLV/dLP)t -1.00405 (dLV/dLT)p 1.1160 Cp, KJ/(KG)(K) 2.4353 GAMMAs 1.1695 SON VEL,M/SEC 891.4 MOLE FRACTIONS *Ar 0.00864 *CO 0.01776 *CO2 0.11144 *H 0.00067 *H2 0.00344 H2O 0.12291 *NO 0.00310 *N2 0.71922 *O 0.00053 *OH 0.00446 *O2 0.00784
prob hp p(atm) = 1 reac name = C2H2,acetylene t(k)= 298.15 end
THERMODYNAMIC PROPERTIES P, BAR 1.0132 T, K 2849.07 RHO, KG/CU M 1.0633-1 : MOLE FRACTIONS : *H 0.03179 *H2 0.30807 C(gr) 0.64748
prob
dataset に uv
(内部エネルギー・体積一定 =
定容断熱) を指定した例を示しています. 燃料と空気(酸素)を反応物に指定した場合,
定容容器中の断熱既燃平衡状態を計算することになります.
体積膨張がないので, (定圧)断熱火炎温度より高温になります.
prob uv phi(eq.ratio) = 1.05 rho(kg/m**3) = 1.1316 reac fuel = CH4 t(k) = 298.15 oxid = Air t(k) = 298.15 end
prob
dataset
に uv
を指定したときは, 密度 または 比体積 (単位質量あたりの体積)
を与える必要があります. ここでは rho
で密度を与えています.
出力の一部を以下に示します.
(定圧)断熱火炎よりも高温になっていることがわかります.
THERMODYNAMIC PROPERTIES P, BAR 9.0164 T, K 2601.97 RHO, KG/CU M 1.1316 0 H, KJ/KG 436.17 U, KJ/KG -360.61入力に必要な密度は, もちろん 理想気体の状態方程式から計算できるわけですが, 複雑な混合気では計算が大変になるので, 例えば, 下に示すような (意味のない) 計算をしてその値を得ることもできます.
prob tp phi(eq.ratio) = 1.05 p(atm) = 1 t(k) = 298.15 reac fuel = CH4 oxid = Air only CH4 N2 O2 Ar CO2 end
tp
) ですが, only
キーワードを使って, 組成が変わらないようにしています.
ここでは, 以下の出力にある密度 (ρ = 1.1316×100 kg/m3)
を上の定容断熱計算の入力に使用しました.
THERMODYNAMIC PROPERTIES P, BAR 1.0132 T, K 298.15 RHO, KG/CU M 1.1316 0 H, KJ/KG -271.07 U, KJ/KG -360.61 G, KJ/KG -2426.74 S, KJ/(KG)(K) 7.2302この出力は, 以下の断熱圧縮の計算でも使用します.
prob
dataset に sv
および sp
を指定した例を示します. 等エントロピー過程である, 断熱圧縮・膨張
などを計算することができます. 多くの場合, 圧縮比 (体積比)
が既知なので sv
を使いますが, 急速圧縮機の実験で用いられる
コアガスモデル の計算では, 観測圧力から温度を計算するので sp
を使います.
prob sv phi(eq.ratio) = 1.05 s/r = 0.869588 rho(kg/m**3) = 1.1316 3.3948 11.316 reac fuel = CH4 oxid = Air only CH4 N2 O2 Ar CO2 end
prob
dataset に sv
を指定したときは,
比エントロピー (単位質量あたりのエントロピー) と 密度
(または 比体積) を与える必要があります.
ここでは上の 「意味のない」 計算の結果を使っています.
比エントロピーは気体定数で割った値
7.2302/8.31451 = 0.869588 を入力します.
この計算では未反応の状態でメタン-空気混合気を, 初期体積,
体積 1/3 (密度 3 倍), 体積 1/10 (密度 10 倍)
に断熱圧縮した状態を計算しています.
以下の結果に示すように, 比エントロピーは確かに保存されています.
THERMODYNAMIC PROPERTIES P, BAR 1.0133 4.6099 23.364 T, K 298.16 452.15 687.46 RHO, KG/CU M 1.1316 0 3.3948 0 1.1316 1 H, KJ/KG -271.05 -102.54 172.01 U, KJ/KG -360.60 -238.34 -34.453 G, KJ/KG -2426.84 -3371.70 -4798.48 S, KJ/(KG)(K) 7.2302 7.2302 7.2302
prob sp phi(eq.ratio) = 1.05 s/r = 0.869588 p(bar) = 1.01325 10 23 reac fuel = CH4 oxid = Air only CH4 N2 O2 Ar CO2 end
prob
dataset に sp
を指定したときは,
比エントロピーと圧力を与える必要があります.
上と同様, 比エントロピーは気体定数で割った値を入力しています.
この計算では未反応の状態でメタン-空気混合気を, 初期圧,
10 bar, 23 bar まで断熱圧縮した状態を計算しています.
THERMODYNAMIC PROPERTIES P, BAR 1.0132 10.000 23.000 T, K 298.16 554.67 684.81 RHO, KG/CU M 1.1315 0 6.0030 0 1.1183 1 H, KJ/KG -271.06 14.197 168.78 U, KJ/KG -360.60 -152.39 -36.887 G, KJ/KG -2426.81 -3996.18 -4782.51 S, KJ/(KG)(K) 7.2302 7.2302 7.2302
only
の入力を取った計算は, 決して,
断熱圧縮気体が燃焼した平衡状態の計算にはならないことに注意して下さい.
そのような計算をしたい場合は, 圧縮状態を初期条件として uv
の計算をします. 例えば, 上の体積 1/10 圧縮状態を初期条件とする計算は
以下のように入力することになります.
prob uv phi(eq.ratio) = 1.05 rho(kg/m**3) = 11.316 reac fuel = CH4 t(k) = 687.46 oxid = Air t(k) = 687.46 end
! CH4 steam reforming problem tp p(atm) = 1 t(k) = 400 450 500 550 600 650 700 750 800 850 900 950 1000 1050 1100 1150 1200 reac name = CH4 moles = 1 name = H2O moles = 1 omit C(gr) outp plot p t CH4 H2O CO H2 CO2 end
problem tp
で計算する場合は, 少量の Ar を初期組成に加えるのが ミソ です.
これがないと低温では, 指定された圧力の解が存在しないために計算は破綻します. また全圧を 10 atm
としているのも, 高温で解がなくなることを避けるためです. 1 atm にしていると
平衡蒸気圧が 1 atm 以上になる温度では計算が破綻します.
! decomposition of CaCO3(cr) problem tp p(atm) = 10 t(k) = 700 750 800 850 900 950 1000 1050 1100 1150 1200 1250 1300 reac name = CaCO3(cr) moles = .9 name = Ar moles = .1 outp plot t p Ar CO2 CaCO3(cr) CaO(cr) end
thermo.inp
を編集することで容易に計算することができます.
thermo.inp
の末尾を以下のように編集して END REACTANTS
の行の前に reactant (反応物) としてフランを追加します.... RP-1 Mehta et.al. AIAA 95-2962 1995. Hcomb(high) = 19923.BTU/# 0 gll/00 C 1.00H 1.95 0.00 0.00 0.00 1 13.9761830 -24717.700 298.150 0.0000 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.000 C4H4O,furan Furan 0 am1209 C 4.00H 4.00O 1.00 0.00 0.00 0 68.0739600 -34700.000 298.150 0.0000 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.000 END REACTANTS
am1209
はオプションの id で空白でもかまいません.
その後ろに化学式を [元素記号, その数] を繰り返して入力します.
68.0739600
の部分には分子量を書きます. 2 行目の最後は
標準生成エンタルピーで
NIST Chemistry WebBook
から得た, フラン(気体)の値 (−34.7 kJ mol−1)
を 入力してあります. 単位は J mol−1 です.
3 行目は, 通常, この例の通りで問題ありません.
いずれの入力欄も FORTRAN の書式付入力で読み込まれますので, 桁位置に注意して下さい.
fcea2.exe
または fcea2m.exe
を起動して thermo
(enter)
と入力します. この名称は特別な予約名称で, この操作により CEA は新しい thermo.lib
を生成します.
以後はこの thermo.lib
をカレントディレクトリに置いて計算して下さい.
problem case=furan-air hp p(atm)=1 phi,eq.ratio=1 reac fuel= C4H4O,furan wt%=100 t(k)= 298.15 oxid= Air wt%=100 t(k)= 298.15 end
thermo.lib
に登録されたら,
後は登録済みの場合と全く同様に計算することができます.