Linear Analyses of Magnetohydrodynamic
Richtmyer-Meshkov Instability in Cylindrical Geometry
Dissertation by
Abeer Habeeb Allah Bakhsh
In Partial Fulfillment of the Requirements
For the Degree of
Doctor of Philosophy
King Abdullah University of Science and Technology
Thuwal, Kingdom of Saudi Arabia
(April, 2018)
2
EXAMINATION COMMITTEE
The dissertation of Abeer Habeeb Allah Bakhsh is approved by the examination
committee
Committee Chairperson: Prof. Ravi Samtaney
Committee Members: Prof. David Keyes, Dr. Aslan Kasimov and Dr. Vincent
Wheatley
3
©April, 2018
Abeer Habeeb Allah Bakhsh
All Rights Reserved
4
ABSTRACT
Linear Analyses of Magnetohydrodynamic Richtmyer-Meshkov
Instability in Converging Geometry
Abeer Habeeb Allah Bakhsh
We investigate the Richtmyer-Meshkov instability (RMI) that occurs when an
incident shock impulsively accelerates the interface between two different fluids. RMI
is important in many technological applications such as Inertial Confinement Fusion
(ICF) and astrophysical phenomena such as supernovae. We consider RMI in the
presence of the magnetic field in converging geometry through both simulations and
analytical means in the framework of ideal magnetohydrodynamics (MHD). In this
thesis, we perform linear stability analyses via simulations in the cylindrical geometry,
which is of relevance to ICF. In converging geometry, RMI is usually followed by the
Rayleigh-Taylor instability (RTI). We show that the presence of a magnetic field
suppresses the instabilities. We study the influence of the strength of the magnetic
field, perturbation wavenumbers and other relevant parameters on the evolution of
the RM and RT instabilities. First, we perform linear stability simulations for a
single interface between two different fluids in which the magnetic field is normal to
the direction of the average motion of the density interface. The suppression of the
instabilities is most evident for large wavenumbers and relatively strong magnetic
fields strengths. The mechanism of suppression is the transport of vorticity away
from the density interface by two Alfvén fronts. Second, we examine the case of an
azimuthal magnetic field at the density interface. The most evident suppression of the
instability at the interface is for large wavenumbers and relatively strong magnetic
fields strengths. After the shock interacts with the interface, the emerging vorticity
5
breaks up into waves traveling parallel and anti-parallel to the magnetic field. The
interference as these waves propagate with alternating phase causing the perturbation
growth rate of the interface to oscillate in time. Finally, we propose incompressible
models for MHD RMI in the presence of normal or azimuthal magnetic field. The
linearized equations are solved numerically using inverse Laplace transform. The
incompressible models show that the magnetic field suppresses the RMI, and the
mechanism of this suppression depends on the orientation of the initially applied
magnetic field. The incompressible model agrees reasonably well with compressible
linear simulations.
6
DEDICATED TO MY PARENTS
This dissertation is dedicated to and in the memory of my mother, Refa’t Jihan
Amir Ahmad. She was passionate about knowledge and science. She tried her best to
make our lives comfortable and taught us how to continue our educational life. Once
she had a dream, and yes dreams can become true. This doctoral of philosophy degree
was her dream. Even during her illness, she pushed me to keep my spirits up and
finish what I started. May her soul rest in peace.
It is also dedicated to my father Habeeb Allah Bakhsh, who has stood by my side since
I decided to join KAUST and untill now. He is a great teacher who spent his youth
teaching many friends and young students. His pure soul gives us the strength to
continue our lives and chase our dreams. I am very proud of them as parents.
7
ACKNOWLEDGEMENTS
Sincere thanks to my supervisor, Prof. Ravi Samtaney, for his profound knowledge,
patience, and enthusiasm. The members of my dissertation committee, Prof. David
Keyes, Dr. Aslan Kasimov and Dr. Vincent Wheatley, have generously given their
time and expertise to better my work. I thank them for their contribution and support. I must acknowledge many friends, colleagues, teachers and others who assisted,
advised, and supported my research over the years. Especially, I am deeply grateful
and appreciate of Dr. Valeriya Ratushna, who was a postdoc in our group, my best
friend, Mrs. Nafisa Bakhsh and others for helping me emotionally and professionally
in this part of my life. Thanks to my group members for their cooperations, kindness,
and supports throughout the past few years. They have helped me keep perspective on
what is essential in life and shown me how to deal with reality.
Thanks to my parents along with my family for their support and prayers. They give
me unconditional love that I cannot forget. My thanks also must go to my friend and
husband Abdulmajid Turkistani who came into my life in the most critical part of my
journey and puts his effort to support me.
This work is supported by KAUST Office of Sponsored Research under Award URF/1/216201.
8
TABLE OF CONTENTS
Examination Committee Page
2
Copyright
3
Abstract
4
Dedicated to my Parents
6
Acknowledgements
7
List of Figures
11
List of Tables
Abbreviation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
List of Symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
16
17
18
1 Introduction
1.1 Background and Motivation . . . . . . . . . . . . . . . . . . . . . .
1.1.1 Importance of RMI and RTI in Inertial Confinement Fusion
1.2 The Richtmyer-Meshkov Instability . . . . . . . . . . . . . . . . . .
1.2.1 Various Stages of RMI . . . . . . . . . . . . . . . . . . . . .
1.2.2 Linear Stability Analysis of RMI . . . . . . . . . . . . . . .
1.2.3 Vortex Dynamical Interpretation of RMI . . . . . . . . . . .
1.2.4 RMI in the Presence of a Magnetic Field . . . . . . . . . . .
1.2.5 Planar Versus Converging Geometries . . . . . . . . . . . . .
1.2.6 Converging Geometry . . . . . . . . . . . . . . . . . . . . .
1.3 Objectives and Contributions . . . . . . . . . . . . . . . . . . . . .
1.3.1 Thesis Outline . . . . . . . . . . . . . . . . . . . . . . . . . .
19
19
20
23
25
25
27
28
35
37
45
46
.
.
.
.
.
.
.
.
.
.
.
2 Linear Simulations of the Cylindrical Richtmyer-Meshkov Instability
in the Presence of Radial Magnetic Field
49
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
2.2 Physical Setup and Governing Equations . . . . . . . . . . . . . . . . 51
9
2.3
2.2.1
2.2.2
2.2.3
Linear
2.3.1
2.3.2
2.3.3
2.3.4
2.3.5
Physical Setup . . . . . . . . . . . . . . . . . . .
Governing Equations . . . . . . . . . . . . . . . .
Linear Stability Analysis . . . . . . . . . . . . . .
Stability Results . . . . . . . . . . . . . . . . . . .
Base State Solution . . . . . . . . . . . . . . . .
Evolution of the Perturbations . . . . . . . . . . .
Effect of Converging Geometry . . . . . . . . . .
Effect of Magnetic Field Strength . . . . . . . . .
Comparison Between Linear and Nonlinear
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
Simulations
51
52
56
59
59
63
77
78
79
3 Linear Analysis of Converging Richtmyer-Meshkov Instability in the
Presence of an Azimuthal Magnetic Field
85
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
3.2 Physical Setup and Governing Equations . . . . . . . . . . . . . . . . 86
3.2.1 Physical Setup . . . . . . . . . . . . . . . . . . . . . . . . . . 86
3.2.2 Initial Conditions . . . . . . . . . . . . . . . . . . . . . . . . . 88
3.2.3 Governing Equations . . . . . . . . . . . . . . . . . . . . . . . 88
3.2.4 Linear Stability Analysis . . . . . . . . . . . . . . . . . . . . . 90
3.3 Results and Discussions . . . . . . . . . . . . . . . . . . . . . . . . . 92
3.3.1 Base State Solution . . . . . . . . . . . . . . . . . . . . . . . . 92
3.3.2 Evolution of the Perturbations . . . . . . . . . . . . . . . . . . 95
3.4 Comparison Between Linear and Nonlinear Simulations . . . 108
4 Incompressible Models of Magnetohydrodynamic Richtmyer-Meshkov
Instability in Cylindrical Geometry
110
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
4.2 Incompressible Model Formulation . . . . . . . . . . . . . . . . . . . 111
4.2.1 Normal Magnetic Field . . . . . . . . . . . . . . . . . . . . . . 111
4.2.2 Azimuthal Magnetic Field . . . . . . . . . . . . . . . . . . . . 112
4.3 Linearized Initial Value Problem . . . . . . . . . . . . . . . . . . . . . 116
4.3.1 Normal Magnetic Field . . . . . . . . . . . . . . . . . . . . . . 117
4.3.2 Azimuthal Magnetic Field . . . . . . . . . . . . . . . . . . . . 120
4.4 Numerical Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
4.5 Results and Discussions . . . . . . . . . . . . . . . . . . . . . . . . . 126
4.5.1 Verification of the Incompressible Model in Planar Geometry . 127
4.5.2 Normal Magnetic Field in Cylindrical Geometry . . . . . . . . 129
4.5.3 Azimuthal Magnetic Field in Cylindrical Geometry . . . . . . 130
10
4.5.4
Hydrodynamic case . . . . . . . . . . . . . . . . . . . . . . . . 132
5 Concluding Remarks
135
5.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
5.2 Future Research Work . . . . . . . . . . . . . . . . . . . . . . . . . . 138
References
139
Appendices
A.0.1 Convergence Test . . . . . . . . . . . . . . . . . . . . . . . . .
B.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
B.1.1 Numerical Method . . . . . . . . . . . . . . . . . . . . . . . .
148
150
151
155
11
LIST OF FIGURES
1.1
1.2
1.3
1.4
1.5
1.6
1.7
1.8
1.9
Schematics of indirect-drive (upper left) and direct-drive (upper right)
ICF. In both cases, a spherical capsule is prepared at t = 0 with a layer
of DT fuel on its inside surface. As the capsule surface absorbs energy
and ablates, pressure accelerates the shell of remaining ablator and DT
fuel inwardsan implosion. By the time the shell is at approximately
one-fifth of its initial radius it is traveling at a speed of many hundreds
of kilometers per second. By the time the implosion reaches minimum
radius, a hotspot of DT has formed, surrounded by colder and denser
DT fuel [1]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Reproduced from [2] (schematic before-and-after shock-interface interaction where the incident plane shock, the corrugated transmitted and
reflected shocks are shown. The materials are at rest before the shock
arrival). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Various Stages of RMI [3]. . . . . . . . . . . . . . . . . . . . . . . . .
Approximation of the asymptotic value of perturbed growth rate ȧ.
The initial parameters of the model: density ratio: ̺ = 1/8, 1/16, the
specific heats: γ = 3/2, velocity of the incident shock: U0 and fluids
initial pressures are equal [2]. . . . . . . . . . . . . . . . . . . . . . .
Vorticity deposition at a light/heavy interface (a) initial configuration [3] (b) clockwise rotation (b) counterclockwise rotation. . . . . .
Setup of the physical domain and boundary conditions for the RM
simulations. Lower end of density interface is initialized at x = 0 [4].
Density field at time t=3.33 for (HD), t=3.37 for (MHD) with nonzero
magnetic field β −1 = 0.5 are reflected about the x-axis. Simulation
parameters: Ms = 2, η = 3, θ = 45◦ [4]. . . . . . . . . . . . . . . . . .
Growth rate for Mach number M = 1.05, 1.25, 2.0 shocks for (a) the
HD case (b) the MHD case in planar geometry [5]. . . . . . . . . . .
Geometry for incompressible model problem, where ρ1 < ρ2 fluid densities, a contact discontinuity has the initial amplitude, wavelength and
acceleration, η0 , λ, V δ(t) respectively [6]. . . . . . . . . . . . . . . . .
22
24
25
26
29
30
31
32
34
12
1.10 Hydrodynamics simulation of the RTI [7]. . . . . . . . . . . . . . . .
Schematic diagram of the physical setup: (a) (r, θ)-section with constant z and (b) (r, z)-section with constant θ. Two fluids of densities
ρ1 and ρ2 are separated by a perturbed interface located at r = R0 . A
Chisnell-type incident shock is initially located at r = Rs , with Mach
number M . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Profiles of base state variables behind the shock front at t = 0 (a)
radial velocity ůr (b) pressure p̊. . . . . . . . . . . . . . . . . . . . . .
2.3 Density profiles of base state at t = 0.0 (initial conditions) and t =
0.552. IS denotes the incident shock, TS denotes the transmitted shock,
RS denotes the reflected shock and CD is the contact discontinuity. .
2.4 Spacetime diagram of density field of the base state. IS denotes the
incident shock, TS denotes the transmitted shock, RS denotes the reflected shock and CD is the contact discontinuity. . . . . . . . . . . .
2.5 Time history of radial velocity ůr at CD. . . . . . . . . . . . . . . . .
2.6 Time history of (a) acceleration of CD, (b) acceleration of CD with
approximate demarcation of RMI and RTI phases. . . . . . . . . . . .
2.7 Time history of scaled growth rate in hydrodynamics for azimuthal
perturbation wavenumber m = 16, 64, 128, 256. . . . . . . . . . . . . .
2.8 Profiles of perturbed variables in hydrodynamics at t = 0.552: (a) rcomponent of velocity, (b) θ-component of velocity, (c) pressure, (d)
z-component of vorticity. The azimuthal perturbation wavenumber is
m = 256. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.9 Time history of scaled growth rate in hydrodynamics for axial perturbation wavenumber k = 16, 64, 128, 256. . . . . . . . . . . . . . . . . .
2.10 Profiles of perturbed variables in hydrodynamics at t = 0.552: (a) rcomponent of velocity, (b) z-component of velocity, (c) pressure, (d)
θ-component of vorticity. The axial perturbation wavenumber is k = 256.
2.11 Time history of scaled growth rate in HD for 3-D perturbation wave
numbers k = m = 16, 64, 128, 256. . . . . . . . . . . . . . . . . . . . .
2.12 Time history of scaled growth rate in MHD for azimuthal perturbation
wavenumber m = 16, 64, 128, 256. The initial plasma beta is β = 4 at
the interface R0 = 1.0. . . . . . . . . . . . . . . . . . . . . . . . . . .
38
2.1
53
61
62
63
64
64
66
68
69
70
71
72
13
2.13 Profiles of perturbed variables in MHD (a) r-component of velocity,
(b) θ-component of velocity, (c) r-component of magnetic field, (d)
θ-component of magnetic field, (e) z-component of vorticity, (f) zcomponent of current density. The perturbation wavenumber is m =
256 and the initial plasma beta is β = 4 at the interface R0 = 1.0. . .
2.14 Time history of scaled growth rate in MHD for axial perturbation
wavenumber k = 16, 64, 128, 256. The initial plasma beta is β = 4
at the interface R0 = 1.0. . . . . . . . . . . . . . . . . . . . . . . . . .
2.15 Profiles of perturbed variables in MHD: (a) θ-component of vorticity,
(b) θ-component of current density. The perturbation wavenumber is
k = 256 and the initial plasma beta is β = 4 at the interface R0 = 1.0.
2.16 Time history of scaled growth rate for 3-D perturbations MHD for
different wavenumbers k = m = 256, m = 16 with k = 64, 128, 256,
k = 16 with m = 64, 128, 256, and k = m = 16 for β = 2. . . . . . . .
2.17 (a) Spacetime diagram of vorticity ω̂z for hydrodynamics. (b) Spacetime diagram of vorticity ω̂z for MHD. The perturbation wavenumber
is m = 256 and the plasma beta is β = 4 for MHD. . . . . . . . . . .
2.18 Time history of scaled growth rate of planar and cylindrical geometries
for k = 256 and m = 256 (a) hydrodynamics case and (b) MHD case
for beta β = 4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.19 Time history of scaled growth rate in MHD for plasma beta β =
∞(HD), 256, 64, 4, 2 at the interface for (a) m = 256, k = 0, (b)
m = 0, k = 256 and (c) β = ∞(HD), 16, 4, 2, 1 for m = 256, k = 256.
2.20 (a) Density and vorticity snapshot in nonlinear simulation. (b-d) Comparison of nonlinear growth rate between linear and nonlinear simulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1
3.2
Schematic diagram of the physical setup: Two fluids of densities ρ1
and ρ2 are separated by a perturbed interface located at r = R0 . A
fast magnetosonic MHD shock is initially located at r = Rs , with Mach
number Ms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Initial conditions of the base flow variables (a) Incident shock position
(b) Modified magnetic field (c) Pressure and (d) Radial velocity. . . .
74
75
76
77
81
82
83
84
87
89
14
3.3
(a) Density profiles of base state at t = 0.0 (initial conditions) and
t = 0.5 for β = 16. IS denotes the incident shock, TF denotes the transmitted fast MHD shock, RF denotes the reflected fast MHD shocks,
and CD is the contact discontinuity (b) Spacetime diagram of density field of the base state. IS denotes the incident shock, TF denotes
the transmitted fast MHD shock, RF denotes the reflected fast MHD
shock, and CD is the contact discontinuity. . . . . . . . . . . . . . . . 94
3.4 (a) Density profiles of base state at t = 0.0 (initial conditions) and
t = 0.2 for β = 4. IS denotes the incident shock, TF denotes the
transmitted fast MHD shock, RR denotes the reflected rarefaction fan,
and CD is the contact discontinuity (b) Spacetime diagram of density
field of the base state. IS denotes the incident shock, TF denotes the
transmitted fast MHD shock, RR denotes the rarefaction fan, and CD
is the contact discontinuity. . . . . . . . . . . . . . . . . . . . . . . . 96
3.5 Scaled ḣ for HD and m = 256 for A < 0 and A > 0. . . . . . . . . . . 98
3.6 Time history of growth rates for different β and fixed azimuthal wavenumber m = 256 (a) A > 0 (b) A < 0. . . . . . . . . . . . . . . . . . . . 100
3.7 2D Vorticity ωz in the interface for β = 16, m = 256 at different scaled
times chosen to correspond to a full oscillation period of the growth
rate. The arrows in the growth rate plot (a) Indicate the times chosen. 103
3.8 Time history of the growth rates for different m and β = 4 (a) A > 0
(b) A < 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
3.9 Time history of growth rates for high β values, and fixed azimuthal
wavenumber m = 256 (a) A > 0 (b) A < 0. . . . . . . . . . . . . . . 107
3.10 Comparison of growth rate between linear and nonlinear simulations
(a) m = 16, β = 16 (b) m = 128, β = 16 (c)m = 16, β = 4 (d)
m = 128, β = 16. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
4.1
4.2
4.3
(a) Initial condition for the shock-driven compressible MHD RMI (b)
Initial condition for the corresponding impulse driven, incompressible
model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
(a) Initial condition for the shock-driven compressible MHD RMI (b)
Initial condition for the corresponding impulse driven incompressible
model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
Profiles of velocity ŵ (a) the analytical solution (b) using the numerical
method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
15
4.4
4.5
4.6
4.7
cylind incomp. model . . . . . . . . . . . . . . . . . . . . . . . . . . .
Comparison between the incompressible model and the linear simulations for growth rate (a) m = 64, β = 16 (b) m = 128, β = 16. . . . .
Comparison between impulsive model and numerical simulations for
scald growth rate (a) m = 32, β = 16 (b) m = 128, β = 16. . . . . . .
Scaled ḣ for HD and m = 128 for the model and numerical simulation.
130
131
133
134
A.1 Time history of unscaled growth rate of the shocked interface for β = 16
wavenumber m = 128 using different mesh sizes: 800, 1600, 3200, 6400, 12800
zones in the radial domain and (b) Relative error of unscaled growth
rate at time t = 0.3. . . . . . . . . . . . . . . . . . . . . . . . . . . . 152
B.1 cylind incomp. model . . . . . . . . . . . . . . . . . . . . . . . . . . . 156
16
LIST OF TABLES
Abbreviation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
List of Symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
17
18
1.1
1.2
1.3
42
43
44
Summary of previous works that studied RMI . . . . . . . . . . . . .
Continue: Summary of previous works that studied RMI . . . . . . .
Continue: Summary of previous works that studied RMI . . . . . . .
17
LIST OF ABBREVIATION
Symbol
Meaning
RMI
Richtmyer-Meshkov instability
RTI
Rayleigh-Taylor instability
MHD
Magnetohydrodynamics
HD
Hydrodynamics
CD
Contact Discontinuity
IS
Incident Shock
RS
Reflected Shock
TS
Transmitted Shock
RR
Reflected Rarefaction
LT
Laplace Transform
ILT
Inverse Laplace Transform
18
LIST OF SYMBOLS
Symbol
Meaning
ρ
Density
Ms
Shock Mach number
ur , uθ , uz
The velocity components in cylindrical geometry
B r , Bθ , Bz
The magnetic filed components in cylindrical geometry
e
Total energy
pt
Total pressure
m, k
Azimuthal and axial wavenumbers
ȧ, ḣ
The growth rate of the interface
a0 , a+
0
Initial amplitude, amplitude after impulse
A, A+
Atwood number, post shock Atwood number
i
Complex number
λθ
Wave length
η0
Initial amplitude
W
Solution vector
F, G, H
Flux vectors
M, N, Q
Jacobian matrices
Br∗ , Bθ∗
Curl free part of magnetic field
Br1 , Bθ1
Corrected magnetic field
Im(x)
Imaginary part of x
H(x)
Heaviside function
U, V, Hr , Hθ , P
The temporal LT of ûr , ûθ , B̂r , B̂θ , p̂t
∆V
The impulse velocity
s
The Laplace variable
N
Number of points
19
Chapter 1
Introduction
In this chapter, we give an adequate introduction to the Richtmyer Meshkov instability and state the objectives of this dissertation. Various results for the hydrodynamic
Richtmyer Meshkov instability are reviewed, along with relevant results from magnetohydrodynamics. Literature specific to each of the subtopics is addressed herein,
but the investigation as a whole is reviewed in the relevant chapters. Finally, the
problems that are addressed in this investigation are outlined and justified.
1.1
Background and Motivation
The Richtmyer-Meshkov instability (RMI) arises when a shock wave interacts with
an interface separating two different fluids. The most significant technological application where RMI is important is in Inertial Confinement Fusion (ICF), which has
been a major motivation for the study of shock-accelerated interfaces [3]. Another
hydrodynamical instability is the Rayleigh-Taylor Instability (RTI) that emerges at
the interface between two fluids of different densities when the lighter fluid pushes
the heavier fluid due to acceleration such as gravity (Sharp [8]). RTI is also important in the context of ICF, where it arises mostly due to the effects of converging
geometrical configurations. RMI is also important in astrophysical phenomena; it
is used to account the lack of stratification in the products of supernova 1987A [9].
RMI is beneficial for supersonic and hypersonic air-breathing engines, where it may
be useful to enhance the mixing of fuel and air [10]. It also arises in many combus-
20
tion systems where shock wave interacts with a flame and the resulting instability
is important in the deflagration-detonation transition [11]. High-speed combustion
applications also are another field where RMI has its importance. For example, in a
scramjet engine, this instability arises due to the interaction of a shock wave with a
flame front, which can greatly enhance the fuel-oxidizer mixing rate [12]. Next, we
discuss the importance of RMI and RTI in ICF followed by a review of both of these
instabilities.
1.1.1
Importance of RMI and RTI in Inertial Confinement
Fusion
With each passing year, the global demand for alternative energy sources to fossil
fuels becomes more intense, particularly for sources that are carbon-free. In the very
long term, mankind’s only hope for energy may be fusion energy. At the heart of
fusion lies Einstein’s famous E = mc2 formula. Fusion of heavy isotopes of hydrogen,
namely deuterium and tritium releases a vast amount of energy (one example of a
fusion reaction is D +T → n +4 He +17.6MeV). Furthermore, there is enough fusion
fuel to last tens of thousands of years. Broadly, fusion energy research is divided
into MFE (magnetic fusion energy) and IFE (inertial fusion energy). In this thesis,
we focus on IFE. In the 1960s, work commenced on the second major approach,
now generally known as inertial confinement fusion, or ICF. The main idea behind
ICF was to attain very high densities, precisely the converse of the conditions sought
by MFE. The goal of ICF is to produce fusion reactions at sufficiently high matter
densities such that a small mass of the fuel maybe confined only by its own inertia.
Each fuel ignition is akin to a small explosion (releasing energy of approximately
hundreds of pounds of high explosive). A fusion reactor would extract pulsed energy,
most of it emerging as neutrons, which would be used to generate steam and breed
tritium fuel. No attempt is made to constrain the burning fuel, and in that sense
21
one might think of ICF as unconfined fusion. By imploding the fuel to very high
densities (approximately two orders of magnitude denser than lead) before ignition,
the time scales for fusion material dispersal are longer than the time scales required
to “burn”the fuel. The only force holding the fuel in place during heating is its own
inertia: hence, inertial confinement fusion [13]. Today in ICF, powerful lasers impinge
on a deuterium-tritium (DT) filled target, either directly (the so-called direct drive
method) or indirectly (the so-called indirect drive). The most prominent facility for
indirect drive is the National Ignition Facility (NIF) commissioned at the Lawrence
Livermore National Laboratory in the United States in 2009. NIF comprises of 192
laser beams with a total power exceeding 500TW. These lasers typically deposit
approximately 1.9MJ of energy on a hollow gold chamber target dubbed the hohlraum
(See Fig. 1.1). Soft X-rays generated by the laser interaction with the hohlraum
impinge radially on a spherical DT pellet ( 2mm in radius) surrounded by an ablative
shell (approximately 150 kJ of energy is absorbed by the outer layers of the target)
generating about 100Mbar of pressure in the ablator. The radiation pressure blows
off the outer ablative shell sending a converging shock into the DT pellet. If the
physical processes were perfectly symmetric with no perturbations, the resulting hot
spot in the center of the DT pellet would exceed 4KeV (1eV = 11,604 K), followed
by fusion of the D-T fuel. The National Ignition Campaign (NIC) followed the NIF
commissioning for the period 2009-2012 with the hope of igniting the target and
creating sustained fusion. While initial experiments demonstrated that conditions in
the facility can meet the simultaneous requirements on laser coupling and capsule
implosion symmetry for future ignition shots [14]. Unfortunately, the NIC did not
achieve all of its stated goals. Lindl et al.[15] provides a comprehensive review of the
NIC. In the final section entitled “The Path Forward”of this paper, the authors write:
“Current evidence points to low-mode asymmetry and hydrodynamic instability as
key areas of research to improve the performance of ignition experiments on the NIF
22
Figure 1.1: Schematics of indirect-drive (upper left) and direct-drive (upper right)
ICF. In both cases, a spherical capsule is prepared at t = 0 with a layer of DT fuel
on its inside surface. As the capsule surface absorbs energy and ablates, pressure
accelerates the shell of remaining ablator and DT fuel inwardsan implosion. By the
time the shell is at approximately one-fifth of its initial radius it is traveling at a
speed of many hundreds of kilometers per second. By the time the implosion reaches
minimum radius, a hotspot of DT has formed, surrounded by colder and denser DT
fuel [1].
23
and are a central focus of the Ignition Program going forward.”It is clear that a
key bottleneck towards achieving fusion is hydrodynamic instabilities. One of the
main objectives of this dissertation is to examine the RMI in a cylindrical geometry
within the context of linear stability analyses. However, before we describe the main
objectives of our investigations, we begin with a review of RMI particularly in the
context of converging geometry and magnetohydrodynamics (MHD). Because of the
secondary instability, i.e., RTI that arises in ICF, we briefly review RTI after the
review of RMI.
1.2
The Richtmyer-Meshkov Instability
The Richtmyer-Meshkov instability (RMI) occurs when a corrugated interface between two fluids of different densities is impulsively accelerated. Typically, the impulse is provided by a shock (aka the “incident shock”) striking the interface. At the
interface the shock undergoes a refraction. Generally, if the density ratio is greater
ρ2 − ρ1
than unity (positive Atwood number, A =
, where ρ1 and ρ2 are the densities
ρ2 + ρ 1
of first and second fluids), then the shock bifurcates into transmitted and reflected
shocks, as shown in Fig. 1.2. A theoretical prediction of Taylor instability in a shock
acceleration of compressible fluids was provided by Richtmyer [2] who developed the
linear theory of RMI. Meshkov [16] performed the first shock-tube experiments with
a sinusoidally-shaped thin membrane a perturbed interface between different kinds
of gases: the experimental results were in agreement with Richtmyer’s theory and
hence this instability came to be known as the Richtmyer-Meshkov instability. We
next discuss the various stages in the development of the instability, followed by a
discussion of linear stability analysis of RMI. We also provide a description of the
physical mechanism that produces this instability using the dynamics of the baroclinic vorticity, which is generated by the misalignment of the pressure gradient of
the shock wave and the density gradient at the perturbed interface [17].
24
Figure 1.2: Reproduced from [2] (schematic before-and-after shock-interface interaction where the incident plane shock, the corrugated transmitted and reflected shocks
are shown. The materials are at rest before the shock arrival).
25
Figure 1.3: Various Stages of RMI [3].
1.2.1
Various Stages of RMI
Linear stability analysis of several flows is considered an important precursor to nonlinear stages in the flow development, which in turn is an important stage before
turbulence ensues. Indeed, hydrodynamic RM flows follow a similar development:
very small amplitude perturbations grow linearly until they are comparable to the
perturbation wavelength, after which nonlinear mode coupling ensues followed by
full-blown turbulent mixing. The evolving of RM from a linear stage to the nonlinear one is shown in Fig.1.3. In the early times and for small initial amplitudes, the
growth of the perturbation on the interface is considered linear in time (Fig.1.3 (a,
b)). However, the increase in amplitude leads to crests and troughs of the perturbations evolve asymmetrically. As a result, spikes of heavy fluid penetrate into the light
one and bubbles of light fluid rise into the heavy one (Fig. 1.3 c). Next, The KelvinHelmholtz instability becomes important, causing the roll-up “mushrooming”of the
spike (Fig.1.3 d) and the appearance of smaller scales. The final stage when a turbulent mixing zone (TMZ) develops between the two fluids (Fig.1.3(e)) [3].
1.2.2
Linear Stability Analysis of RMI
It is pointed out that RMI is not a classical hydrodynamic instability with modal
analysis, because the perturbations grow linearly (or linearly in an averaged sense)
rather than exponentially. The original analysis by Richtmyer was performed on a
26
vertical interface with a single mode sinusoidal corrugation, between two different
ideal gases. A planar incident shock interacts with the interface and bifurcates into
transmitted and reflected shocks. Figure 1.2 shows the schematic of before-andafter interaction between an incident shock and the interface. In this analysis, the
Euler equations of motion were linearized and an initial boundary value problem was
derived. After the interaction of an incident shock with a contact discontinuity a
reflected and a transmitted shock are generated and considered as a boundaries for
the problem. The base state in this analysis is a self-similar solution involving a
reflected shock, a transmitted shock, and a contact discontinuity that translates at
an average constant speed due to the impulse by the incident shock. The linearized
compressible inviscid equations were solved by Richtmyer, in which the perturbation
amplitude growth rate undergoes a short transient before settling into a oscillatory
state around a constant asymptotic value. The perturbations initially increases to
specific value then it increases linearly in time (Fig. 1.4).
Figure 1.4: Approximation of the asymptotic value of perturbed growth rate ȧ. The
initial parameters of the model: density ratio: ̺ = 1/8, 1/16, the specific heats:
γ = 3/2, velocity of the incident shock: U0 and fluids initial pressures are equal [2].
27
Impulse Model
Richtmyer proposed an incompressible model for an impulsively accelerated perturbed
interface. Immediately after the passage of this impulse, the initial amplitude a0 of
the interface is changed to a+
0 which is the amplitude after the impulse. Also, the
fluids densities ρ1 , ρ2 are changed, which affects the density ratio at the interface,
Atwood number. Its new value after the impulse is notated as A+ . The asymptotic
growth rate of the perturbed interface will have a rate of change
da
dt
+
= a+
0 ∆wA ,
(1.1)
∞
where ∆w is the impulsive change in the interface velocity and ∞ denotes the asymptotic value.
1.2.3
Vortex Dynamical Interpretation of RMI
As we mentioned earlier, the interaction between the incident shock and a perturbed
interface generates the instability which are physically can be explained by the generation of the baroclinic vorticity. The general vorticity equation for a compressible,
inviscid flow with variable fluid properties can be written as
∂ω
∇ρ × ∇p
,
+ (u.∇)ω = (ω.∇)u − ω(∇.u) +
∂t
ρ2
(1.2)
where ω is the vorticity, u is the flow velocity, ρ and p are the local density and
pressure [18]. The first term on the right-hand side of the equation is due to vortex
stretching due to the flow velocity gradients, which vanishes in two dimensions, second
term is the bulk dilatation, the third is the baroclinic term. As seen from the equation,
the baroclinic term is generated by the misalignment of the pressure gradient across
the shock with the density gradient across the interface. Figure 1.5 (a) shows light
28
and heavy fluids ρ1 < ρ2 separated by a sinusoidal interface. As the incident shock
interacts with the interface, clockwise rotation is generated. The interface continues
rotating which extend the positive growth of the perturbation (Fig 1.5 (b)). For the
negative Atwood, ratio interface A < 0, where the shock moves from heavy to light
fluids the amplitude will decrease in time goes to zero and continue decreasing as
negative values which generate counterclockwise rotation. This is referred as phase
change of perturbations (Fig 1.5 (c)).
1.2.4
RMI in the Presence of a Magnetic Field
The first work to investigate the RMI in the presence of a magnetic field was by Samtaney [4]. Subsequently, Wheatley and co-workers have carried out simulation and
modeling studies in planar geometry (elaborated further in the next two sub-sections)
[6, 19, 20, 21]. Sano et al. [22] performed a numerical simulation on compressible Ideal
MHD equations. Their work concentrated on the evolution of an ambient magnetic
field affected by the nonlinear motions of the RMI. It was found that due to the fluid
motions associated with the RMI, the ambient magnetic field strength was amplified
efficiently at the perturbed interface between the two fluids. Hohenberger et al. [23]
applied experimentally laser-driven magnetic-flux compression to (ICF) targets to enhance the implosion performance. This technique is used to reach the strong magnetic
fields that are required to magnetize a hot spot in the target.
Non-zero Normal Component
In the planar geometry, Samtaney [4] showed that the presence of the magnetic field
(with a net non-zero component of the field normal to the interface) inhibits the
growth of the RMI. The mathematical model of two different fluids was described
by ideal MHD equations with a shock propagating from left to right and an oblique
interface as shown in Fig. 1.6. The principal parameters of these problems are the
29
(a) Hydrodynamic case
Light
Heavy
(b) Positive growth of perturbation
Heavy
Light
(c) Phase change of perturbation
Figure 1.5: Vorticity deposition at a light/heavy interface (a) initial configuration [3]
(b) clockwise rotation (b) counterclockwise rotation.
30
strength of the incident shock Ms , the Atwood ratio across the interface A = η, the
angle between the incident shock and the density interface θ, and the non-dimensional
strength of the magnetic field written as β = 2p0 /B02 , where p0 is the initial pressure,
B0 is the magnitude of a magnetic field. A density field is shown in Fig. 1.7 for the
hydrodynamic RMI (top) as well as the evolution in the presence of the magnetic field
(bottom). During the early time after the shock-interface interaction, the interface
is compressed by the incident shock, and baroclinic vorticity generation takes place.
As the simulation time increases, the absence of the magnetic field β −1 = 0, the
interface acts like a vortex layer and rolls up, as expected for the standard RMI. On
the other hand, the presence of the magnetic field β −1 = 0.5 remains the interface
smoothly, and no evidence of roll-up is observed. The vortex sheet splits, and two
of the MHD shocks transport the baroclinically generated vorticity away from the
contact discontinuity, which results in the suppression of RMI. Subsequent work by
Figure 1.6: Setup of the physical domain and boundary conditions for the RM simulations. Lower end of density interface is initialized at x = 0 [4].
Wheatley et al.[19] analyzed in detail the early refraction phase of the incident shock
in MHD. Furthermore, Wheatley et al. [6] proposed an incompressible model of the
RMI in MHD, and found that the initial growth rate of the interface is unaffected by
the presence of a magnetic field but the interface amplitude asymptotes to a constant
31
Figure 1.7: Density field at time t=3.33 for (HD), t=3.37 for (MHD) with nonzero
magnetic field β −1 = 0.5 are reflected about the x-axis. Simulation parameters:
Ms = 2, η = 3, θ = 45◦ [4].
value for a finite magnetic field, resulting in the suppression of the interface instability.
In their incompressible model, the vorticity is transported away from the interface by
Alfvén fronts.
Linear simulations for hydrodynamic (HD) and the MHD case were investigated
by Samtaney [5]. The growth rate of the perturbed interface for each case is shown in
Figure 1.8. In these figures, the time t = 0 is when the incident shock has completely
traversed the interface. The growth rate plots in the gas dynamics has the fast initial
growth followed by gradual decaying oscillations about an asymptotic value. While
for MHD case the growth rate increase initially, but eventually decays and oscillates
about zero. The stronger incident shocks shows larger amplitude of these oscillations.
32
(a) Hydrodynamic case
(b) MHD case
Figure 1.8: Growth rate for Mach number M = 1.05, 1.25, 2.0 shocks for (a) the HD
case (b) the MHD case in planar geometry [5].
33
Tangential Magnetic Field Case
The dynamics of the interface instability is quite different if the magnetic field is
initially tangential to the interface. If the direction of the magnetic field is parallel
to the shock surface, the stabilizing mechanism is the Lorentz force [24]. The case
of a tangential field was investigated in detail by Wheatley et al.[21]. They found
that the baroclinic vorticity on the interface immediately after the shock impulse
breaks up into waves traveling parallel and anti-parallel to the magnetic field, which
transport the vorticity. The interference of these waves as they propagate causes the
perturbation amplitude of the interface to oscillate in time.
Linear Analysis of RMI in Incompressible MHD
Richtmyer realized that as the reflected and transmitted shocks move away from the
interface, the flow becomes nearly incompressible. This formed the justification for
his incompressible impulse model which agrees well with linear compressible results
over a large extent of the parameter space. In a similar spirit, Wheatley et al. [6]
examined the behavior of an impulsively accelerated perturbed interface separating
incompressible conducting fluids of different densities, in the presence of a normal
magnetic field. Their incompressible model is characterized by the density ratio at the
interface ρ2 /ρ1 , a magnetic field B parallel to the acceleration with the strength β =
2p0 /B 2 , where p0 is the initial pressure, B is the magnitude of the applied magnetic
field, the ratio of the specific heats γ, and the normalized magnitude of the impulse
p
V ρ1 /p0 , (see Fig. 1.9). They linearized the initial value problem about a base flow
which results from the impulsive acceleration of an unperturbed interface. The initial
conditions are given at t = 0− which is before the impulsive acceleration, and there are
no perturbations in the magnetic field. They applied the temporal Laplace transforms
to the perturbed equations and derived a fourth order ordinary differential equation,
along with appropriate interface conditions. The boundary conditions assumed that
34
as |z|→ ∞ the perturbations must be bounded, there are no incoming waves from
z = ∓∞. Across the interface, the jump condition is obtained by integrating the
momentum equation in z direction. Finally, the inverse of Laplace transforms applied
to find the solution of the velocity components in each fluid. They found that the
solutions do not grow exponentially in time. It was found that the initial growth
rate of the interface at t = 0+ is unaffected by the presence of a magnetic field. But
for the later time, the authors found that there are two fronts of the velocity arise
naturally in the solution, propagating at the local Alfvén speed. This incompressible
Figure 1.9: Geometry for incompressible model problem, where ρ1 < ρ2 fluid densities,
a contact discontinuity has the initial amplitude, wavelength and acceleration, η0 , λ,
V δ(t) respectively [6].
model agreed well with nonlinear simulations as well as linear compressible results.
For the case of the tangential magnetic field, Wheatley et al. [21] proposed an
incompressible model in a manner similar to the one for the normal field component
case. In this case, however, the resulting Laplace transformed linear equations resulted in a second order ODE. Applying appropriate conditions, the solution showed
similar vorticity dynamics as in the nonlinear compressible case. Good agreement was
35
reported between the incompressible model and nonlinear compressible simulations
over a wide range of parameters.
1.2.5
Planar Versus Converging Geometries
RMI has attracted the attention of scientists worldwide because of the key role it
plays in multiple scientific areas. Linear and nonlinear simulations or experimental
studies of RMI have been performed in planar or converging geometries. We have
reviewed some of the basics of RMI above along with a few key studies. Here we
review some others: Tables 1.1, 1.2 and 1.3 contain a summary of these different
studies.
Planar Geometry
In this section, we will present some of the important studies that have been done
previously in the planar geometry. Meyer and Blewett [25] analyze numerically the
case of a rarefaction wave which is reflected after the incident shock interacts with the
material surface. Fraley [26] found that the initially deposition of the vorticity on the
interface causes pressure perturbations across the interface. Zhang and Sohn [27] presented an analytical solution for the nonlinear theory for RMI for compressible fluids.
They discussed the case when the reflected wave was a shock. This method showed
that the decay of the growth rate of the instability is due to the effect of nonlinearity,
rather than the effect of compressibility. The effect of compressibility dominates in
the early stage of development of the instability, while the effect of nonlinearity dominates in the intermediate and late stages of instabilities. Once the shock moves away
from the material interface, the fluids can be treated as approximately incompressible. Yang et al. [10] solved the linearized equations numerically for two cases, when
an incident shock moves from a heavy to a light one. They determined the nature of
the reflected wave, i.e., whether it is a rarefaction or a shock. Velikovich [28] obtained
36
an analytical solution for the reflected rarefaction case using power series expansions.
For strong shocks, Samtaney and Meiron [29] investigated a numerical model for the
compressible Euler equations. The impulse model and Richtmyer’s linear theory are
extended to include equilibrium dissociation chemistry. They discussed the case of
frozen chemistry, where the reaction rates are zero, and found that the growth rate
of impulse model agrees with the late-time asymptotic growth rate from the linear
theory. They also discussed both cases of the reflected wave when it becomes a reflected shock or rarefaction wave. Many experiments in shock tubes have been used
to study RMI. Brouillette [3] reviewed the fundamental physical processes underlying
the onset and development of the RMI in simple geometries, due to theoretical results
along with their experimental and numerical validation. One of these experiments
was carried by Jones and Jacobs [30]. In their work, the growth rate showed similar
results as the linear theory of Richtmyer. A theoretical framework to study linear and
nonlinear (RMI) was presented by Nishihara et al. [31], using an analytical model to
determine the asymptotic linear growth rate in planar and cylindrical geometries for
incompressible and irrotational fluids. In the linear theory of RMI, they described the
deposition of vorticity through the interaction between shock fronts and an interface
and gave an exact formula for the asymptotic growth rate. The effect of transverse
magnetic field, oriented parallel to the interface, recently has been investigated analytically and numerically [32, 24, 33, 21]. In fact, in a super-conductive plasma,
the existence of a transverse magnetic field behaves in a fashion somewhat similar to
surface tension on the interface. The transverse magnetic field inhibits the growth of
surface perturbations and as a result, stabilizes the instability [32]. Cao et al. [24]
studied this case in the framework of incompressible media through an analytical
development of a magnetized impulsive model. They also found an inhibition of the
RMI growth with oscillations of the perturbation amplitude. Wheatley et al. [21]
investigated nonlinear compressible simulations and incompressible model of MHD-
37
RMI for the case of a transverse magnetic field, i.e. there was no normal component
of the magnetic field at the interface and found oscillatory solutions of the growth
rates with the RMI still being suppressed. These studies in planar geometry clarified
some of the RMI mechanisms. Since our interest is to investigate the linear stability
analysis in converging geometry, we will review some of the former works in the next
section.
1.2.6
Converging Geometry
The most important applications of RMI is the ICF which typically occurs in converging geometry. In this geometry, Rayleigh-Taylor instability (RTI) occurs as a
secondary instability due to the acceleration experienced by the interface as it moves
radially inwards. The RTI emerges at the interface between two fluids of different
densities when the lighter fluid is pushing the heavier fluid due to acceleration such
as gravity (Sharp [8]). Figure 1.10 shows a typical hydrodynamic simulation of the
RTI. The figure shows the starting of formation of a “mushroom cap”which are more
visible at a later stage in the third and fourth frame in the sequence [7].
Previous Studies in Converging Geometry
It is of interest to examine the RM and RTI instabilities in a converging geometry,
especially in the presence of a magnetic field, and investigate the physical mechanisms of suppression, if any. In this section, since previous studies in the converging
geometry are discussed. Bell [34] was the first to consider the RTI in curved geometry. An incompressible fluid-void system was used to derive equations that govern
the growth of linearized perturbations. He also generalized the results to include a
uniform compressible fluid. Plesset [35] studied the linear stability of RTI in spherical
symmetry, for the interface between two different incompressible immiscible fluids in
normal motion. Plesset’s linear analysis was extended to the spherical and cylindrical
38
Figure 1.10: Hydrodynamics simulation of the RTI [7].
system composing ˝Nconcentric fluid shells by Mikaelian [36]. Mikaelian treated the
linear regime by assuming that perturbations are small. In his work, a formulated
model to estimate the mixing layer width between different fluids was presented. Also,
for the cylindrical case, Mikaelian discussed the freeze-out phenomena1 which leads
to a zero asymptotic growth rate. Hosseini and Takayama [37] performed an experimental study of RMI, which was induced by cylindrical shock waves propagating
across cylindrical interfaces. The results showed that after shock wave interaction,
the re-shocked cylindrical interfaces had a high growth rate of the turbulent mixing
zone. One of the early investigations of RMI in cylindrical geometry was carried
out by Zhang and Graham[38], in which they examined both the imploding (shock
converging towards cylinder axis) case and the exploding (shock moving radially outwards) case for positive and negative Atwood ratio interfaces. Furthermore, they also
1
Freez-out phenomena of the growth rate: where the pressure perturbations coming from the
stable shock fronts affect the growth, typically slowing it down, even a full cancellation is possible
39
examined the growth rate of small amplitude perturbations via “direct numerical simulations”rather than solving the linearized Euler equations. Zhang & Graham write
in their paper: “The equations and boundary conditions of the linear theory of RM
instability in curved geometry are considerably more complicated than the equations
and boundary conditions of the linear theory of RM instability in plane geometry,
and to date there is no linear theory for the curved case”. It is not uncommon in the
RMI literature to encounter linear stability investigations via small amplitude nonlinear simulations. Our approach is different in that we actually solve the linearized
equations of motion. Lombardini [39] investigated the RMI in converging geometries
analytically and numerically, the linear regime is covered as it is the onset to nonlinear stages of the perturbation growth. Another study of RMI in incompressible flow
was conducted by Lombardini and Pullin [40] for cylindrical and spherical geometries.
They investigated the linear regime of an interface between two fluids following the
passage of imploding or exploding shock waves. This was done by using an extension
of the results of Mikaelian [41] to three-dimensional azimuthal and axial perturbations in cylindrical perturbations and numerical Euler-based simulations. They used
the approximate analytical solution of Chisnell [42] who provided a simple analytic
description of the flow behind collapsing shock waves. They developed a theory for
an asymptotic growth rate of density interface. The growth amplitude and rate of the
spike (heavy fluid penetrating into light fluid) and the bubble (light fluid penetrating
into heavy fluid) had been computed. They showed that in curved geometry, the perturbation growth at the interface is not only due to the initial impulsive deposition
of the vorticity but also due to radial motion of the interface. In the numerical simulations of Yu and Livescu [43], the linear stability analysis of the RTI between two
ideal inviscid immiscible2 compressible fluids was performed in cylindrical geometry.
It was shown that the growth rate increased as the ratio of specific heats decreased,
2
Immiscible fluids: Do not form a homogeneous mixture when added together.
40
but it decreased as the Mach number increased. A Front Tracking method3 to solve
RTI and RMI in cylindrical and spherical geometries is used by Glimm et al. [44].
They assumed that the fluid is axi-symmetric, so the dynamic change of the fluid
interface can be described in a two-dimensional plane. The authors showed that the
azimuthal asymmetry is one of the physical properties of spherical RMI simulations.
Many numerical studies were done to examine the converging shocks which used
to initiate the RMI. Pullin et al. [45] examined the behavior of a converging cylindrical fast magnetosonic MHD shock onto a time-wise constant line current within ideal
MHD using Whitham’s theory of geometrical shock dynamics (GSD) [46], where the
shock is moving into a region with a spatially varying magnetic field. This shock can
weaken in terms of pressure ratio and Mach number as it reaches the origin, contrary
to the expected behavior under similar conditions of a gas-dynamic shock [42]. Recently, the effect of applying an external magnetic field on cylindrical and spherical
implosions in ideal MHD flows was presented by Mostert et al. [47]. The effects of
magnetic fields were studied with the aim of reducing the heat loss caused by RM and
RT instabilities. It was made evident that the azimuthal, and spherical symmetry
may be violated in cylindrical and spherical geometry in the presence of a magnetic
field, respectively. They found that the saddle configuration of magnetic field yields
the best symmetry during the implosion. Mostert et al.[48] investigated the behavior
of a cylindrical fast MHD shock wave that collapses onto an axial line where a timedependent current generates an azimuthal magnetic field that has a power law time
behavior. The time variation of this magnetic field is tuned in a manner such that
it becomes zero at the instant where the shock reaches the axis. In these aforementioned studies, the emphasis was on the behavior of the shock as it collapses on to
the central axis. Mostert et al. [49] carried out an extensive investigation of nonlin3
A Front Tracking method is an adaptive computational method in which a lower dimensional
domain is used. It is similar to a moving grid so that it can follow the dynamical evolution of
distinguished waves in a fluid flow.
41
ear RMI in cylindrical and spherical geometry for various initial seed magnetic field
configurations. All cases showed suppression of the RM and RT instabilities.
Paper
Bell and Taylor [34]
Plesset [35]
Richtmyer [2]
Meshkov [16]
Meyer
and
Blewett [25]
Fraley [26]
Mikaelian [50]
Velikovich [28]
Jones and Jacobs [30]
Samtaney
and
Meiron [29]
Zhang
and
Sohn [27]
Zhang and Graham [38]
Holmes et al. [51]
Geometry
Cylindrical,
spherical
Spherical
Initializing the instability
Steady acceleration
Flow
Incompressible,
compressible
Incompressible
Equations
Linear
Linear
HD
Planar
Planar shock, impulse acceleration
Compressible
Linear, nonlinear
HD
Planar
Planar shock
Compressible
Nonlinear
HD
Planar
Planar shock
Compressible
Linear
HD
Planar
Planar
Compressible
Incompressible/
compressible
Incompressible
Compressible
Linear, nonlinear
Nonlinear
HD
HD
Spherical shells
Planar
Planar shock
Steady acceleration/ planar
shock
Impulse acceleration
Planar shock
Linear, nonlinear
Linear
HD
HD
Planar
Planar
Planar shock
Planar shock
Compressible
Incompressible
Linear
Linear, nonlinear
HD
HD
Planar
Planar shock
Compressible
Linear, nonlinear
HD
Planar
Planar shock
Compressible
Linear, nonlinear
HD
Numerical simulations
Analytical,
experimental study,
numerical simulations
Cylindrical
Exploding/imploding shock
Nonlinear
HD
Planar
Planar shock
Compressible
Linear
Steady acceleration
Table 1.1: Summary of previous works that studied RMI
Fluid
HD
HD
42
Mikaelian [36]
Yang et al. [10]
Kind of study
(RTI) Analytical
study
(RTI) Analytical
study
Analytical study,
numerical simulations
Experimental
study
Numerical simulations
Analytical study
(RTI/ RMI) Analytical study
Analytical study
Numerical simulations
Analytical study
Experimental
study
Numerical simulations
Analytical study
Paper
Glimm et al. [44]
Zabusky
and
Zhang [52]
Samtaney [4]
Hosseini
and
Takayama [37]
Mikaelian [41]
Wheatley et al. [6]
Yu
and
Livescu [43]
Lombardini and
Pullin [40]
Samtaney [5]
Wheatley
al. [20]
Nishihara
al.[31]
et
Hohenberger
al.[23]
et
et
Geometry
Cylindrical,
spherical
Initializing the instability
Steady acceleration/ spherical
shock
Flow
Compressible
Equations
Nonlinear
Fluid
HD
Planar
Impulse model
Compressible
Nonlinear
HD
Planar
Planar shock
Compressible
Nonlinear
HD/MHD
Cylindrical
Cylindrical converging shock
Compressible
Nonlinear
HD
Cylindrical shells
Impulse acceleration
Linear, nonlinear
HD
Analytical study,
numerical simulations
Numerical simulations
Experimental
study
Analytical study
Planar
Impulse model
Incompressible,
irrotational,
viscous
Incompressible,
compressible
Linear
MHD
Planar
Planar shock
Compressible
Nonlinear
MHD
Planar
Impulse model
Linear
HD
Cylindrical,
spherical
Cylindrical
Chisnell-type shock
Incompressible,
miscible
Incompressible
Linear, nonlinear
HD
Linear
HD
Cylindrical,
spherical
Planar
Chisnell-type shock
Compressible, immiscible
Incompressible
Linear, nonlinear
HD
Planar shock
Compressible
Linear
HD/MHD
Planar
Planar shock
Compressible
Nonlinear
HD/MHD
Planar, cylindrical
Planar shock, cylindrical converging shock
Linear, nonlinear
Cylindrical,
spherical
Cylindrical converging shock
Compressible, irrotational, incompressible
Compressible
(RTI) Analytical
study
Numerical simulations
Numerical simulations
Numerical simulations
Analytical study,
numerical simulations
Experimental
study
43
Wheatley
et
al. [19]
Chapman and Jacobs [53]
Lombardini [39]
Kind of study
(RTI/
RMI)
Numerical simulations
Numerical simulations
Numerical simulations
Experimental
study
Analytical study
Steady acceleration
Table 1.2: Continue: Summary of previous works that studied RMI
Linear
HD
MHD
Paper
Sano et al.[22]
Sano et al.[54]
Pullin et al.[45]
Wheatley
al.[21]
et
Mostert et al.[49]
Mostert et al.[48]
Geometry
Planar
Initializing the instability
Planar shock
Flow
Compressible
Equations
Nonlinear
Fluid
MHD
Planar
Planar shock
Compressible
Nonlinear
MHD
Cylindrical
Converging cylindrical shock
Compressible
Nonlinear
MHD
Cylindrical,
spherical
Planar
MHD Converging shock
Compressible
Nonlinear
MHD
Planar shock
Compressible ,incompressible
Linear
MHD
Cylindrical,
spherical
Cylindrical
MHD Converging shock
Compressible
Nonlinear
MHD
MHD Converging shock
Compressible
Nonlinear
MHD
Table 1.3: Continue: Summary of previous works that studied RMI
44
Mostert et al.[47]
Kind of study
Numerical simulations
Numerical simulations
Analytical study ,
numerical simulations
Numerical simulations
Numerical simulation, analytical
study
Numerical simulations
Numerical simulation
45
1.3
Objectives and Contributions
As described earlier, several investigation of RMI have been conducted in planar geometry. Apart from the works of Mostert et al.([47, 49, 48]), in which they examined
nonlinear simulations of RMI in MHD, there are no other studies of RMI in the
literature. Particularly lacking are any linear analyses of RMI in MHD in converging geometry. In the actual context of ICF, the physical processes are exceedingly
complex: these include ablation physics, laser-plasma interaction, phase change of
materials, radiation coupled with hydrodynamics (an exceeding challenging field in
its own right) etc. In the spirit of scientific reductionism, one breaks down these physical processes into their “atoms” and investigates the details in term of simplified and
canonical problems. The whole history of hydrodyamics studies in gas shock-tube
experiments and in simulations stemmed from such a reductionist approach towards
understanding the basic hydrodynamic instabilities in ICF. Here we adopt a similar
philosophy and investigate RMI in converging cylindrical geometry. Furthermore,
linear stability analysis has a long tradition in fluid mechanics. Nonlinear processes
in RMI are complex, and it is natural to begin investigations of any fluid instabilities
with a linear approach. Furthermore, linear analyses is a computationally inexpensive approach, compared with computationally expensive nonlinear simulations, to
elucidate a number of physical processes.
The main topics of our investigations are:
• Linear simulations of RMI in the cylindrical geometry of two or three-dimensional
perturbations with normal (radial) magnetic field.
• Linear simulations of RMI in the cylindrical geometry of two-dimensional perturbations with an azimuthal magnetic field.
46
• Linear analysis of RMI in the framework of incompressible ideal MHD (analytical study).
• A comparison between linear and nonlinear simulations of MHD RMI for cylindrical geometry.
The main mathematical contributions of the thesis are:
• The nonlinear compressible ideal MHD equations are linearized about a time dependent base state, which will be discussed later, so we drive a set of hyperbolic
equations with source terms for the perturbed quantities.
• A numerical approach to investigate linear stability analysis of RMI in MHD in
cylindrical geometry. This includes an approach to deal with singular magnetic
fields at the origin.
• Developing analytical approaches in incompressible models of RMI in MHD. The
analytically derived models are then solved using a numerical inverse Laplace
transform method: this is the first application of these techniques to investigate
linear stability.
• An insight into the physics of RMI (and associated RTI) in cylindrical geometry,
and the suppression of these instabilities by a magnetic field. Furthermore, we
provide a physical interpretation of the suppression of these instabilities.
1.3.1
Thesis Outline
The investigation of MHD RMI in this thesis is performed in the cylindrical geometry.
Mainly, it includes two parts. First, the numerical simulations for compressible MHD
RMI. Second, the analytical analysis for an incompressible model of MHD RMI. The
thesis is divided into five chapters.
47
Chapter 2: The presence of hydrodynamic instabilities in inertial confinement fusion
and suppression by means of a magnetic field motivates us to investigate the RMI via
linear MHD simulations in cylindrical geometry. The physical setup of the ideal MHD
equations are introduced in the presence of a normal magnetic field. Linearizing the
ideal MHD equations about a time dependent base state we derive a set of hyperbolic
equations with source terms for the perturbed quantities. These equations are solved
using an upwind technique (the details of the numerical method are relegated to
Appendix A). The numerical results for both hydrodynamics and MHD for a positive
Atwood ratio interface, with purely azimuthal or axial, and combined azimuthalaxial perturbations, accelerated by a Chisnell-type shock[42] are investigated. In
particular, we observe in converging geometry, the Rayleigh-Taylor instability due to
the continued acceleration of the interface towards the center of convergence. The
presence of magnetic field suppresses the instabilities at the interface.
Chapter 3: The linear stability of both positive and negative Atwood ratio interfaces accelerated either by a fast magnetosonic or hydrodynamic shock in cylindrical
geometry. Since, the physical setup and the ideal MHD equations in the presence
of an azimuthal magnetic field are introduced. The linearized ideal MHD equations
about a time-dependent base state are solved using an upwind technique. An investigation of numerical results for both hydrodynamics and MHD for a positive or
negative Atwood ratio interface, with purely azimuthal perturbations accelerated by
a converging MHD shock is presented. The presence of magnetic field suppresses the
instabilities at the interface. The suppression mechanism is associated with the interference of two waves running parallel and anti-parallel to the interface that transport
vorticity and cause the growth rate to oscillate in time.
Chapter 4: A linearized model problem is studied here to understand the behavior of the interface seen in Chapters 2 and 3. In the framework of ideal incompressible
48
MHD equations, the stability of an impulsively accelerated, sinusoidally perturbed
density interface in the presence of the magnetic field. Here, we extend the earlier
linear stability analysis by Wheatley et al. [6, 21] to cylindrical geometry. The magnetic field first is initialized parallel to the acceleration direction. The second case
where the magnetic field is unperturbed and aligned with position of the interface.
The appropriate linearized initial value problem is solved by Laplace transforms. The
key features of the resulting incompressible linear models for the MHD RMI are then
discussed. We present a comparison between the incompressible models and the linear
compressible simulations that investigated in Chapters 2 and 3.
Chapter 5: The major findings to arise from this thesis and future work are
summarized in Chapter 5.
49
Chapter 2
Linear Simulations of the Cylindrical Richtmyer-Meshkov
Instability in the Presence of Radial Magnetic Field
2.1
Introduction
In the previous chapter, we introduced the linear analysis of RM flows. It is evident
from an examination of the RMI literature that there is a gap in our knowledge
regarding the suppression of RMI in converging geometry in the presence of a magnetic
field. Investigations of the MHD RMI thus far have primarily addressed the case where
the incident shock propagates from a light gas to a heavy one, corresponding to a
positive Atwood number (A > 0). Consequently, the effects of cylindrical convergence
on the MHD RMI should be first investigated for the A > 0 case. Note, however,
that the planar MHD RMI is also suppressed for A < 0, as clearly illustrated in Fig.
1 of Wheatley et al.[20]. This is important even in the A > 0 case because, although
not studied here, the reshock of the interface after shock reflection from the axis of
convergence will correspond to a diverging MHD RMI with A < 0. It is therefore
anticipated that the reshock induced perturbation growth will also be attenuated
by the presence of a magnetic field. In the original analysis by Richtmyer an initial
boundary value problem was derived, where the boundaries were formed by a reflected
and a transmitted shock emerging from the zeroth order interaction of an incident
shock with a contact discontinuity. The base state in this analysis is a self-similar
solution involving a reflected shock, a transmitted shock, and a contact discontinuity
50
that translates at an average constant speed due to the impulse by the incident shock.
The linearized compressible inviscid equations were solved by Richtmyer, in which the
perturbation amplitude growth rate undergoes a short transient before settling into
an oscillatory state around a constant asymptotic value. Clearly, there are a wealth
of configurations and parameters for which the cylindrically converging MHD RMI
could be investigated. A logical first step in such an investigation is linear stability
analysis of an A > 0 interface with azimuthal and/or axial perturbations.
Motivated by the presence of hydrodynamic instabilities in inertial confinement
fusion and suppression by means of a magnetic field, we investigate the RMI via
linear MHD simulations in cylindrical geometry. The physical setup is that of a
Chisnell-type converging shock interacting with a density interface with either axial
or azimuthal (2D) perturbations. The linear stability is examined in the context of
an initial value problem (with a time-varying base state) wherein the linearized ideal
MHD equations are solved with an upwind numerical method. Linear simulations
in the absence of a magnetic field, indicate that RMI growth rate during the early
time period is similar to that observed in Cartesian geometry. However, this RMI
phase is short-lived and followed by a Rayleigh-Taylor instability (RTI) phase with
an accompanied exponential increase in the perturbation amplitude. We examine
several strengths of the magnetic field (characterized by β =
2p
)
Br2
and observe a
significant suppression of the instability for β ≤ 4. The suppression of the instability
is attributed to the transport of vorticity away from the interface by Alfvén fronts.
51
2.2
2.2.1
Physical Setup and Governing Equations
Physical Setup
A schematic of the physical setup is shown in Fig. 2.1, depicting (a) A (r, θ)−section
with a fixed axial coordinate z and (b) A (r, z)−section with constant azimuthal
orientation θ. A cylindrical annulus domain is confined by left and right boundaries at
r = Rl and r = Rr , respectively. A density interface of azimuthal or axial perturbation
with wavenumber m or k is located at r = R0 that separates two fluids of densities
ρ1 and ρ2 . A Chisnell-type converging shock of Mach number Ms travels inwards
(further discussion about the Chisnell-type shock in Section 2.3.1), which is initialized
at r = Rs . Samtaney [4], Wheatley et al.[6, 20] studied the suppression of the RMI
for the case of an initial magnetic field normal to the plane of the driving incident
shock in planar geometry. It was recently made evident by the work of Mostert et
al.[47] that azimuthal, and spherical symmetry may be violated in cylindrical and
spherical geometry, respectively, in the presence of a magnetic field. Our intention
is to examine the effect of a magnetic field on the RMI in cylindrical geometry in
which the base state is only a function of the radial coordinate and employing any of
the configurations in Mostert et al.would make the base state in our linear stability
analysis a multi-dimensional function. To retain simplicity of the base state, our
preference is to examine a magnetic field that is initially oriented along the direction of
the shock motion so that it results in a purely hydrodynamic incident shock impinging
on the density interface. This configuration also implies a somewhat unrealistic radial
magnetic field because the solenoidal constraint on the magnetic field precludes a
converging purely radial field. Nevertheless, we further justify this choice by trimming
the domain of interest and cutting out the origin. At location Rl , the field lines exit
the domain radially and do not violate the solenoidal condition outside the domain. A
52
purely radial field inside the domain of interest may be conceivably created by current
carrying coils. We do not elaborate on such possibilities because this is beyond the
scope of the present study. An additional justification for the choice of a radial field
may be made by invoking the benefits of a reductionist approach: our approach is
a natural extension of previous studies in planar geometry by including cylindrical
geometry effects in a somewhat straightforward manner.
The magnetic field is characterized by its nondimensional strength using the
plasma beta β = 2p/Br2 , evaluated at R0 , where p is the initial pressure ahead of
the shock and Br is the radial component of magnetic field strength.The adopted
boundary conditions ar the with zero gradient of the velocity. Five parameters characterize the problem: the Mach number M of the incident Chisnell shock, the Atwood
number A = (ρ2 − ρ1 )/(ρ2 + ρ1 ), the azimuthal or axial perturbation wavenumber m or k, the plasma β at the interface and the distance ratio Rs /R0 . Our
numerical results are completely determined by these five parameters (Ms , A, m or
k, β, Rs /R0 ). To limit the scope of our investigation, our simulation domain is chosen
to be Rl = 10−5 , Rr = 1.9, and we fix Ms = 2.0, A = 0.667 as the ratio for an air-SF6
interface and Rs /R0 = 1.2/1.0 and mainly study the cases of different perturbation
wavenumbers m or k, and plasma β.
2.2.2
Governing Equations
Consistent with other studies we expect the time-scale for diffusion to be much longer
than the convective time scale in RMI [4, 19, 6, 20, 21]. Hence we neglect viscous,
thermal and plasma resistivity effects. Hence, the governing equations are those of
ideal MHD written below in conservation form in cylindrical coordinates:
∂W 1 ∂(rF(W)) 1 ∂G(W) ∂H(W)
+
+
+
= S,
∂t
r
∂r
r ∂θ
∂z
(2.1)
53
(a) (r, θ)-section
(b) (r, z)-section
Figure 2.1: Schematic diagram of the physical setup: (a) (r, θ)-section with constant
z and (b) (r, z)-section with constant θ. Two fluids of densities ρ1 and ρ2 are separated
by a perturbed interface located at r = R0 . A Chisnell-type incident shock is initially
located at r = Rs , with Mach number M .
54
where W ≡ W(r, θ, z, t) = {ρ, ρur , ρuθ , ρuz , Br , Bθ , Bz , e}T is the solution vector, ρ
is the density, (ρur , ρuθ , ρuz ) is the momentum, (Br , Bθ , Bz ) is the magnetic field,
e is the total energy per unit volume. These equations are already expressed in
non-dimensional form by choosing a reference density, pressure, and length scale (for
example, see Wheatley et al. [21] for the non-dimensionalization). F(W), G(W) and
H(W) are the fluxes of mass, momentum, magnetic flux and energy in the r, θ and
z directions given as
F=
ρur
ρu2r + p̃ − Br2
ρur uθ − Br Bθ
ρur uz − Br Bz
0
ur B θ − uθ B r
ur B z − uz B r
(e + p̃)ur − (B · u)Br
H=
G =
ρuz
ρuθ
ρur uθ − Br Bθ
ρu2θ + p̃ − Bθ2
ρuθ uz − Bθ Bz
uθ B r − ur B θ
0
uθ B z − uz B θ
(e + p̃)uθ − (B · u)Bθ
ρur uz − Br Bz
ρuθ uz − Bθ Bz
ρu2z + p̃ − Bz2
uz B r − ur B z
uz B θ − uθ B z
0
(e + p̃)uz − (B · u)Bz
,
(2.2)
where p̃ = p + 12 (B · B), is the sum of gas and the magnetic pressures. S is a
55
source term arising in cylindrical coordinates, given by
1
S=−
r
0, Bθ2
−
ρu2θ
− p̃, ρur uθ − Br Bθ , 0, 0, uθ Br − ur Bθ , 0, 0
T
To close the system we use the perfect gas equation of state so that e =
.
p
γ−1
(2.3)
+ 21 (B ·
B) + 12 ρ(u · u), where γ is the ratio of specific heats fixed at 5/3. Finally, we have
the solenoidal constraint on the magnetic field, i.e., ∇ · B = 0.
The initial magnetic field has an inverse radial dependence. Close to the origin,
large numerical errors are encountered owing to a lack of analytical cancellation of a
term proportional to r−3 in the radial momentum equation. These errors are eliminated by rewriting the radial component of the magnetic field as Br = Br∗ + Br1 where
β
is a curl-free portion of the magnetic field. The ideal MHD equations are
Br∗ =
r
rewritten as
∂ W̃ 1 ∂(rF̃(W̃)) 1 ∂ G̃(W̃) ∂ H̃(W̃)
+
+
+
= S̃,
∂t
r
∂r
r ∂θ
∂z
(2.4)
where W̃ = {ρ, ρur , ρuθ , ρuz , Br1 , Bθ , Bz , ẽ}T . Also, the modified fluxes and source
terms are
F̃ = {ρur , ρu2r + p̃t − Br Br1 , ρur uθ − Br Bθ , ρur uz − Br Bz , 0, ur Bθ − uθ Br , ur Bz (2.5)
− uz Br , (ẽ + p̃t + Br∗ Br1 )ur − (B.u)Br }T
G̃ = {ρuθ , ρur uθ − Br Bθ , ρu2θ + p̃t − Bθ2 , ρuθ uz − Bθ Bz , uθ Br − ur Bθ , 0, uθ Bz (2.6)
− uz Bθ , (ẽ + p̃t + Br∗ Br1 )uθ − (B.u)Bθ }T
H̃ = {ρuz , ρur uz − Br B, ρuθ uz − Bθ Bz , ρu2z + p̃t − Bz2 , uz Br − ur Bz , uz Bθ (2.7)
− uθ Bz , 0, (ẽ + p̃t + Br∗ Br1 )uz − (B.u)Bz }T ,
1
S̃ = − {0, Bθ2 − ρu2θ − p̃t + Br∗ Br1 , ρur uθ − Br Bθ , 0, 0, uθ Br − ur Bθ , 0, , 0}T ,
r
(2.8)
56
p
where ẽ =
+ 1 ((Br1 )2 +Bθ2 +Bz2 )+ 12 ρ(u2r +u2θ +u2z ) and p̃t = p+ 21 ((Br1 )2 +Bθ2 +Bz2 ).
γ−1 2
2.2.3
Linear Stability Analysis
If a small parameter is present, then the governing equations can be linearized using
an expansion in terms of the small parameter by retaining terms up to first order and
neglecting higher order terms. For the physical situation of a shock interacting with
a perturbed density interface, the ratio of the perturbation amplitude to wavelength
may be considered a small parameter. In general, linear stability can be broadly
classified into normal mode analysis or an initial value problem (IVP) approach.
The normal mode approach is suitable when the instability grows exponentially in
time (as in the Rayleigh-Taylor instability case). The IVP approach is the chosen
method when we have algebraic growth (as in the RMI case) or when transient growth
is investigated. Our approach detailed below is somewhat similar in spirit to the
original analysis of Richtmyer [2] and Yang et al. [10]. Richtmyer derived a set of
linear equations for the case of a shock interacting with a density interface wherein
the reflected and transmitted waves after the shock refraction were both shocks.
Yang et al. on the other hand, derived a set of linear equations for the case where
the reflected wave is a rarefaction. In both analyses, the amplitude of the interface
itself is not present; and in both cases the growth rate (or rather the time derivative
of the interface amplitude) is computed numerically. In this sense, both these were
instances of linear stability analysis as IVPs. In our investigation, an analytical
process similar to Richtmyer and Yang et al. is not tractable owing to the complexity
of the refractions in MHD, and the converging geometry. Instead, we follow the
approach of Samtaney [5] who proposed a numerical method (also an IVP approach)
for solving the linearized equations resulting in stability analysis of shocked interfaces
in gas dynamics and magnetohydrodynamics.
57
To examine the linear stability of a converging shock interacting with a perturbed interface, we linearize the ideal MHD equations by splitting the solution
into a base solution with {˚} notation, and perturbed solution, with {ˆ} notation as
W̃(r, θ, z, t) = W̊(r, t) + ǭŴ(r, t)ei(mθ+kz) , where m and k denote the azimuthal and
√
axial perturbation wavenumbers and i ≡ −1. The time dependent base state is a
function of the radial coordinate, shocks and a contact discontinuity. Details of the
base state are presented later. Substituting the base and perturbed solutions into the
governing equations, we get the partial differential equations governing the base and
perturbed states as follows:
∂ W̊ 1 ∂
+
(rF(W̊)) = S̊,
∂t
r ∂r
(2.9a)
1
∂ Ŵ 1 ∂
+
(rM(W̊)Ŵ) + imN(W̊)Ŵ + ikQ(W̊)Ŵ = Ŝ,
∂t
r ∂r
r
(2.9b)
58
0
1 2
1
2
∗
˚
− [B̊θ − ρ̊ůθ − p̃t + Br B̊r ]
r
1
−
[ρ̊ů
ů
−
B̊
B̊
]
r θ
r θ
r
0
,
S̊ =
0
1
[ů
B̊
−
ů
B̊
]
−
θ
r
r
θ
r
0
0
0
− 1 [2B̊θ B̂θ − 2ρ̊ůθ ûθ − ρ̂ů2θ − p̃ˆt + Br∗ B̂r1 ]
r
− 1 [ρ̊ů û + ρ̂ů ů + ρ̊û ů − B̊ B̂ − B̂ B̊ ]
r θ
r θ
r θ
r θ
r r θ
0
Ŝ =
0
− 1r [ûθ B̊r + ůθ B̂r − ûr B̊θ − ůr B̂θ ]
0
0
(2.10)
,
◦
where M(W̊), N(W̊) and Q(W) are respectively the Jacobians of the fluxes F̃, G̃
and H̃ with respect to the base state. We note here that the above equations, derived
from ideal MHD equations, are scale-free and as such one may choose the radial
domain extent arbitrarily.
Equations (2.9a) and (2.9b), governing the base and perturbed state, respectively,
are both hyperbolic partial differential equations. These are numerically solved with
an explicit time marching upwind procedure. The details of the numerical method
are given in Appendix A.
59
2.3
Linear Stability Results
In this section, we present results from linear stability simulations. We begin by
examining the time dependent radially varying base state solution, followed by the
perturbed quantities for hydrodynamics and MHD: both of these are presented first
for purely azimuthal perturbations (m 6= 0, k = 0), then for purely axial perturbations
(m = 0, k 6= 0) followed by a combination of axial and azimuthal perturbations. In
particular, we focus on vorticity in the domain as we consider the baroclinic vorticity
to be the dominant mechanism for the instability of the interface. We also present
comparisons with nonlinear simulations. In all cases, except when specified otherwise,
the radial domain is discretized with 1600 cells. As briefly shown for a specific case
in Section 2.3.5, differences in growth rates between employing 1600 or 8000 cells are
small.
2.3.1
Base State Solution
The time-dependent base state equations are independent of perturbed state variables;
the azimuthal, m, and axial perturbation, k wavenumbers. Besides, choosing B̊r to
be inversely proportional to radius, the induction equation is satisfied identically and
decoupled from the other equations.
The base state is examined by plotting the density profiles at time t = 0.0 (initial
condition) and t = 0.552 in Fig. 2.3. At t = 0, the r− dependent flow field behind the
initialized circularly converging gas-dynamic shock is computed based on the approximate analysis by Chisnell [42] and is dubbed the ”Chisnell-type” shock. Chisnell’s
analysis provides the radial variation of the density, radial velocity and pressure behind the shock as a function of radius. It turns out that at t = 0 the density profile is
a mildly increasing function of radius just behind the shock (see Fig 2.3), whereas the
60
radial velocity and pressure decrease in magnitude behind the shock front as shown
in Fig. 2.2 (a,b) respectively. Specifying a constant state behind the shock, while
convenient, leads to disturbances which are best avoided. An alternative to this socalled Chisnell-type shock would be to employ a circular Riemann problem to drive
a shock onto the interface. Although this latter method is also convenient, it results
in extra waves, i.e., a contact and rarefaction wave, in addition to a shock [47, 49]
The reflected waves from the shock-density interface may subsequently interact with
the contact and can be re-reflected and affect the results of the linear analysis. Hence
the circular Riemann problem is not our choice to induce the incident shock. The
Chisnell-type incedint shock (IS) is perhaps the cleanest technique for generating a
single radially moving shock. IS moves radially towards the contact discontinuity cd
and refracts into a transmitted shock (TS) and a reflected shock (RS).
The density spacetime diagram illustrates the wave diagram of the base state more
clearly. The density interface (cd) is a vertical line in Fig. 2.4 indicating that it is
stationary until it interacts with the incident Chisnell-type shock (IS). Then, the IS
bifurcates into a transmitted shock (TS) and a reflected shock (RS). The IS provides
an impulse to the CD accelerating it which is seen between the two shocks RS and TS.
Although the CD appears to be a straight line in the spacetime diagram, the second
time derivative is in fact non-zero. We note that at t = 0.552, the transmitted shock is
about to leave the domain and therefore, all the results are generated until this time.
In Fig. 2.5 the base velocity ůr of the CD is plotted, where it shows that as the CD
approaches the origin it decelerates. As discussed later, this deceleration (or negative
acceleration) is responsible for a Rayleigh-Taylor (RT) phase of the instability of
the interface. Hence, the interface is subject to an impulsive acceleration by the IS
(resulting in the RMI) and a continuous acceleration phase due to geometric effects
(resulting in the RTI). Fig. 2.6(a) shows the impulsive acceleration of the interface
ůr
61
0
t=0
−0.2
−0.4
−0.6
−0.8
−1
−1.2
−1.4
−1.6
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
r
p̊
(a) ůr
5
4.5
4
3.5
3
2.5
2
1.5
1
t=0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
r
(b) p̊
Figure 2.2: Profiles of base state variables behind the shock front at t = 0 (a) radial
velocity ůr (b) pressure p̊.
62
20
t = 0.0
t = 0.552
18
16
14
◦
ρ
12
CD
10
TS
8
6
RS
RS
4
2
CD t=0
IS t=0
0
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
r
Figure 2.3: Density profiles of base state at t = 0.0 (initial conditions) and t = 0.552.
IS denotes the incident shock, TS denotes the transmitted shock, RS denotes the
reflected shock and CD is the contact discontinuity.
63
Figure 2.4: Spacetime diagram of density field of the base state. IS denotes the
incident shock, TS denotes the transmitted shock, RS denotes the reflected shock
and CD is the contact discontinuity.
by the interaction with the incident shock; while the magnified portion in Fig. 2.6(b)
shows the acceleration of the interface as it moves towards the origin. We demarcate
the RMI from the RTI phase at t ≈ 0.21 in our simulations.
2.3.2
Evolution of the Perturbations
In this section, the results for hydrodynamics (HD) (i.e., β = ∞), and MHD (i.e.,
finite β) are presented.
For each case we discuss the case of purely azimuthal
(m 6= 0, k = 0 in equation 2.9b) and purely axial perturbations (m = 0, k 6= 0 in
equation. 2.9b). Following Lombardini and Pullin [40], we non-dimensionalize the
time scale using 1/(a0 K) as our reference time scale, where a0 is the unshocked
sound speed ahead of the incident shock, K = m/R0 (k) for purely azimuthal (axial)
modes. The perturbation growth rate ḣ is normalized by the theoretical growth rate
64
-0.84
ůr
-0.86
-0.88
ůr
-0.9
-0.92
RTI phase
-0.94
RMI phase
-0.96
-0.98
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
0.55
t
Figure 2.5: Time history of radial velocity ůr at CD.
Figure 2.6: Time history of (a) acceleration of CD, (b) acceleration of CD with
approximate demarcation of RMI and RTI phases.
65
ḣ∞ expressed as
ḣ∞ =
+
F Ms h +
0 ∆w(1 + KR0 A )
,
R0
(2.11)
where ∆w is the impulsive change in the interface velocity, A+ is the post-shock
Atwood number, FMs is a parameter depending on Mach number of the incident
shock (FMs = 0.85 for our cases) [40]. h+
0 the initial perturbation amplitude, which
in linear stability cancels out, is chosen as unity for convenience .
Solution Details for Hydrodynamics
We begin our discussion with the hydrodynamics case followed by MHD results.
Purely Azimuthal Perturbations: We now examine the case of purely azimuthal perturbations, in which there is no variation of the solution in the z direction. The scaled growth rate history is plotted in Fig. 2.7 for different perturbation
wavenumber m = 16, 64, 128, 256. The growth rate for the largest wave number case
m = 256 in hydrodynamics is similar to the Cartesian geometry case until about
scaled time a0 tm/R0 ≈ 30. The growth rate increases from zero and then oscillates
around unity; this is the asymptotic growth rate one would attain in Cartesian geometry. However, in our cylindrical case, as the interface approaches the origin, the RTI
manifests itself and we observe an exponential increase in magnitude of the growth
rate. The onset of RTI and its distinction from the RMI phase is most evident for
the largest wavenumber (m = 256), while the distinct separation between RMI and
RTI for smaller wavenumber is not very clear and for the lowest perturbation wave
number (m = 16) the RTI lasts only a short time. In all cases, the exponential growth
during the RTI results in a phase reversal of the amplitude. The spatial variation of
the perturbed radial and azimuthal components of the velocity, perturbed pressure,
and z− component of the vorticity are plotted in Fig. 2.8 at t = 0.552. The radial
ḣ/ḣ1
66
2
1.5
1
0.5
0
−0.5
−1
−1.5
−2
0
RMI
20
40
m = 256
m = 128
m = 64
m = 16
RTI
60 80 100 120 140 160
a0tm/R0
Figure 2.7: Time history of scaled growth rate in hydrodynamics for azimuthal perturbation wavenumber m = 16, 64, 128, 256.
velocity component shows a large negative peak at the interface location (the average
position of the interface in the base state solution is superimposed on the perturbed
radial velocity). This magnitude of this peak corresponds to the growth rate of the
interface amplitude. During the refraction of the incident shock at the interface, vorticity is deposited because of the baroclinic source term, and is the essential driving
mechanism of the RMI, as discussed by Hawley and Zabusky [17] and Samtaney and
Zabusky [55]. The z− component of the vorticity ω̂z =
1 ∂(r ûθ )
r ∂r
− 1r imûr is a vortex
sheet (VS) for a sharp interface and visible as a large peak in Fig. 2.8(d). In this
case with zero magnetic field, the position of the VS coincides with the position of
the interface. This can also be seen as the large gradient in ûθ at the interface in
Fig. 2.8(b).
After the primary interaction of the shock at the interface, secondary waves propagate along the azimuthal direction. These secondary interactions are manifested in
67
our linear simulations as the oscillations of the fields between the interface and the
reflected/transmitted shocks. The wavelength of these oscillations decreases with increasing azimuthal wavenumber m. The secondary interactions also leave a signature
on the growth rate in the form of oscillations before the RTI overwhelms the RMI
growth rate (see Fig. 2.7).
Purely Axial Perturbations: We plot the scaled growth rate history for different perturbation wavenumber k = 16, 64, 128, 256 in Fig. 2.9. Similar to the case
of purely azimuthal perturbations, the growth rate for purely axial perturbations exhibits an initial rapid increase, after which the perturbation growth rate oscillates
around a constant value for a short duration followed by the exponential increase in
the magnitude of the growth rate because of the RTI. The onset of RTI is clear for
the largest wavenumber (k = 256), although for a smaller one (k = 16), the demarcation between RMI and RTI is not so apparent. For the highest axial wave number
(k = 256) the perturbed radial and azimuthal components of the velocity, perturbed
pressure, and θ− component of the vorticity ω̂θ = ikûr −
∂ ûz
∂r
are plotted in Fig. 2.10
at t = 0.552. The behavior of these fields is very similar to those observed for the
azimuthal case. The peak vorticity magnitude is coincident with the mean interface
position and the peak magnitude of the radial velocity (also on the interface) corresponds to the growth rate of the perturbation amplitude. The secondary waves are
manifested as oscillations between the interface and the reflected/transmitted shocks.
68
20
40
t = 0.552
CD
0
20
0
−20
10
ûθ
−40
ûr
t = 0.552
CD
30
−60
0
−10
−80
−20
−
−100
−
−30
−120
−40
−
−140
−
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
r
r
(a) ûr
(b) ûθ
8000
120
t = 0.552
CD
100
6000
80
VS
0
ω̂z
40
0
4000
60
p̂
t = 0.552
CD
2000
20
0
0
−
−
−20
−40
−60
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
−2000
−4000
−
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
r
r
(c) pb
(d) ω̂z
Figure 2.8: Profiles of perturbed variables in hydrodynamics at t = 0.552: (a) rcomponent of velocity, (b) θ-component of velocity, (c) pressure, (d) z-component of
vorticity. The azimuthal perturbation wavenumber is m = 256.
69
2
RMI
1.5
ḣ/ḣ1
1
k = 256
k = 128
k = 64
k = 16
RTI
0.5
0
−0.5
−1
−1.5
0
20
40
60
80 100 120 140 160
a0tk
Figure 2.9: Time history of scaled growth rate in hydrodynamics for axial perturbation wavenumber k = 16, 64, 128, 256.
Combined Azimuthal-Axial Perturbations: We now examine a combination
of three dimensional azimuthal and axial perturbations. Following [40], the time
scale is non-dimensionalized by 1/(a0 K), where the combination of two wavenumbers
p
K = m2 + k 2 R02 is used in 3-D modes. The growth rate is non-dimensionalized by
the same factor given in equation 2.11.
The time histories of normalized growth rate in the hydrodynamics case for
different combinations of axial/azimuthal wavenumbers (k, m) = (16, 16), (16, 256),
(256, 16), (256, 256) are plotted in Fig. 2.11. The overall behavior of growth rate is
similar to that seen for purely axial and azimuthal wave numbers. One interesting
feature is that after the onset of RTI, the same total axial/azimuthal combination
with the larger azimuthal contribution has a larger growth of the RTI. This is also
evident by comparing the purely axial (k=256) and azimuthal perturbations (m=256)
in Figures 2.7 and 2.10: the growth rate is larger during RTI for the azimuthal case
70
20
t = 0.552
CD
0
−20
ûz
ûr
−40
−60
−80
−100
−120
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
r
40
t = 0.552
CD
30
20
10
0
−10
−20
−30
−40
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
r
(a) ûr
40
(b) ûz
t = 0.552
CD
30
8000
20
6000
0
VS
4000
10
0
ω̂θ
p̂
t = 0.552
CD
−10
2000
0
−20
−
−30
−2000
−40
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
−4000
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
r
r
(c) pb
(d) ω̂θ
Figure 2.10: Profiles of perturbed variables in hydrodynamics at t = 0.552: (a) rcomponent of velocity, (b) z-component of velocity, (c) pressure, (d) θ-component of
vorticity. The axial perturbation wavenumber is k = 256.
71
ḣ/ḣ1
than the corresponding axial one with the same wave number.
2
1.5
1
0.5
0
−0.5
−1
−1.5
−2
0
RMI
k = m = 256
k = m = 128
k = m = 64
k = m = 16
RTI
50
100
200
150
p
2
2
a0t m + k /R0
250
Figure 2.11: Time history of scaled growth rate in HD for 3-D perturbation wave
numbers k = m = 16, 64, 128, 256.
Solution Details of MHD
Several values of β characterizing the strength of the magnetic field are examined.
Here we present the details of the solution for β = 4 for purely azimuthal and purely
axial perturbations. A discussion of varying β will be presented later in Section 2.3.5.
Purely Azimuthal Perturbations: In Fig. 2.12, we plot the time history of
scaled growth rate in MHD for different perturbation wavenumber m = 16, 64, 128, 256.
For the two largest wavenumber case (m = 256, 128), the growth rate increases initially and then reduces because of the magnetic field. During the time duration of the
linear simulation, the growth rate eventually reduces and oscillates around a small
value below zero. Recall that in Cartesian geometry, the growth rate at the interface
is suppressed by the magnetic field and it oscillates around the asymptotic value zero
72
(See Ref. [5] for examples). We conjecture that, at the late time stage for these cases,
vorticity is being continually generated at the accelerating interface, and MHD waves
remove it continuously from the interface. This probably results in a small residual
growth rate. For the smaller wavenumber cases (m = 64, 16), the time duration of
the simulation is insufficient for any large suppression to occur. In all these cases, the
onset of RTI is not clearly demarcated from the RMI phase.
1
m = 256
m = 128
m = 64
m = 16
0.8
ḣ/ḣ1
0.6
0.4
0.2
0
−0.2
−0.4
0
20 40 60 80 100 120 140 160
a0tm/R0
Figure 2.12: Time history of scaled growth rate in MHD for azimuthal perturbation
wavenumber m = 16, 64, 128, 256. The initial plasma beta is β = 4 at the interface
R0 = 1.0.
Because the suppression is largest for the highest wavenumber, as one may intuitively expect, we now examine the solution details and elucidate the suppression
mechanisms for m = 256. The perturbed r- and θ-components of velocity, r- and θcomponents of magnetic field, z-component of vorticity and z-component of current
density ĵz (defined as ĵz =
1 ∂(r B̂θ )
r ∂r
− 1r imB̂r ) are plotted in Fig. 2.13. Wheatley et
al. [6], in analysis of impulsively accelerated interfaces in the context of incompressible MHD, concluded that Alfvén fronts transport vorticity away from the interface
73
and constitute the main mechanism for the RMI suppression. A similar behavior is
observed in our linear simulations. Fig. 2.13(e) shows two dominant vortex sheets
which coincide with Alfvén fronts evident in the current density in Fig. 2.13(f). As in
the planar case, these Alfvén fronts transport vorticity away from the interface and
hence stabilize the interface.
Purely Axial Perturbations: The time history of scaled growth rate in MHD
is plotted in Fig. 2.14, for different perturbation wavenumber k = 16, 64, 128, 256. For
the two largest wavenumber case (k = 256, 128), the growth rate increases initially
followed by a reduction because of the magnetic field. During the time duration of
the investigation, the growth rate eventually reduces and oscillates around a small
value below zero. The behavior is very similar to that observed for the azimuthal
perturbations. For the smaller wavenumber cases (k = 64, 16), the time duration of
the simulation is insufficient for significant suppression to occur. In all these cases,
the onset of RTI is not clearly demarcated from the RMI phase.
We now examine the solution for k = 256 for which the suppression of the instability is the largest. The θ-component of vorticity and θ-component of current
density ĵθ = ik B̂r −
∂(rB̂z )
∂r
are plotted in Fig. 2.15. The two vortex sheets evident
in the vorticity plot, are coincident with the Alfvén fronts seen in current density.
The interface, devoid of the large initial baroclinic vorticity, loses the main driving
mechanism of the instability and perturbation growth is significantly suppressed.
Combined Azimuthal-Axial Perturbations: The time histories of normalized
growth rate in the MHD case for different combinations of axial/azimuthal wavenumbers (k, m) = (16, 16), (16, 64), (16, 128), (16, 256), (64, 16), (128, 16), (256, 16), (256, 256)
are plotted in Fig. 2.16. For these cases, we chose the plasma beta to be β = 2 at the
interface R0 = 1.0. While the instability is suppressed for these wavenumber combinations, the suppression is still somewhat less than that observed for β = 4 for purely
74
50
20
t = 0.552
CD
40
t = 0.552
CD
15
30
10
ûθ
ûr
20
10
5
0
0
−10
−5
−20
−10
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
r
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
r
(a) ûr
200
(b) ûθ
t = 0.552
CD
150
B̂θ
B̂r
100
50
0
−50
−100
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
r
4000
t = 0.552
CD
VS
3000
2000
VS
1000
0
−1000
−2000
−3000
−4000
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
r
(e) ω̂z
(d) B̂θ
Jˆz
ω̂z
(c) B̂r
40
t = 0.552
CD
30
20
10
0
−10
−20
−30
−40
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
r
8000
t = 0.552
CD
6000
AF
4000
AF
2000
0
−2000
−4000
−6000
−8000
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
r
(f) ĵz
Figure 2.13: Profiles of perturbed variables in MHD (a) r-component of velocity,
(b) θ-component of velocity, (c) r-component of magnetic field, (d) θ-component of
magnetic field, (e) z-component of vorticity, (f) z-component of current density. The
perturbation wavenumber is m = 256 and the initial plasma beta is β = 4 at the
interface R0 = 1.0.
75
1
0.8
ḣ/ḣ1
0.6
k = 256
k = 128
k = 64
k = 16
0.4
0.2
0
−0.2
−0.4
0
20 40 60 80 100 120 140 160
a0tk
Figure 2.14: Time history of scaled growth rate in MHD for axial perturbation
wavenumber k = 16, 64, 128, 256. The initial plasma beta is β = 4 at the interface
R0 = 1.0.
axial or azimuthal perturbations. Another noteworthy feature can be distinguished by
comparing the (k, m) = (16, 128), (128, 16) combinations. Both show similar overall
behavior, but the oscillations associated with transverse reverberations in the higher
azimuthal wavenumber persist while those in the higher axial wavenumber are nearly
absent.
Vorticity spacetime evolution: To further elaborate on the vorticity evolution
and distinguish between the MHD and hydrodynamics cases, we plot the vorticity
space-time diagram for perturbation wavenumber m = 256 in Fig. 2.17. In hydrodynamics, the vorticity is confined to the interface and undergoes a sign change in
keeping with the phase change associated with the RTI phase. Moreover, the vorticity
on the interface increases in magnitude as time progresses. In MHD, on the other
hand, the dominant vorticity is on Alfvén waves which straddle the interface in the
76
4000
VS
3000
2000
t = 0.552
CD
VS
ω̂θ
1000
0
−1000
−2000
−3000
−4000
0
0.2 0.4 0.6 0.8
1
r
1.2 1.4 1.6 1.8
2
(a) ω̂θ
8000
t = 0.552
CD
6000
AF
4000
Jˆθ
2000
AF
0
−2000
−4000
−6000
−8000
0
0.2 0.4 0.6 0.8
1
r
1.2 1.4 1.6 1.8
2
(b) ĵθ
Figure 2.15: Profiles of perturbed variables in MHD: (a) θ-component of vorticity,
(b) θ-component of current density. The perturbation wavenumber is k = 256 and
the initial plasma beta is β = 4 at the interface R0 = 1.0.
77
1
m = 256, k = 256
m = 256, k = 16
m = 16, k = 256
m = 128, k = 16
m = 16, k = 128
m = 64, k = 16
m = 16, k = 64
m = 16, k = 16
0.8
ḣ/ḣ1
0.6
0.4
0.2
0
−0.2
−0.4
0
50
100
200
150
p
2
2
a0t m + k /R0
250
Figure 2.16: Time history of scaled growth rate for 3-D perturbations MHD for
different wavenumbers k = m = 256, m = 16 with k = 64, 128, 256, k = 16 with
m = 64, 128, 256, and k = m = 16 for β = 2.
spacetime diagram. There is also some evidence of weak secondary vorticity in the
spacetime diagram but this is of little consequence for the growth of the instability.
2.3.3
Effect of Converging Geometry
We investigate the effect of converging geometry on the linear simulations of hydrodynamics case and MHD case for plasma β = 4. The wavenumber for 2-D perturbations
in planar geometry is k = 256, and in cylindrical geometry for purely azimuthal (axial) modes are m(k) = 256. We plot the growth rate using 12800 cells, since the planar
case needs high resolution, and it is scaled by ḣ∞ . The time is scaled by (1/a0 K)
where K = k in the planar geometry and for the cylindrical geometry it is presented
in section 2.3.2. For the hydrodynamics instabilities, Fig. 2.18(a) shows an agreement
between the planar and the cylindrical geometries at the early time, t < 20. Then,
in the planar case the growth rate oscillates around an asymptotic value less than
78
one. Whereas, in the cylindrical geometry the RTI starts to play role as presented in
Section 2.3.2. The scaled growth rate of linear simulations of MHD-RMI is plotted in
Fig. 2.18(b). It also shows similarities in the early time since the initial growth rate is
unaffected by the magnetic field. Then as the suppression of RMI occurs, the growth
rate in the planar geometry oscillates around zero while in the cylindrical geometry
they oscillate around values less than zero.
2.3.4
Effect of Magnetic Field Strength
We now examine the effect of initial magnetic field strength characterized by the nondimensional parameter β. The time history of scaled growth rate in MHD for plasma
beta β = ∞ (HD), 256, 64, 4, 2 is shown in Fig. 2.19, 2-D perturbations of azimuthal
wavenumber m = 256, Fig. 2.19(a), and axial wavenumber k = 256, Fig. 2.19(b).
The 3-D perturbations case we chose β = ∞ (HD), 16, 4, 2, 1 and the perturbations
wavenumbers are k = m = 256, Fig. 2.19(c). As expected the smallest β ≤ 4, (largest
magnetic field) has the strongest suppression of the instability for both the azimuthal
and axial perturbation cases. Interestingly, the weakest magnetic field case (β = 256)
halts the exponential increase in the magnitude of the growth rate during the RTI
phase. We observe a second frequency related to the speed of the Alfvén fronts. The
1
Alfvén speed is proportional to β − 2 . Comparing the time histories in Figure 2.19 for
β = 256 to β = 64, we observe the first local minimum shift from t ≈ 128 to t ≈ 64
representing a period halving, and comparing β = 64 to β = 4 we observe roughly
a decrease by a factor of four in the period of the oscillations. We further point out
that even though growth rates for the azimuthal and axial cases seem remarkably
similar, there are minor differences.
79
2.3.5
Comparison Between Linear and Nonlinear Simulations
Nonlinear simulation for RMI are performed using the ideal MHD equations expressed
in strong conservation form. The numerical method is a predictor-corrector unsplit
method in cylindrical coordinates, and is similar to that developed by Samtaney[4],
and Samtaney et al.[56]. A seven-wave linearized Riemann solver method is used to
compute the fluxes; and the solenoidal condition on the magnetic field is maintained
by a projection technique. The details of the numerical method are presented in
Appendix B. The nonlinear simulation is conducted on a grid extending azimuthally
from θ = 0 to 2π/128 with periodic boundary conditions applied at these bounds. In
the radial direction, the grid extended from r = 0.3 to 1.9. The boundary conditions
applied at these bounds are identical to those applied in the linearized simulations.
The Chisnell-type shock is initialized at Rs = 1.2 and the initial interface location is
R0 = 1, as for the linearized simulations. The computational grid has 3000 (8000) cells
in the radial direction for m = 16 (m = 128) and 128 cells in the azimuthal direction.
Examining a snapshot of density and vorticity for wavenumber m = 16 and magnetic
field β = 4 in the nonlinear simulation, shown in Fig. 2.20(a), shows clear evidence
of vorticity on MHD shock fronts and none on the interface. These observations are
similar those in Wheatley et al.[20] for the Cartesian case. For this case, a comparison
of the perturbation growth rate between linear and nonlinear simulations is shown
in Fig. 2.20(b). Further comparisons of the nonlinear growth rate of a single mode
interface in the linear and nonlinear simulations are shown for wavenumber m = 16
and a lower value of the magnetic field β = 16 in Fig. 2.20(c), and for a higher
wavenumber m = 128 in Fig. 2.20(d). The linear simulations performed here had
the same number of radial cells as the nonlinear simulations. In Fig. 2.20(d), we also
80
plot the growth rate for linear simulations using 1600 cells. The difference in linear
simulations at 1600 and 8000 cells is relatively small and one may consider this as a
rudimentary convergence test (further discussion on convergence is in Appendix A).
The initial amplitude of the interface in the nonlinear simulation is chosen to be
5% of the wavelength, i.e.,
1 2π
.
20 m
In the nonlinear case, the perturbation amplitude is
compressed by the incident shock. In linear theory, the amplitude of perturbation is
arbitrary and the compression by the incident shock is not part of the linear theory.
The growth rate comparison is performed after the incident shock passage over the
perturbed interface. The nonlinear and linear growth rates show satisfactory agreement at early times. As nonlinear effects become important, the agreement worsens
at late times for the low wavenumber cases (m = 16) particularly for the weaker
magnetic field case (β = 16).
81
(a) Hydrodynamics
(b) MHD
Figure 2.17: (a) Spacetime diagram of vorticity ω̂z for hydrodynamics. (b) Spacetime
diagram of vorticity ω̂z for MHD. The perturbation wavenumber is m = 256 and the
plasma beta is β = 4 for MHD.
ḣ/ḣ1
82
1.5
1
0.5
0
−0.5
−1
−1.5
−2
Planar k = 256
−2.5
Cyl. axial k = 256
Cyl. azimuthal m = 256
−3
0 20 40 60 80 100 120 140 160 180
Scaled time
ḣ/ḣ1
(a) Hydrodynamics.
0.6
0.5
0.4
0.3
0.2
0.1
0
−0.1
−0.2
0
Planar k = 256
Cyl. axial k = 256
Cyl. azimuthal m = 256
20 40 60 80 100 120 140 160 180
Scaled time
(b) MHD.
Figure 2.18: Time history of scaled growth rate of planar and cylindrical geometries
for k = 256 and m = 256 (a) hydrodynamics case and (b) MHD case for beta β = 4.
83
1.5
β=1
β = 256.0
β = 64.0
RTI β
= 16.0
β = 4.0’
RMI
ḣ/ḣ1
1
0.5
0
−0.5
−1
0
20 40 60 80 100 120 140 160
a0tm/R0
(a) Purely azimuthal perturbations m = 256
1.5
β=1
β = 256.0
RTI β = 64.0
β = 16.0
β = 4.0
RMI
1
ḣ/ḣ1
0.5
0
−0.5
−1
−1.5
0
20 40 60 80 100 120 140 160
a0tk
(b) Purely axial perturbations k = 256
1.5
RMI
1
ḣ/ḣ1
0.5
β=1
β = 16.0
β = 4.0
β = 2.0
RTI β = 1.0
0
−0.5
−1
−1.5
0
50
100
200
150
p
a0t m2 + k2/R0
250
(c) 3-D perturbations m = 256, k = 256
Figure 2.19: Time history of scaled growth rate in MHD for plasma beta β =
∞(HD), 256, 64, 4, 2 at the interface for (a) m = 256, k = 0, (b) m = 0, k = 256 and
(c) β = ∞(HD), 16, 4, 2, 1 for m = 256, k = 256.
84
1
Nonlinear
Linear
0.8
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
0
1
2
3
4
5
6
7
8
9
10
11
(a) Nonlinear
16, β = 4
simulation,
=
m
(b) m = 16, β = 4
1
Nonlinear
Linear, 8000 cells
Linear, 1600 cells
1.5
Nonlinear
Linear
0.5
1
0.5
0
0
-0.5
-0.5
0
1
2
3
4
5
6
(c) m = 16, β = 16
7
8
9
10
0
10
20
30
40
50
60
70
80
(d) m = 128, β = 16
Figure 2.20: (a) Density and vorticity snapshot in nonlinear simulation. (b-d)
Comparison of nonlinear growth rate between linear and nonlinear simulations.
85
Chapter 3
Linear Analysis of Converging Richtmyer-Meshkov
Instability in the Presence of an Azimuthal Magnetic Field
3.1
Introduction
The effect of transverse magnetic field, oriented parallel to the interface, recently
has been investigated analytically and numerically [32, 24, 33, 21]. It is found that
the existence of a transverse magnetic field in a super-conductive plasma, behaves
in a fashion somewhat similar to surface tension on the interface. In this chapter,
we investigate the linear stability of both positive and negative Atwood ratio interfaces accelerated either by a fast magnetosonic or hydrodynamic shock in cylindrical
geometry. For the (MHD) case, we examine the role of an initial seed azimuthal
magnetic field on the growth rate of the perturbation. In Section 3.2, the physical
setup is formulated and the ideal MHD equations are summarized. Here, we extend
the earlier numerical linear stability analysis by Samtaney [5] to MHD in cylindrical geometry. In Section 3.3, we present numerical results for both hydrodynamics
and MHD for positive/negative Atwood ratio interface (A > 0, A < 0, respectively),
with azimuthal perturbations, accelerated by an MHD converging shock. Finally, a
comparison between linear and nonlinear simulations is discussed in Section 3.4.
86
3.2
3.2.1
Physical Setup and Governing Equations
Physical Setup
The schematic of the physical setup is shown in Fig. 3.1 for (r, θ)−section with fixed
axial coordinate z. A cylindrical annulus domain is confined by left and right boundaries at r = Rl and r = Rr , respectively. A density interface of azimuthal perturbation
wavenumber m is located at r = R0 that separates two conducting fluids of different
densities ρ1 and ρ2 . A converging cylindrical fast magnetosonic MHD shock with
MHD Mach number Ms (based on the fast magnetosonic speed) is located at r = Rs ,
where the setup of the initial condition is given in Section 3.2.2. The spatially varying
azimuthal magnetic field is characterized by its non-dimensional strength using the
parameter beta β = 2p/(Bθ,0 )2 , where p is the initial pressure ahead of the shock and
Bθ,0 is the azimuthal component of the magnetic field at r = R0 . Five parameters
characterize the problem: the Mach number Ms of MHD shock, the Atwood number
A = (ρ2 − ρ1 )/(ρ2 + ρ1 ), the azimuthal perturbation wavenumber m, the plasma β
and the distance ratio Rs /R0 . These five parameters (Ms , A, m, β, Rs /R0 ) essentially
span the space of numerical results. To limit the scope of our investigation, our simulation domain is chosen to be Rl = 0.1, Rr = 2.6, and we fix Ms = 2.0, A = ±0.667,
the positive Atwood ratio corresponding to air-SF6 interface (vice verse for negative
Atwood ratio) and Rs /R0 ≃ 1.2. Mainly, we study the cases of different perturbation
wavenumber m and plasma β. It has been found that in the MHD case, the growth
rate of the instability reduces in proportion to the strength of the applied magnetic
field. The suppression mechanism is associated with the interference of two waves
running parallel and anti-parallel to the interface that transport vorticity and cause
the growth rate to oscillate in time with nearly a zero mean value
87
Figure 3.1: Schematic diagram of the physical setup: Two fluids of densities ρ1 and
ρ2 are separated by a perturbed interface located at r = R0 . A fast magnetosonic
MHD shock is initially located at r = Rs , with Mach number Ms .
88
3.2.2
Initial Conditions
The details of the setup of a converging cylindrical ideal fast magnetosonic shock may
be found in Reference [45]. We set up a planar MHD shock at a very large radius
and let this shock collapse onto the axis of convergence. When the shock reaches its
desired location at r = Rs , we extract the profiles and use that as the initial shock
profile. The incident converging shock may not be simply set up as a planar shock
because it is necessary to have correct radial variation of the flow variables behind
the shock. Figure 3.2 shows the initial profile of density, magnetic field, pressure,
and radial component of velocity of the base flow incident shock as a function of
radius. The perturbed interface in our model is affected by an impulse due to the
shock interaction together with the acceleration effects caused by the slow pressure
increase with radius behind the shock.
3.2.3
Governing Equations
It is expected that the diffusion time-scale would be much larger than the convective
time-scale in this converging RMI scenario. Hence we neglect viscous, thermal and
plasma resistivity effects. The flow is considered to be cylindrically symmetric with
non-zero velocity in the radial direction. The governing equations of ideal MHD are
presented in Chapter 2, Section 2.2.2 in the conservative form hence, not repeated
here. The initial magnetic field has an inverse radial dependence. Close to the origin,
large numerical errors are encountered close to the origin. These errors are eliminated
by rewriting the azimuthal component of the magnetic field as Bθ = Bθ∗ + Bθ1 where
β
is a curl-free portion of the magnetic field. The ideal MHD equations are
Bθ∗ =
r
89
(a) ρ̊
(b) B̊θ1
(c) p̊
(d) ůr
Figure 3.2: Initial conditions of the base flow variables (a) Incident shock position
(b) Modified magnetic field (c) Pressure and (d) Radial velocity.
90
rewritten
∂ W̃ 1 ∂(rF̃(W̃)) 1 ∂ G̃(W̃)
+
+
= S̃,
∂t
r
∂r
r ∂θ
(3.1)
where W̃ = {ρ, ρur , ρuθ , Br , Bθ1 , ẽ, ρφ}T . Also, the modified fluxes and source terms
are
F̃ = {ρur , ρu2r + p̃t − Br2 + Bθ∗ Bθ1 , ρur uθ − Br Bθ , 0, ur Bθ − uθ Br , (ẽ + p̃t + Bθ∗ Bθ1 )ur (3.2)
− (B.u)Br , ρur φ}T ,
G̃ = {ρuθ , ρur uθ − Br Bθ , ρu2θ + p̃t − (Bθ1 )2 − Bθ1 Bθ∗ , uθ Br
(3.3)
− ur Bθ , 0, (ẽ + p̃t + Bθ1 Bθ∗ )uθ − (B.u)Bθ , ρuθ φ}T ,
1
S̃ = − {0, (Bθ1 )2 − ρu2θ − p̃t + Bθ∗ Bθ1 , ρur uθ − Br Bθ , 0, uθ Br − ur Bθ , 0, 0}T , (3.4)
r
where ẽ =
3.2.4
p
+ 1 (B 2 + (Bθ1 )2 ) + 21 ρ(u2r + u2θ ) and p̃t = p + 21 (Br2 + (Bθ1 )2 ).
γ−1 2 r
Linear Stability Analysis
The linearization of the ideal MHD equations already has been discussed in Chapter
2. The examination of the linear stability of a converging shock interacting with a
perturbed interface is done by linearizing the ideal MHD equations by splitting the
solution into a background one-dimensional (radial) solution and perturbed solution
as W̃(r, θ, t) = W̊(r, t) + ǭŴ(r, t)eimθ . The time-dependent background state is a
function of the radial coordinate alone and comprises of MHD shocks, and a contact discontinuity. We term this background state as the “base”state (even though,
because it is time dependent, it is strictly speaking not a base state of standard
hydrodynamic stability analysis). Details of the base state are presented later. Substituting the base and perturbed solutions into the governing equations, we get the
91
partial differential equations governing the base and perturbed states as follows:
∂ W̊ 1 ∂
+
(rF(W̊)) = S̊,
∂t
r ∂r
(3.5a)
im
∂ Ŵ 1 ∂
+
(rM(W̊)Ŵ) +
N(W̊)Ŵ = Ŝ,
∂t
r ∂r
r
(3.5b)
where
0
∗ 1
(B̊ 1 )2 − ρ̊ů2 − ˚
p̃t + Bθ B̊θ
θ
θ
ρ̊ůr ůθ − B̊r B̊θ
1
S̊ = −
,
0
r
ů
B̊
−
ů
B̊
θ r
r θ
0
0
0
2B̊ 1 B̂θ − 2ρ̊ůθ ûθ − ρ̂ů2 − p̃ˆt + B ∗ B̂ 1
θ
θ
θ θ
ρ̊ůr ûθ + ρ̂ůr ůθ + ρ̊ûr ůθ − B̊r B̂θ − B̂r B̊θ
1
Ŝ = −
0
r
ûθ B̊r + ůθ B̂r − ûr B̊θ − ůr B̂θ
0
0
(3.6)
,
where M (W̊ ) and N (W̊ ) are the Jacobian matrices of the fluxes F̃ and G̃ with respect
to the base state. The total pressures ˚
p̃t , p̃ˆt are specified in terms of other variables
in equations (3.7) and (3.8).
γ−1
γ−2 2
˚
ρ̊(ů2r + ů2θ ) −
(B̊r + (B̊θ1 )2 ),
p̃t = (γ − 1)e̊ −
2
2
(3.7)
92
γ
−
1
p̃ˆt = (γ − 1)ê −
ρ̂(ů2r + ů2θ ) − (γ − 2)(B̊r B̂r + B̊θ1 B̂θ1 ).
2
(3.8)
We note here that the above equations, derived from ideal MHD equations, are scalefree and the choice of the radial domain can be arbitrary. Equations (3.5a) and (3.5b),
governing the base and perturbed state, respectively, are both hyperbolic partial
differential equations which are solved numerically with an explicit time marching
upwind procedure the details of which may be found in References [45, 57].
3.3
Results and Discussions
In this section, we present results from linear stability simulations. First, a convergence test is performed to quantify the order of accuracy of the numerical method.
Then, we examine the time-dependent radially varying background state solution,
followed by the investigation of azimuthal perturbation quantities for hydrodynamics
and MHD for both A > 0 and A < 0. Adopting the viewpoint that the main driving
mechanism of the RMI is the vorticity generated at the interface due to baroclinic
source term induced by the incident shock, we focus on the vorticity in the domain
and examine the suppression of RMI in terms of vorticity.
3.3.1
Base State Solution
Here we present the solution of the base state equations for MHD RMI. The timedependent base state equation (2.9a) are independent of the azimuthal wavenumber
m. It turns out that at t = 0.0 the density profile is mildly increasing with radius
just behind the shock (Fig. 3.3(a)). A similar behavior is seen for the radial velocity,
ůr , and the pressure, p̊, which are increasing in magnitude behind the shock front.
Next, the density profiles for positive, negative Atwood ratios (A > 0, A < 0) are
examined, respectively.
93
Positive Atwood Number
The base state is examined using the density profiles at time t = 0.0 (initial condition)
and t = 0.5 as shown in Fig. 3.3(a) for β = 16. The density spacetime diagram
illustrates the wave diagram of the base state more clearly. The contact discontinuity
(CD) can be seen as vertical line at initial position R0 in Fig. 3.3(b) indicating that it
is nearly stationary until the incident MHD shock (IS) overruns it. The IS bifurcates
into transmitted fast (TF), and a reflected fast (RF) MHD shocks. In one-dimensional
flow and the presence of azimuthal magnetic field with uθ = Br = 0, the slow and
intermediate speeds collapse to zero [45], and hence the slow-mode shocks are not
seen. The IS provides an impulse to the CD, thereby accelerating it. The CD is
located between the two shocks RF and TF. We note that at t = 0.5, the transmitted
fast shock is about to leave the domain so, all simulations are stopped at that time.
There is also an acceleration effect from the converging geometry, and due to the
pressure magnitude changing behind the shock.
Negative Atwood Number
The density profiles at time t = 0.0 (initial condition) and t = 0.2 are plotted in
Fig. 3.4(a) as solutions for the base state equations for β = 4. At t = 0.0, the CD is
seen as a vertical line at the initial position R0 , indicating that it is nearly stationary
until it is overrun by the IS. Then, the IS bifurcates into TF, and the fast mode
rarefaction wave (RR) moves in the opposite direction of the shock. The IS provides
an impulse to the CD which is located between the TF shock and RR waves. The
density spacetime diagram shows these waves clearly in Fig. 3.4(b). The small blip
in the profile discerned near r = 0.9 in density profile is attributed to a startup error
due to a sharp shock initialization. We note that at t = 0.2, the transmitted fast
94
(a) Density profiles at t = 0.0 and t = 0.5
(b) Spacetime diagram of density
Figure 3.3: (a) Density profiles of base state at t = 0.0 (initial conditions) and t = 0.5
for β = 16. IS denotes the incident shock, TF denotes the transmitted fast MHD
shock, RF denotes the reflected fast MHD shocks, and CD is the contact discontinuity
(b) Spacetime diagram of density field of the base state. IS denotes the incident
shock, TF denotes the transmitted fast MHD shock, RF denotes the reflected fast
MHD shock, and CD is the contact discontinuity.
95
shock is about to leave the domain and therefore, all simulations are stopped at that
time.
3.3.2
Evolution of the Perturbations
In this section, we present the results for both hydrodynamics (HD) (i.e., β = ∞),
and MHD (i.e., finite β) cases with azimuthal perturbations. In principle, both azimuthal and axial perturbation modes should also be examined, either separately and
in combination, for the sake of completeness. However, in Chapter 2, the inclusion of
both axial and azimuthal perturbations showed somewhat similar results. Hence, we
limit our investigation to azimuthal perturbations at this time. Following Lombardini
& Pullin [40], we non-dimensionalize the time scale using 1/(a0 (m/R0 )) as our reference time scale, where a0 is the sound speed in the unshocked medium ahead of the
incident shock. This time scale corresponds to the time taken by an acoustic wave
to traverse a wavelength of the unshocked interface. The perturbation growth rate ḣ
is normalized by the theoretical growth rate ḣ∞ was given in Chapter 2 in equation.
2.11 for the case of A > 0. On the other hand, for the case of A < 0 a modified
ḣ∞∗ is used to scale ḣ, where the post-shock initial amplitude, h+
0 , is replaced by the
+
−
h + h0
, where h−
average 0
0 is the pre-shock initial amplitude [25], and FMs = 1 for
2
this case. We note that there is indeed a small difference in ∆w between the MHD
and HD base states. However, in order to compare results between HD and MHD
cases with different magnetic field strengths, we use the same theoretical growth rate
in equation. ?? to normalize the growth rate in both HD and MHD cases.
Hydrodynamic Case
We begin our discussion with the HD case for both positive and negative Atwood
numbers followed by MHD results for the perturbation wavenumber m = 256. The
96
(a) Density profiles at t = 0.0 and t = 0.2
(b) Spacetime diagram of density
Figure 3.4: (a) Density profiles of base state at t = 0.0 (initial conditions) and t = 0.2
for β = 4. IS denotes the incident shock, TF denotes the transmitted fast MHD
shock, RR denotes the reflected rarefaction fan, and CD is the contact discontinuity
(b) Spacetime diagram of density field of the base state. IS denotes the incident
shock, TF denotes the transmitted fast MHD shock, RR denotes the rarefaction fan,
and CD is the contact discontinuity.
97
growth rate in HD is similar to that obtained in the Cartesian geometry case until
about scaled time a0 tm/R0 ≃ 30 as shown in Fig. 3.5. The growth rate increases
from zero and after a small time lag of the order of a few acoustic time periods
(the acoustic time period is the time taken by a sound wave to travel the distance
equal to the wavelength in the azimuthal direction) reaches a value of unity and then
oscillates around unity; this would correspond to the asymptotic growth rate obtained
in Cartesian planar geometry and this is the growth rate associated with the RMI.
The shock provides an impulsive acceleration to the interface. But in the cylindrical
case, after the shock passage, there is another acceleration due to the pressure rise
behind the shock (see Fig. 3.2(c)). Moreover, the interface continues to accelerate
(resp. decelerate) in the negative (resp. positive) Atwood ratio case for a certain
time duration which the acceleration vector points in the direction from the heavy
to the light fluid (i.e. RT unstable stratification). Therefore, in the cylindrical case,
as the interface approaches the origin, the RTI manifests itself, and we observe an
exponential increase in the magnitude of the growth rate. The RTI is associated
with a phase change for the positive Atwood interface. For the A > 0 interface,
this growth rate is opposite in sign to the model growth rate ḣ∞∗ . For the A < 0
interface, the RMI itself undergoes a phase change and the RTI reinforces this with
the corresponding exponential increase in growth rate as seen in Fig. 3.5. In reality,
the instability of the interface is a combination of both RM and RT instabilities.
However, the RMI phase dominates the early part of the growth rate curve, and the
RTI dominates the later portion of the growth rate curve. Hence the demarcation
between RMI and RTI phases in Fig. 3.5 is somewhat artificial; nonetheless after
the scaled time ≈ 30 the exponential growth of the RTI is evident. Another subtle
feature is that the frequency of oscillation of the growth rate increases with time: this
is because, although the wavenumber m is fixed, the wavelength 2πR/m decreases as
98
Figure 3.5: Scaled ḣ for HD and m = 256 for A < 0 and A > 0.
interface moves radially inwards.
Effect of Azimuthal Seed Magnetic Field
In this section, we investigate the effect of varying β parameter which characterizes
the strength of the seed magnetic field in the azimuthal direction, as well as the effect
of varying the azimuthal perturbation wavenumber m on the time history of scaled
growth rates in MHD. Note that β is a non-dimensional parameter that is proportional
to the square of the ratio of the acoustic speed in the unshocked incident gas and the
Alfvén speed. In our investigation, β is larger than unity so that the Alfvén speed is
always lower than the acoustic sound speed. First, the effect of varying the magnetic
field strength is presented in the left and right panel in Fig. 3.6(a), 3.6(b) for A > 0
and A < 0 interfaces, respectively, for β = ∞, 256, 128, 16, 4 for fixed wavenumber
m = 256. For finite β we note two frequencies: the higher frequency oscillations
99
are attributed to the acoustic waves traveling in the azimuthal direction. These are
clearly discerned in the β = ∞ or hydrodynamic cases discussed earlier and also
present for finite β cases (e.g. see the β = 128, 256 cases where these are more clear
especially at late times). For finite magnetic field strength or finite β, the oscillations
corresponding to the propagation of Alfvén waves in the azimuthal direction occur
at a lower frequency than those associated with acoustic modes. The frequency of
the Alfvén mode oscillations may be approximately predicted by the model proposed
by Wheatley et al. [21] for the planar case. An obvious effect of the azimuthal seed
field is that it suppresses the instability or growth of the interface. The larger the
magnetic field (smaller β) the larger the suppression: the amplitude of the growth
rate decreases with decreasing β and the oscillation frequency of the Alfvén mode
increases. Here, by instability we mean a combination of RMI and RTI and the
magnetic field affects both. Later in this section we will estimate the minimum field
strength required to suppress the RTI. For any given β, as time progresses, we see a
decrease in the amplitude of the growth rate. We note that the reduction of growth
rate amplitude with time is attributed to two reasons: (a) the dominant reason is that
as the density interface moves radially inwards the magnetic field strength increases in
proportion to 1/r and this reduces the growth amplitude; and (b) the second reason is
that our numerical method based on upwind schemes has implicit numerical diffusion
which results in a diffusion of the baroclinic vorticity away from the interface causing
a decrease in the growth rate of the interface perturbations.
As discussed for the hydrodynamic cases, there is a short time lag before the RMI
fully develops (or scaled growth reaches unity). For finite but relatively high β cases
(β ≥ 128) we note that the growth rate reaches unity before reducing. The lower
β cases see a suppression of the growth rate very early so that the scaled growth
rate never even approaches unity. We next elucidate the physical mechanism of the
100
(a) Effect of different β and m = 256 for A > 0
(b) Effect of different β and m = 256 for A < 0
Figure 3.6: Time history of growth rates for different β and fixed azimuthal wavenumber m = 256 (a) A > 0 (b) A < 0.
101
instability suppression. At the CD, for finite β, the vorticity distribution breaks
up into two vorticity carrying waves that propagate parallel and anti-parallel to the
magnetic field, causing the amplitude of interface to oscillate in time. This effect
is similar to that observed by Wheatley et al. [21] in planar nonlinear ideal MHD
simulations. Moreover, simulations presented here in cylindrical geometry show that
as the magnetic field strength increases the interface will oscillate with higher frequencies and smaller amplitudes. This can be due to the fact that vorticity carrying
waves propagate faster at the interface causing more rapid or higher frequency oscillations. Furthermore, as in the classic treatment of the RTI with a magnetic field
parallel to the interface where Chandrasekhar [58] suggests that the magnetic field
acts as a surface tension term, a higher field strength (lower β) is correlated with the
lower amplitude of the growth rate. These results also agree with Levy et al. [33]
who examined the cases of A < 0, A > 0 interfaces in the presence of parallel and
transverse magnetic field. They found that a higher magnetic field corresponds to a
longer oscillation period, and a smaller perturbation amplitude.
To further illustrate the vorticity dynamics in the vicinity of the interface, we com1 ∂(rûθ ) 1
− imûr ) (where Im(x) indicates
pute z-component of vorticity ω̂z = Im(
r ∂r
r
the imaginary value of x) evaluated at the interface for the case of A > 0, β = 16 and
m = 256. Using this, we recreate a two dimensional vorticity field in the vicinity of
the interface e to examine the vorticity as ωz (r, θ, t) = ω̂z (r, t) sin(mθ) in the range of
θ ∈ [−2π/m, 2π/m]. Figure 3.7(b-f) shows this two dimensional vorticity field (two
wavelengths shown) at scaled times t = 18, 23, 27, 33 and 40 which correspond to a
full oscillation cycle of the growth rate (shown in Fig. 3.7(a)). The dashed line in
the 2D vorticity images corresponds to the mean location of the interface. There
is a clear correspondence, somewhat obvious from a kinematics viewpoint, between
the maximum positive growth rate and positive vorticity on the interface (for half-
102
wavelength shown in θ ∈ [0, π/m]), zero growth when interface vorticity is zero on the
interface, and negative growth with negative vorticity on the interface (θ ∈ [0, π/m]).
We note that the dominant effect is that the vorticity on either side of the interface is
transported by waves propagating at the local Alfvén speed, parallel and anti-parallel
to the interface. The interference of these waves, with alternating phase causes constructive and destructive interference and vorticity patterns alternate with changing
signs. The net result is that the growth rate of the interface perturbation oscillates
from positive to negative. In addition, the periodic interference of these waves continues with time and causes the interface growth rate to oscillate. As mentioned earlier,
the process is very similar to that seen in planar nonlinear simulations of Wheatley
et al.(See Figure 6 in Reference [21]).
Next, for a constant magnetic field strength β = 4 we plot the time history
of scaled growth rates by varying the azimuthal perturbation wavenumber m =
16, 64, 128, 256 in Figures 3.8(a), 3.8(b) for A > 0, A < 0 interfaces, respectively.
The growth rates oscillate in time with a frequency that is proportional to m so
that on the scaled time plot the frequency is nearly the same for all wavenumbers.
However, the amplitude of the growth rate decreased by increasing the wavenumber,
the highest wavenumber case (m = 256) with the smallest amplitudes compared with
other cases. This indicates the level of suppression of the instability increases with the
wavenumber. Besides, for the MHD cases, the onset of RTI is not as well demarcated
from the RMI phase as it was for the hydrodynamic cases.
From the discussion earlier, it is evident that a magnetic field strength corresponding to β = 4 effectively suppresses the instability for both positive and negative
Atwood ratio interfaces, for all wavenumbers m ≥ 64 within several Alfvén wave transit times. We now examine the behavior of the instability for magnetic fields which
are considerably weaker (large β cases) than the one examined above. The objective
103
(a)
ḣ
for β = 16 and m = 256
ḣ∞
(b) t = 18
(c) t = 23
(d) t = 27
(e) t = 33
(f) t = 40
Figure 3.7: 2D Vorticity ωz in the interface for β = 16, m = 256 at different scaled
times chosen to correspond to a full oscillation period of the growth rate. The arrows
in the growth rate plot (a) Indicate the times chosen.
104
(a) Effect of different m and β = 4 for A > 0
(b) Effect of different m and β = 4 for A < 0
Figure 3.8:
(b) A < 0.
Time history of the growth rates for different m and β = 4 (a) A > 0
105
is to determine, albeit numerically, a critical value of the field as βcrit that would
effectively have a suppressing effect on the instability. We plot the growth rate in
Fig. 3.9(a), 3.9(b) for both Atwood ratio interface for large values of β for m = 256.
It is clear that β = 20, 000 has virtually no effect (apart from slightly delaying the
onset of RTI for A > 0). However, as the field strength increases, we observe a strong
effect at β = 1000.
The classic dispersion relation based on normal mode analysis (the solution is
proportional to exp(nt)) for magnetic Rayleigh-Taylor from Chandrasekhar [58] for
the case of a magnetic field that is parallel to the interface is expressed below for a
2D planar case:
2
n = gk
2B 2
ρ1 − ρ2
−
ρ1 + ρ2 (ρ1 + ρ2 )gk
,
(3.9)
where k is the perturbation wavenumber, and g is the acceleration in the direction
of the heavy-to-light fluid. This result strictly is not valid for our instability study,
because of the difference in geometry (planar vs. cylindrical), and because ours is
a combination RM/RT instability with no constant acceleration. Nonetheless, the
above result can be used for the A < 0 interface to estimate the strength of the
critical magnetic field Bcrit , which can suppress the instability at least during the RT
phase. For zero growth, we estimate this as:
2
Bcrit
=
+
(ρ+
1 − ρ2 )g(τ )R(τ )
,
2m
βcrit =
2p0
,
2
Bcrit
(3.10)
where the superscript + indicates post-shock density values, and τ is chosen as the
time when the RTI begins to manifest itself, g(τ ) is the acceleration experienced by
the interface located at radius R(τ ) at that time. Choosing scaled time τ = 20, at
which g(τ ) = 0.4, and R(τ ) = 0.8, and using the post-shock values as ρ+
1 = 1.9, and
ρ+
2 = 0.72, we obtain βcrit = 2712. An examination of Fig. 3.9(b) indicates that this is
106
a reasonable estimate. This estimation of the critical value of the magnetic field relies
on the knowledge of the acceleration and the post-shock density values: these can be
obtained from a simple one-dimensional simulation. Alternatively, if we simply use
pre-shock conditions, the only unknown is g(τ ): if this is estimated as the maximum
acceleration after the initial impulsive acceleration, we obtain βcrit = 3200.
107
(a) Effect of different β on RTI for A > 0, m = 256
(b) Effect of different β on RTI for A < 0, m = 256
Figure 3.9: Time history of growth rates for high β values, and fixed azimuthal
wavenumber m = 256 (a) A > 0 (b) A < 0.
108
3.4
Comparison Between Linear and Nonlinear Simulations
The nonlinear simulation for RMI are performed using the ideal MHD equations. The
numerical method is introduced in B. The applied magnetic field is aligned with the
mean interface. The boundary and initial conditions are identical to those applied in
the linearized simulations. The nonlinear simulation is conducted on a grid extending
azimuthally from θ = 0 to 2π/128. In the radial direction, the grid extends from
r = 0.3 to 1.3. An incident fast MHD shock is initialized at Rs = 1.05 and the initial
interface location is R0 = 1, as done for the linearized simulations. The computational
grid has 2048 cells in the radial direction for m = 16 (m = 128) and 512 cells in the
azimuthal direction. A comparison of the perturbation growth rate of a single mode
interface in the linear and nonlinear simulation for the wavenumbers m = 16, 128 and
the magnetic field β = 16 in Fig. 3.10(a,b), and for higher value of the magnetic field,
β = 4, with wavenumbers m = 16, 128 in Fig. 3.10(b-c). The growth rate comparison
is performed after the passage of the incident shock over the perturbed interface.
At early times, the nonlinear and linear growth rates show satisfactory agreement.
As time increases, the nonlinear effects become important which the worsens this
agreement. This behavior is more obvious for the lowest wavenumber case (m = 16)
particularly for the weaker magnetic field case (β = 16).
109
0.1
0.02
Noninear
Linear
0.015
0.05
0.01
0
0.005
0
-0.05
Noninear
Linear
-0.005
-0.1
-0.01
0
1
2
3
4
5
0
6
1
(a) β = 16, m = 16.
2
3
4
5
(b) β = 16, m = 128.
0.01
0.02
Noninear
Linear
Noninear
Linear
0.015
0.005
0.01
0.005
0
0
-0.005
-0.01
-0.005
-0.015
-0.02
-0.01
0
1
2
3
(c) β = 4, m = 16.
4
5
6
0
5
10
15
20
25
30
35
(d) β = 4, m = 128.
Figure 3.10: Comparison of growth rate between linear and nonlinear simulations
(a) m = 16, β = 16 (b) m = 128, β = 16 (c)m = 16, β = 4 (d) m = 128, β = 16.
40
110
Chapter 4
Incompressible Models of Magnetohydrodynamic
Richtmyer-Meshkov Instability in Cylindrical Geometry
4.1
Introduction
In Chapter 4, our goal is to understand the effect of a magnetic field on RMI of
sinusoidally perturbed density interface separating incompressible conducting fluids.
Richtmyer proposed that, by the time the transmitted and reflected waves have traveled a long distance from the interface compared to the wavelength of the perturbation, the subsequent motion at the interface can be considered incompressible [59].
Therefore, we propose an incompressible model which is faster in time of calculations
than the nonlinear and linear simulations that solve a full system of hyperbolic PDE.
It is also useful as a model of the early linear stages of RMI. We consider two cases
as above: (a) the magnetic field is oriented parallel to the motion of the shock and
(b) the case of a magnetic field parallel to the mean interface location. The initial
condition for this flow and the setup for the model problem are illustrated in Fig.
4.1 and 4.2 for each case respectively in section 4.2. We use an approach similar
to the one presented by Wheatley et al. [6, 21], so that the interface motion occurs
only due to the impulsive acceleration. In Section 4.3, the ideal incompressible MHD
equations are linearized which reduced to an initial value problem. The numerical
method presented in Section 4.4 used to find the solutions of the proposed models
. Then, the behavior predicted by this model is compared with the results of the
111
shock-driven simulations of the linearized compressible ideal MHD equations. these
results are shown in Section 4.5.
4.2
Incompressible Model Formulation
In this section, an incompressible model for the MHD RMI when the magnetic field
orientation is either normal or parallel to the interface is presented. The case of a
normal magnetic field is discussed first followed by the case of an azimuthal magnetic
field.
4.2.1
Normal Magnetic Field
In the case the magnetic field is oriented in the direction of the shock motion normal
to the contact discontinuity (CD). Figure 4.1(a) illustrates the initial condition for
a flow characterized by the density ratio across the CD, ρ2 /ρ1 , the ratio of specific
heats γ, the incident shock of Mach number Ms and the nondimensional strength
of the magnetic field, β = 2p0 /B 2 , where p0 is the initial pressure in the unshocked
regions of the flow and B is the magnitude of the magnetic field. The model for this
flow is shown in Fig. 4.1(b), where the shock is replaced by an impulsive acceleration
in r-direction, ∆V δ(t), where δ(t) is the Dirac delta function in time. The equations
are linearized about a base flow (denoted with a subscript 0) that results from the
impulsive acceleration of an unperturbed interface. This base flow has no θ dependence and zero azimuthal velocity uθ . Our choice of reference frame is attached to
the impulsively started contact discontinuity which results in the mean radial velocity
(ur ) at the contact being zero for all time. The base flow is governed by the following
equations
112
ρ̊(r) = ρ1 + H(r)(ρ2 − ρ1 ),
ůr = 0,
B̊r = br ,
ůθ = 0,
B̊θ = 0,
(4.1)
p̊(r, t) = −ρ1 ∆V δ(t)r − H(r − R0 )(ρ2 − ρ1 )∆V δ(t)r,
where ρ is the density, p is the pressure, u is the velocity, and B is the magnetic
field with initial constant magnitude br in the normal direction to the interface and
H(r − R0 ) is the Heaviside function. When the interface is perturbed, the density
becomes ρ̊(r −f ), where f (θ, t) is the location of the interface and f ≪ λθ , where λθ is
the wavelength. We assume the interface moves with a constant impulsive speed of the
shocked contact discontinuity, this will prevent the deceleration effect of converging
geometry on the interface as it converges to the origin. As a result, the RTI phase of
the instability is not captured by the model.
4.2.2
Azimuthal Magnetic Field
This section presents an incompressible model of the MHD RMI for a magnetic field
parallel to the mean interface location. The initial condition for a shocked sinusoidally
perturbed density interface is shown in Fig. 4.2(a), in the presence of an azimuthal
magnetic field B. The corresponding setup for the incompressible model is shown in
Fig. 4.2(b). An impulsive acceleration in r-direction, ∆V δ(t), is used to replace the
shock acceleration of uniform conducting fluids of different densities ρ1 and ρ2 . The
initial magnetic field is aligned with the mean interface location in θ-direction with
the nondimensional strength β = p0 /B 2 . The base flow is governed by the following
113
(a) Shock-driven setup
(b) incompressible model
Figure 4.1: (a) Initial condition for the shock-driven compressible MHD RMI (b)
Initial condition for the corresponding impulse driven, incompressible model.
114
equations:
ρ̊(r) = ρ1 + H(r)(ρ2 − ρ1 ),
B̊r = 0,
B̊θ = bt ,
ůr = 0,
ůθ = 0,
(4.2)
p̊(r, t) = p0 ,
where bt is the initial magnitude of the magnetic field in the azimuthal direction. We
linearize the ideal MHD equations about this flow.
115
(a) Shock-driven setup
(b) Incompressible model
Figure 4.2: (a) Initial condition for the shock-driven compressible MHD RMI (b)
Initial condition for the corresponding impulse driven incompressible model.
116
4.3
Linearized Initial Value Problem
We seek solutions to the model problem that satisfy the linearized equations of ideal
incompressible MHD. The perturbed equations are linearized about the base flow by
setting the density equal to
ρ(r, θ, t) = ρ1 + H[r − f (θ, t)](ρ2 − ρ1 ),
where f (θ, t) = η(t)eimθ ≪ λθ , and by assuming that all other flow quantities have the
form q ′ (r, θt) = q̊(r) + q ′ (r, θ, t), where q ′ are small perturbations to the base flow in
the form q ′ (x, z, t) = q̂(r, t)eimθ . Using these expressions in the ideal MHD equations
and neglecting terms of higher orders in perturbations, the linearized equations of the
perturbed state are given in the non-conservative form as
1 ∂(rûr ) im
+
ûθ = 0,
r ∂r
r
ρ
ρ
(4.3a)
∂ ûr 1 ∂(r[2ρûr ůr + p̂t − 2B̂r B̊r ]) im
+
+
[ρ(ûr ůθ + ůr ûθ ) − (B̂r B̊θ + B̊r B̂θ )] =
∂t
r
∂r
r
(4.3b)
1
− [2B̂θ B̊θ − 2ûθ ůθ − p̂t ] + (ρ2 − ρ1 )[(H(r) − H(r − f )]∆V δ(t),
r
∂ ûθ 1 ∂(r[ρ(ûr ůθ + ůr ûθ ) − (B̂r B̊θ + B̊r B̂θ )]) im
+
+
[2ρûθ ůθ + p̂t − 2B̂θ B̊θ ] =
∂t
r
∂r
r
(4.3c)
1
− [ρ(ûr ůθ + ůr ûθ ) − (B̂r B̊θ + B̊r B̂θ )],
r
∂ B̂r im
+
[ůθ B̂r + ûθ B̊r − (ůr B̂θ + ûr B̊θ )] = 0,
∂t
r
(4.3d)
117
∂ B̂θ 1 ∂r[ûr B̊θ + B̂θ ůr − (B̂r ůθ + B̊r ûθ )]
+
=
∂t
r
∂r
−1
[ůθ B̂r + ûθ B̊r − (ůr B̂θ + ûr B̊θ )],
r
∂(rB̂r )
+ imB̂θ = 0,
∂r
(4.3e)
(4.3f)
where the total pressure is p̂t = p̂ + B̂r B̊r + B̂θ B̊θ .
4.3.1
Normal Magnetic Field
Expanding the derivatives of equations (4.3) and use the base flow quantities (4.1),
the system (4.3) is simplified as follow
r
ρ
ρ
∂ ûr
+ ûr + imûθ = 0,
∂r
(4.4a)
imbr
∂ B̂r 2br
∂ ûr ∂ p̂t
B̂r −
B̂θ =
+
− 2br
−
∂t
∂r
∂r
r
r
(ρ2 − ρ1 )[(H(r) − H(r − f )]∆V δ(t),
(4.4b)
∂ ûθ
im
∂ B̂θ 2br
B̂θ +
− br
−
p̂t = 0,
∂t
∂r
r
r
(4.4c)
∂ B̂r imbr
+
ûθ = 0,
∂t
r
(4.4d)
∂ B̂θ
∂ ûθ
− br
= 0,
∂t
∂r
(4.4e)
118
r
∂ B̂r
+ B̂r + imB̂θ = 0.
∂r
(4.4f)
The initial conditions are imposed at t = 0− , just prior to the impulsive acceleration,
when the velocity and magnetic field perturbations are zero. The forcing of the
perturbations due to the impulse at t = 0+ is non-zero only in the vanishingly small
region r ∈ [R0 , f (θ, t)]. Taking the temporal Laplace transforms (LT) of (4.3) outside
the forcing region in each fluid gives the transformed homogeneous equations as
rDUj + Uj + imVj = 0,
ρj sUj + DPj − 2br DHrj −
ρj sVj − br DHθj −
imbr
2br
Hrj −
Hθj = 0,
r
r
im
2br
Hθj +
Pj = 0
r
r
(4.5a)
(4.5b)
(4.5c)
imbr
Vj = 0,
r
(4.5d)
sHθj − br DVj = 0,
(4.5e)
rDHrj + Hrj + imHθj = 0,
(4.5f)
sHrj +
where U, V, Hr , Hθ , and P are the temporal LT of ûr , ûθ , B̂r , B̂θ and p̂t respectively.
j = 1 and 2 for the two fluids, s is the Laplace variable and D ≡ d/dr. Then,
119
combining equations (4.5 a-f) gives a single fourth order ODE with a differential
1
1
operator L = L0 + L1 + 2 L2 as
r
r
"
2
!
2 2
#
sn
s
+ n2 D2 + 2 Uj
LUj ≡ D4 −
2
CAj
C
"
#
"
! # Aj
2
2
1
s
1
3s
+
4D3 −
+ 2n2 D Uj + 2 3D2 − 2 Uj = 0,
2
r
CAj
r
CAj
(4.6)
where, the differential operator L is a combination of the exact terms of 4th order
ODE of linear analysis of Wheatley et al. [6] in Cartesian geometry, L0 , plus the
1
1
corrections to the ODE stemming from cylindrical geometry, L1 + 2 L2 . Here,
r
r
√
n = m/r refers to the wavenumber and CAj = br / ρj is the Alfvén wave speed
in fluid j. The perturbations must be bounded as |r − R0 |→ ∞ and there are no
incoming waves from r = ±∞, where r is the distance from the interface. The
matching conditions across the interface r = R0 , refer to jump conditions of the
velocity normal component ûr , B̂r and B̂θ which must be continuous (see pages 458495 [58]). Taking Laplace transforms of these variables and using 4.5 to express each
in terms of U these boundary conditions become, to leading order in f ,
[U ]r=R0 = 0, [DU ]r=R0 = 0, and [D2 U ]r=R0 = 0,
(4.7)
where [q]R=R0 = q2 |r=R0 −q1 |r=R0 . Another matching condition is the dynamical condition, thatobtained by integrating equation (4.4 b) with respect to correction of r as
ξ from 0 to f across its inhomogeneous region, where the total pressure is represented
by the terms p̂t − 2br Br . The integration is contained within fluid 1, so the resulting
equation is
ρ1
Z
f
0
∂ ûr1
dξ + p̂1 (f, θ, t) − p̂1 (0, θ, t) + H.O.T =
∂t
Z
f
0
(ρ2 − ρ1 )vδ(t)dξ.
(4.8)
120
∂ ûr1
is considered smooth within the inhomogeneous
∂t
region, so it is approximated by
Another assumption is that
ρ1
Z
f
0
∂ ûr1
∂ ûr1
dξ h ρ1 f
|ξ=0 .
∂t
∂t
(4.9)
This term is O(f 2 ), and thus it can be neglected. Next, the continuity of p̂ across the
contact discontinuity ξ = f is expressed as [p̂] = 0. Thus to the leading order and
assuming p̂1 and p̂2 are smooth then
p̂1 (f, θ, t) = p̂2 (f, θ, t) h p̂2 (0, θ, t).
(4.10)
Combining 4.8 together with 4.9, 4.10 and neglecting the higher order terms we get:
p̂2 (0, θ, t) − p̂1 (0, θ, t) = (ρ2 − ρ1 )∆V δ(t)eimθ
(4.11)
Taking LT of the previous equation gives the transformed pressure magnitude as
follows:
[P ]ξ=0 = (ρ2 − ρ1 )∆V η0 .
(4.12)
The details of the numerical formulation of this condition are given in section (4.4).
4.3.2
Azimuthal Magnetic Field
For this case, we expand the derivatives in equations (4.3) and use the base flow
quantities (4.2). The resulting linearized equations are
r
∂ ûr
+ ûr + imûθ = 0,
∂r
(4.13a)
121
ρ
ρ
∂ ûr ∂ p̂t 2bt
imbt
+
+
B̂θ −
B̂r =
∂t
∂r
r
r
(ρ2 − ρ1 )[(H(r) − H(r − f )]∆V δ(t),
∂ ûθ im
∂ B̂r 2bt
2bt im
+
p̂t − bt
−
B̂r −
B̂θ = 0,
∂t
r
∂r
r
r
r
(4.13b)
(4.13c)
∂ B̂r imbt
+
ûr = 0,
∂t
r
(4.13d)
∂ B̂θ
∂ ûr
+ bt
= 0,
∂t
∂r
(4.13e)
∂ B̂r
+ B̂r + imB̂θ = 0.
∂r
(4.13f)
The initial conditions are imposed at t = 0− , just before the impulsive acceleration,
when the velocity and magnetic field perturbations are zero. Taking the temporal LT
of (4.13) outside the forcing region in each fluid gives
rDUj + Uj + imVj = 0,
ρj sUj + DPj +
ρj sVj − bt DHrj −
imbt
2bt
Hθj −
Hrj = 0,
r
r
2bt
imbt
im
Hrj −
Hθj +
Pj = 0,
r
r
r
(4.14a)
(4.14b)
(4.14c)
122
imbt
sHrj +
Uj = 0,
r
(4.14d)
sHθj + bt DUj = 0,
(4.14e)
rDHrj + Hrj + imHθj = 0.
(4.14f)
The equations (4.14 a-f) are combined which gives a second order ODE of the trans1
1
formed radial velocity magnitude in each fluid, with L = L0 + L1 + L2 + 2 L3
r
r
as
"
#
#
"
2
2
2
s
s
1
3s
LUj ≡ D2 − n Uj +
D2 − 2 Uj +
+ 1 DUj +
CA2 j n2
CAj
r CA2 j n2
"
#
1
s2
+ 1 Uj = 0
r2 CA2 j n2
2
(4.15)
Also in this case, the differential operator L is a combination of the 2nd order ODE of
linear analysis of Wheatley et al. [21], L0 in Cartesian geometry, plus the corrections
1
1
of the ODE that include the cylindrical geometry, L1 + L2 + 2 L3 . This ODE is
r
r
simplified as follows:
r 2 D 2 Uj + r
"
3s2 r2 + CA2 j m2
s2 r 2
+
CA2 j m2
#
DUj + 1 − m2 Uj = 0.
(4.16)
The solution is subjected to boundary and matching conditions which are similar
to the conditions presented in section 4.3.1. The perturbations must be bounded
as |r − R0 |→ ∞ and there are no incoming waves from r = ±∞, where r is the
distance from the interface. Across the interface r = R0 , the jump conditions must
satisfied as the velocity normal component ûr , B̂r and B̂θ are continuous. Taking
123
Laplace transforms of these variables and using 4.5 to express each in terms of U
these boundary conditions become, to leading order in f ,
[U ]r=R0 = 0, [DU ]r=R0 = 0, and [D2 U ]r=R0 = 0,
(4.17)
where [q]r=R0 = q2 |r=R0 −q1 |r=R0 . The last matching condition is obtained by integrating equation (4.13 b) with respect to correction of r as ξ from 0 to f across its
inhomogeneous region, where the total pressure is represented by the terms p̂t −2br Br .
Next, the continuity of p̂ across the contact discontinuity ξ = f is expressed as [p̂] = 0.
the final form of dynamical condition:
p̂2 (0, θ, t) − p̂1 (0, θ, t) = (ρ2 − ρ1 )∆V δ(t)eimθ
(4.18)
Taking LT of the previous equation gives the transformed pressure magnitude as
follows:
[P ]ξ=0 = (ρ2 − ρ1 )∆V η0 .
(4.19)
For the case of azimuthal perturbations, equation 4.15 has an analytical solution in
terms of a hypergeometric functions. This analytical solution has to be transformed
into the physical domain using inverse Laplace Transform, which turns out to be in
tractable so we use numerical means to solve the above system.
4.4
Numerical Method
In the cylindrical geometry the ODEs (4.6) and (4.15) are solved numerically since
they are not separable. We note here, the ODEs in the Cartesian geometry are
separable and analytically solved as seen in [6, 21], while the solution of the ODEs
124
is obtained numerically. Next, the numerical inversion of the Laplace transform,
(ILT) is obtained by Gaver-Stehfest method [60, 61, 62], which is based on deforming
the contour in the Bromwich integral to enclose all the singularities of the operator L
inside it. Then, applying the residual theorem to obtain ILT of U (s). This integration
is performed along the vertical line s = α in the complex plane, s ∈ C. The number
α ∈ ℜ is chosen so that the line s = α lies in the right of all singularities of the
contour. This method approximates the solution of ûr (r, t) by sequence of functions
at specific time t, given as
M
ln 2 X
wl U (sl ),
t l=1
(4.20)
where l = 1, 2, .., M, M is even, the expression of sl =
l ln 2
, and the coefficients wl
t
ûr (r, t) ≈
can be easily computed as
min(l,M/2)
wl = (−1)
(M/2)+l
X
k=⌊ 21 (1+l)⌋
2k! k ( M/2)
.
(M/2 − k)! k! (k − 1)! (l − k)! (2k − l)!
The numerical implementation of the pressure continuity is performed by combining equation (4.5 b) with (4.5 f), so we get:
ρj sUj + DPj = 0.
(4.21)
The Trapezoidal rule is used to approximate the integration of equation (4.21) w.r.t.
ξ at the interface position ξ = 0 in the interval [−a, a] as
[Pj ]|ξ=0 = −ρj s
Z
a
−a
Uj dξ = −ρj sd
l−1
X
Ul,1
U−l,2
+
Ui,j +
2
2
i=−l+1
!
.
(4.22)
125
Then we express this equation in each fluid as follows:
!
−1
X
U−l,2
U0,2
P2 |ξ=0 = −ρ2 sd
,
+
Ui,2 +
2
2
i=−l+1
!
l−1
Ul,1
U0,2 X
+
Ui,1 +
.
P1 |ξ=0 = −ρ1 sd
2
2
i=1
(4.23)
Substituting these equations into equation (4.12), we get the numerical jump condition at the interface. The one-sided central difference approximation is used for the
matching conditions [D2 U ]|ξ=0 = 0 and [DU ]|ξ=0 = 0. A uniform grid of 1024 points
with equal spacing ∆x is utilized in the numerical solution. The numerical discretization of the ODE’s is now describes The mesh has 2N + 1 points with equal spacing
cells, ∆x, so Ui,2 represents the velocity magnitude for i = −N, −N + 1, ..., −1 and
Ui,1 represents the velocity magnitude for i = 1, 2, ..., N − 1, N , j = 1, 2 and the
interface conditions are applied at U0 . The central second order difference is used for
discretization as listed below:
D4 Ui,j =
Ui+2,j − 4Ui+1,j + 6Ui,j − 4Ui−1,j + Ui−2,j
,
∆x4
(4.24a)
Ui+2,j − 2Ui+1,j + 2Ui−1,j − Ui−2,j
,
2∆x3
(4.24b)
Ui+1,j − 2Ui,j + Ui−1,j
,
∆x2
(4.24c)
Ui+1,j − Ui−1,j
.
2∆x
(4.24d)
D3 Ui,j =
D2 Ui,j =
DUi,j =
The set of linear equations are combined together and build the following linear
126
system:
AU = R,
(4.25)
where AN +1×N +1 is the coefficient matrix, the U is the solution vector and R is the
RHS that includes the impulsive velocity at the interface. The jump conditions at
i = 0 are discretized using the one-sided central difference first order accurate, listed
below:
U0− ,2 = U0+ ,1 ,
U1,1 − U0,1
−U−1,2 + U0,2
, DU0+ ,1 =
,
∆x
∆x
(4.26b)
U−2,2 − 2U−1,2 + U0,2 2
U1,1 − 2U0,1 + U2,1
, D U0+ ,1 =
.
2
∆x
∆x2
(4.26c)
DU0− ,2 =
D2 U0− ,2 =
4.5
(4.26a)
Results and Discussions
In this section, we present results of the incompressible model and the compressible
simulations of MHD RMI in planar and cylindrical geometries. First, in a planar
geometry the numerical result of the incompressible model [6] is presented and our
method is verified against the analytical solution of Wheatley et al.[6]. Following
the verification, we present the numerical solution of the incompressible model in
cylindrical geometry, and compare the results with compressible linear simulations of
MHD RMI for both cases of normal and azimuthal magnetic fields.
127
4.5.1
Verification of the Incompressible Model in Planar Geometry
The analytical solution of the incompressible model presented by Wheatley et al.[6]
is used to verify the numerical solution for same model. We use (x, z) as the planar
coordinate for this model. The linearized initial value problem of the incompressible
model was simplified by using LT and they obtained a 4th order ODE
D 4 Wj − (
s2 k 2
s2
2
2
+
k
)D
W
+
Wj = 0,
j
2
2
CAj
CAj
(4.27)
where Wj is the temporal LT of the perturbed horizontal velocity component ŵ , k
d
. The initial conditions correspond to the case of
is the wavenumber and D =
dz
accelerated interface by shock of strength Ms = 1.25 with ρ1 = 1.48372, ρ2 = 4.43149,
∆V = 0.319125, η0 = 0.00799 and β = 16.The boundary conditions used to solve
this equation are clarified in reference [6]. Figure 4.3 shows the profiles of ŵ(z, t)
at various times. At z = 0 the value of ŵ(z, t) represents the growth rate of the
interface, as t increases the maximum velocity peaks coincide with Alfvén fronts and
propagate away from the interface, therefore the growth rate decays to zero. The
analytical solution is shown in Fig. 4.3(a) from reference [6] where there are two
sharp cusps of Alfvén fronts. The comparison between numerical results and the
previous analytical results are quite reasonable although, dispersion errors stemming
from our purely central difference schemes affect the sharpness of the cusps and we
also get some undershoots. We can conclude that the solutions can reproduce the
analytical solution.
Next, the results of the incompressible model in cylindrical geometry are presented
and compared with the numerical simulations of the compressible MHD-RMI which
have been already presented in Chapters 2 and 3.
128
−3
normalized w perturbation amplitude
8
x 10
7
6
t / t* = 0
t / t* = 1
t / t* = 4
CD location
5
4
3
2
1
0
−1
−3
−2
−1
0
1
2
z/λ
ŵz
p"""""""""""""solution[6]
(a) Analytical
0.008
0.007
0.006
0.005
0.004
0.003
0.002
0.001
0
−0.001
t=0
t=1
t=4
interface
−3
−2
−1
0
z
1
2
3
(b) Present numerical solution
Figure 4.3: Profiles of velocity ŵ (a) the analytical solution (b) using the numerical
method.
129
4.5.2
Normal Magnetic Field in Cylindrical Geometry
The geometry of linear simulations is presented in Fig. 4.2(a). In the reference shockdriven linear simulations, the ratio of specific heats is set to γ = 5/3, the initial
density ratio is ρ2 /ρ1 = 3, the pressure upstream of the shock is p0 = 1, the shock
Mach number Ms = 1.25 and the normal magnetic field is characterized by β = 16.
The behavior of the present incompressible model is evaluated using, ρ1 = 1.4837, ρ2 =
4.4315 corresponding to post-shock densities. The magnitude of the impulse is set to
the interface velocity in the shock-driven problem as ∆V = 0.3191 with magnitude of
magnetic field β = 16. The solution of equation (4.6) is obtained numerically using
LU decomposition, then taking the numerical ILT to evaluate the perturbed velocity
in each fluid. The profile of ûr is plotted at different times in Fig. 4.4 where the
values at R0 = 1 are the growth rate of the interface. As in the previous observation
in the incompressible simulations of MHD-RMI, (Chapter 2 Fig. 2.13(a)), we see two
dominant Alfvén fronts coincide with the vortex sheets Fig. 2.13(e). These Alfvén
fronts are seen in the incompressible model solution which transport the vorticity
away from the interface position. They propagate with the local Alfvén speed. The
growth rate of the interface after the initial increase, decays to value less than zero.
For different wavenumbers, m = 64, 128 and β = 16 the results of scaled growth rate
of the interface are compared with the results of linear simulations and presented
in Fig. 4.5, the model matches the simulations at the early time. It is known in
converging geometry that as the interface converges to the origin, it decelerates. This
deceleration effect is not considered in the incompressible analysis because it assumes
a constant impulsive speed of the shocked contact discontinuity. As a result, the
RTI phase of the instability is not captured by the model. Nevertheless, it shows the
decaying growth rate behavior which represents the stabilizing effect of the magnetic
130
t=0
t = 0.23
t = 0.51
interface
0.8
ûr
0.6
0.4
0.2
0
0.6
0.7
0.8
0.9
r
1
1.1
1.2
1.3
Figure 4.4: Profiles of velocity ûr , using numerical ILT at t = 0, 0.23, 0.51, for
ρ1 = 1.48372, ρ2 = 4.4315, ∆V = 0.319125, β = 16.
field on RMI.
4.5.3
Azimuthal Magnetic Field in Cylindrical Geometry
In the presence of an azimuthal magnetic field, Fig. 4.2(a, b) show the geometry
of compressible linear simulations and the incompressible model for the RichtmyerMeshkov instability, respectively. The behavior of the compressible simulations and
the impulsively perturbed interface are examined using the initial values that presented in section 4.5.2. For the case of the incompressible model, the 2nd order ODE
(4.15) is solved numerically. The results are presented in regards to the growth rate
of the perturbed interface and compared with linear simulation as shown in Fig. (4.6)
for wavenumbers m = 32, 128 and β = 16. For large m, the model shows better results
as it is picking the right frequency and magnitude of the amplitude of the interface
ḣ/ḣ1
131
1
0.8
0.6
0.4
0.2
0
−0.2
−0.4
−0.6
model, m = 64
simulation,m = 64
0
10
20
30
40
50
ta0m/R0
(a) m=32
1
model, m = 128
simulation,m = 128
0.8
ḣ/ḣ1
0.6
0.4
0.2
0
−0.2
−0.4
0
10
20
30
40
50
ta0m/R0
(b) m=128
Figure 4.5: Comparison between the incompressible model and the linear simulations
for growth rate (a) m = 64, β = 16 (b) m = 128, β = 16.
132
and matches the linear simulation results at early scaled time < 30. The previous
discussion in chapter 3 (sec. 3.3.2) explained the dominant effect on either side of
the interface, the vorticity is transported by waves propagating at the local Alfvén
speed, parallel and anti-parallel to the interface. The interference of these waves,
with alternating phase causes constructive and destructive interference and vorticity
patterns alternate with changing signs. The net result is that the growth rate of the
interface perturbation oscillates from positive to negative. As mentioned earlier, the
process is very similar to that seen in planar nonlinear simulations of Wheatley et
al.(See Figure 6 in Reference [21]). An obvious effect of the azimuthal seed field is
that it suppresses the instability of the interface.
4.5.4
Hydrodynamic case
Earlier in section 1.2.2 we introduced the impulse model of Richtmyer, which is solved
analytically. Our goal is to introduce and solve the incompressible model for the HD
case in cylindrical geometry, in the absence of the magnetic field (i.e., β = ∞). This
is done by reducing the ODE’s 4.6 and 4.15 into the following equation
ρj s2 (D2 − n2 )Uj +
ρj s2
3ρj s2
DUj + 2 Uj = 0.
r
r
(4.28)
To investigate this model, the initial quantities presented in section 4.5.2 with wavenumber m = 128. The numerical method is used to solve this ODE and we compare the
solution with the relative numerical simulations, see Fig. 4.7. The model shows results that predicted by the impulse model, where the growth rate is a constant value
in time. It is important to emphasize that our model catch only the RMI as we
mentioned in section 4.5.2 not the RTI., while the numerical simulations shows the
growth rate increases from zero and then oscillates around unity. However, as the in-
133
1.5
model, m = 32
simulation, m = 32
1
ḣ/ḣ1
0.5
0
−0.5
−1
−1.5
−2
0
5
10
15
20
ta0m/R0
(a) m=32
model, m = 128
simulation, m = 128
1.5
ḣ/ḣ1
1
0.5
0
−0.5
−1
−1.5
0
10
20
30
40
50
ta0m/R0
(b) m=128
Figure 4.6: Comparison between impulsive model and numerical simulations for
scald growth rate (a) m = 32, β = 16 (b) m = 128, β = 16.
134
1.5
1
ḣ/ḣ1
0.5
0
−0.5
−1
−1.5
model, m = 256
simulation,m = 256
0
20
40
60
80 100 120 140
ta0m/R0
Figure 4.7: Scaled ḣ for HD and m = 128 for the model and numerical simulation.
terface approaches the origin, the RTI manifests itself and we observe an exponential
increase in magnitude of the growth rate.
135
Chapter 5
Concluding Remarks
5.1
Summary
The Richtmyer-Meshkov instability (RMI) occurs when a shock wave impulsively accelerates a perturbed density interface separating two different fluids [2, 16]. RMI
may not be classified as a classical hydrodynamic instability with modal analysis
because the perturbations grow algebraically (linearly in time) rather than exponentially. The RMI occurs in a variety of physical applications. It is important
in astrophysical phenomena such as supernova explosions, supersonic combustion in
hypersonic air-breathing engines, and inertial confinement fusion (ICF) [9, 63, 3].
The RMI in the presence of a magnetic field initially oriented parallel to the incident
shock direction has been investigated numerically and theoretically. Nonlinear ideal
magnetohydrodynamics (MHD) simulations of a shock interaction with an oblique
planar contact discontinuity were discussed by Samtaney [4]. It was demonstrated
that the RMI is suppressed by the existence of such a magnetic field. The primary
suppression mechanism may be attributed to transport of baroclinic vorticity away
from the density interface by slow- or intermediate mode shocks reducing the vorticity
at the interface. This results were reconfirmed in the context of incompressible MHD
linear stability analysis by Wheatley et al. [6, 20]. They demonstrated that vorticity
is transported away from the interface by Alfvén fronts.
As a first step, ideal MHD simulations of the RMI of a sinusoidally perturbed
136
interface were carried out, both in the presence and absence of a magnetic field. Then,
an incompressible model for ideal MHD RMI when the magnetic field is normal or
parallel to the mean interface location is formulated, which is solved using numerical
inverse Laplace transform.
In Chapter 2, linearized ideal MHD equations with a time-dependent base state
with discontinuities have been solved numerically. The numerical linear stability
analysis is performed in cylindrical geometry, with purely azimuthal or axial perturbations at a positive Atwood number interface interacts with the incident shock of
a Chisnell-type that propagates from a light gas to a heavy one. In hydrodynamics, the Richtmyer-Meshkov instability is driven by the baroclinic vorticity deposited
by the incident shock on the interface, followed by Rayleigh-Taylor instability which
dominates the growth rate. The transition from RMI to RTI in cylindrical geometry is accompanied by a phase change in the perturbation amplitude. Moreover,
this transition between RMI and RTI is most evident for the highest wavenumbers
investigated, and not so clearly demarcated for low wavenumbers. The growth rate
during the RMI phase is proportional to wavenumber (m or k) and agrees well, for
large wave numbers, with the asymptotic growth rate ḣ∞ of Lombardini and Pullin.
For the MHD cases of normal magnetic field, the suppression of the instability is
most evident for large wavenumbers and magnetic fields strengths characterized by
β = 4 or less, are sufficient to suppress the growth rate at the interface. There is a
residual small growth rate associated with the RTI. The mechanism of suppression is
the transport of vorticity away from the density interface by two Alfvén fronts.
In Chapter 3, we consider the case of a cylindrically collapsing fast MHD shock
interacting with an azimuthally perturbed interface. The linear simulations are carried out for cases where a fast magnetosonic converging MHD shock propagates from
a light gas to a heavy one or vice versa, corresponding to a positive/negative Atwood
137
number (A > 0/A < 0). In hydrodynamics, the baroclinic vorticity deposited by the
incident shock on the interface is the driving mechanism of the Richtmyer-Meshkov
instability. The RMI phase is somewhat short-lived and followed by Rayleigh-Taylor
instability that dominates the growth rate. For the MHD cases, the suppression of the
instability at the interface is most evident for large wavenumbers and relatively strong
magnetic fields strengths. After the shock interacts with the interface, the emerging
vorticity breaks up into waves traveling parallel and anti-parallel to the magnetic
field. The interference as these waves propagate with alternating phase causing the
perturbation growth rate of the interface to oscillate in time. Moreover, a critical
value of the field that would effectively have a suppressing effect on the instability is
determined numerically.
In Chapter 4, the linear analysis of an incompressible model MHD RMI when the
magnetic field has the normal or parallel direction to the interface is investigated. The
approach of Wheatley et al. [6, 21] is used so that the interface motion occurs only
due to the impulsive acceleration and not to the initial conditions. The linearized
initial value problem is solved by a numerical inverse Laplace transform. This model
is verified in terms of solving the incompressible model in planar geometry. The
predicted behavior of this model is compared with the results of the shock-driven
simulations of the linearized compressible ideal MHD equations. The results show that
the model with high wavenumbers are piking the same frequencies as the numerical
simulations. The model in cylindrical geometry does not capture the RT effect, as we
assumed the interface moves with constant speed.
138
5.2
Future Research Work
The present research can be further extended in many aspects: simulating flows with
multiple interfaces, applying other magnetic field configurations and investigating the
nonlinear development of these flows. We are also planning to extended our numerical
framework to handle two-fluid ideal MHD in which we distinguish between ions and
electrons [64]. Numerical challenges here stem from the finite speed of light, and the
electron-to-ion mass ratio.
139
REFERENCES
[1] R. Betti and O. A. Hurricane, “Inertial-confinement fusion with lasers,” Nature
Physics, vol. 12, no. 5, pp. 435–448, 2016.
[2] R. D. Richtmyer, “Taylor instability in shock acceleration of compressible fluids,” Communications on Pure and Applied Mathematics, vol. 13, no. 2, pp.
297–319, 1960.
[3] M. Brouillette, “The Richtmyer-Meshkov instability,” Annual Review of Fluid
Mechanics, vol. 34, pp. 445–468, 2002.
[4] R. Samtaney, “Suppression of the Richtmyer-Meshkov instability in the presence
of a magnetic field,” Physics of Fluids, vol. 15, no. 8, pp. L53–L56, 2003.
[5] ——, “A method to simulate linear stability of impulsively accelerated density
interfaces in ideal-MHD and gas dynamics,” Journal of Computational Physics,
vol. 228, no. 18, pp. 6773–6783, 2009.
[6] V. Wheatley, D. I. Pullin, and R. Samtaney, “Stability of an impulsively
accelerated density interface in magnetohydrodynamics,” Phys. Rev. Lett.,
vol. 95, p. 125002, 2005. [Online]. Available: http://link.aps.org/doi/10.1103/
PhysRevLett.95.125002
[7] S. Li and H. Li, “Parallel amr code for compressible MHD or HD equations,”
Los Alamos National Laboratory. Retrieved, pp. 09–05, 2006.
[8] D. H. Sharp, “An overview of Rayleigh-Taylor instability,” Physica D: Nonlinear Phenomena, vol. 12, no. 1, pp. 3–18, 1984.
[9] D. Arnett, “The role of mixing in astrophysics,” Astrophysical Journal Supplement Series, vol. 127, no. 2, pp. 213–217, 2000.
[10] Y. Yang, Q. Zhang, and D. H. Sharp, “Small amplitude theory of Richtmyer–
Meshkov instability,” Physics of Fluids, vol. 6, no. 5, pp. 1856–1873, 1994.
[11] A. M. Khokhlov, E. S. Oran, and G. O. Thomas, “Numerical simulation of
deflagration-to-detonation transition: the role of shock–flame interactions in
turbulent flames,” Combustion and Flame, vol. 117, no. 1, pp. 323–339, 1999.
140
[12] I. A. Waitz, F. E. Marble, and E. E. Zukoski, “Investigation of a contoured wall
injector for hypervelocity mixing augmentation,” AIAA journal, vol. 31, no. 6,
pp. 1014–1021, 1993.
[13] T. H. Johnson, “Inertial confinement fusion: Review and perspective,” Proceedings of the IEEE, vol. 72, no. 5, pp. 548–594, 1984.
[14] S. H. Glenzer, B. J. MacGowan, P. Michel, N. B. Meezan, L. J. Suter, S. N.
Dixit, J. L. Kline, G. A. Kyrala, D. K. Bradley, D. A. Callahan et al., “Symmetric inertial confinement fusion implosions at ultra-high laser energies,” Science,
vol. 327, no. 5970, pp. 1228–1231, 2010.
[15] J. Lindl, O. Landen, J. Edwards, E. Moses, N. Team et al., “Review of the
national ignition campaign 2009-2012,” Physics of Plasmas, vol. 21, no. 2, p.
020501, 2014.
[16] E. E. Meshkov, “Instability of the interface of two gases accelerated by a shock
wave,” Pure Applied Mathematics, vol. 4, no. 5, pp. 101–104, 1969.
[17] J. F. Hawley and N. J. Zabusky, “Vortex paradigm for shock-accelerated
density-stratified interfaces,” Physical Review Letters, vol. 63, pp. 1241–1244,
1989. [Online]. Available: http://link.aps.org/doi/10.1103/PhysRevLett.63.
1241
[18] J. Wu, H. Ma, and M. D. Zhou, Vorticity and vortex dynamics.
Science & Business Media, 2007.
Springer
[19] V. Wheatley, D. I. Pullin, and R. Samtaney, “Regular shock refraction at an
oblique planar density interface in magnetohydrodynamics,” Journal of Fluid
Mechanics, vol. 522, pp. 179–214, 2005.
[20] V. Wheatley, R. Samtaney, and D. I. Pullin, “The Richtmyer-Meshkov instability in magnetohydrodynamics,” Physics of Fluids, vol. 21, no. 8, p. 082102,
2009.
[21] V. Wheatley, R. Samtaney, D. I. Pullin, and R. M. Gehre, “The transverse field
Richtmyer-Meshkov instability in magnetohydrodynamics,” Physics of Fluids,
vol. 26, no. 1, p. 016102, 2014.
[22] T. Sano, K. Nishihara, C. Matsuoka, and T. Inoue, “Magnetic field amplification
associated with the Richtmyer-Meshkov instability,” Astrophysical Journal, vol.
758, no. 2, 2012.
[23] M. Hohenberger, P. Y. Chang, G. Fiksel, J. P. Knauer, R. Betti, F. J. Marshall,
D. D. Meyerhofer, F. H. Seguin, and R. D. Petrasso, “Inertial confinement
141
fusion implosions with imposed magnetic field compression using the OMEGAI
Laser,” Physics of Plasmas, vol. 19, no. 5, 2012.
[24] J. Cao, Z. Wu, H. Ren, and D. Li, “Effects of shear flow and transverse magnetic
field on Richtmyer–Meshkov instability,” Physics of Plasmas, vol. 15, no. 4, p.
042102, 2008.
[25] K. A. Meyer and P. J. Blewett, “Numerical investigation of the stability of a
shock-accelerated interface between two fluids,” Physics of Fluids, vol. 15, no. 5,
pp. 753–759, 1972.
[26] G. Fraley, “Rayleigh–taylor stability for a normal shock wave–density discontinuity interaction,” Physics of Fluids, vol. 29, no. 2, pp. 376–386, 1986.
[27] Q. Zhang and S. I. Sohn, “An analytical nonlinear theory of Richtmyer-Meshkov
instability,” Physics Letters A, vol. 212, no. 3, pp. 149–155, 1996.
[28] A. L. Velikovich, “Analytic theory of richtmyer–meshkov instability for the case
of reflected rarefaction wave,” Physics of Fluids, vol. 8, no. 6, pp. 1666–1679,
1996.
[29] R. Samtaney and D. I. Meiron, “Hypervelocity Richtmyer-Meshkov instability,”
Physics of Fluids, vol. 9, no. 6, pp. 1783–1803, 1997.
[30] M. A. Jones and J. W. Jacobs, “A membraneless experiment for the study of
Richtmyer-Meshkov instability of a shock-accelerated gas interface,” Physics of
Fluids, vol. 9, no. 10, pp. 3078–3085, 1997.
[31] K. Nishihara, J. G. Wouchuk, C. Matsuoka, R. Ishizaki, and V. V. Zhakhovsky,
“Richtmyer-Meshkov instability: theory of linear and nonlinear evolution,”
Philosophical Transactions of the Royal Society of London A: Mathematical,
Physical and Engineering Sciences, vol. 368, no. 1916, pp. 1769–1807, 2010.
[32] Z. Qiu, Z. Wu, J. Cao, and D. Li, “Effects of transverse magnetic field and
viscosity on the Richtmyer-Meshkov instability,” Physics of Plasmas, vol. 15,
no. 4, pp. 42 305–42 305, 2008.
[33] Y. Levy, S. Jaouen, and B. Canaud, “Numerical investigation of magnetic
Richtmyer-Meshkov instability,” Laser and Particle Beams, vol. 30, no. 03, pp.
415–419, 2012.
[34] G. I. Bell, “Taylor instability on cylinders and spheres in the small amplitude
approximation,” Los Alamos National Laboratory, Los Alamos, NM, Report
LA-1321, 1951.
142
[35] M. S. Plesset, “On the stability of fluid flows with spherical symmetry,” Journal
of Applied Physics, vol. 25, no. 1, pp. 96–98, 1954.
[36] K. O. Mikaelian, “Rayleigh-Taylor and Richtmyer-Meshkov instabilities and
mixing in stratified spherical shells,” Physical Review A, vol. 42, no. 6, p. 3400,
1990.
[37] S. Hosseini and K. Takayama, “Experimental study of Richtmyer-Meshkov instability induced by cylindrical shock waves,” Physics of Fluids, vol. 17, no. 8,
p. 084101, 2005.
[38] Q. Zhang and M. J. Graham, “A numerical study of Richtmyer-Meshkov instability driven by cylindrical shocks,” Physics of Fluids, vol. 10, no. 4, pp.
974–992, 1998.
[39] M. Lombardini, Richtmyer-Meshkov instability in converging geometries. California Institute of Technology, 2008.
[40] M. Lombardini and D. I. Pullin, “Small-amplitude perturbations in the threedimensional cylindrical Richtmyer-Meshkov instability,” Physics of Fluids,
vol. 21, no. 11, p. 114103, 2009.
[41] K. O. Mikaelian, “Rayleigh-Taylor and Richtmyer-Meshkov instabilities and
mixing in stratified cylindrical shells,” Physics of Fluids, vol. 17, no. 9, p.
094105, 2005.
[42] R. F. Chisnell, “An analytic description of converging shock waves,” Journal of
Fluid Mechanics, vol. 354, pp. 357–375, 1998.
[43] H. Yu and D. Livescu, “Rayleigh–taylor instability in cylindrical geometry with
compressible fluids,” Physics of Fluids, vol. 20, no. 10, p. 104103, 2008.
[44] J. Glimm, J. Grove, and Y. Zhang, “Numerical calculation of Rayleigh-Taylor
and Richtmyer-Meshkov instabilities for three dimensional axi-symmetric flows
in cylindrical and spherical geometries,” Preprint, SUNY at Stony Brook, 1999.
[45] D. I. Pullin, W. Mostert, V. Wheatley, and R. Samtaney, “Converging cylindrical shocks in ideal magnetohydrodynamics,” Physics of Fluids, vol. 26, no. 9,
p. 097103, 2014.
[46] G. B. Whitham, “On the propagation of shock waves through regions of nonuniform area or flow,” Journal of Fluid Mechanics, vol. 4, no. 04, pp. 337–360,
1958.
143
[47] W. Mostert, V. Wheatley, R. Samtaney, and D. I. Pullin, “Effects of seed
magnetic fields on magnetohydrodynamic implosion structure and dynamics,”
Physics of Fluids, vol. 26, no. 12, p. 126102, 2014.
[48] W. Mostert, D. I. Pullin, R. Samtaney, and V. Wheatley, “Converging cylindrical magnetohydrodynamic shock collapse onto a power-law-varying line current,” Journal of Fluid Mechanics, vol. 793, pp. 414–443, 2016.
[49] W. Mostert, V. Wheatley, R. Samtaney, and D. I. Pullin, “Effects of magnetic
fields on magnetohydrodynamic cylindrical and spherical Richtmyer-Meshkov
instability,” Physics of Fluids, vol. 27, no. 10, p. 104102, 2015.
[50] K. O. Mikaelian, “Turbulent mixing generated by Rayleigh-Taylor and
Richtmyer-Meshkov instabilities,” Physica D: Nonlinear Phenomena, vol. 36,
no. 3, pp. 343–357, 1989.
[51] R. L. Holmes, G. Dimonte, B. Fryxell, M. L. Gittings, J. W. Grove, M. Schneider, D. H. Sharp, A. L. Velikovich, R. P. Weaver, and Q. Zhang, “RichtmyerMeshkov instability growth: experiment, simulation and theory,” Journal of
Fluid Mechanics, vol. 389, pp. 55–79, 1999.
[52] N. J. Zabusky and S. Zhang, “Shock–planar curtain interactions in two dimensions: Emergence of vortex double layers, vortex projectiles, and decaying
stratified turbulence,” Physics of Fluids, vol. 14, no. 1, pp. 419–422, 2002.
[53] P. R. Chapman and J. W. Jacobs, “Experiments on the three-dimensional incompressible Richtmyer-Meshkov instability,” Physics of Fluids, vol. 18, no. 7,
p. 074101, 2006.
[54] T. Sano, T. Inoue, and K. Nishihara, “Critical magnetic field strength for suppression of the Richtmyer-Meshkov instability in plasmas,” Physical review letters, vol. 111, no. 20, p. 205001, 2013.
[55] R. Samtaney and N. J. Zabusky, “Circulation deposition on shock-accelerated
planar and curved density-stratified interfaces: models and scaling laws,”
Journal of Fluid Mechanics, vol. 269, pp. 45–78, 1994. [Online]. Available:
http://journals.cambridge.org/article S0022112094001485
[56] R. Samtaney, P. Colella, T. J. Ligocki, D. F. Martin, and S. C. Jardin, “An
adaptive mesh semi-implicit conservative unsplit method for resistive mhd,”
vol. 16, no. 1, p. 40, 2005.
144
[57] A. Bakhsh, S. Gao, R. Samtaney, and V. Wheatley, “Linear simulations of the
cylindrical Richtmyer-Meshkov instability in magnetohydrodynamics,” Physics
of Fluids, vol. 28, no. 3, p. 034106, 2016.
[58] S. Chandrasekhar, Hydrodynamic and hydromagnetic stability, ser. Dover
classics of science and mathematics. New York: Dover Publications, cop.
1961, 1981. [Online]. Available: http://opac.inria.fr/record=b1090797
[59] M. Brouillette and R. Bonazza, “Experiments on the Richtmyer–Meshkov instability: Wall effects and wave phenomena,” Physics of Fluids, vol. 11, no. 5,
pp. 1127–1142, 1999.
[60] A. Davies, “The solution of differential equations using numerical laplace transforms,” International Journal of Mathematical Education in Science and Technology, vol. 30, no. 1, pp. 65–79, 1999.
[61] A. Kuznetsov, “On the convergence of the gaver–stehfest algorithm,” SIAM
Journal on Numerical Analysis, vol. 51, no. 6, pp. 2984–2998, 2013.
[62] W. Abate, J.and Whitt, “A unified framework for numerically inverting laplace
transforms,” INFORMS Journal on Computing, vol. 18, no. 4, pp. 408–421,
2006.
[63] J. Yang, T. Kubota, and E. E. Zukoski, “Applications of shock-induced mixing
to supersonic combustion,” AIAA journal, vol. 31, no. 5, pp. 854–862, 1993.
[64] D. Bond, V. Wheatley, R. Samtaney, and D. I. Pullin, “Richtmyer–meshkov
instability of a thermal interface in a two-fluid plasma,” Journal of Fluid Mechanics, vol. 833, pp. 332–363, 2017.
[65] N. F. Ponchaut, H. G. Hornung, D. I. Pullin, and C. A. Mouton, “On imploding cylindrical and spherical shock waves in a perfect gas,” Journal of Fluid
Mechanics, vol. 560, pp. 103–122, 2006.
[66] A. Bakhsh and R. Samtaney, “Linear analysis of converging richtmyer–meshkov
instability in the presence of an azimuthal magnetic field,” Journal of Fluids
Engineering, vol. 140, no. 5, p. 050901, 2018.
[67] S. Gao, “Linear simulations of the cylindrical Richtmyer-Meshkov instability in
hydrodynamics and MHD,” MS Thesis, King Abdullah University of Science
and Technology, 2013.
[68] D. J. Hill, C. Pantano, and D. I. Pullin, “Large-eddy simulation and multiscale
modelling of a Richtmyer–Meshkov instability with reshock,” Journal of fluid
mechanics, vol. 557, pp. 29–61, 2006.
145
[69] J. McFarland, D. Reilly, S. Creel, C. McDonald, T. Finn, and D. Ranjan, “Experimental investigation of the inclined interface richtmyer–meshkov instability
before and after reshock,” Experiments in fluids, vol. 55, no. 1, pp. 1–14, 2014.
[70] A. Bakhsh and R. Samtaney, “Linear analysis of converging richtmyer-meshkov
instability in the presence of an azimuthal magnetic field,” vol. 140, 11 2017.
[71] J. U. Brackbill and D. C. Barnes, “The effect of nonzero ∇ · b on the numerical
solution of the magnetohydrodynamic equations,” Journal of Computational
Physics, vol. 35, no. 3, pp. 426–430, 1980.
[72] M. Latini, O. Schilling, and W. S. Don, “Effects of weno flux reconstruction
order and spatial resolution on reshocked two-dimensional richtmyer–meshkov
instability,” Journal of Computational Physics, vol. 221, no. 2, pp. 805–836,
2007.
[73] J. Glimm, J. Grove, Y. Zhang, and S. Dutta, “Numerical study of axisymmetric Richtmyer-Meshkov instability and azimuthal effect on spherical mixing,”
Journal of statistical physics, vol. 107, no. 1-2, pp. 241–260, 2002.
[74] A. L. Zachary, A. Malagoli, and P. Colella, “A higher-order godunov method for
multidimensional ideal magnetohydrodynamics,” SIAM Journal on Scientific
Computing, vol. 15, no. 2, pp. 263–284, 1994.
[75] S. K. Godunov, “Symmetric form of the equations of magnetohydrodynamics,”
Numerical Methods for Mechanics of Continuum Medium, vol. 1, pp. 26–34,
1972.
[76] Y. Kuramitsu, Y. Sakawa, T. Morita, S. Dono, H. Aoki, H. Tanji, C. D. Gregory,
J. N. Waugh, B. Loupias, M. Koenig et al., “Formation of density inhomogeneity in laser produced plasmas for a test bed of magnetic field amplification
in supernova remnants,” Astrophysics and Space Science, vol. 336, no. 1, pp.
269–272, 2011.
[77] R. Krechetnikov, “Rayleigh-Taylor and Richtmyer-Meshkov instabilities of flat
and curved interfaces,” Journal of Fluid Mechanics, vol. 625, pp. 387–410, 2009.
[78] R. K. Crockett, P. Colella, R. T. Fisher, R. I. Klein, and C. F. McKee, “An
unsplit, cell-centered godunov method for ideal mhd,” Journal of Computational
Physics, vol. 203, no. 2, pp. 422–448, 2005.
[79] J. D. Lindl, R. L. Mccrory, and E. M. Campbell, “Progress toward ignition and
burn propagation in inertial confinement fusion,” Physics Today, vol. 45, no. 9,
pp. 32–40, 1992.
146
[80] M. R. Mankbadi and S. Balachandar, “Compressible inviscid instability of
rapidly expanding spherical material interfaces,” Physics of Fluids, vol. 24,
no. 3, p. 034106, 2012.
[81] S. Falle, S. S. Komissarov, and P. Joarder, “A multidimensional upwind scheme
for magnetohydrodynamics,” Monthly Notices of the Royal Astronomical Society, vol. 297, no. 1, pp. 265–277, 1998.
[82] J. Kim, D. Ryu, T. W. Jones, and S. S. Hong, “A multidimensional code
for isothermal magnetohydrodynamic flows in astrophysics,” The Astrophysical Journal, vol. 514, no. 1, p. 506, 1999.
[83] J. H. Kim, “Small amplitude theory of Richtmyer-Meshkov instability in cylindrical and spherical geometries,” 1997.
[84] D. Ryu and T. W. Jones, “Numerical magetohydrodynamics in astrophysics:
Algorithm and tests for one-dimensional flow,” The Astrophysical Journal, vol.
442, pp. 228–258, 1995.
[85] M. Lombardini, D. I. Pullin, and D. I. Meiron, “Turbulent mixing driven by
spherical implosions. part 1. flow description and mixing-layer growth,” Journal
of Fluid Mechanics, vol. 748, pp. 85–112, 2014.
[86] S. Gao, “Linear simulations of the cylindrical Richtmyer-Meshkov instability in
hydrodynamics and mhd,” Ph.D. dissertation, 2013.
[87] O. Schilling and J. W. Jacobs, “Richtmyer-Meshkov instability and reaccelerated inhomogeneous flows,” Scholarpedia, vol. 3, no. 7, p. 6090, 2008.
[88] K. G. Powell, P. L. Roe, T. J. Linde, T. I. Gombosi, and D. L. De Zeeuw, “A
solution-adaptive upwind scheme for ideal magnetohydrodynamics,” Journal of
Computational Physics, vol. 154, no. 2, pp. 284–309, 1999.
[89] A. Ghasemizad, H. Zarringhalam, and L. Gholamzadeh, “the investigation of
rayleigh-taylor instability growth rate in inertial confinement fusion,” Fusion
Res Ser, vol. 8, pp. 1234–1238, 2009.
[90] S. Li and H. Li, “A modern code for solving magnetohydrodynamics or hydrodynamic equations,” submitted to ApJS, 2003.
[91] V. Wheatley, R. Samtaney, and D. I. Pullin, “The magnetohydrodynamic
Richtmyer-Meshkov instability: The transverse field case,” in Proceedings of
the 18th Australasian Fluid Mechanics Conference, Launceston, Australia, 2012,
pp. 3–7.
147
[92] V. A. Thomas and R. J. Kares, “Drive asymmetry and the origin of turbulence
in an icf implosion,” Physical review letters, vol. 109, no. 7, p. 075004, 2012.
[93] R. J. Stalker and K. C. A. Crane, “Driver gas contamination in a high-enthalpy
reflected shock tunnel,” AIAA Journal, vol. 16, no. 3, pp. 277–279, 1978.
[94] A. L. Velikovich, J. P. Dahlburg, A. J. Schmitt, J. H. Gardner, L. Phillips, F. L.
Cochran, Y. K. . Chong, G. Dimonte, and N. Metzler, “Richtmyer–Meshkovlike instabilities and early-time perturbation growth in laser targets and z-pinch
loads,” Physics of Plasmas, vol. 7, no. 5, pp. 1662–1671, 2000.
[95] A. L. Velikovich, J. G. Wouchuk, C. H. R. de Lira, N. Metzler, S. Zalesak, and
A. J. Schmitt, “Shock front distortion and Richtmyer-Meshkov-type growth
caused by a small preshock nonuniformity,” Physics of Plasmas, vol. 14, no. 7,
p. 072706, 2007.
[96] L. Massa and P. Jha, “Linear analysis of the Richtmyer-Meshkov instability in
shock-flame interactions,” Physics of Fluids, vol. 24, no. 5, p. 056101, 2012.
[97] A. Shiferaw and R. C. Mittal, “Fast finite difference solutions of the three
dimensional poissons equation in cylindrical coordinates,” American Journal of
Computational Mathematics, vol. 3, no. 04, p. 356, 2013.
[98] J. Lindl, “Development of the indirect-drive approach to inertial confinement
fusion and the target physics basis for ignition and gain,” Physics of Plasmas,
vol. 2, no. 11, pp. 3933–4024, 1995.
[99] R. Samtaney and D. Meiron, “Physics of the strong shock Richtmyer–Meshkov
instability,” in Proceedings of Fifth International Workshop on Compressible
Turbulent Mixing at Stony Brook, USA, 1995, pp. 373–380.
[100] M. Vetter and B. Sturtevant, “Experiments on the Richtmyer-Meshkov instability of an air/sf6 interface,” Shock Waves, vol. 4, no. 5, pp. 247–252, 1995.
[101] L. F. Henderson, “General laws for propagation of shock waves through matter,”
Handbook of Shock Waves, vol. 1, pp. 144–183, 2001.
148
APPENDICES
A
Numerical Method
In this section, we describe the numerical method that used to solve the ideal MHD
equation 2.9, following the method originally developed by Samtaney[5]. We use the
method-of-lines approach for the time discretization by rewriting, equation 2.9a as
∂ Ŵ
∂ W̊
= R(W̊ ),
= R(Ŵ ),
∂t
∂t
and adopting the third order TVD Runge-kutta time discretization approach as
W̊0 = W̊ n ,
W̊1 = W̊0 + ∆tR(W̊0 ),
3
W̊2 = W̊0 +
4
1
W̊3 = W̊0 +
3
1
W̊1 +
4
2
W̊2 +
3
∆t
R(W̊1 ),
4
2∆t
R(W̊2 ),
3
(A.1)
W̊ n+1 = W̊3 ,
where ∆t is the time step and W̊ n is the solution at time step n. The perturbed state
equations 2.9b are similarly updated.
A finite volume method is adopted for the spatial discretization in which we define
the solution vectors at the centroids of the subintervals indexed as j. The fluxes are
computed at the faces indexed as j ±
1
2
and then used to update the solution at the
149
centroids.
The Jacobian matrix M (W̊ ) is singular and hence we adopt the seven-wave method
to compute the fluxes [56]. Denoting the remaining seven variables and their fluxes
by an over-bar, we express the flux of the base state and perturbed quantities as
7
F̄ (W̊ )
j+ 12
¯
[M̄ (W̊ )Ŵ ]
1
1X
¯
¯
α k rk
=
(F̄ (W̊L,j+ 1 , B̊r ) − F̄ (W̊R,j+ 1 , B̊r )) −
2
2
2
2 k=1
(A.2)
1
1
¯
¯
= (M̄j+ 1 ŴL,j+ 1 + M̄j+ 1 ŴR,j+ 1 ) − M̄ j+ 1
2
2
2
2
2
2
2
¯
¯
(ŴR,j+ 1 − ŴL,j+ 1 )
(A.3)
j+ 12
2
2
where M̄ = M̄ + − M̄ − , lk and rk are the left and right eigenvectors of M̄ (M̄ is
¯
¯
M after deleting the row/column corresponding to B̊r ), αk = lk · (W̊R,j+ 1 − W̊L,j+ 1 ).
2
At each side of interface j +
1
,
2
2
both the base state and the perturbed quantities
¯
¯
(W̊L/R,j+ 1 and ŴL/R,j+ 1 ) are obtained by fitting linear profiles and slope limiting.
2
2
In the seven-wave method, the contribution B̂r with other variables is considered by
including the following term to the RHS of Eq. (2.9b)
1 ∂(rB̂r1 )
(A.4)
r ∂r
C(W̊ , B̂r ) = [0, γ B̊r1 + Br∗ , B̊θ , B̊z , ůθ , ůz , (γ B̊r1 ůr + B̊θ ůθ + B̊z ůz )]T
At the final stage of each time step, maintaining the solenoidal property of the
magnetic field requires updating the perturbed variables B̂θ and B̂z as:
• for purely azimuthal perturbation:
B̂θ
j
i ∂
i
(rB̂r1 ) ≈
=
m ∂r
m
j
rB̂r1
j+1
− rB̂r1
2∆r
j−1
.
(A.5)
150
• For purely axial perturbations:
B̂z
j
i
i ∂
(rB̂r1 ) ≈
=
krj ∂r
krj
j
rB̂r1
j+1
− rB̂r1
j−1
2∆r
.
(A.6)
All of B̂r1 , B̂θ and B̂z in the above equations are at the same time step.
A.0.1
Convergence Test
Numerical Convergence:
In Chapter 2, we have not presented detailed results on numerical convergence and order of accuracy of the proposed method. Numerical convergence tests of this method,
performed by Gao [67], show second order accuracy for smooth initial conditions.
The method proposed here is very similar to that developed by Samtaney [5] for
Cartesian slab geometry, which exhibits second order convergence for smooth flows
and convergence rates between one and two for flows with shocks. In Section 2.3.5,
linear simulations for m = 128, β = 16 (Fig.
2.20) show a small difference in the
growth rate at two mesh resolutions of 1600 and 8000 radial cells. Other such convergence tests were also performed for other parameters and showed little difference
in the growth rate at mesh resolutions of 1600 and 3200 radial cells.
In Chapter 3, we examine the convergence of our numerical results, by choosing
a test case with plasma beta β = 16, and perturbation wavenumber m = 128 and
A = 0.667. The solution is computed with different mesh sizes: 800, 1600, 3200, 6400
and 12800 zones in the radial domain. The unscaled growth rate of the perturbed
interface, (essentially the time derivative of the perturbation amplitude) ḣ, computed
as the real value Re(ûr ) at the interface position, is plotted in Fig. A.1(a). The
finest mesh resolution, 12800 cells, is designated as the “exact” numerical solution
for our simulations. We fix the time t = 0.3 to obtain the relative error of the
151
growth rates which is plotted in Fig. A.1(b). It shows that as the mesh resolution
increased the error is decreased which implies that the growth rates are converging.
The order of accuracy is 1.4914, which is between the first and second order: a value
between one and two is generally considered acceptable for flows dominated by shocks
and contact discontinuities. In ideal MHD, with the lack of a physical dissipation,
one cannot expect a point-wise convergence of the solution with grid refinement.
Numerical diffusion decreases as the grid is refined, which causes all discontinuities
to be more sharply resolved: this is also observed in nonlinear ideal MHD simulation
of planar shock interactions with a perturbed interface with an initially transverse
magnetic field [21]. For clearly illustration of the wave structures in the vicinity of
the interface, results in our simulations are obtained on a mesh with 6400 grid cells.
B
B.1
Nonlinear simulation of compressible MHD-RMI
Introduction
The nonlinear simulations were carried out with a nonlinear compressible MHD solver
similar to the method developed by Samtaney [56]. In cylindrical geometry, the ideal
MHD equations expressed in strong conservative form as follows
∂U 1 ∂(rF(U)) 1 ∂G(U)
+
+
= S + S∇·B ,
∂t
r
∂r
r ∂θ
(B.1)
where U ≡ U(r, θ, z, t) = {ρ, ρur , ρuθ , ρuz , Br , Bθ , Bz , e}T is the solution vector, ρ is
the density, (ρur , ρuθ , ρuz ) is the momentum, (Br , Bθ , Bz ) is the magnetic field, e is
the total energy per unit volume. F(U) and G(U) are the fluxes of mass, momentum,
152
(a) Unscaled growth rate of the shocked interface
102
Error
Error at t=0.3
O(∆r)
O(∆r2)
100
10
-2
10
-4
10
-3
-2
10
∆r
(b) Relative error at time t = 0.3
Figure A.1: Time history of unscaled growth rate of the shocked interface for β = 16
wavenumber m = 128 using different mesh sizes: 800, 1600, 3200, 6400, 12800 zones in
the radial domain and (b) Relative error of unscaled growth rate at time t = 0.3.
153
magnetic flux and energy in the r and θ directions given as
F=
ρur
ρu2r
+ p̃ −
Br2
ρur uθ − Br Bθ
ρur uz − Br Bz
0
u r B θ − uθ B r
u r B z − uz B r
(e + p̃)ur − (B · u)Br
,G =
ρuθ
ρur uθ − Br Bθ
ρu2θ + p̃ − Bθ2
ρuθ uz − Bθ Bz
uθ B r − ur B θ
0
uθ B z − uz B θ
(e + p̃)uθ − (B · u)Bθ
,
(B.2)
where p̃ = p + 12 (B · B), is the sum of gas and the magnetic pressures. Equation B.3
presents the source term S arising in cylindrical coordinates, and the stabilization
term S∇·B , it was introduced by Godunov [75]
0
2
Bθ − ρu2θ − p̃
ρu u − B B
r θ
r θ
0
1
S=−
r
0
u B −u B
r θ
θ r
0
0
, S∇·B = −(∇.B)
0
Br
Bθ
Bz
ur
uθ
uz
B.u
.
To close the system we use the perfect gas equation of state so that e =
(B.3)
p
γ−1
+ 21 (B ·
B) + 12 ρ(u · u), where γ is the ratio of specific heats fixed at 5/3. Using the solenoidal
154
constraint on the magnetic field:
∇·B=
⇒
1 ∂(rBr ) 1 ∂(Bθ )
+
r ∂r
r ∂θ
∂Br 1 ∂(Bθ ) Br
+
+
∂r
r ∂θ
r
˜ · B + Br .
∇·B=∇
r
(B.4)
The advecting term can be written as:
˜
+ S̄ ⇒ −(∇.B)
= S∇·B
˜
S∇·B
0
0
Br
Br
B
Bθ
θ
Bz Br Bz
−
r
ur
ur
u
uθ
θ
uz
uz
B.u
B.u
.
(B.5)
Rewriting equation B.1 as:
F
∂U 1 ∂U
∂U
+ S̄ − ,
+ Ar
+ Aθ
= S + S∇·B
˜
∂t
∂r
r
∂θ
r
where the Jacobian singular matrices of the fluxes are Ar =
∂F
∂U
(B.6)
and Aθ =
∂G
.
∂U
to the LHS of equation B.6 we obtain
moving S∇·B
˜
∂U Acθ ∂U
∂U
+ Acr
+
= Sc ,
∂t
∂r
r ∂θ
where Sc = S + S̄ −
F
and Acr , Acθ are nonsingular matrices.
r
(B.7)
155
B.1.1
Numerical Method
In this section, the adopted numerical method is a predictor-corrector unsplit method
in cylindrical coordinates, and is similar to that developed by Samtaney [4], and
Samtaney et al. [56]. Solving equation B.1 by the unsplit method requires nonconservative form for primitive variable W ≡ W (U ). In our implementation we
1
choose W = ρ, ui , Bi , ui , pt T ,where pt = p + Bk Bk is the total pressure. Rewriting
2
the equations using W in quasilinear form, we get for 2-D
∂W Acθ ∂W
∂W
+ Acr
+
= SW ,
∂t
∂r
r ∂θ
(B.8)
where SW = ∇U W Sc . The unsplit algorithm in which face-centered and time-centered
primitive variables are predicted, followed by a corrector step in which a Riemann
problem is solved using the predicted values to compute a second order accurate
estimate of the fluxes and the solenoidal condition on the magnetic field is maintained
by a projection technique. The schematic of the numerical discretization is presented
n+ 1
in Fig. B.1 where i, j is the centroid of each computational cell. We require Wi±,j2
where ± represent the i + 1/2, i − 1/2 (i.e. the right and left faces of cell i, j in the
radial direction).
First in this step, we predict the primitive values at n +
1
2
time step at the faces.
Then, we split the update for primitive variables as follows first for r-direction i.e
n+ 1
i±, j (the same process is applied in θ-direction, i, j±). A Taylor series of Wi±,j2 is:
n+ 1
Wi±,j2 = Wijn +
∆t ∂W n
∆r ∂W n
(
)ij ±
(
) ,
2 ∂t
2 ∂r ij
156
Figure B.1: Schematic of 2-D computational cell.
157
∂W
equation B.8 is used to replace the term
in the previous equation to get:
∂t
n+ 1
Wi±,j2
→
n+ 1
Wi±,j2
=
=
Wijn
Wijn
1
+
2
∆t
+
2
∂W
Aθ ∂W
−Ar
−
+ SW
∂r
r ∂θ
∆t
∂W
±I −
Ar )∆r(
∆r
∂r
n
+
ij
n
ij
±
∆r ∂W n
(
) ,
2 ∂r ij
∆t
∆t Aθ ∂W n
(SW )nij −
(
) . (B.9)
2
2 r ∂θ ij
The predictor step is further divided into a normal and a transverse predictor
steps, which are outlined below.
Normal Predictor: Compute the effect of the normal derivative terms and the
source term on the extrapolation in space and time from cell centers to faces. In this
step we split the primitive variables into everything except the normal field component
(we call it W̄ ) and normal field component as follows
Win =
Then:
W̄i±,j =
W̄in j
1
+
2
W̄in
Bir
∆t
±I −
Ār P± (∆W̄ ),
∆r
(B.10)
(B.11)
where Ār is the matrix obtained from Ar after deleting the row and column corresponding to the normal component of the magnetic field Br ,
r
Bi±
= Bir ,
n
Wi±,j
=
n
W̄i±,j
,
r
Bi±,j
X
P± (W ) =
(lk .W )rk ,
(B.12)
±λ>0
where λk are eigenvalues of Ār , k = 1, 2, .., 7 (for the seven-waves method) and lk
158
and rk are the corresponding left and right eigenvectors. Adding the source term to
Wi±,j , we get
Wi±,j = Wi±,j +
∆t
SW .
2
(B.13)
Stone Correction: In multi-dimensions the derivative of the normal component of
the magnetic field in the θ-direction is not zero. Therefore, a correction originally
attributed to Stone added to the above normal predicted states ([78]) which becomes
W̄i±,j = W̄i±,j −
∆t ∂Bθ
{0, Br , Bθ , Bz , uθ , uz , −(γ − 1)u · B}T .
2∆r ∂θ
(B.14)
Likewise, we compute W̄i,j± and Wi,j± , so we now have Wi±,j and Wi,j± that are used
to compute solution of 1-D Riemann problem (R(.,.)) at i ± 1/2, j and i, j ± 1/2.
The above normal predictor step gives us the left and right states at each cell
interface (Wi+,j , Wi+1−,j ). A seven-wave linearized Riemann solver is applied to obtain
the primitive variables at the cell faces, except the normal component of the magnetic
field, which is taken as the arithmetic mean of the left and right states at i + 1/2 and
j + 1/2. The fluxes are then computed from the solution of the Riemann problems
R(.,.), gives solution at i + 12 , j which is used to compute flux.
F 1D = F (R(Wi+,j , Wi+1−,j ))
(B.15)
G1D = G(R(Wi,j+ , Wi,j+1− )).
Transverse Predictor: Compute final corrections to Wi,j due to the transverse
derivatives,
n+ 1
Wi±,j2 = Wi±,j −
∆t
− G1D
).
∇U W (G1D
i,j− 21
i,j+ 12
2
(B.16)
159
Similarly update
n+ 1
Wi,j±2 .Using
Uijn+1 = Uijn −
the predicted values we calculate the final update:
∆t
∆t
{Fi+ 1 ,j − Fi− 1 ,j } −
{G 1 − Gi,j− 1 } + ∆tSc ,
2
2
2
2∆r
2r∆θ i,j+ 2
(B.17)
where
n+ 1
n+ 1
2
))
Fi+ 1 ,j = F (R(Wi+,j2 , Wi+1−,j
2
Gi,j+ 1 =
2
n+ 1
n+ 12
F (R(Wi,j+2 , Wi,j+1−
)).
(B.18)
Projection: The normal component of the magnetic field at i + 12 , j is used
to compute a cell centered divergence. The magnetic field is projected as Bi+ 1 ,j =
2
Bi+ 1 ,j −∇χ, and replace the corrected magnetic field into
2
n+ 21
,
Wi+1−,j
where the following
Poisson equation is solved using a multigrid technique with a Gauss-Seidel Red-Black
ordering smoother, and a BiCGStab bottom solver
1 ∂
∇ χ=
r ∂r
2
n+ 1
∂Br
∂r
+
1 ∂ 2 Bθ
.
r2 ∂ 2 θ
n+ 1
Then, we recompute the fluxes as Fi+ 12,j = F (Wi+ 12,j ).
2
2
(B.19)
160
C
Papers Submitted and Under Preparation
• A. Bakhsh, S. Gao, R. Samtaney, and V. Wheatley, ˝The cylindrical RichtmyerMeshkov instability]Linear Simulations of the Cylindrical Richtmyer-Meshkov Instability in Magnetohydrodynamics, Physics of Fluids, March, 2016, 28(3), p. 034106.
Also presented in conference +poster.
• A. Bakhsh, and R. Samtaney, ˝Linear Analysis of Converging Richtmyer-Meshkov
Instability in the Presence of an Azimuthal Magnetic Field, Journal of Fluids Engineering, November, 2017, 140(5),p. 050901. Also presented in conference +poster.
Abeer Bakhsh