朧モデルの差分化方法

杉山耕一朗 (北大理・地惑) sugiyama@gfd-dennou.org

小高正嗣 (東大数理) 倉本圭 (北大理・地惑) 林祥介 (北大理・地惑)

平成 15 年 6 月 22 日

\underline{朧モデルの差分化方法 1}

目次

第 1 章大気の平衡状態の記述 2

1.1 ギブス自由エネルギーの極値の与え方 . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.1.1 大気が理想気体のみから成る場合 . . . . . . . . . . . . . . . . . . . . . . . 2 1.1.2 多相系への拡張 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.1.3 負のモル数になった場合には . . . . . . . . . . . . . . . . . . . . . . . . . 7

1.2 RAND 法の適用できないケースとその解決方法 . . . . . . . . . . . . . . . . . . . 8 1.2.1 ある 2 つの相に存在する物質が全く同じな場合 . . . . . . . . . . . . . . . 8 1.2.2 ある相に含まれる物質の存在量が 0 に近付く場合 . . . . . . . . . . . . . . 9

1.3 反復計算手法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 第 2 章大気の断熱変化の記述 14 第 3 章化学ポテンシャル, エントロピーの計算方法 17

3.1 比熱を用いた標準化学ポテンシャル, 標準モルエントロピーの計算 . . . . . . . . . 17 3.1.1 希ガス . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.1.2 スプライン補間 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

3.2 飽和蒸気圧を用いた標準化学ポテンシャル, 標準モルエントロピーの計算 . . . . . 19 3.2.1 Antoine 式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.2.2 AMP 式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

3.3 化学反応の圧平衡定数を用いた標準化学ポテンシャル, 標準モルエントロピーの計算 20 第 4 章雲密度の計算方法 22 第 5 章静的安定度の計算 23 第 6 章参考文献 24 付録 A AMP 式 25 付録 B 相平衡条件を用いた存在可能な化学種の判定 27

B.1 気相と凝縮相の相平衡条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 B.2 液相と固相の相平衡条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 B.3 気相と溶液の相平衡: 溶液について勉強しなおし! . . . . . . . . . . . . . . . . . . 28 B.4 化学反応の平衡 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

付録 C 初期値の作成方法 29 C.1 使わなかった方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

\underline{朧モデルの差分化方法第 1 章大気の平衡状態の記述 2}

第1章大気の平衡状態の記述

大気の平衡状態を熱力学関数の最小化された状態として記述する. 大気の状態量として温度・圧力・組成を選択すると, 大気の平衡状態は大気中の元素数保存の下でギブス自由エネルギー G の最小化された状態として記述することができる.

大気の温度・圧力を与えた場合, 平衡状態は条件付き極値問題として扱うことが可能となる. すなわち, 元素数保存の条件,

∑_{i} a_{ie}n_{i} = B_{e}.  (1.1)
n_{i} ≥ 0  (1.2)

の下, ギブス自由エネルギー G,

G = ∑ μ_{i}n_{i}.  (1.3)

の最小値を求めれば良い. ここで添字 i は大気に含まれる i 番目 (0 ≤ i ≤ N) の物質を意味し, 添字 e は系に含まれる e 番目 (0 ≤ e ≤ M) の元素を意味する. μ_{i} は i 番目の物質の化学ポテンシャルを意味する. また a_{ie} は i 番目の物質に含まれる e 番目の元素の数, B_{e} は系全体に含まれる e 番目の元素の個数を意味する. 例えば元素 H (e = 1), O (e = 2) と物質 H_{2}O (i = 1) を考えると a_{1,1} = 2, a_{1,2} = 1 となる. (1.2) 式はモル数は常に正の数であることを意味する.

本章では数値的に (1.3) 式の条件付き極値を求める方法を議論する. この計算方法の骨子は, 逐次近似法によってギブス自由エネルギー G が最小となる物質存在量を求めることにある. ギブス自由エネルギー G を適当な物質存在量のまわりで 2 次の関数に近似し, その最小値を計算する. 得られた物質存在量が収束するまで同様の操作を繰り返す.

まず第 1.1 節において, ラグランジュの未定乗数法によるギブス自由エネルギーの最小化法を議論する. そして第 1.3 節でギブスの自由エネルギーを反復的に収束させる方法を議論する.

1.1 ギブス自由エネルギーの極値の与え方

1.1.1 大気が理想気体のみから成る場合

大気が理想気体のみから成る場合を考える. この場合, White et al. (1958) らによって開発された「RAND 法」を用いる (Van Zeggerren and Storey, 1970).

適当な物質存在量を m_{i} とする. m_{i} よりも熱化学平衡に近い物質存在度を n_{i}(= m_{i} + δm_{i}) とする. 平衡時の熱力学ポテンシャル G を m_{i} のまわりでテーラー展開し, δm_{i} に関する 2 次の微小量まで近似式 Q(n_{i}) を作る.

Q(n_{i}) = G(m_{i}) + ∑_{i} (\frac{∂G}{∂m_{i}} ) δm_{i} + \frac{1}{2} ∑_{i} ∑_{j} (\frac{∂^{2}G}{∂m_{i}∂m_{l}} ) δm_{i}δm_{j}.  (1.4)

ここで理想気体の化学ポテンシャルの定義式を用いると, 右辺第 2 項のカッコ内は以下のように書ける.

(\frac{∂G}{∂m_{i}} ) = μ_{i} = μ_{i}^{∘}(T) + RT ln (\frac{m_{i}}{∑ m_{i}} ) + RT ln (\frac{p}{p_{0}} )  (1.5)

但し p_{0} は基準圧力で通常 1.0 × 10^{5} Pa である. (1.5) 式を用いて (1.4) 式の右辺第 3 項のカッコ内は以下のように書ける.

(\frac{∂^{2}G}{∂m_{i}∂m_{l}} ) = RT (\frac{δ_{il}}{m_{i}} − \frac{1}{∑ m_{i}} ) .  (1.6)

(1.5), (1.6) 式を (1.4) 式に代入し, RT で割ることによって,

\frac{Q(n_{i})}{RT} = \frac{G(m_{i})}{RT} + ∑_{i} \frac{μ_{i}}{RT} δm_{i} + \frac{1}{2} [∑_{i} \frac{(δm_{i})^{2}}{m_{i}} − \frac{\{∑(δm_{i})\}^{2}}{∑ m_{i}} ] ,  (1.7)

が導かれる.

(1.7) 式の Q(n_{i})/RT と質量保存則から, ラグランジュの未定係数法を用いて極値の条件を与える. L を以下のように定義すると,

L = \frac{G(m_{i})}{RT} + ∑_{i} \frac{μ_{i}}{RT} δm_{i} + \frac{1}{2} [∑_{i} \frac{(δm_{i})^{2}}{m_{i}} − \frac{\{∑(δm_{i})\}^{2}}{∑ m_{i}} ] − ∑_{e} χ_{e} [∑_{i} a_{i,e}n_{i} − B_{e}] ,

(1.7) 式が極値をとる条件は ∂L/∂n_{i}^{ϕ} = 0 なので,

\frac{μ_{i}}{RT} + \frac{n_{i}}{m_{i}} − \frac{∑ n_{i}}{∑ m_{i}} − ∑_{e} χ_{e}a_{ie} = 0,  (1.8)

となる. (1.8) 式と (1.1) 式を連立することで, 変数 n_{i}, n, χ_{e}(ラグランジュの未定係数) に関する N + M 個の線形方程式が得られ, 方程式系が閉じる.

(1.8) 式は式変形を行うことで連立方程式の次数を減らすことができる. まず (1.8) 式を n_{i} に関する式に書き直す.

n_{i} = m_{i} \{−\frac{μ_{i}}{RT} + \frac{∑ n_{i}}{∑ m_{i}} + ∑_{e} χ_{e}a_{ie}\} .  (1.9)

この (1.9) 式を i に関して 1 から N までの和をとると,

∑_{e} B_{e}χ_{e} = ∑_{i} m_{i} \frac{μ_{i}}{RT} ,  (1.10)

となる. また (1.9) 式の両辺に a_{if} をかけて i に関して和を取る.

∑_{e} \{∑_{i} a_{ie}a_{if}m_{i}\} χ_{e} + B_{f} (\frac{∑ n_{i}}{∑ m_{i}} − 1) = ∑_{i} a_{if}m_{i} \frac{μ_{i}}{RT} ,  (1.11)

である. a_{if} の添字 f (0 ≤ f ≤ N) はダミーサフィックスで, 実質的には e と同等である.

式変形の結果得られた (1.10), (1.11) 式は M + 1 個の方程式から成り, それらを連立することで χ_{e} と (∑ n_{i}/ ∑ m_{i} − 1) の値が求まる. 各物質のモル数は得られた値を (1.9) 式に代入することで求まる.

この式変形の利点は, N + M + 1 個の方程式を連立すべきところを M + 1 個の方程式の連立で済ませられる点にある. この変形によって逆変換する行列の大きさが M + 1 となるので, 計算量や計算時間を短縮できるものと期待される.

最終的に解くべき式をまとめる.

[figure]

温度・圧力一定, 初期状態での物質のモル数 m_{i}, 理想気体を仮定した場合, 平衡状態での物質のモル数は以下の式から得られる.

n_{i} = m_{i} \{−\frac{μ_{i}}{RT} + \frac{∑ n_{i}}{∑ m_{i}} + ∑_{e} χ_{e}a_{ie}\} .

但し, 上式に含まれる変数 χ_{e}, ∑ n_{i} は以下の 2 式を連立させることで得られる.

∑_{e} B_{e}χ_{e} = ∑_{i} m_{i} \frac{μ_{i}}{RT} ,

∑_{e} \{∑_{i} a_{ie}a_{if}m_{i}\} χ_{e} + B_{f} (\frac{∑ n_{i}}{∑ m_{i}} − 1) = ∑_{i} a_{if}m_{i} \frac{μ_{i}}{RT} .

[figure]

1.1.2 多相系への拡張

RAND 法を気相系から多相系へ拡張する (Van Zeggerren and Storey, 1970). 各相を独立に扱い, 多相系へ拡張する. 相の区別を添字 ϕ (1 ≤ ϕ ≤ Φ) で表し, 前章で導出した式を書き換える.

まず物質の保存則は以下のように書き換えることができる.

∑_{ϕ} ∑_{i} a_{ie}^{ϕ}n_{i}^{ϕ} = B_{e}  (1.12)

尚, (1.12) 式の B_{e} は各元素の全相における総和に拡張されている. 多相系で保存されているのは, 全相にわたる元素の存在量である.

次に理想系での化学ポテンシャルの表現を拡張する.

μ_{i}^{ϕ} = μ_{i}^{∘}^{ϕ}(T) + RT ln (\frac{m_{i}^{ϕ}}{∑_{i} m_{i}^{ϕ}} ) + αRT ln (\frac{p}{p_{0}} )  (1.13)

但し, 係数 α は気相と凝縮相を区別するために導入した定数で, 気相ならば 1, 液相・固相ならば 0 になる.

前章の方法を多相系に拡張した方程式に対して適用する. 適当な物質存在量 m_{i}^{ϕ} のまわりで 2 次のテーラー展開を行い, その極値を求める. (m_{i}^{ϕ} においても質量保存を満足していると仮定する.)

Q(n_{i}^{ϕ}) = G(m_{i}^{ϕ}) + ∑_{ϕ} ∑_{i} \underbrace{(\frac{∂G}{∂m_{i}^{ϕ}} )}_{A} δm_{i}^{ϕ} + \frac{1}{2} ∑_{ϕ} ∑_{i} ∑_{j} \underbrace{(\frac{∂^{2}G}{∂m_{i}^{ϕ}∂m_{l}^{ϕ}} )}_{B} δm_{i}^{ϕ}δm_{j}^{ϕ},  (1.14)

と書ける. ここで (1.14) 式の A 項は (1.13) 式を用いることによって,

(\frac{∂G}{∂m_{i}^{ϕ}} ) ≡ μ_{i}^{ϕ} = μ_{i}^{∘}^{ϕ}(T) + RT ln (\frac{m_{i}^{ϕ}}{∑_{i} m_{i}^{ϕ}} ) + αRT ln (\frac{p}{p_{0}} ).  (1.15)

と書ける. また (1.14) 式の B 項を (1.13) 式を用いて書き下すと,

(\frac{∂^{2}G}{∂m_{i}^{ϕ}∂m_{l}^{ϕ}} ) = RT (\frac{δ_{il}}{m_{i}^{ϕ}} − \frac{1}{∑_{i} m_{i}^{ϕ}} ) ,  (1.16)

となる. (1.15), (1.16) 式を (1.14) 式に代入し, RT で割ることによって,

\frac{Q(n_{i}^{ϕ})}{RT} = ∑_{ϕ} \frac{G(m_{i}^{ϕ})}{RT} + ∑_{ϕ} ∑_{i} \frac{μ_{i}^{ϕ}}{RT} δm_{i} + \frac{1}{2} ∑_{ϕ} [∑_{i} \frac{(δm_{i}^{ϕ})^{2}}{m_{i}^{ϕ}} − \frac{\{∑_{i}(δm_{i}^{ϕ})\}^{2}}{∑_{i} m_{i}^{ϕ}} ] , (1.17)

が導かれる. (1.17) 式の Q(n_{i})/RT と質量保存則から, ラグランジュの未定係数法を用いて極値の条件を与える. L を以下のように定義すると,

L = ∑_{ϕ} \frac{G(m_{i}^{ϕ})}{RT} + ∑_{ϕ} ∑_{i} \frac{μ_{i}^{ϕ}}{RT} δm_{i} + \frac{1}{2} ∑_{ϕ} [∑_{i} \frac{(δm_{i}^{ϕ})^{2}}{m_{i}^{ϕ}} − \frac{\{∑_{i}(δm^{ϕ})\}^{2}}{∑_{i} m^{ϕ}} ] − ∑_{e} χ_{e} ∑_{ϕ} [∑_{i} a_{ie}^{ϕ}n_{i}^{ϕ} − B_{e}^{ϕ}] ,

極値をとる条件は ∂L/∂n_{i}^{ϕ} = 0 なので,

\frac{μ_{i}^{ϕ}}{RT} + \frac{n_{i}^{ϕ}}{m_{i}^{ϕ}} − \frac{∑_{i} n_{i}^{ϕ}}{∑_{i} m_{i}^{ϕ}} − ∑_{e} χ_{e}a_{ie}^{ϕ} = 0,  (1.18)

である. (1.18) 式と (1.12) 式を連立することで, 変数 n_{i}^{ϕ}, χ_{e} (ラグランジュの未定係数) に関する ΦN + M 個の線形方程式が得られ, 方程式系が閉じる. 単相系と同様に (1.18) 式は式変形を行うことで連立方程式の次数を減らすことができる. まず (1.18) 式を n_{i}^{ϕ} に関する式に書き直す.

n_{i}^{ϕ} = m_{i}^{ϕ} \{−\frac{μ_{i}^{ϕ}}{RT} + \frac{∑_{i} n_{i}^{ϕ}}{∑_{i} m_{i}^{ϕ}} + ∑_{e} χ_{e}a_{ie}^{ϕ}\} .  (1.19)

この (1.19) 式を i に関して 1 から N までの和をとると,

∑_{i} ∑_{e} a_{ie}^{ϕ}m_{i}^{ϕ}χ_{e} = ∑_{i} m_{i}^{ϕ} \frac{μ_{i}^{ϕ}}{RT}  (1.20)

となる. また (1.19) 式の両辺に a_{if}^{ϕ} をかけて i と ϕ に関して和を取る.

∑_{e} \{∑_{i} ∑_{ϕ} a_{ie}^{ϕ}a_{if}^{ϕ}m_{i}^{ϕ}\} χ_{e} + ∑_{ϕ} ∑_{i} a_{if}^{ϕ}m_{i}^{ϕ} (\frac{∑_{i} n_{i}^{ϕ}}{∑_{i} m_{i}^{ϕ}} ) = ∑_{i} ∑_{ϕ} a_{if}^{ϕ}m_{i}^{ϕ} \frac{μ_{i}^{ϕ}}{RT} + B_{f} (1.21)

である. a_{if} の添字 f (0 ≤ f ≤ M) はダミーサフィックスで, 実質的には e と同等である.

(1.20), (1.21) 式は M + Φ の方程式であり, また変数の数は χ_{e} (M 個) と ∑_{i} n_{i}^{ϕ} (Φ 個) なので, この 2 式を連立することによって方程式が閉じる.

最後に解くべき式は以下のようになる.

[figure]

温度・圧力一定, 初期状態での物質のモル数 m_{i}^{ϕ}, 理想気体を仮定した場合, 平衡状態での物質のモル数は以下の式から得られる.

n_{i}^{ϕ} = m_{i}^{ϕ} \{−\frac{μ_{i}^{ϕ}}{RT} + \frac{∑_{i} n_{i}^{ϕ}}{∑_{i} m_{i}^{ϕ}} + ∑_{e} χ_{e}a_{ie}^{ϕ}\} .

但し, 上式に含まれる変数 χ_{e}, n^{ϕ} は以下の 2 式を連立させることで得られる.

∑_{i} ∑_{e} a_{ie}^{ϕ}m_{i}^{ϕ}χ_{e} = ∑_{i} m_{i}^{ϕ} \frac{μ_{i}^{ϕ}}{RT}

∑_{e} \{∑_{i} ∑_{ϕ} a_{ie}^{ϕ}a_{if}^{ϕ}m_{i}^{ϕ}\} χ_{e} + ∑_{ϕ} ∑_{i} a_{if}^{ϕ}m_{i}^{ϕ} (\frac{∑_{i} n_{i}^{ϕ}}{∑_{i} m_{i}^{ϕ}} ) = ∑_{i} ∑_{ϕ} a_{if}^{ϕ}m_{i}^{ϕ} \frac{μ_{i}^{ϕ}}{RT} + B_{f}

[figure]

1.1.3 負のモル数になった場合には

(1.19), (1.20), (1.21) 式を連立することによって, 平衡状態に近い物質存在量 n_{i}^{ϕ} が求まる. しかし得られた n_{i}^{ϕ} は負になる可能性がある. もしも計算から得られた存在量 n_{i}^{ϕ} が負の場合, 以下の式より得られる n_{i}^{ϕ∗} をより平衡状態に近い物質存在量とみなす.

n_{i}^{ϕ∗} = m_{i}^{ϕ} + λ{n_{i}^{ϕ} − m_{i}^{ϕ}}.  (1.22)

ここで λ(> 0) は定数で, 以下のように定義する.

λ ≡ \frac{m_{i}^{ϕ}}{m_{i}^{ϕ} − n_{i}^{ϕ}} × 0.999 (n_{i}^{ϕ} ≤ 0),  (1.23)
λ ≡ 1 (n_{i}^{ϕ} > 0).  (1.24)

なお, (1.22) 式の n_{i}^{ϕ∗} は元素保存の式 ((1.1) 式) を満たす.

1.2 RAND 法の適用できないケースとその解決方法

前節までで議論した RAND 法を適用できないケースがいくつか存在する. そのようなケースは主に行列の非正則性として現れる. 代表的なケースとして,

• ある 2 つの相に存在する化学種が全く同じな場合. 例えば, H_{2}O(l) と H_{2}O(s) がそれぞれ独立に 1 つの相を形成する場合がこれに相当する. この際, 係数行列の 2 つの行の要素が全く同じとなり, 行列式が非正則となる.

• ある相に存在する全ての物質の存在量が 0 に近付いた場合. この際, 行列式が 0x = 0 となり解が不定となる.

が挙げられる. 以下の節では非正則になる場合を示し, さらにその解消方法を述べる.

1.2.1 ある 2 つの相に存在する物質が全く同じな場合

RAND 法では, 熱力学関数 G の関数形は物質量に関する 2 階以上の関数であることを前提としている. しかし, 1 成分 2 相系の熱力学関数 G の関数形は物質量に関する 1 次方程式である. そのため RAND 法では扱うことができない. このケースを避けるためには, RAND 法による計算を始める前 (温度圧力を与えた段階) において存在し得る化学種を予め判定する必要がある. 各化学種の化学ポテンシャルを比較し, 化学ポテンシャルの大きい化学種を削除すれば良い.

具体例としてまず 1 成分 2 相系の場合を考える. 例として H_{2}O(l), H_{2}O(s) からなる 1 成分 2 相系を考える. 温度・圧力を与えると, 熱力学関数 G は以下のように書ける.

G = ∑_{ϕ} μ^{ϕ}n^{ϕ}
= μ^{l}n^{l} + μ^{s}n^{s}  (1.25)

ここで添字 ϕ は相を表し, l は液相, s は固相である. μ は化学ポテンシャル, n はモル数である. さらに物質の保存則は以下のように書ける.

n^{l} + n^{s} = n  (1.26)

物質は保存されるので n は定数である. (1.25) 式に (1.26) 式を代入し, さらに化学ポテンシャルを標準化学ポテンシャル μ^{∘} を用いて書き直すと,

G = μ^{∘}^{l}(T)n^{l} + μ^{∘}^{s}(T)(n − n^{l})
= {μ^{∘}^{l}(T) − μ^{∘}^{s}(T)}n^{l} + μ^{∘}^{s}(T)n  (1.27)

(1.27) 式より, 1 成分 2 相系の G は物質量に関する 1 次関数であることがわかる. (1.27) 式の最小値を与える n^{l} は以下のように決まる:

μ^{∘}^{l}(T) < μ^{∘}^{s}(T) の場合
n^{l} = n
μ^{∘}^{l}(T) > μ^{∘}^{s}(T) の場合
n^{l} = 0
μ^{∘}^{l}(T) = μ^{∘}^{s}(T) の場合 (相平衡の場合)
0 ≤ n^{l} ≤ n

従って標準化学ポテンシャルの値から, どの相に存在するか判定することができる.

次に多成分多相系においても, ある 2 つの相に含まれる物質が全く同じ場合を考える. 結論から言えば, この場合にも 1 成分 2 相系と同様の問題が生じる. 例として多成分系ではあるが, H_{2}O(l), H_{2}O(s) がそれぞれ単独で 1 つの相を形成する場合を考える. (1.20) 式を H_{2}O(l), H_{2}O(s) に関する式に変形すると,

∑_{e} a_{H_{2}O,e}^{l}m_{H_{2}O}^{l}χ_{e} = m_{H_{2}O}^{l}μ_{H_{2}O}^{l}
∑_{e} a_{H_{2}O,e}^{s}m_{H_{2}O}^{s}χ_{e} = m_{H_{2}O}^{s}μ_{H_{2}O}^{s}

両辺を物質量で通分すると,

∑_{e} a_{H_{2}O,e}χ_{e} = μ_{H_{2}O}^{l}
∑_{e} a_{H_{2}O,e}χ_{e} = μ_{H_{2}O}^{s}  (1.28)

ここで, a_{H_{2}O,e}^{l} = a_{H_{2}O,e}^{s} の関係を用いた. 物質が同じなので, 相が異なっても a_{ie} の値は変化しない. (1.28) 式を見ると, 左辺の係数は同じであるが, 右辺の値は異なることがわかる. これでは χ_{e} を求める行列式が非正則となり, 解が求まらない.

1.2.2 ある相に含まれる物質の存在量が 0 に近付く場合

ある相に存在する全ての物質存在量が 0 に近付くと, (1.20) 式は

0 · χ_{e} = 0  (1.29)

となり, χ_{e} は不定となる. このケースを避けるためには, 当該温度圧力条件において存在できない化学種を「推定」し, 系から取り除くという作業が必要となる. 経験上, 単純に相に存在する全ての物質の存在量が 0 に近付いた場合にその相が存在しないものとして計算を進めると, 存在するはずの化学種まで系から取り除いてしまう可能性がある.

平衡条件から存在しない化学種を推定する. 化学種の気相のモル比はだいたい保存されると仮定した時に, 気相と凝縮相の平衡条件を満たすために必要な凝縮相のモル比を計算し, その値と実際の凝縮相のモル比を比較する. そして凝縮相のモル比が実際のモル比よりも著しく小さくなっている化学種を存在しないものとする.

行列が非正則になったときのある化学種の気相のモル比を x^{g}, 凝縮相のモル比を x^{l} とする. モル比を x^{g} の気相と平衡する凝縮相のモル比 x^{l}′ は以下のように求まる.

μ_{0}^{g} + RT ln x^{g} + RT ln p/p_{ref} = μ_{0}^{l} + RT ln x^{l}′  (1.30)
\frac{(μ_{0}^{g} + RT ln p/p_{ref}) − μ_{0}^{l}}{RT} = ln x^{l}′/x^{g}  (1.31)
x^{l}′ = exp (\frac{(μ_{0}^{g} + RT ln p/p_{ref}) − μ_{0}^{l}}{RT} )x^{g}  (1.32)

但し R は気体定数, T は温度, p は圧力, p_{ref} は標準圧力, μ_{0} は標準化学ポテンシャルである. 以上より求めた x^{l}′ と実際の凝縮相のモル比 x^{l} との比を取り, 最も値が小さいもの, すなわち平衡するためにはモル比が著しく減少せねばならないもの, を存在しえない化学種であると「推定」する.

1.3 反復計算手法

(1.19), (1.20), (1.21) 式を連立させることにより得られる n_{i}^{ϕ} を平衡点まで収束させる必要がある (図 1.1 参照). 平衡点に収束させるために反復法を用いる.

第 1.2 節で述べたように, RAND 法は一成分一相系を扱うことができず, また行列式が正則かどうか常にチェックする必要がある. そのため, 以下の手順で反復計算を繰り返す.

[figure]

(1) 温度, 圧力, 適当な組成 m_{i}^{ϕ} を与える.

(2) 気相と凝縮相の平衡条件, 及び液相と固相の平衡条件を用いて, 当該温度, 圧力, 組成に対して存在しえない凝縮性化学種を削除する. 気相と凝縮相の平衡条件は必ずしも必要ないが, 計算時間の短縮に大きく貢献する. 液相と固相の平衡条件は, 第 1.2.1 節で述べたように, RAND 法の破綻を招かないために必要. また各元素量をチェックし, 微小量な元素に関しては, その元素から成る化学種を存在しないものとみなし削除する.

(3) (1.20), (1.21) 式を連立し行列式を作成する. 必要となる各化学種の化学ポテンシャルは, 比熱や飽和蒸気圧より計算する.

(4) 行列式が正則かどうか判定する. 非正則な場合には, 液相と固相の平衡条件を用いた存在可能物質のチェックを行い, 反復の前後での m_{i}^{ϕ} の変化率が最も大きい相を削除し, 手順 (3) に戻る. この段階で再び液相と固相の平衡条件を用いた存在可能物質のチェックを行うのは, 手順 (6) で微小量となった化学種を削除するためである. 削除した結果, 溶液を維持できず, 単体の液相, 固相が計算中に現れるケースがままあるからである.

(5) 行列式を計算し χ_{e}, ∑_{i} n_{i}^{ϕ} を求める. それらを (1.19) 式に代入し, より平衡に近い物質存在量 n_{i}^{ϕ} とその時のギブス自由エネルギー G の値を求める.

(6) 極小量となった化学種 (気相) を存在しないものとみなし削除する. ある化学種が一方的に無限に 0 に近付くことで, ローカルミニマムにおちこむことがある. この場合, 行列は正則のままである.

(7) (6) で得られた G が収束したか判定する. 収束の判定条件としては, 反復の前後でのギブス自由エネルギーの相対誤差がマシンイプシロン程度になったか否かで判断する. 収束していない場合は手順 (5) で得られた n_{i}^{ϕ} を m_{i}^{ϕ} として手順 (4) に戻る. 収束した場合には n_{i}^{ϕ} を平衡組成とみなし反復計算を止める.

[figure]

図 1.1: 逐次近似法の概要. 適当な物質存在度のまわりでギブス自由エネルギー G を 2 次関数に近似し極値を求める. 得られた極値を使って同様の操作を繰り返し, 平衡時の物質存在量を求める. 太線は大気のギブス自由エネルギー, 1, 2, 3, ... と番号を付けた細線は 2 次の微少量まで考慮して近似した時のギブス自由エネルギー.

\underline{朧モデルの差分化方法第 1 章大気の平衡状態の記述 13} {図 1.2: 大気の平衡計算の流れ図. 計算に必要となる化学ポテンシャルは比熱や飽和蒸気圧等の既 \atop 存のデータから予め計算しておく.}

\overline{main.tex 2003/06/22(杉山耕一朗, 小高正嗣, 倉本圭, 林祥介)}

\underline{朧モデルの差分化方法第 2 章大気の断熱変化の記述 14}

第2章大気の断熱変化の記述

擬湿潤断熱線に沿った温度, 圧力, 組成を以下の手順で求める. 通常, 基準温度圧力として大気上端の温度圧力が与えられるので, 以下の手順で得られた断熱線のうち基準温度圧力を通るものを探すこととなる.

[figure]

1. 初期圧力, 初期温度での平衡組成を求め, エントロピーを計算する. 平衡組成は温度・圧力固定, 元素数保存の条件のもとギブスの自由エネルギー G を最小化させることで求まる (第 1 章参照).

2. 圧力を dp だけ変化させる. 前のステップで凝縮が生じた場合, 空気塊から凝縮物を除去し, 凝縮物以外のもつエントロピーが保存するものとみなす.

3. 変化させた圧力面において前のステップでのエントロピーを保存するような平衡組成, 温度を反復的に求める.

4. 対流圏上層まで手順 (2) – (3) のステップを繰り返すことで擬湿潤断熱線が求まる.

[figure]

上記の手順 2 (図 2.1 の Step 2) において, 温度を変化させながら前のステップでのエントロピーにシューティングさせる必要がある. シューティングの手順として以下の方法を採用する (前のステップでのエントロピーを S_{pre} とする).

[figure]

ステップ 1 適当な温度幅 δT の間隔でエントロピーを計算し, S_{pre} をはさむ温度を探す.

ステップ 2 ステップ 1 で得られた 2 点の温度の中間値でのエントロピー S_{now} を計算し, S_{pre} と比較する. S_{now} と S_{pre} の変化率がマシンイプシロン程度になった場合には, S_{now} と S_{pre} は一致したものとみなす. 一致しない場合は, さらに中間値でのエントロピーを計算し S_{pre} と比較する.

[figure]

図 2.1: 断熱曲線 dS = 0 の計算方法. Step 1: 初期温度での平衡組成を計算し, エントロピー S_{0} を求める (左上). Step 2: 圧力を dp だけ変化させる. 温度を変化させた時のエントロピーを順次計算し, 前のステップでのエントロピー S_{0} と一致する温度を求める (右上). Step 3: (T_{1}, p_{1}) において凝縮が生じた場合, 次のステップにおいて保存させるエントロピーは (T_{1}, p_{1}) での気相のエントロピー S_{4}_{gas} (左下). Step 4: Step 1 から Step 3 において得られた温度・圧力を結んで断熱線を引く (右下).

図 2.2: エントロピーを保存するための方法. ステップ 1 では, 前のステップでのエントロピーをはさむ温度を求める. ステップ 2 では温度の中間値でのエントロピーを計算しつづけ, 前のステップでのエントロピーと一致する温度を求める.

\underline{朧モデルの差分化方法第 3 章化学ポテンシャル, エントロピーの計算方法 17}

第3章{化学ポテンシャル, エントロピーの計 \atop 算方法}

oboro モデルにおいて必要とされる物性値は化学ポテンシャル及びエントロピーである. 大気中の物質が全て理想気体/ 液体/ 固体/ 溶液であるとすると, 化学ポテンシャル μ_{i}^{ϕ}(T, p, n_{i}^{ϕ}) は,

μ_{i}^{ϕ}(T, p, n_{i}^{ϕ}) = μ_{i}^{∘}^{ϕ}(T) + RT ln (\frac{m_{i}^{ϕ}}{∑_{i} m_{i}^{ϕ}} ) + αRT ln (\frac{p}{p_{0}} )
= μ_{i}^{∘}^{ϕ}(T) + RT ln (\frac{p_{i}^{ϕ}}{p_{0}} )  (3.1)
μ_{i}^{∘}(T) = h_{i}^{ϕ}(T_{0}) − T s_{i}^{ϕ}(T_{0}, p_{0}) + ∫_{T_{0}}^{T} c_{p}_{i}^{ϕ}dT − T ∫_{T_{0}}^{T} \frac{c_{p}_{i}^{ϕ}}{T} dT  (3.2)

と書け, モルエントロピー s_{i}^{ϕ}(T, p, n_{i}^{ϕ}) は,

s_{i}^{ϕ}(T, p, n_{i}^{ϕ}) = s_{i}^{∘}^{ϕ}(T) − \frac{∂}{∂T} \{RT ln (\frac{m_{i}^{ϕ}}{∑_{i} m_{i}^{ϕ}} ) + αRT ln (\frac{p}{p_{0}} )\}
= s_{i}^{∘}^{ϕ}(T) − \frac{∂}{∂T} \{RT ln (\frac{p_{i}^{ϕ}}{p_{0}} )\} + R ln (p_{0})  (3.3)
s_{i}^{∘}^{ϕ}(T) = s_{i}^{ϕ}(T_{0}, p_{0}) + ∫_{T_{0}}^{T} \frac{c_{p}_{i}^{ϕ}}{T} dT  (3.4)

と書ける. 但し, μ_{i}^{∘}^{ϕ}(T) は標準化学ポテンシャル, s_{i}^{∘}^{ϕ}(T) は標準モルエントロピー, c_{p}_{i}^{ϕ} は比熱, T_{0} = 298.15 K は基準温度, p_{0} = 1.0 × 10^{5} Pa は基準圧力, p_{i}^{ϕ} = α ∗ p ∗ m_{i}^{ϕ}/ ∑_{i} m_{i}^{ϕ} は分圧である. (3.1) — (3.4) 式より, 化学ポテンシャル及びエントロピーは温度, 圧力, 組成, 基準状態のエントロピーとエンタルピーに加え, 各化学種の飽和蒸気圧もしくは比熱から求まることがわかる.

エントロピーと化学ポテンシャルの関係式は s = −∂μ/∂T なので, 化学ポテンシャルを数値微分することでエントロピーを求めることもできる. しかし数値微分による数値誤差は大きいので, エントロピーも飽和蒸気圧もしくは比熱から求める方が無難である.

3.1 比熱を用いた標準化学ポテンシャル, 標準モルエントロピーの計算

比熱は NIST Chemical Web Book, 化学工学便覧, 化学便覧より簡単に得られる. しかし温度に対して離散的にしか与えられていないので, まず比熱を温度に対する方程式として表現する必要がある. oboro モデルでは比熱を 3 次のスプライン補完している. 以下ではスプライン補間した場合の標準化学ポテンシャル, 標準モルエントロピーの計算方法を述べる.

3.1.1 希ガス

理想気体を仮定した場合, 希ガスの比熱は 20.786 J K^{−1} mol^{−1} の一定値である. そのため簡単に計算することができる.

μ_{i}^{∘}(T) = h_{i}^{ϕ}(T_{0}) − T s_{i}^{ϕ}(T_{0}, p_{0}) + 20.786 \{(T − T_{0}) − T ln \frac{T}{T_{0}} \}  (3.5)
s_{i}^{∘}^{ϕ}(T) = s_{i}^{ϕ}(T_{0}, p_{0}) + 20.786 ln \frac{T}{T_{0}}  (3.6)

3.1.2 スプライン補間

oboro では科学計算サブルーチン SSL II の INSPL ルーチンを用いて補間式を作る. INSPL ルーチンの補間係数 C_{j}, D_{j}, E_{j}(1 ≤ j ≤ n) を用いることにより, 比熱 c_{p}_{j}^{ϕ}(T) は以下のように書ける.

c_{p}_{j}^{ϕ}(T) = c_{p}_{i}^{ϕ}_{1} + C_{1}(T − T_{1}) + D_{1}(T − T_{1})^{2} T < T_{1}

c_{p}_{i}^{ϕ}(T) = c_{p}_{i}^{ϕ}_{j} + C_{j}(T − T_{i}) + D_{j}(T − T_{i})^{2} + E_{j}(T − T_{i})^{3} T_{j} ≤ T ≤ T_{j+1} (1 ≤ j ≤ n − 1)

c_{p}_{i}^{ϕ}(T) = c_{p}_{i}^{ϕ}_{n} + C_{n}(T − T_{n}) + D_{n}(T − T_{n})^{2} T > T_{n}

但し c_{p}_{i}^{ϕ}_{j} は物性表から得た温度 T_{j} での比熱とする.

∫_{T_{0}}^{T} c_{p}_{i}^{ϕ}dt = ∫_{T_{0}}^{T_{j}} c_{p}_{i}^{ϕ}dT + ∫_{T_{j}}^{T} c_{p}_{i}^{ϕ}dT
= ∑_{j=a}^{j−1} (c_{p})_{int}^{j}

+C_{j} \{\frac{1}{2}(T^{2} − T_{j}^{2}) + T_{j}(T − T_{j})\}

+D_{j} \{\frac{1}{3}(T^{3} − T_{j}^{3}) − T_{j}(T^{2} − T_{j}^{2}) + T_{j}^{2}(T − T_{j})\}

+E_{j} \{\frac{1}{4}(T^{4} − T_{j}^{4}) − T_{j}(T^{3} − T_{j}^{3}) + \frac{3}{2}T_{j}^{2}(T^{2} − T_{j}^{2}) − T_{j}^{3}(T − T_{j})\} (3.7)

∫_{T_{0}}^{T} \frac{c_{p}_{i}^{ϕ}}{T} dt = ∫_{T_{0}}^{T_{j}} \frac{c_{p}_{i}^{ϕ}}{T} dT + ∫_{T_{j}}^{T} \frac{c_{p}_{i}^{ϕ}}{T} dT
= ∑_{j=a}^{j−1} (\frac{c_{p}}{T} )_{int}^{j} + c_{p}_{i}^{ϕ} ln (T/T_{j})

+C_{j} {(T − T_{j}) + T_{j} ln(T/T_{j})}

+D_{j} \{\frac{1}{2}(T^{2} − T_{j}^{2}) − 2T_{j}(T − T_{j}) + T_{j}^{2} ln(T/T_{j})\}

+E_{j} \{\frac{1}{3}(T^{3} − T_{j}^{3}) − \frac{3}{2}T_{j}(T^{2} − T_{j}^{2}) + 3T_{j}^{2}(T − T_{j}) − T_{j}^{3} ln(T/T_{j})\} (3.8)

但し T_{j} = T_{0} となる j を a = j とした. (c_{p})_{int}^{j}, (c_{p}/T)_{int}^{j} はそれぞれ, c_{p}_{i}^{ϕ}, (c_{p}_{i}^{ϕ}/T) を T_{j−1} から T_{j} まで積分した値である. 従って, (3.7), (3.8) 式を (3.2), (3.4) 式の比熱の積分項に以下の値を代入することで標準化学ポテンシャル, 標準モルエントロピーが求まる.

3.2 飽和蒸気圧を用いた標準化学ポテンシャル, 標準モルエントロピーの計算

凝縮相の場合, 比熱よりも飽和蒸気圧の方がデータとして入手しやすい. そこで飽和蒸気圧を用いた化学ポテンシャル, エントロピーの計算方法をまとめる. 相平衡が成立しているならば凝縮相と気相の化学ポテンシャルは等しいので,

μ_{i}^{c}^{∘} + RT ln x_{i}^{c} = μ_{i}^{g}^{∘} + RT ln{x_{i}^{c}p^{∗}/p_{0}}
μ_{i}^{c}^{∘} = μ_{i}^{g}^{∘} + RT ln{p^{∗}/p_{0}}  (3.9)

となる. ここで μ_{i}^{c}^{∘} は凝縮相の標準化学ポテンシャル, μ_{i}^{g}^{∘} は気相の標準化学ポテンシャル, x_{i} はモル比, p^{∗} は飽和蒸気圧である. 以上から気相の化学ポテンシャルと飽和蒸気圧より, 凝縮相の化学ポテンシャルが計算できる.

NIST Chemical Web Book, 化学工学便覧, 化学便覧より, 飽和蒸気圧の式が得られる. 以下では代表的な 2 つの経験式を用いた場合での化学ポテンシャル, エントロピーの定式化を行う.

3.2.1 Antoine 式

飽和蒸気圧の式として有名なのは Antoine の式である. Antoine の式は,

ln\{p^{∗}(T)\} = A − \frac{B}{C + T}  (3.10)

の形式で与えられる. T は温度, A, B, C は定数であり, その値は表 3.1 に示されている. この式はクラウジウス–クラペイロンの式と同じ温度依存性を持つ.

標準化学ポテンシャルは (3.2) 式, (3.10) 式より,

μ_{i}^{c}^{∘} = μ_{i}^{g}^{∘}(T) + RT (A − \frac{B ln (10)}{C + (T − 273.15)}) − RT ln p_{0} + RT ln (133.322).  (3.11)

標準モルエントロピーは (3.4) 式, (3.10) 式より,

s_{i}^{c}^{∘} = s_{i}^{g}^{∘}(T) − R (A − \frac{B ln (10)}{C + (T − 273.15)}) − RT (\frac{B ln (10)}{(C + (T − 273.15))^{2}} ) + R ln p_{0} (3.12)

表 3.1 の単位系を SI に変換する必要があるために, (3.11), (3.12) 式に補正項を加えてある.

3.2.2 AMP 式

飽和蒸気圧の式として AMP 式がある. この式は分子構造に基づいた蒸気圧の推算式 (グループ寄与法) である. 詳細はアペンデックス参照.

ln p^{∗}(T) = \frac{a_{1}}{T} + a_{2} + a_{3} ln T + a_{4}T + a_{5}T^{2}  (3.13)

表 3.1: Antoine の式の定数 A, B, C の値 (化学工学便覧より). 但し T の単位は^{∘}C, p^{∗} の単位は mmHg とする.

物質名  A  B  C
H_{2}O(l) 7.9186968 1636.909 224.92
H_{2}S(l)  6.99392 768.1302 247.093
H_{2}S(s)  7.61418 885.319 250.25
NH_{3}(l)  7.3605 926.132 240.17
NH_{3}(s)  9.96382 1617.907 272.55
CH_{4}(l)  6.61184  389.93  266.0
CH_{4}(s)  7.6954  532.2  275.0
CO_{2}(s)  9.81062 1347.79 272.99

で与えられる. a_{1}, a_{2}, a_{3}, a_{4}, a_{5} は定数であり, 代表的な定数は表 3.2 に与えられる.

標準化学ポテンシャルは (3.2) 式, (3.10) 式より,

μ_{l}(T, p^{∗}(T)) = μ_{g}^{∘}(T) + RT (\frac{a_{1}}{T} + a_{2} + a_{3} ln T + a_{4}T + a_{5}T^{2}) − RT ln p_{0} + RT ln (0.1) (3.14)

但し, 表 3.2 の単位系は cgs なので, SI 単位系に変換するために係数をかけてある. 標準モルエントロピーは (3.4) 式, (3.10) 式より,

s_{i}^{c}^{∘} = s_{i}^{g}^{∘}(T) − R (\frac{a_{1}}{T} + a_{2} + a_{3} ln T + a_{4}T + a_{5}T^{2})
−RT (−\frac{a_{1}}{T^{2}} + \frac{a_{3}}{T} + a_{4} + 2a_{5}T ) + R ln p_{0}  (3.15)

表 3.2: AMP 式の定数 a_{1}, a_{2}, a_{3}, a_{4}, a_{5} の値 (Table 4 of Briggs and Sacket, 1989, ICARUS, 80, p98.). 但し単位系は cgs.

物質名  a_{1}  a_{2}  a_{3}  a_{4}  a_{5}
NH_{3}(s)  -4122  41.67871  -1.81630  0  0
NH_{3}(l) -4409.3512 76.864252 -8.4598340  5.51029 × 10^{−3} 6.804632 × 10^{−6}
H_{2}O(s) -5631.1206 -8.363602 8.2312 -3.861449 × 10^{−2} 2.77494 × 10^{−5}
H_{2}O(l) -2313.0338 -164.03307 38.053682 -1.3844344 × 10^{−1} 7.4465367 × 10^{−5}

3.3 化学反応の圧平衡定数を用いた標準化学ポテンシャル, 標準モルエントロピーの計算

化学種の比熱, 飽和蒸気圧の知られていなくとも, 化学反応の圧平衡定数を用いて標準化学ポテンシャル, 標準モルエントロピーを計算することができる. Lewis (1979) によって, 木星大気中に NH_{4}SH の雲が生じる可能性が示唆されている. この反応の圧平衡定数を用いて NH_{4}SH の標準化学ポテンシャル, 標準モルエントロピーを求める.

NH_{4}SH の生成反応は以下のように書ける.

NH_{3} + H_{2}S → NH_{4}SH  (3.16)

Briggs and Sackett (1989) によれば, (3.16) 式の圧平衡定数は以下のように書くことができる.

K_{p} = ln (p_{NH3} · p_{H_{2}S}) = 61.781 − \frac{10834}{T}  (3.17)

但し, p_{NH3}, p_{H_{2}S} はそれぞれ, アンモニア気体の分圧, 硫化水素の分圧である. 単位系は cgs 単位系である (圧力の単位は dynes/cm^{2}).

平衡条件より, NH_{4}SH の標準化学ポテンシャルは NH_{3} 及び H_{2}S の標準化学ポテンシャルを用いて以下のように書ける.

μ_{NH_{4}SH}^{s} = μ_{NH_{3}}^{g} + μ_{H_{2}S}^{g}
= \{μ_{NH_{3}}^{g}^{∘} + RT ln (\frac{p_{NH3}^{g}}{p_{0}} )\} + \{μ_{H_{2}S}^{g}^{∘} + RT ln (\frac{p_{H_{2}S}^{g}}{p_{0}} )\}
= μ_{NH_{3}}^{g}^{∘} + μ_{H_{2}S}^{g}^{∘} + RT ln (\frac{p_{NH3}^{g} · p_{H_{2}S}^{g}}{p_{0}^{2}} )
= μ_{NH_{3}}^{g}^{∘} + μ_{H_{2}S}^{g}^{∘} + RT ln (p_{NH3} · p_{H_{2}S}) − RT ln p_{0}^{2}
= μ_{NH_{3}}^{g}^{∘} + μ_{H_{2}S}^{g}^{∘} + RT (61.781 − \frac{10834}{T})
−RT ln (1.0 × 10^{2}) − RT ln p_{0}^{2}.  (3.18)

標準モルエントロピーは以下のようになる.

s_{NH_{4}SH}^{s} = −\frac{∂μ_{NH_{4}SH}^{s}}{∂T}
= s_{NH_{3}}^{g}^{∘} + s_{H_{2}S}^{g}^{∘} − R (61.781 − \frac{10834}{T})
−RT (\frac{10834}{T^{2}}) + R ln (1.0 × 10^{2}) + R ln p_{0}^{2}.  (3.19)

但し, μ_{NH_{3}}^{g}^{∘}(T) はアンモニア気体の標準化学ポテンシャル, μ_{H_{2}S}^{g}^{∘}(T) は硫化水素気体の標準化学ポテンシャルである. また単位を SI 単位系にするために補正項が加えられている.

\underline{朧モデルの差分化方法第 4 章雲密度の計算方法 22}

第4章雲密度の計算方法

静水圧平衡にある大気での雲密度は Weidenschilling and Lewis (1973) の方法によって導くことができる. 静水圧平衡にある大気でのある高度 J において

M^{J} = \frac{p^{J}}{g}  (4.1)

が成り立つ. ここで M^{J} は高度 J より上における大気の単位面積当たりの質量である. 物質 k の高度 J より上での単位面積当たりの大気の質量は

m_{k}^{J} · M^{J},  (4.2)

である. ここで m_{k}^{J} は大気の質量分率である. それゆえ物質 k が凝縮する際, 高度 I から J の変化幅 dz における平均雲密度 D̄ は,

D̄ = \frac{(m_{k}^{I} − m_{k}^{I})M̄}{dz} ,  (4.3)

である. ここで M̄ は考えている領域での単位面積当たりの平均質量である. (4.4) 式はモル分率と圧力の項によって書き換えることができる.

D = \frac{m_{k}(X_{k}^{I} − X_{k}^{J})p̄}{m̄gdz}.  (4.4)

但し m̄ は大気の平均分子量, m_{k} は物質 k の分子量である. 分子量の単位は kg mol^{−1} に換算する必要がある.

さらに高度の変化 dz を圧力変化 dp で書き換える. 静水圧平衡の式と理想気体の状態方程式を用いると,

dp = −ρgdz,
= −\frac{m̄p}{RT} gdz,
dz = −\frac{RT g}{m̄p} dp,  (4.5)

である. (4.4) 式に (4.5) 式を代入することによって,

D = \frac{m_{k}(X_{k}^{I} − X_{k}^{J})p̄^{2}}{RT dp} .  (4.6)

が求まる. もし 1 つ以上の物質が凝縮して雲を作る場合 (NH_{4}SH やアンモニア水溶液), (4.6) 式の右辺は凝縮する物質を足し合わせたものになる.

\underline{朧モデルの差分化方法第 5 章静的安定度の計算 23}

第5章静的安定度の計算^{本章では静的安定度の計算方法を述べる.}

\overline{main.tex 2003/06/22(杉山耕一朗, 小高正嗣, 倉本圭, 林祥介)}

\underline{朧モデルの差分化方法第 6 章参考文献 24}

第6章参考文献

American Chemical Soc, American Institute of Physics, National Bureau of Standards, 1986:

JANAF Thermochemical Tables (3rd. Edition), Part I, II. American Chemical Society, Atreya, S.K., Romani. P.N., 1985: Photochemistry and clouds of Jupiter, Saturn and Uranus.

In: Hunt. G.E. (Ed.), Planetary Meteorology. Cambridge University Press, pp. 17–68. Briggs, F. H., Sackett, P. D., 1998: Radio Observations of Saturn as a Probe of its Atmosphere

and Cloud Structure. ICARUS, 80, p77–103. Lewis, J.S., 1969: Observability of Spectroscopically Active Compounds in the Atmosphere of

Jupiter. ICARUS, 10, 393–409 National Institute of Standards and Technology, 2000: NIST Chemical Web Book. http://webbook.nist.gov/chemistry/. Seiff, A., Kirk, D.B., Knight, T.C.D., Young, R.E., Mihalov, J.D., Young, L.A., Milos. F.S.,

Schubert, G., Blanchard, R.C., Atkinson, D., 1998: Thermal structure of Jupiter’s atmosphere near the edge of a 5-μm hot spot in the north equatorial belt. J. Geophys. Res. 103, 22,857–22,889.

Van Zeggerren, F., Storey, S.H., 1970: The computation of chemical equilibria, Cambridge

Univ. Press. Weidenschilling, S.J., Lewis, J.S., 1973: Atmospheric and cloud structure of the Jovian planets.

ICARUS, 20, 465–476. White, W.B., Johnson, S.M., Dantzig, G.B., 1958: Chemical equilibrium in complex mixture.

J. Chem. Phys. 28, 751–755 化学工学会編, 1999: 化学工学便覧, 改訂 6 版. 丸善. 日本化学会編, 1999: 化学便覧, 改訂 4 版. 丸善. 久保亮五, 1998: 熱学・統計力学 (修正版). 裳華房. プリゴジーヌ, デフェイ, 1966: 化学熱力学 I, II. みすず書房.

\underline{朧モデルの差分化方法付録 A AMP 式 25}

付録A AMP 式

分子構造に基づいた蒸気圧の推算式として AMP 式がある (化学工学便覧改訂 6 版).

ln p^{∗} = A + \frac{B}{T} + C ln T + DT + ET^{2}  (A.1)

但し p^{∗} は飽和蒸気圧, T は温度, A–E までの定数は以下で与えられる.

A = ln \frac{R}{V_{w}} + (s − \frac{1}{2}) ln \frac{E_{0}}{R} − ln \{(s − 1)!\} + ln α.  (A.2)
B = −\frac{E_{0}}{R}  (A.3)
C = \frac{3}{2} − s  (A.4)
D = \frac{s − 1}{(E_{0}/R)}  (A.5)
E = \frac{(s − 3)(s − 1)}{2(E_{0}/R)}  (A.6)

但し,

α = 0.0966  (A.7)
(s − 1)! = exp[−s]s^{s−1/2}(2π)^{1/2} (1 + \frac{1}{12s} + \frac{1}{288s^{2}} − \frac{139}{51480s^{3}} − \frac{571}{2488320s^{4}} ) . (A.8)
s = ∑_{i} ν_{i}s_{i}.  (A.9)
E_{0}/R = ∑_{i} ν_{i}ε_{0}_{i}.  (A.10)
V_{w} = ∑_{i} ν_{i}V_{w}_{i}.  (A.11)

である. ここで ν_{i} はグループ i の数, 例えば aliphatic–CH_{3} 鎖の数である. 代表的なグループパラメタを表 A.1 に示す.

\underline{朧モデルの差分化方法付録 A AMP 式 26}

表 A.1: AMP 式のグループパラメタ (化学工学便覧より)

s_{i}  ε_{0}_{i}  V_{w}_{i}
グループ名  [-]  [K] [cm^{3} mol^{−1}]
aliphatic–CH3  2.359 1162.7  13.67
aromatic AR=C(Ar)H  1.175 939.5  8.06
aromatic AR=C(Ar)R  -0.520 583.0  5.54
condensed aromatic AR=C(Ar)cond  -0.774 432.5  4.74
condensed aromatic AR=C(cond)cond 0.321 623.5  4.74
naphthenic CH_{2} 5–member ring  1.335 985.4  9.30
heterocycloaromatic NH  4.847 4253  8.08
heterocycloaromatic N 6–member ring 2.201 2184  5.2

\underline{朧モデルの差分化方法付録 B 相平衡条件を用いた存在可能な化学種の判定 27}

付録B {相平衡条件を用いた存在可能な化学 \atop 種の判定}

第 1 章で述べたギブスの自由エネルギー最小化法では, 行列式が非正則にならないようにせねばならない. そのため, 行列式を解く前に, 相平衡条件を用いて当該温度圧力において存在可能な化学種を判定しておく必要がある.

温度, 圧力, 組成を与えた場合, 相 ϕ_{1} と相 ϕ_{2} が平衡である条件は以下のように書ける.

μ_{i}^{ϕ_{1}}(T, p, n_{i}^{ϕ_{1}}) = μ_{i}^{ϕ_{2}}(T, p, n_{i}^{ϕ_{2}}).  (B.1)

また化学反応の場合, 平衡条件は化学量論係数 ν_{i} を用いて以下のように書ける.

∑ ν_{i}μ_{i} = 0.  (B.2)

但し, 化学反応の場合添字 ϕ を省略した. それぞれの化学反応によってとりうる相がばらばらであるからである.

知りたい情報は, 温度と圧力を与えた際に存在できない化学種はどれか, ということである. この判定を行うのはギブスの自由エネルギー最小化法を実行する前なので, 組成は定まっていないためである. 以下の節でケース毎に相平衡条件を用いた存在可能な化学種の判定方法を述べる.

B.1 気相と凝縮相の相平衡条件

気相と凝縮相の間での平衡条件は,

μ_{i}^{g}^{∘}(T) + RT ln (\frac{p}{p_{0}} ) + RT ln x_{i}^{g} = μ_{i}^{c}^{∘}(T)  (B.3)

但し気相の添字を g, 凝縮相の添字を c とした. (B.12) 式の RT ln x_{i}^{g} 項の大きさは −∞ < RT ln x_{i}^{g} ≤ 0 なので, 凝縮相の化学ポテンシャルの値は以下の範囲となる.

−∞ < μ_{i}^{c}^{∘}(T) ≤ μ_{i}^{g}^{∘}(T) + RT ln (\frac{p}{p_{0}} )  (B.4)

したがって (B.4) 式より以下が言える.

[figure]

当該温度圧力において,

μ_{i}^{c}^{∘}(T) > μ_{i}^{g}^{∘}(T) + RT ln (\frac{p}{p_{0}} )  (B.5)

の関係式が成立するならば凝縮相は存在できない.

[figure]

B.2 液相と固相の相平衡条件

液相と固相の間での平衡条件は,

μ_{i}^{l}^{∘}(T) = μ_{i}^{s}^{∘}(T)  (B.6)

但し液相の添字を l, 固相の添字を s とした. (B.6) 式は温度のみの関数である. 温度を与えた際に μ_{i}^{l}^{∘}(T) ≠ μ_{i}^{s}^{∘}(T) ならば, より化学ポテンシャルの小さな化学種のみ存在する. したがって (B.6) 式より以下が言える.

[figure]

当該温度において,

μ_{i}^{l}^{∘}(T) < μ_{i}^{s}^{∘}(T)  (B.7)

の関係式が成立するならば液相のみ存在し, 固相は存在できない. その一方,

μ_{i}^{l}^{∘}(T) > μ_{i}^{s}^{∘}(T)  (B.8)

の関係式が成立するならば固相のみ存在し, 液相は存在できない.

[figure]

B.3 気相と溶液の相平衡: 溶液について勉強しなおし!

気相と溶液の間での平衡条件は,

μ_{i}^{g}^{∘}(T) + RT ln (\frac{p}{p_{0}} ) + RT ln x_{i}^{g} = μ_{i}^{s}^{∘}(T) + RT ln x_{i}^{s}

(μ_{i}^{g}^{∘}(T) + RT ln (\frac{p}{p_{0}} )) − μ_{i}^{s}^{∘}(T) = RT ln x_{i}^{g} − RT ln x_{i}^{s}

(B.9)

但し気相の添字を g, 溶液の添字を s とした. (B.9) 式の RT ln x_{i}^{g} 項, RT ln x_{i}^{s} 項の大きさはそれぞれ −∞ < RT ln x_{i}^{g} ≤ 0, −∞ < RT ln x_{i}^{s} ≤ 0 なので, 溶液の化学ポテンシャルの値は以下の範囲となる.

B.4 化学反応の平衡

化学反応の場合を一般的に取り扱うのは難しいので, モデル中で使われている NH_{4}SH の生成反応に関して議論する. NH_{4}SH の生成反応の平衡条件は,

\{μ_{NH_{3}}^{g}^{∘}(T) + RT ln (\frac{p}{p_{0}} ) + RT ln x_{NH_{3}}^{g}\} + \{μ_{H_{2}S}^{g}^{∘}(T) + RT ln (\frac{p}{p_{0}} ) + RT ln x_{H_{2}S}^{g}\}
= μ_{NH_{4}SH}^{s}^{∘}(T) + RT ln x_{NH_{4}SH}^{g}  (B.10)
μ_{NH_{3}}^{g}^{∘}(T) + μ_{H_{2}S}^{g}^{∘}(T) + 2RT ln (\frac{p}{p_{0}} ) + RT ln x_{NH_{3}}^{g} · x_{H_{2}S}^{g} = μ_{NH_{4}SH}^{s}^{∘}(T) + RT ln x_{NH_{4}SH}^{g} (B.11)

(B.12)

\underline{朧モデルの差分化方法付録 C 初期値の作成方法 29}

付録C 初期値の作成方法

C.1 使わなかった方法

• 化学ポテンシャル, 低温と高温

• NH4SH の平衡条件 -¿ どちらかが全然存在しない場合には...

• 前処理