Introduction toMagnetohydrodynamic
Simulations
Seiji ZENITANIDivision of Theoretical Astronomy
SOKENDAI 2014/2015

Outline• Day 1 - Magnetohydrodynamics (MHD)• Introduction & MHD gallery• Basic equations, fundamental notions• MHD waves• Day 2
• MHD Theory (cont.) - Conservation form, shocks• Simulations - Riemann solver, HLL, TVD• Simulations - Recent topics (Higher-order schemes, div B)• Assignments - 1D simulation

What is MagnetoHydrodDynamics (MHD)?
磁気流体力学／電磁流体力学、磁流体力学（中文）• A fluid framework for macroscopic plasma motion• “Plasma” consists of...
• Charged particles (ionized gas)
• Magnetic fields

PlasmaBionAtomAtomicnucleus

Magneto fluidelectronelectron

magnetic
field line

MHDmotion

MHD gallery

MHD gallery - solar corona• Highly filamented structure• Triggering of flares (magnetic eruption)

Kusano et al. 2012 ApJ

NASA / TRACE satellite

MHD process - Magnetic reconnection• Underlying mechanism of solar flares

Jet

Jet
Zenitani & Miyoshi 2011 Phys. Plasmas

Turbulent reconnection

MHD gallery - Heliosphere• Global simulations ofHeliosphere & planetary magnetosphere

Shiota et al. 2014 Space WeatherMatsumoto et al. 2010 IEEE trans.

MHD process - Kelvin-Helmholtz instability

Weathernews.jp 2015

Matsumoto & Hoshino 2004 GRL

MHD gallery - Pulsar M’shpere• Magnetically-dominated environments• Exotic variants of MHD - Force-free / Relativistic MHD

Spitkovsky 2011 ASSP

Tchekovskoy et al. 2013 MNRASMHD gallery - Black hole system• Black hole accretion disks + jet launching

Machida et al. 2003 ApJ

Kato et al. 2004 ApJ

MHD equations

Basic equations - fluid part

• Continuity• Momentum• Energy*
Lorentz force

Lorentz force = magnetic pressure+ magnetic tension• Examples: “magnetic pressure”“magnetic tension”
B = (y,1,0)• B = (0,x,0)

How magnetic tension works ...

• It tends to be straightBasic equations - EM part

• Faraday’s law• Ampere’s law~O(V2/c2)• Of course it follows the constraints in Maxwell eqs.

MHD equations ...?

• Continuity• Momentum• Energy*• Magnetic field
We have 8 equations, but 11 unknowns...

• Constraint

Closure - ideal condition• Ideal condition• No electric field in the moving frame at v• Particle motion in the perpendicular plane
p+vE B =Ee-cE B
B2

Ideal MHD equations• Continuity• Momentum• Energy*• Magnetic field• Ideal condition• Constraint

Scales of MHD• Time:Plasmaci =eB
mi celectronmagnetic
field line

Alfvén speed

Length:Log TB

MHDion1ci⌦ce1Kinetic
model!pe1

rL,e

di

Log L

Flux frozen-in• Magnetic flux in a comoving surface is conserved• Similar to Kelvin’s circulation theorem in hydrodynamics
Bvdl• Magnetic field is “frozen-in” to the flow field VStern 1966 SSR

Line connectivity• Convection of a line element• From the induction equation andthe continuity equation ...

Landau & Lifshitz 1960

• Line element (δl) is parallel to magnetic field (B)δlPriest & Forbes 2000

Field line motion in ideal MHD• All in all, it is convenient to consider that plasmas and“magnetic field lines” move together at the same speed
Bp+Ee-Resistivity• We sometimes consider the resistivity in the MHD frame• Plasma version of Ohm’s law, or simply “Ohm’s law”
• This modifies the induction equation in the following way

Convection term

Diffusion term

• Energy equation needs to be modified accordingly, because
it brings Ohmic heating.

General discussion...

• Ohm’s law• Ideal condition• Frozen-in condition• Line connection• Topology condition ...
• In nonideal case of R≠0, it is not always possible to define
the field line velocity.

• Don’t worry, eventually, magnetic field lines are virtual lines.

MHD waves

Waves in hydrodynamics• Linearization
0: average value• Hydrodynamics• Sound (compressional) wave

1: small oscillation

MHD waves - Alfvén wave

• Incompressible caseCompressional part

• Alfvén wave - transverse wave by magnetic tension• See the movie

Alfvén wave paper (1942)

Hannes Alfvén (1908-1995)

Alfvén 1942 Nature• 1970 Novel Prize in Physics

Alfvén wave - observation• Solar photosphere, observed by Hinode satellite

http://hinode.nao.ac.jp/gallery/Okamoto et al. 2007 Science

MHD waves - magnetosonic waves• Small oscillation• MHD equations• After SOME algebra

MHD waves - magnetosonic waves• Again, linearized MHD equationszB

Compressional* waveby magnetic pressure

Alfvén wave

θ(Parallel propagation)ykSound wave
(Spherical propagation)• Dispersion relation:• Fast/slow magnetosonic waves• Compressional waves due to fluid and magnetic* pressuresx

MHD waves - magnetosonic waves• Alfvén wave
• Fast and slow

magnetosonic waves• Friedrichs diagram (phase-speed)•c > c•c < cs

A

sAcAcosθ

BB

MHD waves - simulation results• MHD waves are no longer spherically symmetric• Hydrodynamic wave

Credit: Dr. Magara• MHD (c > c )sA• MHD (c < c )sA

http://solar.physics.montana.edu/magara/Research/Topics/MHD_waves/mhdw.html

Summary• MHD = Magnetohydrodynamics• Plasmas and magnetic field lines move together• MHD physics highly depends on the direction of B
•Magnetic pressure vs magnetic tension•Alfvén wave
•Magnetosonic waves

Further reading•宇宙流体力学、坂下＆池内 (1996)• Plasma Physics for Astrophysics, R. M. Kulsrud (2004)•Magnetohydrodynamics of the Sun, E. R. Priest (2014)•シリーズ現代の天文学 12. 天体物理学の基礎 II, ２章•シリーズ現代の天文学 14. シミュレーション天文学• W. A. Newcomb, Ann. Phys. (1958)• gnuplot code for Friedrichs diagram
set polar
unset keya=1; s=1.2
plot a*abs(cos(t)) lw 3, \sqrt(0.5*(a**2+s**2+sqrt((a**2+s**2)**2-(2*a*s*cos(t))**2))) lw 3, \sqrt(0.5*(a**2+s**2-sqrt((a**2+s**2)**2-(2*a*s*cos(t))**2))) lw 3

Outline• Day 1 - Magnetohydrodynamics (MHD)• Introduction & MHD gallery• Basic equations, fundamental notions• MHD waves• Day 2
• MHD Theory (cont.) - Conservation form, shocks• Simulations techniques - Riemann solver, HLL, MUSCL• Simulations techniques - div B• Assignments - 1D simulation

Review - ideal MHD equations• Continuity• Momentum• Energy*• Magnetic field• Ideal condition• Constraint

Review - MHD waves• Alfvén wave
• Fast and slow

magnetosonic waves• Friedrichs diagram (phase-speed)•c > c•c < csAs

AcAcosθBB

Ideal MHD equations - conservation form

• Continuity• Momentum• Energy• Mag. field• Constraints

Ideal MHD equations - conservation form

Primitive variables（基本量）

Conserved quantities（保存量）

Numerical flux
（数値流束）

Source term

Characteristic variables（特性変数）

Finite volume method• Conservation form is friendly with finite volume method• One can deal with shocks (“Shock-capturing” code)

Shocks in the Universe

Mach 1887

NASA/CXC/Rutgers/J.Hughes et al. 2005
Jokipii 2008

Shocks• Sudden jump in physical quantities• Nonlinear wave steepens into a shock• Supersonic flow slows down to a subsonic speed• Fast explosion propagates into a medium...

• Shock is a big enemy to numerical simulations!

MHD shocks & discontinuities
Fast shock (FS)Intermediate shock*

Rotational discontinuity

Slow shock (SS)Contact discontinuity

Tangential discontinuity

Upwind schemec• We use numerical flux on the “upwind” side
Generic form

• Time-evolution by the 1st-order Euler method
Diffusion term

Monotonicity（単調性）

Riemann solver (1/2)• Riemann problem

U• Time evolution from two flat states• Nonlinear problem in HD• Riemann solver• Solve Riemann problem at each cellinterface

UUx

Riemann solver (2/2)• Switches to an Upwind schemetx
tx

UUMHD Riemann problem• MHD Riemann problem has no analytic solution• Ex. Brio=Wu problemtx• Left state• (ρ, p, v, B , B , B ) = (1, 1, 0, 0.75, 1, 0)
xyz• Right state:• (ρ, p, v, B , B , B ) = (0.125, 0.1, 0, 0.75, -1, 0)xBrio & Wu 1988yz※ B is renormalized by [4π]-0.5

Approximate Riemann solver - HLL solver (1/2)• Harten–Lan–van Leer (1983)

tUx
Utx

Approximate Riemann solver - HLL solver (2/2)
Fastest signal = fast wave

tx
Jump conditions
for the λ-moving discontinuity

Note that

F hll = F (U hll )Better Riemann solvers• HLLD solver (Miyoshi & Kusano 2005)• Relativistic HLLD solver (Mignone et al. 2009)• Exact solver (Takahashi & Yamada 2014)

Higher accuracy• Time-marching - Higher-order Runge-Kutta method• Spatial profile - MUSCL interpolation (linear version)• ex. minmod limiterRR

L

i-1/2

i+1/2

ii+1

MUSCL interpolation• We should interpolate characteristic variables for higher (≧3 ) order.
• We often interpolate primitive/conservative variables for 2 orderrdndLRlimiterfunction (≧0)

B

B=2r

B=r210r0≦x≦1Multi-dimensional problem• Fractional step• Compute x for Δt/2 (L ) --> Compute y for Δt/2 (L ) --> ...
• U(t+Δt) = L (Δt/2) L (Δt/2) L (Δt) L (Δt/2) L (Δt/2) U(t)xx• All at once

yyzyx

Multi-dimensional problem• MHD scientists worry about the “solenoidal condition”• Non-zero error:• The magnetic topology is distorted• Theorists are just unhappy• Unphysical force grows in the system

Unphysical force

Divergence cleaning methods (1/2)• 1. Projection scheme (Brackbill & Barnes 1983, Tóth 2000)• Solve Poisson equation every timestep• 2. Constraint-Transport (Evans & Hawley 1988)• Set up the staggered grid for B

Divergence cleaning methods (2/2)• 3. Hyperbolic divergence cleaning method (Dedner et al. 2002)• A new virtual potential• Telegraph equation of div B• 4. Artificial diffusivity for finite difference code (Kawai 2013)• And many others

[Option] Assignment• Use an MHD code or write your own code• Athena - https://trac.princeton.edu/Athena/• CANS+ - http://www.astro.phys.s.chiba-u.ac.jp/cans/doc/• PLUTO - http://plutocode.ph.unito.it/• OpenMHD - http://th.nao.ac.jp/MEMBER/zenitani/openmhd-j.html• In these codes, the factor 4π is often set to 1.

• Try EITHER of the following problems, plot something, andexplain the figure in a paragraph by 12 Feb 2015.

• 1. Brio-Wu problem• 2. MHD fast shock problem

MHD fast-shock problem• Consider a Riemann problem of interpenetrating flows• We will see shocks when the flow speed exceeds the fast speed,t

By=B0+V-Vx

t

Shock propagationx

Further reading•流体力学の数値計算法、藤井孝藏（1994）•岩波講座地球惑星科学６ 地球連続体力学、寺澤敏夫（1996）• A multi-state HLL approximate Riemann solver for ideal

magnetohydrodynamics, Miyoshi & Kusano, J. Comput. Phys. (2005)• The div B = 0 Constraint in Shock-Capturing

Magnetohydrodynamics Codes, Toth, J. Comput. Phys. (2000)•数値シミュレーションは div B=0 を守れるか?、三好隆博、物理学会
誌（2014年７月号）