% 表題   DCPAM5  力学過程の支配方程式の導出 -- 基礎方程式系の導出
%
% 履歴 
%\Drireki{1994/04/13 石渡正樹}
%\Drireki{1997/04/15 赤堀浩司}
%\Drireki{2008/06/15 森川靖大}
%

\section{基礎方程式系の導出}

方程式系は 6 本の予報方程式と 1 本の診断方程式からなる. 予報方程式は, 
全質量の連続の式, 水蒸気量の式, 運動方程式(3 成分), 熱力学の式からなる. 
これらは, それぞれ, 全質量保存則, 水蒸気量の保存則, 全質量に関する運動
量保存則, 全質量に関する全エネルギー保存則から導出する. 診断方程式には, 
理想気体の状態方程式を用いる
    \footnote{
	乾燥空気と水蒸気は, 同じ速度と温度をもつことを暗黙のうちに仮定
	している. したがって, 水蒸気に関する運動量保存則および全エネルギー保存
	則および状態方程式を考慮する必要がない.
     }.

！！注意: この付録中では導出の都合上, 乾燥空気の気体定数を $R^d$, 
定圧比熱を $C_p^d$, 全大気の気体定数を $R$ とおく. 
しかし, モデルの実装について記した別紙
『\htmladdnormallinkfoot{支配方程式系とその離散化}
{http://www.gfd-dennou.org/library/dcpam/dcpam5/dcpam5\_current/doc/basic\_equations/htm/basic\_equations.htm}』
の『力学過程』
では, 乾燥空気の気体定数を $R$, 定圧比熱を $C_p$ と表記している
ので留意いただきたい. 

\subsection{状態方程式}

乾燥空気, 水蒸気の状態方程式はそれぞれ
    \begin{align}
       p^d & = \rho^{d} R^d T, \\
       p^v & = \rho^{v} R^v T
    \end{align}
である. ここで$p$ は圧力, $\rho$は密度, $R$は気体定数, $T$は温度であり,
$\bullet^d$, $\bullet^v$ はそれぞれ乾燥空気および水蒸気に
関する量であることを示す. したがって, 全圧 $p=p^d+p^v$ は, 
    \begin{align}
     \begin{split}
       p & =  (\rho^d R^d + \rho^v R^v) T \\
        & = \rho R^d ( 1 + \epsilon_v q ) T
     \end{split}    
    \end{align}
となる. ここで, $q=\rho_v/\rho$ は比湿, であり,
$\epsilon_v \equiv 1/\epsilon -1$,
$\epsilon \equiv R^d/R^v$ である. したがって, 
全大気の状態方程式は,
    \begin{align}
       \Deqlab{状態方程式}
       p = \rho R T.
    \end{align}
ただし, $R \equiv R^d ( 1+\epsilon_v q )$ である. あるいは, 仮温度 
$T_v \equiv T ( 1 + \epsilon_v q )$ を用いれば,
    \begin{align}
       \Deqlab{状態方程式仮温度使用}
       p = \rho R^d T_v
    \end{align}
と表される. 

\subsection{連続の式}

全大気の質量保存則は, 水蒸気の生成消滅を無視すれば
    \footnote{
	次で示すように水蒸気式では生成消滅を含めている. したがって,
	全大気の質量保存則は, 水蒸気の生成消滅が起きても全質量が
	保存するように, 乾燥大気量が変化することを
	要請していることになる.
	}, 
    \begin{align}
        \DP{\rho}{t}
      + \DP{}{x_j}( \rho v_j )
      = 0.
      \Deqlab{全大気質量フラックス}
    \end{align}
ここで, $v$は風速である. 
ラグランジュ形式で記述すれば, 
    \begin{align}
        \DD{\rho}{t}
      + \rho \Ddiv \Dvect{v}
      = 0.
    \end{align}

\subsection{水蒸気の式}

水蒸気密度 $\rho^v$ に対する質量保存則は, 単位時間単位体積あたりの生成
消滅量を $S$ とすれば, 
    \begin{align}
     \DP{\rho^v}{t}
      + \DP{}{x_j} ( \rho^v v_j )
      = S.
      \Deqlab{水蒸気質量フラックス}
    \end{align}
比湿 $q=\rho^v/\rho$ に関する式は, 原理的には\Deqref{全大気質量フラッ
クス} と\Deqref{水蒸気質量フラックス} から得ることができる. しかし, 
今の場合, \Deqref{全大気質量フラックス}で水蒸気の生成消滅を無視したので, 
正しくは得られない. そこで比湿の生成消滅に関する項を改めて $S_q$ と定
義する. 
    \begin{align}
      \DD{q}{t} = S_q.
    \end{align}

\subsection{運動方程式}

運動量保存則は, 水蒸気の生成消滅にともなう運動量変化を無視すれば次のよ
うに書ける.
    \begin{align}
        \DP{}{t}(\rho v_i)
      + \DP{}{x_j}( \rho v_i v_j )
      + \DP{p}{x_i} 
      - \DP{\sigma_{ij}}{x_j}
      + \rho \DP{\Phi^*}{x_i}
     = {\cal F}_i^{\prime}.
     \Deqlab{運動量フラックス}
    \end{align}
%
ここで, $\sigma_{ij}$ は粘性応力テンソル, $\Phi^*$ は惑星
の引力によるポテンシャル
    \footnote{
        これは遠心力を考慮しない惑星の質量にのみ起因
        したポテンシャル.},
${\cal F}_i^{\prime}$ はその他の外力項である. あるいは連続の式
を用いてラグランジュ形式で記述すると
    \begin{align}
        \rho \DD{v_i}{t}
      + \DP{p}{x_i}
      - \DP{\sigma_{ij}}{x_j}
      + \rho \DP{\Phi^*}{x_i}
     = {\cal F}_i^{\prime }
    \end{align}
となる. ここで, 粘性項と外力項を ${\cal F}_i$ とおき, さらにベクトル表示する.
    \begin{align}
        \rho \DD{\Dvect{v}}{t}
      + \Dgrad p
      + \rho \Dgrad \Phi^*
     = \Dvect{\cal F}.
    \end{align}

\subsection{熱力学の式}

単位質量あたりの全エネルギーは, 運動エネルギー $\Dvect{v}^2/2$ と内部エ
ネルギー $\varepsilon$ およびポテンシャルエネルギー $\Phi^*$ の和で表
現される. この時間変化率の式は, 水蒸気の生成消滅による影響を無視すれば,
    \begin{align}
            \DP{}{t} 
             \left[ \rho 
                   \left(   \frac{1}{2} \Dvect{v}^2 
                 + \varepsilon + \Phi^* \right) \right]
         + \DP{}{x_j} \left[ 
             \rho 
               \left(   \frac{1}{2} \Dvect{v}^2 
                    + \varepsilon + \Phi^* \right)v_j
                    + p v_j - \sigma_{ij}v_i  
            \right]
      =  \rho Q + {\cal F}_i^{\prime} v_i
      \Deqlab{全エネルギーフラックス}
    \end{align}
である. ここで, $Q$ は外部からの加熱率である. 一方, 運動エネルギーとポ
テンシャルエネルギーの和の保存式は, 運動量保存式 \Deqref{運動量フラック
ス} に $v_i$ をかけ, 連続の式を用いて変形することで得られる
    \footnote{
        導出の過程を示す. 左辺第1項と第2項は次のように変形される. 
            \begin{align}
              v_i \DP{}{t} ( \rho v_i ) 
                     + v_i \DP{}{x_j} ( \rho v_j v_i )
              & =   \DP{}{t} ( \rho v_i^2 ) 
                     + \DP{}{x_j} ( \rho v_j v_i^2 )
                     - \rho \DP{}{t} \left( \frac{1}{2}  v_i^2 \right) 
                     - \rho v_j \DP{}{x_j} \left( \frac{1}{2} v_i^2 \right)
                      \nonumber \\
              & =   \DP{}{t} ( \rho v_i^2 )
                     + \DP{}{x_j} ( \rho v_j v_i^2 )  
                     - \DP{}{t} \left( \frac{1}{2}  \rho v_i^2 \right)  
                     - \DP{}{x_j} \left( \frac{1}{2} v_i^2 \rho v_j \right)
                       \nonumber \\ 
              & \quad
                     + \frac{1}{2} v_i^2 \DP{\rho}{t}
                     + \frac{1}{2} v_i^2 \DP{}{x_j} ( \rho v_j )  \nonumber \\
              & =   \DP{}{t} \left( \frac{1}{2} \rho v_i^2 \right)  
                     + \DP{}{x_j} ( \frac{1}{2} \rho v_j v_i^2 )
                     + \frac{1}{2} v_i^2  
                       \left\{ \DP{\rho}{t} + \DP{}{x_j} ( \rho v_j ) \right\}
                      \nonumber  \\
              & =   \DP{}{t} \left( \frac{1}{2} \rho v_i^2 \right)
                     + \DP{}{x_j} ( \frac{1}{2} \rho v_j v_i^2 ).  \nonumber 
            \end{align}
        また, 左辺第5項は次のように変形される.
	変形の際には $\DP{\Phi^*}{t}=0$ であるとしている. 
            \begin{align}
              v_i \rho \DP{\Phi^*}{x_i} 
               & = \Phi^* \left\{ \DP{\rho}{t} + \DP{}{x_i}(\rho v_i) \right\} 
                      + \rho \DP{\Phi^*}{t}
                      + v_i \rho \DP{\Phi^*}{x_i} \nonumber \\
               & = \DP{}{t} ( \rho \Phi^* )
                      + \DP{}{x_i} ( \rho \Phi^* v_i ). \nonumber 
            \end{align}
    }.
    \begin{align}
       \Deqlab{運動とポテンシャルエネルギーフラックス}
          \DP{}{t} \left( \frac{1}{2} \rho v_i^2 + \rho \Phi^* \right) 
        + \DP{}{x_j} \left( \frac{1}{2} \rho v_j \Dvect{v}^2
                            + \rho \Phi^* v_j
                            + p v_j - \sigma_{ij} v_i \right)
     =  p \DP{v_j}{x_j} - \sigma_{ij} \DP{v_i}{x_j} + {\cal F}_i^{\prime} v_i .
    \end{align}
ここで, 変形の際には $\DP{\Phi^*}{t}=0$ であるとしている. 
\Deqref{全エネルギーフラックス}と
\Deqref{運動とポテンシャルエネルギーフラックス}
との差をとると, 次のように内部エネルギーの式が得られる.
    \begin{align}
       \DP{}{t} ( \rho \varepsilon )
         + \DP{}{x_j} ( \rho \varepsilon v_j )
       =  - p \DP{v_j}{x_j} + \sigma_{ij} \DP{v_i}{x_j}
         + \rho Q .
    \end{align}
連続の式を用いてラグランジュ形式に書き直せば
    \begin{align}
     \Deqlab{内部エネルギー}
       \rho \DD{\varepsilon}{t} 
      =   \frac{p}{\rho} \left( \DD{\rho}{t} \right)
        + \rho Q.
    \end{align}
以降では, 外部からの加熱の項と粘性による加熱の項を
まとめて $Q^*$ とおくこととする.

内部エネルギーを温度を用いて表現すると $\varepsilon = C_v T$ である.
$C_v$は定圧比熱である. 
さらに状態方程式 \Deqref{状態方程式} を用いて\Deqref{内部エネルギー}
を変形する. $C_p = C_v + R$ であることに注意すれば
    \begin{align}
        \DD{C_p  T}{t} = \frac{1}{\rho} \DD{p}{t} + Q^*,
    \end{align}
となる. ここで, $C_p$ を乾燥空気の定圧比熱 $C_p^d$ 
で近似すると
    \footnote{
        この近似には疑問が残る. 状態方程式においては, 
	気体定数 $R$ を $R^d$ とする近似は
	(仮温度 $T_v$ を導入することで)行なわなかった. 
	$C_p$ についてだけ近似するのは近似のレベルに
	一貫性がないように思われる. 
	
        以下はその主張. 混合比 $r=\rho^v/\rho^d$ を用いている. 全大気の内部
        エネルギーは 
        \begin{align*}
          \rho \varepsilon 
         & = \rho^d \varepsilon^d + \rho^v \varepsilon^v  \\
         & = \rho^d C_v^d T + \rho^v C_v^v T  \\
         & = \rho \left( \frac{ \rho^d C_v^d
             + \rho^v C_v^v}{\rho} \right) T  \\
         & = \rho \left( \frac{ C_v^d + r C_v^v }{ 1+r } \right) T,  
        \end{align*}
        となる. したがって, 
        \[
          C_v \equiv \frac{ C_v^d + r C_v^v }{ 1+r },
        \]
        である. また, 
        \[
          R \equiv \frac{R^d + r R^v}{1+r}
            \left( = R^d \frac{ 1 + r/\epsilon}{1+r} \right) ,
        \]
        であるから, 
        \begin{align*}
         C_p & = C_v + R  \\
             & = \frac{ C_v^d  + r C_v^v + R^d + r R^v }{1+r}  \\
             & = \frac{ C_p^d + r C_p^v}{1+r}  \\
             & = C_p^d \frac{ 1+rC_p^v/C_p^d}{1+r}  \\
             & \sim C_p^d \frac{ 1+ 8 r/7 \epsilon}{1+r}, 
        \end{align*}
        となる. ここで, $C_p^d = ( C_v^d+R^d ) \sim ( \frac{5}{2}R^d +R^d ) =
        \frac{7}{2} R^d $,  および $C_p^v = ( C_v^v + R^v ) \sim ( 3R^v + R^v
        ) = 4 R^v$ を用いた. 熱力学の式では, この状況に対して, $C_p \sim
        C_p^d$ と近似した. しかし, 静力学平衡の式では, たった $8/7$ の違いなの
        に $R$ を $R^d$ に近似せず, 仮温度の導入により厳密に
	取り扱おうとしている. 
    }, 
次の熱力学の式を得る. 
\begin{align}
 \DD{T}{t}  =  \frac{1}{C_p^d \rho} \DD{p}{t} + \frac{Q^*}{C_p^d}.
\end{align}