% 表題  移流方程式の差分解法: MPDATA スキーム  一般形
%
% 履歴  1997/11/25 小高正嗣
%       1998/01/20 小高正嗣
%       1998/05/11 小高正嗣
%
\section{一般的な場合の MPDATA}

     ここではより一般的な場合, すなわち N 次元, 速度場が発散的, 速度場が
     非定常の各場合の MPDATA について記述する.

     \subsection{N 次元の場合}

     第 $I$ 次元方向の座標, 速度をそれぞれ $x^{I},u^{I}$ とする.
     2 次元の場合の結果から考えれば $\Dvect{i}$ での $\psi $ の時間変化は 

     \begin{eqnarray}
       \left.\DP{\psi }{t}\right|_{\Dvect{i}}^{n}
       &=& - \sum _{I=1}^{N}\left.\DP{}{x^{I}}(u^{I}\psi )
       \right|_{\Dvect{i}}^{n} \nonumber \\
       && + \sum _{I=1}^{N}\DP{}{x^{I}}\left.\left[\frac{1}{2}
           (|u^{I}|\Delta x^{I}-\Delta t(u^{I})^{2})
           \DP{\psi }{x^{I}} 
           - \sum _{J\neq I}^{N}\frac{1}{2}\Delta t u^{I}u^{J}
           \DP{\psi }{x^{J}}\right]\right|_{\Dvect{i}}^{n}  \Deqlab{eq: 35}
     \end{eqnarray}

     となる. これより反拡散速度は

     \begin{equation}
       u_{d,\Dvect{i}}^{I} = \left.
           0.5[|u^{I}\Delta x^{I}-\Delta t(u^{I})^{2}]\frac{1}{\psi }
           \DP{\psi }{x^{I}}\right|_{\Dvect{i}}
           - \left.\sum _{J\neq I}^{N}0.5\Delta tu^{I}u^{J}\frac{1}{\psi }
           \DP{\psi }{x^{J}}\right|_{\Dvect{i}},
         \Deqlab{eq: 36}
     \end{equation}

     と与えられる.

     \subsection{速度場が発散場の場合}

     速度場が発散場の場合, 例えば 2 次元では\Deqref{eq: 18}式の右辺第 5 項
     の展開を以下のようにする.

     \begin{eqnarray*}
       \DP[2]{\psi }{t} 
       &=& - \DP{}{x}\left
         [u \left( -\DP{}{x}(u\psi ) -\DP{}{y}(v\psi ) \right) \right]
       - \DP{}{y}\left
         [v \left( -\DP{}{x}(u\psi ) -\DP{}{y}(v\psi ) \right) \right] \\
       &=& \DP{}{x}\left[  u^{2} \DP{\psi }{x} + uv \DP{\psi }{y}\right]
       + \DP{}{y}\left[  v^{2} \DP{\psi }{y} + uv \DP{\psi }{x}\right] \\
       && +\DP{}{x}\left[u\psi \left(\DP{u}{x} + \DP{v}{y}\right)\right]
       + \DP{}{y}\left[ v\psi \left(\DP{u}{y} + \DP{v}{x}\right)\right]
     \end{eqnarray*}

     よって\Deqref{eq: 19}式は以下のように変更される.

     \begin{eqnarray}
       \left.\DP{\psi }{t}\right|_{i,j}^{n}
       &=& -\left.\DP{}{x}(u\psi )\right|_{i,j}^{n} 
           -\left.\DP{}{y}(v\psi )\right|_{i,j}^{n} \nonumber \\
       && + \DP{}{x}\left.\left[\frac{1}{2}(|u|\Delta x-\Delta t u^{2})
           \DP{\psi }{x} - \frac{1}{2}\Delta tuv\DP{\psi }{y}\right]
       \right|_{i,j}^{n} \nonumber \\
       && + \DP{}{y}\left.\left[\frac{1}{2}(|v|\Delta y-\Delta t v^{2})
           \DP{\psi }{y} - \frac{1}{2}\Delta tuv\DP{\psi }{x}\right]
       \right|_{i,j}^{n} \nonumber \\
       && - \frac{1}{2}\left.\DP{}{x}\left[\Delta t u\psi 
         \left(\DP{u}{x} + \DP{v}{y}\right)\right]\right|^{n}_{i,j}
       - \frac{1}{2}\left.\DP{}{y}\left[\Delta t v\psi 
         \left(\DP{u}{x} + \DP{v}{y}\right)\right]\right|^{n}_{i,j}.
       \Deqlab{eq: 37}
     \end{eqnarray}

     反拡散速度は

     \begin{eqnarray}
       \tilde{u} &\equiv &\frac{1}{\psi }
       \left[\frac{1}{2}(|u|\Delta x-\Delta t u^{2})
           \DP{\psi }{x} - \frac{1}{2}uv \DP{\psi }{y}\right]
         - \frac{1}{2}\Delta t u\left(\DP{u}{x} + \DP{v}{y}\right),
         \Deqlab{eq: 38} \\
       \tilde{v} &\equiv &\frac{1}{\psi }
       \left[\frac{1}{2}(|v|\Delta y-\Delta t v^{2})
           \DP{\psi }{y} - \frac{1}{2}uv \DP{\psi }{x}\right]
         - \frac{1}{2}\Delta tv\left(\DP{u}{x} + \DP{v}{y}\right),
         \Deqlab{eq: 39}
     \end{eqnarray}

     となる. N 次元では

     \begin{eqnarray}
       \left.\DP{\psi }{t}\right|_{\Dvect{i}}^{n}
       &=& - \sum _{I=1}^{N}\left.\DP{}{x^{I}}(u^{I}\psi )
       \right|_{\Dvect{i}}^{n} \nonumber \\
       && + \sum _{I=1}^{N}\DP{}{x^{I}}\left.\left[\frac{1}{2}
           (|u^{I}|\Delta x^{I}-\Delta t(u^{I})^{2})
           \DP{\psi }{x^{I}} 
           - \sum _{J\neq I}^{N}\frac{1}{2}\Delta t u^{I}u^{J}
           \DP{\psi }{y}\right]\right|_{\Dvect{i}}^{n} \nonumber \\
       && - \sum _{I=1}^{N}\DP{}{x^{I}}\left. \left[\left(
             \frac{1}{2}\Delta t u^{I}\sum _{J=1}^{N}\DP{u^{J}}{x^{J}}\right)
           \psi \right]\right|_{\Dvect{i}}^{n},  \Deqlab{eq: 40}
     \end{eqnarray}

     となって, 反拡散速度は

     \begin{equation}
       u_{d,\Dvect{i}}^{I} = \left.
           0.5[|u^{I}\Delta x^{I}-\Delta t(u^{I})^{2}]\frac{1}{\psi }
           \DP{\psi }{x^{I}}\right|_{\Dvect{i}}
           - \left.\sum _{J\neq I}^{N}0.5\Delta tu^{I}u^{J}\frac{1}{\psi }
           \DP{\psi }{x^{J}}\right|_{\Dvect{i}}
           - \left. 0.5 \Delta t u^{I}\sum _{J=1}^{N}\DP{u^{J}}{x^{J}}
             \right|_{\Dvect{i}}
         \Deqlab{eq: 41}
     \end{equation}
     
     と表される.

     \subsection{速度場が非定常な場合}

     速度場が非定常な場合, $\psi $ の $t$ に関する 2 階微分をこれまでの
     ように簡単に空間微分に置き換えることができない. たとえば 1 次元の
     場合,

     \begin{equation}
       \DP[2]{\psi }{t} 
       = -\DP{}{t}\left[  \left( \DP{}{x}(u\psi ) \right) \right]
       = -\DP{}{x}\left[ \psi \DP{u}{t}+u\DP{\psi }{t}\right]
       = -\DP{}{x}\left[ \psi \DP{u}{t}\right] + 
       \DP{}{x}\left[  u^{2} \DP{\psi }{x} \right], \Deqlab{eq: 42}
     \end{equation}
     
     となる. スキームには上式の右辺第 1 項が考慮されなければならない. 
     しかし $u$ と $\psi $ が結合していなければ, この項の存在は以下の
     手続きに従い簡単に取り込むことができる. まず 1 次元移流方程式
     \Deqref{eq: 2}を次のように時間方向に差分化する.

     \begin{equation}
       \frac{\psi ^{n+1}-\psi ^{n}}{\Delta t } = 
       -\DP{}{x}(u^{n+\frac{1}{2}}\psi ^{n}), \Deqlab{eq: 43}
     \end{equation}

     $\psi ^{n+1}, u^{n+\frac{1}{2}}$ をそれぞれ時刻 $t$ の周りで展開する.     
     このとき $\psi $ の $t$ に関する 2 階微分は\Deqref{eq: 42}のように表す.

     \begin{eqnarray}
      \psi _{i}^{n+1} &=& \psi _{i}^{n} + \left.\DP{\psi }{t}\right|^{n}
      \Delta t
      + \left.\frac{1}{2}\DP[2]{\psi }{r}\right|^{n}(\Delta t)^{2} 
      + O(\Delta t^{3}) \nonumber \\
      &=&\psi _{i}^{n} + \left.\DP{\psi }{t}\right|^{n}\Delta t
      - \frac{1}{2}(\Delta t )^{2}\DP{}{x}
       \left(\psi \left.\DP{u}{t}\right)\right|^{n}+ O(\Delta t^{3}), \\
      -\DP{}{x}(u^{n+\frac{1}{2}}\psi ^{n}) &=&
      -\DP{}{x}(u^{n}\psi ^{n}) 
      -\frac{1}{2}\Delta t \DP{}{x}\left(\psi
           \left.\DP{u}{t}\right)\right|^{n} + O(\Delta t^{2})
     \end{eqnarray}

     これらを\Deqref{eq: 43}に代入すると,

     \begin{equation}
       \left.\DP{\psi }{t}\right|^{n} =
       -\DP{}{x}(u^{n}\psi^{n})+\left.\DP{}{x}\left(
           \frac{1}{2}\Delta t (u^{n})^{2} 
         \DP{\psi ^{n}}{x}\right)\right|^{n}+ O(\Delta t^{2}), \Deqlab{eq: 45}
     \end{equation}

     となる. 右辺第 1 項を上流差分で空間方向に差分化すると\Deqref{eq: 8}
     が求まる.

     以上より速度場が非定常な場合には時刻 $t+\Delta t/2$ における速度場を用いて
     計算すればよいことがわかる. $u^{n+\frac{1}{2}}$ は,

     \[
       u^{n+\frac{1}{2}}=\frac{u^{n}+u^{n+1}}{2}+O(\Delta t^{2})
     \]
       
     等として求めてやればよい.

