Transcript MHD数値シミュレーションの基礎 ＜講義

MHD数値シミュレーションの基礎
•
•
•
•
ニュートン流体力学方程式
（ニュートン）磁気流体力学方程式
相対論的流体力学方程式
相対論的磁気流体力学方程式
• 千葉大学・花輪知幸
一部のスライドは 高橋博之さんからの借用
保存形式
U Fx Fy Fz



S
t
x
y
z
U UV Fx  Fx V

 

v 

 x
V  vy  , U  

v 

 z
 
  


vx
vy
vz
 v
S SV, 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 0dV    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 自然に空間２次精度、磁束の保存
o 短所
o 補間による拡散、モードの分離
数値天文学でふつうの解き方
• セルの面ごとに流束を計算し、その和により変化さ
せる （方向別分離）
– fractional time step は流行っていない
• セルの境界とベクトルは垂直
– 円筒座標なら
vr r, , z 
– 航空などではセルを斜めにしても、」ベクトルはデカルト
のまま。
• 拡散や加熱冷却は、流体の計算と別ステップ
• 流体は陽解法で時間変化を求める。
• ガスは圧縮性
セル内の物理量の分布
• 高階微分は変数としない。
– CIPなどでは１階微分も計算
• セル幅より短波長の波を考えない
– サブグリッドモデルは採用しない
• 適当な物理量が一定とする
– 角速度、温度などが一定とする
– セル内に圧力勾配を考える（ LeVeque の方
法）
課題
• 原始変数Vから 流束Fを求める
– 空間精度 （２次以上）
• 衝撃波に対する安定性
• 負の密度、負の圧力の排除
• カーバンクル不安定
– 時間精度 （２次以上）
• 磁気モノポール（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




７成分
div B = 0 の
制約があるた
めに自由度
が減る。
(Sod の)衝撃波管問題
Riemann の解
ρ
P
By
セル境界の値は「平均」でない！
t  t


t
t F (x j1/ 2 , t) dt  2 F x j , t  F x j1, 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  R1 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とλｎ は 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   kwk 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 sgnBx 
 f vy   s  y cA sgnBx 
 v    c sgnB  
 v    c sgnB  
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 xx / 2
– 混合によるエントロピー増加
• セル内を一様と仮定した誤差
– 空間１次精度
空間２次精度
• 線形補間ではうまくゆかない
– Godunov の定理
• TVD: total variation dimishing
– 「波の」単調性 (数値振動なし)
空間２次精度 (高次精度化)
• 変数が階段的に変わるとしたことが原因
• 変数が滑らかに変化していることを考慮
• 単純平均はだめ
– 元の木阿弥
– 「風上」性が失われる
• 補間のしかたに工夫が必要
– 面倒な公式が不可避
•
Godunov の定理 線形補間は×
風上性を考慮して外挿
極大・極小で不自然な振動
単調性の保持 (TVD)
total variation diminishing
Minmod 補間
穏健な補間
勾配は少なめに
どの変数を補間するべきか
•
•
•
•
•
原理的には波の振幅(δwk) を補間すべき
流体力学では 原始変数を補間している。
保存量の補間は良くない。
温度より、圧力がまし。
磁気流体では、縦波・横波に分けて補間した
ほうが良い。
– 大振幅アルフヴェン波のテスト (Fukuda & TH)
• 角速度を補間すると良い。（対称軸の付近）
• 座標値やシフトベクトルは補間しない
時間２次精度
• 時刻とともに流束は変化する。
• [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
n1/ 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
方向分離の問題
•
•
•
•
波面の向きによる依存性
斜めの波は２成分に分離
斜めの情報は直接入らない
∇・Bの発生源
Bx,i 1, j ,k  Bx,i 1, j ,k
B 
2x
B
 By,i, j 1,k Bz ,i , j ,k 1  Bz ,i, j ,k ^1
 y,i , j 1,k

2y
2x
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
１ステップ弱で隣に
モノポールを移動。
数ステップで消滅
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.
を独立変数とした１次元の非線形方程式
この式を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
計算コストの削減
• 実効的な分解能の２倍低下を、細かい格子で
補うにと１６倍の計算量増加
• 一般には高次精度のほうが得
• コストを決めているのは２割の要素
– パレトの法則 （２：８の法則）