The Richtmyer–Meshkov instability in magnetohydrodynamics
V. Wheatley,
1
R. Samtaney,
2
and D. I. Pullin
3
1
School of Mechanical Engineering, University of Adelaide, Adelaide, South Australia 5005, Australia
2
Princeton Plasma Physics Laboratory, Princeton University, Princeton, New Jersey 08543, USA
3
Graduate Aeronautical Laboratories, California Institute of Technology, Pasadena, California 91125, USA
Received 8 February 2009; accepted 23 June 2009; published online 25 August 2009
In ideal magnetohydrodynamics
MHD
, the Richtmyer–Meshkov instability can be suppressed by
the presence of a magnetic field. The interface still undergoes some growth, but this is bounded for
a finite magnetic field. A model for this flow has been developed by considering the stability of an
impulsively accelerated, sinusoidally perturbed density interface in the presence of a magnetic field
that is parallel to the acceleration. This was accomplished by analytically solving the linearized
initial value problem in the framework of ideal incompressible MHD. To assess the performance of
the model, its predictions are compared to results obtained from numerical simulation of impulse
driven linearized, shock driven linearized, and nonlinear compressible MHD for a variety of cases.
It is shown that the analytical linear model collapses the data from the simulations well. The
predicted interface behavior well approximates that seen in compressible linearized simulations
when the shock strength, magnetic field strength, and perturbation amplitude are small. For such
cases, the agreement with interface behavior that occurs in nonlinear simulations is also reasonable.
The effects of increasing shock strength, magnetic field strength, and perturbation amplitude on both
the flow and the performance of the model are investigated. This results in a detailed exposition of
the features and behavior of the MHD Richtmyer–Meshkov flow. For strong shocks, large initial
perturbation amplitudes, and strong magnetic fields, the linear model may give a rough estimate of
the interface behavior, but it is not quantitatively accurate. In all cases examined the accuracy of the
model is quantified and the flow physics underlying any discrepancies is examined. ©
2009
American Institute of Physics
.
DOI:
10.1063/1.3194303
I. INTRODUCTION
The Richtmyer–Meshkov instability
RMI
is important
in a wide variety of applications
1
including inertial confine-
ment fusion
2
and astrophysical phenomena.
3
In these appli-
cations, the fluids involved may be ionized and hence be
affected by magnetic fields. Samtaney
4
demonstrated, via nu-
merical simulations, that the growth of the RMI is sup-
pressed in the presence of a magnetic field. The particular
flow studied was that of a shock interacting with an oblique
planar contact discontinuity
CD
separating conducting flu-
ids of different densities within the framework of planar
ideal magnetohydrodynamics
MHD
. It was shown that the
suppression of the instability is caused by changes in the
shock refraction process at the CD with the application of a
magnetic field.
5
These changes prevent the deposition of cir-
culation on the CD.
A more widely studied flow results from a shock wave
accelerating a density interface with a single-mode sinu-
soidal perturbation in amplitude. Our goal is to understand
the effect of a magnetic field on this flow when conducting
fluids are involved. The magnetic field is again aligned with
the motion of the shock. The initial condition for this flow is
illustrated in Fig.
1
a
. It is characterized by the incident
shock sonic Mach number
M
, the density ratio across the
CD,
2
/
1
, the ratio of the CDs initial amplitude to its wave-
length,
0
/
, the ratio of specific heats,
, and the non-
dimensional strength of the applied magnetic field,
−1
=
B
2
/
2
p
0
. Here
B
is the magnitude of the applied mag-
netic field and
p
0
is the initial pressure in the unshocked
regions of the flow. To demonstrate the suppression of the
instability in this geometry, the flow was simulated both in
the presence and absence of a magnetic field. For both simu-
lations,
M
=2,
2
/
1
=3,
0
/
= 0.1, and
=5
/
3. In the simu-
lation in which a magnetic field is present,
−1
= 1. The nu-
merical method used is described in Sec. III D. Time
sequences of density and vorticity fields from the simulations
without,
i
, and with,
ii
, an initial magnetic field are shown
in Fig.
2
.In
i
, the transmitted and reflected shocks are
clearly visible in the density fields. Note the presence of the
transverse waves that are generated downstream of the lead-
ing transmitted and reflected shocks in each simulation. It
can be seen from Fig.
2
i
that the vorticity generated by the
shock refraction process remains at the interface. This causes
the interface to roll up into the mushroom shape characteris-
tic of the hydrodynamic RMI. In
ii
, the transmitted and
reflected fast shocks that are produced by the interaction are
clearly visible in the density fields. The tangential velocity
jumps across these shocks are considerably smaller than
those across the slow shocks; thus they are not as visible in
the vorticity fields. The transmitted and reflected slow shocks
have small density jumps across them and therefore do not
feature prominently in the density fields. It can be seen from
the vorticity fields, however, that the majority of the vorticity
generated during the shock refraction process is transported
away from the interface via the tangential velocity jumps
across the slow shocks. In the ideal case, this leaves the
PHYSICS OF FLUIDS
21
, 082102
2009
1070-6631/2009/21
8
/082102/13/$25.00
© 2009 American Institute of Physics
21
, 082102-1
Downloaded 16 Oct 2009 to 131.215.220.165. Redistribution subject to AIP license or copyright; see http://pof.aip.org/pof/copyright.jsp
interface with zero circulation per unit length, which drasti-
cally alters its evolution. This is evident from a comparison
of Fig.
2
c
i
and Fig.
2
c
ii
, where the interface from the
simulation with a magnetic field present exhibits none of the
roll-up seen in the hydrodynamic case.
A second case of interest occurs when the initial mag-
netic field is parallel to the interface. In this case a jump in
tangential velocity across the interface is permitted, allowing
the majority of the vorticity generated by the shock refrac-
tion process to remain on the interface. Thus the instability
may not be suppressed by transport of vorticity. An incom-
pressible model for this case was formulated by Cao
et al.
6
They determined that when the magnetic field is parallel to
the interface and there is no mean shear, the instability is
suppressed by a Lorentz force that opposes the perturbation
of the interface and no discontinuous waves are generated.
As the characteristics of the flow in this case differ substan-
tially from that shown in Fig.
2
, we do not consider it any
further here.
As a model for this flow in the case where a magnetic
field aligned with the shock motion is present, we examined
the growth of a sinusoidally perturbed interface separating
incompressible conducting fluids that is impulsively acceler-
ated at
t
=0.
7
The setup for the model problem is illustrated
in Fig.
1
b
. This problem is characterized by
1
/
,
2
/
,
0
/
,
, and the normalized magnitude of the impulse,
V
/
p
0
. Appropriate values of
V
and
are computed from
the corresponding shock driven flow;
V
is the change in
mean interface velocity produced by the shock interaction
process, while
is selected to be the post-interaction value
of
1
. We shall subsequently refer to this as the incompress-
ible linear theory
ILT
. The ILT differs from the full MHD
RMI in that it is incompressible, linear, and driven by an
impulse rather than by the impact of a shock wave. Here, the
performance of the model is assessed for a variety of cases
by comparing it to the results of compressible numerical
simulations. In each case, an impulse driven linearized
IDL
simulation, a shock driven linearized
SDL
simulation, and
a nonlinear
NL
simulation were carried out. This allows the
effects on the flow of compressibility, shock acceleration,
and nonlinearity to be assessed systematically: differences
between the ILT and an IDL simulation are mainly due to the
effects of compressibility, differences between IDL and SDL
simulations are due to the effects of shock rather than impul-
sive acceleration, and differences between SDL and NL
simulations are due to NL effects.
It is expected that, for fixed
, the effects of compress-
ibility increase with shock Mach number
M
, while NL ef-
fects increase with the initial amplitude of the interface
0
.
From the ILT, the propagation speeds of the fronts that carry
circulation away from the interface scale like the Alfvén
speed,
C
A
=
B
/
. This must be small compared to the sound
speed, which corresponds to
2
a
2
/
C
A
2
being large, if
these fronts are not to interact with the shocks present in the
compressible case. Thus it is anticipated that the ILT will be
most accurate for a flow characterized by small
M
, small
0
,
and large
. The performance of the model for such a set of
parameters is analyzed as baseline case. We then examine
how the performance of the model is affected as
M
,
0
, and
−1
are increased.
In the sequel, we first present the incompressible model
for the MHD Richtmyer–Meshkov flow, including features
not previously discussed. The numerical methods and simu-
lation setups utilized are described next. Following this, we
analyze the results of all methods for a baseline case from
which we began our investigations. In this section, we
present many finer details of the MHD Richtmyer–Meshkov
flow that were not discussed in previous works. The effects
of increasing shock strength, magnetic field strength, and
perturbation amplitude are then investigated sequentially. Af-
ter examining the performance of the model for a case where
all of these parameters are large, we present the conclusions
that have arisen from this investigation.
FIG. 2.
Color online
Negative vorticity and density fields from compress-
ible simulations with
M
=2,
2
/
1
=3,
0
/
= 0.1,
=5
/
3, and
i
B
=0 or
ii
=1 at
a
t
/
t
= 0.2,
b
t
/
t
= 0.8, and
c
t
/
t
= 3.4. The top half of each plot
shows vorticity while the bottom half shows density. At the time of these
images, the incident shock has interacted with the interface.
(a)
λ
shock
M
B
ρρ
x
z
η
0
12
(b)
λ
ρ
B
ρ
z
x
η
0
12
Vt
δ( )
FIG. 1.
a
Initial condition geometry for compressible RMI.
b
Geometry
for incompressible model problem.
082102-2 Wheatley, Samtaney, and Pullin
Phys. Fluids
21
, 082102
2009
Downloaded 16 Oct 2009 to 131.215.220.165. Redistribution subject to AIP license or copyright; see http://pof.aip.org/pof/copyright.jsp
II. INCOMPRESSIBLE MODEL
The derivation of the incompressible model for the
MHD Richtmyer–Meshkov flow was presented in Wheatley
et al.
7
To derive the model, the linearized equations of ideal,
incompressible MHD are solved in a noninertial reference
frame that has acceleration
V
t
in the
z
-direction. Here,
t
is the Dirac delta function and
V
c
. The equations are
linearized about a base flow that results from the impulsive
acceleration of an unperturbed interface. Our choice of ref-
erence frame results in the horizontal velocity
w
of the base
flow being zero for all time. The resulting linearized equa-
tions are subject to an impulsive forcing that is nonzero only
in a small region between
z
= 0 and the perturbed interface
location
z
=
h
x
,
t
=
t
e
ikx
h
. All perturbations are as-
sumed to have the form
q
x
,
z
,
t
=
q
ˆ
z
,
t
e
ikx
. Subject to
boundary and initial conditions, the linearized equations can
be solved to yield the following solutions for
w
in each fluid
see Ref.
7
for details
:
w
1
=
a
1
t
e
kz
+
H
t
+
z
/
C
A
1
c
1
t
+
z
/
C
A
1
e
ikx
,
1
w
2
=
b
2
t
e
−
kz
+
H
t
−
z
/
C
A
2
d
2
t
−
z
/
C
A
2
e
ikx
,
2
where
H
z
is the Heaviside function. In the above expres-
sions,
a
1
t
=
K
A
2
1
2
e
1
t
1
−
2
+
2
+
R
+
i
1
+
+
i
e
+
i
t
i
+
i
−
1
,
3
b
2
t
=
K
A
2
2
2
e
2
t
2
−
2
+
2
+
R
+
i
2
+
+
i
e
+
i
t
i
+
i
−
2
,
4
c
1
t
=
K
C
1
+
2
e
1
t
1
−
2
+
2
+
R
2
+
+
i
e
+
i
t
i
+
i
−
1
,
5
d
2
t
=
K
D
1
+
2
e
2
t
2
−
2
+
2
+
R
1
+
+
i
e
+
i
t
i
+
i
−
2
,
6
where
1
=
Bk
1
,
2
=
Bk
2
,
=−
Bk
1
+
2
1
+
2
,
7
=
B
2
k
2
1
+
2
−2
1
2
1
/
2
1
+
2
,
8
and
K
A
=
kV
0
A
,
K
C
=−
2
Bk
2
V
0
A
2
1
2
+
2
1
,
K
D
=
1
2
K
C
.
The Atwood number
A
2
−
1
/
2
+
1
. The above ex-
pressions are not valid if
= 0, but this requires that either
B
=0,
k
=0, or
1
=
2
, which corresponds to cases that are not
of interest here.
A. Initial solution and growth rate
Profiles of
w
ˆ
z
,
t
at various times are shown in Fig.
3
for one set of parameters. The initial
t
=0
+
velocity
distribution,
w
x
,
z
,0
+
=
0
kV
A
e
−
k
z
+
ikx
,
9
is identical to the steady velocity distribution that arises from
the hydrodynamic
B
=0
case. This implies that the initial
growth rate of the interface, which to leading order is given
by
w
ˆ
i
0,0
, is unaffected by the presence of a magnetic field.
Indeed, from Eq.
1
or Eq.
2
it can be shown that this
initial growth rate is
t
t
=0
=
0
kV
A
,
10
as in the hydrodynamic case.
8
This is consistent with the fact
that the baroclinic generation of vorticity
is unaffected by
the presence of the magnetic field.
B. Circulation distribution
On any interface with unit tangent
t
ˆ
, the circulation per
unit length
u
is given by
u
=
u
·
t
ˆ
.
For the interfaces in our problem,
u
·
t
ˆ
=
u
to leading order.
Using the fact that
u
=
iDw
/
k
, where
D
denotes a derivative
with respect to
z
,
−3
−2
−1
0
1
2
−1
0
1
2
3
4
5
6
7
8
x10
−
3
z
/
λ
normalized
w
perturbation amplitude
t/t*
=0
t/t*
=1
t/t*
=4
CD location
FIG. 3. Profiles of
w
ˆ
z
,
t
/
p
0
at
t
/
t
=0,
t
/
t
= 1, and
t
/
t
= 4 for
1
/
= 1.483 72,
2
/
= 4.431 59,
V
/
p
0
= 0.319 125,
0
/
= 0.007 992 76, and
= 16. Here
t
/
p
0
. The maxima of
w
ˆ
z
,
t
coincide with the Alfvén
fronts.
082102-3 The Richtmyer–Meshkov instability in MHD
Phys. Fluids
21
, 082102
2009
Downloaded 16 Oct 2009 to 131.215.220.165. Redistribution subject to AIP license or copyright; see http://pof.aip.org/pof/copyright.jsp
u
=
u
=
i
k
Dw
=
i
k
Dw
ˆ
e
ikx
.
Thus the gradient discontinuities in
w
ˆ
seen in Fig.
3
indicate
the presence of interfaces that carry circulation on a half-
period. At
t
=0
+
, Fig.
3
shows that circulation is present on
the density interface at
z
= 0. This was baroclinically gener-
ated during the impulsive acceleration of the interface. Away
from the interface the flow is irrotational at
t
=0
+
; thus the
total circulation in a half-period of the domain must be equal
to
1
/
2
=
i
k
Dw
ˆ
z
=0,
t
=0
0
/
2
e
ikx
dx
=4
0
V
A
.
11
In MHD, the incompressible vorticity equation
t
+
u
·
=
·
u
+
p
2
+
B
B
12
has an additional term involving the magnetic field. The ad-
ditional term implies that, even in the absence of baroclinic
generation, vortex lines are not necessarily material lines as
they are in hydrodynamics. For
t
0
+
,
w
ˆ
is smooth around
z
= 0, indicating that the circulation has been removed from
the density interface. Instead, circulation is carried by two
fronts that propagate at the local Alfvén speed in each fluid.
These fronts correspond to the locations where the Heaviside
functions change magnitude.
In the smooth regions of the flow, the vorticity is given
by
=
u
z
−
w
x
=
i
k
D
2
w
−
k
2
w
.
13
By substituting our solution for
w
into the above equation,
we find that the flow is irrotational upstream of the Alfvén
fronts in each fluid. Downstream of the Alfvén fronts, how-
ever, we find that the vorticity is nonzero. This is illustrated
in Fig.
4
, which shows the vorticity field for one particular
case. Note that the vorticity decays exponentially down-
stream of each Alfvén front.
C. Interface behavior
The value of
w
ˆ
z
,
t
at
z
= 0 is the growth rate of the
interface. From Fig.
3
, it can be seen that as
t
increases and
the Alfvén fronts propagate away from the interface, carrying
away the majority of the vorticity produced by the impulsive
acceleration, the growth rate of the interface decays to zero.
Thus the instability of the interface is suppressed and its
amplitude asymptotes to a constant value. For
t
→
, the in-
terface amplitude tends to
=
0
+
0
w
ˆ
0,
t
dt
=
0
1+
V
B
2
−
1
.
14
This shows that the change in interface amplitude is in-
versely proportional to
B
. Thus for
B
→
0,
→
, which is
in agreement with the result from hydrodynamic linear sta-
bility analysis.
8
Interestingly,
is independent of wave-
number. For finite times the interface amplitude is given by
t
=
0
+
0
t
w
ˆ
0,
T
dT
=
−
−
0
e
t
cos
t
,
15
where
and
are as defined in Eqs.
7
and
8
, respectively.
III. SIMULATION TECHNIQUES
A. Numerical method for linearized simulations
The linearized simulations presented here were carried
out using a method developed by Samtaney
9
for obtaining
numerical solutions to the linearized ideal MHD equations
when the base flow is temporally evolving. In this method,
the equations of compressible ideal MHD are specialized to
two dimensions,
x
and
z
. The solution is then assumed to
have the form
U
x
,
z
,
t
=
U
o
z
,
t
+
U
ˆ
z
,
t
exp
ikx
,
where
1,
U
o
z
,
t
is an unsteady one-dimensional base
flow, and
U
ˆ
z
,
t
exp
ikx
is the perturbation to the base
flow. A finite volume upwind approach is adopted to solve
for both the base flow and the perturbations. The equations
are integrated in time using a third order Total Variation Di-
minishing Runge–Kutta scheme and the fluxes are evaluated
using Roe’s method.
FIG. 4.
Color online
Negative vorticity field at
t
/
t
= 4 for
1
/
= 1.483 72,
2
/
= 4.431 59,
V
/
p
0
= 0.319 125,
0
/
= 0.007 992 76, and
= 16. Here
t
/
p
0
.
082102-4 Wheatley, Samtaney, and Pullin
Phys. Fluids
21
, 082102
2009
Downloaded 16 Oct 2009 to 131.215.220.165. Redistribution subject to AIP license or copyright; see http://pof.aip.org/pof/copyright.jsp
B. Setup for shock driven linearized simulations
Let us first consider the initial conditions for the base
flow of a SDL simulation. Prior to the interaction of the
incident shock with the density interface, which is unper-
turbed in the base flow, we designate the quiescent condi-
tions to the left
z
z
if
and right
z
z
if
of the interface as
states 1 and 2, respectively. The conditions downstream of
the incident shock are referred to as state 3. For given values
of
and the incident shock Mach number
M
, state 3 is
obtained from the normal shock relations for an ideal gas.
For the range of parameters considered herein, the interac-
tion of the incident shock with the interface generates a re-
flected shock and a transmitted shock. The conditions down-
stream of the reflected and transmitted shocks are referred to
as states 4 and 5, respectively. These states are shown in Fig.
5
, which shows the paths of the discontinuities in the base
flow in the
z
-
t
plane.
In the initial condition for the base flow, the incident
shock is represented by a sharp discontinuity located at
z
shock
=
z
if
−
/
5. The flow is initialized to state 3 for
z
z
shock
and state 1 for
z
z
shock
. For
z
z
shock
, the base flow
density is set to
o
x
,0
=
1
2
2
+
1
+
2
−
1
tanh
z
−
z
if
to represent the density interface. With this initial condition,
the base flow for a SDL simulation is the numerical solution
to the Riemann problem that arises from the interaction of a
shock with an unperturbed density interface. The only non-
zero perturbation at
t
= 0 is that in density, which approxi-
mates a delta function as follows:
ˆ
z
,0
=−2
R
−
L
exp
2
z
−
z
if
1 + exp
2
z
−
z
if
2
.
16
Here,
is a parameter that regularizes the sharp density in-
terface. It should be chosen judiciously such that the density
interface remains relatively sharp with its thickness generally
comparable to the numerical smearing of the shock front. As
long as this condition is met, the results do not change sig-
nificantly when
is varied. For our simulations,
=80
worked satisfactorily. Note that the initial perturbation am-
plitude of the interface,
0
, has been scaled out of the prob-
lem. For comparison with the results of NL simulations and
the ILT, the scaled perturbations from the linearized simula-
tions must be multiplied by
0
.
C. Setup for impulse driven linearized simulations
For IDL simulations, the base flow is initialized to state
4 for
z
z
if
and state 5 for
z
z
if
. These are the postshock
states from the Riemann problem described in Sec. III B. The
sharp interface between these two uniform states is approxi-
mated by hyperbolic tangent profile, with all quantities hav-
ing the form
q
z
,
t
=
1
2
q
L
+
q
R
+
q
R
−
q
L
tanh
z
−
z
if
,
where the subscripts
L
and
R
indicate values to the left and
right of the interface, respectively. The perturbations are ini-
tialized as described in Sec. III B. When scaling the pertur-
bations from IDL simulations for comparison with other re-
sults, the initial perturbation amplitude of the interface is
taken from the corresponding NL simulation immediately af-
ter the interface has been compressed by the passage of the
shock wave. These same postshock initial conditions are
used in the ILT.
D. Setup for nonlinear simulations
The NL simulations were carried out with a NL com-
pressible MHD solver developed by Samtaney. It uses the
eight-wave upwinding formulation of Powell
10
within an un-
split upwinding method.
11
The solenoidal property of the
magnetic field is enforced at each time step using a projec-
tion method. A constrained transport step is then used to
remove divergence modes with a centered finite difference
representation. This uses the formulation prescribed by
Toth.
12
The setup for the NL simulations is as described in
Sec. I, except with different parameter values.
E. Characterization of interface behavior
In the linearized simulations, the interface corresponds
to the location where the density perturbation is maximum.
The interfacial growth rate
̇
is approximated by the magni-
tude of the
z
-velocity perturbation at this location. The inter-
face amplitude
lin
is then computed by numerically integrat-
ing
̇
. At time step
N
,
lin
is given by
lin
=
0
+
n
=1
n
=
N
̇
n
t
n
,
where the subscript
n
denotes a quantity is evaluated at the
n
th time step. For simulations where the interface is shock
accelerated, the time origin of the interface amplitude histo-
ries is shifted to the point where the amplitude is minimum.
This is done to allow direct comparison to the results arising
from an impulsive acceleration at
t
= 0. The interface ampli-
tude histories obtained from the simulations must be quanti-
CD
CD
R
I
z
t
State 1
State 3
State 2
State 4
State 5
z
T
if
FIG. 5. Illustration of the base flow for SDL simulations in the
z
-
t
plane.
The lines shown are the paths of the discontinuities in the flow.
I
,
T
, and
R
designate the incident, transmitted, and reflected shocks, respectively, while
CD designates the contact discontinuity.
082102-5 The Richtmyer–Meshkov instability in MHD
Phys. Fluids
21
, 082102
2009
Downloaded 16 Oct 2009 to 131.215.220.165. Redistribution subject to AIP license or copyright; see http://pof.aip.org/pof/copyright.jsp