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

MHD数値シミュレーションの基礎
•
•
•
•
ニュートン流体力学方程式
（ニュートン）磁気流体力学方程式
相対論的流体力学方程式
相対論的磁気流体力学方程式
• 千葉大学・花輪知幸
一部のスライドは高橋博之さんからの借用
保存形式
U Fx Fy Fz



S
t
x
y
z
Fx  Fx V 
U  UV 

 

v 

 x
V  v y  , U  

v 

 z
  
  


vx
v y
vz
 v
S  SV, r, t 
2
vx






2

v

P
x



 , Fx  

vx v y



vx vz



2
2 
    v 2  P vx 


 
重力、遠心力、加熱・冷却
T
相対論的磁気流体力学
ー保存形式ー

g

n u
F



,
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 
t
  F dS dt
0
保存量と流束
保存量はセル平均の値
Fi 1/ 2, j ,k
U i , j ,k
U i , j ,k 
Fi 1/ 2, j ,k 
U
dV

Vi , j ,k
F

d
S
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 j 3
x


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
 v x2  P 


  


8
4
 v 


Bx B y
x


v x v y 


4

 v y 


B
B

U   v z  , F   v x v z  x z


4
 By 


v
B

v
B
x
y
y
x




B
z


v x Bz  v z 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
s 
f 
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
F
F

UR 
U

U6 
U
U1 U2 U3
U4
U 6    f  U R  U 6 
U 5    A U 6  U 5 
U5 U6 UR
衝撃波の
Jump
Condition
..............
Unとλｎは UR
F U  U    U  U 
L
f
1
L
と ULの関数
U 1
HLL
FR
UL
F
UM
HLL
FL
UR
Lt  0
保存則
F
HLL
 FL 
F HLL  FR 
F
HLL
R t  0
L U M  U L 
R U M  U R 
R FL  L FR  R L U R  U L 

R L
λL とλR の
とり方に多
様性
特性速度
速い磁気音波  f   v x  c f
アルフヴェン波 A   v x  cA
s   v x  c s
エントロピー波 0  v x
遅い磁気音波
2
x
B
c 
,
4
2
A
P
P
c 
 , cA2  
 s 
2
s


Bx2  By2  Bz2
4
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
vy
 L v yL   L v yR
L  R
, By 
UL
UR
 R ByL   L ByR
Q  QR  QL ,  
L  R
vx    vx  vx 
R L
保存量 U は W の自乗で表せる

W ,

 vx ,  v y ,  vz ,
P

,
Bz 
,


 
By
固有モード
U 
7
 w
k
k 1
F 
rk : 固有ベクトル
7
  w
k 1
f
k : 特性速度
rk
k
A
k
rk
s
固有モード
c
s
A
f
x
固有モード
P, vx , v// , B//
• 縦波 f, s
• 横波 A
• エントロピー波 c
f
A
s
v , B

c
s
A
f
x
横波とエントロピー波は区別しやすい
 z 
 By 
y 

   B//    B 


B

 z
 z
 y 
 z 
 vy 
y 

   v//    v 


v

y
 z
 z


 By 
y 
1
 
  
B y2  Bz2  Bz 
 z 
  
  c P    s
 s  P
2
s
縦波: Fast & Slow
rf 
s
f












c

v

c

v

s
x
s
f
x
f




 s v y   f  y cs sgn Bx 
 f v y   s  y c A sgn Bx 
 v    c sgn B  
 v    c sgn B  
x
z s
f
x
z A
s
f z

 s z


s  y c f
,
 , rs      f  y c f




4 
4

   c


s z cf
z f
f




4 
4



2
2
2




v

v

v
vx2  v y2  v z2
z
y
x


 f
f
 gs 
 gf 
2
2




f 
c 2f  c A2
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


H 

P
 1 

v
2
2
速度差は実効的音
速を増やす。
CA'* も平均の磁場
からけいさんされる
より大きい
Cargo & Gallice (1998)
一般の状態方程式 Nobuta & TH, Mikami et al.
HLL 系では付加条件が必要
• 波の分解が不十分なので、整合条件だけで
は足りない。
• 付加的な条件が必要
• 特性速度を大きめにとると、密度・圧力が正
に保たれる (cf. Miyoshi et al.)
どこで特性速度を見積もるか
k , L , k , k , R
• 絶対値が大きいものが良い。
• 符号が同じなら平均で十分（とくに滑らか
な変化の場合は平均が良い）
• 符号が変わる場合は、平均より絶対値を
大きめにとる（エントロピーフィックス）
特性速度と数値粘性
F
ROE
F
HLL
1
 FL  FR 
2

k wk rk 

k 1

7


1
R  L 
R  L
FR  FL  
U R  U L  
 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の発生源
B 

Bx ,i 1, j ,k  Bx ,i 1, j ,k
By ,i , j 1,k
2x
 By ,i , j 1,k
2y

Bz ,i , j ,k 1  Bz ,i , j ,k ^1
2x
div B の消去
Hyperbolic Divergence Cleaning
Dedner et al., JCP, 175, 645 (2002)

   v   0
t
2


B


v      vv   P 
t
8





BB
I 
  0

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 


  0
v      vv   P 
I

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  B 0  B1
  B 0  0 のとき有効
強い双極子磁場をもつ場合
球座標で一様磁場をとく場合
Tanaka
計算コストの削減
• 実効的な分解能の２倍低下を、細かい格子で
補うにと１６倍の計算量増加
• 一般には高次精度のほうが得
• コストを決めているのは２割の要素
– パレトの法則（２：８の法則）

Download Report