% 表題 移流方程式の差分解法: MPDATA スキーム 数値実験 % % 履歴 1997/11/25 小高正嗣 % 1998/01/20 小高正嗣 % 1998/05/11 小高正嗣 % \section{2 次元 の計算} 実際に 2 次元の移流方程式, \begin{equation} \DP{\psi }{t} + \DP{}{x}(u\psi ) + \DP{}{y}(v\psi ) = 0, \Deqlab{eq: 29} \end{equation} を MPDATA スキームを用いて計算する. 比較のため上流差分, 中心差分, フッラ クス修正法(Flux Corrected Transpot; 以下 FCT と表記)による計算も行なう. \subsection{設定} 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 $\psi $ の初期値は円錐形分布を与える. 底面の円の中心は $(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: 1}に示す. \begin{figure}[bth] \begin{center} \Depsf[][80mm]{ps-fig/init.ps} \caption{計算領域と初期値の分布.} \Dfiglab{fig: 1} \end{center} \end{figure} \newpage \subsection{計算例: 上流差分} 比較のためまず上流差分スキーム\Deqref{eq: 22}式だけを用いて計算を 行なった. スキームを再掲すると以下のようになる. \begin{eqnarray} \psi _{i,j}^{n+1} &=& \psi _{i,j}^{n} - \{F_{x}(\psi _{i,j}^{n},\psi _{i+1,j}^{n},u_{i+\frac{1}{2},j}^{n}) - F_{x}(\psi _{i-1,j}^{n},\psi _{i,j}^{n},u_{i-\frac{1}{2},j}^{n})\} \nonumber \\ && - \{F_{y}(\psi _{i,j}^{n},\psi _{i,j+1}^{n},v_{i,j+\frac{1}{2}}^{n}) - F_{y}(\psi _{i,j-1}^{n},\psi _{i,j}^{n},v_{i,j-\frac{1}{2}}^{n})\}. \Deqlab{eq: 30} \end{eqnarray} ここで, \begin{eqnarray*} F_{x}(\psi _{i,j}^{n},\psi _{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}|)\psi _{i,j}^{n} + (u_{i+\frac{1}{2},j}^{n}-|u_{i+\frac{1}{2},j}^{n}|)\psi _{i+1,j}^{n}] \frac{\Delta t}{2\Delta x}, \\ F_{y}(\psi _{i.j}^{n},\psi _{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}|)\psi _{i,j}^{n} + (v_{i,j+\frac{1}{2}}^{n}-|v_{i,j+\frac{1}{2}}^{n}|)\psi _{i,j+1}^{n}] \frac{\Delta t}{2\Delta y}, \end{eqnarray*} である. これを時間方向に修正オイラー法(2 次のルンゲクッタ)を用いて 計算した. 1 回転後(628ステップ)と 3 回転後(1884ステップ)後の結果 を\Dfigref{fig: 2}と\Dfigref{fig: 3}にそれぞれ示す. 上流差分では 数値拡散が大きいため 1 回転後 で既に初期分布は大きく損なわれている. 3 回転もすると初期分布はあとかたもなくなってしまう. \begin{figure}[p] \begin{center} \Depsf[][80mm]{ps-fig/upstream1.ps} \caption{上流差分による計算. 1 回転(628ステップ)後の結果.} \Dfiglab{fig: 2} \Depsf[][80mm]{ps-fig/upstream2.ps} \caption{上流差分による計算. 3 回転(1884ステップ)後の結果.} \Dfiglab{fig: 3} \end{center} \end{figure} \newpage \subsection{計算例: 2 次中心差分} 次に 2 次の中心差分スキームを用いて計算を行なった. スキームは, \begin{eqnarray} \psi _{i,j}^{n+1} &=& \psi _{i,j}^{n} - \{F_{x}(\psi _{i,j}^{n},\psi _{i+1,j}^{n},u_{i+\frac{1}{2},j}^{n}) - F_{x}(\psi _{i-1,j}^{n},\psi _{i,j}^{n},u_{i-\frac{1}{2},j}^{n})\} \nonumber \\ && - \{F_{y}(\psi _{i,j}^{n},\psi _{i,j+1}^{n},v_{i,j+\frac{1}{2}}^{n}) - F_{y}(\psi _{i,j-1}^{n},\psi _{i,j}^{n},v_{i,j-\frac{1}{2}}^{n})\}. \Deqlab{eq: 31} \end{eqnarray} ただしここでの $F_{x}, F_{y}$ は \begin{eqnarray} F_{x}(\psi _{i,j}^{n},\psi _{i+1,j}^{n},u_{i+\frac{1}{2},j}^{n})&=& u_{i+\frac{1}{2},j}^{n}(\psi _{i,j}^{n} + \psi _{i+1,j}^{n}) \frac{\Delta t}{\Delta x}, \Deqlab{eq: 32} \\ F_{y}(\psi _{i.j}^{n},\psi _{i,j+1}^{n},v_{i,j+\frac{1}{2}}^{n})&=& v_{i,j+\frac{1}{2}}^{n}(\psi _{i,j}^{n} + \psi _{i,j+1}^{n}) \frac{\Delta t}{\Delta y}, \Deqlab{eq: 33} \end{eqnarray} である. 時間方向には同様に修正オイラー法を用いて計算した. ただし時間 格子間隔を $\Delta t=0.1$ として計算すると計算不安定を起こしてしまっ たので, ここでは $\Delta t=0.05 $ とした. また壁での波の反射を抑える ため 壁際の 1 グリッドでは $\psi $ の 1 次に比例した減衰(減衰係数は 5 )を加えてある. 1 回転後(1256ステップ)と 3 回転後(3768ステップ)後の結果 を\Dfigref{fig: 4}と\Dfigref{fig: 5}にそれぞれ示す. 中心差分では初期 分布の形状はそれなりに維持されるが, 進行方向の後側に波が発生しそれに 伴い負の値が生じていることがわかる. \begin{figure}[p] \begin{center} \Depsf[][80mm]{ps-fig/center2-1.ps} \caption{中心差分による計算. 1 回転(1256ステップ)後の結果.} \Dfiglab{fig: 4} \Depsf[][80mm]{ps-fig/center2-2.ps} \caption{中心差分による計算. 3 回転(3768ステップ)後の結果.} \Dfiglab{fig: 5} \end{center} \end{figure} \subsection{計算例: MPDATA} MPDATA (\Deqref{eq: 22}$\sim $ \Deqref{eq: 25}式)を用いて同様の 計算を行なった. なお $\Delta t=0.1$ に戻してある. 途中必要となる 平均値 計算は\Deqref{eq: 26}$\sim $ \Deqref{eq: 28}に従い求めた. これらの平均操作については Smolarkiewicz(1982) の記述を参考にしてある. 1 回転後(628ステップ)と 3 回転後(1884ステップ)後の結果 を\Dfigref{fig: 6}と\Dfigref{fig: 7}にそれぞれ示す. 初期分布は比較的 よく保たれている. 2 次中心差分で発生したような波も生じていない. \begin{figure}[p] \begin{center} \Depsf[][80mm]{ps-fig/mpdata1.ps} \caption{MPDATA による計算. 1 回転(628ステップ)後の結果.} \Dfiglab{fig: 6} \Depsf[][80mm]{ps-fig/mpdata2.ps} \caption{MPDATA による計算. 3 回転(1884ステップ)後の結果.} \Dfiglab{fig: 7} \end{center} \end{figure} \newpage \subsection{計算例: FCT} FCT を用いて同様の計算を行なった. FCT の一般形は例えば 2 次元の場合 次のように表される. \begin{eqnarray} \psi _{i,j}^{n+1} &=& \psi _{i,j}^{n} - \{[CF_{x}\cdot FH_{x}+(1-CF_{x})\cdot FL_{x}]_{i+\frac{1}{2},j} \nonumber \\ && - [CF_{x}\cdot FH_{x}+(1-CF_{x})\cdot FL_{x}]_{i-\frac{1}{2},j}\} \nonumber \\ && - \{[CF_{y}\cdot FH_{y}+(1-CF_{y})\cdot FL_{y}]_{i,j+\frac{1}{2}} \nonumber \\ && - [CF_{y}\cdot FH_{y}+(1-CF_{y})\cdot FL_{y}]_{i,j-\frac{1}{2}}\} \Deqlab{eq: 34} \end{eqnarray} ここで $FH_{x}, FH_{y}$ は高次精度(例えば 2 次中心差分)で求めたフ ラックス, $FL_{x}, FL_{y}$ は低次精度(例えば上流差分)で求めたフラッ クス, $CF_{x}, CF_{y}$ は補正係数(Corrective Factor)と呼ばれる係 数である. FCT は Boris and Book(1973, 1976), Book et al.(1975)に よって提案され, Zalesak(1979) によって一般形が示された. 補正係数 の決めかたは Zalesak(1979) にその詳細が述べられている. ここでは高 次精度に 2 次中心差分, 低次精度に上流差分を用いた. まず補正係数を $CF_{x}= CF_{y}=0.8$ と固定して計算を行なった. 1 回転後(628ステップ)と 3 回転後(1884ステップ)後の結果を \Dfigref{fig: 8}と\Dfigref{fig: 9}にそれぞれ示す. 上流差分と比べ 初期分布は比較的よく保たれている. 2 次中心差分で発生したような波 も生じない. しかし上流差分による数値拡散が時間とともに無視できなく なる. 続いて補正係数をきちんと毎ステップ計算するようにして同様の計算を 行なった. 1 回転後(628ステップ)と 3 回転後(1884ステップ)後の結果 を\Dfigref{fig: 10}と\Dfigref{fig: 11}にそれぞれ示す. 全体的に数 値拡散は抑えられているが, 山の頂上付近が削られてしまい ``clipping''と呼ばれる現象が回避できていないことを示している. \begin{figure}[p] \begin{center} \Depsf[][80mm]{ps-fig/fct1.ps} \caption{補正係数一定の FCT による計算. 1 回転(628ステップ)後の結果.} \Dfiglab{fig: 8} \Depsf[][80mm]{ps-fig/fct2.ps} \caption{補正係数一定の FCT による計算. 3 回転(1884ステップ)後の結果.} \Dfiglab{fig: 9} \end{center} \end{figure} \begin{figure}[p] \begin{center} \Depsf[][80mm]{ps-fig/fct3.ps} \caption{FCT による計算. 1回転(628ステップ)後の結果.} \Dfiglab{fig: 10} \Depsf[][80mm]{ps-fig/fct4.ps} \caption{FCT による計算. 3回転(1884ステップ)後の結果.} \Dfiglab{fig: 11} \end{center} \end{figure}