% 表題: ECCM の解法(熱力学ポテンシャル G の最小化法)
%
% 履歴: 2000/08/02 杉山耕一朗
% 履歴: 2002/09/23 杉山耕一朗


\chapter{大気の平衡状態の記述}
\label{method}


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


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

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

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

\begin{eqnarray}
G =\sum \mu_{i} n_{i}.
\Deqlab{2001-02-13-1}
\end{eqnarray}

の最小値を求めれば良い. 

ここで添字 $i$ は大気に含まれる $i$ 番目($0 \leq i \leq N$) の
物質を意味し, 添字 $e$ は系に含まれる $e$ 番目($0 \leq e \leq M$)
の元素を意味する. 
$\mu_{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$ となる. 
\Deqref{1-1} 式はモル数は常に正の数であることを意味する. 
\bigskip

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

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





\section{ギブス自由エネルギーの極値の与え方}
\label{minimum}

\subsection{大気が理想気体のみから成る場合}

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


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


\begin{eqnarray}
Q(n_{i}) = G(m_{i}) + \sum_{i} \left( \DP{G}{m_{i}}\right) \delta m_{i}
+ \Dinv{2} \sum_{i} \sum_{j} \left( \frac{\partial^{2}G}{\partial m_{i}\partial m_{l}}\right) \delta m_{i} \delta m_{j}.  
\Deqlab{3.3.24}
\end{eqnarray}

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

\begin{eqnarray}
 \left(\DP{G}{m_{i}}\right) 
= \mu_{i}
= \mu_{i}^{\circ}(T) 
  + RT \ln{\left( \frac{m_{i}}{\sum m_{i}}\right)}
  + RT \ln\left(\frac{p}{p_{0}}\right)
\Deqlab{3.3.25}
\end{eqnarray}

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

\begin{eqnarray}
\left( \frac{\partial^{2}G}{\partial m_{i}\partial m_{l}}\right)
 = RT \left( \frac{\delta_{il}}{m_{i}} - \Dinv{\sum m_{i}}\right). 
\Deqlab{3.3.27}
\end{eqnarray}

\bigskip


\Deqref{3.3.25}, \Deqref{3.3.27} 式を 
\Deqref{3.3.24} 式に代入し, $RT$ で割ることによって, 

\begin{eqnarray}
\frac{Q(n_{i})}{RT} = \frac{G(m_{i})}{RT} 
 + \sum_{i}  \frac{\mu_{i}}{RT} \delta m_{i}
 + \Dinv{2} \left[ \sum_{i}\frac{(\delta m_{i})^{2}}{m_{i}}
 - \frac{\{\sum (\delta m_{i})\}^{2}}{\sum m_{i}}\right], 
\Deqlab{3.3.28}
\end{eqnarray}

が導かれる. 
\bigskip

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

\begin{eqnarray}
L =
\frac{G(m_{i})}{RT} 
+ \sum_{i} \frac{\mu_{i}}{RT} \delta m_{i}
+ \Dinv{2} \left[ \sum_{i}\frac{(\delta m_{i})^{2}}{m_{i}} 
   - \frac{\{\sum (\delta m_{i})\}^{2}}{\sum m_{i}}\right]  
- \sum_{e} \chi_{e} 
 \left[\sum_{i} a_{i,e}n_{i} - B_{e}\right], \nonumber 
\end{eqnarray}

\Deqref{3.3.28} 式が極値をとる条件は $\partial L / \partial n_{i}^{\phi} = 0$ なので, 

\begin{eqnarray}
 \frac{\mu_{i}}{RT}
+  \frac{n_{i}}{m_{i}} 
- \frac{\sum n_{i}}{\sum m_{i}}
- \sum_{e} \chi_{e} a_{ie} = 0,
\Deqlab{3.2.30}
\end{eqnarray}

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


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

\begin{eqnarray}
 n_{i} = 
  m_{i} \left\{
	 - \frac{\mu_{i}}{RT} 
	 +  \frac{\sum n_{i}}{\sum m_{i}} 
	 + \sum_{e} \chi_{e} a_{ie}
  \right\}. 
\Deqlab{3.3.31}
\end{eqnarray}

この \Deqref{3.3.31} 式を $i$ に関して 1 から $N$ までの和をとると, 

\begin{eqnarray}
 \sum_{e}B_{e}\chi_{e} = \sum_{i} m_{i} \frac{\mu_{i}}{RT},
\Deqlab{3.3.33}
\end{eqnarray}

となる. また \Deqref{3.3.31} 式の両辺に $a_{if}$ をかけて $i$ に関して和を取る.

\begin{eqnarray}
\sum_{e} \left\{ \sum_{i} a_{ie} a_{if} m_{i} \right\} \chi_{e}
 + B_{f} \left( \frac{\sum n_{i}}{\sum m_{i}} - 1 \right) 
 = \sum_{i} a_{if} m_{i} \frac{\mu_{i}}{RT},
\Deqlab{3.3.34-2}
\end{eqnarray}

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

式変形の結果得られた \Deqref{3.3.33}, \Deqref{3.3.34-2} 式は $M + 1$ 個
の方程式から成り, それらを連立することで $\chi_{e}$ と 
$(\sum n_{i}/\sum m_{i} - 1)$ の値
が求まる.  各物質のモル数は得られた値を \Deqref{3.3.31} 式に
代入することで求まる.  \bigskip

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


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

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

\begin{eqnarray}
 n_{i} = 
  m_{i} \left\{
	 - \frac{\mu_{i}}{RT} 
	 + \frac{\sum n_{i}}{\sum m_{i}}
	 + \sum_{e} \chi_{e} a_{ie}
  \right\}. \nonumber
\end{eqnarray}

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


\begin{eqnarray}
&& \sum_{e}B_{e}\chi_{e} = \sum_{i} m_{i} \frac{\mu_{i}}{RT},
\nonumber \\
&& \sum_{e} \left\{ \sum_{i} a_{ie} a_{if} m_{i} \right\} \chi_{e}
 + B_{f} \left( \frac{\sum n_{i}}{\sum m_{i}} - 1 \right) 
 = \sum_{i} a_{if} m_{i} \frac{\mu_{i}}{RT}. \nonumber
\end{eqnarray}

\end{screen}



\subsection{多相系への拡張}

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

まず物質の保存則は以下のように書き換えることができる. 
\begin{eqnarray}
\sum_{\phi}\sum_{i} a_{ie}^{\phi} n_{i}^{\phi} 
= B_{e}
\Deqlab{20010114-1}
\end{eqnarray}

尚, \Deqref{20010114-1} 式の $B_{e}$ は各元素の全相における総和
に拡張されている. 多相系で保存されているのは, 
全相にわたる元素の存在量である. 
\bigskip

次に理想系での化学ポテンシャルの表現を拡張する. 
\begin{eqnarray}
\mu_{i}^{\phi}
= {\mu_{i}^{\circ}}^{\phi}(T) 
+ RT \ln\left({\frac{m_{i}^{\phi}}{\sum_{i} m_{i}^{\phi}}} \right)
+ \alpha RT \ln \left( \frac{p}{p_{0}} \right)
\Deqlab{20010114-2}
\end{eqnarray}

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

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

\begin{eqnarray}
Q(n_{i}^{\phi}) 
= G(m_{i}^{\phi}) 
+ \sum_{\phi} \sum_{i} 
  \underbrace{
    \left( \DP{G}{m_{i}^{\phi}}\right)
  }_{A} 
  \delta m_{i}^{\phi}
+ \Dinv{2} \sum_{\phi} \sum_{i} \sum_{j} 
  \underbrace{
  \left( 
     \frac{\partial^{2}G}{\partial m_{i}^{\phi} \partial m_{l}^{\phi}}
  \right)
  }_{B}  
  \delta m_{i}^{\phi} \delta m_{j}^{\phi}, 
\Deqlab{20010114-3}
\end{eqnarray}

と書ける. 
ここで \Deqref{20010114-3} 式の A 項は \Deqref{20010114-2} 式を用いるこ
とによって, 

\begin{eqnarray}
\left(\DP{G}{m_{i}^{\phi}}\right) 
\equiv \mu_{i}^{\phi}
= {\mu_{i}^{\circ}}^{\phi}(T) 
+ RT \ln{\left( \frac{m_{i}^{\phi}}{\sum_{i} m_{i}^{\phi}}\right)}
+ \alpha RT \ln{\left( \frac{p}{p_{0}}\right)}. 
\Deqlab{20010114-4}
\end{eqnarray}

と書ける. また \Deqref{20010114-3} 式の B 項を 
\Deqref{20010114-2} 式を用いて書き下すと,

\begin{eqnarray}
\left( \frac{\partial^{2}G}{\partial m_{i}^{\phi}\partial m_{l}^{\phi}}\right)
 = RT \left( \frac{\delta_{il}}{m_{i}^{\phi}} 
   - \Dinv{\sum_{i} m_{i}^{\phi}}\right), 
\Deqlab{20010114-6}
\end{eqnarray}

となる. 
\bigskip


\Deqref{20010114-4}, \Deqref{20010114-6} 式を 
\Deqref{20010114-3} 式に代入し, $RT$ で割ることによって, 

\begin{eqnarray}
\frac{Q(n_{i}^{\phi})}{RT} 
= \sum_{\phi}\frac{G(m_{i}^{\phi})}{RT} 
+ \sum_{\phi} \sum_{i} \frac{\mu_{i}^{\phi}}{RT} 
  \delta m_{i}
+ \Dinv{2} \sum_{\phi}
  \left[ 
     \sum_{i}\frac{(\delta m_{i}^{\phi})^{2}}{m_{i}^{\phi}}
     - \frac{\left\{\sum_{i} (\delta m_{i}^{\phi})\right\}^{2}}{\sum_{i} m_{i}^{\phi}}
  \right], 
\Deqlab{20010114-7}
\end{eqnarray}

が導かれる. 
\bigskip

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

\begin{eqnarray}
L =
\sum_{\phi}\frac{G(m_{i}^{\phi})}{RT} 
+ \sum_{\phi}\sum_{i} 
      \frac{\mu_{i}^{\phi}}{RT}
  \delta m_{i}
+ \Dinv{2} \sum_{\phi}
    \left[ 
      \sum_{i}\frac{(\delta m_{i}^{\phi})^{2}}{m_{i}^{\phi}} 
      - \frac{\left\{\sum_{i}(\delta m^{\phi})\right\}^{2}}{\sum_{i} m^{\phi}}
    \right]  
- \sum_{e} \chi_{e} \sum_{\phi}
    \left[
       \sum_{i} a_{ie}^{\phi} n_{i}^{\phi} - B_{e}^{\phi} 
    \right], \nonumber 
\end{eqnarray}

極値をとる条件は $\partial L / \partial n_{i}^{\phi} = 0$ なので, 

\begin{eqnarray}
\frac{\mu_{i}^{\phi}}{RT} 
  + \frac{n_{i}^{\phi}}{m_{i}^{\phi}} 
  - \frac{\sum_{i} n_{i}^{\phi}}{\sum_{i} m_{i}^{\phi}}
  - \sum_{e} \chi_{e} a_{ie}^{\phi} = 0,
\Deqlab{20010114-9}
\end{eqnarray}

である. \Deqref{20010114-9} 式と \Deqref{20010114-1} 式を連立することで, 
変数 $n_{i}^{\phi}$, $\chi_{e}$ (ラグランジュの未定係数)に関する $\Phi N
+ M $ 個の線形方程式が得られ, 方程式系が閉じる.  \bigskip


単相系と同様に \Deqref{20010114-9} 式は式変形を行うことで連立方程式の次
数を減らすことができる.  まず \Deqref{20010114-9} 式を $n_{i}^{\phi}$ に
関する式に書き直す.

\begin{eqnarray}
n_{i}^{\phi} = 
m_{i}^{\phi} \left\{
- \frac{\mu_{i}^{\phi}}{RT}
+ \frac{\sum_{i} n_{i}^{\phi}}{\sum_{i} m_{i}^{\phi}}
+ \sum_{e} \chi_{e} a_{ie}^{\phi}
\right\}.
\Deqlab{20010123-2-1}
\end{eqnarray}

この \Deqref{20010123-2-1} 式を $i$ に関して 1 から $N$ までの和をとると, 

\begin{eqnarray}
\sum_{i} \sum_{e} a_{ie}^{\phi} m_{i}^{\phi} \chi_{e} = 
\sum_{i} m_{i}^{\phi} \frac{\mu_{i}^{\phi}}{RT}
\Deqlab{20010123-3}
\end{eqnarray}

となる. また \Deqref{20010123-2-1} 式の両辺に $a_{if}^{\phi}$ をかけて
 $i$ と $\phi$ に関して和を取る.

\begin{eqnarray}
\sum_{e} 
  \left\{ 
    \sum_{i} \sum_{\phi} a_{ie}^{\phi} a_{if}^{\phi} m_{i}^{\phi}
  \right\} 
  \chi_{e}
+ \sum_{\phi} \sum_{i} a_{if}^{\phi} m_{i}^{\phi}
  \left( 
     \frac{\sum_{i} n_{i}^{\phi}}{\sum_{i} m_{i}^{\phi}}
  \right) 
= \sum_{i}\sum_{\phi} a_{if}^{\phi} m_{i}^{\phi} \frac{\mu_{i}^{\phi}}{RT} + B_{f}
\Deqlab{3.3.34}
\end{eqnarray}

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

\Deqref{20010123-3}, \Deqref{3.3.34} 式は $M + \Phi$ の方程式であり, 
また変数の数は $\chi_{e}$ (M 個) と $\sum_{i} n_{i}^{\phi}$ ($\Phi$ 個)なので, 
この 2 式を連立することによって方程式が閉じる. 
\bigskip


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



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

\begin{eqnarray}
n_{i}^{\phi} = 
m_{i}^{\phi} \left\{
- \frac{\mu_{i}^{\phi}}{RT}
+ \frac{\sum_{i} n_{i}^{\phi}}{\sum_{i} m_{i}^{\phi}}
+ \sum_{e} \chi_{e} a_{ie}^{\phi}
\right\}.
\nonumber
\end{eqnarray}

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

\begin{eqnarray}
&&\sum_{i} \sum_{e} a_{ie}^{\phi} m_{i}^{\phi} \chi_{e} = 
\sum_{i} m_{i}^{\phi} \frac{\mu_{i}^{\phi}}{RT}
\nonumber \\ 
&&\sum_{e} 
  \left\{ 
    \sum_{i} \sum_{\phi} a_{ie}^{\phi} a_{if}^{\phi} m_{i}^{\phi}
  \right\} 
  \chi_{e}
+ \sum_{\phi} \sum_{i} a_{if}^{\phi} m_{i}^{\phi}
  \left( 
     \frac{\sum_{i} n_{i}^{\phi}}{\sum_{i} m_{i}^{\phi}}
  \right) 
= \sum_{i}\sum_{\phi} a_{if}^{\phi} m_{i}^{\phi} \frac{\mu_{i}^{\phi}}{RT} + B_{f}
\nonumber
\end{eqnarray}

\end{screen}



\subsection{負のモル数になった場合には}

\Deqref{20010123-2-1}, \Deqref{20010123-3}, \Deqref{3.3.34} 式を
連立することによって, 平衡状態に近い物質存在量 $n^{\phi}_{i}$ が求まる. 
しかし得られた $n^{\phi}_{i}$ は負になる可能性がある. 
もしも計算から得られた存在量 $n^{\phi}_{i}$ が負の場合, 
以下の式より得られる $n_{i}^{\phi *}$ をより平衡状態に近い物質存在量とみなす. 

\begin{eqnarray}
n_{i}^{\phi *} = m^{\phi}_{i} + \lambda \{n^{\phi}_{i} - m^{\phi}_{i}\}.
\Deqlab{20010129-51}
\end{eqnarray}

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

\begin{eqnarray}
&&\lambda \equiv 
 \frac{m^{\phi}_{i}}{m^{\phi}_{i} - n^{\phi}_{i}} \times 0.999 
\Deqlab{20010214-10}   \;\;\;\; (n_{i}^{\phi} \leq 0),\\
&&\lambda \equiv 1 \;\;\;\; (n_{i}^{\phi} > 0).
\Deqlab{20010214-10-1} 
\end{eqnarray}


なお, \Deqref{20010129-51} 式の $n_{i}^{\phi *}$ は 
元素保存の式 (\Deqref{2} 式) を満たす. 



