計算例

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

• 計算領域に $0\leq x \leq 1, 0\leq y \leq 1$ をとる. それぞれの方向に 100 個の格子に分割する. すなわち $\Delta x = \Delta y = 0.01$ である.

• 速度場は $(x_{0}, y_{0})=(0.5, 0.5)$ を中心に角速度 $\omega = 0.1$ で剛体回転する場を与える. すなわち $(u, v) = ( - \omega (y-y_{0}), \omega (x-x_{0}))$ である.

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

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

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

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

図 7: 計算領域と初期値の分布.

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

$$\rho _{i,j}^{n+1} = \rho _{i,j}^{n} - \{F_{x}(\rho _{i,j}^{n},\rho _{i+1,j}^{n},u_{i+... ...2},j}^{n}) - F_{x}(\rho _{i-1,j}^{n},\rho _{i,j}^{n},u_{i-\frac{1}{2},j}^{n})\}$$

$$- \{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})\}.$$ (34')

ここで,

$$\begin{eqnarray*} F_{x}(\rho _{i,j}^{n},\rho _{i+1,j}^{n},u_{i+\frac{1}{2},j}^{... ...}{2}}^{n}\vert)\rho _{i,j+1}^{n}] \frac{\Delta t}{2\Delta y}, \end{eqnarray*}$$

である. これを時間方向に修正オイラー法(2 次のルンゲクッタ)を用いて 計算した.

1回転後(628ステップ)と3回転後(1884ステップ)後の結果をFig.とFig.9にそれぞれ示す. 上流差分では数値拡散が大きいた め1回転後で既に初期分布は大きく損なわれている. 3回転もすると初期 分布はあとかたもなくなってしまう.

図 8: 上流差分による計算. 1回転(628ステップ)後の結果.
図 9: 上流差分による計算. 3回転(1884ステップ)後の結果.

同様の計算を FCT を用いて行なった. 低次のスキームに上流差分, 高次 のスキームに2次中心差分を用いている. 1回転後(628ステップ)と3回転後 (1884ステップ)後の結果をFig.10とFig.11にそれ ぞれ示す. 全体的に数値拡散は抑えられているが, 山の頂上付近が削られ てしまい clipping が回避できていないことを示している.

図 10: FCT による計算. 1回転(628ステップ)後の結果.
図 11: FCT による計算. 3回転(1884ステップ)後の結果.

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

図 12: 計算領域と初期値の分布. 円筒分布を与えた場合.
図 13: FCT による計算. 1回転(628ステップ)後の結果. 初期値に円筒分布 を与えた場合