第 15 回数値流体力学シンポジウム
D08-4
重合格子法を用いた球面座標系での流体計算
Fluid simulation on Spherical coordinate system by using Overset grid
吉田 浩, 東京工業大学院, 〒 152-8550 東京都目黒区大岡山 2-12-1, E-mail:[email protected]
青木 尊之, 東京工業大学 学術国際情報センター, 〒 152-8550 東京都目黒区大岡山 2-12-1, E-mail:[email protected]
Hiroshi Yoshida, Univ.of Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-ku, Tokyo 152-8550 Japan
Takayuki Aoki, GSIC.Univ.of Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-ku,Tokyo 152-8550 Japan
We develop a numerical scheme of fluid dynamics on spherical coordinate system suitable for meteorological and
ocean simulation of the global model. In order to keep higher-order accuracy, Interpolated Differential Operator
(IDO) scheme with Hermite interpolation is used. Introduction of the overset grid method avoids the problems
appearing at the pole, where the singular point exists and the mesh concentrates. The overset grid requires
communications with high accurate interpolation at the boundary between grid systems. We estimate the case I,
II and V of Williamson’s benchmark tests, and obtain quite good results.
1.
はじめに
従来、全球モデルの気象・海洋シミュレーションでは
座標系として球面座標系が使われてきた。しかし、球面
座標系では解像度が上がるにつれ北極、および南極付近
の格子が過剰に集中してしまう。極近傍に格子が集中す
ることは、CFL 条件による制限からタイムステップを非
常に小さくしなければならないという困難を生じる。そ
の結果、無駄な計算が増え、計算時間が大幅に増大する。
また、一般に球面座標系の支配方程式中には 1/ cos θ(θ
は緯度)の項が含まれるが、極点 θ = π/2, −π/2 のとき
無限大になるので計算が行えない。
計算スキームでは差分法や球面調和関数を用いた擬ス
ペクトル法が広く使われてきた。しかし、差分法は高速
に計算できるが一般に精度が低く、長時間スケールなグ
ローバルシミュレーションには向かない。また、球面調
和関数を用いた疑スペクトル法は精度が高いが計算時間
が非常にかかる。
そこで、本研究では従来の方法とは異なる計算スキー
ム、座標系を用いて、全球モデルでの新しい計算手法の開
発を行う。計算スキームにはエルミート補間を用いた「局
所微分オペレータ法 (IDO 法)」(1−4) を採用する。IDO
法では CIP 法 (5−6) と同様に、物理量とその空間 1 階微
係数を従属変数として持つ。補間関数の特徴を活かして、
移流項、非移流項の両者とも高精度で近似ができる。新た
な格子モデルの構築には、元の球面座標系とそれに直交
する補助的な球面座標系を用意し、元の座標系の極付近
の格子を補助座標系の赤道付近の格子で置き換える「重
合格子法」を用いる。(Fig.1)座標系の相互参照には物
理量とその 1 階空間微係数を従属変数として用いた CIP
(Cubic Interporated Propagation)法を利用して 3 次精
度補間で行う。
そして、これらを用いて Williamson らの提唱する球
面上浅水波方程式のテストケース (7−8) 1,2,5 を解きその
有用性を確かめる。
1
Fig. 1: Schematic view of the overset grid. The gridpoints of the neighborhood of a pole in the original
coordinate system are substituted with the equatorial
gridpoints in the auxiliary coordinate system.
2.
2.1
局所微分補間オペレータ法(IDO)
補間関数
IDO 法では各グリッド点毎に物理量とその空間 1 階微
係数を定義し、それらの値を用いて局所的にエルミート
補間を形成する。ただし、移流項、非移流項で補間関数
の形成方法が異なる。計算の際に必要な高次空間微係数
は補間関数を微分することで求める。
移流相に関してはグリッド点における移流速度から風
上方向を検出し、風上方向に対して補間関数を張る。こ
c 2001 by JSCFD
Copyright °
れを風上補間と呼ぶ。補間関数は以下のようになる。
F (x) = ax3 + bx2 + fx,i x + fi
偏微分方程式に含まれる微分項は移流項と非移流項に
分かれるのでそれぞれの項に対し上で述べた補間を適用
する。
(1)
移流速度 u が正のとき i − 1 番目から i 番目のグリッドに
かけて補間関数を張る。ここで i は格子点の番号、fi , fx,i
はそれぞれ格子点 i における物理量と空間 1 階微係数を
表す。高次空間微係数は fxxx = Fxxx (0) = 6a、fxx =
Fxx (0) = 2b、と補間関数を微分することで求める。係数
a, b は適合条件 F (−∆x) = fi−1 、Fx (−∆x) = fx,i−1 か
ら求める。格子点 i における高次空間微係数は
µ
¶
fi − fi−1
fx,i + fx,i−1
fxxx,i = 6 −2
(2)
+
∆x3
∆x2
µ
¶
fi − fi−1
2fx,i + fx,i−1
fxx,i = 2 −3
(3)
+
∆x2
∆x
2.2
IDO 法では物理量とその空間 1 階微分を従属変数とし
て持っているので、時間前進もこれらについて行う必要
がある。時間前進スキームは時間に関するテーラー展開
が用いられる。
∆t2 n
f
2 tt
∆t3 n ∆t4 n
∆t5 n
(11)
+
fttt +
ftttt +
f
6
24
120 ttttt
∆t2 n
n
fxn+1 = fxn +ftx
∆t+
f
2 ttx
∆t4 n
∆t5 n
∆t3 n
+
ftttx +
fttttx +
f
(12)
6
24
120 tttttx
ここで上付きの添え字 n, n + 1 はそれぞれ時刻 t, t + ∆t
での値、ft は時間に対する微係数を表している。f n , fxn
の時間微分は支配方程式を用いて、すべて空間微分に置
き換える。例えば定常速度場における 1 次元の移流方程
式 ft + ufx = 0 の場合、ftt = u2 fxx , fttx = u2 fxxx ,
fttt = u3 fxxx , · · · となり時間微分はすべて空間微分で置
き換わる。
2.3 2 次元の場合
f n+1 = f n +ftn ∆t+
となる。移流速度が負のときは i 番目から i+1 番目のグリ
ッドにかけて補間関数を張る。適合条件は F (∆x) = fi+1 、
Fx (∆x) = fx,i+1 となり、高次空間微係数は
µ
¶
fi − fi+1
fx,i + fx,i+1
fxxx,i = 6 2
+
(4)
∆x3
∆x2
¶
µ
fi − fi+1
2fx,i + fx,i+1
fxx,i = 2 −3
−
(5)
∆x2
∆x
となる。
非移流項に関してはグリッド点を中心として 5 次多項
式で補間関数を張る。これを 5 次中心補間と呼ぶ。補間
関数は以下のようになる。
F (x) = ax5 + bx4 + cx3 + dx2 + fx,i x + fi
時間前進
球面座標系は 2 次元座標系なので IDO 法も 2 次元化す
る必要がある。2 次元問題として以下のような方程式を
考える。
(6)
5 次中心補間では i − 1 番目から i + 1 番目のグリッド
にかけて補間関数を張る。適合条件は F (−∆x) = fi−1 、
Fx (−∆x) = fx,i−1 、F (∆x) = fi+1 、Fx (∆x) = fx,i+1
となる。格子点 i において fxxxxx = Fxxxxx (0) = 120a、
fxxxx = Fxxxx (0) = 24b、fxxx = Fxxx (0) = 6c、fxx =
Fxx (0) = 2d から高次空間微係数は
µ
3(fi+1 − fi−1 )
fxxxxx,i = 120 −
4∆x5
¶
(fx,i+1 + 4fx,i + fx,i−1 )
+
(7)
4∆x4
µ
3(fi+1 − fi + fi−1 )
fxxxx,i = 24 −
2∆x4
¶
(fx,i+1 − fx,i−1 )
+
(8)
4∆x3
µ
5(fi+1 − fi−1 )
fxxx,i = 6
4∆x3
¶
(fx,i+1 + 8fx,i + fx,i−1 )
−
(9)
4∆x2
µ
(fi+1 − 2fi + fi−1 )
fxx,i = 2
∆x2
¶
(fx,i+1 − fx,i−1 )
−
(10)
4∆x
となる。
2
ft = fx + fy
(13)
1 次元の場合と同様に与えられた方程式を空間微分する
ことで不足分の方程式を生成すると以下のようになる。
ftx
=
fxx + fxy
(14)
fty
=
fxy + fyy
(15)
fxx , fyy は 1 次元のケースと同様に求めることができる
が fxy は補間関数が作れないため求めることができない。
そこで、fxy もまた従属変数として定義する。従属変数
を増やしたことで ftxy = fxxy + fxyy を解く必要がある。
fxx が f の x 方向の補間関数から求められるのと同様
に、fxxy 、fxyy はそれぞれ fy の x 方向、fx の y 方向の
補間関数から求めることができる。すなわち、1 次元の
エルミート補間関数を f の x および y 方向、fx の y 方
向、fy の x 方向と 4 回適用することで必要な微係数を得
ることができる。また、fxxyy のような微係数は計算精度
にほとんど影響しないので 0 とみなして計算する。
3. 重合格子法
極付近でおこる問題を回避するため重合格子法を用い
る。元の球面座標系 (メイングリッド) の他にもう 1 つ、
緯度方向に π/2 だけ回転した球面座標系 (サブグリッド)
を用意し、メイングリッドの極付近をサブグリッドの赤
道付近で置き換える。(Fig.1)
c 2001 by JSCFD
Copyright °
3.2
サブグリッドはメイングリッドを回転しただけで形状
は変化しないので解く式の形はメイングリッドのものと
同じである。しかし、サブグリッドでの初期条件の与え
方とそれぞれのグリッドでの境界の時間前進に関しては
注意が必要である。
3.1
境界点の補間
メイングリッド、サブグリッドとも毎ステップごとに
境界点での物理量、その空間 1 階微係数を求める必要が
ある。ここでは境界の時間前進を CIP-Type C(6) の手法
で補間して求める。Fig.2 はサブグリッドの境界点(●)
での物理量、その空間 1 階微係数をそれを囲むメイン
グリッドの物理量、その空間 1 階微係数(■)を用いて
補間する例である。(i, j)、(i + 1, j) の間に 3 次のエル
(i, j + 1)
一般に座標変換によって密度や温度などのスカラー量
の値は変わらないが、微分演算子や速度などのベクトル
量は値が変わってしまう。よって、サブグリッドでの初
期条件や補間で求めた境界点での値がベクトル量の場合、
座標変換を施す必要がある。サブグリッドの点における
物理量またその空間 1 階微係数を f 0 , fλ0 , fθ0 , fλ0 θ0 とする
と chain rule を用いることで以下のようになる。
f0 = f
(16)
∂λ
∂θ
fλ0 =
fλ +
fθ
(17)
∂λ0
∂λ0
∂λ
∂θ
fθ0 = 0 fλ + 0 fθ
(18)
∂θ
∂θ
∂λ ∂λ
∂θ ∂θ
∂λ ∂θ
fλ0 θ0 =
fλλ +
fθθ +
fλθ
0
0
0
0
∂λ ∂θ
∂λ ∂θ
∂λ0 ∂θ0
∂θ ∂λ
∂2λ
∂2θ
(19)
+ 0 0 fλθ +
fλ +
fθ
0
0
∂λ ∂θ
∂λ ∂θ
∂λ0 ∂θ0
(i + 1, j + 1)
B
ここで、(λ, θ)、(λ0 , θ0 ) はそれぞれメイングリッド、サブ
グリッドでの経度、緯度を表す。両者の関係は以下のよ
うになる。
·
¸
cos θ sin λ
0
λ = arctan
(20)
sin θ
T
C
座標変換
D
θ0
= − arcsin [cos λ cos θ]
·
η
ζ
(i, j )
A
λ
=
θ
=
cos θ0 sin λ0
− arctan
sin θ0
0
arcsin [cos λ cos θ0 ]
(21)
¸
(22)
(23)
ところで、式 (19) では fλ0 θ0 を求めるのに fλλ , fθθ といっ
た空間 2 階微分の値が必要となる。fθθ は Fig.2 の点 A,
点 B 間の補間関数から、
(i + 1, j )
fθθ = Fθθ (η) = 6aη + 2b
Fig. 2: Illustration of the boundary interpolation.
(24)
と求め fλλ は点 C、点 D 間の補間関数から、
ミート補間関数 (F (λ) = aλ3 + bλ2 + fλ(i,j) λ + f(i,j) 、
Fθ (λ) = aλ3 + bλ2 + fλθ(i,j) λ + fθ(i,j) ) を形成する。適合
条件より a, b を求める。得られた補間関数と補間関数の微
分に λ = ξ を代入することにより点 A での f, fλ , fθ , fλθ
を求める。(i, j +1)、(i+1, j +1) 間も同様にして、点 B で
の f, fλ , fθ , fλθ を求める。さらに点 A、点 B 間に補間関数
(F (θ) = aθ3 +bθ2 +fθ θ+f 、Fλ (θ) = aθ3 +bθ2 +fλθ θ+fλ )
を形成し、適合条件より a, b を求める。得られた補間関
数またその微分に θ = η を代入することにより点 T での
f, fλ , fθ , fλθ が求まる。
メイングリッドの境界点での物理量、その空間 1 階微
係数も同様の手順でサブグリッドの物理量、その空間 1
階微係数を用いて補間して求める。
fλλ = Fλλ (ξ) = 6aξ + 2b
と求める。
メイングリッド、サブグリッド上における移流速度を
(u, v)、(u0 , v 0 ) とすると両者の関係は以下のようになる。
u0
=
v0
=
∂λ
∂λ
+v
∂λ
∂θ
∂θ
∂θ
u
+v
∂λ
∂θ
u
(26)
(27)
∂u0
∂v 0 ∂v 0
∂v 0
∂u0 ∂u0
、
、
、
、
、
も求
∂λ0 ∂θ0 ∂λ0 θ0 ∂λ0 ∂θ0 ∂λ0 θ0
める必要がある。
サブグリッドの初期条件や境界点での値は式 (16)-(19)、
(26)、(27) から求めることができる。
IDO 法では
4.
4.1
以上の作業で●における f, fλ , fθ , fλθ が求まるが、こ
れらはメイングリッドでの値のためサブグリッドでの値
に変換する必要がある。
(25)
数値実験
IDO 法での精度検証
最初に、座標系は通常の球面座標系として IDO 法の精
度を検証する。
3
c 2001 by JSCFD
Copyright °
離散化は球面座標系を直交座標系にマッピングすること
で行う (Fig.3)。経度方向は周期境界、緯度方向は極を固
定境界とし、極での値は時間的に変化しないものとする。
θ
λ
Fig. 3: Mapping of the spherical coordinate system into
the cartesian coordinate system
検証問題として Williamson のテストケース 1、2、5 を
解く。支配方程式は、
u
v
uφ + uθ
a cos θ
a
u tan θ
g
−(f +
)v +
hφ = 0
a
a cos θ
v
u
vt +
vφ + vθ
a cos θ
a
u tan θ
g
+(f +
)u + hθ = 0
a
a
v
u
hφ + hθ
ht +
a cos θ
a
h
+
(uφ + (v cos θ)θ ) = 0
a cos θ
ut +
I [|h∗ − h∗T |]
(36)
I [|h∗T |]
PP
1
ただし I(h∗) = 4π
h cos θdθdλ である。ここで h∗
は数値解、h∗T は解析解を表す。
Fig.4 は 12 日後におけるプロファイル(格子点数 180 ×
135)の θ = 0 での断面図を 1 次精度風上差分(FDM1)、
3 次精度風上差分(FDM3)、IDO 法で比較した結果であ
る。FDM と比較して、IDO は高精度でプロファイルを
維持できていることがわかる。
また、Fig.5 は各スキームの空間精度の比較した結果
である。解像度を高くすることで、FDM ではそれほど
精度が上がらないのに対し IDO 法では精度がよくなって
いる。
E(h∗ )
(28)
(29)
u =
u0 (cos θ cos α + sin θ cos λ sin α)
(37)
v
=
(38)
h∗
=
−u0 sin λ sin α
µ
¶ µ ¶2
1
u20
f
h0 −
aΩu0 +
g
2
2Ω
(30)
で与えられる。ここで、h は自由表面の高さ、h は流体の
深さを表す。もし、hs を山の高さとすると h = h∗ +hs で
ある。u, v はそれぞれ経度(λ)、緯度(θ)方向の速度、a(=
6.37122 × 106 [m]) は地球の半径、g(= 9.80616[m/s−2 ])
は重力加速度、f はコリオリパラメータである。
u0 (cos θ cos α + sin θ cos λ sin α)
(31)
v
−u0 sin λ sin α
(32)
4.1.3 テストケース 5 テストケース5では中緯度に高
い山を置き、それを過ぎる流れの時間発展を追う問題で
ある。このテストの目的は、質量やエネルギーといった
保存量がどの程度保存できるかを検証することである。
初期条件はケース2のときと同じだが、山を配置するこ
とにより hs がゼロでなくなるので、支配方程式がケース
2のときと異なる。保存誤差は以下の式で評価する。
で与える。ここで α は回転軸と座標軸の角度を表す。式
(30) に代入すると
∂h∗
u ∂h∗
v ∂h∗
=−
−
∂t
a cos θ ∂λ
a ∂θ
と簡単な移流方程式になる。初期プロファイルは

 h0 (1 + cos( πr )) if r < R
h=
2
R

0 if r ≥ R
(39)
各パラメータは f = 2Ω(− cos λ cos θ sin α + sin θ cos α)、
Ω = 7.292×10−5 、u0 = 2πa/(12day)、h0 = 2.998×103 、
α = 0.0 とする。誤差の評価は式 (36) で行う。
Fig.6 は格子点数 (a)60 × 45、(b)120 × 90、(c)180 ×
135 での誤差、Fig.7 は各格子点での誤差の平均を IDO
と FDM(移流項:1次風上差分、非移流項:2 次中心差
分)で比較した結果である。FDM に比べ IDO は非常に
よい結果を得ることができた。これは IDO 法が非移流項
を高精度に近似できることと、時間、空間の 1 階微分を
解析的に持っていることが要因だと考えられる。
4.1.1 テストケース 1 テストケース 1 は移流の問題で
ある。12 日で地球上を一周する速度場を与え、cosine 型
のプロファイルをどれほど維持できるかを検証する問題
である。定常速度場を
u =
=
4.1.2 テストケース 2 テストケース 2 はケース 1 で用
いた定常流のテストである。時間発展しないような初期
プロファイルを与え、どの程度の誤差で定常解を保つこ
とができるかを検証する問題である。このテストでは3
つの式すべてを解く。ただし、山は存在しない(hs = 0)
ので h = h∗ となる。初期条件 (かつ解析解) は以下の式
で与えられる。
∗
=
また誤差の評価は以下の式で行う。
(33)
(34)
r = a arccos [sin θc sin θ + cos θc cos θ cos(λ − λc )] (35)
で与える。ここで各パラメータは λc = 3π/2、θc = 0、
h0 = 1000、R = a/3、2πa/(12days)、α = π/6 である。
4
I(ξ(t)) =
I [ξ (t)] − I [ξ (0)]
I [ξ (0)]
(40)
ここで、
ξ
=
h∗ (質量)
(41)
c 2001 by JSCFD
Copyright °
ξ
=
1 ∗ 2
h (u + v 2 )
2
1
+ g(h2 − h2s )(エネルギー)
2
参考文献
(42)
である。
Fig.8 は格子点数 128 × 64 での (a) 初期、(b)5 日後、
(c)10 日後、(d)15 日後の自由表面(h)の等高線図(50m
間隔)である。また、Fig.9 は保存誤差を FDM(移流項:
3 次風上差分、非移流項:2 次中心差分)と IDO で比較
した結果である。FDM に比べ IDO は保存性が高いこと
がわかる。
4.2
格子モデルの検証
新 格 子 モ デ ル で は メ イ ン グ リッド の 緯 度 の 範 囲 を
−π/4 ≤ θ ≤ π/4 としそれ以外の範囲をサブグリッド
で置き換えた。
検証問題としてテストケース1を用い球面座標系の場
合と比較する。ただし、新格子モデルで2つの格子を通
過するようにパラメータ条件は α = π/3 とした。
Fig.10 は計算時間と精度 (式 (36) で評価) を比較した
結果である。精度に関しては2つの格子間を補間してい
るにもかかわらず、球面座標系の場合とほぼ変わらない
ことがわかる。また、極付近の格子集中を回避したこと
により、計算時間を大幅に減らすことができた。
Fig.11 は極を通過する条件 (α = π/2) で解いたときの
プロファイルの様子である。見てわかるように格子間を
またぐときプロファイルをくずすことなく、極も問題な
く通過できることがわかる。
5.
結論
本研究では、球面上で数値計算を行う上での問題点を
回避すべく新しい計算モデルの構築を行い検証した。ま
ず、計算スキームには IDO 法を用い、検証問題として
Williamson のテストケース 1、2、5 を解いた。いずれの
テストケースも非常によい結果を得ることができた。
また、座標系の改善のため重合格子法を用いて極の特
異点と極近傍の格子点集中を回避する格子モデルを作成
した。格子間の相互作用は CIP-Type C を用いて行った。
テストケース 1 を解いた結果、球面座標系の場合と精度
はほぼ変わらずに計算時間を大幅に短縮することができ
た。さらに、極を通過する条件でも問題なく解くことが
できた。
5
(1) T.Aoki:Interpolated Differential Operator (IDO)
scheme for solving partial differential equations.Comput.Phys.Commun,Vol.102(1997) 132146.
(2) K.Sakurai,T.Aoki:Implicit IDO (Interpolated Differential Operator) Scheme, Comput. Fluid Dynamics Journal, Vol.8, No. 1 (1999) 6-12.
(3) K. Sakurai, T. Aoki: A Numerical Procedure to
Evolving Contact Discontinuity by using Moving
Cut-Cell Method, Comput. Fluid Dynamics Journal, Vol.10, No.1 (2001), 76-84.
(4) T. Aoki, S. Nishita, K. Sakurai: Interpolated Differential Operator (IDO) Scheme and Application
to Level Set Method, Comput. Fluid Dynamics
Journal , Vol.9, No.4 (2001) 406-417.
(5) 青木尊之, 肖 鋒: 複雑形状の複雑でない計算法, 情
報処理, Vol. 42, No.6 (2001), 550-556.
(6) T.Aoki:Multi-dimensional Advection of CIP
(Cubic-Interpolated Propagation) Scheme, Comput. Fluid Dynamics Journal , Vol.4, No. 3 (1995)
279-291.
(7) D.L.Williamson et.al: A Standard Test Set for
Numerical Approximations to the Shallow Water
Equations in Spherical Geometry , J.Comput.Phys
, Vol.102(1992)211-224.
(8) R.Jakob-Chen, J.J.Hanck, and D.L.Williamson:
Spectral Transform Solutions to the Shallow Water
Test Set,J.Comput.Phys,Vol.119(1995) 164-187.
c 2001 by JSCFD
Copyright °
(a)
1000
h
α=30[deg]
FDM1
after 12day
analysis solution
500
0
λ
(b)
1000
0
10
h
(60
(180
E
α=30[deg]
135)
(120
90)
45)
FDM3
after 12day
analysis solution
FDM1
FDM3
500
IDO
-2
10
-1
10
0
λ
(c)
1000
Fig. 5: The deviation from the exact solution versus the grid interval of the longitudinal(λ) in the
Case I. ●:1st-order upwind scheme, ▲:3rd-order upwind scheme, ■:IDO scheme.
h
α=30[deg]
after 12day
IDO
analysis solution
500
0
λ
Fig. 4: The cross section of the profile(grid number
180 × 135) in θ = 0 in the Case I. (a)1st-order upwind
scheme,(b)3rd-order upwind scheme,(c)IDO scheme. 6
c 2001 by JSCFD
Copyright °
(a)
-2
10
-4
10
E
-6
10
-8
10
-10
10
FDM Grid.No 60
45
IDO Grid.No 60
45
-12
10
0
1
2
3
4
5
time[day]
(b)
-2
10
-2
10
-4
10
-4
10
(60
E
-6
E
-6
10
45)
(120
90)
10
(180
-8
10
135)
-8
10
-10
10
-10
10
-12
10
0
1
FDM Grid.No 120
90
IDO Grid.No 120
90
2
3
4
FDM
IDO
-12
10
5
100
200
Grid number(longitudinal)
time[day]
Fig. 7: The average of the error(Ē) versus the grid
interval of the longitudinal(λ)in the Case II.
(c)
-2
10
-4
10
E
-6
10
FDM Grid.No 180
135
IDO Grid.No 180
135
-8
10
-10
10
-12
10
0
1
2
3
4
5
time[day]
Fig. 6: Time evolution of the error in the Case II.
Grid.No.(a)60 × 45(b)120 × 90(c)180 × 135.
7
c 2001 by JSCFD
Copyright °
(a)
(b)
(a)
(c)
(b)
(d)
Fig. 9: Time evolution of conservation error (a)mass
(b)energy in the Case IV.
Fig.
8: Contour maps on a rectangular latitude/longitude projection of the height of the free
surface (h) in the Case IV.(a)day 0;(b)day 5;(c)day
10;(d)day 15.The contour interval is 50m
8
c 2001 by JSCFD
Copyright °
(a)
1
(60
(120
E
45)
90)
spherical coordinate system
overset grid system
(180
135)
-2
10
-1
10
(b)
2000
spherical coordinate system
Time[sec]
overset grid system
(180
(60
0
50
45)
(120
100
90)
90)
150
200
Grid number(longitudinal)
Fig.
10:
Comparison of (a)spatial accuracy and
(b)computer performance between spherical coordinate
system and overset grid system.
Fig. 11: Advection of Cosine Bell over the Pole on the
Overset grid system in the Case I.
9
c 2001 by JSCFD
Copyright °
ダウンロード

重合格子法を用いた球面座標系での流体計算