% 表題 移流方程式の差分解法: FCT 法 Zalesak(1979)によるFCT % % 履歴 1998/01/19 小高正嗣 % \section{Zalesak(1979) による表現: 1次元} Zalesak(1979) はより一般化された FCT の表現を提案した. 現在 FCT と いうときには Zalesak(1979) による定義を指すのが普通である. この定 義に従うと全節で挙げたいくつかの問題点をいくらか解決することができ る. 以下ではこの Zalesak(1979) によって与えられた FCT について見て いくことにする. 簡単のため解く式は1次元の移流方程式\Deqref{eq: 1}とする. これを保 存形で差分化すると, \begin{equation} \rho _{i}^{n+1}=\rho _{i}^{n} - \frac{1}{\Delta x}[F_{i+\frac{1}{2}} -F_{i-\frac{1}{2}}], \Deqlab{eq: 13} \end{equation} となる. これに FCT を適用することを考える. \subsection{計算手順} Zalesak(1979) に従うとその計算手順は以下のようになる. \begin{enumerate} \item 精度の低い低次のスキームを用いてフラックスを計算する. これを $F_{i+\frac{1}{2}}^{L}$ と表す. このとき用いるスキームは上流差分 のような「単調な」(数値波を作らない)スキームであることが要請され る. \item 精度の高い高次のスキームを用いてフラックスを計算する. これを $F_{i+\frac{1}{2}}^{H}$ と表す. \item antidiffusive フラックスを以下のようにして定義する. \begin{equation} A_{i+\frac{1}{2}}\equiv F_{i+\frac{1}{2}}^{H}-F_{i+\frac{1}{2}}^{L}. \Deqlab{eq: 14} \end{equation} \item 低次のスキームで計算する. \begin{equation} \rho _{i}^{td}=\rho _{i}^{n} - \frac{1}{\Delta x}[F_{i+\frac{1}{2}}^{L} -F_{i-\frac{1}{2}}^{L}]. \Deqlab{eq: 15} \end{equation} \item 補正フラックスを適当な方法で計算する. 補正の与える際の条件は Boris and Book(1973) と同様である. \begin{equation} A_{i+\frac{1}{2}}^{c} = C_{i+\frac{1}{2}}A_{i+\frac{1}{2}}, \; 0 \le C_{i+\frac{1}{2}}\le 1. \Deqlab{eq: 16} \end{equation} なお補正係数 $C_{i+\frac{1}{2}}$ の与え方については後述する. \item 補正フラックスを与える. \begin{equation} \rho _{i}^{td}=\rho _{i}^{n} - \frac{1}{\Delta x}[A_{i+\frac{1}{2}}^{c} -A_{i-\frac{1}{2}}^{c}]. \Deqlab{eq: 17} \end{equation} \end{enumerate} このように定義することにより応用範囲が大きく広がる. 低次のスキーム には条件がつくが, 高次のスキームに対しては何も制約がないのであらゆ るスキームを用いることができる. 組み合わせるスキームを決めてしまえ ばあとは上記の手順に沿って機械的に計算を進めることができる. 実際計 算をする場合には $F_{i+\frac{1}{2}}^{L}, F_{i+\frac{1}{2}}^{H}$ を 計算する部分のコードだけを書き換えればよいことになるので, 種々のス キームを用いてそのパフォーマンスを比較検討することが容易になる. なお Harten and Zwas(1972) による自己調節混合スキーム (Self-Adjusting Hybrid Schemes; 以下 SAHS とする)は, 上記の 5 番 目に相当する所の計算方法が異なるだけで, その他の手順は同様である. \Deqref{eq: 14}$\sim $\Deqref{eq: 17}をまとめて書くこともできる. この場合 \begin{screen} \begin{equation} \rho _{i}^{n+1}=\rho _{i}^{n} - \frac{1}{\Delta x}\left\{[C_{i+\frac{1}{2}}F_{i+\frac{1}{2}}^{H}+ (1-C_{i+\frac{1}{2}})F_{i+\frac{1}{2}}^{L}] - [C_{i-\frac{1}{2}}F_{i-\frac{1}{2}}^{H}+ (1-C_{i-\frac{1}{2}})F_{i-\frac{1}{2}}^{L}]\right\}, \Deqlab{eq: 18} \end{equation} \end{screen} となる. FCT とは低次と高次のスキームの重ね合わせであり, その比率 $C_{i+\frac{1}{2}}$ の与え方に工夫を凝らしたものと捉えることができ る. $C_{i+\frac{1}{2}}$ の与え方を変えると SAHS にすることもできる. 一般に FCT として普通目にするのは\Deqref{eq: 18}式の表現である. ここまで来ると Boris and Book(1973), Book {\it et al.}(1975) の表 現とは一見異なるが, やっていることは同じである. \subsection{補正係数の決め方} 補正係数 $C_{i+\frac{1}{2}}$ は補正フラックスによって補正された値 が極値を取らないように, すなわち $\rho _{i}$ がある最大値 $\rho _ {i}^{max}$ 及び $\rho _{i}^{min}$ を超えないように与える. その際対 象となる格子点に関わる全ての補正フラックスの存在を考慮すれば, 全節 で挙げたような多次元に拡張した場合に発生する問題を回避することがで きるだろう. これは次のような手順で $C_{i+\frac{1}{2}}$ を与えるこ とで達成することができる. \begin{enumerate} \item antidifuusive フラックスに対しあらかじめ以下の条件を課しておく. \begin{eqnarray} A_{i+\frac{1}{2}} =0 & \mbox{if} & A_{i+\frac{1}{2}}(\rho _ {i+1}^{td} - \rho _{i}^{td}) < 0, \nonumber \\ & \mbox{and either} & A_{i+\frac{1}{2}}(\rho _ {i+2}^{td} - \rho _{i+1}^{td}) < 0, \nonumber \\ & \mbox{or} & A_{i+\frac{1}{2}}(\rho _ {i}^{td} - \rho _{i-1}^{td}) < 0. \Deqlab{eq: 19} \end{eqnarray} これは Book {\it et al.}(1975) で言うところの ``terrace'' を作 らないようにするためのものである(Book {\it et al.} 1975, 図9参 照). 実際は $ A_{i+\frac{1}{2}}(\rho _ {i+1}^{td} - \rho _ {i}^{td}) < 0 $ となることは稀なので, この条件を与えなくてもさ ほど問題にはならない. \item $\rho _{i}^{max}, \rho _{i}^{min}$ を以下のようにして与える. \begin{eqnarray} \rho _{i}^{max} &=& \mbox{max}(\rho _{i-1}^{td}, \; \rho _ {i}^{td}, \; \rho _{i+1}^{td}), \Deqlab{eq: 20} \\ \rho _{i}^{min} &=& \mbox{min}(\rho _{i-1}^{td}, \; \rho _ {i}^{td}, \; \rho _{i+1}^{td}). \Deqlab{eq: 21} \end{eqnarray} または clipping を抑えるために, \begin{eqnarray} \rho _{i}^{a} &=& \mbox{max}(\rho _{i}^{n}, \; \rho _{i}^{td}), \Deqlab{eq: 22} \\ \rho _{i}^{max} &=& \mbox{max}(\rho _{i-1}^{a}, \; \rho _ {i}^{a}, \; \rho _{i+1}^{a}), \Deqlab{eq: 23} \\ \rho _{i}^{b} &=& \mbox{min}(\rho _{i}^{n}, \; \rho _{i}^{td}), \Deqlab{eq: 24} \\ \rho _{i}^{min} &=& \mbox{min}(\rho _{i-1}^{b}, \; \rho _ {i}^{b}, \; \rho _{i+1}^{b}), \Deqlab{eq: 25} \end{eqnarray} とする. \item 以下の3つの量を用意する. \begin{eqnarray} P_{i}^{+} &=& \mbox{格子点 $i$ に入る全ての antidiffusive フ ラックスの和} \nonumber \\ &=& \mbox{max}(0, A_{i-\frac{1}{2}}) - \mbox{min}(0, A_{i+\frac{1}{2}}), \Deqlab{eq: 26} \\ Q_{i}^{+} &=& (\rho _{i}^{max}-\rho _{i})\Delta x, \Deqlab{eq: 27} \\ R_{i}^{+} &=& \left\{ \begin{array}{lcl} \mbox{min}(1, Q_{i}^{+}/P_{i}^{+}) & \mbox{if} & P_{i}^{+}>0, \\ 0 & \mbox{if} & P_{i}^{+}=0. \end{array} \right. \Deqlab{eq: 28} \end{eqnarray} $R_{i}^{+}$ は $i$ 点に流入しうる antidifuusive フラックスの上 限を決める. \item 同様に以下の3つの量を用意する. \begin{eqnarray} P_{i}^{-} &=& \mbox{格子点 $i$ から出て行く全ての antidiffusive フ ラックスの和} \nonumber \\ &=& \mbox{max}(0, A_{i+\frac{1}{2}}) - \mbox{min}(0, A_{i-\frac{1}{2}}), \Deqlab{eq: 29} \\ Q_{i}^{-} &=& (\rho _{i}-\rho _{i}^{min})\Delta x, \Deqlab{eq: 30} \\ R_{i}^{-} &=& \left\{ \begin{array}{lcl} \mbox{min}(1, Q_{i}^{-}/P_{i}^{-}) & \mbox{if} & P_{i}^{-}>0, \\ 0 & \mbox{if} & P_{i}^{-}=0. \end{array} \right. \Deqlab{eq: 31} \end{eqnarray} $R_{i}^{-}$ は $i$ 点から流出しうる antidifuusive フラックスの上 限を決める. \item 上で与えた $R_{i}^{+}, R_{i}^{-}$ を用いて補正係数を以下のよ うに与える. \begin{equation} C_{i+\frac{1}{2}} = \left\{ \begin{array}{lcl} \mbox{min}(R_{i+1}^{+}, R_{i}^{-}) & \mbox{if} & A_{i+\frac{1}{2}}\ge 0, \\ \mbox{min}(R_{i}^{+}, R_{i+1}^{-}) & \mbox{if} & A_{i+\frac{1}{2}} < 0. \end{array} \right. \Deqlab{eq: 32} \end{equation} \end{enumerate}