近似リーマン解法による磁気流体方程式の差分解法

2015年8月3日（月）-7日（金）千葉大学統合情報センター
宇宙磁気流体・プラズマシミュレーションサマースクール
近似リーマン解法による
磁気流体方程式の差分解法
三好隆博
広島大学大学院理学研究科
目的
風上？
それどこの風？
目的
開始
初期設定
CFL条件の評価
SSP-RK法
高次精度補間
数値流束の評価
時間積分
t = tout
Yes
データ出力
Yes
t < tend
終了
No
No
目的
簑島先生の
講義では、
開始
初期設定
CFL条件の評価
SSP-RK法
高次精度補間
数値流束の評価
時間積分
t = tout
No
Yes
データ出力
Yes
t < tend
終了
No
＋多次元化
おぉぉ、相当大変・・・
内容
 はじめに
 MHD方程式
 風上型解法の基礎
 スカラー方程式の解法
 システム方程式の解法
 MHD方程式の近似リーマン解法
 線形近似リーマン解法
 HLL型近似リーマン解法
 フローチャート
2015年8月3日（月）-7日（金）千葉大学統合情報センター
宇宙磁気流体・プラズマシミュレーションサマースクール
はじめに
流体
 身近な流れは非圧縮性流体として非常によく近似
 縮まない（密度変動が少ない）流れ
 低マッハ数（流速/音速 < 0.3）流れ
 ただし、音響学分野では音波考慮
 非日常的で極限的な流れでは圧縮性が本質的に重要
 高マッハ数（流速/音速 > 1）流れ
 衝撃波
地上での衝撃波
Vapor cone (F-14)
Schlieren image (.30-06 Springfield)
Afterburner (SR-71)
Meteor (Chelyabinsk)
Volcanic eruption (Sarychev Peak：芙蓉山)
宇宙流体
円谷プロ
野望
 宇宙流体は圧縮性流体としての取り扱いが不可欠
 密度・圧力変動が大きい流れ
 高マッハ数（流速/音速 >> 1）流れ
 磁場、放射などが流体と非線形相互作用
 高エネルギー天体周りでは相対論的効果
偏在する磁気流体力学（MHD）的な不連続！
本講義では、不連続解をも恐れぬ
数値磁気流体力学の猛者を養成。
2015年8月3日（月）-7日（金）千葉大学統合情報センター
宇宙磁気流体・プラズマシミュレーションサマースクール
MHD方程式
MHD方程式
 Euler方程式＋ローレンツ力＋磁場の誘導方程式

   v   0
t
 v

   v  v   p  J  B
 t

 p 
 p 
    v      0
t   
 
B
   v  B   0 ,   B  0
t

連続の式

運動方程式

断熱の式

磁場の誘導方程式
 流体と磁場の相互作用を支配
速進磁気音波 > アルフェン波 > 遅進磁気音波
MHD方程式
 保存形式


   v   0
t
 v 
   vv  pT I  BB   0 
t
e
   e  pT v  B v  B   0 
t
B

   vB  Bv   0
t

v 2 B 2 
 ,
  B  0 , p    1 e 

2
2 

質量保存則
運動量保存則
エネルギー保存則
磁束保存則
B2
pT  p 
2
MHD方程式
 1次元MHD方程式
U F

 0 , Bx  const. ,
t
x


u
  


 
2
uu  pT  Bx


 u 


 v 
 vu  B x B y


 

U   w  , F  
wu  Bx Bz


B 
B y u  Bx v


 y


 Bz 
Bz u  Bx w


 
 e 
 e  pT u  Bx uB x  vB y  wB z 
MHD方程式
 1次元MHD方程式
U F U
U
F


A
0, A
U
t
x
t
x
 Aの固有値（特性速度）
AR  RΛ , Λ  diag 1 , 2 ,  , 7 
1  u  c f  2  u  ca  3  u  cs  4  u
 5  u  cs  6  u  ca  7  u  c f
bB
c
2
f ,s
 , ca2  bx2 , a 2  p  ,
1 2
  a  b2 
2
a  b 
2 2
 4a 2bx2 

MHD方程式
 MHD方程式の不連続解
 U  F 
  t  x dxdt  x U L  U R   t FR  FL   0
S  x t
t
 S U   F  , X   X R  X L
不連続にのった系で、
t
UL
F   0, u   0 （速進/遅進衝撃波）
    p   B y2  Bz2   0 ,   v   B y , 
v   w  B y   Bz    p   0 , Bx  0
p  B
2
y

UR
x
x
 w  Bz 
（回転不連続）
（接触不連続）
 Bz2  2  0 , Bx  0 （接線不連続）
2015年8月3日（月）-7日（金）千葉大学統合情報センター
宇宙磁気流体・プラズマシミュレーションサマースクール
風上型解法の基礎
話を進めるその前に
 時間・空間座標、変数の離散表現
t
i, n  1
n2
t
n 1
i  1, n 
n
n 1
t
n2
i  2 i 1 i
i 1 i  2
x
i  1, n 
i, n 
i, n  1
x
x

xi , t n , uin  u xi , t n

流れの基礎方程式
 双曲型保存則
U F U
U
F


A
 0, A 
U
x
t
x
t
Aが独立な実固有値をもつとき双曲型
 線形移流方程式
 非粘性Burgers方程式
 Euler方程式
 MHD方程式
など
スカラー方程式
システム方程式
スカラー方程式の解法
 線形移流方程式
u
u
a
 0 , a  const.
x
t
u
a
 u  x, t   u  x  at ,0 
x
 風上差分法
uin 1  uin
uin  uin1
a
 0, a  0
t
x
スカラー方程式の解法
 線形移流方程式の保存型解法
u f

 0 , f  au
t x
u
n 1
i
f
u

t
n
i

i 1 / 2
f
x

i 1 / 2
f i 1/ 2
f i 1/ 2
ui 1
ui
ui 1
xi 1/ 2
0
xi 1/ 2
f i 1/ 2 ：数値流束
 流束ベクトル分離法（FVS法）



f

f
f  f  f,
 0,
0
u
u
 f i 1/ 2  f i   f i 1
f i 1 / 2
u
fi 

f i 1
uin
uin1
a | a |
a n
|a| n

n
f 
u  f i 1/ 2  ui 1  ui 
ui 1  uin
2
2
2





x
スカラー方程式の解法
 線形移流方程式の保存型解法
f i 1/ 2
u f

 0 , f  au
t x
u
n 1
i
u
f

t
n
i

i 1 / 2
f
x

i 1 / 2
ui 1
ui
xi 1/ 2
0
f i 1/ 2
ui 1
xi 1/ 2
f i 1/ 2 ：数値流束
 流束ベクトル分離法（FVS法）



f

f
f  f  f,
 0,
0
u
u
 f i 1/ 2  f i   f i 1
Lax法
a  t x
a n
t n

n
f 
u  f i 1/ 2  ui 1  ui 
ui 1  uin
2
2
2 x





スカラー方程式の解法
 線形移流方程式の保存型解法
 リーマン解法（Godunov法）
リーマン問題＝衝撃波管問題を利用
u
リーマン解法の手順：
1. 区分定数分布を仮定
uin1 uin
uin1
x
スカラー方程式の解法
 線形移流方程式の保存型解法
 リーマン解法（Godunov法）
リーマン問題＝衝撃波管問題を利用
at
u
リーマン解法の手順：
at
1. 区分定数分布を仮定
2. 移流方程式の厳密解
at
x
スカラー方程式の解法
 線形移流方程式の保存型解法
 リーマン解法（Godunov法）
リーマン問題＝衝撃波管問題を利用
u
リーマン解法の手順：
uin11 uin 1 uin11
x
1. 区分定数分布を仮定
2. 移流方程式の厳密解
3. 厳密解の空間平均値
スカラー方程式の解法
 線形移流方程式の保存型解法
 リーマン解法（Godunov法）
リーマン問題＝衝撃波管問題を利用
u
リーマン解法の手順：
uin11 uin 1 uin11
x
1. 区分定数分布を仮定
2. 移流方程式の厳密解
3. 厳密解の空間平均値
数値流束
 u f 
  t  x dxdt  0
 at uin  uin1  t f i n1  f i 1/ 2  0 , a  0




at
t
f i 1/ 2
uin
xi 1/ 2
uin1
f i n1
x
スカラー方程式の解法
 線形移流方程式の解法
 風上差分法
t n

ui  uin1 , a  0
u u a
x
 流束ベクトル分離法
n 1
i
u
n 1
i
n
i
t 
u 
f i 1/ 2  f i 1/ 2 , f i 1/ 2  au in , a  0
x
n
i


 リーマン解法
u
n 1
i
a t n x  a t n

ui 1 
ui , a  0
x
x
全て同じ解法に帰着（でも、思想は異なる）
スカラー方程式の解法
 非線形移流方程式
u f
u
u
f

0 
 a u   0, a u  
t x
t
x
u
非保存形式
保存形式
dx
 a u  に沿って du  0  u x, t   u  x  a u t ,0 
dt
x  x0  u  x 0 , 0  t
特性曲線
t
t
非線形移流方程式
x
線形移流方程式
x
スカラー方程式の解法
 非粘性Burgers方程式の解法
u
u
u   u 2 
u
u
    0 
u
0
t x  2 
t
x
 風上差分法（非保存形式）
n
n
uin 1  uin
u
u

i 1
 uin i
 0 for uin  0
t
x
n
n
u
u
uin 1  uin

 uin i 1 i  0 for uin  0
t
x
n
n
n

uin 1  uin
u
u
|
u
n i 1
i 1
i |


 ui

uin1  2uin  uin1 
t
2 x
2 x
x
スカラー方程式の解法
 非粘性Burgers方程式の解法
u   u 2 
u
u
    0 
u
0
t x  2 
t
x
 流速ベクトル分離法
uin 1  uin f i 1/ 2  f i 1/ 2

0
t
x
u u  | u |




f i 1/ 2  f i  f i 1 , f 
4
n
n
f

f
1 n

i 1
i
 f i 1/ 2 
 | ui 1 | uin1  | uin | uin
2
4


スカラー方程式の解法
 非粘性Burgers方程式の解法
 衝撃波管問題
衝撃波 t
S
x
0  uin1  uin
t
S t
S
uin1  uin  0
x
S t
S
f i 1  f i
ui 1  ui
ui 1  ui

2
uin1  0  uin
x
膨張波 t
0  uin  uin1
uin1  0  uin
x
t
x
uin  0  uin1
t
x
uin  uin1  0
x
スカラー方程式の解法
 非粘性Burgers方程式の解法
u   u 2 
u
u
    0 
u
0
t x  2 
t
x
 リーマン解法
uin 1  uin f i 1/ 2  f i 1/ 2

0
t
x
n
n
 un 2 2
if
u
0
,
u

i
i
i 1  0
 n 2
n
n
u
2
if
u
0
,
u

 i 1
i
i 1  0


f i 1/ 2  0
if uin  0  uin1
 n 2
n
n
u
2
if
u
0
u


i
i 1 , S i 1 / 2  0
 i 2
 uin1 2 if uin  0  uin1 , S i 1/ 2  0
 
 
 
 
スカラー方程式の解法
 非粘性Burgers方程式の解法
u   u 2 
u
u
    0 
u
0
t x  2 
t
x
 線形近似リーマン解法（Roe法）
局所的に線形化
⇒ 膨張波を無視
uin 1  uin f i 1 / 2  f i 1 / 2

0
t
x
n
n

|
f

f
a

i 1 / 2 |
i 1
i 1

 f i 1 / 2 

uin1  uin 
2
2
n
n
n
n
f

f
u

u
i 1
i
i
ai1 / 2  in1

2
ui 1  uin
ai1 / 2 t
t
f i 1/ 2
uin
xi 1/ 2
uin1
f i n1
x
スカラー方程式の解法
 非粘性Burgers方程式の解法
非保存型解法
1.5 for x  0.1
u  x ,0   
0.5 for x  0.1
保存型解法
スカラー方程式の解法
 非粘性Burgers方程式の解法
非保存型解法
0.5 for x  0.1
u  x ,0   
1.5 for x  0.1
保存型解法
スカラー方程式の解法
 保存型解法
 Lax-Wendroffの定理[1960]
数値解が収束すれば、その解は保存則の弱
解に収束する
 Hartenのエントロピー条件[1980]
数値解がエントロピー条件を満足し、収束す
れば、その解は保存則の物理解に収束する
 非保存型解法
 Hou-LeFlochの定理[1994]
数値解が収束したとしても、衝撃波を含むそ
の解は非物理解に収束する
スカラー方程式の解法
 非粘性Burgers方程式の解法
FVS法/Godunov法
Roe法
  0.5 for x  0.5
u  x ,0   
 0.5 for x  0.5
局所線形化により膨張衝撃波
⇒エントロピー補正
スカラー方程式の解法
 非粘性Burgers方程式の解法
FVS法/Godunov法
Roe法＋エントロピー補正
様々な補正法があるが、ここでは、

n
n
n
n



|
a
|
a
a
2
if
a
0
a




i 1
i
i
i 1
| ai1 / 2 |  i1 / 2
otherwise
| ai 1 / 2 |
スカラー方程式の解法
 Nonconvexな非線形移流方程式
u f
u3

 0, f 
t x
3
a u   u 2
uL2  uLuR  uR2
S
3
aL  aR  S
混合波
（衝撃波+膨張波）
aL  S  aR
aL  S  aR
t
t
x
t
x
S
x
システム方程式の解法
 双曲型保存則系
U F

0
t
x
U ：保存変数ベクトル F ：流束ベクトル
U
F
U
A
 0 , AU  
t
x
U
W
W
1 U
1
1 U
R
 R ARR

Λ
 0 , AR  ΛR
t
x
t
x
A：流束ヤコビアン行列 R：右固有行列 Λ：固有値行列
独立な実固有値・固有ベクトルをもつとき
⇒ 双曲型
システム方程式の解法
 双曲型保存則系
 線形2×2保存則系
 v 
u
U F
 0 , U    , F   2  , a  const.  0

x
t
a u
v
 u  0
    2
t  v   a
1
R  
a
1  u 
    0
0  x  v 
1 
 ,
 a
 w1 
1  au  v 
1  u 
   R   


 v  2a  au  v 
 w2 
  w1   a 0    w1 
   

    0
t  w2   0  a  x  w2 
システム方程式の解法
 双曲型保存則系
 等温Euler方程式
u 
 

U F
 0 , U    , F   2

 , a  const.  0
2
x
t
 u 
 u  a 
    u
    2
t  u   a 

R  
a
    
    0
u  x  u 
   w1 
  1 a  u
 ,    R   


 a   w2 
 u  2a  a  u 
1
0    w1 
  w1   u  a
   

    0
u  a  x  w2 
t  w2   0
システム方程式の解法
 双曲型保存則系
W
W
 0 , dW  R 1dU
Λ
x
t
W ：特性変数ベクトル
 w1   1 0  0   w1 
  
  
    w2 
  w2   0 2

0






 0 x 
t 
  
  
 wm   0  0 m   wm 
dx dt  k に沿って dwk  0
連立移流方程式 ⇒ スカラー方程式の解法の拡張
システム方程式の解法
 線形双曲型保存則系の解法
U F U
U

A
 0 , F  AU , A  const.

x
t
x
t
 FVS法
1
1 1
R F  R F  | Λ | R 1U 
2
R 1Fi 1 / 2  R 1Fi   R 1Fi 1
1

1
Fi 1 / 2
1
R U
R 1Fi 1 / 2
R 1Fi 
R 1Fi 1
R 1U i R 1U i 1
x
R Fi 1  R Fi 1

 | Λ | R 1U i 1  | Λ | R 1U i 
2
2
Fi 1  Fi 1

 | A | U i 1  | A | U i  , | A | R | Λ | R 1
2
2
もう時間方向のnはやめますね・・・
システム方程式の解法
 線形双曲型保存則系の解法
U F U
U

A
 0 , F  AU , A  const.

x
t
x
t
 リーマン解法
R 1Fi 1 / 2  R 1Fi 1  Λ R 1U i 1  R 1U i 
1

i 1 / 2
R F
 R Fi  Λ R U i 1  R U i 
1

1
1
ここで Λ   Λ | Λ | 2
1

i 1 / 2
R F

i 1 / 2
F
1 1 
 R Fi 1 / 2  R 1Fi 1 / 2 
2
Λ
t
R 1Fi 1 / 2
R 1U i
R 1U i 1
R 1Fi 1
x
Λ
R 1Fi
t
R 1U i 1
R 1U i
Fi 1  Fi | A |
U i 1  U i  , | A | R | Λ | R 1


2
2
R 1Fi 1 / 2
x
システム方程式の解法
 非線形双曲型保存則系の解法
U F U
U
F

A
 0 , AU  

U
x
t
x
t
Conservative vs Non-conservative
 非保存型解法
たとえ数値解が収束しても、
不連続解を含む問題では、
非物理的な解
“Computational Tutorial: MHD” by G. Toth
http://www.lorentzcenter.nl/lc/web/2011/441/presentations/Advanced_Toth.pdf
システム方程式の解法
 非線形双曲型保存則系の解法
U F U
U
F

A
 0 , AU  

U
x
t
x
t
 FVS法 [e.g., Steger+, 1981; van Leer, 1991]
一般にはFVS法が適用できるとは限らない
でも、たまたまEuler方程式では、
F  AU
なので、

i 1 / 2
F
Fi 1  Fi 1

 | A |i 1 U i 1  | A |i U i 
2
2
でも、MHD方程式ではだめ
システム方程式の解法
 非線形双曲型保存則系の解法
U F U
U
F

A
 0 , AU  

U
x
t
x
t
 FVS法 [e.g., Steger+, 1981; van Leer, 1991]
あ、でも、なんでも使える簡単なFVS法が！

F  F   F  , F   1 2  F   max U
だって、
1

2 R F
なので、

i 1 / 2
F


U R  R
1
A  
max


I R  Λ   max I
local Lax-Friedrichs法
Fi 1  Fi 1

 |  |max i 1 U i 1  |  |max i U i 
2
2
システム方程式の解法
 非線形双曲型保存則系の解法
U F U
U
F

A
 0 , AU  

U
x
t
x
t
 リーマン解法（Godunov法） [Godunov, 1959]
これは中々たいへん
衝撃波管問題の厳密解（衝撃波、膨張波、混合
波の組合せ）を求めるには繰り返し計算が必要
しかも方程式ごとにアルゴリズムを考えねば・・・
という労力にも関わらず、FVS法や近似リーマン
解法と精度は変わらなかったでしょ？
システム方程式の解法
 非線形双曲型保存則系の解法
U F U
U
F

A
 0 , AU  

U
x
t
x
t
 線形近似リーマン解法（Roe法） [Roe, 1981]
局所的に線形化 ⇒ 局所的にAを凍結

i 1 / 2
F
Fi 1  Fi A U i 1 , U i 
U i 1  U i 


2
2
Property U
1． U i , U i 1  U  A U i , U i 1   A
2． Fi 1  Fi  A U i 1 , U i U i 1  U i 
3．実固有値、線形独立な固有ベクトル
システム方程式の解法
 非線形双曲型保存則系の解法
U F U
U
F

A
 0 , AU  

U
x
t
x
t
 HLL型近似リーマン解法 [Harten+, 1983]
HLL法、HLLC法、HLLD法など
このあとすぐ！お楽しみに！
2015年8月3日（月）-7日（金）千葉大学統合情報センター
宇宙磁気流体・プラズマシミュレーションサマースクール
MHD方程式の
近似リーマン解法
近似リーマン解法リターンズ
 近似リーマン解法
 U F 
  t  x dxdt   Udx  Fdt   0
U
 物理量を区分定数分布
x
近似リーマン解法リターンズ
 近似リーマン解法
 U F 
  t  x dxdt   Udx  Fdt   0
U
 物理量を区分定数分布
 衝撃波管問題の近似解
x
近似リーマン解法リターンズ
 MHD衝撃波管問題
t
U  U x t ;U R , U L 
FS / FR RD SS / SR CD SS / SR RD FS / FR
UL
UR
他にも複合波が・・・
x
近似リーマン解法リターンズ
 近似リーマン解法
 U F 
  t  x dxdt   Udx  Fdt   0
U
 物理量を区分定数分布
 衝撃波管問題の近似解
 近似解の空間積分
x
近似リーマン解法リターンズ
 近似リーマン解法
 U F 
  t  x dxdt   Udx  Fdt   0
U
 物理量を区分定数分布
 衝撃波管問題の近似解
 近似解の空間積分
x
 数値流束による形式（時空間保存則から評価）

xi 1 / 2
xi 1 / 2

xi 1 / 2
xi
U n 1dx  
xi 1 / 2
xi 1 / 2
t n 1
t n 1
t
t
U n dx   n Fi 1 / 2dt   n Fi 1 / 2dt  0
 x  xi 1 / 2 n n 
U
;Ui , Ui 1 dx   xi 1 / 2  xi Uin  t Fi 1 / 2  Fi n   0
 t

近似リーマン解法リターンズ
 近似リーマン解法
 U F 
  t  x dxdt   Udx  Fdt   0
 x  xi 1 / 2 n n 
U
;Ui , Ui 1 
 t

t n  t
U
U
n
i 1
xi 1 / 2
U in1
n
i
知ってる！
xi 1 / 2
tn
近似リーマン解法リターンズ
 近似リーマン解法
 U F 
  t  x dxdt   Udx  Fdt   0
 x  xi 1 / 2 n n 
U
;Ui , Ui 1 
xi 1 / 2
n 1
n
Ui   U x, t  t dx
 t

x
i 1 / 2
答えでる！
U
U
n
i 1
xi 1 / 2
U in1
n
i
知ってる！
t n  t
xi 1 / 2
tn
近似リーマン解法リターンズ
 近似リーマン解法
 U F 
  t  x dxdt   Udx  Fdt   0
U
n 1
i

U
t n  t
U in1
n
i
…
…
xi 1 / 2
知らない
知ってる！
知らない
U
xi 1 / 2
n
i 1
U x, t n  t dx
xi 1 / 2
知ってる！
t n  t
Fi 1 / 2   n
t
xi 1 / 2
F  xi 1 / 2 , t dt
tn
近似リーマン解法リターンズ
 近似リーマン解法
 U F 
  t  x dxdt   Udx  Fdt   0
 x  xi 1 / 2 n n 
U
;Ui , Ui 1 
 t

U
答えでる！
答えでる！
U
n
i 1
t n  t
U in1
n
i
xi 1 / 2
xi 1 / 2
t n  t
Fi 1 / 2   n
t
F  xi 1 / 2 , t dt
tn
近似リーマン解法リターンズ
 近似リーマン解法
 U F 
  t  x dxdt   Udx  Fdt   0
 x  xi 1 / 2 n n 
U
;Ui , Ui 1 
 t

U
答えでる！
答えでる！
U
n
i 1
t n  t
U in1
n
i
xi 1 / 2
xi 1 / 2
t n  t
Fi 1 / 2   n
t
F  xi 1 / 2 , t dt
tn
線形近似リーマン解法
 線形近似リーマン解法 [Roe, 1981; Brio+, 1988]
 局所的に線形化（ヤコビアンを凍結）
Fi 1  Fi A U i 1 , U i 
U i 1  U i 
F


2
2
1
|
A
|

R
|
Λ
|
R
ここで

i 1 / 2
 さてと・・・
MHDの固有ベクトルと固有値がわかればよい
A U i 1 , U i  を求めるのは一般には面倒だけど
線形近似リーマン解法
 線形近似リーマン解法 [Roe, 1981; Brio+, 1988]
 MHD方程式の固有ベクトル






1
0










1








u  c f ,s
0






u





Bx B y c f , s 





 v 2
 Bz sgn Bx 
2




v

c
B



f ,s
x





B x Bz c f , s 






 B y sgn Bx 
Ru 
, Ru  ca  
w
 , Ru  c f , s   w  c 2  B 2 


f ,s
x




2


Bz
B y c f ,s




0


2
2





c f , s  B x






2


B
0
Bz c f , s




 y


2
2



 u 2  v 2  w2 
c f , s  B x







 2
u  v 2  w2
  vBz  wB y sgn Bx 
2


 



 h f ,s




2


c 2f ,s
B c vB  wB    2 2


c f ,s  a 2 
h f ,s 
 c f ,s u  x f ,s 2 y 2 z 
固有ベクトルに特異点
 1
c f , s  B x
 1
線形近似リーマン解法
 線形近似リーマン解法 [Roe, 1981; Brio+, 1988]
 固有ベクトルの再規格化
2
2
2
c
a

c
f a
2
2
s
f  2
, s  2
2
c f  cs
c f  cs2
[Roe+, 1996]
 s2   2f  1,  s2 cs2   2f c 2f  a 2 ,  f  s 
by
y 
b b
2
y
2
z
, z 
b  b  0   y ,z
2
y
2
z
bz
b b
2
y
1

2
2
z
a by2  bz2
c 2f  cs2
bx , y , z 
Bx , y , z

線形近似リーマン解法
 線形近似リーマン解法 [Roe, 1981; Brio+, 1988]
 固有ベクトルの再規格化
・・・で色々と計算するわけですが、ここに書いても
あれなので参考文献をご参照ください
[1] Brio, Wu, JCP 75, 400 (1988)
[2] Ryu, Jones, ApJ 442, 228 (1995)
[3] Roe, Balsara, SIAM J. Appl. Math. 56, 57 (1996)
[4] Balsara, ApJS 116, 119 (1998)
[5] Powell, et al., JCP 154, 284 (1999)
ヒント
U
1
1 U
1 VP
1
1 VP
 R ARR
 RP
 RP AP RP RP
0
R
t
x
t
x
U
dW  R 1dU  R 1
dVP  RP1dVP
VP
1
HLL型近似リーマン解法
 HLL近似リーマン解法 [Harten+, 1983]
 衝撃波近似
 2-wave近似
t
SL
SR
FL
FR
UL
S R  max uL  cL , uR  cR ,0
S L  minuL  cL , uR  cR ,0
UR
i 1/ 2
S R , L：最大/最小情報伝播速度
x
 U F 
  t  x dxdt   Udx  Fdt   0
HLL型近似リーマン解法
 HLL近似リーマン解法 [Harten+, 1983]
 衝撃波近似
 2-wave近似
t
SL
SR
U
FL
FR
UL
S R  max uL  cL , uR  cR ,0
S L  minuL  cL , uR  cR ,0
UR
i 1/ 2
S R , L：最大/最小情報伝播速度
x
 U F 
  t  x dxdt   Udx  Fdt   0
 S R  S L U *  S RU R  S LU L  FR  FL  0
HLL型近似リーマン解法
 HLL近似リーマン解法 [Harten+, 1983]
 衝撃波近似
 2-wave近似
t
SL
SR
U
F
FL
UL
FR
S R  max uL  cL , uR  cR ,0
S L  minuL  cL , uR  cR ,0
UR
i 1/ 2
S R , L：最大/最小情報伝播速度
x
 U F 
  t  x dxdt   Udx  Fdt   0
S R FL  S L FR  S R S L U R  U L 
*
*
 F  FR ,L  SR ,L U  U R ,L  
SR  SL
HLL型近似リーマン解法
 HLL近似リーマン解法 [Harten+, 1983]
 衝撃波近似
 2-wave近似
S RU R  S LU L  FR  FL
U 
SR  SL
S R FL  S L FR  S R S L U R  U L 
*
F 
 F U
SR  SL

 



固有ベクトルの計算不要
正値性保存を保証
 HD [Einfeldt+, 1991] / MHD [Miyoshi+, 2005]
接触不連続の分解不可能
HLL型近似リーマン解法
 HLLD近似リーマン解法 [Miyoshi+, 2005]
 衝撃波近似
 5-wave近似
速進磁気音波×2、アルフェン波 ×2、遅進磁気音
波×2、エントロピー波×1のどの波を残すべき？
SR
t
SL
S R , L ：速進磁気音波
FL
FR
UL
UR
i 1/ 2
x
HLL型近似リーマン解法
 HLLD近似リーマン解法 [Miyoshi+, 2005]
 衝撃波近似
 5-wave近似


リーマンファンで法線方向速度一定
リーマンファンで全圧力一定
SR
t
SL
S R , L ：速進磁気音波
S M , pT
FL
UL
FR
UR
i 1/ 2
x
S M ：エントロピー波
HLL型近似リーマン解法
 HLLD近似リーマン解法 [Miyoshi+, 2005]
 衝撃波近似
 5-wave近似


リーマンファンで法線方向速度一定
リーマンファンで全圧力一定
SL
S L* t
SM
S R*
U L U L U R U R
S M , pT
FL
UL
SR
S R , L ：速進磁気音波
FR
UR
i 1/ 2
x
S M ：エントロピー波
S R* ,L ：アルフェン波
HLL型近似リーマン解法
 HLLD近似リーマン解法 [Miyoshi+, 2005]
 エントロピー波の評価 [Batten+, 1997]

u SR  uR RuR  SL  uL LuL  pTR  pTL
SM   
SR  uR R  SL  uL L

 全圧力の評価
HLL型近似リーマン解法
 HLLD近似リーマン解法 [Miyoshi+, 2005]
 速進磁気音波に対するジャンプ条件


   
   
 SM
u


 

   
2
2
 2

 SM  pT  Bx
u  pT  Bx


 u  
  SM  

 


 v  
  v  
u
B
B
S
B
B


v
v


x y
  
x y
  M
   



   

 
 
S   w     w SM  Bx Bz

  S   w   
 wu  Bx Bz


 

   


B
B
B
u
B
B
S
B
v
v


y 
x 
y M
x 


 y  
 y  




 Bz  
 Bz  
Bzu  Bx w
Bz SM  Bx w

 

     

 
 e   e  pT u  Bx v  B 
 e   e  pT SM  Bx v  B 


v  S M , v , w , B  B x , B y  , B z  ,   R , L


HLL型近似リーマン解法
 HLLD近似リーマン解法 [Miyoshi+, 2005]
 エントロピー波の評価 [Batten+, 1997]

u SR  uR RuR  SL  uL LuL  pTR  pTL
SM   
SR  uR R  SL  uL L

 全圧力の評価
pT  pTL  L SL  uL SM  uL 
 pTR  R SR  uR SM  uR 

SR  uR R pTL  SL  uL L pTR  LR SR  uR uR  uL 

SR  uR R  SL  uL L
HLL型近似リーマン解法
 HLLD近似リーマン解法 [Miyoshi+, 2005]
 速進磁気音波に対するジャンプ条件


   
   
 SM
u


 

   
2
2
 2

 SM  pT  Bx
u  pT  Bx


 u  
  SM  

 


 v  
  v  
u
B
B
S
B
B


v
v


x y
  
x y
  M
   



   

 
 
S   w     w SM  Bx Bz

  S   w   
 wu  Bx Bz


 

   


B
B
B
u
B
B
S
B
v
v


y 
x 
y M
x 


 y  
 y  




 Bz  
 Bz  
Bzu  Bx w
Bz SM  Bx w

 

     

 
 e   e  pT u  Bx v  B 
 e   e  pT SM  Bx v  B 


v  S M , v , w , B  B x , B y  , B z  ,   R , L


HLL型近似リーマン解法
 HLLD近似リーマン解法 [Miyoshi+, 2005]

 HLLD解： U
  
S  u
S  SM
HLL型近似リーマン解法
 HLLD近似リーマン解法 [Miyoshi+, 2005]
 速進磁気音波に対するジャンプ条件


   
   
 SM
u


 

   
2
2
 2

 SM  pT  Bx
u  pT  Bx


 u  
  SM  

 


 v  
  v  
u
B
B
S
B
B


v
v


x y
  
x y
  M
   



   

 
 
S   w     w SM  Bx Bz

  S   w   
 wu  Bx Bz


 

   


B
B
B
u
B
B
S
B
v
v


y 
x 
y M
x 


 y  
 y  




 Bz  
 Bz  
Bzu  Bx w
Bz SM  Bx w

 

     

 
 e   e  pT u  Bx v  B 
 e   e  pT SM  Bx v  B 


v  S M , v , w , B  B x , B y  , B z  ,   R , L


HLL型近似リーマン解法
 HLLD近似リーマン解法 [Miyoshi+, 2005]

 HLLD解： U
  
S  u
S  SM
SM  u
 
vt  vt  Bx Bt  S  u S  S   B2
M
x
 




vt  0, v, w, Bt  0, By , Bz 
2
2
 S  u   Bx
B  B
t
 t
 S  u S  SM   Bx2
HLL型近似リーマン解法
 HLLD近似リーマン解法 [Miyoshi+, 2005]
 速進磁気音波に対するジャンプ条件


   
   
 SM
u


 

   
2
2
 2

 SM  pT  Bx
u  pT  Bx


 u  
  SM  

 


 v  
  v  
u
B
B
S
B
B


v
v


x y
  
x y
  M
   



   

 
 
S   w     w SM  Bx Bz

  S   w   
 wu  Bx Bz


 

   


B
B
B
u
B
B
S
B
v
v


y 
x 
y M
x 


 y  
 y  




 Bz  
 Bz  
Bzu  Bx w
Bz SM  Bx w

 

     

 
 e   e  pT u  Bx v  B 
 e   e  pT SM  Bx v  B 


v  S M , v , w , B  B x , B y  , B z  ,   R , L


HLL型近似リーマン解法
 HLLD近似リーマン解法 [Miyoshi+, 2005]

 HLLD解： U
  
S  u
S  SM
SM  u
 
vt  vt  Bx Bt  S  u S  S   B2
 


M
x


vt  0, v, w, Bt  0, By , Bz 
2
2
 S  u   Bx
B  B
t
 t
 S  u S  SM   Bx2

e

S  u e  pT  pT  Bx v  B  v  B 

S  SM
HLL型近似リーマン解法
 HLLD近似リーマン解法 [Miyoshi+, 2005]
 アルフェン波に対するジャンプ条件


   
   
SM
 SM




  
  
2
2
 2

 2

 SM  pT  Bx
 SM  pT  Bx


  SM  
  SM  


 
 


  v  
 v  
S
B
B
S
B
B
v
v




  M
x y
x y
  M
   
   











 
S   w   

  S   w     w SM  Bx Bz
wSM  Bx Bz


   
   




B
B
B
S
B
B
S
B
v
v


y M
x 
y M
x 


 y  
 y  






 Bz  
 Bz  
Bz SM  Bx w
Bz SM  Bx w
     
     

 

 
 e   e  pT SM  Bx v  B 
 e   e  pT SM  Bx v  B 
HLL型近似リーマン解法
 HLLD近似リーマン解法 [Miyoshi+, 2005]

 HLLD解： U
  

R
S  SM 
Bx
R
t
SL

L
, S  SM 
Bx
L
SM
 L
L
SR
 R
R
i 1/ 2
x
HLL型近似リーマン解法
 HLLD近似リーマン解法 [Miyoshi+, 2005]
 アルフェン波に対するジャンプ条件


   
   
SM
 SM




  
  
2
2
 2

 2

 SM  pT  Bx
 SM  pT  Bx


  SM  
  SM  


 
 


  v  
 v  
S
B
B
S
B
B
v
v




  M
x y
x y
  M
   
   











 
S   w   

  S   w     w SM  Bx Bz
wSM  Bx Bz


   
   




B
B
B
S
B
B
S
B
v
v


y M
x 
y M
x 


 y  
 y  






 Bz  
 Bz  
Bz SM  Bx w
Bz SM  Bx w
     
     

 

 
 e   e  pT SM  Bx v  B 
 e   e  pT SM  Bx v  B 
HLL型近似リーマン解法
 HLLD近似リーマン解法 [Miyoshi+, 2005]
 アルフェン波に対するジャンプ条件
det M v , B

t

t
  0
 エントロピー波に対するジャンプ条件

  LvtL    LvtLSM  Bx BtL 
  R vtR    R vtR
SM  Bx BtR 
  SM      

SM      
 
 



 BtL   BtL SM  Bx vtL 
 BtR   BtR SM  Bx vtR 

tL

tR

t

tL

tR

t
v  v  v ,B  B  B
for Bx  0
HLL型近似リーマン解法
 HLLD近似リーマン解法 [Miyoshi+, 2005]
 
 
 
 







v
v
v








R tR
R t
L t
L vtL 
SR  SR     SR  SM     SM  SL     SL  SRL   
 BtR 
 Bt 
 Bt 
 BtL 
  R vtR 
  LvtL    R vtRuR  Bx BtR    LvtLuL  Bx BtL 
  SL 


0
 SR 








B
B
B
v
B
v

u
B
u
B
x tR  
tL L
x tL 
 tR 
 tL   tR R
SL
S L* t
S R*
SR
vtL , BtL vt , Bt vtR , BtR
vtL , BtL
vtR , BtR
i 1/ 2
x
HLL型近似リーマン解法
 HLLD近似リーマン解法 [Miyoshi+, 2005]

 HLLD解： U
  
 
 L vtL  R vtR  BtR  BtL sgnBx 
vt 
L  R


L BtR   R BtL   L R vtR  vtL sgnBx 
 
Bt 





L
R

HLL型近似リーマン解法
 HLLD近似リーマン解法 [Miyoshi+, 2005]
 アルフェン波に対するジャンプ条件


   
   
SM
 SM




  
  
2
2
 2

 2

 SM  pT  Bx
 SM  pT  Bx


  SM  
  SM  


 
 


  v  
 v  
S
B
B
S
B
B
v
v




  M
x y
x y
  M
   
   











 
S   w   

  S   w     w SM  Bx Bz
wSM  Bx Bz


   
   




B
B
B
S
B
B
S
B
v
v


y M
x 
y M
x 


 y  
 y  






 Bz  
 Bz  
Bz SM  Bx w
Bz SM  Bx w
     
     

 

 
 e   e  pT SM  Bx v  B 
 e   e  pT SM  Bx v  B 
HLL型近似リーマン解法
 HLLD近似リーマン解法 [Miyoshi+, 2005]

 HLLD解： U
  
 
 L vtL  R vtR  BtR  BtL sgnBx 
vt 
L  R


L BtR   R BtL   L R vtR  vtL sgnBx 
 
Bt 





L
R

e  e   v  B  v  B sgnBx 
 : R,  : L
HLL型近似リーマン解法
 HLLD近似リーマン解法 [Miyoshi+, 2005]
 衝撃波近似
 5-wave近似
*
t SM
S
SL
L
S R*
U L U L U R U R
S M , pT
FL
UL
SR
S R , L ：速進磁気音波
FR
UR
i 1/ 2
S M ：エントロピー波
S R* ,L ：アルフェン波
x
S R ,L U R* ,L  U R ,L   FR*,L  FR ,L , S R* ,L U R**,L  U R* ,L   FR*,*L  FR*,L ,
1 S R t
n 1


S M U  U   F  F ,
U
x
,
t
dx  S RU R  S LU L  FR  FL  0

S
t

t L
**
R
**
L
**
R
**
L
HLL型近似リーマン解法
 HLLD近似リーマン解法 [Miyoshi+, 2005]
 数値流束
F1/ 2  FL if SL  0
F1/ 2  FL  SLUL  SLUL  FL if SL  0  SL
t 0
SL
UL
FL
F1/ 2
UL
x
HLL型近似リーマン解法
 HLLD近似リーマン解法 [Miyoshi+, 2005]
 数値流束
F1/ 2  FL if SL  0
F1/ 2  FL  SLUL  SLUL  FL if SL  0  SL
F1/ 2  FL  SLUL  SL  SL UL  SLUL
 FL  SLUL  SLUL  FL if SL  0  SM
SL
SL

L
U
FL
t 0
UL
F1/ 2
UL
x
HLL型近似リーマン解法
 HLLD近似リーマン解法 [Miyoshi+, 2005]
 数値流束
FL
F
 L
FL
F1/ 2   
FR
FR

FR
if SL  0
if SL  0  SL
if SL  0  SM
if SM  0  SR
if SR  0  SR
if SR  0
F/  F/, SM , vt/, Bx , Bt/, e/, pT 
HLL型近似リーマン解法
 HLLD近似リーマン解法 [Miyoshi+, 2005]
 孤立した接線不連続（TD）の分解
t
SM
U L  U L
S M  u , Bx  0
U R  U R
x
i 1/ 2
S M U   F 
HLL型近似リーマン解法
 HLLD近似リーマン解法 [Miyoshi+, 2005]
 孤立した接線不連続（TD）の分解
 孤立した接触不連続（CD）の分解
t
U L  U L  U L
SM
U R  U R  U R
x
i 1/ 2
S M  u , Bx  0
S M U   F 
HLL型近似リーマン解法
 HLLD近似リーマン解法 [Miyoshi+, 2005]
 孤立した接線不連続（TD）の分解
 孤立した接触不連続（CD）の分解
 孤立した回転不連続（RD）の分解
S R
t

L

L
UL  U  U  U

R
S u
*
R
U R  U R
x
i 1/ 2
Bx

S R* U   F 
, Bx  0
HLL型近似リーマン解法
 HLLD近似リーマン解法 [Miyoshi+, 2005]
 孤立した接線不連続（TD）の分解
 孤立した接触不連続（CD）の分解
 孤立した回転不連続（RD）の分解
 孤立した速進衝撃波（FS）の分解
t
U L  U L  U L  U R  U R
SR
UR
x
i 1/ 2
S R U   F 
HLL型近似リーマン解法
 MHDの正値性
 物理的な解の集合

G  U |   0, e   v
2
2 B
2

20
 物理的な解の重み付き平均値
U1, 2  G  U  1   U1  U 2  G 0    1
  1   1   2  0
p  1    p1  p2

  1     1 v 1 2   B
2
2
 20
HLL型近似リーマン解法
 HLLD近似リーマン解法 [Miyoshi+, 2005]
 HLLD解の正値性
   0
 
   0

2
2






 p    1 e   v 2  B 2   0



 

  2
 2

  0


p

1
e

v
2
B
2





 







HLL型近似リーマン解法
 HLLD近似リーマン解法 [Miyoshi+, 2005]
  S R  uR  0,  S R  S M  0,   S M  uR
 密度の正値性

    R  0


R

R
HLL型近似リーマン解法
 HLLD近似リーマン解法 [Miyoshi+, 2005]
  S R  uR  0,  S R  S M  0,   S M  uR
 圧力の正値性

  2
 2

    eR   R vR 2  BR 2 


2


B
 R 
pR
tR

2
  pR 
1

2 

 1
2   R  Bx 


2


BtR
 R 
pR

2



p




1


R
2
2
 1
2   R c fR  Bx 


一生懸命テキストの方に書きました。
HLL型近似リーマン解法
 HLLD近似リーマン解法 [Miyoshi+, 2005]
  S R  uR  0,  S R  S M  0,   S M  uR
 圧力の正値性
D    0     0



HLL型近似リーマン解法
 HLLD近似リーマン解法 [Miyoshi+, 2005]
  S R  uR  0,  S R  S M  0,   S M  uR
 圧力の正値性
2

 2
BtR
2  R pR 
2
  0
1
D    pR 
  1   R c 2fR  Bx2 

  1 pR 
 
1

2
2 R

1

BtR
    1 c 2fR
2
 R c 2fR  Bx2 
 1
S R  uR 
c fR
2
2
HLL型近似リーマン解法
 HLLD近似リーマン解法 [Miyoshi+, 2005]
  S R  uR  0,  S R  S M  0,   S M  uR
 圧力の正値性

  2
 2

p    1 eR   R vR 2  BR 2 


 pR  0

R
 正値性保存の条件
 1
 1
S R  uR 
c fR , S L  uL 
c fL
2
2
HLL型近似リーマン解法
 HLL型近似リーマン解法
 HLLD近似解の重み付き平均値
 正値性保存
(MHD HLL-type)
HLL型近似リーマン解法
 HLLD近似リーマン解法 [Miyoshi+, 2005]
 固有ベクトルの計算不要
 正値性保存を保証
 接線不連続を厳密に分解
 接触不連続を厳密に分解
 回転不連続を厳密に分解
 速進衝撃波を厳密に分解
 支配方程式に依存
MHD方程式の近似リーマン解法
 精度・計算速度の検証
 ロバスト性の検証
[Mignone+, 2007]
MHD方程式の近似リーマン解法
 その他の数値実験は元論文などを参照してください
[1] Miyoshi, Kusano, JCP 208, 315 (2005)
[2] Mignone, et al., ApJS 170, 228 (2007)
[3] Kritsuk, et al., ApJ 737:13 (2011)
わかりました？
 どこの風を感じるのか違いはわかりましたか？
 流束ベクトル分離法
F   A U
セル内の風を感じる
 近似リーマン解法
F   A  U
セル境界の風を感じる
２つのlocal Lax-Friedrichs法

i 1 / 2
F
Fi 1/ 2

Fi 1  Fi 1

  max i 1U i 1  |  |max i U i
2
2
Fi 1  Fi  max i 1/ 2
U i 1  U i 


2
2

何かご利益があるかどうかは別として、
流束ベクトル分離法的HLL法も作れますね！
流束ベクトル分離法的HLL法
 作れますね！と書いたけどちょっと不安
1m
1m
 1

F 
F
U, F 
F
U
m  1
m  1
m  1
m  1

m
でできるよね？
 1

m 
Λ 
m  1 



 1

 1 
Λ 
m  1 



2
2

1




 1m  1


0





m  1





m 
1 



1



 1m  1


   
0


m
1





m 
1  できた！

かな？
2015年8月3日（月）-7日（金）千葉大学統合情報センター
宇宙磁気流体・プラズマシミュレーションサマースクール
フローチャート
フローチャート
 基本手順
開始
初期設定
CFL条件の評価
数値流束の評価
時間積分
t = tout
Yes
データ出力
Yes
t < tend
No
終了
No
フローチャート
 基本手順
パラメータ
 ,,
開始
初期設定
グリッド
xi ,  x , 
CFL条件の評価
初期条件
数値流束の評価
U  xi , t  0 
時間積分
t = tout
Yes
データ出力
Yes
t < tend
No
終了
No
フローチャート
 基本手順
開始
初期設定
i  u i  c f i
CFL条件の評価
t  cCFL
数値流束の評価
時間積分
t = tout
Yes
データ出力
Yes
t < tend
No
終了
No
x
max i
フローチャート
 基本手順
開始
初期設定
for i  1 ,  , N  1
U L  U i 1
CFL条件の評価
数値流束の評価
U R  Ui
F1 / 2 U L , U R 
時間積分
t = tout
Yes
データ出力
Yes
t < tend
No
終了
No
フローチャート
 基本手順
開始
初期設定
CFL条件の評価
数値流束の評価
時間積分
t = tout
Yes
データ出力
Yes
t < tend
No
終了
No
for i  1 ,  , N
t
Fi
Ui  Ui 
x
フローチャート
 基本手順
開始
初期設定
CFL条件の評価
数値流束の評価
時間積分
t = tout
Yes
データ出力
Yes
t < tend
No
終了
No
if t  t out then
Data Output
t out  t out  t out
フローチャート
 高次精度化
開始
初期設定
CFL条件の評価
SSP-RK法
高次精度補間法
数値流束の評価
時間積分
t = tout
Yes
データ出力
Yes
t < tend
No
終了
No
フローチャート
 MUSCL法
CFL条件の評価
n-step SSP-RK法
MUSCL法
数値流束の評価
for i  0 ,  , N  1
x
 U
Lim 
UR  Ui 
2
 x
x
 U
Lim 
UL  Ui 
2
 x
UL
時間積分
UR
x 2
x 2
i






フローチャート
 MUSCL法
CFL条件の評価
n-step SSP-RK法
MUSCL法
数値流束の評価
時間積分
for i  1 ,  , N  1
F1 / 2 U L , U R 
フローチャート
 MUSCL法
（2nd-order SSP（TVD）-RK）
CFL条件の評価
n-step SSP-RK法
MUSCL法
数値流束の評価
時間積分
U 1  U n   tL U n 
1 n 1 1 1
n 1
U  U  U   tL U 1 
2
2
2
（3rd-order SSP（TVD）-RK）
U 1  U n   tL U n 
3
1
1
U 2   U n  U 1   tL U 1 
4
4
4
1 n 2 2  2
n 1
U  U  U   tL U 2  
3
3
3
フローチャート
 FV-WENO法など変数補間で高次精度化を実現する
数値解法は同様のフローチャートです。
 数値流束を高次精度補間する数値解法もあります。
 補間の実装は結構面倒です。簑島先生よろしく！
 多次元化では、ざっくりと言って、Split法もUnsplit法も
単に1次元の拡張です。
 ただし、数値的な磁場発散の処理が必要です。これ
も大変面倒。簑島先生よろしく！
2015年8月3日（月）-7日（金）千葉大学統合情報センター
宇宙磁気流体・プラズマシミュレーションサマースクール
おしまい
お疲れ様でした。

Download Report