Upwind MHD Code メモ

風上系のMHDコードを書いたのでそのメモとギャラリー。

最近業界ではHLLD法が流行りのようだけれど、 状態方程式が変わるたびにスキームを構築しなおすのは面倒なので、 拡張が簡単なものをとりあえず実装。

Reconstruction Step

Riemann Solverでセル境界のRiemann問題を解くためには、 セル境界での物理量の値を補間するステップが必要で、これをReconstuction Stepといいます。 有限差分法ではセル中心における物理量から、 有限体積法ではセル内の物理量の平均値から、それぞれ境界値を再構築します。 現時点で実装完了したスキームは以下の表のとおり。

Piecewise Constant Method セル内一定値補間。つまり何もしない。1次精度。[Godunov(1959)]
Piecewise Linear Method セル内1次多項式補間。流束制限関数でTVD化。2次精度。[van Leer(1979)]
Piecewise Parabolic Method セル内2次多項式補間。流束制限関数でTVD化。3次精度。[Kermani et al.(2003)]
3rd-order WENO Weighted Essentially Non-Oscillatory scheme。異なるステンシルで1次多項式補間し、不連続を含まないように滑らかさに応じた重みで足し合わせることで、スムーズな領域で3次精度を達成。[Liu et al.(1994)]
5th-order WENO WENOの5次精度版。2次多項式補間を足し合わせる。[Jiang and Shu(1996)]

Riemann Solver

各セル境界に右側と左側からそれぞれ補間された物理量の値は、一般には異なる。 この２つの状態を左右の初期状態として（近似的に）Riemann問題を解くことで、 次の時間ステップの物理量を求めよう、 というのが（近似）Riemann Solverの基本的なコンセプト。 実装済みのものは以下のとおり。

HLL Harten-Lax-van Leer scheme。特性速度最大の波だけを解く近似Riemann Solver。接触不連続面やシアー成分で数値拡散が大きい。HLLC、HLLDなどのHLL系スキームの原型。[Harten et al.(1983)]
Kurganov and Tadmor central scheme HLLの背景移流を考えない版。Rusanovスキーム、Local Lax-Friedrichsスキームとも。[Kurganov and Tedmor(2000)]

時間積分

今回のコードは、一貫して3rd-order TVD Runge-Kutta schemeを使用しました。

div(B)の取扱い

MHD計算では、数値的なエラーによって本来ないはずの磁場の発散が生じるので、 これをケアすることがとても大事。 方法はいくつかあるが、今回はConstrained Transportを採用。 実装の詳細については、Balsara and Spicer(1999)を参照。

ギャラリー

テスト計算をした中から、いくつか見栄えのいいものを。

Orszag-Tang vortex problem

2次元MHDの有名なテスト問題。3次元版もあるらしい。 密度、圧力一様のプラズマ中を超音速の渦が流れていて、 まわりの磁場と相互作用して何枚もの衝撃波ができる。 解像度は512×512。

左：Piecewise Parabolic Method + van Albada's flux limiter + HLL
右：5th-order WENO Reconstruction + HLL

Gas Pressure @ t = 3.14

両者とも精度よく衝撃波を捉えており、見た目には違いは分からない。

Gas Pressure @ t = 6.28

渦の中心に発達した一枚の大きな電流シートに注目すると、 左の3次精度PPMの計算では電流シートは安定なままだが、 右の5次精度WENOの計算では、電流シートがグリッドスケールまで薄くなると 数値拡散によりPetschekライクな磁気リコネクションを起こす （陽に電気抵抗を入れていないので、すべて数値拡散によるもの）。 低次精度の計算では、 数値拡散が効きすぎるために電流シートがグリッドスケールまで局在化できず、 Sweet-Parker的なリコネクションを起こすのかも。

Balsara's Rotor problem

こちらも有名なテスト問題。 背景の静止プラズマの中心に回転する円柱（円盤？）を置き、 そこからAlfven波が伝搬する様子を見るもの。 解像度は512×512。

左：Piecewise Parabolic Method + van Albada's flux limiter + HLL
右：5th-order WENO Reconstruction + HLL

Gas Pressure @ t = 0.15