% 表題  移流方程式の差分解法: FCT 法  Zalesak(1979)によるFCT(2D)
%
% 履歴  1998/01/20 小高正嗣
%       1998/05/11 小高正嗣
%       1998/05/19 小高正嗣
%       2006/02/13 小高正嗣: ps ファイルへの相対パスを修正
%
\section{Zalesak(1979) による表現: 2次元}

    \subsection{計算手順と補正係数の決め方}

    2次元の FCT は1次元のものをそのまま拡張することで得ることができる.
    \Deqref{eq: 18}を2次元形に拡張すると,

    \begin{screen}
      \begin{eqnarray}
        \rho _{i.j}^{n+1}&=&\rho _{i.j}^{n} - \frac{1}{\Delta x\Delta
        y}\{[C_{i+\frac{1}{2},j}F_{i+\frac{1}{2},j}^{H} +
          (1-C_{i+\frac{1}{2},j})F_{i+\frac{1}{2},j}^{L}] \nonumber \\ 
          && - [C_{i-\frac{1}{2},j}F_{i-\frac{1}{2},j}^{H}+
          (1-C_{i-\frac{1}{2},j})F_{i-\frac{1}{2},j}^{L}], \nonumber
          \\ && + [C_{i,j+\frac{1}{2}}G_{i,j+\frac{1}{2}}^{H}+
          (1-C_{i,j+\frac{1}{2}})G_{i,j+\frac{1}{2}}^{L}] \nonumber \\ 
          && - [C_{i,j-\frac{1}{2}}G_{i,j-\frac{1}{2}}^{H}+
          (1-C_{i,j-\frac{1}{2}})G_{i,j-\frac{1}{2}}^{L}]\},
        \Deqlab{eq: 33}
    \end{eqnarray}
    \end{screen}
    
    となる. 補正係数 $C_{i\pm \frac{1}{2},j}, C_{i,j\pm \frac{1}{2}}$
    を決める手順は\Deqref{eq: 19}$\sim $\Deqref{eq: 32}と同様である.

    {
    \setcounter{equation}{18}
    \renewcommand{\theequation}{\arabic{equation}'}
    \begin{enumerate}
      \item antidifuusive フラックスに対しあらかじめ以下の条件を課しておく.
        \begin{eqnarray}
          A_{i+\frac{1}{2},j} =0 & \mbox{if} & A_{i+\frac{1}{2},j}(\rho _
          {i+1,j}^{td} - \rho _{i,j}^{td}) < 0, \nonumber \\
          & \mbox{and either} & A_{i+\frac{1}{2},j}(\rho _
          {i+2,j}^{td} - \rho _{i+1,j}^{td}) < 0, \nonumber \\
          & \mbox{or} & A_{i+\frac{1}{2},j}(\rho _
          {i,j}^{td} - \rho _{i-1,j}^{td}) < 0, \nonumber \\
          A_{i,j+\frac{1}{2}} =0 & \mbox{if} & A_{i,j+\frac{1}{2}}(\rho _
          {i,j+1}^{td} - \rho _{i,j}^{td}) < 0, \nonumber \\
          & \mbox{and either} & A_{i,j+\frac{1}{2}}(\rho _
          {i,j+2}^{td} - \rho _{i,j+1}^{td}) < 0, \nonumber \\
          & \mbox{or} & A_{i,j+\frac{1}{2}}(\rho _
          {i,j}^{td} - \rho _{i,j-1}^{td}) < 0. \Deqlab{eq: 19'}
        \end{eqnarray}

      \item $\rho _{i}^{max}, \rho _{i}^{min}$ を以下のようにして与える.
        \begin{eqnarray}
          \rho _{i,j}^{max} &=& \mbox{max}(\rho _{i-1,j}^{td}, \; \rho _
          {i,j}^{td}, \; \rho _{i+1,j}^{td}, \; \rho _{i,j-1}^{td}, \; 
          \rho _{i,j+1}^{td}), \Deqlab{eq: 20'} \\
          \rho _{i,j}^{min} &=& \mbox{min}(\rho _{i-1,j}^{td}, \; \rho _
          {i,j}^{td}, \; \rho _{i+1,j}^{td}, \; \rho _{i,j-1}^{td}, \; 
          \rho _{i,j+1}^{td}). \Deqlab{eq: 21'}
        \end{eqnarray}

        または clipping を抑えるために, 
        \begin{eqnarray}
          \rho _{i,j}^{a} &=& \mbox{max}(\rho _{i,j}^{n}, \; \rho _
          {i,j}^{td}), \Deqlab{eq: 22'} \\ 
          \rho _{i}^{max} &=& \mbox{max}(\rho _{i-1,j}^{a}, \; \rho _ {i,j}^{a}
          , \; \rho _{i+1,j}^{a}, \; \rho _{i,j-1}^{a}, \; 
          \rho _{i,j+1}^{a}), \Deqlab{eq: 23'} \\ 
          \rho _{i,j}^{b} &=&  \mbox{min}(\rho _{i,j}^{n}, \; \rho _
          {i,j}^{td}), \Deqlab{eq: 24'} \\ 
          \rho _{i,j}^{min} &=& \mbox{min}(\rho _{i-1,j}^{b}, \;
          \rho _ {i,j}^{b}, \; \rho _{i+1,j}^{b}, \; \rho _{i,j-1}^{b}, \;
          \rho _{i,j-1}^{b}), \Deqlab{eq: 25'}
        \end{eqnarray}

        とする.

      \item 以下の3つの量を用意する.
        \begin{eqnarray}
          P_{i,j}^{+} &=& \mbox{格子点 $i,j$ に入る全ての antidiffusive フ
          ラックスの和} \nonumber \\ 
          &=& \mbox{max}(0, A_{i-\frac{1}{2},j}) - \mbox{min}(0,
          A_{i+\frac{1}{2},j}) \nonumber \\
          && + \mbox{max}(0, A_{i,j-\frac{1}{2}}) - 
          \mbox{min}(0, A_{i,j+\frac{1}{2}}), \Deqlab{eq: 26'} \\
          Q_{i,j}^{+} &=& (\rho _{i,j}^{max}-\rho _{i,j})\Delta x \Delta y, 
          \Deqlab{eq: 27'} \\
          R_{i,j}^{+} &=& \left\{
          \begin{array}{lcl}
            \mbox{min}(1, Q_{i.j}^{+}/P_{i,j}^{+}) & \mbox{if} & P_{i,j}^{+}>0,
            \\
            0 & \mbox{if} & P_{i,j}^{+}=0.
          \end{array}
          \right. \Deqlab{eq: 28'} 
        \end{eqnarray}

      \item 同様に以下の3つの量を用意する.
        \begin{eqnarray}
          P_{i,j}^{-} &=& \mbox{格子点 $i,j$ から出て行く全ての antidiffusive
          フラックスの和} \nonumber \\ 
          &=& \mbox{max}(0, A_{i+\frac{1}{2},j}) - \mbox{min}(0,
          A_{i-\frac{1}{2},j}) \nonumber \\
          && + \mbox{max}(0, A_{i,j+\frac{1}{2}}) - 
          \mbox{min}(0,A_{i,j-\frac{1}{2}}), \Deqlab{eq: 29'} \\
          Q_{i,j}^{-} &=& (\rho _{i,j}-\rho _{i,j}^{min})\Delta x \Delta y
          , \Deqlab{eq: 30'} \\
          R_{i,j}^{-} &=& \left\{
          \begin{array}{lcl}
            \mbox{min}(1, Q_{i,j}^{-}/P_{i,j}^{-}) & \mbox{if} & P_{i,j}^{-}>0,
            \\
            0 & \mbox{if} & P_{i,j}^{-}=0. 
          \end{array}
          \right. \Deqlab{eq: 31'} 
        \end{eqnarray}

      \item 補正係数を以下のように与える.
        \begin{eqnarray}
          C_{i+\frac{1}{2},j} &=& \left\{
            \begin{array}{lcl}
               \mbox{min}(R_{i+1,j}^{+}, R_{i,j}^{-}) & \mbox{if} &
              A_{i+\frac{1}{2},j}\ge 0, \\
               \mbox{min}(R_{i,j}^{+}, R_{i+1,j}^{-}) & \mbox{if} &
              A_{i+\frac{1}{2},j} < 0.
            \end{array}
            \right. \nonumber \\
            && \Deqlab{eq: 32'} \\
          C_{i,j+\frac{1}{2}} &=& \left\{
            \begin{array}{lcl}
               \mbox{min}(R_{i,j+1}^{+}, R_{i,j}^{-}) & \mbox{if} &
              A_{i,j+\frac{1}{2}}\ge 0, \\
               \mbox{min}(R_{i,j}^{+}, R_{i,j+1}^{-}) & \mbox{if} &
              A_{i,j+\frac{1}{2}} < 0.
            \end{array}
            \right. \nonumber 
       \end{eqnarray}
    \end{enumerate}
    }
    \addtocounter{equation}{1}
    \newpage
    \subsection{計算例}

    実際に式\Deqref{eq: 33}を\Deqref{eq: 19'}$\sim $\Deqref{eq: 32'}を
    用いて計算する. 計算設定は Zalesak(1979) と同様のもの, 初期値の与
    え方は Smolarkiewicz(1983) で用いられたもの同様とする.

    \begin{itemize}
      \item 計算領域に $0\leq x \leq 1, 0\leq y \leq 1$ をとる.
        それぞれの方向に 100 個の格子に分割する. すなわち $\Delta x = 
        \Delta y = 0.01$ である.
      
      \item 速度場は $(x_{0}, y_{0})=(0.5, 0.5)$ を中心に角速度 $\omega = 0.1$
        で剛体回転する場を与える. すなわち $(u, v) = ( - \omega (y-y_{0}), 
        \omega (x-x_{0}))$ である.

      \item $\rho $ の初期値は円錐形分布を与える. 底面の円の中心は 
        $(x_{m}, y_{m})=(0.75, 0.5)$ ,半径は 0.15 , 円錐の高さは 4 とする.

      \item 境界条件は全ての変数が壁で対称となるように与える.

      \item 時間格子間隔 $\Delta t$ は 0.1 とした. これより与えられる
       クーラン数 MAX$(u\Delta t/\Delta x)$ は約 0.7 である.

    \end{itemize}
    以上の設定で計算すると, 時間方向に 628 ステップ計算でほぼ一回転する.
    計算領域と初期値の分布を\Dfigref{fig: 7}に示す.

    \begin{figure}[bth]
      \begin{center}
        \Depsf[][80mm]{ps-fig/init-cone.ps}
        \caption{計算領域と初期値の分布.} \Dfiglab{fig: 7}
      \end{center}
    \end{figure}

    比較のためまず上流差分スキームを用いて計算を行なった. スキームを
    書くと以下のようになる.

    \begin{eqnarray}
      \rho _{i,j}^{n+1} &=& \rho _{i,j}^{n} - 
      \{F_{x}(\rho _{i,j}^{n},\rho _{i+1,j}^{n},u_{i+\frac{1}{2},j}^{n}) - 
      F_{x}(\rho _{i-1,j}^{n},\rho _{i,j}^{n},u_{i-\frac{1}{2},j}^{n})\}
      \nonumber \\
      && - \{F_{y}(\rho _{i,j}^{n},\rho _{i,j+1}^{n},v_{i,j+\frac{1}{2}}^{n})
      - F_{y}(\rho _{i,j-1}^{n},\rho _{i,j}^{n},v_{i,j-\frac{1}{2}}^{n})\}. 
      \Deqlab{eq: 34}
    \end{eqnarray}
    ここで,
    \begin{eqnarray*}
      F_{x}(\rho _{i,j}^{n},\rho _{i+1,j}^{n},u_{i+\frac{1}{2},j}^{n})&=& 
      [(u_{i+\frac{1}{2},j}^{n}+|u_{i+\frac{1}{2},j}^{n}|)\rho _{i,j}^{n} + 
      (u_{i+\frac{1}{2},j}^{n}-|u_{i+\frac{1}{2},j}^{n}|)\rho _{i+1,j}^{n}]
      \frac{\Delta t}{2\Delta x}, \\
      F_{y}(\rho _{i.j}^{n},\rho _{i,j+1}^{n},v_{i,j+\frac{1}{2}}^{n})&=&
      [(v_{i,j+\frac{1}{2}}^{n}+|v_{i,j+\frac{1}{2}}^{n}|)\rho _{i,j}^{n} + 
      (v_{i,j+\frac{1}{2}}^{n}-|v_{i,j+\frac{1}{2}}^{n}|)\rho _{i,j+1}^{n}]
      \frac{\Delta t}{2\Delta y}, 
    \end{eqnarray*}
    である. これを時間方向に修正オイラー法(2 次のルンゲクッタ)を用いて
    計算した. 

    1回転後(628ステップ)と3回転後(1884ステップ)後の結果を\Dfigref{fig:
    8}と\Dfigref{fig: 9}にそれぞれ示す. 上流差分では数値拡散が大きいた
    め1回転後で既に初期分布は大きく損なわれている.  3回転もすると初期
    分布はあとかたもなくなってしまう.
     
    \begin{figure}[p]
      \begin{center}
        \Depsf[][80mm]{ps-fig/upstream1.ps}
        \caption{上流差分による計算. 1回転(628ステップ)後の結果.}
        \Dfiglab{fig: 8}
        \Depsf[][80mm]{ps-fig/upstream2.ps}
        \caption{上流差分による計算. 3回転(1884ステップ)後の結果.}
        \Dfiglab{fig: 9}
      \end{center}
    \end{figure}
    
    同様の計算を FCT を用いて行なった. 低次のスキームに上流差分, 高次
    のスキームに2次中心差分を用いている. 1回転後(628ステップ)と3回転後
    (1884ステップ)後の結果を\Dfigref{fig: 10}と\Dfigref{fig: 11}にそれ
    ぞれ示す. 全体的に数値拡散は抑えられているが, 山の頂上付近が削られ
    てしまい clipping が回避できていないことを示している.

    \begin{figure}[p]
      \begin{center}
        \Depsf[][80mm]{ps-fig/fct-exp1.ps}
        \caption{FCT による計算. 1回転(628ステップ)後の結果.}
        \Dfiglab{fig: 10}
        \Depsf[][80mm]{ps-fig/fct-exp2.ps}
        \caption{FCT による計算. 3回転(1884ステップ)後の結果.}
        \Dfiglab{fig: 11}
      \end{center}
    \end{figure}

    続いて初期値分布を変えて同様の計算を行なった. \Dfigref{fig: 7}で円
    錐を置いた所に同じ高さを持つ円筒を置く(\Dfigref{fig: 12}参照). こ
    こでも低次のスキームに上流差分,高次のスキームに2次中心差分を用いた.
    1回転後(628ステップ)の結果を\Dfigref{fig: 13}に示す. 初期値に円錐
    分布を与えた場合と比べ clipping はあまり目立たない. 

    \begin{figure}[p]
      \begin{center}
        \Depsf[][80mm]{ps-fig/init-cylind.ps}
        \caption{計算領域と初期値の分布. 円筒分布を与えた場合.}
        \Dfiglab{fig: 12}
        \Depsf[][80mm]{ps-fig/fct-exp3.ps}
        \caption{FCT による計算. 1回転(628ステップ)後の結果. 初期値に円筒分布
        を与えた場合} \Dfiglab{fig: 13}
      \end{center}
    \end{figure}


     



