マントル対流シミュレーション手法に関する補足情報
亀山 真典
平成 22 年 8 月 21 日

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 1

基礎方程式系 (無次元、非弾性近似)➢基礎方程式系➢格子分割と変数配置➢格子分割と座標系の設定➢球座標系での基礎方程式➢SIMPLE 系解法➢べき乗法➢1 次元定常移流–拡散方程式の解➢各種「上流化」法
➢Gershgorin の定理➢行列演算からみた
CFL 条件➢Gershgorin の定理
による CFL➢格子界面での熱伝導率➢格子界面での剪断粘性率➢参考文献

以下、˜ で無次元化された量を示す❏ 運動方程式 (Navier-Stokes 方程式 + 連続の式)˜ +∇˜ · τ̃ − Ra T̃ g̃,
0 = −∇p̃˜ · (ρ̃ṽ) = 0∇

❏ 粘性流体の構成方程式


1 ˙
˙
τ̃ = 2η̃ ε̃ − tr(ε̃I )I32 ˜˜˜∇ · ṽ I= η̃ ∇ ⊗ ṽ + ṽ ⊗ ∇ −3α∆T gd3: Rayleigh 数
ただし、 Ra ≡ηκ2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 2

基礎方程式系 (無次元、非弾性近似)➢基礎方程式系➢格子分割と変数配置➢格子分割と座標系の設定➢球座標系での基礎方程式➢SIMPLE 系解法➢べき乗法➢1 次元定常移流–拡散方程式の解➢各種「上流化」法
➢Gershgorin の定理➢行列演算からみたCFL 条件➢Gershgorin の定理
による CFL➢格子界面での熱伝導率➢格子界面での剪断粘性率➢参考文献❏ 熱輸送方程式

∂ T̃
˜ · ṽ T̃ − k̃ ∇˜ T̃ += −∇∂ t̃
|{z}
移流+拡散Di α̃T̃ (ṽ · g̃)|{z}

静水圧下での断熱温度変化

2ρHdDi2 ˜2˙++∇ · ṽ
η̃ 2(ε̃II ) −3{z } |Ra|k∆T{z}2

内部発熱

粘性散逸αgd
ただし、 Di ≡: 散逸数 (運動から熱への変換効率に相当)
Cppε̇II = ε̇211 + ε̇222 + ε̇233 + 2ε̇212 + 2ε̇223 + 2ε̇213: 歪速度テンソルの第 2 不変量

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 3

格子分割と変数配置➢基礎方程式系➢格子分割と変数配置➢格子分割と座標系の設定➢球座標系での基礎方程式➢SIMPLE 系解法➢べき乗法➢1 次元定常移流–拡散方程式の解➢各種「上流化」法
➢Gershgorin の定理➢行列演算からみたCFL 条件➢Gershgorin の定理
による CFL➢格子界面での熱伝導率➢格子界面での剪断
粘性率
➢参考文献

Staggered grid に基づく有限体積法
δx(I)ydx(i-1)

dx(i)

J+1
yp(j)
yc(J)

jdy(j)
δy(J)J

解くべき変数のうち、❏ 点 (i, j) でスカラー量j-1dy(j-1)J-1(温度、圧力、· · · )❏ 点 (I, j) で x 方向流速I-1i-1Ii

I+1xxc(I)❏ 点 (i, J) で y 方向流速を定義する。

xp(i)

点 (i, j) は有限体積の中心にとる。このため
dx(i−1) + dx(i)のようにとっている。
δx(I) =2
2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 4

格子分割と座標系の設定

yp(j=2)yp(j=Ny+1)

yc (J=Ny+2)(3 次元の場合も同様)

yc(J=2)➢基礎方程式系➢格子分割と変数配置➢格子分割と座標系の設定➢球座標系での基礎方程式➢SIMPLE 系解法➢べき乗法➢1 次元定常移流–拡散方程式の解➢各種「上流化」法
➢Gershgorin の定理➢行列演算からみた
CFL 条件➢Gershgorin の定理
による CFL➢格子界面での熱伝導率➢格子界面での剪断粘性率
➢参考文献

yx

xc(I=2)
2010 年夏の GFD セミナー

xp(i=2)xp(i=Nx+1)

xc(I=Nx+2)

平成 22 年 8 月 21 日 – slide 5

球座標系での基礎方程式系 (1)➢基礎方程式系➢格子分割と変数配置➢格子分割と座標系の設定➢球座標系での基礎方程式➢SIMPLE 系解法➢べき乗法➢1 次元定常移流–拡散方程式の解➢各種「上流化」法
➢Gershgorin の定理➢行列演算からみたCFL 条件➢Gershgorin の定理
による CFL➢格子界面での熱伝導率➢格子界面での剪断粘性率➢参考文献

有限体積法に則って離散化をしたいので、基礎方程式はなるべく保存形で書くようにする。
熱輸送方程式をあらわに書くと以下の通り。

∂Tfr = vr T − κ∂r
1 ∂Tfθ = vθ T − κr ∂θ1 ∂Tfφ = vφ T − κr sin θ ∂φ1 ∂1 ∂ 2 1 ∂fφ∂Tr fr += 2(sin θ fθ ) +∂tr ∂rr sin θ ∂θr sin θ ∂φ+qただし簡単のため、Boussinesq 近似 (断熱温度変化なし、粘性散逸なし) を仮定してある。

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 6

球座標系での基礎方程式系 (2)➢基礎方程式系➢格子分割と変数配置➢格子分割と座標系の設定➢球座標系での基礎方程式➢SIMPLE 系解法➢べき乗法➢1 次元定常移流–拡散方程式の解➢各種「上流化」法
➢Gershgorin の定理➢行列演算からみた
CFL 条件➢Gershgorin の定理
による CFL➢格子界面での熱伝導率➢格子界面での剪断粘性率➢参考文献

歪速度テンソルの対角成分をあらわに書くと以下の通り。

∂vrε̇rr =∂r
∂vrvr 1 ∂vθ
1 1 ∂ 21 ∂vθ
+=(r vr ) −ε̇θθ =+2rr ∂θ2 r ∂r∂rr ∂θvr cot θ1 ∂vφ
ε̇φφ =+vθ +rrr sin θ ∂φ
1 1 ∂ 2∂vr=(r vr ) −22 r ∂r
∂r


1 ∂vφ
11 ∂∂vθ++(sin θvθ ) −r sin θ ∂θ∂θr sin θ ∂φ
粗い空間解像度 (多重格子法で登場する) であっても∇ · v = ε̇rr + ε̇θθ + ε̇φφ を厳密に満たすよう工夫している。

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 7

球座標系での基礎方程式系 (3)➢基礎方程式系➢格子分割と変数配置➢格子分割と座標系の設定➢球座標系での基礎方程式➢SIMPLE 系解法➢べき乗法➢1 次元定常移流–拡散方程式の解➢各種「上流化」法
➢Gershgorin の定理➢行列演算からみたCFL 条件➢Gershgorin の定理
による CFL➢格子界面での熱伝導率➢格子界面での剪断
粘性率
➢参考文献

歪速度テンソルの非対角成分をあらわに書くと以下の通り。

∂  vθ  1 ∂vr2ε̇rθ = 2ε̇θr = r+∂r rr ∂θ1 ∂vr∂  vφ +2ε̇rφ = 2ε̇φr = r∂r rr sin θ ∂φ1 ∂vθ
sin θ ∂  vφ +2ε̇θφ = 2ε̇φθ =r ∂θ sin θr sin θ ∂φ
これらと粘性率の積をとれば、応力テンソルの 9 成分の表式を得る。

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 8

球座標系での基礎方程式系 (4)➢基礎方程式系➢格子分割と変数配置➢格子分割と座標系の設定➢球座標系での基礎方程式➢SIMPLE 系解法➢べき乗法➢1 次元定常移流–拡散方程式の解➢各種「上流化」法
➢Gershgorin の定理➢行列演算からみたCFL 条件➢Gershgorin の定理
による CFL➢格子界面での熱伝導率➢格子界面での剪断
粘性率
➢参考文献

各方向の運動方程式をあらわに書くと以下の通り。
”∂ “ 21τθθ
τφφ
∂
∂τφr
1r τrr +(sin θτθr ) +−−r 2 ∂rr sin θ ∂θr sin θ ∂φrr”
1∂∂τφr
11 ∂ “ 2r τrr +(sin θτθr ) +=r 2 ∂rr sin θ ∂θr sin θ ∂φﬀi
11 ∂ h 2∂−r (τθθ + τφφ ) −(τθθ + τφφ )2 r 2 ∂r
∂r

−br =

11”
∂ “ 211τθr
cot θτφφ∂∂τφθ
r τrθ +(sin θτθθ ) ++−r 2 ∂rr sin θ ∂θr sin θ ∂φrr”
∂τφθ
11
∂
1 ∂ “ 2r τrθ +(sin θτθθ ) +=r 2 ∂rr sin θ ∂θr sin θ ∂φ»–»–
11 ∂∂τθr
11∂τφφ∂2+(r τθr ) −−(sin θτφφ ) −2 r 2 ∂r
∂rr sin θ ∂θ
∂θ

−bθ =

2010 年夏の GFD セミナー1”∂ “ 211τφrcot θτφθ∂ `∂τφφ´r τrφ +sin θτθφ +++r 2 ∂rr sin θ ∂θr sin θ ∂φrr”“∂τφφ11∂ `1 ∂´
2sin θτθφ +r τrφ +=r 2 ∂rr sin θ ∂θr sin θ ∂φ
»
–
»–
1 ∂
111
∂τφr
∂τφθ
∂
2(r τφr ) −+(sin θτφθ ) −+2 r 2 ∂r
∂rr sin θ ∂θ∂θ−bφ =平成 22 年 8 月 21 日 – slide 9

SIMPLE 系解法 (Patankar, 1980)➢基礎方程式系➢格子分割と変数配置➢格子分割と座標系の設定➢球座標系での基礎方程式➢SIMPLE 系解法➢べき乗法➢1 次元定常移流–拡散方程式の解➢各種「上流化」法
➢Gershgorin の定理➢行列演算からみた
CFL 条件➢Gershgorin の定理
による CFL➢格子界面での熱伝導率➢格子界面での剪断
粘性率
➢参考文献

有限体積法ベースのものではまだまだ一般的。
0
 −G10= −G2−G3

D1
C11C21C31D2C12C22C32







0pD3


 C13   v1  +  b1 C23   v2   b2 b3v3
C33

1. 圧力場の推定値 p∗ を用意する。
2. p∗ に対応する速度場 v1∗ 、v2∗ 、v3∗ を以下より求める。
Ckk vk∗ = Gk p∗ − bk · · · (k = 1, 2, 3)′3. 
以下の方程式により、圧力の補正量pを求める。


−1
∗
D1 (diagC11 ) G1
D 1 v1
 +D2 (diagC22 )−1 G2  p′ = −  +D2 v2∗ +D3 v3∗+D3 (diagC33 )−1 G3

4. 速度の補正量 v1′ 、v2′ 、v3′ を以下より求める。
vk′ = (diagCkk )−1 Gk p′ (k = 1, 2, 3)
2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 10

べき乗法 (Patankar, 1980)➢基礎方程式系➢格子分割と変数配置➢格子分割と座標系の設定➢球座標系での基礎方程式➢SIMPLE 系解法➢べき乗法➢1 次元定常移流–拡散方程式の解➢各種「上流化」法
➢Gershgorin の定理➢行列演算からみたCFL 条件➢Gershgorin の定理
による CFL➢格子界面での熱伝導率➢格子界面での剪断粘性率➢参考文献

数値的安定性のための移流–拡散流束の「上流化」の一種❏ ある点でのある方向の熱流束 F を、その両側の温度 T0と T1 を用いた適切な内挿により求める❏ 各点での各方向の熱流束 F を、そこでの 1 次元定常熱流束と (なるべく) 一致させるように定める
advection velocity u
T=T0heat flux
F

T=T1

diffusion coefficient kx=0

2010 年夏の GFD セミナー

x=Ld0=dx

dTuT − κdx

境界条件:xT = T0 at x = 0
T = T1 at x = L

平成 22 年 8 月 21 日 – slide 11

1 次元定常移流–拡散方程式の解➢基礎方程式系➢格子分割と変数配置➢格子分割と座標系の設定➢球座標系での基礎方程式➢SIMPLE 系解法➢べき乗法➢1 次元定常移流–拡散方程式の解➢各種「上流化」法
➢Gershgorin の定理➢行列演算からみた
CFL 条件➢Gershgorin の定理
による CFL➢格子界面での熱伝導率➢格子界面での剪断
粘性率
➢参考文献

advection velocity uT1

T=T0

P=-10
P=-2

heat flux
F

P=0

T=T1T(x)P=2diffusion coefficient kx=0

x=Lx

P=10T00x

L0 ≤ x ≤ L での定常温度分布は
eP (x/L) − 1T (x) = T0 + (T1 − T0 ) P(右図)e −1uLκと表わされる。ただし、P ≡は cell Péclet 数、D ≡κ
Lはコンダクタンスと呼ばれる。
dT
これより熱流束 F ≡ uT − κはdx
PPPeP
e T0 − T1=D P
T0 − P
T1
F =u Pe −1
e −1e −1

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 12

各種「上流化」法 (1)➢基礎方程式系➢格子分割と変数配置➢格子分割と座標系の設定➢球座標系での基礎方程式➢SIMPLE 系解法➢べき乗法➢1 次元定常移流–拡散方程式の解➢各種「上流化」法
➢Gershgorin の定理➢行列演算からみた
CFL 条件➢Gershgorin の定理による CFL➢格子界面での熱伝導率➢格子界面での剪断
粘性率
➢参考文献

熱流束の表式として F = D[B(P ) T0 − A(P ) T1 ] を採用すると、各種「上流化」法は以下のようにまとめられる。
スキーム

A(P )

厳密解(指数法)

P
eP − 111− P
2
1 + max(0, −P )max[0, (1 − 0.5|P |)]+ max(0, −P )max[0, (1 − 0.1|P |)5 ]+ max(0, −P )中心差分

1 次の風上差分
ハイブリッド法
べき乗法

B(P )(= A(P ) + P )P ePeP − 1
11+ P
2
1 + max(0, P )(A(P ) + P )(A(P ) + P )ただし熱流束 F が温度 T0 と T1 の「内挿」で決まるようにするには、A(P ) 、B(P ) は共に非負であってほしい。
2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 13

各種「上流化」法 (2)熱流束の表式として F = D[B(P ) T0 − A(P ) T1 ] を採用すると、各種「上流化」法は以下のようにまとめられる。

A(P )

B(P )
151010B(P)

15A(P)➢基礎方程式系➢格子分割と変数配置➢格子分割と座標系の設定➢球座標系での基礎方程式➢SIMPLE 系解法➢べき乗法➢1 次元定常移流–拡散方程式の解➢各種「上流化」法
➢Gershgorin の定理➢行列演算からみた
CFL 条件➢Gershgorin の定理
による CFL➢格子界面での熱伝導率➢格子界面での剪断粘性率➢参考文献

5centered difference05

1st-order upwinding
exact / power-law-5

1st-order upwinding
0

centered difference

exact / power-law-5-15-10-505
Peclet number P

1015-15-10-505
Peclet number P

1015

ただし熱流束 F が温度 T0 と T1 の「内挿」で決まるようにするには、A(P ) 、B(P ) は共に非負であってほしい。
❏ 上流化をかけない中心差分では |P | < 2 でのみ OK❏ 上流化をかけることにより、大きな P でも (一応) OK

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 14

Gershgorin の定理➢基礎方程式系➢格子分割と変数配置➢格子分割と座標系の設定➢球座標系での基礎方程式➢SIMPLE 系解法➢べき乗法➢1 次元定常移流–拡散方程式の解➢各種「上流化」法
➢Gershgorin の定理➢行列演算からみたCFL 条件➢Gershgorin の定理
による CFL➢格子界面での熱伝導率➢格子界面での剪断粘性率➢参考文献

中心が aii 、半径が ri =n
X|aij | の円で囲まれた複素j=1,j6=i平面内の領域を Si とする。このとき行列 A = (aij ) のn[全ての固有値 λk は和集合
Si = S1 ∪ S2 ∪ · · · ∪ Sn のi=1内部に存在する。

証明: 行列 A の固有値 λk に対応する 固有ベクトルを x とする。x の成分のうちで絶対値が最大のものを xi とすると、全ての j 6= i に対して当然 |xi | ≥ |xj | を満たす。このとき、固有方程式 Ax = λk x の第 i 行に注目すると、

ai1 x1 + · · · + aii xi + · · · + ain xn = λk xi
これをさらに書き直すと、˛˛˛ X
˛ X˛˛ X˛˛˛˛xxjj
˛aij ˛ ≤|aii − λk | = ˛˛−|aij |aij ˛˛ ≤˛˛xxi˛i
˛ j6=ij6=ij6=iを得る。即ち、λk を含む領域 Si が必ず存在する。

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 15

行列演算からみた CFL 条件 (1)➢基礎方程式系➢格子分割と変数配置➢格子分割と座標系の設定➢球座標系での基礎方程式➢SIMPLE 系解法➢べき乗法➢1 次元定常移流–拡散方程式の解➢各種「上流化」法
➢Gershgorin の定理➢行列演算からみたCFL 条件➢Gershgorin の定理
による CFL➢格子界面での熱伝導率➢格子界面での剪断粘性率➢参考文献

係数一定の時間発展方程式 (拡散方程式とか) を離散化して得られる方程式は形式的に

xn+1 − xn= d − Cxn∆tと書ける。ただし、空間離散化を表わす行列 C は正則な疎行列とし、時間 t の方向には陽的な離散化を用いることにする。
これを書き直すと、

xn+1 − C −1 d = (I − ∆tC)(xn − C −1 d)を得る。即ち、このスキームでは定常解 C −1 d との差が反復の度に (I − ∆tC) 倍に「縮小」されていくことになる。

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 16

行列演算からみた CFL 条件 (2)➢基礎方程式系➢格子分割と変数配置➢格子分割と座標系の設定➢球座標系での基礎方程式➢SIMPLE 系解法➢べき乗法➢1 次元定常移流–拡散方程式の解➢各種「上流化」法
➢Gershgorin の定理➢行列演算からみたCFL 条件➢Gershgorin の定理
による CFL➢格子界面での熱伝導率➢格子界面での剪断
粘性率
➢参考文献

時間方向に陽的な離散化を行って得られたスキーム
xn+1 − C −1 d = (I − ∆tC)(xn − C −1 d)
このスキームが収束するための必要十分条件は、行列I − ∆tC のベキが 0 に収束すること、言い換えれば行列I − ∆tC のスペクトル半径 (絶対値最大の固有値) が 1 より小さいことである。
さらに振動解の発生を抑えるために、行列 I − ∆tC の固有値 (の実部) が非負でなければならない。
ちなみに、時間方向に完全に陰的に離散化した場合には

xn+1 − xn= d−Cxn+1 ,∆t2010 年夏の GFD セミナー

xn+1 −C −1 d = (I+∆tC)−1 (xn −C −1 d)

平成 22 年 8 月 21 日 – slide 17

Gershgorin の定理による CFL (1)➢基礎方程式系➢格子分割と変数配置➢格子分割と座標系の設定➢球座標系での基礎方程式➢SIMPLE 系解法➢べき乗法➢1 次元定常移流–拡散方程式の解➢各種「上流化」法
➢Gershgorin の定理➢行列演算からみた
CFL 条件➢Gershgorin の定理
による CFL➢格子界面での熱伝導率➢格子界面での剪断
粘性率
➢参考文献

反復の際に現われる実行列 I −∆tCの全ての固有値 λ が 0 < λ < 1 を満たすようにしたい。
そのために、Gershgorin の円が右図のような位置にあるようにしたい。

Im(λ)

0Re(λ)
1具体例として、3 次元熱輸送方程式を考える。
newT(i,j,k)− T(i,j,k)∆t

−c1 T(i,j,k−1) − c2 T(i,j−1,k) − c3 T(i−1,j,k)
+c4 T(i,j,k)= d(i,j,k) −−c5 T(i+1,j,k) − c6 T(i,j+1,k) − c7 T(i,j,k+1)

この場合、点 (i, j, k) に対応する Gershgorin の円は、中心 (1 − c4(i,j,k) ∆t, 0) の実軸上
半径 ∆t(c1 + c2 + c3 + c5 + c6 + c7 )(i,j,k)で特徴づけられる。

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 18

Gershgorin の定理による CFL (2)➢基礎方程式系➢格子分割と変数配置➢格子分割と座標系の設定➢球座標系での基礎方程式➢SIMPLE 系解法➢べき乗法➢1 次元定常移流–拡散方程式の解➢各種「上流化」法
➢Gershgorin の定理➢行列演算からみたCFL 条件➢Gershgorin の定理
による CFL➢格子界面での熱伝導率➢格子界面での剪断粘性率➢参考文献

Gershgorin の円が満たす条件はIm(λ)

1. 中心が実軸上 0 以上 1 未満
2. 左端が実軸上 0 以上

0Re(λ)
1

3. 右端が実軸上 1 未満
❏ 条件 1 より 0 < 1 − c4(i,j,k) ∆t < 1❏ 条件 2 より0 < (1 − c4(i,j,k) ∆t) − ∆t(c1 + c2 + c3 + c5 + c6 + c7 )(i,j,k)
このうち、より厳しい条件をとると1を得る。
(0 <)∆t <(c1 + c2 + c3 + c4 + c5 + c6 + c7 )(i,j,k)
ただし、Gershgorin の定理が示す固有値の存在範囲はかなり広いので、ここで求まった∆t の閾値も必要以上に厳しくなっていることに注意2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 19

格子界面での熱伝導率の設定 (Patankar, 1980)➢基礎方程式系➢格子分割と変数配置➢格子分割と座標系の設定➢球座標系での基礎方程式➢SIMPLE 系解法➢べき乗法➢1 次元定常移流–拡散方程式の解➢各種「上流化」法
➢Gershgorin の定理➢行列演算からみた
CFL 条件➢Gershgorin の定理
による CFL➢格子界面での熱伝導率➢格子界面での剪断粘性率➢参考文献

T(i-1)k(I) ?

k(i-1) >> k(i)xp(i-1)

xc(I)

T(i)
xp(i)

x

体積要素の境界面 xc(I) での熱流束が連続になるように、この点での熱伝導率 k(I) を定める。
ただし、各体積要素の内部では k(i) は一定であるとする。

x = xc (I) での温度を Tc とおく。熱流束の連続性より、T(i) − Tc
Tc − T(i−1)
T(i) − T(i−1)= k(i) 1= k(i−1) 1
k(I)
δx(I)
dx(i−1)dx(i)22
これを満たす k(I) は2δx(I)
dx(i−1) dx(i)=+k(I)k(i−1)k(i)により与えられる。特に格子が等間隔 (dx(i−1) = dx(i) ) の2k(i−1) k(i)の如き調和平均となる。
場合には、k(I) =k(i−1) + k(i)2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 20

格子界面での剪断粘性率の設定 (Ogawa et al., 1991)➢基礎方程式系➢格子分割と変数配置➢格子分割と座標系の設定➢球座標系での基礎方程式➢SIMPLE 系解法➢べき乗法➢1 次元定常移流–拡散方程式の解➢各種「上流化」法
➢Gershgorin の定理➢行列演算からみた
CFL 条件➢Gershgorin の定理
による CFL➢格子界面での熱伝導率➢格子界面での剪断粘性率➢参考文献

熱伝導率の定義と同様、格子界面での剪断応力の連続性(に相当するもの) を考える。
η2(I,j)y j
界面での熱伝導率の定義に倣い、
dy(j−1) dy(j)2δy(J)
ηyx=+J
η1(i,J)η(i,j−1)η(i,j)
ηxyη
(i,J)
η1(i-1,J)
12δx(I)dx(i−1) dx(i)


=+
j-1
η2(I,j)η(i−1,j)η(i,j)
η2(I,j-1)
で与えられる η1(i,J) 及び η2(I,j)
ii-1
Ixを考える。
x 方向に作用する剪断応力 τyx(I,J) の計算に用いる粘性は
η1(i−1,J) dx(i−1) + η1(i,J) dx(i)ηyx(I,J) ≡2δx(I)y 方向に作用する剪断応力 τxy(I,J) の計算に用いる粘性は
η2(I,j−1) dy(j−1) + η2(I,j) dx(j)ηxy(I,J) ≡2δy(J)2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 21

参考文献➢基礎方程式系➢格子分割と変数配置➢格子分割と座標系の設定➢球座標系での基礎方程式➢SIMPLE 系解法➢べき乗法➢1 次元定常移流–拡散方程式の解➢各種「上流化」法
➢Gershgorin の定理➢行列演算からみたCFL 条件➢Gershgorin の定理
による CFL➢格子界面での熱伝導率➢格子界面での剪断
粘性率
➢参考文献

❏ M. Kameyama, A. Kageyama, and T. Sato. Multigrid

iterative algorithm using pseudo-compressibility forthree-dimensional mantle convection with strongly variableviscosity. J. Comput. Phys., 206, 162-181, 2005.
❏ M. Kameyama. ACuTEMan: A multigrid-based mantleconvection simulation code and its optimization to theEarth Simulator. J. Earth Simulator, 4, 2-10, 2005.
❏ M. Ogawa, G. Schubert, and A. Zebib. Numericalsimulations of three-dimensional thermal convection in a
fluid with strongly temperature-dependent viscosity.
J. Fluid Mech., 233, 299-328, 1991.
❏ S. V. Patankar. Numerical Heat Transfer and Fluid Flow.

Hemisphere, Washington DC, 197 pp., 1980.
(邦訳「コンピュータによる熱移動と流れの数値解析」
森北出版)2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 22