% 表題: oboro モデルの概要
%
% 履歴: 2001/01/28 杉山耕一朗


\chapter{計算手法の検討}


本章では大気の断熱温度減率と凝縮物質の鉛直分布を計算するための
熱力学計算手法に関して検討する. 
大気中で生じる反応式を考えずに済み, 
さらに大気組成を簡単に変更することのできる計算手法を開発するには, 

\begin{itemize}
\item 大気の熱力学変数として何を選択し, そのとき大気の状態はどの熱力学
      関数によって表現するのが適切か. 
\item 大気の断熱変化をどのように表現するか
\end{itemize}

の 2 点を検討することが必要である. 
以下の節ではこの 2 点に関して検討する. 


\section{全体の枠組}

本研究では


本モデルの全体の枠組は先に紹介した Weidenschilling and
Lewis (1973) の鉛直 1 次元の雲凝結平衡モデルに準じる:
\Dfigref{20010128-10} のように, 圧力 $p$ において温度 $T$, 圧力 $p$, 物
質量 $n_{i}$ の空気塊を考える. ここで添字 $i$ は空気塊に含まれるそれぞれ
の物質を意味する. 空気塊の熱化学平衡状態とそれに含まれる元素の保存を仮定
し, 空気塊を断熱的に上昇させた(圧力を変化させた)際に実現する熱力学的状態を
計算する. 熱力学的な状態量を計算することによって, ある圧力面高度における
空気塊の温度とそこに含まれる物質の存在量が求まる. 
\bigskip

\begin{figure}[h]
\begin{center}
 \Depsf{./images/kei.ps}
\end{center}
\caption{対流熱力学モデルの系の設定. 圧力 $p$ の地点に存在する空気塊を
断熱的に圧力 $p + dp$ の地点まで上昇させる. }
\Dfiglab{20010128-10}
\end{figure}



このモデルで考えている空気塊において, 以下の事項が成立しているものと仮定
している. 

\begin{screen}
\begin{description}
\item[1.] 考える空気塊は熱的・化学的に平衡状態にある. 
\item[2.] 空気塊の温度は断熱的に変化する. 
\item[3.] 空気塊は孤立系であり, 外界から空気塊へ物質が輸送されない. 
\item[4.] 空気塊に含まれる物質は全て理想的であるとする. 気相においては
      理想気体の方程式が成立し, 溶液は理想溶液であるとする. 
\end{description}
\end{screen}
\smallskip

 1 番目の仮定は空気塊の熱力学的状態に関する仮定である. 
 2 番目の仮定は空気塊の断熱上昇を考えていることを意味する. 
3 番目の仮定は空気塊に含まれる物質が外界の影響を受けないことを意味してい
る. 本モデルにおいて凝結した物質は空気塊から離脱すること, つまり降水が生
じることを考えているが, 外界から系に物質が輸送されることを考えていない.
4 番目の仮定は物質の状態方程式に関する仮定である.  現在のところ物質の非
理想性は全く考慮されていないが, 将来的に非理想性を考慮することは可能であ
る. 
\bigskip


本モデルの計算の核をなすのは, 

\begin{itemize}
\item 空気塊の平衡状態をいかに計算するか.
\item 断熱の条件をどのようにモデルに組み込むか. 
\end{itemize}

である. 以下の節ではこの 2 点をどのように扱うか議論する. 




\section{大気の平衡状態の記述}

大気の状態を表現するための熱力学変数として温度 $T$, 圧力 $p$, 物質存在量
$n_{i}$ を選択する. 大気を考える上でこれらの熱力学変数が最も扱い易いから
である. このとき大気の状態を与える適切な熱力学関数は $G$ である. 
\bigskip

大気の状態を熱力学関数 $G$ で表現することに決めたので, 大気の平衡状態は
熱力学関数 $G$ の最小化された状態である. 
温度・圧力を与えると熱力学関数 $G$ は, 

\begin{eqnarray}
G =\sum \mu_{i} n_{i}, 
\Deqlab{1}
\end{eqnarray}

と書ける. ここで添字 $i$ は各物質を意味する添字, $\mu_{i}$ は化学ポテンシャ
ルである. \Deqref{1} 式より熱力学関数 $G$ の最小値は, 物質存在量 $n_{i}$
をパラメータとして変化させることによって求めることができる. 
したがって大気の平衡状態は

\begin{itemize}
 \item 予め与えられた大気の温度 $T$ ・圧力 $p$
 \item \Deqref{1} 式の $G$ を最小化する物質存在量 $n_{i}$
\end{itemize}

によって記述される. また平衡時の熱力学関数 $G$ の値も求まる. 
\bigskip


大気の状態を熱力学関数 $G$ によって表現することにしたので, 潜熱・反応熱
の代わりに化学ポテンシャル $\mu_{i}$ のリストを与えればよい.  化学ポテン
シャルは既存の物性データから求めることができる. 今や Atreya and Romani
(1985) のように, 大気中で生じる全ての化学反応を知っている必要がなくなっ
た.
\bigskip



\section{大気の断熱変化の記述}


大気の断熱変化は熱力学関数 $S$ の保存として記述することができる. 
熱力学関数 $S$ は $G$ を用いて以下のように与えられる. 

\begin{eqnarray}
\left( \DP{G}{T} \right)_{p, n_{i}} = -S
\end{eqnarray}


である. 
温度・圧力空間で $dS = 0$ の曲線の通る軌跡を順にたどれば, 
断熱温度減率と凝縮物質の存在量を求めることができる. 

\section{過去の熱力学モデルとの比較}

外惑星大気対流圏の相変化や化学反応を計算した例として, 
 Weidenschilling and Lewis (1973), Atreya and Romani (1985) がある. 
彼らは熱化学平衡と物質の保存, 物質の理想性を仮定し, 
空気塊を断熱上昇させた際に凝縮する物質の種類とその量を求めた. 
彼らは空気塊の平衡状態を内部エネルギーの保存 dU = 0 から
計算している. 従って自動的に断熱の条件は満たされている. 
\Dfigref{20010128-19} に彼らの方法の模式図を示す. 
\bigskip

例として Atreya and Romani (1985) の定式化を挙げる. 
 1 モルの理想気体が断熱膨張する際のエネルギー保存則は, 

\begin{eqnarray}
\bar{c}_{\rm p} dT - v dP + \sum_{\rm i} L_{\rm i} dn_{\rm i} 
 + \sum L_{\rm RX} dn_{\rm RX} = 0 
\Deqlab{e12.1}
\end{eqnarray}

と書ける. 
ここで $\bar{c}_{\rm p}$ は定圧モル比熱の平均値,
 $dT$ は絶対温度の変化, $v$ は気体分子 1 モル当たりの体積, 
 $dP$ は全圧の変化, $L_{\rm i}$ は i${}_{\rm th}$ 番目の気体 1 モルの
凝縮のエンタルピー, $dX_{i}$ は凝縮による i${}_{\rm th}$ 番目の
気体のモル数の変化, $L_{\rm RX}$ は化学反応の反応熱, $dX_{\rm RX}$ は
化学反応に伴うモル数の変化である. 
\bigskip

 \Deqref{e12.1} 式から断熱温度減率や凝縮量を求めるためには, 
 \Deqref{e12.1} 式中の $dn_{\rm i}$ や $dn_{\rm RX}$ を
$dT$ と $dp$ (静水圧平衡の成立を仮定するならば, $dT$ と $dZ$) 
で書き換える必要がある. $dn_{\rm i}$ はクラジウス--クラペイロンの式, 
$dn_{\rm RX}$ は化学反応の圧平衡定数を用いて書き換えることが
可能である. 例えば化学反応として, 

\begin{eqnarray}
&&{\rm NH_{3}(g)} + {\rm H_{2}S(g)}  \rightarrow {\rm NH_{4}SH(s)}, 
\Deqlab{e12.2} \\
&& \left\{
\ln(K_{\rm p}) = \ln (P_{\rm NH_{3}} \cdot P_{\rm H_{2}S}) 
= 6.467 - \frac{10834}{T} 
\right\} \nonumber
\end{eqnarray}

だけを考えた場合でも, \Deqref{e12.1} 式は


\begin{eqnarray}
\frac{dT}{dZ} = - \frac{\bar{m}g}{\bar{C_{\rm p}}}
\frac
{\left[
 1 + \Dinv{RT} 
\left\{ 
       \sum_{\rm i} L_{\rm i} n_{\rm i} + 
       \frac{2 (n_{\rm H_{2}S} \cdot n_{\rm NH_{3}}) L_{\rm RX}}
       {(n_{\rm H_{2}S} + n_{\rm NH_{3}})}
\right\}
\right]}
{\left[ 
1 + \Dinv{\bar{c}_{\rm p}T^{2}}
\left\{ 
       \sum_{\rm i} \frac{L_{\rm i}^{2} n_{\rm i}}{R} + 
       \frac{ (n_{\rm H_{2}S} \cdot n_{\rm NH_{3}}) }
       {(n_{\rm H_{2}S} + n_{\rm NH_{3}})} 
       \cdot L_{\rm RX} \cdot 10834 
\right\} 
\right] }.
\Deqlab{e12.16}
\end{eqnarray}

となり, 非常に繁雑になる. 
\bigskip


Atreya and Romani (1985) の \Deqref{e12.16} 式を解くためには, 

\begin{screen}
\begin{description}
\item[1.] 空気塊の温度
\item[2.] 物質の凝縮のエントロピー(潜熱).
\item[3.] 化学反応の反応熱と平衡定数.
\item[4.] 物質の定圧モル比熱.
\end{description}
\end{screen}

を知っている必要がある. 
\bigskip

対流熱力学モデルと Atreya and Romani (1985) の方法の違いの 1 つは, 
その定式化において利用する物性値が異なることである. 前者は
化学ポテンシャルであり, 後者は潜熱と比熱と反応熱である. 
しかし物質の化学ポテンシャルを計算する際には潜熱や比熱の
データを使うので, その意味ではこの違いは本質的でない. 
最も大きな違いは式の形式にある. \Deqref{2}, \Deqref{2.3.3} 式において
空気塊に含まれる物質を取り換える場合, 定数である化学ポテンシャルの項だけ
必要に応じて取り換えれば良い. その一方で
\Deqref{e12.1} 式は化学反応式を陽に取り込んでいるので, 
空気塊に含まれる物質を取り換える場合は空気塊内部で生じうる
化学反応式を全て洗い出し, 再度 \Deqref{e12.1} 式を書き下して
 \Deqref{e12.16} 式のような形式に直さねばならない. 
それはとりもなおさず空気塊に含まれる物質を変える時には
数値計算のコードを大幅に書き換える必要があることを意味する. 
\bigskip

確かに Atreya and Romani (1985) のモデルの方が, 大気組成やそこで
生じる化学反応が全て分かっているならば, 
大気の断熱温度減率と物質存在量変化の計算としてはより直接的である. 
しかし外惑星大気では, そもそも元素組成が明らかにされていないので, 
\Deqref{e12.1} 式を利用するのは好ましくないと考える. 
 \Deqref{2}, \Deqref{2.3.3} 式のように, 物質自体を簡単に変化させる
ことのできるモデルの方がより好ましいと考える. 
\bigskip


本モデルは今までの平衡雲凝縮モデルに比べて, 空気塊内の物質自体を簡単に変
更することができるという特徴を持つ. またこの特徴ゆえに, 外惑星以外の大気, 
例えば火星, 金星のような他の惑星にも簡単に応用することが可能であると期待
される.







\begin{figure}
\begin{center}
 \Depsf{./images/map.ps}
\end{center}
\caption{Atreya and Romani (1985) の計算方法の概要}
\Dfiglab{20010128-19}
\end{figure}




数学的に熱力学関数 $G$ を最小化する操作を考える. 
熱力学ポテンシャル $G$ の自然な変数は温度 $T$, 圧力 $p$, 物質の存在量
$n_{i}$ である.  そのうち $T$, $p$ を一定とすると, 

\begin{eqnarray}
G(n_{i}) =\sum \mu_{i} n_{i} 
\Deqlab{1}\\
(但し T, p 一定) \nonumber
\end{eqnarray}

と書ける. ここで添字 $i$ は系に含まれる $i$ 番目($0 \leq i \leq N$)
の物質を意味し, また $\mu_{i}$ は $i$ 番目の物質の化学ポテンシャルを
意味する. 
また系に含まれる元素の保存という条件は, 

\begin{eqnarray}
\sum_{i} a_{ie} n_{i} = B_{e},
\Deqlab{2} \\
n_{i} \geq 0 
\Deqlab{1-1} 
\end{eqnarray}

と書くことができる. ここで添字 $e$ は系に含まれる $e$ 番目($0 \leq e \leq M$)
の元素を意味し, 
また $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$ となる. 
尚, \Deqref{2} 式は全部で $M$ 本の式から成る. 
また \Deqref{1-1} 式は, モル数は常に正の数であることを意味する. 
\bigskip


$G$ の最小値を求めるためには, \Deqref{2} 式の条件下で
 \Deqref{1} 式の極値を求めれば良いので, 
ラグランジュの未定係数法を用いて極値の条件を求めると,  

\begin{eqnarray}
&&\frac{\partial}{\partial n_{i}} 
\left\{
   \sum_{i} \mu_{i} n_{i} - \sum_{e} \chi_{e} 
   \left( \sum_{i}a_{i e}n_{i} - B_{e} \right) 
\right\} = 0, \nonumber \\
&&\mu_{i} - \sum_{e} \chi_{e} a_{ie} = 0.
\Deqlab{2.3.3}
\end{eqnarray}

と書ける. 
最終的に \Deqref{2}, \Deqref{2.3.3} 式を連立することで, 変数 $n_{i}$,
$\chi_{e}$ について解くことができる.  それらの変数が求まることによって
熱力学関数 $G$ の値も求めることができる. 
 \bigskip

