マントル対流の「数値」「流体」「力学」
亀山 真典 @ 愛媛大学地球深部ダイナミクス研究センター
平成 22 年 8 月 21 日

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 1

「固体」地球内部の「流体」現象➢地球内部流体現象➢地震学的層構造➢温度・圧力状態➢上昇プリューム
➢滞留スラブ➢アイソスタシー➢後氷期回復➢マントルの粘性
マントル対流の数値「流体」「力学」「固体」地球内部の「流体」現象

マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3
マントル対流の「数値」流体「力学」

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 2

「固体」地球内部の「流体」現象の例「固体」地球内部の「流体」現象➢地球内部流体現象➢地震学的層構造➢温度・圧力状態➢上昇プリューム
➢滞留スラブ➢アイソスタシー➢後氷期回復➢マントルの粘性

6370km6346kmマントル(岩石)

マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2

マントル内
地殻 の対流

下図はKageyama, A. et al., (2004)より

温度の水平平均からのずれ

3480km外核(液体鉄)「マグマの海」の中での金属鉄の分離1220kmマントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3
マントル対流の「数値」流体「力学」

2010 年夏の GFD セミナー

内核
(固体鉄)
0

右図は

外核内の
対流

Ichikawa, H., Labrosse, S.,
Kurita, K., (2010) より

渦度の自転軸方向成分
平成 22 年 8 月 21 日 – slide 3

地震学でみた地球の層構造

マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3

地震波の伝わり方の違いから、マントルも深さによって 地震波伝播速度 [km/s] または密度 [g/cc]
0 2 4 6 8 10 12 140
区別されている。
地殻❏ 地表面から深さ約 410km〜660km の部分を「マントル遷移層」と呼ぶ。
❏ 深さ約 660km より

下を「下部マントル」、上を「上部マントル」、と呼ぶ。
❏ 最下部の約 200km の部分は

1000

2000

マントル対流の「数値」流体「力学」

2010 年夏の GFD セミナー

P波
(縦波)

マントル
(岩石)

3000

4000

外核
(流体鉄)

5000

「D” 層」と呼ばれる。(ディーツーダッシュまたはディーダブルプライム)

S波
(横波)
密度

深さ [km]「固体」地球内部の「流体」現象➢地球内部流体現象➢地震学的層構造➢温度・圧力状態➢上昇プリューム
➢滞留スラブ➢アイソスタシー➢後氷期回復➢マントルの粘性

6000

内核
(固体鉄)

右図は Dziewonski and Anderson (1981) のPREM (Preliminary Reference Earth Model) より
平成 22 年 8 月 21 日 – slide 4

地球内部の温度・圧力状態

マントル対流の数値「流体」「力学」

地球内部の圧力状態は、ほぼ正確に求めることができる❏ 深さ 660km で約 23 万気圧(ここより下が下部マントル)

マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3
マントル対流の「数値」流体「力学」

温度 [K]
2000

4000

6000

200

400

600

1000

2000

❏ マントルの底で約 135 万気圧❏ 地球中心で約 360 万気圧

3000

4000

マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1

00

深さ [km]「固体」地球内部の「流体」現象➢地球内部流体現象➢地震学的層構造➢温度・圧力状態➢上昇プリューム
➢滞留スラブ➢アイソスタシー➢後氷期回復➢マントルの粘性

地球内部の温度状態の正確な見積りは難しいのだが、

5000

6000

❏ 深さ 660km で約 1900K❏ マントルの底で約 3000〜4000K❏ 地球中心で約 5000〜6000K

2010 年夏の GFD セミナー0

圧力 [万気圧]唐戸 (2000) による地
球内部の温度・圧力分
布の推定値
平成 22 年 8 月 21 日 – slide 5

地震波トモグラフィーでみた地球の内部 (1)「固体」地球内部の「流体」現象➢地球内部流体現象➢地震学的層構造➢温度・圧力状態➢上昇プリューム
➢滞留スラブ➢アイソスタシー➢後氷期回復➢マントルの粘性

地震の波を使ってやると、地球内部を「CT スキャン」で撮影することができる。
タヒチ〜 ハワイ

マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3
マントル対流の「数値」流体「力学」

地震波 (P 波) の伝播速度異常 (平均からのずれ) の分布 (Zhao, 2004)
その深さでの平均伝播速度からの微小 (∼ ±%) なずれに注目
青 (高速度異常) → 固い → 低温
基本的な解釈は · · ·赤 (低速度異常) → 柔らかい → 高温
タヒチ島やハワイ諸島は「ホットスポット」と呼ばれ、これらを作った火山活動の源の位置はマントルに固定されていると思われている。

特にタヒチ島の地下には、下部マントルに根を持つ高温の上昇流があると解釈できる。

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 6

地震波トモグラフィーでみた地球の内部 (2)「固体」地球内部の「流体」現象➢地球内部流体現象➢地震学的層構造➢温度・圧力状態➢上昇プリューム
➢滞留スラブ➢アイソスタシー➢後氷期回復➢マントルの粘性

マントル対流の下降流にあたる、沈み込んだプレート(「スラブ」と呼ばれる) はどう見えるか?

中国北西部―中部日本

マントル対流の数値「流体」「力学」

左図はZhao, D., (2004) より

マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3
マントル対流の「数値」流体「力学」

日本列島下に沈み込んだ太平洋プレートを想像させる高速度異常が、地表面から深部に向かって延びている。
❏ 太平洋スラブは深さ 660km (上部マントルと下部マントルの境界) 周辺で「横たわっている」ように見える「滞留スラブ」または「スタグナントスラブ」

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 7

「固体」のマントルが「流れる」例 (1)「固体」地球内部の「流体」現象➢地球内部流体現象➢地震学的層構造➢温度・圧力状態➢上昇プリューム
➢滞留スラブ➢アイソスタシー➢後氷期回復➢マントルの粘性
マントル対流の数値「流体」「力学」

高校の地学で習った例は「アイソスタシー」(地殻均衡説)❏ 地殻の厚さと標高の間には

標高の「高い」ところでは地殻が「厚い」
標高の「低い」ところでは地殻が「薄い」という関係がある。

❏ 「アルキメデスの原理」により、軽い地殻がマントルに浮かんでいると考えるとうまく説明できる

マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2

ρ=ρw
ρ=ρc

H

マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3

アイソスタシーのしくみ

h

ρ=ρm

マントル対流の「数値」流体「力学」

2010 年夏の GFD セミナーb

ρmb = ρc(b+h)

海面に浮かぶ氷山のように、密度の小さい地殻は、密度の大きいマント
ルの浮力により支えられて、標高の高いところほどマントルの中に大きく食い込んでいる。
平成 22 年 8 月 21 日 – slide 8

「固体」のマントルが「流れる」例 (2)「固体」地球内部の「流体」現象➢地球内部流体現象➢地震学的層構造➢温度・圧力状態➢上昇プリューム
➢滞留スラブ➢アイソスタシー➢後氷期回復➢マントルの粘性
マントル対流の数値「流体」「力学」

後氷期回復 (Postglacial Rebound)
スカンジナビア半島では、最終氷期に陸地をおおっていた
厚さ数千 m の氷河がとけてなくなった後も、数千年にわたって地表面の隆起が続いている。
❏ 氷河がとけて崩れたアイソスタ
シーを回復するために、地殻が浮き上がっている、と考えれば説明できる。

マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3
マントル対流の「数値」流体「力学」

2010 年夏の GFD セミナー

❏ マントルがとても「ねばねば」

しているために、地殻が完全に浮き上がるまでに時間がかかっている。
図は BIFROST (Baseline Inferences from
Fennoscandian
Rebound
Observations,
Sealevel and Tectonics) Project よ り(Scherneck et al., 2001)
平成 22 年 8 月 21 日 – slide 9

「固体」のマントルが「流れる」例 (3)「固体」地球内部の「流体」現象➢地球内部流体現象➢地震学的層構造➢温度・圧力状態➢上昇プリューム
➢滞留スラブ➢アイソスタシー➢後氷期回復➢マントルの粘性

後氷期回復に伴う地表面の隆起の速度から、マントルの粘性率を見積ることができる。
粘性率 η の非圧縮性流体からなる半無限体のマントルに、荷重によって表面形状に微小で周期的な変化が加えられたとする。

マントル対流の数値「流体」「力学」

w

マントル対流の「数値」「流体」力学 1λマントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3
マントル対流の「数値」流体「力学」

表面形状の変化が内部の流れによって緩和される時定4πη
数は τr =とかける。
ρgλ

典型的な後氷期回復のデータを用いると、

ρgλτr
η=
≃ 1021 [Pa s]
4π
3300[kg/m3] × 9.8[m/s2 ] × 3000[km] × 4400[yr]=4π2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 10

「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」➢物性値➢基礎方程式系➢状態方程式簡略化➢無次元化➢無次元方程式➢非弾性流体近似➢よくある近似
➢Ra
➢Re
➢Ek
➢P r

マントル対流の数値「流体」「力学」

マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3
マントル対流の「数2010 年夏の GFD セミナー
値」流体「力学」

平成 22 年 8 月 21 日 – slide 11

外核・マントルの「流体」的な物性の比較「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」➢物性値➢基礎方程式系➢状態方程式簡略化➢無次元化➢無次元方程式➢非弾性流体近似➢よくある近似
➢Ra
➢Re
➢Ek
➢P r
マントル対流の「数値」「流体」力学 1

記号

マントル

α
κ
Cp
ρ
ν
σ
V
L

マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3

意味
熱膨張率
熱拡散率
定圧比熱
密度

動粘性率

η
=
ρ

電気伝導度
典型的流速

典型的長さ
(対流層の厚さ)



値
外核

10−5
10−6
103
3.3 ∼ 5.6
9.9 ∼ 12.23

単位
K−1
m2 /s
J/kg K
×103 kg/m3

1016 ∼ 1020

10−6

m2 /s

10−2 ∼ 101
3 × 10−9
10−1

4 × 105
4 × 10−4
104

S/m
m/s
m/yr

2.9 × 106

2.3 × 106m

マントルの物性は本多 (1997)、外核の物性は Kono and Roberts (2002) より

動粘性率 ν 、電気伝導度 σ が桁外れに異なる。

マントル対流の「数2010 年夏の GFD セミナー
値」流体「力学」

平成 22 年 8 月 21 日 – slide 12

基礎方程式系 (1 of 3)「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」➢物性値➢基礎方程式系➢状態方程式簡略化➢無次元化➢無次元方程式➢非弾性流体近似➢よくある近似
➢Ra
➢Re
➢Ek
➢P r

回転系での流体の熱対流を記述する方程式系から、マントルの流れを扱う方程式を導出しよう❏ 質量保存則∂ρ+ ∇ · (ρv) = 0∂t❏ 運動量保存則 (Navier-Stokes 方程式の手前)

マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3
マントル対流の「数2010 年夏の GFD セミナー
値」流体「力学」

粘性抵抗

重力
z
}|
{
z}|{
Dv
ρ
= − ∇p + ∇ · τ + ρg
Dt

回転する座標系でみたときの「みかけの力」}|{
z+ 2ρv × Ω − ρΩ × Ω × x{z}| {z } |
コリオリ力

遠心力

平成 22 年 8 月 21 日 – slide 13

基礎方程式系 (2 of 3)「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」➢物性値➢基礎方程式系➢状態方程式簡略化➢無次元化➢無次元方程式➢非弾性流体近似➢よくある近似
➢Ra
➢Re
➢Ek
➢P r

回転系での流体の熱対流を記述する方程式系から、マントルの流れを扱う方程式を導出しよう❏ 構成方程式

1「歪速度テンソル」
ε̇ = (∇ ⊗ v + v ⊗ ∇)「変形速度テンソル」
2
1「スピンテンソル」Ω = (∇ ⊗ v − v ⊗ ∇)「回転速度テンソル」
2
1
τ = 2η ε̇ − tr(ε̇)I「( 偏差応力テンソル」)
3

マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2❏ 状態方程式

マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3
マントル対流の「数2010 年夏の GFD セミナー
値」流体「力学」

ρ = ρ(T, p, · · · )平成 22 年 8 月 21 日 – slide 14

基礎方程式系 (3 of 3)「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」➢物性値➢基礎方程式系➢状態方程式簡略化➢無次元化➢無次元方程式➢非弾性流体近似➢よくある近似
➢Ra
➢Re
➢Ek
➢P r
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1

回転系での流体の熱対流を記述する方程式系から、マントルの流れを扱う方程式を導出しよう❏ 熱輸送方程式

DT
Dp
ρCp
− αT
= ∇ · (k∇T ) + ρH + |{z}
Φ
|{z}
Dt
Dt
内部発熱

粘性散逸

このうち、粘性散逸項は以下で与えられる。


1 2
∂vi
2
Φ=
τji = ε̇ : τ = 2η ε̇II − ε̇I∂xj
3

マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3
マントル対流の「数2010 年夏の GFD セミナー
値」流体「力学」

ε̇I = tr(ε̇) = ε̇11 + ε̇22 + ε̇33qε̇II = ε̇211 + ε̇222 + ε̇233 + 2ε̇212 + 2ε̇223 + 2ε̇213

平成 22 年 8 月 21 日 – slide 15

マントル対流向け状態方程式の簡略化

マントル対流の数値「流体」「力学」➢物性値➢基礎方程式系➢状態方程式簡略化➢無次元化➢無次元方程式➢非弾性流体近似➢よくある近似
➢Ra
➢Re
➢Ek
➢P r
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2

マントル内の密度変化 (上面と下地震波伝播速度 [km/s] または密度 [g/cc]
0 2 4 6 8 10 12 14
面で約 65%) のほとんどは、静水
0
地殻
圧下の圧縮に起因している。
そこで、分布の変化を❏ 静水圧状態の「基本場」❏ 流れに関係する「擾乱」に分けて考える。
′

ρ(T, p)= ρ + ρ= ρ[1 + χp′ − αT ′ ]ρ= ρ(T , p)
∇p= ρgマントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3
マントル対流の「数2010 年夏の GFD セミナー
値」流体「力学」

1000

S波
(横波)
密度

2000

深さ [km]「固体」地球内部の「流体」現象

P波
(縦波)

マントル
(岩石)

3000

4000

外核
(流体鉄)

5000

6000

内核
(固体鉄)

右図は Dziewonski and Anderson (1981) のPREM (Preliminary Reference Earth Model) より
平成 22 年 8 月 21 日 – slide 16

基礎方程式系の無次元化「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」➢物性値➢基礎方程式系➢状態方程式簡略化➢無次元化➢無次元方程式➢非弾性流体近似➢よくある近似
➢Ra
➢Re
➢Ek
➢P r

スケール
長さ
時間
速度
圧力
温度

定義
L
L2
κ
κ
L
ηκ
L2
∆T

意味
対流層の厚さ
熱拡散時間
熱拡散速度
粘性応力
対流を駆動する特徴的な温度差

マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3
マントル対流の「数2010 年夏の GFD セミナー
値」流体「力学」

平成 22 年 8 月 21 日 – slide 17

基礎方程式系 (無次元) 1 of 2
「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」➢物性値➢基礎方程式系➢状態方程式簡略化➢無次元化➢無次元方程式➢非弾性流体近似➢よくある近似
➢Ra
➢Re
➢Ek
➢P r
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3

以下特に断らない限り、変数は全て無次元化された量を示す。また、「擾乱」を表す ’ は省略する。
❏ 連続の式∇ · (ρv) = 0
ただし非弾性流体近似を用いて ∂ρ/∂t = 0 とし、マントル内を伝わる音波 (地震波; v = O(km/s)) を除去した。

なお静水圧による断熱的密度変化のスケールハイト h は

−1
1 dρ
1
h=
=
≃ O(106 )[m]ρ dz
χa ρ gであり、マントルの厚さ (2900km) と同程度。

マントル対流の「数2010 年夏の GFD セミナー
値」流体「力学」

平成 22 年 8 月 21 日 – slide 18

基礎方程式系 (無次元) 2 of 2
「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」➢物性値➢基礎方程式系➢状態方程式簡略化➢無次元化➢無次元方程式➢非弾性流体近似➢よくある近似
➢Ra
➢Re
➢Ek
➢P r

❏ 運動方程式 (遠心力項は省略)"∇⊗v+v⊗∇
2− (∇ · v)I
3
1− Ra ρ α T g ĝ +
ρ v × Ω̂| {z } Ek| {z }熱的浮力

1 Dv
ρ= − ∇p + ∇ · η
P r Dt!#

コリオリ力❏ 熱輸送方程式

マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3
マントル対流の「数2010 年夏の GFD セミナー
値」流体「力学」



DT
ρCp
= ∇ · k∇ T + T + Diρ g α T (v · ĝ)
|
{z
}
Dt
|
{z
}
熱伝導

Di
H
+
Φ + |{z}{z }|Ra

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

内部発熱

粘性散逸

平成 22 年 8 月 21 日 – slide 19

非弾性流体近似 (Anelastic Liquid Approximation)「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」➢物性値➢基礎方程式系➢状態方程式簡略化➢無次元化➢無次元方程式➢非弾性流体近似➢よくある近似
➢Ra
➢Re
➢Ek
➢P r
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1

流体の圧縮性の効果を (ある程度) 正しく取り入れる(King, 6 authors, and Kameyama, 2010)❏ 流体の静的な圧縮性を考慮した質

量保存則 (∇ · (ρv) = 0)
(ただし、動圧による密度変化を無視する「打ち切り版」のほうが一般的)❏ 力学的な仕事から熱的なエネル

Di = 0ギーへの変換✏ 静水圧下の断熱圧縮に伴う温度変化 (断熱温度勾配)

Di = 0.25✏ 粘性散逸 (摩擦) による加熱

マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3

なお、圧縮性の効果の強さはαgd
散逸数 Di =ではかられる。
Cp

マントル対流の「数2010 年夏の GFD セミナー
値」流体「力学」

Di = 0.5
平成 22 年 8 月 21 日 – slide 20

マントル対流研究でよく用いられる近似「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」➢物性値➢基礎方程式系➢状態方程式簡略化➢無次元化➢無次元方程式➢非弾性流体近似➢よくある近似
➢Ra
➢Re
➢Ek
➢P r

❏ 非弾性流体近似 (Di 6= 0)

静的な圧縮による密度変化の影響を残す。
❏ ブシネスク近似 (Di = 0)運動方程式の浮力項を除いて、密度変化の影響を無視する。
❏ 拡張ブシネスク近似 (extended Boussinesq Approximation)

熱輸送方程式は非弾性流体近似のもの
運動方程式+連続の式はブシネスク近似のもの

マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3

世の中のマントル対流シミュレーションはだいたい、ブシ
ネスク近似あるいは拡張ブシネスク近似を使っている。
非弾性近似に基づくシミュレーションは現在でもあまり盛んではない。その原因の 1 つは、地球内部物質の状態方程式の不確定さがまだ大きいこと。

マントル対流の「数2010 年夏の GFD セミナー
値」流体「力学」

平成 22 年 8 月 21 日 – slide 21

外核・マントルの流体としての性質 1「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」➢物性値➢基礎方程式系➢状態方程式簡略化➢無次元化➢無次元方程式➢非弾性流体近似➢よくある近似
➢Ra
➢Re
➢Ek
➢P r
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2

Rayleigh 数
α∆T gL3
熱的な浮力の大きさ
Ra ≡=νκ
粘性抵抗の大きさマントルでは · · ·ν = 1018 m2 /s∆T = 4000Kとしたら、Ra ∼ 107

外核では · · ·ν = 10−6 m2 /s∆T = 2000Kとしたら、Ra ∼ 2 × 1030

❏ マントルでも外核でも臨界 Ra よりも十分大きいので、熱対流は起こっている

マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3
マントル対流の「数2010 年夏の GFD セミナー
値」流体「力学」

平成 22 年 8 月 21 日 – slide 22

外核・マントルの流体としての性質 2「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」➢物性値➢基礎方程式系➢状態方程式簡略化➢無次元化➢無次元方程式➢非弾性流体近似➢よくある近似
➢Ra
➢Re
➢Ek
➢P r
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2

Reynolds 数
LV
代表的長さ × 流れの速さ
Re ≡
=
ν
動粘性率
マントルでは · · ·ν = 1018 m2 /s
V = 3 × 10−9 m/s
L = 2.9 × 106 mとしたら、Re ∼ 10−20

外核では · · ·ν = 10−6 m2 /s
V = 4 × 10−4 m/s
L = 2.3 × 106 mとしたら、Ra ∼ 109❏ マントル内は Re が非常に小さい「遅い流れ」のため、

速度の非線型項は無視できる。
❏ 外核内は激しい乱流になっている。

マントル対流の「数値」流体力学 3
マントル対流の「数2010 年夏の GFD セミナー
値」流体「力学」

平成 22 年 8 月 21 日 – slide 23

外核・マントルの流体としての性質 3「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」➢物性値➢基礎方程式系➢状態方程式簡略化➢無次元化➢無次元方程式➢非弾性流体近似➢よくある近似
➢Ra
➢Re
➢Ek
➢P r
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2

Ekman 数ν
粘性による散逸
Ek ≡=2
2L |Ω|
コリオリ力による加速マントルでは · · ·ν = 1018 m2 /s
L = 2.9 × 106 mとしたら、Ek ∼ 2 × 1015

外核では · · ·ν = 10−6 m2 /s
L = 2.3 × 106 mとしたら、Ek ∼ 10−15

地球の自転によるコリオリ力による影響を❏ マントルは全く感じないで流れている❏ 外核は強く感じながら流れている(地球の磁場が南北方向を向くのもこのせい)

マントル対流の「数値」流体力学 3
マントル対流の「数2010 年夏の GFD セミナー
値」流体「力学」

上図はKageyama and Sato (1995) より

平成 22 年 8 月 21 日 – slide 24

外核・マントルの流体としての性質 4「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」➢物性値➢基礎方程式系➢状態方程式簡略化➢無次元化➢無次元方程式➢非弾性流体近似➢よくある近似
➢Ra
➢Re
➢Ek
➢P r
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3

Prandtl 数
ν
運動量の拡散の速さ
Pr ≡ =
κ
熱の拡散の速さマントルでは · · ·ν = 1018 m2 /s
κ = 10−6 m2 /sとしたら、P r ∼ 1024

外核では · · ·ν = 10−6 m2 /s
κ = 10−6 m2 /sとしたら、P r ∼ 100

マントルの流れの異常さをよく特徴づける量の 1 つ❏ 運動量の拡散に比べて熱の拡散が圧倒的に遅い(運動量の輸送は瞬時に起こるとみなしてよい)❏ 熱境界層に比べて速度境界層が圧倒的に厚い(マントル全体が速度境界層)

マントル対流の「数2010 年夏の GFD セミナー
値」流体「力学」

平成 22 年 8 月 21 日 – slide 25

「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1➢数値解法骨格➢数値解法難点➢特異な流動特性➢外核対流の困難
マントル対流の「数値」流体力学 1

マントル対流の「数値」「流体」力学 1

マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3
マントル対流の「数値」流体「力学」

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 26

マントル対流シミュレーションの手順「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1➢数値解法骨格➢数値解法難点➢特異な流動特性➢外核対流の困難
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3
マントル対流の「数値」流体「力学」

最も簡単な、非圧縮性流体 + ブシネスク近似 (Di = 0) のもとでの基礎方程式を例に用いる。
do 時間発展ループ
solve 熱輸送方程式 (温度場 T を更新)

DT= ∇ · (κ∇T ) + q
Dtsolve 運動方程式+連続の式 (速度場 v と圧力場 p を更新)
0=∇·v
1 Dv
= −∇p + ∇ · [η(∇ ⊗ v + v ⊗ ∇)] + Ra T ez
P r Dtend do 時間発展ループ

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 27

マントル対流シミュレーションの困難「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1➢数値解法骨格➢数値解法難点➢特異な流動特性➢外核対流の困難
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3
マントル対流の「数値」流体「力学」

各時刻での流れ場を解くのが非常に厄介 (時間がかかる)

0 = ∇·v
1 Dv
0≃= −∇p + ∇ · [η(∇ ⊗ v + v ⊗ ∇)] + Ra T ez
P r Dt
主にマントル物質の特異な流動特性に起因する。
❏ 圧力 p を直接求める式がない (非圧縮性流体の宿命)❏ 粘性率 η が非常に大きい (P r ∼ O(1024 ) → ∞)✏ 慣性項が無視できるほど小さい✏ 浮力と粘性抵抗の 釣り合い (力学的平衡) が成り立つ流

れ場を各時間ステップで求める 必要がある❏ 粘性率 η の空間変化が非常に大きい✏ 速度場 v と圧力場 p を解く式の性質が非常に悪い
2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 28

マントルの特異な流動特性

マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1➢数値解法骨格➢数値解法難点➢特異な流動特性➢外核対流の困難
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3

マントル物質の流動則は一般的に以下の表式で書かれる

∗∗
E + pVn−mε̇ ∝ |{z}
τ
d
exp −
|{z}
|{z}
RT
歪速度
応力依存性 結晶粒径依存性 |{z}
温度・圧力依存性

1024

拡散クリープマントル物質の粘性率の典型的な値はおおよそ 1021 Pa s であるが、変形させる条件 (温度・圧力・変形速度・変形機構など ) に
よって大きく異なる。

1022

粘性率 [Pa s]「固体」地球内部の「流体」現象20

10 マントル対流の典型的な条件1018

マントル対流の「数値」流体「力学」

10

16

τ = 105 Pa
τ = 106 Pa

転位クリープ1000 温度 [K]2010 年夏の GFD セミナー

2000

特に温度依存性に注目してみれば、およそ 100K の温度変化で粘性率は 1 桁低下する。
上部マントル (地表から深さ約 660km まで) 条件での粘性率の推定値。図は Schubert et al., 2001 より
平成 22 年 8 月 21 日 – slide 29

外核内対流シミュレーションの困難「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1➢数値解法骨格➢数値解法難点➢特異な流動特性➢外核対流の困難
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2

主として粘性率が低いことに起因する❏ 電磁流体としての扱いが必要❏ 回転の効果が非常に強い

Ek ∼ O(10−15 )❏ Reynolds 数が非常に大きいRe ∼ O(109 )❏ Rayleigh 数が非常に大きいRa ∼ O(1030 )マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3
マントル対流の「数値」流体「力学」

現実的なパラメータ値を用いた
外核内対流シミュレーションは到底不可能
最近の外核内対流シミュレーションで採用された Ek と Ra と、現実に予想される値との比較
図は Kono and Roberts (2002) より (ここでは粘性を数桁大きめに見積っている)

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 30

「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1➢ここまでのまとめ➢流線関数➢原始変数解法➢反復解法➢ACuTE その 1
➢擬似圧縮性法とは?
➢ACuTE その 2
➢反復解法 2
➢Gershgorin の定理➢反復解法 3➢多重格子法➢多重格子法並列化➢箱型計算例

マントル対流の「数値」流体力学 1

マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2
マントル対流の「数2010 年夏のGFD セミナー
3値」流体力学

平成 22 年 8 月 21 日 – slide 31ここまでのまとめ「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1➢ここまでのまとめ➢流線関数➢原始変数解法➢反復解法➢ACuTE その 1
➢擬似圧縮性法とは?
➢ACuTE その 2
➢反復解法 2
➢Gershgorin の定理➢反復解法 3➢多重格子法➢多重格子法並列化➢箱型計算例
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2❏ 流体力学的にみれば、マントルの流れは極めて特異✏ 粘性率が非常に高い流体 のゆっくりとした対流✏ 粘性率が大きく変化 する流体

どちらかといえば、外核内の対流のほうが、「普通の」流体現象に近い(Prandtl 数 P r で比べれば明らか)❏ マントル対流の数値シミュレーションでは、計算時間の9 割以上が流れ場の求解に消費されている✏ (ほぼ) 非圧縮性流体の定常流れを頻繁に求める❏ マントル対流の流れ場の数値解法には、極めて特異な手法が求められる✏ メジャーな他の数値流体力学的手法をそのまま利用することができない

マントル対流の「数2010 年夏のGFD セミナー
3値」流体力学

平成 22 年 8 月 21 日 – slide 32

2 次元マントル対流問題でよくある方法「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1➢ここまでのまとめ➢流線関数➢原始変数解法➢反復解法➢ACuTE その 1
➢擬似圧縮性法とは?
➢ACuTE その 2
➢反復解法 2
➢Gershgorin の定理➢反復解法 3➢多重格子法➢多重格子法並列化➢箱型計算例
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2

流線関数 Ψ を用いるのが最も簡単
非圧縮 (∇ · v = 0) の流れ場 v は、あるベクトルポテン
シャル Ψ = (Ψx , Ψy , Ψz ) を用いて v = ∇ × Ψ と書ける。
特に 2 次元問題 (∂/∂y = 0) の場合には、∂Ψz ∂Ψy
∂Ψy
vx =
−
=−
,
∂y
∂z
∂z

∂Ψy ∂Ψx
∂Ψy
vz =
−
=
∂x
∂y
∂x

となる。さらに圧力 p も消去すれば、結果的に運動方程式は Ψy のみを含んだ 4 階の偏微分方程式に変形される。
(ただし 3 次元ではそうはいかない)Ψy の等値線分布により流れ場は一目瞭然❏ 接線方向がその点での流速ベクトルの向き❏ 等値線の間隔が狭いほど流れが速い

マントル対流の「数2010 年夏のGFD セミナー
3値」流体力学

平成 22 年 8 月 21 日 – slide 33

原始変数による流れ場の解法: 原理「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1➢ここまでのまとめ➢流線関数➢原始変数解法➢反復解法➢ACuTE その 1
➢擬似圧縮性法とは?
➢ACuTE その 2
➢反復解法 2
➢Gershgorin の定理➢反復解法 3➢多重格子法➢多重格子法並列化➢箱型計算例
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2

0
0

=
=

∇·v
−∇p + ∇ · [η(∇ ⊗ v + v ⊗ ∇)] + Ra T er

速度 v と圧力 p をしかるべく空間離散化してやると、こ
の式は結局


 

0
D1 D2 D3
p0
 −G1 C11 C12 C13   v1   0 
+


0=
 −G2 C21 C22 C23   v2   0 
v3
RaT−G3 C31 C32 C33の如き連立一次方程式に書き直せる。
マントル対流の数値シミュレーションではこの連立一次方
程式を各時間ステップごとに解く必要がある。
言い換えれば、この (大規模) 連立一次方程式を高速に解けるようにすることが、マントル対流シミュレーションの高速化にとって極めて重要である。

マントル対流の「数2010 年夏のGFD セミナー
3値」流体力学

平成 22 年 8 月 21 日 – slide 34

原始変数による流れ場の解法: 難点「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1➢ここまでのまとめ➢流線関数➢原始変数解法➢反復解法➢ACuTE その 1
➢擬似圧縮性法とは?
➢ACuTE その 2
➢反復解法 2
➢Gershgorin の定理➢反復解法 3➢多重格子法➢多重格子法並列化➢箱型計算例
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2



0
 −G1
0=
 −G2
−G3

D1
C11
C21
C31

D2
C12
C22
C32









D3
p
0
 v1   0 
C13 

+

C23   v2   0 
C33
v3
RaT

この連立一次方程式を数値的に解くのも、けっこう大変❏ ベクトル変数 (速度場の 3 成分) を解く必要がある✏ 巷間広く使われている数値計算ライブラリの求解ルーチンが使えない (スカラー変数の求解が前提)❏ 係数行列の対角成分に 0 がある (圧力 p の離散化方程式に由来)✏ メジャーな連立一次方程式の数値解法をそのまま使うことができない (対角成分の逆数を使おうとするから)❏ 係数行列は悪条件かつ特異 (一意な解が存在しない)

マントル対流の「数2010 年夏のGFD セミナー
3値」流体力学

平成 22 年 8 月 21 日 – slide 35

連立一次方程式の反復解法: 基本の「き」「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1

連立一次方程式 Ax = b の解を、反復回数 k として

xk = Hxk−1 + c (k = 1, 2, · · · )なる漸化式を用いて逐次的に求める方法を「定常的反復法」と呼ぶ。

マントル対流の「数値」流体力学 1➢ここまでのまとめ➢流線関数➢原始変数解法➢反復解法➢ACuTE その 1
➢擬似圧縮性法とは?
➢ACuTE その 2
➢反復解法 2
➢Gershgorin の定理➢反復解法 3➢多重格子法➢多重格子法並列化➢箱型計算例

適切な H と c を選ぶ例の一つとして、係数行列 A を
A = M − N のように分解する。特に、M が正則かつ Nの全非零成分が正であるようにとることを「行列 A の正則分離」という。これによりマントル対流の「数値」流体力学 2とすればよい。

xk = M −1 N xk−1 + M −1 b (k = 1, 2, · · · )マントル対流の「数値」「流体」力学 2
マントル対流の「数2010 年夏のGFD セミナー
3値」流体力学

平成 22 年 8 月 21 日 – slide 36

連立一次方程式の反復解法: 基本の「ほ」「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1➢ここまでのまとめ➢流線関数➢原始変数解法➢反復解法➢ACuTE その 1
➢擬似圧縮性法とは?
➢ACuTE その 2
➢反復解法 2
➢Gershgorin の定理➢反復解法 3➢多重格子法➢多重格子法並列化➢箱型計算例
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2

行列 A を A = M − N と正則分離し、これらを用いて以下のような漸化式を考える。

xk = M −1 N xk−1 + M −1 b (k = 1, 2, · · · )ベクトル列 {xk } が収束するとしたら、その極限値 x∞ は
当然 Ax∞ = b を満たす。
さらに実用上は、M −1 が簡単に求まることが望ましい。
よく知られた Jacobi 法もその一例であり、M として Aの対角部分 D をとった場合に相当する。

Jacobi 法をはじめ多くの反復解法は D −1 を用いて近似解を更新する。巷間よくある問題 (スカラー変数のポアソン方程式の求解) にはこれで十分なのだが、A の対角成分に零を含むような場合には、これらの方法を用いることができない。

マントル対流の「数2010 年夏のGFD セミナー
3値」流体力学

平成 22 年 8 月 21 日 – slide 37

ACuTE 法 (Kameyama et al., 2005): 基本原理
「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1➢ここまでのまとめ➢流線関数➢原始変数解法➢反復解法➢ACuTE その 1
➢擬似圧縮性法とは?
➢ACuTE その 2
➢反復解法 2
➢Gershgorin の定理➢反復解法 3➢多重格子法➢多重格子法並列化➢箱型計算例
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2

構成要素その 1: 擬似圧縮性法 (Chorin, 1967)
与えられた温度 T 、粘性率 η の分布のもとで、高粘性・非圧縮性流体の定常流れ場を求める方程式
0 = −∇p + ∇ · [η(∇ ⊗ v + v ⊗ ∇)] + Ra T ez
0 = ∇·v
を直接解く代わりに、擬似的な時間発展方程式∂v
M= −∇p + ∇ · [η (∇ ⊗ v + v ⊗ ∇)] + Ra T ez
∂τ
∂p
−K
=∇·v
∂τ
を定常になるまで時間積分してやる。
ただし τ : 擬似時間
M : 擬似密度K : 擬似圧縮率 は「本物」とは無関係。

マントル対流の「数2010 年夏のGFD セミナー
3値」流体力学

平成 22 年 8 月 21 日 – slide 38

擬似圧縮性法とは?
「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1➢ここまでのまとめ➢流線関数➢原始変数解法➢反復解法➢ACuTE その 1
➢擬似圧縮性法とは?
➢ACuTE その 2
➢反復解法 2
➢Gershgorin の定理➢反復解法 3➢多重格子法➢多重格子法並列化➢箱型計算例
マントル対流の「数値」流体力学 2

長所❏ 非圧縮の速度場が必ず得られる (owing to 粘性による散逸)∂22η ∂  2
1
2(∇ · v) =∇ (∇ · v) +∇ (∇ · v) + · · ·2∂τ|KM {z} |M ∂τ {z}
擬似音波の伝播

粘性による散逸 (拡散)❏ 原始変数 (速度場 v と圧力場 p) をそのまま使って解く✏ 2 次元問題でも 3 次元問題でも OK✏ プログラムの構造が非常に直感的で簡単

短所❏ そのまま使うと非常に遅い✏ 特に長波長成分の収束が遅い (拡散方程式の宿命)

マントル対流の「数値」「流体」力学 2
マントル対流の「数2010 年夏のGFD セミナー
3値」流体力学⇒ 多重格子法との併用が不可欠
平成 22 年 8 月 21 日 – slide 39

ACuTE 法 (Kameyama et al., 2005): 粘性変化対策「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1➢ここまでのまとめ➢流線関数➢原始変数解法➢反復解法➢ACuTE その 1
➢擬似圧縮性法とは?
➢ACuTE その 2
➢反復解法 2
➢Gershgorin の定理➢反復解法 3➢多重格子法➢多重格子法並列化➢箱型計算例
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2

対策: 「局所時間刻み法」 の援用❏ 粘性率 η に応じて「密度」M と「圧縮率」K の大きさ

を「場所ごと」に変える✏ 実効的な時間刻み ∆τ /M 、∆τ /K を空間変化させる

ことに対応❏ 粘性率の空間変化の影響をなるべく打ち消すように、

M と K を空間変化させたい∂2
1
2η ∂  2
2(∇ · v) =∇ (∇ · v) +∇ (∇ · v) + · · ·2∂τ|KM {z} |M ∂τ {z}
擬似音波の伝播

粘性による散逸 (拡散)

✏ 実効的な拡散係数を一様にしたい ⇒M ∝ η✏ 擬似的な「音速」を一様にしたい ⇒K ∝ η −1
(スムーズな変化に対しては) この方法は非常に効果的

マントル対流の「数2010 年夏のGFD セミナー
3値」流体力学

平成 22 年 8 月 21 日 – slide 40

連立一次方程式の反復解法: 基本の「ん」「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1➢ここまでのまとめ➢流線関数➢原始変数解法➢反復解法➢ACuTE その 1
➢擬似圧縮性法とは?
➢ACuTE その 2
➢反復解法 2
➢Gershgorin の定理➢反復解法 3➢多重格子法➢多重格子法並列化➢箱型計算例
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2

連立一次方程式 Ax = b の解は、仮想的な時間発展方程式dx= b − Axdtの定常解とみなすこともできる。ならば、この方程式を定常になるまで時間積分してやることで、求める解が得られるはずである。
例えば、時間方向に陽的なスキームを用いてこの微分方程式を離散化すると、

xk − xk−1= b − Axk−1∆t(k = 1, 2, · · · )のような漸化式を得る。行列 A が正則で、かつ ∆t が十分小さい場合には、ベクトル列 {xk } は真の解に収束する(x∞ = A−1 b) ことが示せる。

マントル対流の「数2010 年夏のGFD セミナー
3値」流体力学

平成 22 年 8 月 21 日 – slide 41

連立一次方程式の反復解法としての ACuTE「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1➢ここまでのまとめ➢流線関数➢原始変数解法➢反復解法➢ACuTE その 1
➢擬似圧縮性法とは?
➢ACuTE その 2
➢反復解法 2
➢Gershgorin の定理➢反復解法 3➢多重格子法➢多重格子法並列化➢箱型計算例
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2

ところで、先の式中のスカラー ∆t を正則な行列 T に入れ換えて得られる式

T −1 (xk − xk−1 ) = b − Axk−1
も同じ極限値に収束する。実際、この式を変形すると

k
−1
k−1
−1
x − A b = (I − T A) x
−A bとなることからも確認できる。その際、行列 T をうまく選んでやれば、このベクトル列の収束を加速することができるはずである。
とはいえ最適な T (= A−1 ) を事前に知ることは不可能だ
から、何らかの形で近似してやる必要がある。ACuTE 法
では、T を適当な対角行列にとる ことによって、収束の加速を試みていることに相当する。

マントル対流の「数2010 年夏のGFD セミナー
3値」流体力学

平成 22 年 8 月 21 日 – slide 42

Gershgorin の定理「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1➢ここまでのまとめ➢流線関数➢原始変数解法➢反復解法➢ACuTE その 1
➢擬似圧縮性法とは?
➢ACuTE その 2
➢反復解法 2
➢Gershgorin の定理➢反復解法 3➢多重格子法➢多重格子法並列化➢箱型計算例

中心が aii 、半径が ri =|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
これをさらに書き直すと、|aii − λk | = −マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2

n
X

X
j6=i

aij

X
X
xj
xj
≤
aij
≤
|aij |
xi
x
i
j6=i
j6=i

を得る。即ち、λk を含む領域 Si が必ず存在する。

マントル対流の「数2010 年夏のGFD セミナー
3値」流体力学

平成 22 年 8 月 21 日 – slide 43

擬似時間刻みのとり方「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1➢ここまでのまとめ➢流線関数➢原始変数解法➢反復解法➢ACuTE その 1
➢擬似圧縮性法とは?
➢ACuTE その 2
➢反復解法 2
➢Gershgorin の定理➢反復解法 3➢多重格子法➢多重格子法並列化➢箱型計算例
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2

ACuTE の反復で出てくる 2 つの行列を
A=

n X
n
X

aij fi ⊗ f j ,i=1 j=1

T =nX

tl fl ⊗ f ll=1と書く。ただし fi 及び f i は空間離散化を表わすベクトルとする。反復行列 I − T A とその成分を書き表わすと、

I − TA =n XnX(δij − ti aij )fi ⊗ f j

i=1 j=1となる。この行列の第 i 行目に注目して Gershgorin の定理
を用いると、複素平面内の領域
Si は中心 (1 − ti aii , 0)、X
半径 ti|aij | の円となる。この円が実軸上で 0 から 1 の
j6=i

範囲に収まるように ti を選ぶとよい。

マントル対流の「数2010 年夏のGFD セミナー
3値」流体力学

平成 22 年 8 月 21 日 – slide 44

多重格子法 (マルチグリッド法)(Brandt, 1977; 以降多数の文献あり)「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1➢ここまでのまとめ➢流線関数➢原始変数解法➢反復解法➢ACuTE その 1
➢擬似圧縮性法とは?
➢ACuTE その 2
➢反復解法 2
➢Gershgorin の定理➢反復解法 3➢多重格子法➢多重格子法並列化➢箱型計算例
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2❏ 楕円型偏微分方程式を数値的に解く最高速の解法✏ 大規模な問題で威力を発揮する

計算量 O(N ): 他の手法では O(N log N ) や O(N 2 ) など
❏ 解像度の異なる格子での計算をうまく組み合わせることにより、細かい格子系での解を高速に求める⇔⇔❏ マントル対流問題以外にも広く適用可能 (事例多数)❏ ただし、大規模並列計算時にはもう一工夫が必要(例えば Kameyama, 2005 など)

マントル対流の「数2010 年夏のGFD セミナー
3値」流体力学

平成 22 年 8 月 21 日 – slide 45

多重格子法計算の並列化「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」

多重格子法の並列実行は agglomeration により効率化❏ 細かい 格子レベル (計算量大) では

全 PE

マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1➢ここまでのまとめ➢流線関数➢原始変数解法➢反復解法➢ACuTE その 1
➢擬似圧縮性法とは?
➢ACuTE その 2
➢反復解法 2
➢Gershgorin の定理➢反復解法 3➢多重格子法➢多重格子法並列化➢箱型計算例
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2❏ 粗いを使って 並列で実行

格子レベル (計算量小) では
1PE のみ を使って 非並列で 実行✍ 通信を「少量多数」ではなく「大量少数」にする
multigrid V-cyclefor 1024x1024x256 mesh problem

domain decompositionby 2x32x8=512

level 8(1024x1024x256)

512x32x32/PE

level 7(512x512x128)

256x16x16/PE

level 6
(256x256x64)

128x8x8/PE

level 5(128x128x32)
level 4(64x64x16)
level 3(32x32x8)level 2
(16x16x4)
level 1(8x8x2)

マントル対流の「数2010 年夏のGFD セミナー
3値」流体力学

64x4x4/PE
restrictionof residuals

gatherto rank 0

prolongationof errors

scatterfrom rank 0

decomposable=> activate all PEs

32x2x2/PE
32x32x8/PE
16x16x4/PE

not decomposable=> activate one PE only

8x8x2/PE

平成 22 年 8 月 21 日 – slide 46

粘性率の温度依存性がある場合の計算例
Case 1(Rt = 105 , r = 1)

「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1➢ここまでのまとめ➢流線関数➢原始変数解法➢反復解法➢ACuTE その 1
➢擬似圧縮性法とは?
➢ACuTE その 2
➢反復解法 2
➢Gershgorin の定理➢反復解法 3➢多重格子法➢多重格子法並列化➢箱型計算例
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2

Ogawa et al. (1991) との比較

1.7×0.5×1 box
64×32×64 mesh
η = ηt exp[E(Tt − T )]
ρgα(Tb − Tt )d3
Rt =
κηt
ηmax
r = exp[E(Tb − Tt )] =
ηmin

T=0.3

T=0.7

Case 4(Rt = 103 , r = 102 )Case 16(Rt = 103 , r = 3.2 × 103 )
T=0.5

T=0.3

T=0.7

T=0.92

Case 18(Rt = 32, r = 105 )T=0.5

T=0.92

計算時間: Case 1 (粘性コントラストなし) で 1 時間ステップあたり約 3 秒
Case 18 (粘性コントラスト5 桁 ) で 1 時間ステップあたり約18 秒
with Pentium IV 2.20GHz (ただし初期条件が異なるので、単純な比較は不可)

粘性コントラストが最大 10 桁ある場合まで動作確認済

マントル対流の「数2010 年夏のGFD セミナー
3値」流体力学

平成 22 年 8 月 21 日 – slide 47

3 次元箱型モデルでの計算例

マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1660km 相転移 + post-perovskite (PPv) 相転移を含めた計算(Kameyama and Yuen, 2006)(a) δT ≡ T − hT i(b) hT i, Tmax , Tmin
(2)case H02; Tbot = 3800K, Rb /Rasurf = 0.25; t = 1.12383057 × 10−31
0.8

マントル対流の「数値」流体力学 1➢ここまでのまとめ➢流線関数➢原始変数解法➢反復解法➢ACuTE その 1
➢擬似圧縮性法とは?
➢ACuTE その 2
➢反復解法 2
➢Gershgorin の定理➢反復解法 3➢多重格子法➢多重格子法並列化➢箱型計算例
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2

0.6
0.4
0.2
0
300
3800
temperature [K]dimensionless height「固体」地球内部の「流体」現象計算の諸元❏ メッシュ分割 512×512×128 (たぶん業界最高解像度)❏ 約 3.2 秒/ステップ (初代 ES の 128CPU 使用時) (たぶん業界最速)❏ 粘性率の温度・深さ依存性 (+熱拡散率の温度依存性)❏ 相転移はやや強めにしてある (相転移の効果を強調したかったので)
クラペイロン勾配 −4.3MPa/K、 密度ジャンプ 10% for 660km 相転移
13MPa/K、
2% for PPv 相転移

マントル対流の「数2010 年夏のGFD セミナー
3値」流体力学

平成 22 年 8 月 21 日 – slide 48

「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2➢Yin-Yang grid
➢球殻計算例 1➢インヤン並列化➢インヤン MG
➢インヤンの困難➢特異問題難点➢特異問題解法➢球殻問題➢特異問題 ACuTE

マントル対流の「数値」流体力学 2

マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3
マントル対流の「数値」流体「力学」

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 49

インヤン格子 (Kageyama and Sato, 2004)「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2➢Yin-Yang grid
➢球殻計算例 1➢インヤン並列化➢インヤン MG
➢インヤンの困難➢特異問題難点➢特異問題解法➢球殻問題➢特異問題 ACuTE
マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3

図はKageyama, A., Sato, T., 2004 より❏ 球座標系の「低緯度」部分のみを 2 つ組み合わせる❏ 2 つの合同な要素格子の重ね合わせで、球殻を覆う❏ 極での特異性なし、極付近での格子間隔の激減なし

マントル対流の「数値」流体「力学」

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 50

インヤン格子を用いたマントル対流「固体」地球内部の「流体」現象(Kameyama et al., 2008)

3 次元箱型モデルと同じ高速解法を球殻モデルにも適用

マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2➢Yin-Yang grid
➢球殻計算例 1➢インヤン並列化➢インヤン MG
➢インヤンの困難➢特異問題難点➢特異問題解法➢球殻問題➢特異問題 ACuTE
マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3
マントル対流の「数値」流体「力学」

figure and movie by courtesy of Dr. Nobuaki Ohno @ ESC/JAMSTEC

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 51

インヤン格子計算の並列化「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」

領域分割法による並列化❏ 空間 3 方向 (r, θ, φ) ×2 (「イン」と「ヤン」) に分割

マントル対流の「数値」「流体」力学 1

部分領域間の通信のうち、マントル対流の「数値」流体力学 1❏ 「イン」または「ヤン」内での通マントル対流の「数値」流体力学 2➢Yin-Yang grid
➢球殻計算例 1➢インヤン並列化➢インヤン MG
➢インヤンの困難➢特異問題難点➢特異問題解法➢球殻問題➢特異問題 ACuTE
マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3
マントル対流の「数値」流体「力学」

図はKageyama, A., Sato, T., 2004より

信は、ごく普通にやれば OK❏ 「イン」と「ヤン」をまたぐ通信には工夫が必要✏ ある方向へのシフト通信の相手が 1 つだけとは限らない✏ 多重格子法ではさらに複雑 (通信相手の数が格子レベルにより変化)

ただし、通信テーブルは「イン」と「ヤン」で使い回し OK

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 52

インヤン格子の多重格子計算「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1

Nθ = 4(Nφ = 3Nθ )マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2➢Yin-Yang grid
➢球殻計算例 1➢インヤン並列化➢インヤン MG
➢インヤンの困難➢特異問題難点➢特異問題解法➢球殻問題➢特異問題 ACuTE

Nθ = 8Nθ = 16

マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3
マントル対流の「数値」流体「力学」

Nθ = 32

2010 年夏の GFD セミナー

Kageyama, A., Sato, T., 2004より

平成 22 年 8 月 21 日 – slide 53

インヤン格子を用いたマントル対流の特殊な事情「固体」地球内部の「流体」現象流れ場を解く連立一次方程式の特異性を確実に除去すべし

マントル対流の数値「流体」「力学」

非圧縮性流体の流れ場を解く場合には、マントル対流の「数値」「流体」力学 1❏ 圧力 p に定数分の不定性

マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2➢Yin-Yang grid
➢球殻計算例 1➢インヤン並列化➢インヤン MG
➢インヤンの困難➢特異問題難点➢特異問題解法➢球殻問題➢特異問題 ACuTE
マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3✏ 非圧縮性流体では、圧力は「周囲との差」(∇p) の形

でしか現れない❏ 速度場 v に剛体運動 (並進・回転) 分の不定性✏ 剛体運動は歪を生まないので、粘性的な抵抗力の発生に寄与しない.... 即ち、「一意な解が得られる」保証が (もともと) ない。
箱型モデルの場合にはさほど気にしなくても大丈夫だが、インヤン格子で球殻形状モデルを扱う場合には要注意。

マントル対流の「数値」流体「力学」

2010 年夏の GFD セミナー(「補間」=「新たな誤差を混入させる」)
平成 22 年 8 月 21 日 – slide 54

特異な連立一次方程式の難点「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2➢Yin-Yang grid
➢球殻計算例 1➢インヤン並列化➢インヤン MG
➢インヤンの困難➢特異問題難点➢特異問題解法➢球殻問題➢特異問題 ACuTE
マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3
マントル対流の「数値」流体「力学」

ある正則な行列 T を用いて、連立一次方程式 Ax = b を

xk = xk−1 + T (b − Axk−1 ) (k = 1, 2, 3, · · · )のような反復で解くにあたり、残差 r k ≡ b − Axk の変化に注目してみよう。上式を書き直すと、

rk = (I − AT )rk−1⇔

rk − rk−1 = −AT rk−1

ベクトル列 {r k } が収束するとしたら、その極限値 r ∞ はAT r∞ = 0 を満たすはずである。
❏ A が正則であれば、r ∞ = T −1 A−1 0 = 0、即ち

b − Ax∞ = 0 を満たす真の解に収束する。
❏ A が特異であれば、Ac = 0 を満たす非零な c が存在する。もし T r ∞ = c であれば、r ∞ = T −1 c 6= 0 となり、真の解に収束しない。

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 55

特異な連立一次方程式を解く: 原理「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2➢Yin-Yang grid
➢球殻計算例 1➢インヤン並列化➢インヤン MG
➢インヤンの困難➢特異問題難点➢特異問題解法➢球殻問題➢特異問題 ACuTE
マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3
マントル対流の「数値」流体「力学」

特異な係数行列 A を持つ n 元連立一次方程式 Ax = b が解を持つ条件を考える。
A の固有値を λi (i = 1, · · · , n)、対応する一次独立な固有ベクトルを ei と書く。n 個の固有値のうち、λ1 , · · · , λmの m (< n) 個が零、λm+1 , · · · , λn の (n − m) 個が非零とする。
b と x を ei の一次結合で書いてみよう。

b=

n
X

cj ej ,

x=

j=1

n
X

yk ek

k=1

これをもとの方程式に代入すると、
!
X
X
X
X
cj ej = Ayk ek =
yk Aek =yk λk ekj2010 年夏の GFD セミナー

k

k

k

平成 22 年 8 月 21 日 – slide 56

特異な連立一次方程式を解く: 原理「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2➢Yin-Yang grid
➢球殻計算例 1➢インヤン並列化➢インヤン MG
➢インヤンの困難➢特異問題難点➢特異問題解法➢球殻問題➢特異問題 ACuTE

固有ベクトルの一次独立性より、全ての j について

cj = λj yj
でなければならない。ここで λ1 = · · · = λm = 0 であるから、右辺項 b の展開係数は

c1 = · · · = cm = 0を満たさなければならない。これが、与式に解が存在するための条件 (「適合条件」) である。
また、与式だけから y1 , · · · , ym の値を決定できないため、解を一意に定めるためには m 個の付加条件が必要である。

マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3
マントル対流の「数値」流体「力学」

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 57

特異な連立一次方程式を解く: 方法「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2➢Yin-Yang grid
➢球殻計算例 1➢インヤン並列化➢インヤン MG
➢インヤンの困難➢特異問題難点➢特異問題解法➢球殻問題➢特異問題 ACuTE
マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3

先の表記を使うと、階数 n − m (< n) の n × n 行列 A は

A=

m
X

λi ei ⊗ ei +

i=1

n
X

λj ej ⊗ ej =

j=m+1

n
X

λj ej ⊗ ej

j=m+1と表わされる。ただし ei は ei · ej = δij を満たすとする。
ここで、A の代わりに、以下で定義される

Ã ≡

m
X
i=1

νi ei ⊗ ei +

n
X

j=m+1

λj ej ⊗ ej = A +mXνi ei ⊗ ei

i=1

を考えよう。ただし ν1 , · · · , νm 6= 0 とする。定義より Ã
は正則であるから、n 元連立一次方程式 Ãx̃ = b は一意な解を持つ。

マントル対流の「数値」流体「力学」

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 58

特異な連立一次方程式を解く: 方法「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1

x̃ =

n
X

y˜k ek と展開してもとの方程式に代入し、

k=1

ν1 , · · · , νm , λm+1 , · · · , λn 6= 0 及び A の固有ベクトルの一次独立性を用いればマントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2➢Yin-Yang grid
➢球殻計算例 1➢インヤン並列化➢インヤン MG
➢インヤンの困難➢特異問題難点➢特異問題解法➢球殻問題➢特異問題 ACuTE

k=1, · · · ,m について y˜k = ck /νkl = m + 1, · · · ,n について ỹl = cl /λl = ylが得られる。
さらにx = x̃ −マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3mX
k=1

y˜k ek = x̃ −mX

(x̃ · ek )ek

k=1などと付加条件を課すことで、一意な解 x が得られる。

マントル対流の「数値」流体「力学」

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 59

3 次元球殻内マントル対流を解くために「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2➢Yin-Yang grid
➢球殻計算例 1➢インヤン並列化➢インヤン MG
➢インヤンの困難➢特異問題難点➢特異問題解法➢球殻問題➢特異問題 ACuTE
マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3
マントル対流の「数値」流体「力学」

球殻内のマントルの流れ場を記述する連立一次方程式を


 

Crr Crθ Crφ Gr
vrbr Cθr Cθθ Cθφ Gθ   vθ   bθ 


=

 Cφr Cφθ Cφφ Gφ   vφ   bφ 
Dr Dθ Dφ 0pbp{z
} | {z } | {z }
|
≡A

≡x

≡b

のように書いておく。
A の固有値 0 に対する固有ベクトルを与える状態は 剛体
回転 3 成分 と 圧力勾配 0 の 4 つが考えられる。これら 4つの「固有状態」に相当する不定性を取り除けばよい。
なお、剛体的な並進運動は、動径方向速度の境界条件により自然に除去される。

またこれらの拘束条件は、インとヤンの補間の際に混入する誤差の影響を除去する上でも重要。

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 60

特異な連立一次方程式用の ACuTE 1 of 4
「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2➢Yin-Yang grid
➢球殻計算例 1➢インヤン並列化➢インヤン MG
➢インヤンの困難➢特異問題難点➢特異問題解法➢球殻問題➢特異問題 ACuTE

先の考察に従い、特異な係数行列 A から派生した正則な係数行列 Ã を持つ連立一次方程式 Ãx̃ = b を、ACuTE に倣った反復法で解くことを考える。

−1
k
k−1
T
x̃ − x̃
= r̃k−1 = b − Ãx̃k−1 (k = 1, 2, 3, · · · )ここで、残差 r̃ k−1 の表式を考えると、"#mX
νi ei ⊗ ei x̃k−1
r̃k−1 = b − A +
i=1

k−1

r k−1

マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3
マントル対流の「数値」流体「力学」= |b − A{zx̃ } −

m
X
i=1

k−1

νi x̃



· ei eiとなる。即ち、残差の定義を上のように修正しておけば、近似解を更新する手順は ACuTE 法と全く同じで OK。

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 61

特異な連立一次方程式用の ACuTE 2 of 4
「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2➢Yin-Yang grid
➢球殻計算例 1➢インヤン並列化➢インヤン MG
➢インヤンの困難➢特異問題難点➢特異問題解法➢球殻問題➢特異問題 ACuTE
マントル対流の「数値」「流体」力学 2

r̃

k−1

=r

k−1

−

m
X

k−1

νi x̃

i=1



· ei ei

このうち、
❏

k−1



x̃
· ei ei は、近似解 x̃k−1 に含まれる、i 番目の特異性に関する成分 (ただし i = 1, · · · , m)。
球殻内熱対流の場合には、以下の 4 つを指す。
✏ 速度場 v のうち、ある角速度 ω = (ω1 , ω2 , ω3 ) で与えられる剛体回転 3 成分
✏ 圧力場 p に含まれる「ゲタ」1 成分❏ 非零なパラメータ νi は、大きすぎても小さすぎても×。

マントル対流の「数値」流体力学 3✏ 小さすぎると、特異性が全く除去されない。

マントル対流の「数値」流体「力学」✏ 大きすぎると、反復の安定性が破れる。

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 62

特異な連立一次方程式用の ACuTE 3 of 4
「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2➢Yin-Yang grid
➢球殻計算例 1➢インヤン並列化➢インヤン MG
➢インヤンの困難➢特異問題難点➢特異問題解法➢球殻問題➢特異問題 ACuTE
マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3

特異性の除去を含める場合には、残差 r の更新に使われる実効的な反復行列の表式を考えよう。
(
)
)( n
m
X
X
I − ÃT = I − A +νl el ⊗ el
tk fk ⊗ f k= I − AT −

l=1
"
n X
n
m
X
X
i=1 j=1

l=1

k=1

#

νl tj (f i · el )(el · fj ) fi ⊗ f j

"
#
n
n
m
XX
X
=
νl tj (f i · el )(el · fj ) fi ⊗ f j
δij − aij tj −i=1 j=1

l=1

このうち、tj (f i · el )(el · fj ) は、j 番目の変数の値が tj であるとき、l 番目の特異性の大きさによって生じる、i 番
目の変数の値 という意味をもっている。

マントル対流の「数値」流体「力学」

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 63

特異な連立一次方程式用の ACuTE 4 of 4
「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2➢Yin-Yang grid
➢球殻計算例 1➢インヤン並列化➢インヤン MG
➢インヤンの困難➢特異問題難点➢特異問題解法➢球殻問題➢特異問題 ACuTE

特異性を除去する効果の強さを表わす νl は、実効的な反
復行列に
n X
n X
m
X
i=1 j=1 l=1という新たな行列要素を足し合わせるが故に収束性に影響し、絶対値が大きいほど反復計算の安定性を損う。
実際の計算では、

νl = ǫ ×

マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3
マントル対流の「数値」流体「力学」

νl tj (f i · el )(el · fj )1nX

(ǫは小さい定数)tj (el · fj )j=1のようにとっている。

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 64

「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2

マントル対流の「数値」「流体」力学 2

マントル対流の「数値」「流体」力学 2➢粘弾性流体➢CIP-CSLR
マントル対流の「数値」流体力学 3
マントル対流の「数値」流体「力学」

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 65

粘弾性流体向けの新しい解法の開発「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2➢粘弾性流体➢CIP-CSLR

粘弾性流体の大変形を扱うための技術的課題を克服(Furuichi, Kameyama and Kageyama, 2008)

粘性+弾性レオロジー⇒ 共回転マクスウェルモデル
物性 (粘性率) の急激な変化⇒ ACuTE 法(境界面) の存在
境界面の追跡と大きな時間幅⇒ 半解析的時間積分法 + CIP 法
独自開発の
Jaumann (共回転) 微分の取り扱い ⇒「共回転セミラグランジュ法」

粘弾性流体の Rayleigh-Taylor 不安定

マントル対流の「数値」流体力学 3
マントル対流の「数値」流体「力学」

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 66

高精度移流スキームとその応用「固体」地球内部の「流体」現象Fluid Rope Coiling

マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1(Furuichi, Kameyama and Kageyama, 2009)

粘性率と密度の異なる 2 種類の流体のストークス流を同時に計算

マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2➢粘弾性流体➢CIP-CSLR
マントル対流の「数値」流体力学 3
マントル対流の「数値」流体「力学」

数値拡散の少ない
保存性移流スキーム の 開 発 に よ り、
(さほど高くない空間解像度

シャープ
な流体界面や質量保存を破ることなく追跡できているであっても)(あまり空間解像度を上げられない 3次元マントル対流計算には有難い)

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 67

「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2

マントル対流の「数値」流体力学 3

マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3➢GMRES
マントル対流の「数値」流体「力学」

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 68

GMRES 法あるいは残差切除法 (その 1)「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3➢GMRES
マントル対流の「数値」流体「力学」

線型な連立一次方程式 f = Au を解くことを考える。
何らかの手段 (例えば過去 n 回のマルチグリッド法のサイクル) によって用意された近似解を u1 , u2 , · · · , un とする。un とその直近の m 個 (m < n) の近似解の線形結合により、新しい近似解 uacc を作ることを考える。
nuacc ≡ u +

m
X
i=1

αi u

n−i

−u

n



本当は un−m , · · · , un が直交しているとなおよろしい (それを保証すると「正統的な」GMRES 法になる?) のだが、とりあえず現時点ではそこまでは要求しないことにする。
この新しい近似解 uacc が、可能な限りよい近似解になる
ように、係数 αi (i = 1, · · · , m) を決めてやろう。

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 69

GMRES 法あるいは残差切除法 (その 2)「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1

新しい近似解 uacc によって与えられる残差は、
"
#
m
X

n
n−i
n
racc ≡ f − Auacc = f − A u +αi u − ui=1n= f − Au −| {z }

マントル対流の「数値」流体力学 2=r n

マントル対流の「数値」「流体」力学 2nマントル対流の「数値」流体力学 3➢GMRES
マントル対流の「数値」流体「力学」

=r +

m
X
i=1

m
X

n−i

n



αi A u − u
|
{z
}
i=1

αi r

n−i

=−r n−i +r n

−r

n



となる。racc の 2 乗ノルムが小さければ小さいほどよい近似解だと言えるので、これを最小化するような係数 αi の組を探そう。

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 70

GMRES 法あるいは残差切除法 (その 3)「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」"||racc ||2 = rn +マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2"· rn +

i=1mX

αk rn−k − rn= (r · r ) + 2

m
X
i=1

+

m
m X
X
i=1 k=1

2010 年夏の GFD セミナー

αi rn−i − rk=1n

マントル対流の「数値」流体力学 3➢GMRES
マントル対流の「数値」流体「力学」

m
X

αi αk





#

n

#

n
n

αi (r · r
n−i

n−k

n−i

n

n

) − (r · r )
n

n−i

(r
· r ) − (r · r )−(rn · rn−k ) + (rn · rn )

平成 22 年 8 月 21 日 – slide 71GMRES 法あるいは残差切除法 (その 4)「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3➢GMRES||racc ||2 の極小値 (より正確には停留値) を探すために、αiで微分すると、
 n n−i∂||racc ||2nn= 2 (r · r ) − (r · r )∂αi
 n−i n−j

m
n
n−i
X
(r
· r ) − (r · r )αj+2−(rn · rn−j ) + (rn · rn )j=1

i = 1, · · · , m の全てについて ∂||racc ||2 /∂αi = 0 とおくと、最小化問題を記述する m 元連立一次方程式

マントル対流の「数値」流体「力学」

Hij αj = βi(i, j = 1, · · · , m)

(1)を得る。これを解いてやることで、最適な α1 , · · · , αm が定まり、その結果最適な uacc も定まる。
2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 72GMRES 法あるいは残差切除法 (その 5)「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3➢GMRES
マントル対流の「数値」流体「力学」❏ リスタートつきの一般化最小残差法 (GMRES(m) 法) の前処理にマルチグリッド法を使ったことに相当✏ 厳密に言えば GMRES 法そのものではない (基底の直交化をしてない)✏ 残差切除法 (Residual Cutting Method) に近い?
❏ m は大きすぎても小さすぎても×✏ m が大きいと、メモリ使用量が増える + 最小化問題

が解きにくくなる✏ m が小さいと、最適解を探索する線型空間の次元が小さく、よい解に到達しにくい
ただしモノの本によれば、多重格子法と併用する場合には m = 5 程度で十分らしく、亀山が試した限りでもそのような感じがする。

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 73

「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2

マントル対流の「数値」流体「力学」

マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3
マントル対流の「数値」流体「力学」➢これから➢superplumes?
➢プレート?
➢「粗視化」妄想2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 74

これからのマントル対流シミュレーション私案「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」

問: マントル対流シミュレーションで得られる描像は、その他の (例えば地震学的観測による) 手法で得られた現在のマントル内部の描像と調和的か否か?

マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3
マントル対流の「数値」流体「力学」➢これから➢superplumes?
➢プレート?
➢「粗視化」妄想

数値シミュレーション

中国北西部 ∼ 中部日本

タヒチ ∼ ハワイ

答はおそらく....
❏ “Yes” in some aspects.
❏ But in (many) other aspects, possibly “No”.

そもそも、最も特徴的な構造ですら、十分に再現されているとも言い難い2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 75How to reconcile with “Superplumes”?
「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」

少数・大規模な上昇流が卓越するようになるためには ....
下部マントルの物性モデル がカギ

マントル対流の「数値」「流体」力学 1


熱膨張率 (or 状態方程式)



粘性率
熱伝導率



組成・相状態による違い

マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3
マントル対流の「数値」流体「力学」➢これから➢superplumes?
➢プレート?
➢「粗視化」妄想

下図は Yuen, Monnereau, Hansen, Kameyama and Matyska (2006) より。
この計算では熱膨張率・粘性率の深さ (圧力) 依存性、熱伝導率の温度依存性を強くしてある。

❏ 高温高圧下でのマントル物質の性質 (状態方程式など) に関する物質科学的知見もまだまだ十分とはいえない❏ マントル対流シミュレーションから (ある程度) 拘束をかけてやることができるはず2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 76

プレート運動とマントル対流を同時に「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1

プレートの挙動を再現する上では、局所的な粘性率の低下が非常に必要例えば、「滞留スラブ」の再現に重要なものは❏ スラブ内の粘性を下げる

マントル対流の「数値」流体力学 1(スラブを曲げやすくする)

マントル対流の「数値」流体力学 2❏ プレート境界に薄い低粘性層を置く

マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3
マントル対流の「数値」流体「力学」➢これから➢superplumes?
➢プレート?
➢「粗視化」妄想(上盤プレートとの結合を弱める)
図は Tagawa et al. (2007)より拝借

しかし、3 次元モデルで正確に解くことは非常に困難❏ 解くべき問題が非常に「悪条件」になる❏ 多重格子法 (3 次元計算の高速化に不可欠) では、粗い格子上で「見えない」微細構造を扱うことが本質的に不可能
2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 77

プログラム改良計画妄想「固体」地球内部の「流体」現象マントル対流の数値「流体」「力学」
マントル対流の「数値」「流体」力学 1
マントル対流の「数値」流体力学 1
マントル対流の「数値」流体力学 2
マントル対流の「数値」「流体」力学 2
マントル対流の「数値」流体力学 3
マントル対流の「数値」流体「力学」➢これから➢superplumes?
➢プレート?
➢「粗視化」妄想「粗視化」のような考え方を多重格子法に入れられるか?
❏ 多重格子法で出てくる relative truncation error τh2h

L2h (u2h ) = f2h|{z}

粗い格子系での離散化方程式+τh2h

τh2h には、粗い格子系での解を細かい格子系での解と一致させるという意味があるということは、❏ 粗い格子系で粘性の微細構造が「見えなく」なったツケを τh2h に押しつけるような定式化ができればよい???
❏ よくある「見かけの力」と同じようなもの?

モデルより細かい時空間スケール
いま考えている
モデルでは含まれない素過程の現象の効果を取り込むようにできる?
(e.g.) 震源時間関数、等価物体力、dislocation theory

2010 年夏の GFD セミナー

平成 22 年 8 月 21 日 – slide 78

参考文献• Kageyama, A. et al. Kageyama, A., Kameyama, M., Fujihara, S., Yoshida,M., Hyodo, M., Tsuda, Y., 2004: A 15.2 TFlops simulation of geodynamo on
the Earth Simulator, Proceedings of the ACM/IEEE Supercomputing
SC2004 conference, p.35-35• Ichikawa, H., Labrosse, S., Kurita, K., 2010: Direct numerical simulation of
an iron rain in the magma ocean, Journal of Geophysical Research, 115• Dziewonski, A. M. , Anderson, D. L., 1981: Preliminary reference Earthmodel, Phys. Earth Planet. Inter., 25, 297-356• Zhao, D., 2004: Global tomographic images of mantle plumes andsubducting slabs: insight into deep Earth dynamics, Phys. Earth Planet.
Inter., 146, 3-34
• Scherneck, H. G., Johansson, J. M., Vermeer, M., Davis, J. L., Milne, G. A.,
Mitrovica, J. X., 2001: BIFROST project: 3-D crustal deformation rates
derived from GPS confirm postglacial rebound in Fennoscandia, Earth,Planets and Space, 53, 703-708

• Kono, M., Roberts, P. H., 2002: Recent geodynamo simulations andobservations of the Geomagnetic Field, Rev. Geophys., 40, 1• King, S. D., Lee, C., Keken, P. E., Leng, W., Zhong, S., Tan, E., Tosi, N.,
Kameyama, M., 2010: A community benchmark for 2-D Cartesiancompressible convection in the Earth's mantle, Geophys. J. Int., 180, 73-87• Kageyama, A. and Sato, T., 1995: Computer simulation of amagnetohydrodynamic dynamo. II, Phys. Plasmas 2, 1421
• Schubert, G., Turcotte, D. L., Olson. P., 2001: Mantle convection in theearth and planets, Cambridge University Press, 940• Kameyama, M., 2005: ACuTEMan: A multigrid-based mantle convection
simulation code and its optimization to the Earth Simulator, Journal of the
Earth Simulator, 4, 2-10
• Kageyama, A., Sato, T., 2004: "Yin-Yang grid" An overset grid in sphericalgeometry, Geochem. Geophys. Geosyst., 5, 15• Furuichi, M., Kameyama, M., Kageyama, A., 2009: Validity test of a Stokes
flow solver by fluid rope coiling: Toward plate-mantle simulation, Physics
of the Earth and Planetary Interiors, 176, Issue 1-2, p. 44-53

• Tagawa, M., Nakakuki, T., Kameyama, M., Tajima, F. 2007: The role ofhistory-dependent rheology in plate boundary lubrication for generating
one-sided subduction, Pure and Applied Geophysics, 164, 879-907• Yuen, D. A., Monnereau, M., Hansen, U., Kameyama, M., Matyska, C.,
2007: Dynamics of superplumes in the lower mantle• 唐戸 俊一郎, 2000: レオロジーと地球科学, 東京大学出版会, 251 pp• 鳥海 光弘, 谷本 俊郎, 高橋 栄一, 本蔵 義守, 玉木 賢策, 本多 了, 巽 好幸 1997: 岩波講座 地球惑星科学〈10〉地球内部ダイナミクス, 岩波書店,
268 pp