\chapter{モデル計算でのテクニック}

モデルを安定して稼働させるために, いくつかテクニックを使っている. 
本章ではそのテクニックを紹介する. 



\section{負のモル数の処理}

空気塊の平衡状態は 
\Deqref{20010123-2}, \Deqref{20010123-3}, \Deqref{3.3.34} 式を
連立して求めるが, 求まる解 $n_{i}$ は負になる可能性がある. 
そのため, 以下のような処理をして正の数に直している. 
\bigskip

もしも計算から得られた存在量 $n_{i}$ が負である場合, 
平衡時の物質存在量 $n_{i}^{*}$ は, 計算前の存在量 $m_{i}$ と計算から
得られた存在量 $n_{i}$ を用いて以下の式から求める. 

\begin{eqnarray}
n_{i}^{*} = m_{i} + \lambda \{n_{i} - m_{i}\} \hspace{5mm} (>0). 
\Deqlab{20010129-51}
\end{eqnarray}

ここで $\lambda$ は定数である.  \Deqref{20010129-51} 式の $n_{i}^{*}$ は 
元素保存の式 (\Deqref{2} 式) を満たすようにつじつまを合わせてい
る. $n_{i}^{*}$ が正であるという条件を$\lambda$ を用いて書き直すと,


\begin{eqnarray}
\lambda < \frac{m_{i}}{m_{i} - n_{i}}, 
\Deqlab{20010129-52}
\end{eqnarray}


\begin{figure}[h]
\begin{center}
\Depsf{./images/ni.ps}
\end{center}
\caption{負のモル数回避の手順. 次の計算でモル数が負になった場合は, 
値が正になるように}
\end{figure}

となる. 計算によって負のモル数が生じた場合, \Deqref{20010129-52} 式より
それぞれの物質に関する $\lambda_{i}$ を計算し, $\lambda_{i}$ の
最小値を定数 $\lambda$ として \Deqref{20010129-51} 式に代入する. 
それによって平衡時の物質存在量 $n_{i}^{*}$ を求める. 
但し, $\lambda$ をそのまま代入すると $n_{i}^{*} = 0$ となってしまう
物質が現れれるので, 実際には $\lambda$ を 0.99 倍したものを 
\Deqref{20010129-51} 式に代入している. 
\bigskip

もしも $n_{i}$ が負になる物質が存在しないならば, $\lambda$ に 1 を
セットして \Deqref{20010129-51} 式に代入すれば良い. 
そうすることで, $n_{i}^{*} = n_{i}$ となる. 



\section{計算の収束条件}

本モデルでは \Deqref{20010123-2}, \Deqref{20010123-3}, \Deqref{3.3.34} 
式の数値解法として解の反復改良を用いているので, 
計算が収束したか, 即ち空気塊が平衡に達したかどうかの判断を行い, 
計算を止める必要がある. 
本モデルでの収束判定は熱力学関数 $G$ の変化率から行う. 

\begin{itemize}
  \item ループの前と後の熱力学関数 G を比較し, 前の値が後の値よりも小さ
	      い場合. 
  \item ループの前と後の熱力学関数 G の変化率が $10^{-15}$ 以下になった
	      場合. 
\end{itemize}

以上の条件のどちらか一方が満たされた場合, 空気塊は平衡状態にあると
みなし, 反復計算を止める. 
\bigskip


\section{相の削除}


ある相に存在する物質存在量が全て 0 に近付くと, \Deqref{20010123-3} 式は

\begin{eqnarray}
0 \cdot \chi_{e} = 0
\end{eqnarray}

となり, $\chi_{e}$ は不定となる. それを避けるために, 
相に存在する全ての物質が 0 に近付いた場合は, その相が存在しないものとし
て計算を進める. 例えば空気塊に気相・液相が存在すると場合を考える. 
温度・圧力条件によっては液相の物質がどれも不安定になり, どれも
存在しない可能性も十分にある. そのような場合は, 液相が存在しないものと
みなし, つまり変数 $\phi$ を 1 つ減らし, 気相だけで計算を行う. 


\section{$\ln(m_{i}^{\phi}/m^{\phi})$ の扱い}

\Deqref{20010114-10} 式
\begin{eqnarray}
\Gamma_{i}^{\phi}
= m_{i}^{\phi} \left( \frac{{\mu_{i}^{\circ}}^{\phi}}{RT} 
	      + \ln{\frac{m_{i}^{\phi}}{m^{\phi}}} + \alpha \ln{p} \right)
\nonumber
\end{eqnarray}

の対数の扱いにも気をつける必要がある. 物質量 $m_{i}$ が 0 に近付いても, 
$\Gamma_{i}^{\phi}$ が $-\infty$ に発散してしまわないようにする必要が
ある. そのため実際の計算では, 
\Deqref{20010114-10}, \Deqref{20010123-2}, \Deqref{20010123-3}, 
\Deqref{3.3.34} 式の $\Gamma_{i}^{\phi}$ に関する項を, 
$m^{\phi}$ で割った形式で置き換えている. 

\begin{eqnarray}
\frac{\Gamma_{i}^{\phi}}{m^{\phi}}
= \frac{m_{i}^{\phi}}{m^{\phi}}  
  \left( \frac{{\mu_{i}^{\circ}}^{\phi}}{RT} 
  + \ln{\frac{m_{i}^{\phi}}{m^{\phi}}} + \alpha \ln{p} \right). 
\Deqlab{200010131-3}
\end{eqnarray}

\Deqref{200010131-3} 式の利点は, 対数に関する項が
\begin{eqnarray}
\frac{m_{i}^{\phi}}{m^{\phi}}  \cdot
 \ln{\frac{m_{i}^{\phi}}{m^{\phi}}}
\end{eqnarray}

となっているので, 対数の中身 $m_{i}^{\phi}/m^{\phi}$ が 0 に近付いても
全体としては 0 に近付くだけで, $-\infty$ に発散しない. 


