MHD数値シミュレーションの基礎 <講義
Download
Report
Transcript MHD数値シミュレーションの基礎 <講義
MHD数値シミュレーションの基礎
•
•
•
•
ニュートン流体力学方程式
(ニュートン)磁気流体力学方程式
相対論的流体力学方程式
相対論的磁気流体力学方程式
• 千葉大学・花輪知幸
一部のスライドは 高橋博之さんからの借用
保存形式
U Fx Fy Fz
S
t
x
y
z
U UV Fx Fx V
v
x
V vy , U
v
z
vx
vy
vz
v
S SV, r, t
2
, Fx
2
vx
2
v
P
x
vxvy
v x vz
2
v 2 P vx
重力、遠心力、加熱・冷却
相対論的磁気流体力学
ー保存形式ー
T
n u
F
g
,
1
g,T
g
2
g
g
,
,
0
0
原始変数
n, , u , u , u ,
B1, B2 , B3
1
F
2
u 0
3
Gauss の法則 (Stokes の定理)
U
F 0
t
U
t dV F dS
U
t
dV U 0dV F dS dt
t
0
保存量と流束
保存量はセル平均の値
Fi 1/ 2, j ,k
U i , j ,k
U i , j ,k
U
dV
Vi , j ,k
Fi 1/ 2, j ,k F dS dt
評価点を変える
, P
Staggered Mesh
Constrained Transport
v, B
は表面で評価
o 長所
o 自然に空間2次精度、磁束の保存
o 短所
o 補間による拡散、モードの分離
数値天文学でふつうの解き方
• セルの面ごとに流束を計算し、その和により変化さ
せる (方向別分離)
– fractional time step は流行っていない
• セルの境界とベクトルは垂直
– 円筒座標なら
vr r, , z
– 航空などではセルを斜めにしても、」ベクトルはデカルト
のまま。
• 拡散や加熱冷却は、流体の計算と別ステップ
• 流体は陽解法で時間変化を求める。
• ガスは圧縮性
セル内の物理量の分布
• 高階微分は変数としない。
– CIPなどでは1階微分も計算
• セル幅より短波長の波を考えない
– サブグリッドモデルは採用しない
• 適当な物理量が一定とする
– 角速度、温度などが一定とする
– セル内に圧力勾配を考える( LeVeque の方
法)
課題
• 原始変数Vから 流束Fを求める
– 空間精度 (2次以上)
• 衝撃波に対する安定性
• 負の密度、負の圧力の排除
• カーバンクル不安定
– 時間精度 (2次以上)
• 磁気モノポール(div B) の消去
• 保存量U から原始変数V を求める
V → F → U 更新 → V
数値流束計算で考えること
• 流束を求める点での物理量を求める
– 物理量はセルでの平均値
• 流束も時間変化する
– 時間方向に積分された流束が必要
LL
L
R
RR
単純平均は失敗
ー奇妙な数値振動の発生ー
x j 1/ 2
x j 1 x j
x j 1 x j 2
x
x j 3
1
F x j 1/ 2 F U x j F U x j 1
2
1
F x j 1/ 2 F U x j U x j 1
2
方向別分離:
1次元磁気流体力学方程式
vx
2
B
B
B
x x
vx2 P
8
4
v
Bx By
x
vxvy
4
vy
BB
U vz , F vxvz x z
4
By
v
B
v
B
x
y
y
x
B
z
vx Bz vz Bx
e
2
B
e P
v B v B
x
x
8
7成分
div B = 0 の
制約があるた
めに自由度
が減る。
(Sod の)衝撃波管問題
Riemann の解
ρ
P
By
セル境界の値は「平均」でない!
t t
t
t F (x j1/ 2 , t) dt 2 F x j , t F x j1, t
データはBalsara (1998)
Godunov type solvers
波の伝播を考慮する
精密
面倒• Godunov (厳密解)
– Fast×2, Slow ×2, Alfven ×2, Entropy
• Roe (振幅の滑らかな変化を無視)
• HLLD (Fast と Slow を区別しない)
• HLLC (Fast, Slow, Alfven を区別しない )
• HLL (中間状態をひとつだけ考える)
簡単
粘性大
Riemann 解と特性速度
f
f
f
s
ρ
A
A
ξ = t/x
0
f
s
f
特性速度
U F
U F U
0
0
t
x
t U x
対角化
1 0
0 2
1 F
R
R
.. ..
U
0 0
.. 0
.. 0
Λ
.. ..
0 n
U
U
Λ
0
t
x
dU R1 dU
U
U
c
0 U ( x, t ) f x ct
t
x
Roe
UL
U1 U2 U3
U4
U
U
U
U
R
6
f
R
6
U
F U 6 U5 A U 6 U5
U
F
U5 U6 UR
衝撃波の
Jump
Condition
..............
Unとλn は UR
F U1 U L f U1 U L
と ULの関数
U
HLL
FR
UL
F
UM
HLL
FL
UR
Lt 0
保存則
F
Rt 0
FHLL FL L UM UL
FHLL FR R UM UR
HLL
R FL LFR R L UR UL
R L
λL とλR の
とり方に多
様性
特性速度
速い磁気音波 f vx c f
アルフヴェン波
遅い磁気音波
エントロピー波
A vx cA
s vx cs
0 vx
2
2
2
B
B
B
B
P P 2
x
y
z
2
2
cA
, cs
, cA
4
s
4
2
x
c f とcs は 4 cs2 cA2 2 a 2 cA2 0 の解
セル境界での特性速度
Roe 平均
L vxL L vxR
R PL L PR
vx
, P
L R
L R
UL
UR
L vyL L vyR
R ByL L ByR
vy
, By
L R
L R
Q QR QL , R L
vx vx vx
保存量 U は W の自乗で表せる
P By Bz
W , vx , vy , vz ,
,
,
固有モード
7
U wk rk
k : 特性速度
k 1
rk : 固有ベクトル
固有モード
7
F kwk rk
k 1
f
A
s
c
s
A
f
x
固有モード
P, vx , v// , B//
v , B
• 縦波 f, s
• 横波 A
• エントロピー波 c
f
A
s
c
s
A
f
x
横波とエントロピー波は区別しやすい
z
By
y
B// B
Bz
z
y
z
vy
y
v// v
vz
z
y
y
z
By
By2 Bz2 Bz
1
c P s
s P
2
s
縦波: Fast & Slow
s
f
c
v
c
v
s
x
s
f
x
f
s vy f y cs sgnBx
f vy s y cA sgnBx
v c sgnB
v c sgnB
x
z s
f
x
s z A
f z
s z
s y c f
,
, rs f y c f
rf
4
4
c
c
z f
f
s z f
4
4
vx2 vy2 vz2
vx2 vy2 vz2
f
f
gs
gf
2
2
f
c 2f cA2
c c
2
f
2
s
, s
c 2f a 2
c 2f cs2
Ryu & Jones (1995)
特性速度は大きめに評価
2
v
2
a 1 H
2
P
P v
H
1
2
2
速度差は実効的音
速を増やす。
CA'* も平均の磁場
からけいさんされる
より大きい
Cargo & Gallice (1998)
一般の状態方程式 Nobuta & TH, Mikami et al.
HLL 系では付加条件が必要
• 波の分解が不十分なので、整合条件だけで
は足りない。
• 付加的な条件が必要
• 特性速度を大きめにとると、密度・圧力が正
に保たれる (cf. Miyoshi et al.)
どこで特性速度を見積もるか
k , L , k , k , R
• 絶対値が大きいものが良い。
• 符号が同じなら平均で十分(とくに滑らか
な変化の場合は平均が良い)
• 符号が変わる場合は、平均より絶対値を
大きめにとる(エントロピーフィックス)
特性速度と数値粘性
FROE
F
HLL
7
1
FL FR k wk rk
2
k 1
1
R L
R L
FR FL
UR UL
FL FR
2
R L
R L
F U x U
D x
数値粘性は必要悪
• 数値粘性は、波を熱に変える。
– 下げないと熱発生源
– 実効的な分解能を低下させる
– 計算量を増加させる
• 数値粘性をなくすと、余計な波を発生
Courant 条件
•
•
•
•
•
•
波の伝播速度による制限
Courant数 |λ| Δt / Δx < 1
隣の隣に波が伝播しない
多次元では条件をきつくする
磁気力優勢の領域ではきつめに
許される範囲内でΔt は大きめに
境界条件
• 対称境界は容易 (鏡像をつくる)
– 無駄な領域(袖)をつくる
– 袖の値は、元の領域から対称性により求める。
– 流束の計算よりさきに境界を処理
• 反射のない境界(消極的な境界)は難しい
Godunov 法でも誤差は残る
• 平均化による誤差
x x / 2
1
U x, t t
U x, t t dx
x xx / 2
– 混合によるエントロピー増加
• セル内を一様と仮定した誤差
– 空間1次精度
空間2次精度
• 線形補間ではうまくゆかない
– Godunov の定理
• TVD: total variation dimishing
– 「波の」単調性 (数値振動なし)
空間2次精度 (高次精度化)
• 変数が階段的に変わるとしたことが原因
• 変数が滑らかに変化していることを考慮
• 単純平均はだめ
– 元の木阿弥
– 「風上」性が失われる
• 補間のしかたに工夫が必要
– 面倒な公式が不可避
•
Godunov の定理 線形補間は×
風上性を考慮して外挿
極大・極小で不自然な振動
単調性の保持 (TVD)
total variation diminishing
Minmod 補間
穏健な補間
勾配は少なめに
どの変数を補間するべきか
•
•
•
•
•
原理的には波の振幅(δwk) を補間すべき
流体力学では 原始変数を補間している。
保存量の補間は良くない。
温度より、圧力がまし。
磁気流体では、縦波・横波に分けて補間した
ほうが良い。
– 大振幅アルフヴェン波のテスト (Fukuda & TH)
• 角速度を補間すると良い。(対称軸の付近)
• 座標値やシフトベクトルは補間しない
時間2次精度
• 時刻とともに流束は変化する。
• [t, t+Δt] で積分した値を求める。
dx
f t の場合
dt
x t t x t f t t O t 2
x t t x t f t t 2 t O t 3
t
x t t x t f t f t t O t 3
2
進んだ時刻での物理量
• Δt/2 進んだセル境界での原始変数を推定
– MUSCL-Hancock
n1/ 2
, S
V
1
V
2
n
t n
n
I
A
V
x
• Δt/2 あるいはΔt 進んだセル中心の解をもとめる。
– Δt 進んだ解を使うほうが衝撃波に強い
U
n
U
U t
t
n
n
方向分離の問題
•
•
•
•
波面の向きによる依存性
斜めの波は2成分に分離
斜めの情報は直接入らない
∇・Bの発生源
Bx,i 1, j ,k Bx,i 1, j ,k
B
2x
B
By,i, j 1,k Bz ,i , j ,k 1 Bz ,i, j ,k ^1
y,i , j 1,k
2y
2x
div B の消去
Hyperbolic Divergence Cleaning
Dedner et al., JCP, 175, 645 (2002)
v 0
t
2
B
BB
I
0
v vv P
t
8
4
2
B
e
B
v
B
v
0
e P
t
8
4
B
vB Bv 0
t
B 0
古典的な ∇・B消去
B 0
B B
B 0
流体と一緒に流す
拡散させる
双極子磁場を付加
GLM形式
拡散しながら伝播
v 0
t
2
B
BB
I
0
v vv P
t
8
4
2
B
e
B
v
B
v
0
e P
t
8
4
B
vB Bv 0
t
ch2
2
ch B 2
t
cp
付加関数 ψ
付加関数の電信方程式
2
c
2
2
2
h
ch 2
v B
2
t
c p t
1ステップ弱で隣に
モノポールを移動。
数ステップで消滅
GLM形式の成分表示と長所
• 陽解法
• Poisson 方程式は解かない
• 遠方への影響は限定的
RMHD近似リーマン解法の発展
•
1983 HLL スキーム(Harten & Lax & van Leer ‘83)
•
2003 HLLスキームをRMHDへ応用(Gammie et al. ’03)
•
2005 exact Riemann solver (Bx=0) (Romero et al. ’05)
•
2006 exact Riemann solver (Bx\=0) (Giacomazzo et al. ’06)
•
2006 HLLC スキーム (Mignone & Bodo ’06, Honkkila Janhunen ’06)
•
2009 HLLD スキーム (Mignone et al. ’08)
•
2010 Roe スキーム (Anton et al. ’10)
ideal RMHDのスキームはだいぶ
確立されている
保存量Uから原始変数Vへ
•Balsara ’01 5x5 Matrix Jacobian
•Komissarov ’99: 3x3 nonlinear eqns.(密度、圧力、速度の絶対値)
•del Zanna et al. ’03: 2x2 nonlinear eqns. (a equation?)(圧力、速度の絶対値)
•Mignone & Bodo ’06: a nonlinear equation.
を独立変数とした1次元の非線形方程式
この式をNewton-Raphson法で解く
(Mignone & Bodo ’06の論文には式の間違いが多いので注意)
γ>100を超えるにはまだ難しい、、、
知っておくと良い知識
•
•
•
•
•
危険箇所
カーバンクル不安定
偽の加熱
既知の磁場と摂動
計算コストの削減
– パレトの法則
カーバンクル不安定
等速度面 (TH, Mikami, Matsumoto)
危険箇所
• 伝播速度λが 0 に近いところ
– 特性曲線が集まるところ 衝撃波
– 特性曲線が開くところ 膨張衝撃波の危険
– 定在衝撃波は解きづらい
• 磁気力が優勢なところ
– わずかな磁場の変化が流体に大きな変化をあ
たえる。
– force free は計算しづらい
偽の加熱
• Godunov 型で静水圧平衡を普通に解くと偽の
加熱が発生
2
B
e
v v g
e P
t
8
ρ v は「数値流束」を用
いると良い。セル中心の
値は良くない。 LeVeque
(2003), Mikami et al. (2008)
既知の磁場を差し引く
B B0 B1
B0 0 のとき有効
強い双極子磁場をもつ場合
球座標で一様磁場をとく場合
Tanaka
計算コストの削減
• 実効的な分解能の2倍低下を、細かい格子で
補うにと16倍の計算量増加
• 一般には高次精度のほうが得
• コストを決めているのは2割の要素
– パレトの法則 (2:8の法則)