U - 千葉大学

簑島敬（海洋研究開発機構）
共同研究者：松本洋介（千葉大学）, 三好隆博（広島大学）
MHDシミュレーションの基本（昨日の復習）
Fi +1/ 2 − Fi −1/ 2
∂U ∂F
n
n +1
.
= 0, ⇒ U i = U i − ∆t
+
∆x
∂t ∂x
Fi +1/ 2
U: 保存量ベクトル（密度、運動量、磁場、エネルギー）
F: 流束ベクトル
Fi +1/ 2 = Riemann(U i , U i +1 ).
（↑Roe法, HLLD法など）
U i U i +1
 今日の話は……実用化
 高次精度化
∂U ∂F ∂G ∂H
 多次元化
+
+
+
=0
∂t
∂x
∂y
∂z
2
高次精度化とは
∂f
∂f
+v
= 0, ( v = const )
∂t
∂x
 物理量プロファイルを高
次関数で近似し、高解像
1次
精度度の解を得る
f ( ∆x ) ~ f (0) + f ' (0)∆x , ←1次
3次
5次
f ( ∆x ) ~ f (0) + f ' (0)∆x
+ f ' ' (0)∆x 2 / 2 + f ' ' ' (0)∆x 3 / 6,
↑3次
厳密解
3
高次精度化の必要性
N次精度スキームの誤差
err ∝ ∆x
∆x
 低次スキーム
で高解像度計
算結果を得る
には、計算負
荷甚大
N
4
高次精度化の必要性
3点補間では実線と点
線の区別がつかない
5点補間では
区別可能
 ステンシルを広げることでプロファイルをより正確に再現
 極値の保持
5
流体シミュレーションの高次精度化例
Fi +1/ 2 − Fi −1/ 2
dU i
∂U ∂F
=−
.
+
= 0, ⇒
dt
∂t ∂x
∆x
高次精度
Fi +1/ 2 = Riemann(U iL+1/ 2 , U iR+1/ 2 ).
例：線形3次精度
U: 保存量ベクトル
F: 流束ベクトル
Fi +1/ 2
U i −1 U i U i +1 U i + 2
U iL+1/ 2 = (− U i −1 + 5U i + 2U i +1 )/ 6
U iR+1/ 2 = (− U i + 2 + 5U i +1 + 2U i )/ 6
線形N次精度化は、流体シミュレーションでは
あまり用いられない
6
高次精度化の欠点
∂f
∂f
+v
= 0, ( v = const )
∂t
∂x
1次
精度
3次
5次
 不連続面近傍の
数値振動
 「線形2次以上ス
キームは解の単
調性を維持できな
い」（Godunovの
定理）
 非線形スキーム
の導入
厳密解
7
非線形スキーム
 低精度@不連続面付近（衝撃波）
 高精度@滑らかな領域（乱流）
 Monotone Upstream-centered Schemes for Conservation
Laws（MUSCL; van Leer 1979）
 Weighted Essentially Non-Oscillatory scheme（WENO;
Jiang+ 1996）
 Monotonicity Preserving scheme（MP; Suresh+ 1997,
CANS+に実装）
http://www.astro.phys.s.chiba-u.ac.jp/cans/doc/riemann.html#id17
8
ui-1
MUSCL法
 セル平均値（既知）
ui
ui+1
x
i-1/2 i+1/2
1 xi+1/ 2
ui =
u( x )dx ,
∫
x
∆x i−1/ 2
 テーラー展開に基づくセルi内の分布
u( x ) = u( x i ) + ui ' ( x i )(x − x i ) + ui ' ' ( x i )(x − x i ) / 2 + 
2
x − xi
∆x 2 
3κ 
2
(ui +1 − ui −1 ) + 2 (x − x i ) − (ui +1 − 2ui + ui −1 ),
≈ ui +
2∆x
2∆x 
12 
 セル境界値
1− κ
(ui − ui −1 ) + 1 + κ (ui +1 − ui ),
4
4
1− κ
(ui +1 − ui ) − 1 + κ (ui − ui −1 ),
= u( x i −1/ 2 ) = ui −
4
4
uiL+1/ 2 = u( x i +1/ 2 ) = ui +
uiR−1/ 2
9
MUSCL法
 流束制限関数Φの導入
 極値の発生を禁止
ui-1
ui
ui+1
x
i-1/2 i+1/2
ui −1 < uiR−1/ 2 < ui < uiL+1/ 2 < ui +1
1− κ
1+ κ
= ui +
Φ(r )(ui − ui −1 ) +
Φ(1/ r )(ui +1 − ui ),
u
4
4
1− κ
1+ κ
R
Φ(1/ r )(ui +1 − ui ) −
Φ(r )(ui − ui −1 ),
ui −1/ 2 = ui −
4
4
r = (ui +1 − ui )/(ui − ui −1 ),
L
i +1/ 2
 代表的な場合： Φ(r)=r Φ(1/r)
uiL+1/ 2 = ui + Φ(r )(ui − ui −1 )/ 2,
uiR−1/ 2 = ui − Φ(1/ r )(ui +1 − ui )/ 2.
10
ui-1
MUSCL法
 流束制限関数
 minmod
Φ(r ) = max[0, min(1, r )],
uiL+1/ 2 = ui + du / 2, uiR−1/ 2 = ui − du / 2,
du = MM (ui +1 − ui , ui − ui −1 ),
ui
ui+1
x
i-1/2 i+1/2
du: セルi内の勾配の
候補値の内、最も滑
らかなものorゼロ
MM ( x , y) = sgn( x ) max[0, min(| x |, sgn( x )y)].
sgn(x): xの符号
 monotonized central
Φ(r ) = max[0, min( 2r ,0.5(1 + r ),2)],
uiL+1/ 2 = ui + du / 2, uiR−1/ 2 = ui − du / 2,
du = MC ( 2(ui +1 − ui ),2(ui − ui −1 ),(ui +1 − ui −1 ) / 2),
MC ( x , y, z ) = MM ( x , MM ( y, z )).
11
U iL+1/ 2
MP法（CANS+に実装）
U i − 2 U i −1 U i U i +1 U i + 2
x
 線形5点補完
U iL+1/ 2 = (2U i − 2 − 13U i −1 + 47U i + 27U i +1 − 3U i + 2 )/ 60,
 このままでは数値振動発生
 周囲の幾つかの候補値と比較し、中央値を取る
U i , U i +1 , UUL , U MD , U LC
 単調性維持（数値振動の抑制）
 極値の保持（MUSCL等のTVDスキームの改良）
12
MP法（CANS+に実装）
Upper limit
UUL = U i + α (U i − U i −1 ),α ≥ 2
Large curvature
Median
U LC = U i + (U i − U i −1 )/ 2 + 4d iMM
−1/ 2 / 3,
d iMM
−1/ 2 = MM (U i + U i − 2 − 2U i −1 ,
U MD = median[(U i + U i +1 ) / 2,
U i + (U i −U i −1 ) / 2, U i +1 − (U i + 2 −U i +1 ) / 2]
U i +1 + U i −1 − 2U i ).
13
MP法（CANS+に実装）
U
Umax
Umin
 単調性維持かつ極値の保持
U iL+1/ 2 ∈ I [U i , UUL , U LC ] and I [U i , U i +1 , U MD ] ,
U iL+1/ 2 = (2U i − 2 − 13U i −1 + 47U i + 27U i +1 − 3U i + 2 )/ 60,
[
= max[min(U , U
]
)],
U max = min max(U i , U i +1 , U MD ), max(U i , UUL , U LC ) ,
U min
i
MD
UL
LC
,
),
min(
,
,
U
U
U
U
i +1
i
U iL+1/ 2 ← median(U iL+1/ 2 , U max , U min ).
14
時間更新
 空間精度を上げた際は、時間精度も上げないと数値的
に不安定
dy
= f (y)
 ルンゲ・クッタ法:
dt
1次
3次（MPで使用）
y n +1 = y n + f ( y n )∆t .
y * = y n + f ( y n )∆t ,
3 n 1 *
**
2次（MUSCLで使用）
y = y + y + f ( y * )∆t ,
4
4
y * = y n + f ( y n )∆t ,
n
2 **
y
n
n +1
**
y
1
=
+
+
y
y
f
(
y
)∆t .
n +1
*
*
y =
+ y + f ( y )∆t .
3 3
2 2
[
[
]
]
[
]
15
線形移流
∂f
∂f
+v
= 0, ( v = const )
∂t
∂x
1次
MUSCLminmod
MP
厳密解
16
1次元理想MHD方程式（保存形式）
∂U ∂F
+
= 0,
∂t ∂x
U = (ρ , ρu, ρv, ρw, By , Bz , e ) , …保存変数（conserved variables）
T

 ρu


2
 ρuu − Bx Bx + p + B / 2 


ρ
vu
B
B
−
x y


,
F =  ρwu − Bx Bz



 By u − B x v



 Bz u − Bx w
 e + p + B 2 / 2 u − B (u ⋅ B )
x


(
(
)
p = (γ − 1) e − ρu 2 / 2 − B 2 / 2 ,
∇ ⋅ B = 0 ⇒ Bx = const,
)
V = (ρ , u, v, w, By , Bz , P ) .
T
…基本変数（primitive variables）
17
1次元MHDシミュレーション
Fi +1/ 2 − Fi −1/ 2
dU i
∂U ∂F
=−
.
+
= 0, ⇒
dt
∂t ∂x
∆x
Fi +1/ 2 = Riemann(U iL+1/ 2 , U iR+1/ 2 ).
 1次精度
U iL+1/ 2 = U i ,
 2次精度
U: 保存量ベクトル
F: 流束ベクトル
Fi +1/ 2
U i − 2 U i −1 U i U i +1 U i + 2
U iL+1/ 2 = MUSCL (U i −1 , U i , U i +1 ),
 5次精度
U iL+1/ 2 = MP(U i − 2 , U i −1 , U i , U i +1 , U i + 2 ).
18
衝撃波管問題（Brio+ 1988）
1次
2次
5次
密度
縦速度
横磁場
HLLD近似リーマン解法使用（Miyoshi+ 2005）
19
双曲型保存則の性質
∂U ∂F
+
= 0,
∂t ∂x
∂U
∂U ∂V
∂V
∂F 

⇒
+A
=
+ Ap
= 0,  A =
,
U
∂t
∂x
∂t
∂x
∂


∂U
∂U
∂V
∂V
⇒L
+ LARL
= Lp
+ L p A pR pL p
= 0, (LR = L p R p = 1),
∂t
∂x
∂t
∂x
∂W
∂W
⇒
+Λ
= 0, (dW = LdU = L p dV, Λ = LAR ). ←非線形移流
∂t
∂x
方程式
 方程式の対角化
 L,R: Aの左,右固有ベクトルを成分とする左,右固有行列
 Λ: 実数固有値を成分とする対角行列
 特性量dWは特性線Λに沿って一定（リーマン不変量）
 風上法に基づく補間は特性量に対して行いたい
20
特性変数変換
 左右固有行列Li,Ri(e.g., Ryu+ 1995)
 特性変数変換：Wi+s=LiUi+s, (s=-2,…,2 for MP)
L
R
 Wi +1/ 2 , Wi −1/ 2 , を補間で求める
 逆変換：U iL+1/ 2 = R i WiL+1/ 2 , U iR−1/ 2 = R i WiR−1/ 2
 固有行列を計算するのが面倒
 プログラムミス多発!!
 固有行列が明らかでない方程式系では使えない
 CANS+に実装済み
21
衝撃波管問題（Brio+ 1988）
1次
2次
5次
密度
縦速度
横磁場
HLLD, 特性変数変換あり
22
2次元理想MHD方程式（保存形式）
∂U ∂F ∂G
+
+
= 0,
∂t ∂x ∂y
U = (ρ , ρu, ρv, ρw, Bx , By , Bz , e ) ,
T

 ρv

 ρu




2
 ρuv − By Bx

 ρuu − Bx Bx + p + B / 2 


2


 ρvv − By By + p + B / 2 

 ρvu − Bx By
 ρwv − B B

 ρwu − B B

y
z
x z
.
, G = 
F=
 B x v − By u


0


B u− B v

x
0

 y

B v− B w

 Bz u − Bx w

y
 z



 e + p + B 2 / 2 u − B (u ⋅ B )
 e + p + B 2 / 2 v − B (u ⋅ B )
x


y


(
)
(
)
23
MHDコードの多次元化
U: 保存量ベクトル
F: X方向の流束ベクトル
G: Y方向の流束ベクトル
Fi +1/ 2, j − Fi −1/ 2, j G i , j +1/ 2 − G i , j −1/ 2
dU i , j
∂U ∂F ∂G
+
+
= 0, ⇒
=−
−
.
∂t ∂x ∂y
∆x
∆y
dt
Gi,j+1/2
 1次元数値解法の利用
Fi +1/ 2, j = Riemann(U iL+1/ 2, j , U iR+1/ 2, j ),
G i , j +1/ 2 = Riemann(U iL, j +1/ 2 , U iR, j +1/ 2 ),
 同じタイミングで求めるUnsplit法
y
(i,j)
Fi+1/2,j
x
 磁場発散の処理
 数値的なdivB≠0で計算破綻
三好、日本物理学会誌（2014）
24
磁場発散の処理
 Constrained Transport (CT) 法（Evans+ 1988）
 DivB=0を離散化レベルで保証する特殊なグリッド配置
 プロジェクション法（Brackbill+ 1980）
 ポアソン方程式を解いて磁場修正
 移流拡散法（Powell 1999, Dedner+ 2002）
 有限のdivBを移流拡散方程式や電信方程式に従って処理
http://www.astro.phys.s.chiba-u.ac.jp/cans/doc/riemann.html#id23
25
プロジェクション法/CT法
 プロジェクション法
B n +1 = B* + ∇φ , ⇒ ∆φ = −∇ ⋅ B* .
 B*: MHDシミュレーションで得た解（divB誤差含む）
 Bn+1: divB=0を満たす真の解
 ポアソン方程式解いてφ得る => B*をBn+1へ更新
 CT法
 セル境界面の中心に面に
垂直な磁場成分を配置
 セル角の電場を求める
（e.g., Balsara+ 1999）
 ∂B
= −∇ × E.
∂t
Bx(i-1/2,j)
Bz(i,j)
Ez(i-1/2,j-1/2) By(i,j-1/2)
26
移流拡散法（CANS+に実装）
 9-wave法（Dedner+ 2002）
 divBを掃除するため方程式を変形、追加
2
c
∂B
∂φ
+ ∇ × E + ∇φ = 0,
+ ch2∇ ⋅ B = − h2 φ
∂t
∂t
cp
 Bとφについて整理すると…
 ∂ 2 ch2 ∂ 
2 2
(
)
c
+
∇
⋅
B
=
 2

h ∇ (∇ ⋅ B ),
2
c p ∂t 
 ∂t
電信方程式（波動+拡散方程式）
 実装しやすく計算負荷も少ない
27
移流拡散法
 1次元の場合
∂  Bx   0
  +  2
∂t  φ   ch

1  ∂  Bx   0
,
   = 
2
2
0  ∂x  φ   − (ch / c p )φ 
演算子分離
∂  Bx   0
  +  2
∂t  φ   ch
1  ∂  Bx 
   = 0,
0  ∂x  φ 
固有値（伝播速度） ±ch
ch: CFL条件から決める最大値
cp: (0.18ch)(1/2) （経験的な値）
∂φ
= −(ch2 / c p2 )φ ,
∂t
φ * = exp[−(ch2 / c p2 )∆t ]φ n ,
28
移流拡散法
∂  Bx   0
  +  2
∂t  φ   ch
1  ∂  Bx 
   = 0,
0  ∂x  φ 
 Lax-Friedrichs flux splitting
(
)
(
)
ch R
1 L
R
φi +1/ 2 = φi +1/ 2 + φi +1/ 2 −
Bx ,i +1/ 2 − BxL,i +1/ 2 ,
2
2
2
c
ch R
L
R
2
h
φi +1/ 2 − φiL+1/ 2 .
ch Bx ,i +1/ 2 =
Bx ,i +1/ 2 + Bx ,i +1/ 2 −
2
2
(
)
(
)
29
2D MHD with 9-wave
 Orszag-Tang 渦問題
温度
|divB|
http://www.astro.phys.s.chiba-u.ac.jp/cans/doc/riemann.html#id23
30
高次精度MHDコードの実用例
 ケルビン・ヘルムホルツ不安定
 磁気リコネクション
 1次, 2次, 5次精度で特性変数補間
 HLLD近似リーマン解法
 3次ルンゲ・クッタ法
31
ケルビン・ヘルムホルツ不安定
 速度差がある流体の間で発生する不安定
 2次的不安定で乱流発展（Matsumoto+2006）
1次
2次
5次
32
KH乱流発展
1次
2次
5次
密度層幅λ2
速度のパワースペクトル
k-2
時間
dλ
dt
= 0.26 (5th)
= 0.15 (2nd)
= 0.02 (1st)
 乱流混合 Dyy =
2
k-3
1次
2次
5次
k-7
波数
33
磁気リコネクション
 宇宙プラズマにおける磁
気エネルギー解放機構
 恒星フレア、磁気圏サブ
ストーム、パルサー磁気
圏…
5次
2次
磁場
電流
1次
流れ
34
磁気リコネクション
1次
2次
5次
磁気
熱
運動x10
 エネルギー時間変化
 1次スキームでは数値散
逸が大きく、MRX失敗
 5次スキームではより効
率的なエネルギー変換
 電気抵抗
 1,2次スキームでは、数値
的な抵抗で散逸（実線 >>
点線）
 5次スキームでは手で与
えた抵抗で散逸（実線 ~
点線）
35
まとめ
 実用的なMHDシミュレーションのためには、高次精度化
と多次元化は必須!!
 高次精度化
 セル境界の左右の値を非線形補間で求める
 次数が上がるほど特性変数変換が重要
 多次元化
 1次元数値解法の利用
 磁場発散の処理（CT, プロジェクション、移流拡散法）
36

Download Report