% 表題 移流方程式の差分解法: FCT 法 1D-SHASTA % % 履歴 1998/01/19 小高正嗣 % 1998/05/11 小高正嗣 % 1998/05/19 小高正嗣 % 2006/02/13 小高正嗣: ps ファイルへの相対パスを修正 % \section{1次元 FCT: SHASTA} ここでは Boris and Book(1973) によって提案されたオリジナルの FCT, 及び そこで示されている1次元移流方程式 \begin{equation} \DP{\rho }{t} + \DP{}{x}(\rho u) = 0, \Deqlab{eq: 1} \end{equation} の差分スキーム SHASTA(Sharp And Smooth Transport Algolithm) につい て解説する. FCTは後述する Zalesak(1979) による定義をもって与えられ るのが一般的であるが, FCT の基本的な考え方を知るためにまずはオリジ ナルの FCT について見ていくことにする. Boris and Book(1973) は FCT には一般的に 2 つの段階(stage)があると している. その第一段階を輸送過程(transport stage), 第二段階を補正 過程(corrective stage)と呼んでいる. 輸送過程で生じた数値拡散を補正 過程で補うというのが FCT の基本的な考え方である. \subsection{輸送過程} \begin{figure}[p] \begin{center} \Depsf[80mm][]{ps-fig/fig1.ps} \caption{1D SHASTA の輸送過程(Boris and Book(1973),図1を基に作成)} \Dfiglab{fig: 1} \end{center} \end{figure} SHASTA の輸送過程は以下のように与えられる. (\Dfigref{fig: 1}参照). \begin{enumerate} \item 隣接する格子点間で台形を作り(\Dfigref{fig: 1}(a)), その面積が 保存するように質量を輸送する. 時刻 $t_{0} (\equiv n\Delta t)$ のときそれぞれ格子点上で密度は $\rho _ {i}^{n}, \rho _{i+1}^{n}$, 速度は $u_{i}^{n}, u_{i+1}^{n}$ とする. 格子間隔は $\Delta x$, 時間間隔は $\Delta t$ とする. 台形の面積を一定に保つので, その高さは底辺の長さに半 比例して与えられる. よって時刻 $t_{0}+\Delta t$ の台形の高さ (\Dfigref{fig: 1}(b))はそれぞれ, \begin{eqnarray} \rho _{1} &=& \frac{\rho _{i+1}^{n}\Delta x}{\Delta x + \Delta t (u_{i+1}^{n} - u_{i}^{n})}, \nonumber \\ && \Deqlab{eq: 2} \\ \rho _{2} &=& \frac{\rho _{i}^{n}\Delta x}{\Delta x + \Delta t (u_{i+1}^{n} - u_{i}^{n})}, \nonumber \end{eqnarray} と与えられる. $\rho _{1}, \rho _{2}$ は $\rho _{i}^{n}, \rho _ {i+1}^{n}$ が正のとき及び \begin{equation} \left|u \frac{\Delta t}{\Delta x}\right| < \frac{1}{2}, \Deqlab{eq: 3} \end{equation} の場合に常に正の値をとる. 以下ではこの条件が満たされるものとして 手続きを進める. \item 輸送された台形を格子を含むセル上に分割し, それぞれのセルに含 まれる面積を格子点上に割り付ける. 具体的には\Dfigref{fig: 1}(c) の A は $i$に, B は $i+1$ に割り当てる. 割り当てた面積を格子間隔 で割った値が次ステップの値になる. セルの境界上の値は線形補間で求 める. \Dfigref{fig: 1}(c)A の面積は \[ \frac{1}{2}\left[\rho _{1} + \left(\rho _{2}+\frac{\rho _{1}-\rho _{2}} {\Delta x + \Delta t(u_{i+1}^{n} - u_{i}^{n})}\cdot (\frac{1}{2}\Delta x + u_{i+1}^{n}\Delta t)\right)\right]\cdot (\frac{1}{2}\Delta x + u_{i+1}^{n}\Delta t), \] B の面積は \[ \frac{1}{2}\left[\rho _{2} + \left(\rho _{2}+\frac{\rho _{1}-\rho _{2}} {\Delta x + \Delta t(u_{i+1}^{n} - u_{i}^{n})}\cdot (\frac{1}{2}\Delta x - u_{i}^{n}\Delta t)\right)\right]\cdot (\frac{1}{2}\Delta x - u_{i}^{n}\Delta t), \] と与えられる. $i$ 格子に割り付けられる面積は\Dfigref{fig: 1}(c)B と $i-1$ からの寄与 \[ \frac{1}{2}\left[\rho _{2} + \left(\rho _{3}+\frac{\rho _{2}-\rho _{3}} {\Delta x + \Delta t(u_{i}^{n} - u_{i-1}^{n})}\cdot (\frac{1}{2}\Delta x + u_{i}^{n}\Delta t)\right)\right]\cdot (\frac{1}{2}\Delta x + u_{i}^{n}\Delta t), \] である. ただし, \[ \rho _{3} = \frac{\rho _{i-1}^{n}\Delta x}{\Delta x + \Delta t (u_{i}^{n} - u_{i-1}^{n})}, \] である. 割り付けられた面積を格子間隔で割って時刻 $t_{0}+\Delta t$ の$i$ 格子点での密度 $\rho _{i}^{n+1}$ を得る. \end{enumerate} 以上をまとめると, \begin{screen} \begin{eqnarray} \rho _{i}^{n+1} &=& \frac{1}{2}Q_{-}^{2}(\rho _{i-1}^{n}-\rho _{i}^{n}) + \frac{1}{2}Q_{+}^{2}(\rho _{i+1}^{n}-\rho _{i}^{n}) + (Q_{+}+Q_{-})\rho _{i}^{n},\nonumber \\ Q_{\pm } &=& \left(\frac{1}{2} \mp u_{i}^{n}\frac{\Delta t}{\Delta x} \right)\left/\left[1 \pm (u_{i\pm 1}^{n} - u_{i}^{n}) \frac{\Delta t}{\Delta x}\right.\right], \Deqlab{eq: 4} \end{eqnarray} \end{screen} となる. 速度一定の場合は $\epsilon \equiv u_{i}^{n}\Delta t/\Delta x$ として, \begin{screen} \begin{equation} \rho _{i}^{n+1} = \rho _{i}^{n} - \frac{\epsilon }{2} (\rho _{i+1}^{n}-\rho _{i-1}^{n}) + \left(\frac{1}{8} + \frac{\epsilon ^{2}}{2}\right) (\rho _{i+1}^{n}-2\rho _{i}^{n}+\rho _{i-1}^{n}), \Deqlab{eq: 5} \end{equation} \end{screen} となる. これは中心差分で表現した移流に拡散が加わった形になっていて, 速度に依存しない拡散を除けば 2 段の Lax-Wendroff スキームと同じである \footnote{2 段の Lax-Wendroff スキームとは以下のような方法である. \begin{eqnarray*} \rho _{i+\frac{1}{2}}^{n+\frac{1}{2}} &=& \frac{1}{2}(\rho _{i}^{n}+ \rho _{i+1}^{n}) - \frac{\Delta t}{2\Delta x}(f_{i+1}^{n}-f_{i}^{n}) \\ \rho _{i}^{n+1} &=& \rho _{i}^{n} - \frac{\Delta t}{\Delta x} (f_{i+\frac{1}{2}}^{n+\frac{1}{2}}-f_{i-\frac{1}{2}}^{n+\frac{1}{2}}) \end{eqnarray*} ただし $f_{i}^{n}=u_{i}^{n}\rho _{i}^{n}$ 等である. まず $\rho _ {i+\frac{1}{2}}^{n+\frac{1}{2}}$ を求め, その値を用いて$\rho _ {i+1}^{n+1}$ を求める. $u_{i}^{n}$ が一定の場合は $\rho _ {i+\frac{1}{2}}^{n+\frac{1}{2}}$ を消去して, \[ \rho _{i}^{n+1} = \rho _{i}^{n} - \frac{1}{2}\epsilon (\rho _{i+1}^{n} -\rho _{i-1}^{n}) + \frac{1}{2}\epsilon ^{2}(\rho _{i+1}^{n}-2\rho _{i} ^{n}+\rho _{i-1}^{n}) \] となる.}. \subsection{補正過程} \begin{figure}[b] \begin{center} \Depsf[80mm][]{ps-fig/fig2.ps} \caption{1D の不連続面の伝播と antidiffusion の影響 (Boris and Book(1973),図2を基に作成)} \Dfiglab{fig: 2} \end{center} \end{figure} \Deqref{eq: 5}式から明らかなように SHASTA には速度に依存した拡散が 含まれている. $u=0$ の場合にも \begin{equation} \rho _{i}^{n+1} = \rho _{i}^{n} + \frac{1}{8}(\rho _{i+1}^{n}-2\rho _{i}^{n} +\rho _{i-1}^{n}), \Deqlab{eq: 6} \end{equation} となって拡散が残ってしまう. これを除いた修正解 $\overline{\rho }_ {i}^{n+1}$ を求めるために, \Deqref{eq: 6}式を陽解法で逆に解く. \begin{equation} \overline{\rho }_{i}^{n+1} = \rho _{i}^{n+1} - \frac{1}{8}(\rho _{i+1}^{n+1}-2\rho _{i}^{n+1} +\rho _{i-1}^{n+1}), \Deqlab{eq: 7} \end{equation} \Deqref{eq: 7}式右辺第二項を ``antidiffusion'' と呼ぶ. \Deqref{eq: 7}は形式的に移流形にすることができる. すなわち, \begin{eqnarray} \overline{\rho }_{i}^{n+1} &=& \rho _{i}^{n+1} - (f_{i+\frac{1}{2}} - f_{i-\frac{1}{2}}), \Deqlab{eq: 8} \\ f_{i\pm \frac{1}{2}} &\equiv & \pm \frac{1}{8}(\rho _{i\pm 1}^{n+1}-\rho _{i}^{n+1}), \Deqlab{eq: 9} \end{eqnarray} となる. $f_{i\pm \frac{1}{2}}$ は antidiffusive フラックス と呼ぶ. この antidiffusion または antidiffusive フラックスを加えれば数値拡散 を除去することができるので,いつでも正しい解が得られるかというと実 はそうではない. 例えば\Dfigref{fig: 2}のような不連続面が一定の速度 で移流されるような場合を考えてみよう. ある時刻で図のように正しい解 が得られているとする. これに antidiffusion を加えると図に示した矢 印のように $i$ での値は増加し $i+1$ の値は減少してしまう. よって新 たに極値を作り出すことになる. とくに非負の物理量を計算している場 合は $i+1$ に負の値が生じてしまう. 明らかにこれは非物理的な解であ る. antidiffusion を加えるときは以上のような非物理的な解を生じないよう に注意しなければならない. すなわち, \begin{itemize} \item 新たに極大, 極小を作らない. \item すでに存在する極大(極小)をさらに大きく(小さく)はしない. \end{itemize} という条件を満たすように $f_{i\pm \frac{1}{2}}$ を加えなければなら ない. そこで以下のような補正フラックス( corrected flux ) $f^{c}_{i\pm \frac{1}{2}}$ を新たに定義し, これを\Deqref{eq: 8}式の $f_{i\pm \frac{1}{2}}$ と置き換えて用いることにする. \begin{screen} \begin{equation} f^{c}_{i+\frac{1}{2}} \equiv S_{i+\frac{1}{2}} \mbox{max} \{0, \mbox{min}(S_{i+\frac{1}{2}}\Delta _{i-\frac{1}{2}}, \; |f_{i+\frac{1}{2}}|, \; S_{i+\frac{1}{2}}\Delta _{i+\frac{3}{2}})\}, \Deqlab{eq: 10} \end{equation} \begin{equation} \Delta _{i+\frac{1}{2}} = \rho _{i+1}^{n+1}-\rho _{i}^{n+1}, \Deqlab{eq: 11} \end{equation} \begin{equation} S = \left\{ \begin{array}{lcl} +1 & \mbox{if} & f_{i+\frac{1}{2}} \geq 0 \\ -1 & \mbox{if} & f_{i+\frac{1}{2}} < 0 \end{array} \right. \Deqlab{eq: 12} \end{equation} \end{screen} \begin{figure}[b] \begin{center} \Depsf[120mm][]{ps-fig/fig3.ps} \caption{変数の取り得る空間配置と補正流束の作用範囲 (Book {\it et al.}(1973),図6を基に作成)} \Dfiglab{fig: 3} \end{center} \end{figure} この補正フラックスはどのような場合に作用するよう定義されたのか具体 的な例を挙げて調べてみる. \Dfigref{fig: 3}には $f_{i+\frac{1}{2}}>0$ のときに取り得る配置を示してある. (a) のよう な場合 $f^{c}_{i\pm \frac{1}{2}}$ は図に示した矢印の範囲まで補正し うる大きさをとることができる. それ以上の大きさになると新たに極値が 生じてしまう. これに対し (b)$\sim $(d) の場合$f^{c}_{i\pm \frac{1}{2}}$ は 0 にしなくてはならない. そうしないと既に存在して いる極大(極小)がさらに成長してしまい, 先に述べたような非物理的な解 が生じるようになってしまうからである. \Deqref{eq: 10}$\sim $\Deqref{eq: 12}式で与えた補正フラックスは\Dfigref{fig: 3}(a)のよ うな場合に極値を生じないよう作用するように巧妙に定義したものなので ある. \newpage \subsection{まとめ: 1D-SHASTA} 1D-SHASTA の実行手続きを以下にまとめる \footnote{ ここで示した補正過程は\Deqref{eq: 5}式に基づいたものなので, \Deqref{eq: 4}式と\Deqref{eq: 5}式のずれに伴う部分は補正されない. $O((u_{i\pm 1}-u_{i})\Delta t/\Delta x) \sim O(\epsilon ) < 1$ として \Deqref{eq: 4}式を展開すると, \Deqref{eq: 5}式とのずれは, \begin{eqnarray*} && \left(\frac{1}{2}+\epsilon \right) ^{2}(u_{i-1}^{n}-u_{i}^{n})(\rho _ {i-1}^{n} - \rho _{i}^{n})\frac{ \Delta t}{\Delta x}- \left(\frac{1}{2}-\epsilon \right) ^{2}(u_{i+1}^{n}-u_{i}^{n})(\rho _ {i+1}^{n} - \rho _{i}^{n})\frac {\delta t}{\Delta x}\\ &+& \left[\left(\frac{1}{2}+\epsilon \right) (u_{i-1}^{n}-u_{i}^{n}) - \left(\frac{1}{2}-\epsilon \right)(u_{i+1}^{n}-u_{i}^{n}) \right]\frac{\Delta t}{\Delta x}\rho _{i}^{n}, \end{eqnarray*} となる. これは $O(\epsilon )$ 以上のオーダーなので補正にとって最も 重要な項は 0 次の拡散項であることがわかる. なおこれから $O(\epsilon )$ の項のみを抜き出すと \[ \frac{1}{4}(u_{i-1}^{n}-u_{i}^{n})(\rho _{i-1}^{n} + \rho _ {i}^{n})\frac{\Delta t}{\Delta x} - \frac{1}{4}(u_{i+1}^{n}-u_{i}^{n}) (\rho _{i+1}^{n} + \rho _ {i}^{n})\frac{\Delta t}{\Delta x} \] となり, 拡散ではない. 速度場が滑らかで $O((u_{i\pm 1}-u_{i})\Delta t/\Delta x) < O(\epsilon )$ とみなせるような場合には, このずれを無視することができる. }. \begin{screen} \begin{enumerate} \item 輸送過程: \Deqref{eq: 4}式を解く. \begin{eqnarray*} \rho _{i}^{n+1} &=& \frac{1}{2}Q_{-}^{2}(\rho _ {i-1}^{n}-\rho _{i}^{n}) + \frac{1}{2}Q_{+}^{2}(\rho _ {i+1}^{n}-\rho _{i}^{n}) + (Q_{+}+Q_{-})\rho _ {i}^{n},\nonumber \\ Q_{\pm } &=& \left(\frac{1}{2} \mp u_{i}^{n}\frac{\Delta t}{\Delta x} \right)\left/\left[1 \pm (u_{i\pm 1}^{n} - u_{i}^{n}) \frac{\Delta t}{\Delta x}\right.\right]. \end{eqnarray*} \item 補正過程(1): 補正フラックスを計算する\footnotemark \[ f^{c}_{i+\frac{1}{2}} = S_{i+\frac{1}{2}} \mbox{max} \{0, \mbox{min}(S_{i+\frac{1}{2}}\Delta _{i-\frac{1}{2}}, \; |f_{i+\frac{1}{2}}|, \; S_{i+\frac{1}{2}}\Delta _ {i+\frac{3}{2}})\}, \] \[ f_{i+\frac{1}{2}}= \frac{1}{8}(\rho _{i+1}^{n+1}-\rho _{i}^{n+1}),\; \Delta _{i+\frac{1}{2}} = \rho _{i+1}^{n+1}-\rho _ {i}^{n+1}, \] \[ S = \left\{ \begin{array}{lcl} +1 & \mbox{if} & f_{i+\frac{1}{2}} \geq 0, \\ -1 & \mbox{if} & f_{i+\frac{1}{2}} < 0. \end{array} \right. \] \item 補正過程(2): 補フラックスによる移流を計算する. \[ \overline{\rho }_{i}^{n+1} = \rho _{i}^{n+1} - (f^{c}_{i+\frac{1}{2}} - f^{c}_{i-\frac{1}{2}}). \] \end{enumerate} \end{screen} \footnotetext{ 実際に計算を行なう場合 Boris and Book(1973) では antidiffusion 係 数(仮に $\mu $ とおく)を以下のようにして精度を上げている. \[ \mu = \left\{ \begin{array}{ll} 0.125[1+(1-0.25\overline{\rho }/|D_{i+\frac{1}{2}}|)^{2}], & |D_{i+\frac{1}{2}}|>0.25 \overline{\rho } \\ 0.125 & \mbox{otherwise}. \end{array} \right. \] ただし $\overline{\rho }=\sqrt{\rho _{min}\rho_{max}}$ であるが, $D_{i+\frac{1}{2}}$ が何であるかは Boris and Book(1973) には明示さ れていない.} \subsection{計算例} Boris and Book(1973) で行なわれたような一次元の計算を行なう. 解く方程式は1次元移流方程式 \[ \DP{\rho }{t} + \DP{}{x}(\rho u) = 0, \] である. 計算設定を以下にまとめる. \begin{itemize} \item 計算領域は $0\le x \le 1$ をとり, 100個の格子点に分割する. よって格子間隔は $\Delta x = 0.01$ である. \item 速度場は $u=1.0$ で一定とする. \item $\rho $ の初期値は $0\le x \le 0.2$ で $\rho =2.0$ とし, それ以外では $0.5$ とする. \item 境界条件は周期境界条件を与える. \item 時間格子間隔は $\epsilon = u\Delta t/\Delta x=0.2$ となる ようとる. \end{itemize} 時間方向にオイラー法を用いて得られた結果を\Dfigref{fig: 4}に示 す. 上段から順に上流差分, フラックス補正を行なわない SHASTA, SHASTA の各スキームを用いた場合であり,左から順に初期値, 100ステッ プ目, 800ステップ目の図である. 上流差分では拡散が効いてしまい時間 とともに初期値分布からのずれは大きくなる. 補正を行なわない SHASTA では\Deqref{eq: 5}式からも明らかであるが波数に依存しない拡散を持つ ため, 上流差分よりも初期値とのずれが大きくなる. 補正を行なった FCT-SHASTA スキームでは初期分布が比較的よく保たれていることがわか る\footnote{Boris and Book(1973) では脚注2に示したようにさらに細か い補正を行なうことでより高精度の結果を得ている(Boris and Book(1973), fig.3, fig.4 参照).}. \begin{figure}[p] \begin{center} \Depsf[][20cm]{ps-fig/fig4.ps} \caption{一次元の計算例. 上から順に上流差分, フラックス補正を行 なわないSHASTA, SHASTA の各スキーム. 左から順に初期値, 100ステッ プ目, 800ステップ目の結果を示す. 各図の点線は解析解を示す.} \Dfiglab{fig: 4} \end{center} \end{figure}