% 表題 移流方程式の差分解法: 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}