1. a7

3.
H1

:nnmﬁu.
. II}!!! '7
3...... 3L

. .. n.»

 

a. A.
w». 1:94..

 

, .

Hut-HIP. : 2.13!

i. .i. . D sun
2....

.

r
at!

imam

if!
L

.. :3! .
ti); in»

w t... .114 I .Ili V”

5... {ii-ﬁt Ithnuwu.“

have 4 Km.

:rs a. .
I... ll
7}

.35... . .Q?
MK...‘ :7»: . 3 .41". I". .
.2 . $.31; 3...“. ..

 

.3... he v.33 «.3 3?. 313....

WES” LIBRARY
N? Michigan State
I University

This is to certify that the
dissertation entitled

LATTICE BOLTZMANN SIMULATION OF LASER
INTERACTION WITH WEAKLY IONIZED PLASMAS

presented by
Huayu Li

has been accepted towards fulﬁllment
of the requirements for the

PhD. degree in Mechanical EmineerinL

 

 

ﬂ [90157011

Majo’r Profé’ssor’s Signature
/’0//? /,2006D

Date

 

MSU is an aﬁnnative-actrbn, equal-opportunity employer

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

5/08 K.IProj/Acca-PrelelRCIDateDue indd

 

LATTICE BOLTZMANN SIMULATION OF LASER INTERACTION WITH
WEAKLY IONIZED PLASMAS

By

Huayu Li

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Mechanical Engineering

2008

ABSTRACT

LATTICE BOLTZMANN SIMULATION OF LASER INTERACTION WITH
WEAKLY IONIZED PLASMAS

By

Huayu Li

Laser-plasma interaction (LPI) is an important subject to a variety of disciplines
in engineering and science, such as laser welding, pulsed laser deposition (PLD), laser—
generated x-rays and laser-aided ignition of inertial confinement fusion (ICF). In
particular, laser interaction with weakly ionized plasmas has invoked a great deal of
interest to the laser manufacturing community because plasmas naturally appear and
interact with a laser beam in such high energy manufacturing processes. Due to the
complexity and richness of physics, numerical model studies have been pivotal in the
understanding of LPI. A number of numerical models have been created to study LPI and
help design LPI equipment, and there are basically two kinds of numerical models: the
kinetic-based model and the hydrodynamic model. Although kinetic models (e.g.,
particle-in-cell model) have been very successful, they are computationally expensive in
most cases and their application is rather limited. Hydrodynamic models are also a
powerful tool for LPI simulations, but they fail in some circumstances because they are
based on the continuum assumption.

In this study, a new numerical model based on the lattice Boltzmann method

(LBM) is introduced to simulate laser interaction with weakly-ionized plasmas. The LBM

is a kinetic theory based method, where the distribution functions of the individual
species of particles in the plasma are solved and thus the macroscopic variables (such as
number density and momentum) are obtained. In this study, the Boltzmann equation with
ionization and recombination terms is solved. Since only number density and momentum
can be correctly retrieved from the two—dimensional nine-bit (D2Q9) discretization
scheme, a set of energy equations is derived from the Boltzmann equation and solved
separately to calculate temperature fields. The electromagnetic ﬁeld from both laser and
plasma is updated by solving Maxwell’s equations using the ﬁnite-difference time-
domain (FDTD) method. In the implementation of the present model, a rescaling scheme
is introduced to select the appropriate simulation parameters for the LBM, so that the
physical properties of the plasma can be used. This rescaling scheme has been validated
by hydrodynamic flow problems and the electron diffusion problem. In this study, a two-
dimensional weakly-ionized helium plasma interaction with a continuous wave C0; laser

beam is simulated.

This model is a mesoscale approach based on the kinetic theory and the LBM, so
it has a number of inherent advantages over previous models. Because the LBM solver is
employed, this approach is computationally efﬁcient and easy to parallelize. In addition,
this model is capable of predicting time-dependent number densities, velocities, and
temperatures of all particle species for a fairly large scale problem without employing the
continuum assumption. It is believed that this model has a lot of potential for the studies

of weakly ionized plasmas in a wide spectrum of applications.

DEDICATED TO MY WIFE AND MY PARENTS

ACKNOWLEDGEMENT

I would like to deeply acknowledge my advisor, committee members, colleagues
and families, because this thesis work could not have been done without their support
and/or contribution.

First of all, I wish to express my special appreciation to my advisor, Professor
Hyungson Ki, whose guidance and directions have been greatly supportive throughout
my entire Ph.D. studies. In the summer of 2004, I joined Dr. Ki’s team and started my
research on modeling laser interactions with matters by using the lattice Boltzmann
method under his supervision. He spent enormous time and energy to help me to start as
a self-motivated graduate student. I am especially grateful not only for his patience and
encouragement on the academic side, but also for his generous supports in my personal
life.

I also would like to thank my committee members, Dr. Andre Benard, Dr. Patrick
Kwon and Dr. Timothy Grotjohn, whose advice and comments were greatly valuable and
helpful for my research and dissertation.

Working with my previous and current colleagues and ﬁ'iends, Kuiyan Xu, Mehdi
Abarham, Taylor Young, Brian Powell, Ankit Aggarwal and Admir ijanin, makes the
past ﬁve-year unforgettable. Their friendship, help and encouragement will remain in my
memory.

Finally, I would thank my families, my wife Lei Jin, my parents and my in-laws.
They have spent extensive time and efforts to support me. Their immeasurable loves and

understandings are overwhelming and go beyond words.

TABLE OF CONTENTS

LIST OF TABLES ........................................................................................................... viii

LIST OF FIGURES ........................................................................................................... ix

NOMENCLATURE ........................................................................................................ xiv

CHAPTER 1

BRIEF REVIEW OF LASER-PLASMA INTERACTIONS ............................................. 1
1.1 Physics oflaser-plasma interactions ................................................................. I
1.2 Applications oflaser-plasma interaction .......................................................... 6
1.3 Numerical Models ............................................................................................. 9
1.4 Main Content of This Research ...................................................................... 23

CHAPTER 2

LATTICE BOLTZMANN METHOD FOR INCOMPRESSIBLE HYDRODYNAMICS

FLOWS ............................................................................................................................. 26
2.1 Boltzmann Transport Equation and Its Collision Term .................................. 26
2.2 From Boltzmann Equation to Lattice Boltzmann Equation ............................ 31
2.3 From Lattice Botzamnn Equation to Navier-Stokes Equation ....................... 36
2.4 LBM Simulation of Two-dimensional Driven Cavity Flow ........................... 40

CHAPTER 3

EVALUATION OF THE EXTERNAL FORCE TERM IN LATTICE BOLTZMANN

METHOD ......................................................................................................................... 43
3.1 Different External Force Models .................................................................... 43
3.2 LBM Simulation of Magnetohydrodynamics Flows ...................................... 56

CHAPTER 4

LATTICE BOLTZMANN SIMULATION OF ISOTHERMAL WEAKLY IONIZED

PLASMAS ........................................................................................................................ 73
4.1 Mathematical Model ....................................................................................... 74
4.2 Rescaling Scheme ........................................................................................... 80
4.3 Simulation Results .......................................................................................... 98

CHAPTER 5 .

LATTICE BOLTZMANN SIMULATION OF WEAKLY IONIZED HELIUM

PLASMA WITH IONIZATION AND RECOMBINATION ........................................ 108
5.1 Mathematical Model ..................................................................................... 108
5.2 F inite—Difference Time-Domain Method ...................................................... 115
5.3 Numerical Implementation ........................................................................... 122
5.4 Simulation Results ........................................................................................ I31

vi

CHAPTER6
CONCLUSKMV ............................................................................................................... 150

REFERENCE .................................................................................................................. 154

vﬁ

LIST OF TABLES

Table 1 Fa and u* in different external force models .................................................... 47
Table 2 Relative errors (X104) ofPoiseuille flow ............................................................. 49
Table 3 Relative errors (X105) of free diffusion problem ................................................. 54
Table 4 Quantitative comparison of simulation results with Dellar’s work ..................... 67

Table 5 Simulation parameters and simulation errors for different values of 2' in the
simulation of two—dimensional driven cavity ﬂows .......................................................... 91

Table 6 Two examples of the air state that give invalid dimensionless relaxation time .. 94

Table 7 Rescaled simulation parameters and # of iterations until convergence ............... 94
Table 8 Simulation parameters for air at 300 K and 1 atm (Case 1) ................................ 95
Table 9 Simulation parameters used for the electron diffusion problem .......................... 99
Table 10 Macroscopic Quantities from Moments of Distribution Function ................... l 13
Table 1 1 Primary simulation variables ........................................................................... 126
Table 12 Secondary simulation variables ....................................................................... 127

viii

LIST OF FIGURES

Figure 1 Schematic of D2Q9 square lattice. The two-dimensional continuous phase space
is discretized by nine microscopic velocity components in a square lattice. .................... 34

Figure 2 Streamlines of driven cavity flow with different Reynolds numbers. The flow
pattern is strongly dependent on the Reynolds number. If the Reynolds number is
sufﬁciently large, secondary vortices appear at the right and left lower corners. ............ 41

Figure 3 (3) Velocity proﬁles for u through the geometric center of the cavity; (b)
Velocity proﬁles for v through the geometric center of the cavity ................................... 42

Figure 4 Proﬁle of x-component velocity obtained by He’s model (g = 5.774X10'2). The
simulation results agree very well with the analytical solution ........................................ 50

Figure 5 Evolution of relative errors of simulation of Taylor vortex ﬂow obtained by the
ﬁve external force models. All the errors are in the same order of values. ...................... 51

Figure 6 Proﬁles of x-component velocity along the center line of x-axis at time steps of
1500 and 3000 by using He’s model. The simulation results are in good agreement with
the analytical solutions ...................................................................................................... 52

Figure 7 Fluid density distribution obtained by He’s model and Guo’s model at 200 time

steps with different magnitudes of the external force: (a) 3x = 4.0 (b) aX = 10.0. Both
Guo‘s model and He’s model can produce good results if the external force is small (as
shown in (a)). However, if the external force is relatively large, Guo’s model induces
larger error than dose He’s model (as shown in (b)). ....................................................... 55

Figure 8 Proﬁles of x-component velocity with different Hartmann numbers: M=O
(squares), M=1 (circles), M=2 (upper triangles), M=5 (lower triangles), M=10
(diamonds), M=20 (stars). The solid lines show the analytical solutions. The simulation
results agree very well with the analytical solutions. ....................................................... 62

Figure 9 Proﬁles of x-component magnetic ﬁeld for different Hartmann numbers: M=0
(squares), M=1 (circles), M=2 (upper triangles), M=5 (lower triangles), M=10
(diamonds), M=20 (stars). The solid lines show the analytical solutions. The simulation
results agree very well with the analytical solutions ........................................................ 62

Figure 10 Evolution of vorticity for the 2-D Orszag-Tang vortex. The flow pattern
becomes complicated due to the nonlinear interactions between the velocity ﬁeld and the
magnetic ﬁeld. A flat quadrupole-like conﬁguration emerges with time lapse ................ 64

ix

Figure 11 Evolution of current density for the 2-D Orszag-Tang vortex. The existing
current sheet at the center of the ﬁgure is enhanced and eventually, a thin elliptic
structure establishes due to the magnetic reconnection occurring there. .......................... 65

Figure 12 Iso-surface contours of magnitudes of vorticity and current density at t=0 and
t=0.598 sec for the 3-D Orszag-Tang vortex. The patterns of the current sheet and the
corresponding vorticity resemble Figure 11 and Figure 10 closely because initially, the
variations of current density and vorticity in z-directions are relatively smaller in values
than those in the x-and y-direction .................................................................................... 69

Figure 13 Evolution of magnetic flux function for doubly periodic coalescence instability.
The position of the magnetic current sheets does not change with time. The magnetic
islands move towards each other. The two comers coalesce into one and the original two
square cells become two adjacent pentagons. A current sheet forms between the two cells
and the intensity of the current sheet increases. Eventually, the neighboring square cells
merge together, simplifying the topology structure of the magnetic ﬁeld to four square-
like islands. ....................................................................................................................... 71

Figure 14 Maximum value of current density vs. magnetic diffusivity. The temporal

maximum of current density can be approximated as jmax oc 77'1/2. .............................. 72

Figure 15 Schematic of second order interpolation method for ion (neutral) distribution
function. The value of distribution function of the heavy particles at p can be found by
using the post-collision values at 0 ', p’ and q’. ................................................................ 79

Figure 16 Distribution of electron number density at t = 3.5 ns. Non-physical peaks
appear if the physical variables of electrons are used ....................................................... 84

Figure 17 Electron number density distribution at (a) t = 1.116 ns and (b) t = 3.347 ns
with different ionization degrees. Non-physical sub-peaks caused by using the physical
variables travel with the constant lattice speed. A higher ionization degree leads to a
higher lattice speed. The ratio of the magnitude of the sub-peaks to that of the primary
peak corresponds to the value ofthe ratio ofthe weight coefﬁcients ............................... 85

Figure 18 Contour plot of electron number density at t = 1.116 ns. The distribution of
electrons is along with the discretized velocity directions if the effects of collisions are
not included in LBM. ........................................................................................................ 86

Figure 19 Velocity proﬁles along the central lines of the computational domain obtained
by using the physical properties of air without the rescaling when the dimensionless
relaxation time is in the valid range. (a) dimensionless x~velocity (b) dimensionless y-
velocity .............................................................................................................................. 90

Figure 20 Number ofiterations until convergence vs. I (Re = 200). A larger r leads to a
quicker convergence ......................................................................................................... 92

Figure 21 (a) r as a function of Re and Ma (N=128). It is apparent that 1 needs to be
close to 0.5 if a high Ma and/or a high Re are desired. (b) Re vs. I at different grid
densities. Increasing grid density is an effective way to increase Re if the other simulation
parameters are kept unchanged. ........................................................................................ 93

Figure 22 Locations of the primary vortices at different Reynolds numbers. The
simulation results obtained by the rescaling scheme agree well with other researchers’
results. ............................................................................................................................... 95

Figure 23 Velocity proﬁles along the central lines of the computational domain
(a) dimensionless x-velocity (b) dimensionless y-velocity. The simulation errors are
relatively large if the Reynolds number is high. It is due to the high Mach number
resulted by the high Reynolds number according to Eq. (119) ......................................... 97

Figure 24 Electron number density distributions with different initial number densities of
neutrals. The simulation results obtained by the rescaling scheme are in good agreement
with the analytical solutions ............................................................................................ 100

Figure 25 Electron number density distribution (Circle: without internal electric ﬁeld,
Squares: with internal electric ﬁeld, both at 0.223 ns. Solid line shows the initial
distribution). The internal electric ﬁeld generated by the charge density largely enhances
the electron diffusion process. ........................................................................................ 102

Figure 26 Snapshot of electron number density (solid lines) and ion number density
(dashed lines) at t = 1.116 ns. The electrons and ions move along opposite directions
because they have opposite charges. The ions travel much more slowly than the electrons
since the ions have much heavier mass ........................................................................... 103

Figure 27 Electron number density distribution (at t = 6.69 ns) obtained by using different
equilibrium distribution functions for the external force term. Apparent simulation errors
would be resulted by using the cross-collision equilibrium distribution function in He’s
external force model. ...................................................................................................... 104

Figure 28 Relative errors in electron number density vs. time for three different grids.
The errors are calculated at the center of the domain with comparison to the analytical
solution. Second order accuracy in space discretization can be seen from the ﬁgure. 105

Figure 29 Time evolution of number density and x-component velocity of electrons at the
center of the computational domain. (Solid lines: simulation results obtained without the
rescaling scheme, circles and squares: simulation results obtained with the rescaling
scheme) ........................................................................................................................... 106

Figure 30 Position of the electric and magnetic ﬁeld vector components about a cubic
unit cell ofthe Yee space lattice ..................................................................................... 1 17

Figure 31 Schematic of computational domain for simulation of a femtosecond laser
pulse interaction with a silicon substrate. ....................................................................... 1 18

xi

Figure 32 Time history of laser pulse propagation. (In these plots, the radial component
ofthe electric ﬁeld (E) is used.) ..................................................................................... 120

Figure 33 Time history of laser energy absorbed by electrons ....................................... 121

Figure 34 Schematic of the rescaling scheme. The value of distribution function at p can
be found by using the post-collision values at 0 ', p ' and q’. .......................................... 125

Figure 35 Positions of the simulation variables. HZ, n3, uS and TS are calculated at the

white node points while Ex and E), are staggered in space and evaluated at the black node

points ............................................................................................................................... 126
Figure 36 Flowchart of the present model ...................................................................... 130
Figure 37 Correlation among three equations ................................................................. 131

Figure 38 Schematic of the computational domain for simulation of interaction between a
continuous C02 laser beam and a weakly ionized helium plasma .................................. 132

Figure 39 Evolution of maximum electron number density with different laser intensities.
At the early stage, the electron number density increases exponentially because of the
laser-induced ionization. However, it shows saturation later due to the equilibrium
between ionization and recombination. .......................................................................... 133

Figure 40 Evolution of maximum electron temperature with different laser intensities.
The electron temperature increases very quickly due to the rapid heating. It decreases
after the ionization takes place because the generation of new free electrons consumes the
laser energy. .................................................................................................................... 134

Figure 41 Evolution of ionization coefﬁcient with different laser intensities. The
ionization coefﬁcient shows a similar pattern with the electron temperature because it
mainly dependent on the electron energy. ...................................................................... 135

Figure 42 Evolution of recombination coefﬁcient under different laser intensities. The
recombination coefﬁcient increases due to the generation of new free electrons and the
decrease of the electron temperature ............................................................................... 136

Figure 43 Evolution of maximum ion temperature with different laser intensities.
Compared to the electron temperature, the ion temperature increases with a lower rate
because of the heavy mass. It decreases due to the attenuation of the electric ﬁeld. ..... 137

Figure 44 Evolution of neutral temperature with different laser intensities. Unlike the
temperatures of electrons and ions, the neutral temperature keeps increasing because the
neutrals are heated by the cross-heat transfer from the charged particles. ..................... 138

Figure 45 Distribution of electron temperature along beam propagation axis at different
time moments. The highest temperature always appears at the interface. ...................... 139

xii

Figure 46 Distribution of electron number density along the beam propagation axis at

different time moments. Most of the new free electrons are generated near the interface.
......................................................................................................................................... 139

Figure 47 Snapshots of electron and neutral number densities at (a) t = 49.93 ps and (b) t
= 149.80 ps ...................................................................................................................... 140

Figure 48 Contours of electron temperature at (a) t = 4.86 ps, (b) t = 13.26 ps, (0) t =
49.93 ps and (CI) 14980 ps .............................................................................................. 142

Figure 49 Electron number density distribution along x-axis atj = 195 at different time
moments .......................................................................................................................... 143

Figure 50 Snapshot of (a) impact ionization coefﬁcient and (b) recombination coefﬁcient
at t = 149.80 ps ................................................................................................................ 143

Figure 51 Snapshots of Joule heating at (a) t = 4.86 ps, (b) t = 13.26 ps, (c) t = 49.93 ps
and (d) t = 149.80 ps ....................................................................................................... 145

Figure 52 Snapshots of x-component of electric ﬁeld at (a) t = 4.86 ps, (b) t = 13.26 ps, (0)
t = 49.93 ps and (d) t = 149.80 ps ................................................................................... 146

Figure 53 Evolution of (a) number densities (b) temperatures (0) ionization and
recombination coefﬁcients and (d) magnitude of electric ﬁeld at Point A ..................... 148

xiii

()

tr,

f5
f val

[5:7

NOMENCLATURE

Acceleration of charged particles caused by the external force

x-component of the acceleration

y-component of the acceleration

Magnetic induction (magnetic ﬂux density)

x-component of magnetic induction

y-component ofmagnetic induction

Lattice speed of particle species 5

Diffusivity
Electric ﬂux density
Electric ﬁeld

x-component ofthe electric ﬁeld

y-component ofthe electric ﬁeld

Unit electronic charge

Discretizcd microscopic velocity ofparticle species 5 along the a -th
direction

Single particle distribution function ofparticle species 5

Self-collision equilibrium distribution function of particle species S

Cross-collision equilibrium distribution function of particle species 3
due to the binary collisions with particle species 3'

xiv

"I

II

([5
q .s-
R 5.
T

S

7:41)):

Discretizcd external force term along the a -th direction

Pressure gradient
Magnetic ﬁeld

z-component ofthe magnetic ﬁeld

Electric current density

Boltzmann constant

Characteristic length of the flow problem
Mass ofparticle species 3

Number density of particle species 5
Number of grids along x-direction

Number of grids along y-direction

Kinetic pressure dyad ofparticle species 3

Electronic charge ofparticle species 5

Heat flux of particle species 5

Ionization and recombination coefﬁcient of particle species 5
Temperature of particle species 5

Post-collision temperature of electrons (ions) due to binary collisions
with neutrals

Macroscopic velocity of particle species 5

Barycentric velocity ofelectrons (ions) due to binary collisions with
neutrals

Equilibrium macroscopic velocity ofparticle species 3

XV

VI!

A

‘90

5 T5

75

.9:

Ps

p,.

Characteristic velocity ofthe flow

First ionization potential of helium
Microscopic velocity of particles species 5
Drift velocity

Interpolation parameter of particle species 5
Electrical permittivity in free space

Total energy ofparticle species 5

Themral energy of particle species 3
Rescaling parameter of particles species 3

Magnetic diffusivity

Thennal conductivity of particles species 5

Relaxation time ofparticle species 3 due to the binary collisions with
particle species 5'

Magnetic permeability
Fluid viscosity
Lattice viscosity

Sound speed of particle species 5

Lattice sound speed ofparticle species 5

Mass density of particle species 3
Charge density

Electrical conductivity

xvi

5)!

Cross section ofbinary collisions between particle species 3 and n
Cross section of electron impact ionization
Dimensionless relaxation time of particle species 5

Time step
Space step along x-direction

Space step along y-direction

xvii

Chapter 1

Brief Review of Laser-Plasma Interactions

Laser-plasma interactions are of great interests to researchers from many different
ﬁelds of expertise. The characteristics of laser—plasma interactions are high nonlinearity,
strong coupling, large spans of spatial and temporal scales, and multiple physics. A
complete description oflaser-plasma interactions should include laser-induced ionization,
absorption of laser energy. plasma responses and electromagnetic wave propagation
inside the plasma, all of which are strongly coupled and dependent upon the incident
laser intensity. The higher the laser intensity, the more complicated the physics involved
will be. In this chapter, the physics, applications and numerical models of laser-plasma
interactions will be brieﬂy reviewed. The main content of this dissertation is also

introduced.

1.] Physics of laser-plasma interactions

When a laser beam impinges on a target which is initially in non-plasma state, the
neutral atoms or molecules will be excited and free electrons will be generated ifthe laser
irradiance (or intensity) exceeds some threshold. For metals. in which there are already
“free" electrons in the conduction band, the electron impact ionization caused by the
thermal collisions between electrons and heavy particles (ions and neutral particles)
dominates the ionization process. The impact ionization mainly depends on the electron
energy that is obtained from the incident laser beam through inverse bremsstrahlung

heating. For dielectrics. there is a well-deﬁned threshold for the laser-induced breakdown

which depends on the laser pulse duration [1, 2], the wavelength of the incident laser
pulse [3] and the band gap or ionization potential of the material. Besides the impact
ionization. the multi—photon ionization also plays an important role in laser-induced
breakdown of dielectric. In the multi-photon ionization, the laser energy carried by the
photons shifts the electrons from the valence band to the conduction band and thus
produces the “feed electrons” for the following impact ionization. Ifthe laser intensity is
high enough, the strong electric ﬁeld ofthe laser beam will distort the potential barrier of
the atom (or molecule) drastically. The electrons can pass the ban‘ier easily to become
free electrons. This ionization mechanism is referred as to tunneling ionization or ﬁeld
ionization. Laser-induced ionization of gases has similar mechanisms as the dielectrics.
Parameters of the incident laser beam, ionization potential of the gas and the gas state
(such as temperature and pressure) are the main impact factors of the ionization process.
Besides the ionization of neutral particle (atoms. molecules, etc.), other atomic physics,
such as excitation of neutral particles or ions, three-body recombination of charged
particles, charge attachment and de-excitation, also happen simultaneously inside the
plasma. Three-body recombination is the inverse process of ionization and lead to losses
of free electrons. It depends on the number density of charged particles and the electrons
temperature. Thus it is only prominent in dense, low-temperature plasmas.

Once the laser-induced plasma emerges, the laser energy will be absorbed into the
plasma through different mechanisms. At relatively low laser intensity, the inverse
bremsstrahlung heating [4, 5] (also referred to as collisional heating), absorption of
photons by charged particles undergoing collisions. dominates the absorption of laser

energy. In dense plasma. where the plasma frequency is higher than the laser ﬁeld

frequency, the stimulated Compton effects [6] prevails the inverse bremsstrahlung for
heating ofplasma electrons [7]. In a nonunifomt plasma, the ponderomotive effect due to
the electron acceleration in the nonuniform electric ﬁeld from the incident laser beam
may be efﬁcient for plasma heating if the nonuniforrnity of the plasma is sufﬁciently
strong so that the characteristic gradient scale length is shorter than the electron
oscillation amplitude [7]. At high laser intensities (above 10'5 V’Wcm2 or so [8]), the
collisions among particles become ineffective [9] and the quiver velocity of electron is
comparable to the thermal velocity, which reduces the effective collision frequency
further [10]. Instead of collisional heating, some collisionless absorption mechanisms
become prominent in this regime of laser intensities. These collisionless mechanisms are
based upon the parametric instabilities that are caused by the plasma parameters
oscillation with the strong laser ﬁled: it is possible for the parameters to become resonant
with the ﬁled ﬂuctuations just as the case of mechanical vibration systems. As the
representative collisionless mechanisms for laser energy absorption, two-plasmon photon
decay, stimulated Raman scattering, stimulated Brillouin scattering and plasma resonance
are strongly depend upon the plasma nonuniformity. The two-plasmon decay [11, 12]
involves the decay ofa laser photon into the quantum ofa Langmuir plasma wave and an
ion acoustic wave. This mechanism is remarkable in sufﬁciently dense plasmas where the

plasma density 11> n(.,./4 where rim. is the critical electron density of the plasma. The

stimulated Raman scattering (SRS) [13, 14] usually happens in underdense plasma where

n < nor/4. ln SRS, the laser energy is absorbed by the plasma Langmuir wave and the

transverse electromagnetic ﬁeld which is shifted in frequency. The stimulated Brillouin

scattering (88$) [15, 16], which is signiﬁcant for longer laser pulses, absorbs the laser

energy by exciting an ion acoustic wave. If the nonuniformity of the plasma is
sufﬁciently strong, the plasma resonance absorption presents and excites a plasma wave.
This plasma wave grows over several laser periods and is eventually damped either by
collisions at low intensity area or by particle trapping and wave breaking at high intensity
area [15]. If the laser intensity increases further to relativistic regime (usually higher than
1018 Wr’cmz), two more parametric instabilities emerge: laser ﬁlamentation [17, 18] and
laser modulation [19. 20], both of which are closely related with the electron mass
oscillation due to the relativistic effect. Another heating mechanism of electrons exists
near the interface between vacuum and plasma, which is called vacuum heating [21, 22].
In this scenario. electrons are dragged away from the plasma surface into the vacuum,
turned around and accelerated back into the solid all within halfa laser cycle and heated
due to the strongly inelastic inverse bremsstrahlung or the ponderornotive scattering [7].
Vacuum heating occurs when a p-polarized laser pulse is obliquely incident on the
vacLiam—plasma surface while another similar heating mechanism, jx B heating [23, 24],
happens when the laser pulse impinges the plasma surface normally. In jx B heating, the
electrons oscillating near the interface of vacuum-plasma are forced into the plasma by
the magnetic ﬁeld ofthe laser beam every halfperiod ofthe laser.

Besides absorptions of laser pulse energy, the laser pulse propagation inside the
plasma is also an important aspect in laser-plasma interactions. A basic phenomenon in
this category is the ionization-induced defocusing [25, 26] of the laser beam. When a
laser pulse whose center intensity is close or above the ionization potential propagates
inside the plasma, more free electrons due to the ionization will appear around the center

of the beam, resulting in a steep radial gradient of electron number density. Since the

refractive index ofthe plasma is dependent on the distribution ofelectron number density,
this radial gradient of electron number density will lead to a smaller refractive index on
the propagation axis, acting as a defocusing lens for the following laser beam. However,
when the laser intensity is in the relativistic range, the self-focusing mechanism of the
laser pulse propagation emerges due to the interplay between relativistic modiﬁcation of
the electron mass and the ponderomotive effect [27]. By considering the relativistic effect,
the distribution of refractive index is not only dependent upon the electron number
density, but also reliant on the relativistic factor that can be modiﬁed by the electron mass.
Meanwhile. the electrons at the center position are expelled by the laser ponderomotive
force, leading to a reduction ofthe electron number density on axis. Owing to above two
effects. the refractive index of the plasma decreases on axis, acting as a focusing lens for
the following portion of the laser beam. Based upon the dynamical equilibrium between
the defocusing and the self-focusing [28], a plasma channel can be established in the
plasma, through which the laser pulse can propagate a relatively long distance which is
signiﬁcant in many applications.

Generation of fast electrons is a representative physics of plasma responses to the
laser-plasma interaction. In the interaction between underdense plasmas and short-pulsed
high-intensity lasers, electrons can be heated to a very high temperature during the laser
pulse period while the ions, due to their large inertia, remain “cold” at the same time.
Thus. a large local charge separation is resulted and the electrons oscillate at the plasma
frequency, creating alternating regions ofnet positive and negative charge. Then a plasma
wave wake is generated whose phase velocity is roughly equal to the velocity ofthe laser

pulse, which is approximately equal to the speed of light: in the tenuous plasmas. The

plasma wave wake has a longitudinal electric ﬁeld so that it can efﬁciently accelerate the
charged particles [29], potentially to GeV level [8]. Those highly energetic electrons can
lead to a large space charge potential and ions can be accelerated in this electrostatic
sheath to MeV level under short-pulse high-intensity laser-plasma interaction [29].
Another direct result of those hot electrons is generation of hard X-rays. A fast electron
can penetrate into the cold target due to its long mean free path where it either emits
bremsstrahlung via collisions with ions, or produces line radiation by knocking out a

bound K-shell electron [8].

1.2 Applications of laser-plasma interaction

Thanks to the robust physics involved in the laser-plasma interaction, there are
many applications of laser-induced plasmas and laser-plasma interactions.

In the ﬁeld of laser material processing. laser-induced breakdown is considered to
be the main reason for the high ablation precision of metals and dielectrics by ultrashort
lasers [2, 30]. In the application of pulsed laser deposition (PLD) [31, 32], the
interactions of incident laser pulses and the resulted plasma plume, such as absorption of
the laser energy within the expanding plume, expansion of the plasma plume and
evolution of the plasma temperature, have signiﬁcant effects on the uniformity and
quality of the deposited ﬁlm. In laser welding, the incident laser beam will ionize the
target material so that partially ionized plasma is generated in the keyhole and near the
top of the workpiece surface. The plasma inside the keyhole, as a good conductive
medium. will assist the energy transport from the laser beam to the workpiece [33]. On

the other hand, the inhomogeneous distribution of the free electrons in the plasma on the

top surface ofthe workpiece will affect the propagation ofthe laser beam, defocusing the
beam and leading to undesirable laser energy dissipation.

Laser-plasma interaction is even more important in the context of inertial
conﬁnement fusion. In the indirect-drive approach [14], bunch of laser beams heat the
inner wall of a gold cavity called holhraum that contains the deuteriurn-tritiurn (DT)
pellet. creating a superhot plasma which radiates a uniform “bath” ofsoft X-rays. The X-
rays rapidly heat the outer surface ofthe fuel pellet. causing a high-speed ablation ofthe
surface material and imploding the fuel capsule. Symmetrically compressing the capsule
with radiation forms a central hot spot where fusion process could take place. Comparing
to the indirect-drive approach, the direct-drive approach [34] uses powerful beams of
laser to directly focus on the pellet. The rapid heating caused by the laser “driver” makes
the outer layer ofthe target explode. The remaining portion ofthe target is driven inwards
in a rocket—like implosion, causing compression of the fuel inside the capsule and the
formation ofa shock wave, which further heats the fuel in the very center and results in a
self-sustaining burn. A new ignition scheme, fast ignition [35], has gained more and more
interests recently. The main feature of fast ignition is that the compression stage is
separated from the ignition phase. Fast ignition uses the same hardware as the direct-
drive approach but adds a high-intensity, ultrafast pulsed laser as the “spark” that
achieves ignition. The DT target is ﬁrst compressed to high density by laser beams, and
then the short-pulsed laser delivers energy to ignite the compressed core. An advantage of
fast ignition is that the density and pressure of the DT target are less than those in
indirect- or direct-drive approaches. In addition, much less mass ofthe DT fuel is used in

the fast ignition, resulting in reduced input energy. From the above description of the

three schemes of ignition for ICF, one can conclude that the laser-plasma interaction
plays a very important role and deep understanding ofthe laser-plasma interaction in the
context oflCF is highly desired.

Based upon the physics of laser-generated fast electrons and ions as described
above, the concepts of laser wakeﬁeld accelerator (LWFA) [36] and self-modulated laser
wakeﬁeld accelerator (SMLWFA) [20, 37] have been studied for the potential realization
of table—top particle accelerator. In LWFA, an electron plasma wave is driven resonantly
by a short laser pulse whose pulse length is approximately equal to the plasma wave
length through the laser ponderomotive force [38, 39]. The electron energy up to 200
MeV has been reported by Malka‘s group at LOA [40] by using LWFA. SMLWFA is a
hybrid scheme which combines the elements of forward stimulated Raman scattering and
the concept of LWFA, in which an electromagnetic wave decays into one plasma wave
and another forward propagating light wave via the stimulated Raman forward scattering
instability. The most impressive results by using this scheme are reported by a group
working with the Vulcan laser at RAL, UK [40, 41] with observations of electrons at
energies as high as 120MeV.

Laser-plasma interaction can also be used for atomic emission spectroscopy. For
example, laser-induced breakdown spectroscopy (LIBS) [42] utilizes a highly energetic
laser pulse as the excitation source. By focusing the laser onto a small area at the surface
of the specimen. a very small amount of material is ablated. The resulted plasma plume
expands and cools within a very small timeframe. At this point the characteristic atomic
emission lines of the elements can be observed. The main feature of LIBS is that it can

analyze any matter regardless of its physical state, be it solid, liquid or gas, limited only

by the power of the laser as well as the sensitivity and wavelength range of the
spectrograph and detector. Another advantage of LIBS is that it is considered essentially
a non-destructive or minimally-destructive technique because only a small amount of
material is ablated during the LIBS process and there is almost no specimen heating
surromrding the ablation site. In addition, since LIBS is purely an optical technique, the
remote analysis can be achieved easily, which is especially signiﬁcant for use in
hazardous emt'ironments. The drawback of LIBS is that it is subject to variation in the
laser spark and resultant plasma which often limits reproducibility. In order to enhance
the signal from the plasma emission, double-pulse LIBS [43] was developed. Depending
on pulse separation, the second pulse is more or less absorbed by the plasma plume
caused by the previous pulse, resulting in a reheating of the laser plasma and leading to
signal enhancement.

Another application of laser-plasma interaction can be seen in the ﬁeld ofmedical
imaging, which is mainly dependent upon the laser-plasma-generated X-rays. The
advantages of using laser-plasma X-ray source include the possibility of substantial dose
reduction [44], enhancement of image contrast and reduction ofunwanted information by
using differential imaging with rapid simultaneous exposure [45], and the possibility to

place the X-ray source inside the object ofinvestigation [46].

1.3 Numerical Models
Besides experimental studies and theoretical analysis. numerical modeling is also
a powerful tool to obtain deep understanding of the laser—plasma interaction physics and

to provide guidelines and optimism methods to applications of laser-plasma interactions.

9

There are mainly two categories of numerical methods for modeling laser—plasma
interactions: hydrodynamics models and kinetic based models.
1.3.1 Hydrodynamics models

The common feature of the hydrodynamics models is that the individual particle
species of the plasma or the plasma itself is treated as a continuum ﬂuid. The
conservations of mass. momentum and energy are described by hydrodynamic equations.
The speciﬁc formulations of the conservation equations, generation of free electrons and
laser energy deposition vary in different models.

Hydrodynamics simulation of laser-plasma interaction are widely applied in the
context of laser-driven ICF, providing guidance for the design of hohlraum targets on the
ignition facility and detailed understanding of the interaction physics between incident
laser beam and the plasma. The very ﬁrst hydrodynamic model for intense laser
interaction with underdense plasma is LASNEX [47] that was developed by the
Lawrence Livermore Laboratory in the late of 1980’s. LASNEX solves hydrodynamic
equations in Lagrangian coordinates which include flux—dependent electron and ion
thermal conductivities. radiative transport for the plasma radiation, laser-radiation
absorption and some other physical processes [48]. In the past 30 years, many physical
packages such as magnetic ﬁeld and ray trace have been added into the code of LASNEX,
making it a powerful simulation tool in many applications [49-51]. Another early stage
hydrodynamics model is SAGE [52] which was developed for simulation ofthemral self-
focusing in laser plasmas. SAGE is a two—dimensional Eulerian code including a self-
consistent ray-tracing algorithm for laser beam propagation. Thus, the dependence of

self-focusing on the laser wavelength and the laser energy absorption can be simulated

more realistically. This model uses perfect-gas equation of state and a one-temperature
fluid assumption. Laser beam diffraction, radiation emission, superthermal-electron
generation and magnetic effects are neglected. However, the above assumptions can be
well justiﬁed ifonly low-atomic-number targets are considered [52].

In the past decade, motivated by the fast development of the new concepts and
designs ofthe ignition facility for ICF, many new hydrodynamics models emerge, such
as HYDRA [53, 54]. PONHFZD [55], PF3D [13, 56-58], HARMONHY [59] and
PARAX [60-63]. HYDRA can be regarded as the successor of LASNEX for three-
dimensional simulation of laser-plasma interaction in the hohlraum. HYDRA is based
upon arbitrary Lagrange Eulerian (ALE) hydrodynamics, employing modem monotonic
artiﬁcial viscosities to stabilize shocks. HYDRA includes a ray tracing package for laser
beam propagation and takes into account ponderomotive effects, heavy ion deposition
and radiation transport. In addition. local themrodynamics equilibrium (LTE) and non-
LTE opacity models are also included in HYDRA. HYDRA has been utilized in design
and simulation of laser-driven experiments [64, 65]. In a later simulation [54], a kernel-
based nonlocal electron thermal conduction package based upon the model in [66] was
installed in HYDRA . enabling the code to handle the effects of long-range electrons.
However. there are still some physics missing in HYDRA [53]. including the power loss
from the beam due to laser power instability, beam steering due to the ponderomotive
force and magnetic fields that may affect the electron temperature proﬁle in the low
density regions inside the hohlraum.

In order to simulate the time-dependent ﬁlamentation and stimulated Brillouin

forward scattering in ICF plasmas, Schmitt and Afeyan developed a two-dimensional

ll

code PONHFZD [55]. The laser beam propagation in this model is achieved by a scalar
paraxial wave equation which is coupled to the plasma through the dielectric constant.
Two different approaches are used to model the plasma response. One is to use a two-
dirnensional Flux-Corrected-Transport (FCT) algorithm [67] to solve the continuity
equations for the plasma density and momentum, assuming a quasineutral one-fluid
plasma. The other approach uses a spectral FFT based algorithm to solve the
ponderomotive driven ion-acoustic wave (IAW) equation. The PONHF2D model
includes background flow and can treat nonlinear plasma behavior. Furthermore. the
derivatives ofthe fluids variables in the direction ofthe beam propagation are not ignored.
which enables correct description of both laser-plasma and plasma-plasma coupling in
the time evolution of ﬁlamentation. However, this model neglects the effects of thermal
ﬁlamentation and nonuniform inverse bremsstrahlung heating. The ionization level ofthe
plasma and the plasma temperature are also assumed to be known independently.

In contrast, interplay between strong inverse bremsstrahlung absorption, laser
pulse temporal proﬁle and target expansion are included in the model of HARMONHY
[59]. A similar three—dimensional model with HARMONHY is PF3D [13], which
includes a nonlinear hydrodynamic package as described in [68] coupled with a paraxial
solver for the laser propagation. PF3D can describe the (nonlinear) evolution of
ponderomotive and thermal ﬁlamentatin, forward Brillouin scattering, and such
phenomena as beam deflection in the presence oftransverse flows. Instead ofprescribing
the plasma parameters (number density, temperature. etc.) as does in PONHFZD,
simulation results of plasma variables obtained from other hydrodynamic plasma models

are input to PF3D as the initial conditions. For example, the results from SAGE and

HYDRA are used as the initial conditions for the work in [56] and [57], respectively. The
signiﬁcant advantage of PF3D over HARMONHY is its prominent parallel computation
performance using message passing.

Another hydrodynamic model that also uses the paraxial wave equation for the
laser beam propagation is PARAX [60]. In this model. nonlocal transport effects [66, 69]
are included by using the nonlocal Navier-Stokes equations [62]. Although PARAX can
describe various physics emerging in the interactions between laser and plasma, the
backseattering processes (SRS and SBS), the expansion of the plasma in the parallel
direction and a nonlocal transport model that would be valid in the strongly nonlinear
regime are still missing. In addition, the present PARAX simulations are implemented
with preformed plasmas. It is expected that the plasmas parameters obtained from other
plasma models need to be input to PARAX as the initial conditions.

The common feature of PONHF2D, HARMONHY, PF3D and PARAX is that all
ofthcm are using a paraxial approximation for the laser beam propagation in the plasma.
However, this approximation is invalid near the critical surface of the plasma inside the
hohlraum [56]. In order to overcome this problem, a hydrodynamic model called
KOLIBRI [70] was developed without the priori ofthe paraxial approximation. Instead, a
scalar model which corresponds to s-polarization in the two-dimensional case is applied
for the laser propagation, allowing for SBS in backward and all sideward directions as
well as its coupling with ﬁlamentation.

The hydrodynamic models are also widely utilized to simulate the laser-produced
plasma and the laser-plasma interaction in the ﬁeld of laser material processing. By

assuming that the electrons in the laser-induced plasma in C02 laser welding are in the

13

local thermal equilibrium, Finke er al. [71] proposed a one-dirnensional steady state
model to describe the laser energy deposition into the workpiece. Based upon a steady
state one-dimensional Boltzmann equation of the ions, the mass and momentum
conservations of the ions and the quasineutral plasma can be established through the ﬁrst
and second moment of the distribution function. Three mechanisms ofthe energy transfer
inside the keyhole: laser energy absorption through inverse bremsstrahlung, generation of
free electrons and recombination of charged particles on the wall, are considered in the
model and the energy balance equation can be established. A similar model was proposed
by Tix and Simon later [72], in which a system of transport equations for all species
particles inside the keyhole plasma is setup also based on stationary kinetic equations.
Both the convective and conductive transports and inverse bremsstrahlung heating are
included in the model. The model assumes cylindrical symmetry and no dependence of
the quantities on the direction of the keyhole axis; that is, the simulation was taken on a
slice of the keyhole, far enough from the surface. The simulation results suggest that the
main part of the plasma (except for the part near the wall) can be considered to be a fully
ionized plasma and can be described by the heat conduction. A model describing the
process of deep-penetration laser welding was developed by calculating the keyhole
proﬁle using a point-by-point energy balance analysis along the wall [73]. The absorption
of laser energy by the plasma is described by Beer-Lambert’s law where the absorption
coefficient is from the dominant inverse bremsstrahlung. The degree of ionization as a
function of temperature is described by Saha’s equation. In the above three models, no
laser beam propagation is considered. However, the laser beam propagation can be

largely affected by the plasma. As a consequence, the plasma absorption and the laser

14

intensity distribution can be altered a lot. In order to deal with this problem. a model
proposed in [74] uses the paraxial wave equation to calculate the laser beam propagation
through the plasma plume whose properties are pre-set within the calculation. The
simulation results indicate that the focusing of the laser beam is strongly deteriorated in
the plasma plume due to refraction and absorption. This plasma defocusing is believed
the main reason of the so-called “plasma shield” in laser welding. By keeping the
treatment of laser beam propagation with the paraxial equation, Kim and Farson [75]
proposed a axisymmetric model with compressible Navier-Stikes equations to simulate a
C0: laser beam impinging on a flat iron surface with helium and argon as the shielding
gas. Besides the refraction and absorption in the plasma plume. the effects of the
shielding gas and the reflection at the workpiece surface are also included in the
simulation. Chen and Wang developed a model [33] with a simpliﬁed three-dimensional
keyhole model to study both the keyhole plasma and the plasma plume. Radiation and
scalar laser energy absorption are included in the governing equations as the source temt.
But the laser beam propagation is not coupled in their model. Another three—dimensional
keyhole model for laser welding was proposed by Ki 6! a]. [76, 77]. The main feature of
this model is that the evolution of the liquid/vapor interface is consistent with the fluid
flow and the heat transfer in the model. The level set method is adopted to incorporate the
liquid/vapor interface boundary condition into the Navier-Stokes equation and the energy
equation. The ray tracing method is applied with the level set method to take into account
the multi-rcflection effects inside the keyhole. The laser energy dissipated into the wall is

obtained by the energy balance analysis on the wall surface.

Except for the applications in laser welding, the hydrodynamic models are also
exploited in the ﬁelds of pulsed laser ablation of metals and dielectrics, pulsed laser
deposition (PLD) [78-80]. and even laser ablation of organic materials [81, 82]. In
modeling of UV pulsed laser ablation of metallic targets, an integral model based on the
energy balance analysis of the different processes involved in the laser-solid interaction
was proposed in [83]. Laser energy absorptions due to inverse bremsstrahlung and
photoionization are considered in this model. Moreover, electron-impact excitation and
ionization. three-body recombination process, and energy transfer from hot-plasma
electrons to the heavy particles are included in the vapor/plasma kinetics. However, no
laser propagation is solved in this model. Another hydrodynamics model, MULTI—FS,
was developed in [84] to simulate subpicosecond laser interaction with solid density
matter. This method is based upon its previous version MULTI [85], which is a one-
dirnensional. one-fluid, two—temperature model including electronic heat conduction and
multigroup radiation transport. Comparing with MULTI, MULTl-FS is improved in the
following three aspects: 1) Maxwell’s equations are solved by the matrix method on a
high resolution mesh whereas a WKB approximation was used in MULTI; 2) electron
collision frequency is modeled to include its effects on the collisional absorption of the
laser, the energy transfer between electron and ion and the themral conductivity; 3)
SESAME equation ofstate are used for electron and ion separately. The code of MULTI-
FS is validated by several absorption measurement performed with aluminum targets
irradiated by a subpicosecond pulsed laser. In order to resolve the multiple timing scales
in the laser—induced plasma plume dynamics. a combined one-dimensional Lagrangian —

two-dimensional Large Particle Model was developed in [86]. This model consists oftwo

16

parts. The ﬁrst part describes the rapid laser-matter interaction within 4 ns. A one-
dimensional Lagrangian hydrodynamic model is used in the ﬁrst part where the two-
temperature approximation is applied. The laser energy absorption is treated by
numerically solving the Helmholtz equation coupled with equations of motion. The
second part models the two-dimensional plasma plume expansion until several hundreds
of nanoseconds with the results from the ﬁrst part as the initial condition. The Large
Particle Method. which includes Euler and Lagrange steps, is used to solve the system of
hydrodynamic equations in the second part. The most recent hydrodynamic model for
diagnostics of the plasma created on a surface of Ag target irradiated by intense
femtosecond laser pulses is proposed in [87]. This is a serniempirical model where the
numerical coefﬁcients are chosen so as to ensure the best accordance ofthe simulation to
the experimental measurements. The wave equations are solved for the laser beam
propagation and inverse bremsstrahlung heating is included in the hydrodynamic
equations as the energy source term.

As for the laser ablation of dielectrics, Jiang and Tsai [88] proposed a semi-
hydrodynarnic model because, in spite that the electron number density is govemed by a
Fokker—Plank equation, the temperature equation is still in the hydrodynamic form.
Another difference with other hydrodynamic models is that the heating term in this model
is .loule heating term. The ablations of barium aluminum borosilicate and fused silica by
a 220 fs pulsed laser with different fluences are simulated and the free electron number
density, the ablation depth and the reflectivity are obtained as the simulation results. In
order to interpret the effects of ionization on the interaction between femtosecond laser

and silicon, Li and Ki [89] proposed a model which solves the full Maxwell’s equations

17

by the ﬁnite-difference time-domain (FDTD) method. The generation of free electrons is
obtained by a local energy balance analysis. The simulation results show that the
ionization of silicon by the incident laser beam makes the silicon behave metallically.
The laser energy is absorbed mostly in the thin layer below the target surface and the
laser beam propagation is even blocked. A threshold characteristic in the regime of
femtosecond laser ablation is clearly seen from the results.
1.3.2 Kinetic Based Models

Although hydrodynamics simulations of laser-plasma interactions have gain great
success in many ﬁelds, they tend to be invalid at high laser intensities where collisionless
and relativistic characterize the laser-plasma interaction. A relevant research indicates

. . . . . . . . 17

that the validity range of laser mtensrty for hydrodynamic Simulations is less than 10
W/cm2 [84]. Moreover, even with low laser intensity, hydrodynamic models are still
invalid ifthe Knudsen number of the problem is large, e.g., when a laser-induced plasma
plume expands in vacuum [90].

From the ﬁrst, simplest one-dimensional collisionless model developed in 1950’s
[91] until today’s sophisticated parallelized three-dimensional collisional code. particle-
in-cell (PIC) method [92, 93] has become the most popular methodology for plasma
modeling, especially in simulating laser-plasma interactions with high laser intensities.
The principle of PIC methods is to solve the equations of motion ofcharged particles in
the self-cor‘rsistent electromagnetic ﬁelds that are calculated from the full set of
Maxwell’s equations. The particles are weighted to the spatial grids so that they have an
effective size on the order of the grid, allowing one to model a plasma problem using

fewer particles than in an actual experiment. Except for the numerical errors generated by

the discretization of space and time, PIC methods retain the physics to the maximum
extent theoretically because there are no physical assumptions (such as continuum
approximation) presented. However, due to the limitation of computational expanse and
the diversity of physical characteristics for different problems, there are more or less
approximations applied on PIC codes in the practical implementations.

Nowadays. the representative PIC codes for laser-plasma interactions include
WAKE [94, 95], OOPIC [96-98], QUICKPIC [99], CALDER [100—102], and etc. The
algorithm of WAKE consists of three approximations [94]. The ﬁrst approximation is
referred to as ponderomotive guiding center motion ofthe particles whereas only the low-
frequency motion of electrons is retained. However, this approximation breaks if there
are self-trapped electrons in the plasma or the axial speed of electron approaches to the
speed of light. The quasi-static approximation is the second approximation in WAKE. It
assumes that the proﬁle ofthe laser pulse dose not change signiﬁcantly during the time of
interaction with the electrons. As the third approximation, extended paraxial
approximation is applied in WAKE to simulate the laser pulse prOpagation. Based upon
the above three approximations. Chessa et a]. [94] applied WAKE including a tunneling
ionization model to simulate short laser pulse self-focusing in underdense plasma.
Besides the tunneling ionization, the impact ionization due to the binary collisions
between electrons and neutrals are also included in OOPIC, which is a two-dimensional
relativistic object-oriented PIC code [96, 97]. The main feature of OOPIC is that a
moving window scheme is adopted because only the small vicinity around the laser pulse
in the plasma is of interest. The moving window scheme in which the computational

window moves with the laser pulse was ﬁrstly adopted in the two-dimensional cylindrical

code of ISIS [103]. In a later work, Decker er a]. [104] extended this scheme to the two-
dirnensional Cartesian geometry and simulated forward Raman scattering induced by
short pulse high intensity lasers. The above introduced PIC codes have a common
drawback: highly computationally expansive even with today’s fastest computers [99].
Encouragingly. a newly developed, full parallelized, full three—dimensional PIC code.
QUICKPIC [99]. can largely save the computational time. The main assumption in
QUICKPIC is that there are two distinct time scales for the laser pulse and the plasma
evolution. Besides. the ponderomotive guiding center motion of the particles and quasi-
static approximations are also included. Simulation of a benchmark problem by
QUICKPIC and a conventional PIC code OSIRIS [105] shows that this new code can
reduce the computational time by 2—3 orders of magnitude while reproduce satisfactory
results. However, the ionization mechanisms are not considered in QUICKPIC. Another
representative two- or three-dimensional, parallelized, full relativistic PIC code is
CALDER [100], which has been applied to simulate the fast electron energy deposition in
fast ignition [102]. proton acceleration [101] and laser-accelerated ion beams [106].

As another particle simulation method of plasma. molecular dynamics (MD)
models are also proactive in simulation of laser—plasma interactions. Unlike PIC, which is
based on particle-mesh (PM) scheme, MD relies on particle-particle (PP) or particle-

particle-particle-mesh (P3M) schemes [93]. Although PM is the most widely applied in

particle simulation, the PP is preferred in some special circumstances, for example, in the
many-body collisions dominated plasmas. In addition, PP can handle three-dimensional
problems just as easily as two-dimensional ones, with insigniﬁcant increase in

computational expanse [107]. However, since solving the individual equation ofmotion

. . 2 . . . .
for each ofN particles would requrre N /2 rndrvrdual force calculations to be undertaken

every integration time step, the main drawback of PP scheme is apparent: the

computational expanse is much higher than that of the PIC approach. In order to utilize
the advantages of both PM and PP to the greatest extent, most MD codes applies P3M as

the underlying scheme, where the force on a particle is divided into a collective, long-
range term from the majority of the particles (PM) and a short range term arising from
particles close to the particle of interest [108]. The main feature of MD models is that the
collisions between particles are inherently treated as many-body collisions instead of
binary collisions in most of PICS. Thus, MD models have been applied to simulate
collective collisional phenomena in plasmas, such as inverse bremsstrahlung [108, 109]
and interactions between KrF laser pulses and plasma cluster with ionization [107].
Gibbon er al. [1 10] also uses MD to simulate the proton acceleration from laser-irradiated
targets in which their calculation is based upon tree algorithm and uses massively parallel
computers.

Comparing to PIC and MD, both of which are particle simulation ofplasrnas, the
Boltzmann equation with Fokker-Planck collision term governing the electron (ion)
collisions can be directly solved to simulate laser-plasma interactions. The ﬁrst numerical
solution of Fokker-Planck-type Boltzmann equation [111] was obtained in the study of
the effects of steep temperature gradients appearing near the critical surface on the
electron heat transport, where the Boltzmann equation was solved by an expansion ofthe
distribution function in spherical harmonics. In a later work [112], a one-dimensional
Fokker-Planck code, FPI (Fokker-Planck International), was developed and inverse

bremsstrahlung and ion motion were included in the model. In the past two decades, FPI

21

was developed and improved by adding more and more physics into it, such as Coulomb
collision, cold ion hydrodynamics, transport for each energy group and self-consistent
electric field to ensure the quasi charge neutrality [113, 114]. In a newer version of FPI
[115], collisional ionization and excitation, three body recombination and de—excitation,
line emission, and the high density effects of continuum lowering and pressure ionization
were included. In this code, the electromagnetic wave was calculated by a Helmholtz
equation solver while the system of equations was solved by ﬁnite difference and time-
splitting methods. In order to save the computational expanse. an average ion model
[I 16], which is often used in hydrodynamics code, was adopted in this code. Besides FPI,
there are also other codes to solve the Fokker-Planck-type Boltzmann equation. For
example, Town er al. [117] used their own Fokker—Planck solver to simulate the short-
pulsed laser interaction with solid but no ionization mechanism was included. In a
subsequent work [118], they added collisional ionization and recombination and also
applied the self-consistent average ion model.

Except for the above mentioned kinetic models, there are other approaches to
simulate laser-plasma interaction on the kinetic base. For example, a cellular automaton
(CA) model was developed to simulate laser-plasma interactions [119]. Although there
are many assumptions applied upon that model, such as no ionization and recombination,
ﬁxed ion background and quasi-static electrons, the preliminary results show an
encouraging potential to apply CA on laser-plasma simulations. In Ref. [90], Itina er (1!.
introduced a hybrid model to simulate the laser-induced plasma plume expansion. In the
first stage of plume expansion, a one-fluid two-temperature gas dynamics model was

used to describe the plume and the Large Particle Method [120] was adapted to solve the

IQ
Ix.)

system numerically. When the plume expands into a region where the gas dynamics
model becomes invalid due to the steep density gradient and fails to describe the
diffusion process. the direct simulation Monte Carlo (DSMC) method [121] was applied.
Relativistic Vlasov simulation [122] and propagator ofdistribution function [7] are other

examples to simulate laser-plasma interactions by kinetic based methods.

1.4 Main Content of This Research

In the present research, a new model is proposed to simulate interactions between
laser and weakly ionized plasmas. In this model, Boltzarnnn equation is used to describe
the transports ofdifferent species ofparticles in the plasma. The contributions ofelectron
impact ionization and three body recombination are treated as part of the collision terms
in the Boltzmann equation. Besides, the elastic collisions between particles are also
included by using the Bhatnagar-Gross-Krook (BGK) approximation. The laser beam
propagation inside the bulk of the plasma is simulated by solving the full set of
Maxwell‘s equations directly. The Maxwell’s equation and the Boltzmann equation are
coupled by the external force term and the current density. In addition, energy equations
of all particle species in the plasma are derived directly from the Boltzmann equations
and are solved by finite volume method. The advantage of this present model over other
hydrodynamic models is that less assumption are made so that the physics of laser-
plasma interactions can be retained to a great extent, especially for the wave-related
phenomena since they are described by the full set of Maxwell’s equations here.
Comparing to the kinetic based models, this model can be applied in much longer length

and time scales with relatively low computational expanse.

Ix)
L’J

In Chapter 2, basic theory of lattice Boltzmann method (LBM) with BGK
approximation of the collision term is described, including derivation of the lattice
Boltzmann equation from the continuous BGK—type Boltzmann equation and the
recovery ofthe macroscopic equations by Chapman-Enskog expansion. As an application
example of LBM for incompressible hydrodynamics flow, the two-dimensional driven
cavity flow is simulated to show the capabilities of the standard LBM. Different
evaluation methods ofthe external force term are introduced in Chapter 3. Poiseuille flow,
Taylor vortex flow and free diffusion problem under a uniform applied external force are
simulated to compare the different external forcing models. Based upon the comparison
results, one external force model is selected for all the plasma simulations in this research.
In this chapter, magnetohydrodynamic (MHD) flows, including Hartmann flow, Orszag-
Tang vortex system (both in 2D and 3D) and magnetic reconnection problem, are
simulated by a new developed hybrid LBM model, where the induction equation is
solved by the ﬁnite difference method and it is coupled with the LBM through the
external force term. Although the standard LBM are successful in many applications of
fluid flow problems, it may produce large simulation errors if the physical properties of
fluid are used in the calculation. It gets more deteriorated in plasma simulations. In order
to overcome this problem, a rescaling scheme is introduced in Chapter 4 to detemrine the
simulation parameters ofLBM. Isothermal weakly ionized helium plasma is simulated by
using LBM with the rescaling scheme in this chapter. The new model as well as the
rescaling scheme is validated by the simulation results of electrostatic wave in the plasma
and electron diffusion problem. In Chapter 5, the LBM model developed in Chapter 4 is

extended to simulation of laser interaction with weakly ionized plasma. This model,

24

which is essentially a thermal LBM model, takes into account for the temperature
evolutions of individual species ofplasma particles so that the electron impact ionization
and its reverse process. three-body recombination, can be included. Preliminary results of
interaction between a continuous laser beam and a weakly ionized helium plasma are
presented in this chapter. The capabilities, limitations and future work necessary for

improvement ofthe present model are concluded in Chapter 6.

25

Chapter 2

Lattice Boltzmann Method for Incompressible Hydrodynamics Flows

In recent years, lattice Boltzmann method (LBM) has been widely applied as a
successful alternative scheme to simulate fluids flows [123]. Unlike the traditional
numerical methods which solve the Navier-Stokes equations for macroscopic variables,
LBM is based on kinetic theory and solves the Boltzmann transport equation for the
evolution of particle distribution functions. The macroscopic averaged properties which
obey the desired macroscopic equations will be obtained from proper moments of the
distribution function. LBM first originated from its Boolean counterparts, the lattice gas
automata (LGA), but it has been proved that it can be derived directly from the
continuous Boltzmann transport equation by descretization in both time and phase space
[124, 125], which gives LBM a more rigorous theoretical foundation. In this chapter,
some basic theories of LBM for incompressible hydrodynamics flows, including the
BGK approximation of the collision term, derivation of lattice Boltzmann equation ﬁ‘om
the continuous Boltzmann equation and recovery of macroscopic conservation equations

by Chapmann-Enskog expansion, are introduced. As a simulation example, a 2D driven

cavity flow is simrrlated by the standard D2Q9 LBM to illustrate its capabilities.

2.1 Boltzmann Transport Equation and Its Collision Term
The Boltzmann equation, a nonlinear seven-dimensional (in phase space) partial
differential equation, describes the evolution of the phase space distribution functions of

different species of particles in various types ofgascs or gas mixtures, including plasmas.

The continuous Boltzmann transport equation for the particle species of s in the fluid has

the following form:

i4" V5. 'Vx/Is' +35 'vas 2 Q5 (ft/’) (1)

(71
where f, = f, (x,vs ,1) is the single particle distribution function in the phase space, v3.
is the microscopic velocity and a, is the acceleration of fluid particles caused by the
external force. Based upon the assumption that only binary collisions in the fluid are
taken into account [126], the collision term Q,(f.f’) representing the change of rate of

the distribution function due to the collisions can be expressed as:

472'
Q_.-(f»f')= ldv.» I (190...»(Qilv...-v.»|[frv;)frvf,«)—frv.)f(v.«)] (2)
I 0

where the subscript s and 5' represent particle s and particle s' which are involved in
the binary collision. the superscript prime stand for the post collision state, Q is the solid
angle in the collision and 05.540) is the differential collision cross section of the
collision between particle s and particle 5'. When s'=s , the above formation of
collision term describes the self-collisions; otherwise, it represents the cross-collisions
between different species of particles in the fluid.

After we obtain the distribution function by solving the Boltzmann equation. the
desirable macroscopic quantities can be deduced by proper moments of the distribution
function. For example, the number density, macroscopic velocity and thermal energy are

given by:

11,: [ftg(x,vs,r)rlv5 (3)

If.

 

"Sub' : I stsix’ V5,!)t/V5 (4)
1 y"
m = *2- I (vs «oz/xx. v.-,Iidv. (5)
—30
D If T. . . .
where 55 2 ~31 8 5 and DO rs the degree of freedom ofthe fluid particle.
m

S

Since the Boltzmann equation is essentially a statistical treatment of the
microscopic behaviors of the particles inside the gas or plasma, it is more fundamental
than Navier-Stokes equation which is based on the assumption that the fluid of interest is
continuum macroscopically; unlike Navier-Stokes equation which is only applicable for
small Knudsen number (usually less than 0.1), the Boltzmann equation can be applied in
a much larger range of Knudsen number [127] and thereby it is especially suitable to
describe the complex fluid flows, such as non-ideal fluids, multiphase or multicomponent
fluids and those which cannot be described by Navier-Stokes equation accurately and
effectively.

Unfortunately, the nature that Boltzmann equation expressed by Eq. (I) with the
collision term (Eq. (2)) is an integral-differential equation makes it almost impossible to
ﬁnd analytical or even numerical solutions [126]. An alternative way to evaluate the
collision term is provided by the relaxation model among which the BGK model [128] is
the most commonly used. In this model it is assumed that a situation initially not in
equilibrium reaches a local equilibrium condition exponentially with time, as a result of

collisions, within a relaxation time 4y The local equilibrium state of the particles is

characterized by a local equilibrium distribution function [SW/(xxx) which often takes

the form ofMaxwell—Boltzmann distribution:

7
.U n. (v, -u\ )-
.I...rI————II-——--—I «II
(““705 ) i “()5

where D is the dimension of the space and 65 =,//cBTS/ms is the sound speed of the
fluid particle s with It); as the Boltzmann constant, Ts as the temperature and ”’3 as the

mass of the fluid particle s. If there is only one type of particle existing in the gas, the

collision term can be mathematically expressed as

_ (f3 _ freq )

7
A ( )

Q5 ([3 ’ ff) :

.8'
Finally, by neglecting the subscript s, the single-particle Boltzmann equation with BGK

approximation for the collision term can be written as:

T+v~fo+a-va:————— (8)
Cl zl

c7 .f—./“"’

Note that the simpliﬁcation of collision term by BGK model expressed by Eq. (7)
is only valid for short-range elastic collisions. Physically, there are also long-range
collisions inside the fluid, such as the Coulomb collision between charged particles in
plasmas. By considering that the charged particle encounters a series ofconsecutive weak
(small deflection angle) binary collisions [129], one can derive the Fokker-Planck
collision form to describe the long-range collisions from Boltzmann equation. The

general form of Fokker-Planck collision term for the charged particle s can be written as

[130]:

29

 

 

,
5 AV 1 8“ AvAv

0...=——- — ‘. +— : (9)
" ’ av <Az>fs 2 am < >f‘

The vector <Av/AI>a and the tensor (AvAv/AQC! describe the dynamical friction and

diffusion in velocity space. respectively. They are given by

 

 

 

 

. 4 3
' ’ 41 . 1 -' 3+ ' "
<:§L> : 73.» 1,, A :2 LIL ”’ If" (V5,) ch" (10)
A! .x m: C" x, (I.\ m , lV - V l
4 7
AvAv 4m . 6“ c
< > : £3 lnA J13 Ifsw(‘s')I‘—VI(IV (II)
A, .5' HI; 6‘ C’VZ

where f,«(v,.I) :f,I(x.v5I,t) is the distribution function of particle species 5', and (I,
and (/_v' are the charges of particles 5 and s', respectively. The summation is over all
charged components including particle s itself. A is the plasma parameter which is the
number ofcharged particles in the Debye sphere. In practical calculation, many simpler
forms of Fokker-Planck equation as well as numerical algorithms are adopted to solve Eq.
(9) [I3I-133].

In weakly ionized plasmas. since the number density of neutral particles is much
greater than that of charged particles, the collisions with neutral particles dominates the
whole process. Thus we neglect the long-range Coulomb interaction between charged
particles in this research and the collision term approximated by BGK model as
illustrated in Eq. (7) is applied to evaluate the elastic collisions between charged particle

and neutral particle in all simulations ofthis research.

30

2.2 From Boltzmann Equation to Lattice Boltzmann Equation
By neglecting the external force temi in Eq. (8), the single-particle Boltzmann

equation can be rewritten as [125]:

 

cq
5L1: f (17)
(ll 21. A.

(I 5 . . . . . . .
where T = :— + v - ‘7‘ is the Lagrangian derivative along the microscopic velocrty v and
( I CI ‘

the equilibrium distribution function f6" takes the form illustrated by Eq. (6). By
integrating Eq. (12) over a small time interval AI , we can get:
—Ar/’/1 A!

‘ e ./ ., _, /
/(x+vAt,v,t+Ar)= ,2 Ie‘S/ﬂft(l(x+vs,v,t+s)ds+e A""1f(x,v,r) (13)
. 0

 

where 0353A! . Assuming that the time step Al is sufficiently small and f?" is
smooth enough locally. f“"’(x+ vs,v,t +5) can be expanded by Taylor expansion with

neglecting the term of order O(Ar2):
f"‘/(x+ vs,v,t +3) :ES—f"‘/(x+vAt,v,r+AI)+[1——As¥]_f‘i‘/(x,v,t) (14)
I 1

Substituting Eq. (14) into Eq. (13) and after some manipulations. we can get:

f(x + V’Al, v,r + At) — f(x, v.t) = (e—A’M — l)[f(x, v,t) — f"‘/ (x, v,1)]

. I, 1 ) (15)
+(1~—§—+-§-6’A”’1)1f‘q(x+ vAt,v,t +Ar)—f‘"(x,v,t)]
t r

—;\I//i -

Again, expanding e in its Taylor expansion and neglecting the terms of order

()(AIZ) , we obtain the Boltzmann equation discretized in time as:

_/'(x+ vAI,v,I +At)—f(x,v,t) = ——A—t—[f(x,v,I)—f"q(x,v.t)] (16)
r

31

where r : xl/Ar is the dimensionless relaxation time. Note that Eq. (16) has first order
accuracy in time [125].
Eq. (16) is still continuous in the phase space and v e (—w,+oc) for 2D problems.

As illustrated in Eqs. (3) to (5), the hydrodynamics variables are proper moments of the

distribution function f(x,v,l) (or equilibrium distribution function f"q(x,v,t)) with

respect to the microscopic velocity. This can be generalized as:
I 7:
ram: 1/‘(x,v,z)z(v)dv = 1 f ""(x. v,r>z(v>dv (17)
-—I —’f‘

where ¢(x,z) is the hydrodynamics variable, which is dependent upon both position and
time, and ;_/(v) is a function of the microscopic velocity v. Assuming that the

temperature is constant and the macroscopic velocity ofthe particle is small (low Mach
number approximation), the equilibrium distribution function can be expanded with

Taylor series up to the order of (u-u) [124, 125]:

 

 

,2 . . 2 2
“"l =————————’2’ 0,2 exp ——-V 2 1+(v 2")+(v "4) — " 2 (18)
(276 ) 219 6 20 20
Substituting Eq.(18) into Eq. (17), ¢(x,t) can be expressed as:
’J.‘ V2
¢(x,!)= [exp ———7— (v’")dv (19)
_x 20"

where t//(V"’) is a m-th order polynomial of v. The integration in Eq. (19) has the

2
structure of c7" and the integration limitation is from 00 to —-oc. Enlightened by the
above observation. it is intuitive to apply Hermite quadrature to find the value of Eq. (19)

exactly to the order ofm:

32

7
e".
¢(.r.r)=ZWa exp -——“3 /(eg’) (20)

C!
where 11"“ and ea are the weights and abscissas (or discrete velocities) of m-th order
Hermite polynomial. respectively. Then the hydrodynamics variables can be computed by

the quadrature as well:

n = if“ (21)
Cl
nu = Zea/a (22)
a
us 2 %Z(ea -u)2fa (23)
a

According to the dimensions ofthe physical problem (1D, 2D or 3D) and the number of
macroscopic variables desired, different number of abscissas of the Hemiite polynomial
may be needed. and accordingly. resulting in different sets of discrete velocities. For
example. the D2Q9 discrete velocity model may be applied for 2D isothermal problems.

To derive this nine-bit lattice Boltzmann equation, a Cartesian coordinate system is used

and w(v’”) is set as w(v"') = i)_¥’i",’. where 1 S n S m . Then Eq. (19) can be rewritten as:
(00”,!) = (ﬁglm'HZ/m/n (24)
where
x ,2
1m = J‘e—b Qym‘ E: (25)
—3c

33

and c; = xix/J25 or t: =v‘./\/20. In the derivation of the nine-bit model in 2D, it is

natural to use the third order Hermite quadrature to evaluate Eq. (25); that is,

3
,VIII , . , . ' . . ' . . ‘ .
I," = Z (015/ , where the three abscrssas and the corresponding three weights are.
i=1

Then Eq. (24) becomes:

4 8
7 7 7
¢(x,l) : 20‘wju/(O)+ Z (alwzu/(ea)+ Z mfg/(ea) (27)
a=1 0:5
where the discrete velocities ea are:
(Or 0) a = 0
ea 2 (cosrpmsin (12a )c' a 21.2.3.

4
(cos (pa ,sin (pa )x/2c‘ a = 5.6, 7,8

 

 

 

 

 

 

6 2 5
A
0
3< >1
Y
7 4 8

Figure 1 Schematic of D2Q9 square lattice. The two-dimensional continuous phase
space is discretized by nine microscopic velocity components in a square lattice.

34

with

{ma =(a—1)7I/2 a=1.2.3,4 (29)

(pa = (a—5)7r/2+7r/4 a = 5,6,7,8
where ('zx/30 is the lattice speed and is deﬁned as ('.=. Arr/AI. The schematic of the
discrete velocities is shown in Figure 1.

Comparing Eq. (20) and Eq. (27), the weights W0 is Eq. (20) can be found as:

 

2
113a = 2,7612 exp ea) (ya (30)
26”
where
'4/9 (1 =0
ma 2 1/9 a=l,2,3,4 (31)

1/36 a = 5.6.7.8

Then the discretized equilibrium distribution function is:

 

-

7
, 3 . c . -
((l _ I ‘L’(] _ i _
fa —Ilaj (x.ea,l)—a)an 1+ 2 + 4 2 (32)
(' 2(' 26

Note that in the above derivation. the phase space and the configuration space are
discretized simultaneously through the relationship between the lattice speed and the

spacing. Finally, the discretized Boltzmann equation. or lattice Boltzmann equation.

takes the follwing form:

/"(X~I)- .IJ‘HM)
Z'

 

_/"(x+eaAI,I+AI)=f(X.1)—‘ (33)

In the practical implementation of standard LBM, there are always two steps
involved: collision step and streaming step. In the collision step, the fluid particles collide

with each other on the local node point to finish the momentum exchange and energy

35

exchange (if it exists in the physical problem). The particles then stream to the
neighboring nodes along the directions of the dcscretized velocity vectors in the
streaming step. These two steps can be illustrated mathematically by:

Collision step:

fa(x9’)—,/;)(l(xvt)

 

f~U(x.t)'-'fa(X,I)- (34)

Streaming step:
‘/‘a(x+eaAt,t+At)=fa(x,t) (35)
where fa(x.t) represents the post-collision distribution function along the direction a.

From Eqs. (34) and (35). we can see that the collision step is purely local and the
streaming step is a uniform data shifting and only demands little computational effort.
The two-step scheme is explicit, easy to implement and straightforward for parallel

computation.

2.3 From Lattice Botzamnn Equation to Navier-Stokes Equation

The desired macroscopic conservation equations, such as the Navier-Stokes
equation, can be recovered by the Chapman—Enskog expansion of the lattice Boltzmann
equation as shown in Eq. (33). The significance of this derivation lies on the fact that
some fluid properties, such as the kinematic viscosity. can be retrieved from this
derivation.

Introduce the following multiscale expansions to the a -component distribution
function and the time derivative:

-() ~ 1 ,2 ‘ 2 , ~
ta =12. Hat; We 1.: ’+--- = Z bra/.5,“ (36)
1:20

36

 

a (3 a 5 a
-ﬁ—=-—-+8—-+82:-—+--°=Z£”ﬁ (37)
(I Clo 8!] C212 ":0 Cl”

where 8 is the time step At.

By taking the Taylor expansion to /;1(x + ea At,t + At) , we can have

fa(x +eaAt,t + At) = fa(x,t) + 0,]? (x,t) + D,2fa(x.t) + ()(At3) (38)

'l
where D, = —:’—+ea 57. Substituting Eqs. (36), (37) and (38) into Eq. (33) and sorting the

0!

resulted equation by the order of 8, we can get:

 

0(50):
((0) =f5" (39)
, ‘1 .
0(1, ).
1
Dzoféo’ =7- 3.” (40)
C(52):
1) IO) _ .
(fa +(2r DD, /;1)=_lfc((2) (41)
all 22' 0 T

q
7

c

where D, 5(—
II 7.

" II

+ea-V). The distribution function in the above expansion is
constrained by:
(0) l 3 n
Zfa : _
0 ea‘ nu

1 (42)
r

m~>[ :0 ">0

0 ea—v

 

Taking the zeroth and first order moments oqu. (40). we get:

37

~——+V- -(nu)=0 (43)
alo

and

C(nu)+ 1V "(0): (44)

 

C to +m

where “(OlznzZeaeafém is the zeroth order momentum flux dyad and

Cl
"(0’ = Inna/“l + uu) according to the kinetic theory. Similarly. by taking the zeroth and

first order moments of Eq. (41), we can get:

—‘—_—.0 (45)

and

(".‘(nu)+ (2 —l) I

811 21' ”I

 

—V [1‘1 =0 (46)
where [1(1) = mZeaea/‘én is the first order momentum flux dyad. From Eq. (40), we
CI

()0 END) (0)
have” )z—rZeae ea D,0/'(1) =—r ——(ﬁl—[—+V-eaeaea/a , where the second term
" O

on the right hand side reads V-eaeaea. “(0) 21211102 [(V u)l + 2Vu] [134]. By neglecting

(1

7.
the terms greater than the order of O(u”) [134], we can get:

7 _
ELEV-Hm: —(2T———1—)tnnr‘)2 V“u (47)
2?
Letting the kinematic viscosity v
v = (r —().5)(}2At (48)

38

Eq. (46) can be rewritten as:

 

 

(3‘ nu t' -.7
(. )— ”\r’ 'u = 0 (49)
CI] A!
. a 0 8 .
Recalling that ::—+At:— (as shown in Eq. (37) and truncated at the order of
CI Eto Ctl

()(A12)) and using Eqs. (43), (44), (45) and (49), we can finally recover the continuity
equation and the Navier-Stokes equation from the lattice Boltzmann equation (Eq. (33))
as:

Continuity equation.“

 

0n
—+V-(nu)=0 (50)
(2
Mtvic’r-Sto/tcs equation:
“.‘ 1
dim)+V-(nuu)=——Vp+vnV2u (51)
o! m

where p = lI/CBT is the scalar pressure.

From the above derivation, we can see that there is a modification of the fluid
viscosity by 0.5 as illustrated in Eq. (48), which leads that the viscosity is not only
dependent on the physics but also dependent on the lattice. This means that ifthe physical
properties of fluid are used in the simulation. the resulted viscosity may not be equal to
the physical one. This problem can be solved by using a rescaling scheme developed in

this research. It will be explained in details in Chapter 4.

2.4 LBM Simulation of Two-dimensional Driven Cavity Flow
In order to illustrate the capability ofthe standard D2Q9 LBM as described above.

a 2D incompressible driven cavity flow is simulated on a uniform 256><256 square grid

39

system. The driven cavity flow which is bounded by a square enclosure and driven by a
uniform outer ﬂow on the top shows rich vortex phenomena at many scales depending on
the Reynolds number [135-137]. In this simulation, the no-slip conditions are applied at
the left. right and bottom walls; that is, the macroscopic velocities at the boundary node
points on the wall are equal to zero. Viewing from the microscopic point. the distribution
functions just reverse their direction but keep their values on the wall. For the top
boundary. since the outer flow remains unchanged physically. the equilibrium values
with the outer ﬂow speed ur are set to the top nodes as the updated distribution function
after each streaming step in the calculation. The simulation parameters are set as: Ar = 1,

At =1 and po = 1.0. The driven velocity on the top of the cavity is u3C 20.1 and the

V
3”“f. A .\'

Reynolds number is defined as Re 2 where NY is the grid number along the top

(I — 0.5)
boundary and the viscosity used here is calculated from Eq. (48).
Figure 2 shows the streamlines of the cavity flow at different Reynolds numbers.
The dependence of the vortex on the Reynolds number can be seen very clearly from the
figure. Figure 3 illustrates the x- and y—component velocities through the geometric
center of the cavity. Their patterns agree with those in the literatures very well. The
detailed quantitative comparison with results in the literatures will be presented in

Chapter 4.

4O

.oy/Lo

.0

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Re = 1000 Re = 2000
Figure 2 Streamlines of driven cavity flow with different Reynolds numbers. The

flow pattern is strongly dependent on the Reynolds number. Ifthe Reynolds number
is sufﬁciently large, secondary vortices appear at the right and left lower corners.

41

 

1.0 . 1

     
 
  

 

 

 

 

   

 

 

 

0.8 3
_r 0-6 ‘ a Re = 100
3‘ b Re = 200 r
c Re = 500
0'4 t d Re = 1000
e Re = 2000
0.2 - .
O-O ' T i I ' I I I l '
-0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0
Ux/UO
(a)
04 r l ‘ l ' l ' l
0.2 -
o 0.0 -
2
5 a Re = 100
-0'2 4 b Re = 200
c Re = 500
(1 Re = 1000
-o_4 - e Re = 2000
'06 ' l V r ' r v 1
0.0 0.2 0.4 0.6 0.8 1.0
x/L
(b)

Figure 3 Velocity proﬁles through the geometric center of the cavity (a) Proﬁles of
x-component velocity; (b) Proﬁles of y-component velocity

Chapter 3

Evaluation of the External Force Term in Lattice Boltzmann Method

In the simulation of the driven cavity flow presented in the previous chapter, no
external force exerts on the fluid particles. Accordingly, the lattice Boltzmann equation
takes the form illustrated by Eq. (33). But in most practical cases, especially in the
simulation of plasmas. there are always various types of external forces acting on the
particles. affecting the rates ofchange of momentum and energy. It is very important to

evaluate the external force term in the lattice Boltzmann equation accurately. However,

the external force term a-va in the Boltzmann equation cannot be calculated directly

because the dependence of the distribution function on the microsc0pic velocity v is
unknown. Thus, it is necessary to develop alternative methods to evaluate the external
force term in the Boltzmann equation. In this chapter, five different external forcing
models are compared through simulations of Poiseuille flow, Taylor vortex flow and free
diffusion problem under a unifomi externally applied force. Based upon the comparisons,
He’s model is selected to evaluate the external force term in this research. By using He’s
model, a new hybrid LBM model is developed to simulate magnetohydrodynamics
(MHD) flows. This model is validated by Hartmenn flow. Orszag-Tang vortex system

and magnetic reconnection problem in this chapter.

3.1 Different External Force Models
The lattice Boltzmann equation with presence of an external force term can be

written as [138]:

43

fa(x+eaAt.t+At) =fa(x,t)—

 

_ er]
fu(x’{) f“ (x°{)+AiF(1 (52)
T

where Fa is the discretized external force term along the a -th direction of the discrete
velocity vector. The corresponding equilibrium distribution function in Eq. (52) is

defined as (in D2Q9 scheme):

:1: a: 2 *2
3(e. -u) 9(e -u) 3n
“7 + a -— 7 (53)

c~ 2C'4 26"

 

 

er] _ .
fa —(.)an 1+

where u* is called the equilibrium velocity. Different formations of Fa and u* lead to
different treatments of the external force term in the lattice Boltzmann equation.

The very first and simplest external force model was developed by He at a]. [139]
and is referred as to Method 1 in this research. In this model. the external force term Fa

(Ha/lea -a
7

0..

takes the form as Fa = where a is the acceleration of the fluid particle due to

. . . . . It . . . .
the external force. The equilibrium velocrty u IS equal to the fluid macroscopic velocrty

u; that rs, u* = u = :z’zea/a . Thrs method rs suitable for the case where the spatral and

a
temporal gradient of the external force is negligible. Thus. it is always applied for the
flows under a constant body force, such as constant pressure gradient or gravity.

The second method to evaluate the external force term (which is referred as to
Method 2 here) was originally developed to simulate nonideal fluids by the LBM [140-
142]. It was originated independent by two groups within the frame of Hermite expansion,
which is consistent with the derivation of the lattice Boltzmann equation. and with some

constraints ofdifferent order moments. In this model, the external force term is written as

44

(ea —u*) (ea -u*)ea
+
02 04

 

Fazwana- and the equilibrium velocity has the same

- . . . . . . a 1 .
definition as the fluid rnacroscoprc velocrty. 1.e. u = u = —Zeafu . In this model. the
n
(1

contribution ofthe external force to the momentum flux is considered. But like Method 1.
it is still only applicable for the case where the spatial and temporal gradient of the

external force is sufficiently small. In a later research [143], this model is updated by

- . . . . 1 I . .. .
defrnrng the fluid rnacroscoprc velocrty as u =—Zeafa +—7—Ata while the equrlrbrrum
n , -

(1
velocity keeps the same definition as before. This approach eliminates the temporal
gradient constraint of the external force. but still cannot be applicable for the large
spatially varying external force. This updated version of Method 2 was applied by
Brcyiannis and Valougeorgis [144] in a recent research to simulate 3D
rnagnetohydrodynarnics flows where the magnetic force J x B is treated as the external
force.

Method 3 was proposed by Buick and Created [145]. The external force term in

 

their model is Fa

:[l—ij (0‘12”(ea-a). The equilibrium velocity and the fluid

0

. . - . . a: 1 l
macroscopic velocrty have the definition as u =u=— E eafa+3Ata. In a later
it ..
Ct

research, Guo et a/. developed a model (which is referred as Guo‘s model) where the

 

xi: :1:
(ea—u )+(ea-u )ea
7

4 and 170th
0“ 0

- . 1
external force term rs defined as Fa :[1—2— mana-
r

the equilibrium velocity and the fluid velocity have the same definition as Method 3.

The most popular external force model was developed by He er a]. [146] (which

is referred as to He’s model here) and it has been widely adopted in many researches
[147-149]. This model is based upon the fact that fa] is the leading part of the

expansion of the distribution function and the gradient of fa] has the most important
contribution to the gradient of f. Then the external force term can be written as:

5...va 2. a.V‘-I./()q :_E.iv7;l£2

‘—

.f°“/ (54)

Accordingly. the Boltzmann equation with the external force term becomes:

 

(3f f‘rfm a'(V-—l|) a .
—+v-V =—‘ + I (55)
at "f 2 93 f

By following the similar discretization approach illustrated in Chapter 2, the lattice

Boltzamnn equation reads:

 

fa(X+eaAl,t+At):fa(x,,)_fa( ) fa ( )+ (a )
T

62 f;‘/(x,z) (56)

Note that the assumption that the momentum of fluid particles conserves at each collision
[150-152] has been applied in derivation ofthe above equation.

In summary. Table 1 lists the extemal force temr F and the equilibrium velocity

(1

u* in different models. The fluid macroscopic velocity u has the same form as the
equilibrium velocity in all the five models and thus has not been listed in the table.

In order to select an accurate external force model for simulation of plasmas, the
above described five models have been applied to simulate the same physical problems
and the results are compared. In all the simulations, the external force term takes the form

listed in Table 1. The discretized equilibrium distribution function follows Eq. (53) where

46

the equilibrium velocity is evaluated according to the form shown in Table 1.

Correspondingly. the lattice Boltzmann equation shown in Eq. (52) is used.

 

 

 

 

 

 

 

 

Model Fa u*
(one -a 1
Method 1 —“—,a— ‘29afa
0~ n a
_ — * . * 1
MethodZ (outta. (ea u )+(ea u )ea ‘Zeafa
62 94 It a
1 (can lZe f +—]—Ata
NICIllOd 3 I“; 07 ( (I ) I] a a a 2
-' 1 ‘ _ It _. . IF 1 1
Guo’s Model [l—— mana- (6“ )u )+(e” u )8“ —Zeafa+—Ata
2r ,9- 64 n a 2
a. e —u , 1
HE‘S lV’lOClCl —(—0aT——)—féq(x,l) ZZeafa
a

 

 

 

 

 

Table 1 Fa and u'" in different external force models

The first simulation example is the Poiseuille flow between two parallel placed
plates which are located at y = O and y = 2H respectively. The analytical solution ofthe

fluid velocity is:

u"‘ {inn—g). 0] (57)

,or'
where ,0 is the mass density of the fluid, v is the kinematic viscosity which can be

determined by Eq. (48) and g = —dp/dx is the constant pressure gradient from the inlet

47

to the outlet. In this simulation. g is treated as the external force term in the lattice

Boltzmann equation and thus, the acceleration term in Fa can be written as a = [5. 0 .
,o

The simulation was taken on a 256><32 grid system. The length and height ofthe
computational domain are 1.0 and 0.125 respectively. The lattice speed C is set as 1.0
and the dimensionless relaxation time 2' is 0.6. No gradient boundary condition was
applied at the inlet and the outlet while no-slip boundary condition was implemented on
the top and bottom walls. The fluid is initially stationary and is accelerated by the
pressure gradient when 1 >0. The calculation is terminated ifthe flow reaches its steady

state. The convergence criterion is:

 

~n+l - - "I ' '
(1.1.7 a (13].)

where the superscript n+1 denotes the values at the next time step. Different values of g
were used in the calculation as long as the low Mach number approximation is retained.
The relative errors between the simulation results and the analytical solution are listed in

Table 2 where the error is defined as:

 

Jinn-58028.1)—tt_‘\.4(.rc.,jA_t‘)]2
error = I (59)

\[Z [113.4 (.rc. , jAv)]2
j

 

 

[,3 . . . . . . ,.
where 1!". represents the x—component velocrty obtained from the srmulatron while it:
is the analytical solution from Eq. (57). The error is only calculated at the center cross
section along the x-axis (i=128) because the flow is in the steady state and it, is not

dependent upon -r. The relative errors listed in Table 2 show that. although the

48

performances of Method 3 and Guo’s model are slightly better, the five external force
models almost have the same accuracy for simulation of Poiseuille flow which is driven
by a constant external force. Without surprise, the error grows with the increase of the
magnitude of the force because a larger Mach number is resulted by a bigger external

force. Figure 4 shows the profile of x-component velocity obtained by He’s model with
g = 5.774x10‘2 . It is apparent that even with the largest error in Table 2. the simulation

result is still in a good agreement with the analytical solution.

 

 

 

 

 

 

 

 

 

 

 

 

g 1x103)
1.155 2.309 3.464 4.619 5.774
Models

Method 1 4.499 4.888 6.314 8.414 11.043
Method 2 4.500 4.888 6.315 8.416 11.046
Method 3 3.672 3.931 5.116 7.057 9.596
Guo’s Moder 3.672 3.931 5.116 7.058 9.597
He‘s Model 4.500 4.888 6.315 8.416 11.046

 

Table 2 Relative errors (x104) of Poiseuille flow

49

 

 

       
    

 

 

 

0.53 7
0,44
’9.
370.1%
:3
o He'sModel
0'2 7 Analytical Solution
0.14
0.04:1) ' r ' r ' r . 1 ' kf
0.000 0.025 0.050 0.075 0.100 0.125

Y

. . . -2
F rgure 4 Proﬁle of x-component velocity obtained by He’s model (g = 5.774X10 ).
The simulation results agree very well with the analytical solution.
As the second simulation example, an unsteady 2D Taylor vortex flow was
simulated by all the five external force models aforementioned. The simulation was
carried on a 129><129 grid system with the dimensions of the computation domain of

—77 S .t‘,_l‘ S Jr. The space- and time-dependent acceleration due the external force is:

P" _

 

 

 

7
It 11“ .
— 1 0 srn(2/q.r)exp[—2|’(/\'12 +/\’,7?)’I
a = A2- 2 ((70)
' 11 .
——l—Q3111(2/(2y)(“pl—214412 + 4'22 III
_ ._ 1'2 J

where no is the magnitude of the flow velocity and it is set as 113 = 0.001 so that the

compressibility of the flow can be neglected, v is the fluid viscosity, k1 and k2 are the

wave numbers and k1 2 k2 21.0 in the calculation. With the external force shown in Eq.

(60). the analytical solution ofthe velocity field can be found as [138]:

u- ’“0 cosfkr-Y)Sin(k2.\’)€XPI-V(/<12 + 4% )rr

— no {1 sin(k1x) cos(k2y) exp[—v(kl2 + k7? )t]
2

In the practical simulation, the time step At is set equal to the spacing step Ar. Then the
dimensionless relaxation time I can be retrieved from Eq. (48) with the value of
v = 0.005 . The initial velocity field is set as Eq. (61) with t: 0. The periodic boundary

condition [153] is applied for all the four boundaries ofthe computational domain.

 

 

 

 

 

8
4
VA
0 a
E: l .l ['1‘ If.
C ;:'. t. .. '1 1
9
L. tlI,["|.|“i‘r
Lg \IC‘l'TtIi'i'T'At/'r’i‘i“\>""""l\'\.!’
. Z 3 _ 319:1“??th ﬁdﬁkw’iﬂﬁgﬁ—ﬁ $33.56" 35:33.51
*5 q — Methodt
g 2_ ---- Meth0d2
. ---'—- Method3
1 - 9 (300's Model
4 0 He's Model
0 ' 7 ' 1 ' 1 ‘ l ' l
0 300 600 900 1200 1500
Time Step

Figure 5 Evolution of relative errors of simulation of Taylor vortex ﬁow obtained by
the ﬁve external force models. All the errors are in the same order of values.

51

Figure 5 shows the relative errors of the simulation versus the time steps where

the error is defined as:

where u U”

 

JZNLBM) — U"'(1'A\3./'A.rllz

t'./'
JZUAII'AY../Ar)2

i. j

 

C’I‘I'OI' =

 

(1.1) represents the velocity at the node point of (i,_/') obtained by the LBM

and ll‘ii(1'A\‘.‘/'A.l') is the analytical solution according to Eq. (61) at the corresponding

position (iAt‘,jAr'). It is clear from Figure 5 that all the relative errors are in the same

order ofrnagnitude. Although Method 3 and Guo’s model present better performance in

Ux(0,y)/U0

Figure 6 Proﬁles of x-component velocity along the center line of x-axis

 

 

 

 

 

0.8
1 Analytical Solutions
06 - 0 He's Model at 1500M
1 ,I \ A He's Model at 3000M
0.4'7 / "3'1”": \
l \
‘ / ( ‘ \
0.23 /’ ‘\
.1 {7 x}
0.0 J L, ,
" i‘:\\ 7";
-0.2 - X . x,"
4 \ i‘ Jill]
x ~, .. /
'O.4-l \\ if a, .L. /’
-o.e~ i
-0'8'1'1'171’1'1'1'1'1‘
-1.0-0.8-0.6 -0.4-0.2 0.0 0.2 0.4 0.6 0.8 1.0

y/rt

at time steps of 1500 and 3000 by using He’s model. The simulation results are in
good agreement with the analytical solutions.

simulate temporally and spatially varying force as illustrated in Figure 5, other methods
also can give highly satisfactory results. For example, Figure 6 shows the profiles of x-
component velocity along the center line of x-axis by using He’s model at different time
moments. Comparing with the analytical solutions which are represented by the solid
lines in the plot, He’s model can reproduce excellent results with respect to the analytical
ones.

Finally. let‘s check the external force models by simulating a simple free
diffusion problem with a constant force pointing to the positive x-direction. Initially. the
fluid density has a Gaussian distribution centered in the computational domain with the

maximum perturbation of the density as 0.01/70. With a constant external force acting

along the positive x-direction, the analytical solution ofthe ﬂuid density is:

2 2
0001 1‘-.. .— ’ I ‘0_ 3.
p'l(x,t».t)=p0+ p0 ex _I‘ ’Q M1) +(.1 .1.)

—-——.— (93)
(1+t/lt0) I'2(I+I/[0)

 

where (.\‘(.,_l-’C.) represents the center point of the computational domain, r=20Ar and
10:7'2/41) where D is the diffusivity which can be recovered from the lattice
Boltzmann equation as:

D=(r——0.5)192Al (64)
In Eq. (63), 14,, is the magnitude ofthe drift velocity vd due to the external force, which

is calculated as follows [154]:

vd = —a 7A! (65)

53

The simulation was taken on a 129><l29 grid system with the dimension of the
computational as 1.0. The time step is set equal to the Spacing and consequently the
lattice speed c: = 1.0. Different amplitudes ofthe acceleration were used in the calculation.

Table 3 lists the relative errors ofthe diffusion problem with different amplitudes
of the external force. The relative error in this case has a similar definition as shown in

Eq. (62) where the velocity is replaced by the fluid density:

 

ZIP/“B (I) j) — pA (1413143)]2
1.] (66)

t/Zp‘tz'mut'Ayiz

 

error =
Iii
We can see from the table that the errors increase with the enlargement of the external
force as expected. However. Method 3 and Guo’s model show larger errors with respect
to Method 1, Method 2 and He’s model.
Figure 7 shows the fluid density distribution at 200 time steps with different

amplitudes ofthe external force. The solid lines represent the analytical solutions. It can

 

 

 

 

 

 

a" 0.5 2.0 4.0 6.0 8.0 10.0

Models
Method 1 1 223 1.240 1.297 1.389 1.515 1.673
Method 2 1 223 1.250 1.335 1.471 1.650 1.864
Method 3 1.261 1.755 2.799 3.968 5.172 6.385
Guo’s Model 1.271 1.862 3.060 4.371 5.705 7.033
He‘s Moder 1.223 1.251 1.336 1.469 1.641 1.840

 

 

 

 

 

 

 

 

 

Table 3 Relative errors (X105) of free diffusion problem
with different amplitudes of the external force

54

1.0030

 

 

 

 

 

 

 

 

 

 

 

 

. Analytical Solution
1.0025 a 0 He's Model
] A Guo's Model ‘
1.0020 ~ I”
O J ‘\
3 . J \
9 "
3 1.0015— \\
9 . J .x
.\
1 .001 O 4 .7, "ﬂ
1.0005~
1.0000 1:11;”,‘55’ v v v v v 1 1 4.1.1:]
-0.5-0.4-0.3-0.2-0.1 0.0 0.1 0.2 0.3 0.4 0.5
X
(a)
1.0030
Analytical Solution
1.00254 0 He's Model
. A Guo’s Model
0 1.0020—
:3
LO- .4
E 10015
1.0010~
1.0005a
1.0000 "f~;’.1.'114;14:’ I. v V r ' Y 1 r j
-0.5-0.4-0.3-0.2-0.1 0.0 0.1 0.2 0.3 0.4 0.5
X
(b)

Figure 7 Fluid density distribution obtained by He’s model and Guo’s model at 200

time steps with different magnitudes of the external force: (a) ax = 4.0 (b) ax = 10.0.
Both Guo’s model and He’s model can produce good results if the external force is
small (as shown in (a)). However, if the external force is relatively large, Guo’s
model induces larger error than dose He’s model (as shown in (b)).

55

be seen that the Guo’s model can give relatively good results ifthe external force is small

( Figure 7 (a)). However, if the external force is large to some level, e.g. the case of

(1, =10.0. Guo’s model cannot reproduce the satisfactory results, as shown in Figure 7

(b). Meanwhile, as illustrated by Table 3 and Figure 7, the results obtained by He’s model
have good agreement with the analytical solutions even with a relatively large external
force.

In conclusion, although Method 3 and Guo’s model behave better in the
simulations of Poiseuille flow and Taylor vortex flow, they show poor performance in the
simulation free diffusion problem, especially when the amplitude ofthe external force is
large. Method I possess higher accuracy for Poiseuille flow and the diffusion problem,
but it gives large errors in simulation of Taylor vortex flow. Method 2 and He’s model
have the same accuracies for the Poiseuille flow and the Taylor vortex ﬂow, but He’s
model is better in diffusion problem. Comprehensively considering the above observation
from the results obtained by the five external force models, He’s model is chosen in this

research for all the simulations thereafter.

3.2 LBM Simulation of Magnetohydrodynamics Flows

By using He’s model to evaluate the external force temr. a hybrid lattice-
Boltzmann finite-difference method is proposed in this research to simulate
incompressible, resistive magnetohydrodynamic (MHD) flows. In recent years. many
attempts have been made to develop LBM algorithms for MHD problems. Chen et al.
[155] and Martinez et a]. [156] employed the bidirectional streaming for 2-D MHD
problems where the distribution is propagated into two different directions associated

with the velocity and magnetic ﬁeld. The former used 37 discrete velocities while the

56

latter reduced it to 13. Schaffenberger and Hanslrneier [157] later reduced the number of
velocities even further to nine by employing the standard streaming rule on a 2-D square
lattice. Dellar [158] developed a new method, where two distribution functions are
utilized to represent the hydrodynamic momentum and the magnetic induction. The
hydrodynamic part is simulated using the conventional low Mach number LBM, and the
magnetic field is represented by a separate vector-valued magnetic distribution function,
which obeys the vector Boltzmann-BGK equation. This method has been later extended
to 3-D by Breyiannis and Valougeorgis [144].

In all the afbre-mentioned methods. both the magnetic induction problem and the
flow problem are dealt with by a lattice kinetic approach. However, while employing
LBM for the Boltzmann equation is natural, the use ofa kinetic approach for solving the
magnetic induction equation is not quite intuitive because, after all, the Boltzmann
equation and the magnetic induction equation constitute a set ofgoverning equations for
MHD. In other words, a lattice kinetic approach does not need to be used for the
magnetic induction problem even though it is a more consistent approach and has many
advantages in many cases. In fact, other numerical methods, such as the finite difference
method. can be easily employed to solve the magnetic induction equation with equal or
better accuracy because those methods are well established.

In this proposed hybrid model. the lattice Boltzmann equation with the Lorentz
force term is solved to update the flow field while the magnetic induction equation is
solved using the finite difference method to calculate the magnetic field. This approach is
methodologically intuitive because the governing equations for MHD are solved in their

respective original forms. In addition, the extension to 3D is straightforward. To validate

57

this newly developed model, the Hartmann flow, the Orszag-Tang vortex system (2D and
3D) and the magnetic reconnection driven by doubly periodic coalescence instability are
simulated. The obtained results agree well with analytical solutions and simulation results
available in the literature.

In the present model, the governing equation is the BGK-type Boltzmann
equation as shown in Eq. (8), where the equilibrium distribution function f” is

described by the Maxwell-Boltzmann distribution as follows:

_ 2
f”l -——’2——exp _____(v u) (67)

_ 2 o2 ‘ 2
(2776’ ) 1‘ 20
where ,o and D are mass density and dimension ofspace. respectively. For MHD flows.

the acceleration a can be written as:

a=-—l—(V><B)><B ((78)
p’u

Here, B is magnetic induction and ,u is magnetic permeability. The evolution of
magnetic field is obtained by solving the magnetic induction equation, which is derived

from Maxwell’s equations with the assumption of 6E/6t = 0 [129]:

,
i—thyV2B+(B-V)u-(u-V)B (69)
CI

where 17 is magnetic diffusivity and is expressed as 1] 20101—1, where o is the
electrical conductivity.

For 2D problems, the Boltzmann equation can be discretized by D2Q9 scheme as
described in Chapter 2. The external force term is evaluated by He’s model. Then the

resulted lattice Boltzmann equation follows the form of Eq. (52) where the discretized

58

equilibrium distribution function takes the form Eq. (32). The mass density and

momentum then can be calculated through:

p=Zfa (70)
pu=Zeafa (71)

The magnetic induction equation (Eq. (69)) is solved by the conventional finite
difference method. In two dimensions, for example, the x-component of Eq. (69) is

written as:

 

 

 

 

 

.. 7 .,7 ’ .,
C‘B . 1 5‘8 0‘8. on . 811 . GB. 58.
" = —- 2x + 2" + Bx “ + 8,. " -—ux ﬁ " —ll,v. —" (72)
(‘1 yo m a}, at " ﬂy or ’ 6y

Then the following discretized equation is obtained if the central difference scheme is

employed.

ail/3.x >19]. _ 1 (8.x‘):i+1yj V218.\‘),nj +(Bx):’_19/

 

 

 

 

 

 

(ff 1110’ ( At‘)2
If n H
+ (B-")1'./‘+l — 2(B.\‘),'_j +(B.r),~.j_1]
(At-)3 (73)
II n II II
+(B )n (”x )i+l._j _(u.t’ )i—1.j +(B )n (“xx-+1.]- ‘(llx )i.j-—1
I 1" I 2 Ar y 1.]. 2A).)
)1 n N N
-(11 )n (Bx )t'+l.j — (Bx )i—l.j -(u )n (Bx ),"j+| —(8.r by-)
.t‘ 1' , j- 2 Ar '1‘ in]. 2A)?

Note that all the terms on the right hand side are evaluated at time step it using the fluid

velocity calculated by the lattice Boltzmann solver. Eq. (73) can be re—written as follows:

——’l= R((B,.)”) (74)

59

In this study, the time derivative on the left hand side is discretized by the second-order

Runge-Kutta method:

1 t A I’
(B )’+1/2=(B_\.)’+—2:R((B,.)’)

.I'

(75)
(BX)H+I : (Bx)" + A!R((Bx)”+l/2)

For validation purposes, the 2D Hartmann flow is simulated because its analytical
solution can be easily obtained. Hartman flow is a channel flow induced by a uniform

magnetic ﬁeld ( 30) applied perpendicular to the flow direction. For 2D Hartmann flows,
velocity has only one component in the channel direction u=(u,..0) and this flow

induces additional magnetic field in the flow direction. Therefore, magnetic field can be

written as B=(B',.,BO). The magnetic induction equation and the hydrodynamic
momentum equation can be simplified as follows:

2
i B. *1 .
tio dy“ d}:

 

 

2
v d 11,. +531 (1'83. :

77
2 # (1y ( )

(11'

(I) . . .
where 51:7]— rs the constant pressure gradrent used to drive the flow along the x-
ix

direction. By using the following boundary conditions

 

11,. =0 at y=iL
' (78)
BI = 0 at y = i].
the analytical solutions to Eqs. (76) and (77) are obtained as:
0' . r '
(7,1)) =—g—Lcoth(M) “WM/L) —1 (79)
' 03g cosh(M)

6O

 

1L sinhMtx-ll. 1‘
B.\‘(.\,):_gr [ ( ., )_,_] (80)

30 sinh(M) L
where M = BOL o/pv is the Hartmann number. If there is no external magnetic field,
the Hartmann number becomes zero and the flow reduces to the Poiseuille flow, which
. , 2 2

has solutrons of it_,.(_1') =g(L -,v )/2,ov and 3, =0.

For the LBM simulation of the Hartmann flow, the initial distribution function is
set to the equilibrium distribution function with constant density, p =1.0; 11,. and B, are
set to zero. A uniform magnetic field BO is applied in the _t~'-direction. The bounce-back

condition is used for wall boundaries while the periodic boundary condition is used at the
inlet and the outlet. The simulation is temiinated until a steady state is reached. To drive

the flow, the driving pressure gx is added to the external force as follows:

paz—l—(VXB)XB+géx (81)
It!

In Eq. (81). e, is the unit vector in the x-direction.

Figure 8 and Figure 9 present simulation results of it, and Bx respectively,
where solid lines denote analytical solutions given by Eqs. (79) and (80). As clearly seen,
the simulation results agree well with the analytical solutions for all Hartmann numbers
considered. In this simulation, a 200><20 grid is used and the parameters of r =0.635,
o :10 and g = 2.5><10‘5 are adopted. The Hartmann numbers of 0. l. 2, 5, 10 and 20
are considered. Note that the Hartmann number can be changed by varying the magnitude

of the applied magnetic field 80. In Figure 8, it is shown that as the Hartmann number

increases the velocity profile becomes flatter, which can be explained by Eq. (80).

61

0.03 .,.,..

 

 

 

 

 

0.02-
X
:7
0.014
0.00 I ' I f I ' l ‘ i
-1.0 -0.5 0.0 0.5 1.0
y/L

Figure 8 Proﬁles of x-component velocity with different Hartmann numbers: M=0
(squares), M=l (circles), M=2 (upper triangles), M=5 (lower triangles), M=10
(diamonds), M=20 (stars). The solid lines show the analytical solutions. The
simulation results agree very well with the analytical solutions.

0.0004 . . + r . .

 

0.0002 —

 

 

 

 

83 0.0000
000024
00004 2 , . A . , .
-10 -0.5 0.0 0.5 1.0
y/L

Figure 9 Proﬁles ofx-component magnetic ﬁeld for different Hartmann numbers:
M=0 (squares), M=l (circles), M=2 (upper triangles), M=5 (lower triangles), M=10
(diamonds), M=20 (stars). The solid lines show the analytical solutions. The
simulation results agree very well with the analytical solutions.

62

Clearly, the applied magnetic field (Bo) tends to reduce the magnitude ofx-component

velocity. Ifthe Hartmann number is finite, the velocity profile cannot be entirely flat even
for very large Hartmann numbers and there must be a region with a large velocity
gradient because velocity is forced to be zero due to the no slip boundary conditions.
Figure 9 presents the magnetic field induced by the flow for a variety of Hartmann
numbers. the simulation results agree with the analytical solutions very well.

As the second test problem, the Orszag-Tang vortex system (which is an unsteady,
nonlinear MHD flow problem) was chosen. Since Orszag and Tang [159] ﬁrst studied
this problem, it has become a popular benchmark problem because many aspects of MHD
turbulent flows appear in this problem, such as the dynamic alignment, selective decay
and magnetic reconnection [160]. In this study, a 2D Orszag-Tang vortex problem is
simulated with the following simple nonrandom detemtinistic initial conditions:

uo = -—no(sin y, —sin x)

. . (82)
Bo = —Bo(srn _v.-—srn(2x))

where no 2 2.0 and Bo = 2.0 are the initial velocity and magnetic induction, respectively.

The simulation is performed on a square domain of ()Sx._1'S 2,7. A 512><512 uniform
grid is used and the periodic boundary conditions are applied on all boundaries. Kinetic
fluid viscosity and magnetic diffusivity are assumed to be the same (v = 77 = 0.02 ), which
leads to the same Reynolds number (Re) and magnetic Reynolds number (Rem) at the
initial stage. With the initial conditions shown in Eq. (82), the evolution ofthe vorticity
of fluid (a) = Vx u ) and the current density (j = Vx B/tio ) are demonstrated in Figure 10
and Figure l 1, respectively. As shown in Figure 10 and Figure l I. initially both velocity

field and magnetic field have symmetric structures. As time elapses. the initial flow

63

 

 

.\‘\‘>
4“ \-
. \. s,

 

 

 

 

 

 

 

 

 

 

t= 0.738

Figure 10 Evolution of vorticity for the 2-D Orszag-Tang vortex. The ﬂow pattern
becomes complicated due to the nonlinear interactions between the velocity ﬁeld
and the magnetic ﬁeld. A ﬂat quadrupole-like conﬁguration emerges with time lapse.

64

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

t= 0.738 t= 1.00

Figure 11 Evolution of current density for the 2-D Orszag-Tang vortex. The existing
current sheet at the center of the ﬁgure is enhanced and eventually, a thin elliptic
structure establishes due to the magnetic reconnection occurring there.

65

pattern becomes complicated due to the nonlinear interactions between the velocity ﬁeld
and the magnetic ﬁeld. In Figure l I, the existing current sheet at the center ofthe ﬁgure
is enhanced and eventually, a thin elliptic structure establishes due to the magnetic
reconnection occurring there. At the same time, a region of sheared flow coexists with the
current sheet. which is shown as the flat quadrupole-like configuration in Figure 10. The
contours of the vorticity and the current density agree well qualitatively with the
simulation results available in the literature [160].

In order to validate this model quantitatively, this problem was simulated by using
the same parameters used in the paper by Dellar [158] (no = 2.0, v =27 =0.02, and a
512><512 grid). The maximum vorticity value at t = 0.5 predicted by this model is 6.764
while Dellar’s result gives a values of 6.758. At t = 1.0, the values are 14.457 and 14.20
respectively. The maximum current densities are also compared and the values are given

in Table 4. In order to check if these two models satisfy the divergence-free property of
the magnetic ﬁeld, the values of [V-Bl are calculated. The presented model gives
0.00463 at t = 0.5 and 0.00922 at t = 1.0 while Dellar’s results indicate 0.0062 and 0.0415

at the corresponding times. Therefore, the presented model is validated not only

qualitatively, but also quantitatively.

66

 

 

 

 

 

Time (sec) Max vorticity Max current IV-Bl
0.50 6.764 18.129 0.00463
Present results
1.0 14.457 45.963 0.00922
0.50 6.758 18.24 0.0062
Dellar’s results
1.0 14.20 46.59 0.0415

 

 

 

 

 

 

 

Table 4 Quantitative comparison of simulation results with Dellar’s work

To show that this model can be easily extended to 3D, the D3Q19 lattice model
[161] has been employed to solve the 3D Orszag-Tang vortex problem with the following
initial conditions [162]:

uo = (—251ny,25inx,0) (83)
BO = 0,8(——2 sin(2y) + sin 2, 2 sin x + sin 2, sin x + sin y)

A cubic domain of 27rx 227x 277 is used and simulations are conducted on a 64x64x64
grid. Periodic boundary conditions are used for all boundary surfaces and edges. Figure

12 (a-b) shows the initial contours ofthe magnitudes of vorticity (Idol) and current density
(ljl ). In Figure 12 (c-d), [or] and [j] at t = 0.598 are shown and the slices of [to] and [j] at

z = 0 are presented. The patterns of the current sheet and the corresponding vorticity
resemble Figure 11 and Figure 10 closely because initially, the variations of current
density and vorticity in :-directions are relatively smaller in values than those in the x-
and y-direction.

Magnetic reconnection is the third example simulated by the present model. It is
the process where magnetic ﬁeld lines from different magnetic domains merge into one

another, changing the overall topology of the magnetic ﬁeld. Meanwhile, the stored

67

magnetic energy is released in heat and kinetic energy forms. The magnetic reconnection
can be driven by different forms of coalescence instability, for example, by the merging
ofa chain ofmagnetic islands [163] or by doubly periodic coalescence instability [164].
In this study, a magnetic reconnection problem driven by the doubly periodic coalescence
instability is simulated using the present model. As the initial distribution of magnetic

fltrx, a checkerboard pattern [165] represented by Eq. (84) is employed.
(you. ,t’) = Bo sin (77(x + y))sin (77(x — ,r')) (84)

The symmetric initial perturbation ofthe kinetic stream function is:
. . _ ,2 2
€00(-‘~.‘)—“0 exp(—IO(.\ +y )) (85)

where Bo = 0.5/rt, no = 0.05 and the initial magnetic and velocity ﬁelds are obtained as

follows:

6 6/

Bo-e~XVV/0=(--ﬂ)—, (”)1
6y 6
a 6 (86)
(PO €00

u —e-xV — —,

0 (P ( 5‘ 8x)

68

 

      

O 7 o
y/pl ' 1

(c) [to] at t = 0.598 with a 2D contour at z=0 (d) [j] at t = 0.598 with a 2-D contour at z=0

Figure 12 Isa-surface contours of magnitudes of vorticity and current density at t=0
and t=0.598 sec for the 3-D Orszag-Tang vortex. The patterns of the current sheet
and the corresponding vorticity resemble Figure 1] and Figure 10 closely because

initially, the variations of current density and vorticity in z-directions are relatively

smaller in values than those in the x-and y-direction.

69

The simulation is conducted on a square domain of —l Sx,y S1 and a 256><256

grid is used. Periodic boundary conditions are used on all boundaries. A fixed viscosity of

v=4><10'3 is used; ﬁve values of magnetic diffusivity (77:05le3 0.8><10_3 ,

w 1 1
3

1x10”. 2x10“ and 4x107") are considered. Figure 13 presents the evolution of

magnetic flux at different times with the magnetic diffusivity of 1x1073. Note that the
position of the magnetic current sheets does not change with time because the initial
perturbation of the kinetic stream function has a mirror symmetry with respect to x and y
directions. It can be seen from Figure 13 that the magnetic islands with currents ofthe
same sign move towards each other. The two comers coalesce into one and the original
two square cells become two adjacent pentagons. A current sheet fomis between the two
cells and the intensity of the current sheet increases. Eventually, the neighboring square
cells merge together. simplifying the topology structure of the magnetic ﬁeld to four
square-like islands. Figure 14. presenting the dependence ofmaximum current density on
magnetic diffusivity, is a quantitative evidence of the present model. As seen from the

-I/2

plot, the temporal maximum of current density can be approximated as jmax at I] as

illustrated by the dashed line in Figure 14, which can be compared with Fig. 3 in [165].

70

 

 

 

 

 

 

 

 

 

t=3.27
. + t
. l
i
;~ ‘1
r «2
. t
.[

t=5.10
' ‘ (l
'l
7.1
l
l

t=7.53

 

 

 

 

t = 4.49

 

 

 

 

 

t49.75

Figure 13 Evolution of magnetic ﬂux function for doubly periodic coalescence
instability. The position of the magnetic current sheets does not change with time.
The magnetic islands move towards each other. The two corners coalesce into one
and the original two square cells become two adjacent pentagons. A current sheet

forms between the two cells and the intensity of the current sheet increases.
Eventually, the neighboring square cells merge together, simplifying the topology
structure of the magnetic ﬁeld to four square-like islands.

 

5O

4 1
.\ TI
\
.\
30- \\ ~
\.
X \
E "
0\
ﬁ 20‘ \\\ "1
~““
-““‘~‘-~—‘.‘
10.4 _.
0

 

 

 

0.000 0.001 0.002 0.003 0.004

Tl

Figure 14 Maximum value of current density vs. magnetic diffusivity. The temporal

maximum of current density can be approximated as jmaX oc 17-le .

In conclusion. three classic problems in MHD flows were solved in this study and
the obtained results agreed well with the data available in the literature. This newly
developed hybrid LBM model has the following advantages: ﬁrst, its implementation is
relatively simple compared to other LBMs and Navier-Stokes equation based methods;
second, the extension to 3-D is straightforward. We believe that this approach can be a
good alternative to other MHD-LBMs that are fully based on the lattice kinetic

algorithms.

Chapter 4

Lattice Boltzmann Simulation of Isothermal Weakly Ionized Plasmas

As the governing equation for all transport phenomena, the Boltzmann transport
equation describes the evolution ofthe distribution function of each species of particles in
the plasma. All macroscopic variables of the plasma, such as number density and
macroscopic velocity, can be retrieved through proper moments of the distribution
function. Many achievements of LBM development (especially multi-component models
[166-170]) can be inherited in the plasma simulation, since plasmas are mixtures of
different types of particles. Among those models, the ﬁnite difference lattice Boltzmann
(FDLB) models [168, 169] can be used for asymmetric system, which is the system in
which the compositive particles have different properties. However, the direct use ofthe
FDLB models for plasma simulation is not sufficient due to some exclusive
characteristics of plasmas. For example. if the physical plasma parameters are used, a
very large dimensionless relaxation time will result. This relaxation time will
significantly reduce the effects of the collision term on the evolution of the distribution
function and thus lead to an ill-favored transport behavior of the electrons. Therefore, to
overcome this problem, the LBM simulation parameters should be selected in such a way
that the dimensionless relaxation time is in a valid value range. In this chapter, a new
LBM model is proposed to simulate weakly ionized isothermal helium plasmas. The
simulation parameters used in LBM are selected according to a rescaling scheme

developed in this research to match the physical properties ofthe plasma. For the purpose

73

of validation, the electron diffusion problem and the electrostatic wave phenomena are
simulated by this new model and the results show good agreement with the analytical

solutions.

4.1 Mathematical Model
The proposed LBM model in this chapter is based upon the following
assumptions:
1) Inelastic collisions. such as ionization and recombination. are not considered.
2) The plasma is isothermal. but different species can have different temperatures.
3) The plasma consists of electrons. neutrals, and singly ionized ions (three species).
Then the Boltzmann transport equation for the three species of particles in a

weakly ionized helium plasma can be written:

 

8 . c3 .
—{; + vs 'st +35 'V\'(.f.s‘ :( (is) (87)
C’ ‘ 0’ coll

The subscript 5' denotes the type of particles and can take 6, i, and n for electrons, ions.

and neutrals. respectively. In this equation, v5 is the microscopic velocity,
1; = ﬁ(x,v,t) is the number density distribution function, and as. is the acceleration due

to the Lorentz force. which is expressed as

E
as : (lb (88)

Ill

 

.S'
if electrostatic behaviors ofplasrnas are considered. Here, E is electric ﬁeld, q, and ”’3-

are charge and mass of species 3 . respectively.

74

By only considering the binary collisions in the plasma and applying the similar
splitting technique adopted for binary gas mixture [166], the collision term for species 5'

can be written as:

[6A] =JSB+JMI+JSH (89)
0’ coll

where J5", J“. and J” are the terms that represent the collisions with electrons, ions,
and neutrals, respectively. It is well known that if the plasma is weakly ionized, the
elastic collision with neutral particles is the dominant collision mechanism for all species.
Therefore. J3“ and J"”are negligible in the weakly—ionized plasmas. For the collisions
with neutral particles. the Bhatnagar-Gross-Krook (BGK) model is used, which assumes
that the particles relax to their equilibrium states during the characteristic time period,
which is called the relaxation time 2.5. Then, the Boltzmann equations for electrons, ions,

and neutrals can be written as

67, f. - ff”
:5 + V. W + 3.) 'Vv..fe = --—“—”- (90)
CI it’ll
c7- fr - f-""
:’-+v.~-Vf.-+a.--V.,.f.~=-——.—’3— (91)
(xi m
a. ("I
(,1 —
41+ vn 'an = _ f”) fun (92)
('1

”NH

5

where do", Am, 2”,, are the relaxation times for electron-neutral. ion-neutral, and

fa], fﬁ,’ are the equilibrium distribution

neutral-neutral collisions. respectively; ft)" m

(’11 ’ .

functions of electrons, ions. neutrals, respectively, due to the collisions with neutrals and

can be written as

 

C, n . (v . - u . )
fs'ii/ (usn) : .s 2 exp __s%_ (93)
7:05 26’;

Note that the collisions with the species other than neutral particles can be easily added in

the model. In Eq. (93), ()3. = kart/"1..- is the sound speed of species 5 where T, is the

temperature of species 5; um is the barycentric velocity of the binary collision with the

neutral particle:

’71 ll . '1' ”I ll
“5" = .S .5 N n (94)

”1.5 + ”In

 

where U, is macroscopic velocity of species 5 (u,,: macroscopic velocity of neutrals).
Note that u,,,, = u,, btrt u,,,, :t no and u," at “t- That 15 due to the fact that the frequency

of self-collisions between charged particles is very low in weakly ionized plasmas and
the charged particles cannot relax to their macroscopic velocity during the relaxation time

period ( 2,,” or xi,” ) when they collide with neutral particles.

Similar with the concept ofdensity dependent relaxation time [171], xi“, in Eqs.

(90) — (92) are written as [129]:

/1. =——1—— (95)

5n _’.
0.4-11”" <l s>

Here, 0' is the cross section ofthe elastic collision between species 8 and neutrals and

All

calculated as 0'," =7r(rS+r,,)2, where rs and r” are the radii of species 3 and n

. . . __ 8 k T.
respectrvely; <13.) rs the average speed ofspecres s, and (its) : [——§——-‘-

1,12
] [129]. Note
7: m

.S'

76

that the relaxation time presented in Eq. (95) is independent oftemperature because only
isothermal plasmas are considered in this chapter.
Now the Boltzmann equations shown in Eqs. (90) — (92) can be discretized by a

traditional .D2Q9 scheme and the external force term follows He’s model:

 

 

.- . l , Ata,-ec,'—u. ,
jj‘(x+effAt,r+At)=/f(x,t)———(ff(x,z)—f,f,;/~a(x,t))+ ‘- (a; 3) ,f‘l‘a(96)
(’11 (3
- . 1 -V .. Ata- (1-. .,
.t.”rx +4141.” 42) = /,A“(x.z) ——t/,“ (x.t) —./,-§,"“‘ 01.1)) + ’ (e; "’ ’ 1f” (97)
7711 (73
. 1 )
/,,C‘(x+ef,‘At,1+251)=f,§"(x,t)——T—-(f,f‘(x,t)-f,§,‘,/~a(x,t)) (98)
H)?

where At is time step; I r rm, are dimensionless relaxation times; superscript a

(311 7 iii 7

[1

denotes the a’ component in the phase space; eff is the am component of discretized

microscopic velocity of species 5 as deﬁned by Eq. (28). The discretized equilibrium
distribution functions foil! as well as the weights (0a are in similar forms deﬁned in

Chapter 2. Note that the self-collision distribution functions of the charged particles are
used in evaluation ofthe external force temi by using He’s model. This selection will be
justified later.

In multi-component LBM, if the same grid is used for all species, the time steps

for different species are all different because time step is determined by At =Ar/cls ,

where Ar is the spacing and (7,5. is the lattice speed of particle species 5. Using several

different time steps are very undesirable for a number of reasons. so in this model we use

a single time step based on the lattice speed of electrons. Then. during a time step At.

77

electrons travel a distance of efilAt to the neighboring node while ions and neutrals travel

a distance of eszt and eSAt , respectively. Since the lattice speeds of ions and neutrals

are signiﬁcantly smaller than that of electrons, the travel distance ofthose heavy particles
will be very small compared to the electron travel distance. As a result, if the same grid is
used for all species, ions and neutrals cannot reach the same node point as the electrons
do and an interpolation scheme [124, 172-175] needs to be used for heavier particles. In
this study, similar to the interpolation scheme used in [173, 174], a second-order
interpolation method is introduced to ﬁnd the on-node values of the discretized

distribution functions for ions and neutrals. In Figure 15, q and 0 denote two

neighboring nodes of [7 along the am direction (pointing from o to p). The ions

(neutrals) that are originally located at o, p, q arrive at 0', p' , (1’ after a streaming step.

The distribution function at p can be obtained by using [3(0'), ff(p') and [3((1') as

 

 

 

 

 

 

follows:
a t (1 I e? a .1 a t o I €6,112
e leo

where s can be e. 1' or 11 and

ff‘td) =f7f’to)
f." (p') = 7:11p) (100)
ftaftt') = 711(4)

78

Ax = effAt ef’émAt

 

 

 

Figure 15 Schematic of second order interpolation method for ion (neutral)
distribution function. The value of distribution function of the heavy particles atp
can be found by using the post-collision values at o’,p’ and q’.

where is“ is the post—collision distribution function after the collision step according to

the lattice Boltzmann equations shown in Eqs. (96) — (98). If we deﬁne
,6: eff /eff =,/7;.mL,/T(,ms , then:

 

 

 

 

[sail/3) = ('1- 431/;“(17'1—05/30 ~ﬂ)fsa(q')+ 0.530 + meaW) (101)

Note that if 5' =6, [3:1 and foa(p)=f,fz(o')=fca(o) which is just the result of the
ordinary streaming step.

Once distribution functions are updated, the number density and velocity of each

species and charge density can be obtained as follows:

125(x) = Eff“) (102)
n,(x)u,(x) = founfm (103)
mix) = etn,-(x)-n..(x)) (104)

The electric field E is updated by solving the following equation:

v3¢=—f’i (105)
50

79

where 8,, is electric permittivity of vacuum; (ti is the electric potential, and V¢ = —E. Eq.

(105) is solved by a Poisson solver.

4.2 Rescaling Scheme
The governing Boltzmann equation of the present model is discretized by the
traditional D2Q9 scheme. In this standard procedure, there are three parameters that can

be chosen freely: dimensionless relaxation time T, grid spacing Ar and time step At.
Then, the sound speed 63, in the D2Q9 model, is determined as 9 =Ar/\/3At and the
viscosity is recovered by the Chapman-Enskog expansion of the lattice Boltzmann
equation as follows [134, 137]:

(r—O.5)9Ax

73

Note that a tilde is used for the viscosity and the sound speed to distinguish themselves

(106)

 

t7=(r-0.5)é3At=

from the physical kinematic viscosity and the physical sound speed. We will call 17 and

~

(9 the lattice kinematic viscosity and the lattice sound speed respectively in this
dissertation. For some flow problems, such as driven cavity flow, the fluid flow is
characterized solely by the Reynolds number and the value of viscosity is not important
as long as the same flow pattern is obtained. In such a case, the lattice kinematic viscosity
need not be identical to the physical one, and the three lattice parameters can be selected
freely as just mentioned. In some cases, other dimensionless numbers, such as Peclet
number and Froude number, can be used to scale LB simulations. But for certain flow
problems, actual physical properties need to be used for the simulation because a
particular dimensionless number (such as Reynolds number) does not characterize the

flow andt’Or the flow depends strongly on some properties of a fluid. For example, in the

80

plasma diffusion problem. the Reynolds number is not important and the diffusion
characteristic is governed by the diffusivity. Therefore, the diffusivity (in this case
kinematic viscosity) must be matched explicitly. If the physical properties are used, the
number of parameters that can be chosen freely decreases due to the restrictions imposed
when matching the physical properties and the lattice properties. For example, only two
properties (grid spacing and time step) were selected freely in the fluctuation lattice
Boltzmann model [176], where the authors matched mass density, temperature, and
viscosity.

In the implementation of the parameter selection procedure based on physical
properties. some problems may arise. Let’s consider the following procedure to choose
the lattice parameters in standard LBM without an external force. First, grid spacing Av

can be determined from the spatial resolution requirement, and then time step is obtained

 

 

. . . k T
from At=Ar/\/30 , where (9 rs the physrcal sound speed (9: B After that,
m
dimensionless relaxation time I is calculated by equating {it with v:
3
.=1CZ+65 non
61M:
where the physical viscosity is
_ 2
v = 8? A (108)
.577

according to kinetic theory [126], and /i is the physical relaxation time, 2.. =l/on<v>

[126], where o is the collision cross section, and n is the number density of the fluid
particles. Apparently, only the grid spacing can be selected freely in this parameter

selection procedure by using physical properties of the fluid. Now, two problems are

81

expected depending on the value of t that is calculated from Eq. (107) because, as well
known, I must be in the range of (0.5.3.0) [147]. Firstly. if the viscosity is very small
(which is the case for most gases under normal conditions), then r is very close to 0.5

and the standard LBM schemes diverge. For example. for air at 1 atm and 300 K, v is

1.8245x10—51112/s and 6 is 293.36 m/s. Ifthe domain size is Im><lm and a 128><128

grid is used, the grid spacing Ar is 0.0078125 m. In this case, r is calculated as
0.500014 from Eq. (107), which causes LBM to diverge. Note that this lower limitation
of r is only applicable to the standard LBM. In the past 5-7 years, a newly developed
entropic LBM [177, 178] eliminates the instability caused by small values of r and the
algorithm is stable with arbitrary small viscosities. Secondly, ifthe viscosity is very large
(which is the case for a gas with low density and high temperature), LBM becomes very
inaccurate. For instance, for air at 0.00001 atm and 1500 K, v and 0 are 20.3991112/5
and 655.97 III/S from kinetic theory, respectively, which lead to r =7.394 with the same
grid described above. This value falls outside the valid range of r and generates a large
simulation error.

The problem of having too large 7 becomes aggravated in the simulation of
weakly-ionized plasmas because a low degree of ionization and the corresponding low
electron temperature of weakly-ionized plasmas lead to a large relaxation time and

subsequently result in a large value of I. For an example, an electron free diffusion

problem is simulated by using the LBM model described in the previous section. In this

simulation, a helium plasma with 1% degree ofionization in a 3.71X3.71 mm2 domain is

considered. and a 256x 256 grid is used. The number densities ofthrce species are set as

no =11, = 1016(11173) and 11,, = 101801173) . Then the temperature of electrons is calculated

as 0.8008 eV according to the Saha equation [129]:

E: e 3.00x1037TE3/3 iexp[—9Ti) (109)

II” "(3

where U, is the first ionization energy of helium (U, =24.59 eV) and T is electron

'6'

temperature in eV. By using the above physical parameters of the plasma, the
dimensionless relaxation time t is calculated as about 1.034x106 according to Eq. (107).
The initial distribution of the electron number density is in the following Gaussian
distribution:

._ 2 ,t— .’ 2
(.t x.) +(.1 .lc) (110)

2
r

 

11,,(x,t=0)=n,.o l+0.01exp —

where r = 0.290 mm and (.x,.,y,. ) represents the center point ofthe computational domain.

This free diffusion problem can be described by the Fick’s law:

7.
one _

“E— — DBNVZNC, (11])

where Do” is the diffusivity of electrons due to the elastic collisions between electron

and neutral and it is calculated as [154]:

 

 

8632-
Den 2 L (.11 (112)
37!
It is easy to ﬁnd the analytical solution of Eq. (111):
2 2
0.01 , .—-.'. r— ._
11,,(x,_v,t) 2 "e0 +————’,i‘iex -(Y i‘z) +0 '1‘) (113)
(1 +1770) r (i +t/to)

where to = 1'2/4Dm .

83

1.012

 

 

 

 

 

- - - - Initial Distribution
-—- Analytical Solution
1.010 n ,\ 0 Result by Using
T I, ‘\ Physcial Variables
1.008 ~
0
E
3, 1.006 -
Z
1.004 —
1.002 ~
1 7.;-., i7: . . "
1.000~—.—,‘*";'4-; "47"”, . I , ' ; ,"‘.—,—‘l
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5

X (mm)

Figure 16 Distribution of electron number density at t = 3.5 ns. Non-physical peaks
appear if the physical variables of electrons are used.

Figure 16 shows the distribution of electron number density at t = 3.5 ns obtained
form Eq. (I I3) and from the simulation. Clearly, by using the physical variables, the
result is very inaccurate and non-physical peaks appear in the solution. The appearance of
these non-physical peaks can be explained as follows. In the case of electrons, the large
non-dimensional relaxation time leads to a near-zero collision term in the lattice
Boltzmann equation, and therefore. the effect of collision step is almost negligible. In
other words, virtually only the streaming step is left in the implementation of LBM.

Figure 17 shows the electron number densities at two time moments with different
ionization degrees. It can be seen that the sub-peaks move with the constant lattice speeds
and the magnitudes do not change. We can also see that the sub-peaks of higher

ionization degree move faster than those oflower ionization degree because the sound

84

 

—— 1% Ionization Degree
1.005 n - - —- 2% Ionization Degree
—-—- 3% Ionization Degree

1.004%

1,003.

Ne/NeO

1.002~

1.0014

 

 

1.000~

 

 

 

— 1% Ionization Degree
1005 s - - - 2% Ionization Degree
—--- 3% Ionization Degree

Ne/NeO

 

 

 

 

Figure 17 Electron number density distribution at (a) t = 1.116 ns and (b) t = 3.347
ns with different ionization degrees. Non-physical sub-peaks caused by using the
physical variables travel with the constant lattice speed. A higher ionization degree
leads to a higher lattice speed. The ratio of the magnitude of the sub-peaks to that of
the primary peak corresponds to the value of the ratio of the weight coefficients.

85

speed is higher in the former case. The magnitudes of all three sub-peaks are also
checked: the ratio of the magnitude of a sub-peak to that of the primary peak is 1:4,
which is the ratio of weight coefficients in the D2Q9 mo‘del employed in this study.
Figure 18 is the contour plot of electron number density at t = 1.1 16 ns. The dashed lines
in the plot represent the discretized phase space and the particles stream to their
neighboring nodes along the directions denoted by the numbers. The numbers in the
brackets are the weight coefficients used in the discretization of the equilibrium
distribution function. Figure 18 clearly illustrates the evolution of the electron number
density along with the vector directions in the discretized phase space if the LBM

simulation was conducted without the collision step.

 

Y (mm)
I
.A “—i—Fl (A)

:5

tr

0 r
:9
”.2

Q

I" .

\
1' ii‘i' ’ i _”
I 'x ‘ ‘ ~ -

u
t
I

., ..__-.--"";.,c-..=
( l‘4(,1/9)\\V8(1/36)

 

 

 

X (mm) 3.71

Figure 18 Contour plot of electron number density at t = 1.116 ns. The distribution
of electrons is along with the discretized velocity directions if the effects of collisions
are not included in LBM.

86

In conclusion, the usage of physical properties of the ﬂuid will probably deduce
an invalid value of the dimensionless relaxation time and thus lead to large simulation
errors. To overcome this problem, a rescaling scheme is presented in this research that
can be used with physical properties of fluids for a wide range of particle number density
in the simulation of fluid flows and weakly-ionized plasmas. The key idea is to design a
scheme in such a way that the kinematic viscosity of the fluid and the characteristic
velocity due to the external force are matched in the implementation of LBM because,
after all, the Boltzmann transport equation solves convection diffusion problems.

First of all, to preserve the diffusion characteristic of the problem, it is assumed

that the lattice kinematic viscosity 17 (Eq. (48)) is identical to the real kinematic viscosity

v (Eq. (108)). Then the lattice sound speed 61~ can be calculated as:

 

 

2
a = 86 ’1 (114)
ﬁerxU—Oj)
From Eq. (1 14) the rescaling parameter 7 is obtained as:
a ﬂux: — 0.5)
y E —. = (115)
(9 86/1

Then the lattice sound speed is expressed as a function of the rescaling parameter as

follows:

67:9 (116)
7

Note that the rescaling parameter y is a function of Ax, /1 and 6, and therefore is ﬁxed
once the type and condition of fluid (1 and 61) and grid spacing (Axiare determined.

Next, the lattice relaxation time x1 is determined from the dimensionless relaxation time,

87

/1 = rd! , where the time step is At = Arc/£61. Using Eq. (1 15), the lattice relaxation time

2: is calculated as follows:

,i:_822/;__,1 (117)
371'(Z'—0.5)

Secondly, the convective property of the ﬂow when an external force exists must not be
altered by the LBM scheme. If the ﬂow is affected by an external force (such as an
electromagnetic force), depending on the nature of the force we may need another
rescaling to obtain accurate velocity fields. So, we propose the second rescaling rule: the
characteristic velocity of the flow should not be changed by the rescaling, i.e.,

u0 = 00 (1 18)
where U,) is the characteristic velocity of the flow. From this rule, we can obtain the

rescaled acceleration 5 , which we call the lattice acceleration in this article.

As mentioned earlier, dimensionless relaxation time t and grid spacing Ar are
the two parameters that can be chosen freely in this rescaling scheme. Once Ax is
determined, 1 can be selected considering the Mach number (Ma) and the Reynolds
number (Re) of the problem. The following equation is the relationship between the Re

and Ma obtained from the definition of Re:

__ x/gMa/V
r—0.5

Re

 

(119)

where N is the number of grid points along the characteristic length L0; Ma =U0/61,

where U0 is the magnitude ofthe characteristic velocity. From Eq. (1 19), we obtain:

JEMaN
Re

 

+0.5 (120)

88

which must lie in the valid range. Note that according to Hou et al. [135], there exists a
critical value of 1 below which simulation results show unphysical patterns or the code
diverges even though I is greater than 0.5. Eq. (120) is especially useful when the ﬂuid
ﬂow is characterized by the Reynolds number, such as driven cavity ﬂow. Given Ma and
Re, I can be determined from Eq. (120). Also, the Mach number of the problem must be
small enough to meet the low Ma requirement. In the mean time, as will be shown later,
LBM converges faster with a higher 2' .
In cases where Re is not important or is hard to be defined, the definition of Ma

can be used to choose I :

Ma: ﬁzzAxUO(r—0.5) (121)

8621

 

from which the following is obtained:

sazxMa

r=——+0.5 (122)
ﬁﬂ'AXUO

In the practical implementation of this rescaling scheme, the ﬁrst step is to
determine the grid spacing Ax and the dimensionless relaxation time 2' (following the

method described above) and then the rescaling parameter y according to Eq. (I 15) with

physical properties of a ﬂuid. Then the lattice sound speed 0 and the lattice acceleration
a are determined from Eqs. (1 l6) and (118), respectively. Once these lattice parameters
are obtained, the remaining steps are same as the standard LBM.

To validate the present rescaling scheme, the two-dimensional driven cavity ﬂow
is simulated with and without the rescaling. In this ﬂow problem, no external force

appears in the lattice Boltzmann equation, and therefore, only the first rescaling rule is

89

 

1.0

   

 

 

 

 

 

   
 

    

 

 

 

 

0.8 — -
Present Simulation
0’6 ’ Re = 100 -
Ej - - -- Re = 200
>- 1 , ----- Re = 400
0.4 t. J Data from Hou et al. (1995) -
;w 0 Re = 200
1' ‘ A Re = 400
A .
0.2 _ \. Data from Ghia et al. (1982) _
A‘ 0 Re = 100
, A Re = 400
0.0 . ' - 1 e 1
-O.5 0.0 0.5 1.0
Ux/UO (Dimensionless x-velocitvl
(a)
0.4 . . - . . . . .
s 0-3: i
:3 0.2- ," ~
g _ " . .
'>. 0.1
a “5’
9 0.0 Present Simulation
C Re = 100
0
g —0.1i --_- Re=200
<1) 1 ----- Re = 400
E '0-2 “ Data from Hou et al. (1995) , N
9, i 0 Re = 200 i), [I i,- 1
O ’0.3 l" A Re=400 ‘\\\\.’1 I] ‘i
g 1 Data from Ghia et al. (1982) '._ r 1
3 '0-4[ 0 Re: 100 ‘1’ ‘
A Re = 400
-0.5 ~ ~
0.0 0.2 0.4 0.6 0.8 1.0
X/LO
(b)

Figure 19 Velocity proﬁles along the central lines of the computational domain
obtained by using the physical properties of air without the rescaling when the
dimensionless relaxation time is in the valid range. (a) dimensionless x-velocity (b)
dimensionless y-velocity

90

employed to obtain the lattice sound speed and the acceleration is not rescaled. The
cavity is Im><lm in size and a l28><128 grid is used. We temiinated the computation

when the maximum relative error of the distribution function between two successive

time steps is less than 1x1076.
Firstly, air at 300 K and 0.0001 atm is considered. At this state, kinematic

. . . , 2, . . . . .
VlSCOSIly v IS 0.18245 m /s and the d1mens1onless relaxation time from Eq. (107) IS

I = 0.638, which happens to fall in the valid range of r . Therefore, this problem can be

simulated with this physical r value. Note that if I calculated from Eq. (107) is used, y

becomes 1 and no rescaling takes place for the sound speed. In Figure 19, the velocity
profiles along the central lines of the computational domain are compared with the data
from [135] and [179]. Apparently, the simulation results agree with the results in [135]
and [179] very well. As expected, ifthe dimensionless relaxation time is valid, LBM can

be used with physical properties ofair without the rescaling.

 

 

 

 

 

 

r=0.6 r=0.638 r=0.8 r=1.0

7 0.7252 1.0 2.1757 3.6262
Re 400 400 400 400

Ma 0.1804 0.2490 0.5413 0.9021
Error in 11,. (0/0) 0.27 1.53 3.08 6.58
Error in u... (”/u) 0.99 1.42 2.56 7.28

 

 

 

 

 

 

 

Table 5 Simulation parameters and simulation errors for different values of r in
the simulation of two-dimensional driven cavity flows

91

 

 

 

 

20000
O
\
\.
160004 \
(0
O. \
g \
U) \0
_§ 12000—
'9
a) -1
z o
8000- \\\\\\\\\\
-1 O \\ O
4000 . . . T . .
(15 1.0 1.5 21)

Dimensionless Relaxation Time T

Figure 20 Number of iterations until convergence vs. 7 (Re = 200). A larger 1 leads
to a quicker convergence.

Now, the same problem (Re = 400) is simulated by using the presented rescaling scheme.
In this example, we select four different values of I using Eq. (122) and the

corresponding values of y are presented in Table 5. The errors in 1.1,. and u,. listed in

Table 5 are the maximum relative errors with respect to the results in [179]. Note that as
mentioned earlier no rescaling takes place when r = 0.638. As seen from Table 5, this
scheme gives reasonably good results for small values of 2' but as 2' increases this
method becomes more inaccurate. This increase in error with increasing r can be
explained by Eq. (1 19): at fixed Re, a higher value of r leads to a higher Ma. Because
the standard LBM is formulated based on the low Ma assumption, a higher Ma tends to

make the algorithm less accurate. On the other hand, as seen in Figure 20, the code

92

 

 

 

 

 

 

 

 

 

: N=128
5- . N=256
, - ----- N=512
a 41 I
a: 1.
5 1 in
3 3- \\\'\‘
2— 0\ 750:7: """ ‘-O~‘_____
1— A1 B: c:
0.5 1.0 1.5 2.0 2.5 3.0 3.5

Dimensionless Relaxation Time T
(b)

Figure 21 (a) r as a function of Re and Ma (N=128). It is apparent that 2' needs to
be close to 0.5 if a high Ma and/or a high Re are desired. (b) Re vs. 2' at different
grid densities. Increasing grid density is an effective way to increase Re if the other
simulation parameters are kept unchanged.

93

converges faster as r increases. In Figure 21 (a), the dimensionless relaxation time I is

plotted as a function of Re and Ma (when N = 128) following Eq. (120); Figure 21 (b)

presents the relationship between 7 and Re at three different grid densities when Ma is

chosen to be 0.2. The points marked by A, B and C are the dimensionless relaxation

times calculated from Eq. (120) when Re=100. At such condition, the dimensionless

relaxation times for N=128, 256 and 512 are found to be 0.943, 1.387 and 2.274,

respectively.

 

 

 

 

Til” Pia/”1) Ming/s) T Ma Re Result
Case 1 300 1 I.8245><10‘5 0.500014 1.24><10'4 200 Unstable
Case 2 1500 1x105 20.399 7.394 6.22 200 Unstable

 

 

 

 

 

 

 

 

Table 6 Two examples of the air state that give invalid dimensionless relaxation time

 

 

 

 

 

 

 

 

 

 

Re Ma r 7 # of iterations
Case 1 200 9.021 ><10’2 0.6 7252.43 20380
Case 2 200 9.021><10'2 0.6 1.45><10'2 20380

 

Table 7 Rescaled simulation parameters and # of iterations until convergence

94

 

 

Re Ma r 7

 

100 0.1 0.72 15955.35

 

200 0.0902 0.6 7252.43

 

500 0.113 0.55 3626.22

 

1000 0.113 0.525 1813.11

 

2000 0.225 0.525 1813.11

 

5000 0.564 0.525 1813.11

 

 

 

 

 

 

Table 8 Simulation parameters for air at 300 K and 1 atm (Case 1)

.0
oo
o

 

 

 

 

+>< Re = 10
E075 %Re=100a
E X Present Simulation i
g 0 70 r + Data from Hou et al. (1995) _ ‘
a - x Data from Ghia et al. (1982) +>< Re ‘ 200
m —I
g 0 65 e
f} Re = 400 1
g 0.60 ~ *8 2
% sex Re = 1000
E 0.55 _ % +>< Re = 2000 —
> Re = 5000

0.50 . 1 1 1 .
0.50 0.55 0.60 0.65

X location of primary vortex (m)

Figure 22 Locations of the primary vortices at different Reynolds numbers. The
simulation results obtained by the rescaling scheme agree well with other
researchers’ results.

95

To show the capabilities of the rescaling scheme, the driven cavity flow is now
simulated for two different conditions of air, under which the dimensionless relaxation
times that satisfy Eq. (107) are outside the valid range. Table 6 shows the properties ofair,
and Re is selected to be 200 for both cases. As expected, LBM fails in both cases without
the rescaling. However, both cases are simulated correctly by using the rescaling scheme.
Table6 lists the parameters used in the simulations with the rescaling and the numbers of
iterations until convergence. Note that Ma and 7 are in the reasonable range. Next, we
simulate the driven cavity ﬂow problem at several different Re using the properties used
in Case 1 (in Table 6) with the rescaling. The main simulation parameters are listed in
Table 7. Figure 22 presents the locations of the primary vortices at different Re, and
Figure 23 shows the velocity profiles along the central lines at two different Reynolds
numbers (200, 5000). As seen, the ﬂow patterns are in good agreement with the ones
shown in [135] and [179]. In Figure 23, we can see that the error at Re = 5000 is
relatively large, which can be explained by the high Ma (Ma = 0.564) resulting from Re =
5000 and r = 0.525 as shown in Table 7. According to Eq. (1 19), there are two ways to
reduce the Ma when Re is large. One way is to decrease 2'. But unfortunately, there is a
lower limit of 7 due to the nature of LBM. The closer the value of 2' is to 0.5, the more
likely the code will be unstable. The other way is to increase the grid density for the
simulation. This is the most straightforward way to solve this problem in a moderate
range of Re but of course will increase the computational expense. Besides these two
methods, there is another way to reduce Ma, i.e., the use of an interpolation scheme after

the collision-streaming steps in LBM [180].

96

 

    
 

 

 

 

 

 

1.0 .
1
0.8 ~ -
3 0'6 0 Present Simulation ‘
; Re = 200 4
O 4 - - — - Re = 5000
' ” Data from Hou et al. (1995)
. 0 Re = 200
O 2 _ A Re = 5000
' Data from Ghia et al. (1982)
. ,1 \ 0 Re = 5000
“ ‘ ‘ -.
00 1 1 1‘ ‘ 1 - _ 1 1 A 1 A A a L
050 -0.25 0.00 0.25 0.50 0.75 1.00
Ux/UO (Dimensionless x-velocity)
(a)
0.6 I I I I I I ' I
0.4 ~ "3’. 4
,’ \3
LI! \ “.
0.2 "I /‘—‘~\ q

I F \ ‘ i
0.0 Present Simulation 9\

 

 

Uy/UO (Dimensionless y-velocity)

 

 

 

t Re = 200 ~ ,
-02 a Re=5000 ‘
r Data from Hou et al. (1995)
04 p 0 Re = 200
_ A Re = 5000
Data from Ghia et al. (1982)
'0-6 ‘ . Re = 5000 ‘
0.0 0.2 0.4 0.6 0.8 1.0
X/LO
(b)

Figure 23 Velocity proﬁles along the central lines of the computational domain
(a) dimensionless x-velocity (b) dimensionless y-velocity. The simulation errors are
relatively large if the Reynolds number is high. It is due to the high Mach number

resulted by the high Reynolds number according to Eq. (119).

97

4.3 Simulation Results

To illustrate the capability of the present LBM model with the rescaling scheme,
two problems of weakly ionized helium plasma are simulated here. First, we simulate the
electron diffusion problem under an externally applied electric ﬁeld by neglecting the
internally generated electric ﬁeld because, in this case, the analytical solution is available.
The details on the problem, such as domain size and grid density, are identical with those
in Section 4.1. If the initial electron distribution is Gaussian and the maximum degree of
ionization is 1%, the temporal and spatial distribution of the electron number density is

calculated analytically as:

2 2
nj(x,y,t) :_____0_01n,,0 exp — (x—xc. — vdt) +0) y“) (123)

,
0+dm) V0+dm1

 

where r and ’0 have the same definitions and values as in Eqs. (110) and (113). The
electron diffusivity can be calculated by Eq. (112). In Eq. (123), vd is the magnitude of
the electron drift velocity vd due to the external electric ﬁeld, which is calculated as

follows [154]:

=a x1 (124)

In this problem, the second rule must be applied to rescale the acceleration term. Because
the drift velocity is the characteristic velocity due to the electric ﬁeld in this case,

i.e., U() = vd, from Eq. (1 l8) and Eq. (124) the following equation can be obtained:

~

ae’ien : ﬁe en (125)

from which the acceleration due to the external force is rescaled as follows:

98

 

 

 

 

 

 

 

,, 8772 . (126)
where Eq. (1 17) is used to replace x16”.
”1200774) Te(€V) r 7 50(V /m) error
1><1016 0.7005 1.5 9.686x10'9 14226.72 2.815><10'8
1><1018 0.8008 1.5 9.686x10'7 16270.90 3.465><10'8
1><1020 0.9343 1.5 9.686x10'S 18975.84 2.815><10'8
1><1022 1.1186 1.5 9.686x10‘3 22731.20 2.817><10‘8
1><1024 1.3891 1.5 0.9686 28288.26 2.825x10'8
1><1026 1.820 1.5 96.86 58822.30 5.360><10'8

 

 

 

 

 

 

 

 

Table 9 Simulation parameters used for the electron diffusion problem

Table 9 lists neutral number density (HMO), electron temperature (T0) and other
parameters used for the simulation. Here, Te is obtained from the Saha equation [129]
corresponding to the 1% ionization degree and the electric ﬁeld E0 directs from left to
right. The dimensionless relaxation time 7 is set 1.5 for all the calculations. In this

. . 16 26 -3 .
reseach, a Wide range of neutral number densrty (from 1X10 to 1X10 m ) IS

considered. In Figure 24, electron number density distributions at later times are plotted

together with the analytical solution (Eq. (123)) for two different neutral number density
values (l><10l6 and IMO26 m'3). The simulation errors for six different initial number

density values are listed in Table 9 using the following formula [181]:

99

 

1.0x1014i — - - - initial Distribution ,,
—-—- Net;B
8.0x10‘3— 0 Ne:

“PE 6.0x10‘3a

 

 

 

 

 

 

 

Q) 1
Z 13 I
4.0x10 7 :
l
13 .'
2.0)(10 '1 7, ‘s i \
g" ‘1 1‘ I.
O O m U ‘Iv‘ﬁtr- ”L.” [A], r \\ r A
, main-urn e111:I‘rig-pituuzz:Liixisfi;|~-‘- I- - — - — — — - — - "f‘"JIHIYLHILLIHIU1UUKUTIZUIIHLI-k
0 1 2 3 4 5 6 7
X (mm)
_ -16
(a) 11,,0 =1x1016m 3 att=2.31><10 s
24 . . . . .
1 .0x10 - - - - - Initial Distribution 1.
LB ,' '
—— Nee , i,
23 o A ‘ ‘
8.0x10 7 Nee : 1‘
.' ‘.
00A ‘ l
‘ 23 1 1
E 6.0)(10 7 l i
v 1 I
(D l i
Z l I
23 ; .
4.0x10 ‘ , i
t I '.
I i
2.0x1023~ . ‘.
I \ I 1
4 rt" fit. ,' “
,1“! 11“.} I \
O_O —li'.l.XIIXIY-"IIIIHIIYH-li'l\ltli,¥5L'iu:[— .. _'_ — .1 _ _ ? - —7:11:11};-"llt.{l.f-IIII."Zl|'.iili‘i:.'{\-SI.Vllfiilli'jllylillil
0 1 2 3 4 5 6 7
X (mm)

(b) 11,10 =1><1026m"3 at t = l.433><10-6 s

Figure 24 Electron number density distributions with different initial number
densities of neutrals. The simulation results obtained by the rescaling scheme are in
good agreement with the analytical solutions.

100

 

2

L8 . . A
2 n6 (1,],1)’”e (xli’yjit) mn (127)

 

error = A
i,j ne (1,},t)

where m and ii are numbers of grid points in x and 1 directions, respectively; 11:80,],1)
is the electron number density at node point (1', j) obtained by the simulation; and
iif(.rl-,_)ij,t) is the electron number density at the corresponding space point by Eq. (123).

As clearly shown from Figure 24 and Table 9, the presented LBM model and the
rescaling scheme give good results for the electron diffusion problem for a wide range of
electron number density.

Next, the electron diffusion problem is simulated again by taking into account the
internal electric ﬁeld generated by the charge density. The computational domain as well
as the grid density is the same as before. The initial number densities of the heavy
particles are uniform in the space while the initial distribution of the electron number

density is in the Gaussian form as shown in Eq. (110). The initial neutral number density

of 11,,0 = 1x101811173 and 1% ionization degree are considered in this calculation. Figure

25 presents the effect of internal electric ﬁeld on the diffusion process. In this case, the
external electrical ﬁeld is not applied and only the internally generated electric ﬁeld is
considered during the diffusion process. It can be seen that the diffusion process is
enhanced since the internal electric ﬁeld generates outward forces which make the
electrons escape from the center quickly. The dynamics of electrons under both external
and internal electric ﬁelds is also simulated. Figure 26 shows electron number density
and ion number density at t = 1.1 16 ns. The dashed line in the plot represents the number

density contours of ions while the solid line stands for electrons. Under the effects ofthe

101

 

1.010-

 

, \
/
1.008 — I” \
O 1.006 -1 7 1 \
g 1 I fﬂtx.\
a f I
2 1.0044 1 \\
If \‘\
1 ll \.
1.0021 1 ‘11
[I \l
l ,5.’ ‘1
10004~~~ * 7..

 

 

Figure 25 Electron number density distribution (Circle: without internal electric
field, Squares: with internal electric field, both at 0.223 ns. Solid line shows the
initial distribution). The internal electric field generated by the charge density

largely enhances the electron diffusion process.
external applied electric ﬁeld, the electrons will move along the inverse direction ofthe
electric ﬁeld. Note that the ions respond to the electric ﬁelds very slowly because they are

much heavier than electrons and their temperature is much lower in this simulation. That

is why the contour of ion number density seems to stay at its original position in the

ﬁgure.

102

 

 

 

 

Figure 26 Snapshot of electron number density (solid lines) and ion number density
(dashed lines) at t = 1.116 ns. The electrons and ions move along opposite directions
because they have opposite charges. The ions travel much more slowly than the
electrons since the ions have much heavier mass.

One thing to recall at this point is the fact that the external force in lattice
Boltzmann equation (Eqs. (96) and (97)) is evaluated by using self—collision equilibrium
distribution function (SCEDF) rather than cross—collision distribution function (SCEDF).
To justify this we used both CCEDF and SCEDF for the same simulation and the results
are shown in Figure 27 (without considering of the internal electric ﬁeld). Since it is a

simple convection-diffusion problem, the ratio 719/7160 cannot be smaller than 1.

However, it is seen that when CCDEF is used the number density away from the peak
can take values smaller than the initial value. As shown in Figure 27, the minimum value

is roughly 0.9997, which means that the number density drops by about 3% ofthe initial

103

1.012

 

   
  

 

 

« — Initial Distribution
1010 , - - - Self-Collision
, -'--- Cross-Collision
1.008 a
o .
g 1.006 4
F) .
Z
1.004 4
1 x’
1.002 - ,/
1 /.?'
10004;.23' ' '\._._,_- _________________

 

 

 

Figure 27 Electron number density distribution (at t = 6.69 ns) obtained by using
different equilibrium distribution functions for the external force term. Apparent
simulation errors would be resulted by using the cross-collision equilibrium
distribution function in He’s external force model.
perturbation of electron number density. On the other hand, physically realistic results

were obtained with the used of SCEDF.

The grid independency and computational efﬁciency ofthe model are also studied
by conducting the same simulation on three different grids (64X64, 128><128 and
256><256). The error at time t was calculated by Eq. (127) (no internal electric ﬁeld is
considered). Figure 28 shows that, as expected, second order convergence is observed
from the simulation results. In order to test the computational efﬁciency of the multi-
component model, CPU times per time step are measured on a single-CPU PC for the
present three-component model with the interpolation and a simple LBE model for
electrons only. The electrostatic equation is not considered in the test. Test results show

that the three-component model takes 4.07 times more CPU time than the electron-only

model.

104

 

2.0x10‘3

 

 

 

 

 

 

 

 

 

//O-_———“ O O
1 ,/
-3 0
1.6x10 4
.I
_ (2)1163, o —<>—— 256*256 Grid Density
9 —D—— 128*128 Grid Density
0‘] —o—— 64*64 Grid Density
8.0x10'4-
/D ELMO X0
4.0x10'44
of” O O o
0.0 , , , , . . . .
1 2 3 4 5
Time (ns)

Figure 28 Relative errors in electron number density vs. time for three different
grids. The errors are calculated at the center of the domain with comparison to the
analytical solution. Second order accuracy in space discretization can be seen from

the figure.

As a second validation of the model, we consider the electrostatic wave problem
by neglecting all collision terms in Eqs. (96) to (98). If the collision term in the
Boltzmann equation is neglected, the Boltzmann equation becomes the Vlasov equation
for collisionless plasmas. The initial spatial distribution of electron number density is
perturbed slightly as follows:

27rx
1.1,,(x,t =0)=n(,0[1—0.01cos{T-]] (128)
.1"

where 1,. is the length of the physical domain (3.71 mm). The periodic boundary

condition is employed for both the streaming step for the LBM and the Poisson equation

105

for the electric potential. For this problem, we used both with and without the rescaling. It
is possible because all the collision terms are neglected. Figure 29 shows the evolution of
electron number density and x-component velocity at the center point of the domain. It is
clear that the results obtained with the rescaling scheme agree very well with the ones

obtained using original variables. In addition, the wave period measured from the ﬁgure

agrees well with the theoretical value of electron oscillation period ( 277 / type [129]). The

maximum relative error between the theoretical and simulation results is 0.22 0/o.

Therefore, both the standard LBM and the rescaled LBM can be used for the collisionless

Vlasov equation.

 

 

 

 

0 Number densit El Velocit

16.005 - y y 4 4.0x104

. 7". 'v IT‘T‘ 1:41:23 17" ’ ”(F’ﬂz-sk-i X; 77".:

'1 -.- .43 2 g 9.1 715 1
1 . 1' 8: ’1' 2‘ ‘13 4

\\ ’7 1" q 1‘ (I r" ‘ 2.0X10

A 16.0024 \ I): 1» 1‘- , y 1’ it. 4 C
a) \ [I I t \ 1’ i 1 (>13
a _-__-.\ ...... l ______ I ..... L! ..... \ ..---___1: ..... l ----- .‘. --00 A
c» 1 \ I I, \ I I 3
o l \ I \ \
‘J ‘1 \ I I \ \ I I 1 ~ 3

15.999 - ’1, \\ 1! [I I \\ (J I] f
‘t t I" ,- ‘1, t (I , ~ -2 0x104

1‘ 'I" 2 \r
'l t] :12 [i til“ 5 ,1
”(LIL-3* ‘ 1 ‘11m \ 1
15.996 — 2,3,?“ « -4.0x104
' I ‘ I ' I ' I
0.0 0.5 1.0 1.5 2.0
Time (ns)

Figure 29 Time evolution of number density and x-component velocity of electrons
at the center of the computational domain. (Solid lines: simulation results obtained
without the rescaling scheme, circles and squares: simulation results obtained with

the rescaling scheme)

106

In summary, a lattice Boltzmann method for weakly-ionized isothemial plasmas
has been presented. A rescaling scheme has been proposed by which the LBM with an
external force term can be used with physical properties of fluids, so that effects of
collisions are taken into account correctly in LBM. This scheme only rescales the sound
speed and the acceleration term due to external forces based on the following two rules:
(1) the physical viscosity is equal to the lattice viscosity, and (2) the characteristic
velocity due to the extemal force is not affected by the rescaling scheme. The 2D driven
cavity flow is simulated to validate the rescaling scheme. Finally, the newly developed
LBM model for weakly ionized plasmas with the rescaling scheme has been applied to
simulate the electron diffusion problem and the electrostatic wave problem. Simulation
results agree well with the data available in literature and the analytical solutions.
Especially, this scheme can be applied to a wide range of number densities in plasma

simulations.

107

Chapter 5
Lattice Boltzmann Simulation of Weakly Ionized Helium Plasma with Ionization

and Recombination

A new LBM model is developed in this chapter to simulate laser interaction with
weakly ionized plasmas. Comparing to the model proposed in the previous chapter, the
present mode has the following features. First, ionization of neutral particles in the
plasma and its inverse process, three-body recombination, have been taken into account
in this model, which is achieved by adding an additional term in the Boltzmann equation;
Secondly, the evolution of temperatures of the three species of particles are resolved by
solving a set of energy equations derived directly from the Boltzmann transport equation;
and thirdly, the Lorentz force is evaluated in this model as the external force term and the
ﬁnite-difference time-domain (FDTD) method is adopted to solve the full set of
Maxwell’s equation to obtain the time- and space-dependent electromagnetic ﬁeld. The
interaction between a continuous laser beam and a weakly ionized helium plasma is
simulated by using this newly developed model in this chapter. The preliminary results
that show ionization and recombination dynamics of the weakly ionized plasma are

presented following the detailed description ofthe model.

5.] Mathematical Model
As the governing equations of the present model, the Boltzmann equations are

listed based upon the following assumptions:

108

l) The plasma is weakly ionized. Therefore, the long range Coulomb interactions
are not considered.

2) Only electrons at the outmost shell are ionized. Consequently, there are only
three species of particles in the plasma: electrons, ions and neutrals.

3) Only electron impact ionization and three-body recombination are included in
the present model. Other inelastic collisions, such as excitation, de-excitation, and charge
attachment are neglected.

Based upon the above assumptions, the Boltzmann transport equations for

electron, ions and neutrals in a weakly ionized plasma can be written as:

 

‘\ t?(]
(.1 ) ) ) ac). v)—u) )
-‘—[‘—+vL,-Vf,=— f ﬂ————’—’l’ +Reffq+ ( ‘ ‘) cf" (129)
C}! (311 I!
C)
. (?(/
(‘7. ‘ 1...". )
4+V,V/,— —f‘ /——~—'——" +Rf“’+———( ’ ’)f,.“’ (130)
("I Am 0,2
i? J/c __ ‘(3(] ,
.— +.-W,,=——"——— (131)

7171

where the external force term a -)-va.(,-) has been evaluated by He’s model as

14(1
described in Chapter 3.

There are two main differences between the above Boltzmann equations and Eqs.
(90) to (92). First, one term is added on the right hand side of the above equations to
represent the rate of change of the particle distribution due to ionization of neutral
particles and recombination of charged particles. In this model, only electron impact
ionization and three-body recombination are considered. Then the physical ionization and

recombination (1R) coefﬁcient for electrons can be expressed as:

109

&=@—@ on)

where R: and R: represent the electron impact ionization rate and the recombination
rate, respectively. The electron impact ionization rate can be expressed as:

R; = 6,11,, (17,.) (133) ,
where (7,.) is the average value of the thermal speed of electrons and
<E.>=./8ltBT(,/miic according to kinetic theory [126] and 0,- is the cross section of
electron impact ionization which is dependent on the electron energy nicvg/Z. The data

of 0,- for helium as a function of electron energy can be found in many literatures [182-

184]. For the three-body recombination, the formula in [185] is adopted in this simulation:

 

. _ 0.822x10‘33

R], T45 nan,- (134)
e

L

where the electron temperature Te is in the unit of eV. Once RC, is calculated by Eqs.

(I33) and (134), the IR coefﬁcients of ions and neutrals can be found by R,- = 112196 and
I.

’I

,7 1 0
— _( ') r:
R — Re, respectively.
)1

The second difference is the deﬁnition of the cross-collision equilibrium

distribution function. f5? is deﬁned in this model as:

 

2
f5? = ""7 exp -(—v—u5—)— (135)
2770;” 26,,"

110

Comparing to Eq. (93), a new deﬁned sound speed 68,, replaces the original electron

sound speed:

R T.
8 L12 (I36)

’71

en 2
e

where T8,, is the electron temperature after elastic collisions with neutrals in the time
period of 21m. From the kinetic theory [126] and considering the ionization and

recombination process, T can be found as:

(’71

 

 

2
+ 2mcii(Tn " Tc) + mnmen(un '— ue) _ Z’iwiRcUie

(137)
int, +171" 3k3(me + m”) 3kg

U)! L’

HIJII. . .
where mm =-—‘—1— IS the reduced mass of electron and neutral, U,- 15 the ﬁrst
NI '1'”!
L‘ H

ionization potential of helium (U,- =24.59eV ) and e is the unit electronic charge

(Ie =1.6X10_19(’). The T6,, is adopted here because, to correctly describe the ionization

and recombination dynamics, it is necessary to describe the post-collision states of the
particles in the plasma, including both the self-collisions and cross-collisions. It can be

seen later that the adoption of T in the equilibrium distribution function just describes

(’11
the heat transfer between electrons and neutrals due to the collisions, which is signiﬁcant

for the energy evolutions between different species of particles in a system where the

composing particles have far different properties. Similarly, the cross-collision

equilibrium distribution function fljlq has the following deﬁnition:

-2
fig] 2 l 2 €Xp _( I 2111) (138)
2776- 26,-”

IN

111

/' T- .
where 0m = ‘8 ’” and T- IS calculated as:

[I]
”11’

 

 

2
2171-(T —T-) m m: (u -u)
_ III I? I H In I? I
Tm — T,- + V , + 3k (139)
111,- +111" 8("11‘ +171”)
[HI/H”

where the reduced mass of ions and neutrals is ”1m: . Note that there is no

INI- +171"

ionization term in the above equation. That is natural because there is no ionization
happening during the collisions between ions and neutrals in the plasma.

The Boltzmann equations shown in Eqs. (129) to (131) can be discretized by the
traditional D2Q9 scheme. However, in such a case, only particle number densities and
momentums can be retrieved correctly from the lattice Boltzmann equation. Thus, a
supplemental set of energy equations are included in this model to describe the heating
patterns ofthe particles. Starting from the Boltzmann equation, the corresponding energy
equation can be obtained by taking proper moments of the distribution functions with
respect to the microscopic velocity. Taking electrons as the example, the distribution

function f, has the following property:

71(,<Z(v€)> = (1(ve)f,,dv, (140)

where ﬁve) is a function only dependent on ve. By setting different expressions of
1(vc), we can have the following macroscopic quantities listed in Table 10. In Table 10,

CL, = Vt, — u(, is the random velocity ofelectron which is closely related to the temperature.
lV'llt'l"'1-'-l ’2 ‘E 179‘d't0‘t' lt "tl
u ip ying [hey—3111,.” on every term in q. ( .. ) an in egra mg eacr erm 111 1

respect to v", from the negative inﬁnity to positive inﬁnity, we can get:

5' 1 1
II<vy)—g;—dve + IZ<V€)VL’ - ijedve : __
V v

 

[1(Ve)(fe - f5." )dve

 

 

 

 

 

 

 

 

[it’ll
, v ,
(141)
er]! at) ' “(10'
+ Re viz/(”)1”), (. ve +37 J-Z(v(,)(ve —u£,)fe ve
v e v
Macrosco ic .
Quantity}? [(vp) Expressron
. n, = ,(lv ,
Number Densuy 1 ‘ ift ‘-
V.
Macroscopic v ) I1e<Ve> = Jivefedvc : neue
Velocrty ‘ v
Kinetic Pressure c,c mene (Cece) = me Icecefedve = P8
Dyad ‘ e v
Momenéum Veve ”Idle <veve> = me JVeVefedVe = PL, + merieueue
Dya v
1 l 2 _ l 2 I _ 3
Temperature — mace -ce ”e <_mece> — "me Icej‘edv e “ _ ”ekBTe
2 2 2 v 2
l l ,2 _ l . 2 _ 3 , l 2
Total Energy Emeve -ve ”e <‘2'me‘e> — Eme Jive/ed“) — E’QJI‘BTe + S’nt’nc’ut’
v

 

 

 

Table 10 Macroscopic Quantities from Moments of Distribution Function

113

 

Using the expressions listed in Table 10 and following the similar derivation

procedure in [129], we can finally ﬁnd the energy equation ofelectrons as:

q

 

 

('8,
‘3‘ + v'(€eue) : ‘V'U’e 'll(,)—V°qe + mencae 'ue
CI
(142)
3n ,k T, —T, m,n,
_ L 32:; (I?) _ 2:1 L (“g—ug)1)+Re£e
en en

3 l . .
where 6, =—n,k T, +—m)n,u% is the total ener of electrons; and , 18 the electron
L 7 L B L 7 L L c L

‘—

l _ .
heat ﬂux and qc,=—KL,VTL, where A’e=§Iz(,kB<v(,>/lm 18 the electron thermal

_e(E+uL, xB)

m

where E and B are the electric ﬁeld

 

conductivity. Recalling that at, =

L,
and magnetic field respectively, the third term on the right hand side of Eq. (142) can be
rewritten as:

mcneae -ue =—en(,(E+ue xB)-ue (143)
Noting that (u(, x B) - up = O and defining the electron current density as J6 = —en£,ue , the
above equation becomes:

meneae cue = Je -E (144)

which can be interpreted as the Joule heating term for electrons. Further, substituting Ten
(Eq. (137)) into Eq. (14.2) and after some manipulations, the energy equation ofelectrons
can be expressed as:

“:1" + V-(ccue) = —V-(P€ -ut,)+ v.(;c(,vr(,)+J(, - E

 

(145)

3kg]? m
L)
—' UL,” ' ("UN - He) — Renal/,6 + R688

(’H

 

<7; — m +

2.8,, (me + m”) [ten

Similarly, the energy equations of ions and neutrals are:

114

ﬂ

5?+V-(a,u,) = —v-(P,- -u,-)+V-(K,.VT,)+J,--E
Cl

 

(146)
3k -- - 7‘. —T -
— BHIWM( I H) + ”11”” "m '(“m Tui)+ Rig!
xt,l,,(m,~ +112”) A,”
and
68
#+V.(Elluzi)=—V'(Prz 'un)+V'(KnVTe)
3k mm, T,—T m n,
+ B . wt. "H N me, .(u,,,, —u,,) (147)

 

Aer: (”18 + ”In ) 2'6”

+ 3k8’11'nl' (7} — T”) + mun,

H1

 

2 2 “in '("m -ll”)+Rn£”
rm (m,- +122”) em

Eqs. (145) to (147) constitute a full set of the energy equations of the three species of
particles in a weakly ionized plasma. The temperatures can be solved from the above
equations by using the hydrodynamics quantities calculated from the Boltzmann
equations. However, evaluation of the Joule heating term in the energy equation as well
as the external force term in the Boltzmann equation requires an accurate solution of the
Maxwell‘s equations. In this model, the ﬁnite-difference time-domain (FDTD) is adopted

to directly solve the full set ofMaxwell’s equations.

5.2 F inite-Difference Time-Domain Method

The extemal force term in the Boltzmann equation is generally the Lorentz force
in plasma simulation. This extemal force essentially govems the dynamics of the
individual species ofcharged particles and thereby inﬂuences the macroscopic quantities
of the plasma. Inversely, the variation of charge density in the plasma affects the
electromagnetic field which determines the Lorentz force together with the macroscopic
velocities. Meanwhile, the Joule heating term in the energy equation can also be

interpreted as the work done by the Lorentz force on the charged particle. Conclusively,

115

the most detailed information about the evolution of the electromagnetic ﬁeld is needed
so that the Lorentz force in LBM solver can be evaluated accurately. As the most exact
governing equations for all electromagnetic phenomena, Maxwell’s equations are solved
directly in this research to obtain the ﬁeld quantities. Maxwell’s equations consist of the

following four equations:

VWng. (M&
VRB=O 04%
Vsz—QE—ht 05m
a
VxH=J+§B USU
a

where ,0, is the free charge density, pv =e(n,- -ne); J is the electric current density,
J =e(n,-u,- —n(,u(,); M is the equivalent magnetic current density which is neglected in

this study; E is the electric ﬁeld; H is the magnetic ﬁeld; D is the electric ﬂux density;
and B is the magnetic flux density.

To solve the Maxwell’s equations, the ﬁnite-difference time-domain (FDTD)
method [186] is adopted in this research. Figure 30 shows the positions ofthe electric and
magnetic ﬁeld vectors about a cubic unit cell ofthe Yee space lattice. This conﬁguration
ensures Eq. ("149) to hold. Accordingly, the Maxwell’s equations can be discretized by
using the central-difference scheme upon both space and the time domain. The concrete
expressions of descretized Maxwell’s equations for different cases can be found in
Taflove and Hagness’s book [186]. Then the new value of an electromagnetic ﬁeld vector
component at any lattice point can be calculated from its previous value and the previous

values ofthe components ofthe other ﬁeld vector at adjacent points.

116

Ez(i-1/2, j+1/2, k+1)

 

’ Hy(i-1,j+1/2. k+1)
/ | / : /.l .
g ,21 ; ,Hx(I-1/2, 1+1, k+1)

 

 

Hx(i-1/2. j. k+1) f Hy“, j+1/2, k+1)

Ey(i-1/2. j, k+‘l/2) HZ(i-_1,j. k+1/2) Ex(i-1. j+1/2. k+1/2)

.. —--:—--l——:f—-—+- _ _
Ex(i,j+1/2,k+1/2) // EytI-1/2. 1+11k+1/2)
d ‘ i Hz(i,j+1,k+‘l/2)

Z __.;._J.’_- -_-
l ' / Hy(i-‘l,j+‘1/2, k)

i
Hx(i-1/2, j, k) ------- #7:"— ----- ----
Lp/
l

 

*l

/

qu, j. k+1/2)

 

—-qr——-——

 

 

 

Hx(i-1/2. j+1, k)

 

1. j, k Hy(i, j+1/2, k)
y Ez(i-1/2, j+1/2, k)

 

X.

Figure 30 Position of the electric and magnetic field vector components about a

cubic unit cell of the Yee space lattice

Since the electromagnetic waves in nature can reflect back into the computational

domain at the computational boundaries, a proper boundary condition must be applied to
terminate the computational domain if a physical problem with Open boundaries is of
interest, such as an electromagnetic wave propagation in an inﬁnite medium. In order to
eliminate the fake waves caused by the reflection ofthe incoming wave at the boundaries,
Berenger introduced a highly effective absorbing boundary condition called the perfectly
matched layer (PM L) in 1994 [187]. However. this PML boundary condition is a
hypothetical medium based on a mathematical model. Indeed, a physical model based on
an anisotropic perfectly matched medium can be formulated. This was ﬁrst discussed by

Sacks er al. [188]. For a single interface, the anisotropic medium is uniaxial and is

117

composed of both electric and magnetic permittivity tensors. This kind of boundary
condition is called UPML (uniaxial PML).

To illustrate the capabilities of FDTD for simulation of laser-plasma interaction, a
femtosecond laser pulse interaction with a target silicon substrate is simulated here. The
laser-induced ionization is included in this calculation by using an ionization model based
upon local energy balance analysis of electrons at successive time steps. The current
density is modeled by Ohm’s law, i.e. J = O'E where 0 is the conductivity of silicon at a
given point, in using the FDTD. A two-temperature model developed by Kaganov et al.
[189] is used to calculate the heating pattern of the electrons and lattices in this

simulation.

 

15 um

 

 

50 um

ital

 

 

 

 

OI

 

Figure 3] Schematic of computational domain for simulation ofa femtosecond laser
pulse interaction with a silicon substrate.

In this simulation, a femtosecond laser with the wavelength of 800 nm and the

pulse width of 80 fs is used. The beam is assumed radially polarized (TEMOO mode) and

118

is focused through a numerical lens located at the top of the computational domain. Four

pulse energies (7.5 [U], 15 yJ, 3O ,uJ and 60 ,uJ) are used for the calculation. The

target silicon with the electrical conductivity of 1.4 9‘1 m is assumed. Figure 31

illustrates the computational domain used in this study.

Figure 32 and Figure 33 are the time histories of the radial component of the

electric ﬁeld (5,.) and the corresponding energy absorption patterns in the silicon
substrate respectively. The pulse energy is 30 ;JJ ; and the pulse width is 80 fs. The top

row in each ﬁgure shows the results obtained without the ionization model; and the
bottom rows are the results after incorporating the developed ionization model. As shown
clearly, the major difference is that the laser beam is blocked signiﬁcantly near the
interface when the ionization process is accounted for while the laser beam propagates
through the silicon freely when the ionization is not considered. In fact, the laser pulse is
blocked almost completely except near the rim of the laser pulse. This is because the
ionization process leads to a sudden increase in the electron number density. As a result,
the electrical conductivity increases accordingly and the silicon substrate can behave
almost like a metal. This can be explained more clearly from the energy balance
viewpoint. As illustrated in the bottom row of Figure 33, almost all the energy carried by
the laser pulse is dissipated in a very thin layer near the interface due to the very high
electron number density. The thickness of this energy absorption layer is called the skin
depth (for metals) and it is believed to be closely related to the ablation thickness. As

clearly seen from these ﬁgures, the thickness is very well deﬁned even for a

119

 

 

 

 

(B) With the ionization model

(1)t=-11.6 fs (2)t=65.2 fs (3)t=l40.2 fs (4)t=218.8 fs

Figure 32 Time history of laser pulse propagation. (In these plots, the radial
component of the electric ﬁeld (Er) is used.)

120

 

 

 

 

 

(B) With the ionization model

(I) t = 65.2 fs (2) t = 103.6 fs (3)1: 180.4 fs (4) t = 257.2 fs

Figure 33 Time history of laser energy absorbed by electrons

121

semiconductor like silicon and this explains the superior ablation quality of femtosecond
laser pulses .for various materials. In the mean time, it is also shown that near the outer
rim of the laser pulse, the EM wave penetrates into the substrate. This “leaking”
phenomenon occurs because, for a Gaussian laser beam, the laser intensity at the outer
edge is not high enough to induce strong ionization.

Other simulations results, such as the evolutions of the electron temperature,
electron number density, electrical conductivity and the absorption coefﬁcient of the
target silicon are also obtained through the simulation and they are discussed in detailed
in [89]. In conclusion, by using the FDTD method, the wave nature ofthe laser beam can
be retrieved to the greatest content. Thus, the FDTD method is adopted in the present
model to calculate the electromagnetic ﬁeld which is essential in evaluation of the
external force term in the Boltzmann equation and the Joule heating term in the energy

equation.

5.3 Numerical Implementation

In the practical implementation of the present LBM model to simulate laser
interactions with weakly ionized plasmas, there are several points worthy to be noted
numerically.

1) Selection of the time step. Because the present model is actually a hybrid

model in which three sub-models involve, the time step should be selected as:
A! = min(AtLBM ’AIFDTD’AIENG) (152)
where A’LBM is the time step determined merely by lattice Boltzmann method,

Ar
A1 LBM = 73:61——

‘ 11111.1(

(where 6m,“ denotes the maximum sound speed ofthe three different

species of particles in the weakly ionized plasma and 6max is usually equal Be in the

. . . A
most cases); A’FDTD IS time step merely dec1ded by FDTD method, A’FDTD S—l
c

(where c is the speed of light) and A’ENG is the time step required by the solution of the

,
. 0.5C Ar“ . .
energy equation, A’ENG S ———e———, where C6 and KL, are the volumetric heat capac1ty
K

e
and thermal conductivity of electrons respectively. Through preliminary quantitative

probe, we can ﬁnd that A’FDTD is the usually the smallest one among the above three

mentioned time steps. Thus, we choose the time step as:

Ax
AletFDTD =5 (153)

where the constant 6 2 2 in our simulation without loss of generality.

2) Rescaling scheme. Although the basic idea ofthe rescaling scheme used in this
model is the same as the one described in Chapter 4, there is a little difference in the
practical implementation. Besides the two free parameters, the spacing Ax and the
dimensionless relaxation time I , there is another free parameter which can be prescribed
here; that is, the time step At as described above. The lattice viscosity recovered from the

LBM has the same form as before (taking the electrons as the example):
0, =(re—o.5)é§m (154)
where the time step is determined by Eq. (153). By equating the physical viscosity to the

lattice viscosity, we can easily have:

37rr—0.5At
ye=\[ (e ) (155)

 

 

8/1

6}?

123

following the same deﬁnition ofthe rescaling parameter y as shown in Eq. (I 15). It can
be seen that once the physical properties of the electron, the dimensionless relaxation

time rt, and the time step At are determined, the lattice sound speed of electrons can be

determined accordingly as:

‘ (156)

 

The rescaling of the relaxation time xi is straightforward, i.e. 2,," = reAt. Then the

U!“ I

rescaled acceleration can be found as:

V" a, (157)

 

according to the second rule of the rescaling scheme introduced in Chapter 4. The
rescaled parameters for ions and neutrals can be found in a similar way.

3) Interpolation scheme. The distance traveled by electrons in a time step At is

At',.=\/30~L.At . Comparing to the spacing Ax , it is apparent that the ratio

Are/Ar:x/391.At/Ar=\/3—6~c/§c is less than one because the lattice sound speed is

impossible to be greater than the speed of light. In another word, during a time step At,
the electrons cannot stream to the adjacent node points, so do the ions and neutrals. Thus,
it is necessary to apply the interpolation scheme to ﬁnd the on-node value of the
distribution functions after the streaming step in LBM. The interpolation scheme applied

here is similar to the one described in Chpter 4. However, since the lattice sound speed

~

0,, is dependent on the temperature, and consequently, dependent on both time and space,

the interpolation fomiula used here is different with the one shown in Eq. (101). Instead,

by referring Figure 34, the formula takes the following form:

124

 

f<P>= ”WP” _ "mate? , kmerl')
o (II~I7I)(k—m) (k._.m)(n__k) (nmenh-k)

(158)
where m = 5.4/7), 11 = ﬂc.(q)+1 and k = ﬂ£.(0)—I . Note that the interpolation parameter
now is dependent on both time and space and it is deﬁned as [3,.(x,r) = J3é.(x,t)At/Ar.

It is easy to see that if [3,, is uniform throughout the whole domain, Eq. (158) becomes

Eq. (101). Further, if A. =1, few) 2 fe(0') which means no interpolation exists.

 

 

 

Ax
1‘ +1
I I H7—
0 0 p p q 6!
k—A l<—>l l+>l
New ﬂew/5x ﬂe(q)Ax

Figure 34 Schematic of the rescaling scheme. The value of distribution function atp
can be found by using the post-collision values at 0’, p’ and q’.

4) Positions of simulation variables in the grid system. Unlike the isothemial
plasma model developed in Chapter 4, the present model also has the temperature as a
variable. The primary variables used in the calculation are listed in Table 11. All
thevariables listed in Table 11 are dependent on both space and time. But in the actual
simulation, not all the variables are located on the same position in the grid system. That
is due to the staggered scheme of the ﬁeld vectors used in FDTD method. In order to
satisfy the divergence free law V~B=O , the electric ﬁeld and magnetic ﬁeld are
discretized on staggered node points as shown in Figure 30. The grid system as well as

the positions ofthe primary variables is shown in Figure 35. As illustrated, magnetic ﬁeld,

number density, macroscopic velocity and temperature are located at the primary nodes
which are denoted by the white knots while the x- and y—component of electric ﬁeld are

located at the staggered grid nodes which are denoted by the solid knots respectively.

 

 

 

 

 

 

 

 

Number density n5.
Macroscopic velocity u 3
Temperature of electron TS
Electric ﬁeld E
Magnetic ﬁeld H

 

Table 11 Primary simulation variables

 

 

Hz. nsi Us. Ts f

l
I l
l - . I
I, +1
ah + o” )+ o
I l
_.__. ____:__-_T_EX..L___.._ _____
E E
gym/ii y <(tr) y 4941.1)
1
l
|

——-—o-———l———-6—§-’§—l——--$— —————

j\
\I
/\

 

 

 

*1

Figure 35 Positions of the simulation variables. Hz, ns, ns and T5 are calculated at

the white node points while Ex and Ey are staggered in space and evaluated at the
black node points.

126

Besides the above mentioned primary variables, there are also many secondary
variables which depend on the primary variables used in the simulation. They are listed in

Table 12.

 

 

 

 

 

 

 

 

 

Speed of sound 63
Current density J

Relaxation times 25”
Field ionization rates Rs-f
Impact ionization rates R:
Recombination rates R;
Volumetric heat capacity CS
Thermal conductivity KS
Interpolation factors ﬂs

 

 

 

 

Table 12 Secondary simulation variables

Note that in the evaluation of the variables which involve electric ﬁeld, such as
current density and Lorentz force, the value of electric ﬁeld should be calculated at the

primary node point (1',j) . In the practical simulation, wejust use the average value ofthe
electric ﬁeld neighboring the node point (i,j) :

E, = 0.5(Ex(i,j) + E_r(z',j + 1))

E, = O.5(E_},(i,_/)+ E.,.(,'+ 1.1)) (159)

127

Besides the space discretization, the electric ﬁeld is also staggered in time. Thus, the
value ofelectric ﬁeld should be used at the same time level in the time marching process:

Exit—l : 0...)!5(E\H—32 + Exn—liilz)
Evﬁ‘i : 0.5(E).”_3/2 + Eyn—l/Z)

(160)

In conclusion, the Boltzmann equations, Maxwell’s equations and energy
equations are coupled to establish the present model. The Boltzmann equations are solved
by using the traditional D2Q9 lattice Boltzmann method. In order to match the physical
properties of the plasma, the rescaling scheme is applied to select the simulation
parameters used in LBM. Besides the standard collision step and streaming step, an
interpolation step is needed in this LBM model to ﬁnd the on-node values of the
distribution functions ofthe plasma particles. As the governing equation ofthe evolution
of the temperatures, a set of particle energy equations are derived directly from the
Boltzmann equations by taking the proper moments of the distribution functions. The
energy equations are solved by the ﬁnite volume method. The FDTD method is adopted
here to solve the Maxwell’s equations where the UPML boundary condition is applied to
terminate the outgoing electromagnetic waves. Figure 36 shows the ﬂowchart of the
present model. It can be seen that there are complicated correlations among those three
equations. The external force in the Boltzmann equation is calculated by using the
solutions of Maxwell’s equation. The sound speeds, relaxation times, ionization and
recombination rates are evaluated by using the results ofthe energy equation. When using
FDTD method to solve the Maxwell’s equations, the current density needs to be attained

from the Boltzmann equation. In solving the energy equation, the Joule heating tenn is

calculated by using the results from both Boltzmann equation and Maxwell’s equations

I28

while the other transport coefﬁcients are evaluated from the results of Boltzmann
equation. Figure 37 shows the correlations among the three equations. The elements

through which one equation affects another are illustrated in the ﬁgure.

Start I

ll
Initialization

 

 

 

 

 

 

 

¥

 

 

 

’1
Store lain-3V2
ll """"""" 'E By using 11”“1

 

 

 

 

 

 

------------------

 

Compute and Store En—l/2

 

 

l

 

 

 

 

Find E"—l
l

Compute macroscopic quantities:
J, FS =qS(E+uS x8), (93, 23”, RS
at time level n—1

Rescaling Scheme: (93, 23”, F3, ,BS

 

 

 

 

 

 

 

 

 

 

 

 

 

 

n+1

 

 

f

 

Equilibrium distribution functions: f5?”

 

 

 

 

ll

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Collision Step
I!
Streaming Step \
Compute H" 1 Compute 7}"
Interpolation Step

 

 

ll

 

n n
Compute n, and us

 

 

 

 

 

 

 

 

 

 

Figure 36 Flowchart of the present model

130

 

Boltzmann Equations

 

 

 

   
   
     
 
  
 

Ionization
Ionization

Recombination

    
 

Recombination
External Force

Current Heat Relaxation Times
Density Transport
Coefficients

Sound Speeds

Joule Heating

 

 

Joule Heating
Maxwell’s Equations > Energy Equations

 

 

 

 

 

 

 

Figure 37 Correlation among three equations

5.4 Simulation Results
The interaction between a continuous C02 laser beam and a weakly ionized

helium plasma is simulated by using the model developed here. The computational
domain is a I35.68><135.68 umz square which is discretized by a unifomi 256><256 grid

system. As illustrated in Figure 38, the computational domain is surrounded by a PML
layer used as the boundary condition in the FDTD solver. The PML thickness is 30 grids.
The origin point of the computational domain is located at the left lower comer. The
interface between helium and vacuum is atj = 200. A continuous C0; laser beam whose
wavelength is 10.6 am is incident from the top of the computational domain and
propagates along the negative y—direction. The no-gradient boundary condition is applied

on all the boundaries enclosing the helium. Initially, the helium is assumed to have a

131

0. 019/ uionization degree. The initial number densities of neutrals, ions and electrons are

2.087xio‘9m‘3 , 2.687x10‘5cm‘3 and 2.687x10‘5cm‘3 , respectively. The

temperatures of all the three species of particles are set identical. According to Saha’s

equation, T 0 — -7}0- — T110 =1 .05.8€V

 

 

 

 

 

 

 

 

 

 

Figure 38 Schematic of the computational domain for simulation of interaction
between a continuous C0; laser beam and a weakly ionized helium plasma.

After the laser beam crosses the interface, the free electros in the weakly ionized
helium plasma will be accelerated and heated by the laser ﬁeld. As a result, the electron
energy, including both the macroscopic kinetic energy and the thermal energy
(temperature) will be increased. If the electron energy is higher than the ionization
potential of helium, new free electrons will be generated through the electron impact

ionization. Figure 39 shows the evolution ofthe maximum electron number densities

18

 

 

 

 

. — 1 x 107W/cm2 - - - 2 x 107W/cm:2
16- ----- 4 x107W/em2 e x 107W/cm2
< ------- 8 x 107W/cm2
14 q - I" -
12 — ‘"
\O T f." I", .JJT ,_ .. r- '
E/ 10 -« I: I: ,2! ' 2- A ' r
x '. 1' t ,
w 1 ' » ,- '
E _ I.’ .‘I [V /
a) 8 :. '71 I" ’1
C ‘ .: .I’ II, II,
6 “ ’i I, l”
‘ J I [I] ’1,
4 .4 I ,’ / ’
. I ,’
2 q I I, 4;:/
0 7‘ ‘5’ ‘7'. T I ' I a 1 ' r ' r Y I
0 50 100 150 200 250 300

Time (ps)

Figure 39 Evolution of maximum electron number density with different laser
intensities. At the early stage, the electron number density increases exponentially
because of the laser-induced ionization. However, it shows saturation later due to

the equilibrium between ionization and recombination.
under exposure to different laser intensities. Note that the values shown in the ﬁgure is
actually the maximum ratio of the electron number density to the initial neutral number
density. Thus, it also can be interpreted as the maximum ionization degree throughout the
whole computational domain. From the ﬁgure, the exponential increase of the electron
number density due to the electron impact ionization can be seen clearly. We can also see
that the electrons are generated with a higher rate ifthe laser intensity is higher. That is
natural because higher laser intensity indicates that more laser energy is input into the
domain in a unit time and thus leads to a quicker heating of the electrons. As long as the

electron energy is above the ionization potential, the ionization will keep taking place and

more and more new free electrons will be generated. However, ionization itselfis a

133

 

— 1 x107w/cm2
— — - 2 x 107w/em2
4 x 107W/‘cm2
e x 107W/cm2
8 x 107 W/cm2

o‘—-—..—. ""’_'— _—'-i-.-.-i
—td"\.-._oo.-...... ..— — - '-

Temax (eV)

 

 

 

 

O 1 r 1 ' 1 ' r ' 1
0 50 100 150 200

250 300
Time (ps)

Figure 40 Evolution of maximum electron temperature with different laser
intensities. The electron temperature increases very quickly due to the rapid heating.
It decreases after the ionization takes place because the generation of new free
electrons consumes the laser energy.
cooling scheme for electrons because, after all, the electron energy obtained form the
laser beam is used to overcome the ionization potential and thus less laser energy is used
to heat the electrons. This is evident in the energy equation where the ionization is just
the sink term. Consequently, with more electrons generated, the electron temperature
begins to decrease as shown in Figure 40. In this plot, the electron temperature increases

very quickly at the early stage when the ionization hasn’t taken place effectively. This

can be veriﬁed by referring to Figure 39. For example, under the laser intensity of

8x 107 W/cm2 , the maximum electron number density starts to increase around 18 ps,

which corresponds to the temperature decrease at about the same time moment as shown

134

 

—¢-— 1 x 107W/cm2
—-°— 2 x 107W/cm2
—o—- 4 X 107 W/cm2
—e-— 6 x 107W/em2
—9—— 8 x 107 W/cm2

Ionization Rate (8'1)

1:)

.4 J 1?.
$1.3“ '» s v. 1

'7" 93t\.§ﬂ'0¢ lu¢w(,?lll‘}

 

r
35,-.“
l; r—_.‘,‘
-T-'m1,eﬁn;.‘.‘.a.~_.‘ \
Inf?“ "_—'.4:.y1"r“:r_
in] 1 “V11 1‘", . .. .
'DL'vtv _‘,,IC.. 1 [‘-’.3 ‘ '

”MILD:

 

 

 

 

 

T

—

r

" ' 1 1 r 1 1 1 ' 1 1
O 50 100 150 200 250 300

Time (ps)

Figure 41 Evolution of ionization coefﬁcient with different laser intensities. The
ionization coefficient shows a similar pattern with the electron temperature because
it mainly dependent on the electron energy.
in Figure 40. As the electron temperature decreases, the ionization rate is reduced
accordingly. At the same time, since the recombination rate is proportional to the number
densities of charge particle and inversely proportional to the electron temperature, the
recombination becomes more and more prominent in the plasma. Reduction of ionization
and enhancement of recombination lead that the generation rate of new free electrons gets
lower and thus the electron number density shows a “saturation” phenomenon in Figure
39. This can be explained more evidently by Figure 41 and Figure 42 which are the
evolutions of the ionization coefﬁcient and the recombination coefﬁcient respectively.
Referring to Eq. (133), the ionization coefﬁcient is mainly dependent on the electron

energy. That is why the evolution ofthe ionization coefﬁcient shown in Figure 41 has the

 

3,0x1o‘°~ — 1 x107W/cm2 — - - 2 x107wxem2
~»--- 4 x 107W/cm2 ~--— 6 ><1o7W/cm2
,z“ 10 , ------ 8 X 107W/cm2
'u) 2.5X10 ‘ ./’I
(U 10 "I ,o I.” ,- r .
m 2 .0X1 0 _, ‘," fir"
C 'II Iii] .'/.'
.Q ’ ; .- /.»’
H 10 .' :I ‘ ".
Q ~—1 0 I .’
.E 1 5x10 5 f X
D . I' I; I" (I
g 10 10
O . X10 d 1‘ I 'l I,
Q) : I. -’ I
tr ,
9 .‘ / I
SCMIO « ,
‘l
00 T -—-;-'-"

 

 

 

T

 

100 150 200 250 300

Time (ps)

Figure 42 Evolution of recombination coefﬁcient under different laser intensities.
The recombination coefﬁcient increases due to the generation of new free electrons
and the decrease of the electron temperature.
similar pattern as the electron temperature. The pattern of evolution of the recombination
coefﬁcient can be explained by Eq. (134). Note that there are also “saturation”
phenomena for the recombination coefﬁcient. That happens when the electron number
density and the electron temperature approach the “static” state as shown in Figure 39

and Figure 40.

Figure 43 and Figure 44 show the temperature evolutions of ion and neutral
respectively. The ions are also heated by the laser ﬁeld directly. But the heating rate is
much slower because the ions are much heavier that the electrons. The ion temperature
also decreases after some time as do the electrons. However, the cause of this temperature
decrease is mainly due to the attenuation of the electric ﬁeld instead of the cooling

scheme of the ionization. That is why the decreasing points of the ion temperature are

136

much later than those of electron temperature. Ofcourse, the attenuation of the electric
ﬁeld, which will be shown later, also affects the heating of electrons. But it is believed to
be the minor factor comparing to the effects of the ionization for electrons. Unlike the
charged particles, the neutrals are only heated by the heat transferred from the charged
particles. That explains the latest increasing time of the neutral temperature as shown in
Figure 44 and its continuous increases with time: as long as there is temperature
difference between charged particles and neutrals, the cross heat transfer will be

continued and the neutrals will be kept heated.

 

 

 

 

 

4.0
.- —— 1 x 107W/cm2
3.5~ '1 --- 2 x107W/cm2
V, t - - 4 x 107 Wi’cm2
- - 6 x 107W/cm2
3‘0 “ .-',- _. ------- 8 x 107W/01112
A i. ,’ \\ “_
a i ..
2 2 5 — 1 I" \ ‘qiiie
.E " x‘l‘t
II I ‘3‘. " ..
'- 2.0 4 'I ’ ‘~ ‘- """
1-0 1 1 1 r 1 r 1 r T ' T m j
0 50 100 150 200 250 300
Time (ps)

Figure 43 Evolution of maximum ion temperature with different laser intensities.
Compared to the electron temperature, the ion temperature increases with a lower
rate because of the heavy mass. It decreases due to the attenuation of the electric
ﬁeld.

137

 

2.2

 

 

 

2 1 ‘ ——- 1 x 107W1’cm2 - - — 2 x 107W/cm2
- T 7 2 7 , 2
. - - 4 x 10 W/cm ------ 6 >110 W1cm ‘
2’0: ------ 8 *1 107W/cm2 , r'
1.9 "‘ z I
A 1.8 - , ,
E 17 "" '1", I”, I",
v ‘1 ‘," z I
a; 1.6 . . ,
E i i. I” "i ”
C 1.5 ‘1 "‘I I" .’v’ I I I
1.4 "‘ 3’ ,’ I 2 ’ I
+ I, "I f, I”
1.3" ‘3. I. (I! Iz”
1.2 " ’4'. 1" j", I I ’ I I
11 _i 414;, .r - ' ”_’;___———-/
I ' I ' 1 if I ' I T T I
O 50 100 150 200 250 300
Time (ps)

Figure 44 Evolution of neutral temperature with different laser intensities. Unlike
the temperatures of electrons and ions, the neutral temperature keeps increasing
because the neutrals are heated by the cross-heat transfer from the charged
particles.

Figure 45 and Figure 46 show the distributions of electron temperature and

electron number density along the beam propagation axis respectively where the incident

laser intensity is 2x 107 W/cm2 . It can be seen from the ﬁgures that, very naturally, the

electrons are heated ﬁrst near the interface and thus the largest electron number density
always emerges at the interface. As analyzed above, after some time, the temperature
near the interface decreases due to the ionization. The temperature distribution on the axis
becomes more and more uniform with time lapse. Finally, an almost uniform temperature
is established along the axis and it decreases very slowly which is believed to be due to
the temperature diffusion process. The saturation of electron number density can also be
seen from Figure 46 where the number density near the interface increase very slowly

with time in the later stage.

138

 

20«
l

1-

161

-l

124
S
3’.
(D

l— 8~

 

y...

." _." '
(-—-_--_- —.q‘.‘--vv-‘.-———o~‘---rc v.-.-‘,‘_o 3’;

I

‘ . v

F- .

/'
/'
I.
/'

   

——486ps

Incident
Laser Beam

 

   

1326ps

 

 

 

4 .1 — — - 3O 05 ps
- - 49 93 ps
------ 9987 ps
149 80 ps
0 I T T I I T I ' I '
O 20 4O 60 80 100 120 140
y (11m)

Figure 45 Distribution of electron temperature along beam propagation axis at
different time moments. The highest temperature always appears at the interface.

 

 

 

 

 

 

10 3 — 30.05 ps 3
‘ - - _, 49.93 ps
----- 99.87 ps :1
- ----- 149.80 ps i:
------- 250.12 ps 1*
1 1 5': Incident
A i j : Laser Beam
(1) 4 / l
C / I
0.1 1 // , ,:
j / I I l
. ,/ , , |
0.01 .
I I I r I T I T I
0 20 4O 60 80 100 120 140
y (11m)

Figure 46 Distribution of electron number density along the beam propagation axis
at different time moments. Most of the new free electrons are generated near the

interface.

139

 

(a) t = 49.93 ps

 

(b) t = 149.80 ps

Figure 47 Snapshots of electron and neutral number densities at (a) t = 49.93 ps and
(b) t = 149.80 ps

140

Figure 47 shows the contours of electron and neutral number densities at two time
moments where the unit ofthe number density in the plot is CHI—3. The laser intensity in

this illustration is 2x107 W/cm2 . The generation of free electrons and loss of neutral

particles are matched very well from the plots, which is evident by the ionization terms
introduced in the Botzmann equations. Further, it is clearly that more electrons are
generated along the axis. That is understandable because the laser intensity is higher
along the axis due to its Gaussian distribution in space. It also can be seen that most of
the free electrons are generated in the ﬂat thin layer below the interface as shown in
Figure 47 (b). This can be explained as follows. The electrons on the axis are heated at
the earliest and when their energy is high enough, they will ionize the neutrals due to the

electron impact ionization. On the other hand, the electrons at the outer rim ofthe laser
ﬁeld will not be heated so quickly because the laser intensity there is not as high as in the
center. This leads to the highest temperature on the axis at the early stage as shown in
Figure 48 (a), (b) and (c) where the temperature is in unit of eV. However, because the
ionization rate is higher at the center, the temperature decreases also quicker there due to
the ionization. In such a case, it is possible that the temperature at the outer rim of the
laser ﬁeld is higher after some time as shown in Figure 48 ((1). As a result, more electrons
are generated at those locations than on the axis. Then a ﬂat pattem of the electron
number density below the interface emerges. This can be seen more clearly from Figure
49 which is the electron number density distribution along the x-axis atj = 195, i.e. 5 grid

points below the interface. The above analysis can be veriﬁed more evidently by Figure

50 which shows the snapshots of ionization coefﬁcient Rf. (s_1 ) and recombination

141

   

(a)t=4.86 ps (b)t=13.26ps

   

(c) t = 49.93 ps (b) t = 149.80 ps
Figure 48 Contours of electron temperature at (a) t = 4.86 ps, (b) t = 13.26 ps, (c) t=
49.93 ps and (d) 149.80 ps
coefﬁcient R; (5‘1) at t = 149.80 ps. It can be seen that the ionization becomes stronger
at the outer rim of the laser beam below the interface. At the same positions, the
recombination is weaker while it is more prominent in the center. This leads to relatively
less electron generated in the center and more generation of electrons at the outer rim of

the laser beam.

142

ne <%>

 

 

 

 

1.6
‘ ,- .2 — 3005 ps
1-4‘ y“. --- 4993ps
1 —-—~- 99 87 ps
1.2 — ,‘ 1‘- — - 149,80 ps
. ,' i ------- 25012 ps
1.0 -* i I-\ '.
_ 1 ' \ .
0.8- ' 1 ‘, '.
. i i '._ |
0.6“ I _’ l“ 1
- J 1" 1 l
0.4- ‘- :1, l":
' ‘ 1' ‘ '.
0.2 ‘ 3' [I ‘. ‘I
‘ II, I _ \ 1. '.
4' ‘3»,
0.0 , . . “a”? . 3 . .
-60 -40 -20 0 20 40 60
X (11m)

Figure 49 Electron number density distribution along x-axis atj = 195 at different
time moments

 

(a)

 

(b)

Figure 50 Snapshot of (a) impact ionization coefﬁcient and (b) recombination
coefﬁcient at t = 149.80 ps

I43

The temperature evolution as shown in Figure 48 can be explained from the point
view ofJoule heating, which is essentially the laser energy deposited into the plasma. The
snapshots ofJoule heating at four time moments are shown in Figure 51. It is clearly seen
that the Joule heating is attenuated along the laser beam propagation and most ofthe laser
energy is deposited in the thin layer below the interface eventually. According to the
deﬁnition, Joule heating is dependent on two factors: the current density and the electric
ﬁeld. After a further analysis, one can ﬁnd that the Joule heating is essentially
proportional to the number densities ofthe charged particles and the square ofthe electric
ﬁeld while inversely proportional to the electron temperature. Thus, it is natural to
conclude that the Joule heating decreases in the region where the number densities of
charged particles are low and the electron temperature is high, just as shown in Figure 51.
For the electric ﬁeld, with the increase of the current density (mainly due to the increase
of number densities of charged particles), the electric ﬁeld is reduced according to the
Maxwell’s equation. It is apparent that the current density is the highest on the axis below
the interface because the strongest ionization takes place there. Thus, the electric ﬁeld is
attenuated most in that region as shown in Figure 52. Comprehensively considering the
effects of the current density and the electric ﬁeld, the Joule heat shows the “leaking”
phenomenon as shown in Figure 51. This also explains why there is similar phenomena

in temperature evolution.

144

   

(a) t = 4.86 ps

   

(c) t = 49.93 ps (d) t = 149.80 ps

Figure 51 Snapshots of Joule heating at (a) t= 4.86 ps, (b) t = 13.26 ps, (c) t = 49.93
ps and (d) t = 149.80 ps

I45

 

 

 

 

(a) t = 4.86 ps

 

 

 

 

 

 

(c) t = 49.93 ps (b) t = 149.80 ps

Figure 52 Snapshots of x-component of electric ﬁeld at (a) t= 4.86 ps, (b) t = 13.26
ps, (c) t = 49.93 ps and (d) t= 149.80 ps

146

Finally, the evolutions of number densities, temperatures, ionization and
recombination coefﬁcients, and magnitude of electric ﬁeld at Point A (the center point

just below the interface whose coordinate is (128, 199) in the grid system) are shown in

Figure 53. Those results are obtained by using the laser intensity of 2x107 W/cm2 .

Generation of electrons and loss of neutrals happen correspondingly as illustrated by
Figure 53 (a) where the saturation phenomenon can be shown very clearly. This
saturation can be explained directly by Figure 53 (c). In that plot, the ionization
coefﬁcient reaches a “static” state due to the decrease ofthe electron energy. As analyzed
above, this energy loss is caused by the ionization itself because ionization here also
plays a role for cooling the electrons. Meanwhile, the recombination coefﬁcient increases
thanks to the combined effects of electron number density and temperature. Moreover,
recombination of charged particles also releases heat which can contribute to increase the
electron temperature. Eventually, the recombination coefﬁcient approaches to the
ionization coefﬁcient as shown in the plot. Then the equilibrium between the ionization
and recombination is expected to be established. The electron temperature and ion
temperature are shown in Figure 53 (b). The different heating rate of electrons and ions
caused by their far different masses are apparent from the plot. The electron temperature
increases very quickly at the early stage mainly due to the rapid heating from the laser
energy through Joule heating in this model. After the effective ionization takes place, the
electron temperature starts to decrease due to the dynamics of ionization and
recombination as analyzed before. Like the electrons, the ions are also heated by the

incident laser beam but with a much lower heating rate. Decrease of the ion temperature

147

is mainly caused by the attenuation of the electric ﬁeld as illustrated in Figure 53 (d).
This attenuation of electric ﬁeld is deduced by the increased current density which is
dependent on the temperatures and number densities of charged particles. Since both the
temperatures and number densities illustrate the saturation behavior, the attenuation of

the electric ﬁeld is also presenting the similar pattern as expected.

 

 

 

 

 

 

 

 

 

 

 

 

25 .
iool-—-- , I 12.7
95' F‘
g? 90.1 ..... _ 9 201/1 \\\ 1'24
V ‘ (D r, ‘ _____
2%“ 85 —— ne(°/o) E 15]" ' _____ Te (eV) 121
C ‘.
w 80? ~~~~~~~ n % i 3 l
o ' “( l e l .1 ------- T1(eV) 11.8
5 151 8.101 .
D 10: M" E i1 :. M_..‘ l
E 1 .9 """"" ‘1 5
0+——-~ 12
ya“..— . - a r- f T A Oi - r . Y? - ~——w—-——-—.——.—.—-——¥
0 50 100 150 200 250 300 0 50 100 150 200 250 300
Time (ps) Time (ps)
(:1) (b)
8X10101 1.2x107“
exio‘Ol R|(( 4)) :5 ”mo 1 ' 0 IE] (V/m)
.— —— S > 5
ml T v 8.0x10
if); 5X101O‘. % 61
0t 4x10 . 1: some)
cf 3x10mi .g s
. 1’5 4.0x10
2x10101 ‘ .......2..4_‘_m_~Tuning; f, l
”1010,. “7” 2.0x106i
l
0L:.—,—-.-—7/’{ - . . . 0.0l . - - . . . - .
0 50 100 150 200 250 300 0 50 100 150 200 250 300
Time (ps) Time (ps)
(C) (d)

Figure 53 Evolution of (a) number densities (b) temperatures (c) ionization and
recombination coefﬁcients and (d) magnitude of electric ﬁeld at Point A

148

In conclusion, the interaction between laser and weakly ionized plasma shows
very complicated behaviors due to the dynamics of ionization and recombination. The
propagation of the laser beam and the laser energy deposition are signiﬁcantly affected
by the generation of the free electrons. Meanwhile, the heating pattern of free electrons
by the laser beam also determines the ionization and recombination dynamics of the
plasma. The spatial distribution of the laser ﬁeld makes the problem more complicated.
The simulation results presented above show the great potential of the current model to
simulate laser interaction with weakly ionized plasma. Not only the wave characteristics
of the laser beam are retained, but also the physical heating patterns for all the plasma
particles are illustrated by using this model. Unlike the kinetic-based model, such as PIC,
this model can be applied for physical problems with much larger timing and spatial
scales. Comparing to the hydrodynamic models, the present model has less assumptions
(such as continuum assumption) and thus can be used to simulate a wide range of

problems.

149

Chapter 6

Conclusion

A new LBM based model is developed in this dissertation, for the ﬁrst time, to
simulate laser interaction with weakly ionized plasmas. By introducing an additional tem1
in the Boltzmann equations, the effects of electron impact ionization and its inverse
process, three-body recombination, on the change ofrate ofthe distribution functions are
included so that the ionization and recombination dynamics can be taken into account in
this model. The Boltzmann equations, particle energy equations and Maxwell’s equations
are coupled and constitute the governing equations ofthe present model. The Boltzmann
equations are solved by the traditional D2Q9 LBM on a uniform square grid system. The
external force tem1 in the Boltzmann equation is evaluated by using He’s model after an
elaborate comparison of the results of some hydrodynamics ﬁows obtained by the other
four models. In the evaluation ofthe extemal force term and the Joule heating term in the
energy equation, the FDTD method is adopted to solve the full set of Maxwell’s equation
so that the wave nature of the electromagnetic ﬁeld carried by the laser beam and
generated by the plasma dynamics can be retained to the greatest content. The energy
(temperature) evolutions of the individual species of particles in the plasma are solved
from the energy equations which are directly derived from the Boltzmann equations by
taking proper moments ofthe distribution functions.

To match the physical properties ofthe plasma, a rescaling scheme is preposed in

this study to select the proper simulation parameters used in LBM. This rescaling scheme

150

is validated by the two-dimensional driven cavity ﬂow and the isothermal weakly ionized
plasma simulations. In the practical implementation ofthe present model, the spacing Ax
is usually determined by the resolution requirements of the physical problem, the time
step AI is selected by the stability requirement of the F DTD solver and the dimensionless

relaxation time r can be prescribed arbitrarily as long as it is in the valid range. Once the

above three parameters are determined, the lattice sound speed 61. , the lattice relaxation

time Jim (5 = 6,1, 11 ) and the lattice acceleration term is. (s = 6,1 ) can be found according

to the rescaling rules and would be used in the LBM calculations. After the standard
collision step and streaming step, it is necessary to introduce an interpolation step to ﬁnd
the on-node values ofthe distribution functions. This is due to the fact that the particles in
the plasma cannot travel to their neighboring nodes in the period of a time step At.
Conclusively, the Boltzmann equations with the ionization and recombination term are
solved by the standard D2Q9 LBM with supplement of the rescaling scheme and an
interpolation step following the streaming step. The number densities, macroscopic
velocities and current densities of the plasma particles can be found as the results of the

LBM solver and are used in the solvers of the FDTD and the energy equations.
As a simulation example, the interaction between a continuous C02 laser beam

and a weakly ionized helium plasma is simulated by using the proposed model in this
dissertation. The results show the physical heating pattems ofthe particles in the plasma
due to the laser energy deposition. The ionization and recombination dynamics and its
effects 011 the heating and generation of free electrons are also illustrated from the
simulation results. Although the simulation conducted in this dissertation is only a

preliminary application, it still could show the capabilities and limitations of the present

151

model. Compared to the kinetic-based methods, such as PIC models and MD simulations,
our model is much less computationally expansive because there is no need to solve the
Newton-Lorentz equation for every single particle. Thus it is possible to use the “real”
plasma particles in the simulation instead ofthose “super-particles” used in PIC and MD
simulations. The advantage of the present model over the hydrodynamics models is that
there is no continuum assumption employed in this model. As a consequence, this model
can be applied in a wide range ofthe physical problems. Another prominent feature ofthe
present model is that Maxwell’s equation is directly solved here to retain the most
detailed wave behaviors of the electromagnetic ﬁeld. Thus, this model is applicable for
even the most complicated wave related problems, such as complex geometries, multiple
incident laser beams and ablation of gradient materials. Finally, the present model is very
easy to implement because each of the three solvers is well developed and many
achievements in the individual ﬁeld can be inherited in development ofthis model.

Of course, no “general” model which can be applied for all physical problems
exists. The limitations ofthe present model include the following aspects. First of all, this
model is only applicable for relatively low laser intensity. This limitation is mainly due to
the low Mach number assumption made in the derivation of the lattice Boltzmann
equations. High laser intensity, as expected, will induce high electron velocities which
might break the low Mach number assumption. To overcome this problem, it is intuitive
to abandon the low Mach number assumption and adopt more velocity components in
discretization of the phase space. But unfortunately, this multiple-speed approach for
high Mach ﬂow is still under development and always suffers from high instability.

Therefore, a high Mach number LBM scheme is highly desired for plasma simulation.

152

Second, the physical spatial and temporal scales that can be simulated by this model are
largely limited by the FDTD method. In order to capture the wave nature of the
electromagnetic ﬁeld accurately, the spacing Ar in FDTD is required to be much less
than the wave length of the laser beam. Further, due to the stability requirement, the time
step At also needs to be sufﬁciently small. Thus, an ultra-stable FDTD solver is needed
for simulation of real size physical problems. Third, an interface tracking technique may
be needed to simulate the plasma plume expansion in vacuum. Plasma plume expansion
is very important in many applications, such as laser welding and pulsed laser deposition.
However, one difﬁculty in simulating this phenomenon is how to track the evolution of
the interface between the plasma and the vacuum. One possible approach is to use the
level set method. But apparently, there is a long way to go to combine the level set
method with the LBM.

In conclusion, although the LBM model developed in this dissertation is only in
its preliminary stage, it has already shown great potential for simulation of laser
interaction with weakly ionized plasmas. It is hoped that those aforementioned limitations
could be eliminated in the future and as a result, this new model could be applied widely

and successfully.

153

10.

11.

l3.

l4.

REFERENCE

Du, D., et al., Laser-Induced Breakdown by Impact Ionization in Si02 with Pulse
Widths from 7 Ns to 150 Fs. Applied Physics Letters, 1994. 64(23): p. 3071-3073.

Perry, M.D., et al., Ultrashort-pulse laser machining of dielectric materials.
Joumal of Applied Physics, 1999. 85(9): p. 6803-6810.

Fan, C .H. and J .P. Longtin, Modeling optical breakdown in dielectrics during
ultrafast laser processing. Applied Optics, 2001. 40(18): p. 3124-3131.

Schlessinger, L. and J. Wright, Inverse—Bremsstrahlung Absorption Rate in an
Intense Laser Field. Physical Review A, 1979. 20(5): p. 1934-1945.

Wierling, A., et al., Inverse bremsstrahlung of hot, weakly coupled plasmas.
Physics of Plasmas, 2001. 8(8): p. 3810-3819.

Leemans, W.P., et al., Stimulated C ompton-Scattering from Preformed
Underdense Plasmas. Physical Review Letters, 1991. 67(1 1): p. 1434-1437.

Kosarev, I.N., Kinetic theory of plasma and gas. Interaction of high-intensity
laser pulses with plasmas. Physics-Uspekhi, 2006. 49(12): p. 1239-1252.

Gibbon, P. and E. Forster, Short-pulse laser-plasma interactions. Plasma Physics
and Controlled Fusion, 1996. 38(6): p. 769-793.

Rozmus, W. and VT Tikhonchuk, Skin Effect and Interaction of Short Laser-
Pulses with Dense-Plasmas. Physical Review A, 1990. 42(12): p. 7401-7412.

Pert, G.J., Inverse Bremsstrahlung in Strong Radiation-Fields at Low-
Temperatures. Physical Review E, 1995. 51(5): p. 4778-4789.

Offenberger, A.A., et al., 2-Plasmon Decay in a CoZ-Laser-Plasma Interaction
Experiment. Physical Review A, 1978. 18(2): p. 746-749.

Machacek, AC. and J .S. Wark, A ﬂuid—kinetic model for the two plasmon decay
instability. Physics of Plasmas, 2001. 8(10): p. 4357-4366.

Berger, R.L., et al., On the dominant and subdominant behavior of stimulated
Raman and Brillouin scattering driven by nonuniform laser beams. Physics of
Plasmas, 1998. 5(12): p. 4337-4356.

Lindl, J .D., et al., The physics basis for ignition using indirect-drive targets on the
National Ignition Facility. Physics of Plasmas, 2004. 11(2): p. 339-491.

154

15.

l6.

17.

18.

19.

20.

21.

[\J
[\J

23.

24.

25.

26.

27.

Kruer, W.L., The physics of laser plasma interactions. 1988, Redwood City,
Calif: Addison-Wesley. xviii, 182 p.

Maximov, A.V., et al., Modeling of stimulated Brillouin scattering near the
critical-density surface in the plasmas of direct-drive inertial confinement fusion
targets. Physics of Plasmas, 2004. 11(6): p. 2994-3000.

Shukla, P.K., et al., Relativistic Nonlinear Effects in Plasmas. Physics Reports-
Review Section ofPhysics Letters, 1986. 138(1-2): p. 1-149.

Shanna, P., et al., Laser beam filamentation and stochastic electron heating at
upper hybrid layer. Physics of Plasmas, 2008. 15(4): p. 042105.

Max, C .E., J. Arons, and AB. Langdon, Self-Modulation and Self-Focusing of
Electromagnetic- Waves in Plasmas. Physical Review Letters, 1974. 33(4): p. 209-
212.

Antonsen, T.M. and P. Mora, Self-Focusing and Raman-Scattering of Laser-
Pulses in Tenuous Plasmas. Physics of Fluids B-Plasma Physics, 1993. 5(5): p.
1440-1452.

Brunel, F., Not-So-Resonant, Resonant Absorption. Physical Review Letters, 1987.
59(1): p. 52-55.

Bauer, D. and P. Mulser, Vacuum heating versus skin layer absorption of intense
femtosecond laser pulses. Physics of Plasmas, 2007. 14(2): p. 023301.

Kruer, W.L. and K. Estabrook, be Heating by Very Intense Laser-Light. Physics
ofFluids, 1985. 28(1): p. 430-432.

Denavit, J., Absorption of High-Intensity Subpicosecond Lasers on Solid Density
Targets. Physical Review Letters, 1992. 69(21): p. 3052-3055.

Leemans, W.P., et al., Experiments and Simulations of T unnel-Ionizea’ Plasmas.
Physical Review A, 1992. 46(2): p. 1091-1105.

Auguste, T., et al., Defocusing Effects of a Picosecond T erawatt Laser-Pulse in
an Underdense Plasma. Optics Communications, 1992. 89(2-4): p. 145-148.

Chen, S.Y., et al., Evolution of a plasma waveguide created during relativistic-
ponderomotive self-channeling of an intense laser pulse. Physical Review Letters,
1998. 80(12): p. 2610-2613.

Vujicic, N., et al., Low-density plasma channels generated byfemtosecond pulses.
Applied Physics B-Lasers and Optics, 2006. 82(3): p. 377-382.

Umstadter, D., Relativistic laser-plasma interactions. Journal of Physics D-
Applied Physics, 2003. 36(8): p. R151-R165.

155

30.

31.

33.

34.

35.

36.

37.

38.

39.

40.

41.

42.

43.

44.

Bogaerts, A. and Z.Y. Chen, Effect of laser parameters on laser ablation and
laser-induced plasma formation: A numerical modeling investigation.
Spectrochimica Acta Part B-Atomic Spectroscopy, 2005. 60(9-10): p. 1280-1307.

Christen, HM. and G. Eres, Recent advances in pulsed-laser deposition of
complex oxides. Journal of Physics-Condensed Matter, 2008. 20(26): p. 264005.

Dwivedi, A., Recent advances in pulsed laser ablated plasma plumes: A review.
Surface review and letters, 2007. 14(1): p. 57-69.

Chen, X. and H.X. Wang, Prediction of the laser-induced plasma characteristics
in laser welding: a new modelling approach including a simpliﬁed keyhole model.
Joumal of Physics D-Applied Physics, 2003. 36(13): p. 1634-1643.

McCrory, R.L., et al., Progress in direct-drive inertial conﬁnement fusion.
Physics of Plasmas, 2008. 15(5): p. 055503.

Badziak, J ., S. J ablonski, and J. Wolowski, Progress and prospect of fast ignition
of ICF targets. Plasma Physics and Controlled Fusion, 2007. 49(128): p. B651-
B666.

Tajima, T. and J.M. Dawson, Laser Electron—Accelerator. Physical Review
Letters, 1979. 43(4): p. 267-270.

Joshi, C., et al., Forward Raman Instability and Electron Acceleration. Physical
Review Letters, 1981. 47(18): p. 1285-1288.

Esarey, E., et al., Overview of plasma-based accelerator concepts. leee
Transactions on Plasma Science, 1996. 24(2): p. 252-288.

Bingham, R., J .T. Mendonca, and PK. Shukla, Plasma based charged-particle
accelerators. Plasma Physics and Controlled Fusion, 2004. 46(1): p. R1-R23.

Wang, PC, et al., Electron acceleration by a propagating laser pulse in vacuum.
Physics of Plasmas, 2007. 14(8): p. 083102.

Modena, A., et al., Electron Acceleration from the Breaking of Relativistic
Plasma- Waves. Nature, 1995. 377(6550): p. 606-608.

Pasquini, C., et al., Laser induced breakdown spectroscopy. Journal of the
Brazilian Chemical Society, 2007. 18(3): p. 463-512.

Babushok, V.I., et al., Double pulse laser ablation and plasma: Laser induced
breakdown spectroscopy signal enhancement. Spectrochimica Acta Part B-
Atomic Spectroscopy, 2006. 61(9): p. 999-1014.

Gordon, C.L., et al., T ime-Gated Imaging with an Ultrashort-Pulse, Laser-
Produced—Plasma X—Ray Source. Optics Letters, 1995. 20(9): p. 1056-1058.

156

45.

46.

47.

48.

49.

50.

51.

52.

53.

54.

55.

56.

57.

Tillman, C ., et al., Elemental biological imaging by differential absorption with a
laser-produced x-ray source. Joumal ofthe Optical Society of America B-Optical
Physics, 1996. 13(1): p. 209-215.

Tillman, C ., et al., Imaging Using Hard X-Rays from a Laser-Produced Plasma.
Applied Physics B-Lasers and Optics, 1995. 61(4): p. 333-338.

Zimmerman G, et al., LASNEX Codefor Inertial Confinement Fusion. Joumal of
the Optical Society of America, 1978. 68(4): p. 549.

The Lawrence Livermore laser fusion program: A status report. Annals of the
New York Academy of Sciences, 2006. 267(Third Conference on the Laser): p.
93-1 10.

Clark, D.S., SW. Haan, and JD. Salmonson, Robustness studies of ignition
targets for the National Ignition Facility in two dimensions. Physics of Plasmas,
2008. 15: p. 056305.

Li, C .K., et al., Measuring E and 8 ﬁelds in laser-produced plasmas with
nzonoenergetic proton radiography. Physical Review Letters, 2006. 97: p. 135003.

Gus'kov, S.Y., et al., Numerical modeling of a 2D wave of electron thermal
conductivity produced by a laser beam in targets of subcritical density. J oumal of
Russian Laser Research, 2000. 21(2): p. 157-167.

Craxton, RS. and R.L. Mccrory, Hydrodynamics of Thermal Self-Focusing in
Laser Plasmas. Journal of Applied Physics, 1984. 56(1): p. 108-117.

Marinak, M.M., et al., Three-dimensional HYDRA simulations of National
Ignition Facility targets. Physics of Plasmas, 2001. 8(5): p. 2275-2280.

Meezan, N.B., et al., Hydrodynamics simulations of 2 omega laser propagation in
underdense gasbag plasmas. Physics of Plasmas, 2004. 11(12): p. 5573-5579.

Schmitt, A.J. and BB. Afeyan, T ime-dependent filamentation and stimulated
Brillouin forward scattering in inertial confinement fusion plasmas. Physics of
Plasmas, 1998. 5(2): p. 503-517.

Myatt, J., et al., Modeling stimulated Brillouin scattering in the underdense
corona of a direct drive inertial conﬁnement fusion target. Physics of Plasmas,
2004. 11(7): p. 3394-3403.

Divol, L., et al., Three-ditnensional modeling of laser-plasma interaction:

Benchmarking our predictive modeling tools versus experiments. Physics of
Plasmas, 2008. 15(5): p. 056313.

Divol, L., et al., T hree-dimensional modeling of stimulated Brillouin scattering in
ignition-scale experiments. Physical Review Letters, 2008. 100(25): p. 255001.

157

59.

60.

61.

62.

63.

64.

65.

66.

67.

68.

69.

70.

71.

Myatt, J ., et al., Nonlinear propagation of a randomized laser beam through an
expanding plasma. Physical Review Letters, 2001. 87(25): p. 255003.

Riazuelo, G. and G. Bonnaud, Coherence properties of a smoothed laser beam in
a hot plasma. Physics of Plasmas, 2000. 7(10): p. 3841-3844.

Walraet, F., G. Bonnaud, and G. Riazuelo, Velocities of speck/es in a smoothed
laser beam propagating in a plasma. Physics of Plasmas, 2001. 8(11): p. 4717-
4720.

Weber, 8., et al., Modeling of laser—plasma interaction on hydrodynamic scales:
Physics development and comparison with experiments. Laser and Particle Beams,
2004. 22(2): p. 189-195.

Lewis, K., G. Riazuelo, and C. Labaune, Modeling of imaging diagnostics for
laser plasma interaction experiments with the code PARAX. Review of Scientiﬁc
Instruments, 2005. 76(9): p. 093502.

Marinak, M.M., et al., Nonlinear Rayleigh-Taylor evolution of a three-
dimensional multimode perturbation. Physical Review Letters, 1998. 80(20): p.
4426-4429.

Marinak, M .M., et al., T hree-dimensional simulations of Nova high growth factor
capsule implosion experiments. Physics of Plasmas, 1996. 3(5): p. 2070-2076.

Schunz, G.P., P.D. Nicolai, and M. Busquet, A nonlocal electron conduction

model for multidimensional radiation hydrodynamics codes. Physics of Plasmas,
2000. 7(10): p. 4238-4249.

Devore, C.R., Flux-Corrected Transport Techniques for Multidimensional
Compressible Magnetohydrodynamics. Journal of Computational Physics, 1991.
92(1): p. 142-160.

Hittinger, J .A.F. and MR. Dorr, Improving the capabilities ofa continuum laser
plasma interaction code. Journal of Physics: Conference Series, 2006. 46: p. 422-
432.

Bibi, F.A. and J .P. Matte, Inﬂuence of the electron distribution function shape on
nonlocal electron heat transport in laser-heated plasmas. Physical Review E,
2002. 66(6): p. 066414.

Huller, S., P. Mounaix, and D. Pesme, Numerical simulation of ﬁlamentation and
its interplay with SBS in underdense plasmas. Physica Scripta, 1996. T63: p. 151-
157.

Finke, B.R., P.D. Kapadia, and J.M. Dowden, A Fundamental Plasma Based
Model for Energy- Transfer in Laser Material Processing. Journal of Physics D-
Applied Physics, 1990. 23(6): p. 643-654.

158

72.

74.

75.

76.

77.

78.

79.

80.

81.

82.

83.

Tix, C. and G. Simon, A Transport T heoretical-Model of the Keyhole Plasma in
Penetration Laser- Welding. Journal of Physics D-Applied Physics, 1993. 26(11):
p. 2066—2074.

Kaplan, A., A Model of Deep Penetration Laser- Welding Based on Calculation of
the Keyhole Proﬁle. Joumal of Physics D-Applied Physics, 1994. 27(9): p. 1805-
1814.

Beck, M., P. Berger, and H. Hugel, The effect of plasma formation on beam
focusing in deep penetration welding with C02 lasers. Journal of Physics D-
Applied Physics, 1995. 28(12): p. 2430-2442.

Kim, K.R. and DE. F arson, C02 laser-plume interaction in materials processing.
Journal oprplied Physics, 2001. 89(1): p. 681-688.

Ki, 11., PS. Mohanty, and J. Mazumder, Modeling of laser keyhole welding: Part
I. Mathematical modeling, numerical methodology, role of recoil pressure,
multiple reﬂections, and free surface evolution. Metallurgical and Materials
Transactions a-Physical Metallurgy and Materials Science, 2002. 33(6): p. 1817-
1830.

Ki, 11., PS. Mohanty, and J. Mazumder, Modeling of laser keyhole welding: Part
II. Simulation of keyhole evolution, velocity, temperature profile, and

experimental verification. Metallurgical and Materials Transactions a-Physical
Metallurgy and Materials Science, 2002. 33(6): p. 1831-1842.

Zhang, D.M., et al., A new model of pulsed laser ablation and plasma shielding.
Physica B-Condensed Matter, 2005. 362(1-4): p. 82-87.

Tan, X.Y., et al., Ionization effect to plasma expansion study during nanosecond
pulsed laser deposition. Physics Letters A, 2007. 370(1): p. 64-69.

Mazhukin, V.I., V.V. Nossov, and l. Smurov, Modeling of plasma-controlled
sit/face evaporation and condensation of Al target under pulsed laser irradiation.

Applied Surface Science, 2007. 253(19): p. 7686-7691.

Fang, Q.Y. and X.1-1. Hu, Modeling of skin tissue ablation by nanosecond pulses

front ultraviolet to near-infrared and comparison with experimental results. leee

Journal of Quantum Electronics, 2004. 40(1): p. 69-77.

Kennedy, P.K., A F irst-Order Model for Computation of Laser-Induced
Breakdown Thresholds in Ocular and Aqueous-Media .1 . Theory. leee Journal of
Quantum Electronics, 1995. 31(12): p. 2241-2249.

Amoruso, 8., Modeling of UV pulsed-laser ablation of metallic targets. Applied
Physics a-Materials Science & Processing, 1999. 69(3): p. 323-332.

159

84.

85.

86.

87.

88.

89.

90.

91.

92.

93.

94.

95.

96.

97.

Eidmann, K., et al., Hvdrodvnamic simulation of subpicosecond laser interaction
with solid-density matter. Physical Review E, 2000. 62(1): p. 1202-1214.

Ramis, R., R. Schmalz, and J. Meyertervehn, Multi - a Computer Code for One-
Dimensional Multigronp Radiation Hydrodynamics. Computer Physics
Communications, 1988. 49(3): p. 475-505.

Itina, T.E., et al., Numerical study of ultra-short laser ablation of metals and of
laser plume dynamics. Applied Physics a-Materials Science & Processing, 2004.
79(4-6): p. 1089-1092.

Veysman, M.E., et al., Femtosecond optical diagnostics and hydrodynamic
simulation of Ag plasma created by laser irradiation of a solid target. Journal of
Physics B—Atomic Molecular and Optical Physics, 2008. 41(12): p. 125704.

J iang, L. and H.L. Tsai, Plasma modeling for ultrashort pulse laser ablation of
dielectrics. Journal of Applied Physics, 2006. 100(2): p. 0231 16.

Li, H.Y. and HS. Ki, Effect of ionization on femtosecond laser pulse interaction
with silicon. Joumal of Applied Physics, 2006. 100(10): p. 104907.

Itina, T.E., et al., Laser-generated plasma plume expansion: Combined
continuous-microscopic modeling. Physical Review E, 2002. 66(6): p. 066406.

Birdsall, C.K., Particle-in-Cell Charged-Particle Simulations, Plus Monte-Carlo
Collisions with Neutral Atoms, Pic-Mcc. leee Transactions on Plasma Science,
1991. 19(2): p. 65—85.

Birdsall, GK. and AB. Langdon, Plasma physics via computer simulation. 1985,
New York: McGraw-Hill. xxiii, 479 p.

Hockney, R.W. and J .W. Eastwood, Computer simulation using particles. 1994,
Bristol [England] ; Philadelphia: Institute of Physics Publishing. xxi, 540 p.

C hessa, P., P. Mora, and T.M. Antonsen, Numerical simulation of short laser
pulse relativistic self-focusing in underdense plasma. Physics of Plasmas, 1998.
5(9): p. 3451-3458.

Mora, P. and T.M. Antonsen, Kinetic modeling of intense, short laser pulses
propagating in tenuous plasmas. Physics of Plasmas, 1997. 4(1): p. 217-229.

Verboncoeur, J.P., A.B. Langdon, and NT. Gladd, An Object-Oriented
Electromagnetic Pic Code. Computer Physics Communications, 1995. 87(1-2): p.
199-211.

Bruhwiler, D.L., et al., Particle-in-cell simulations of plasma accelerators and
electron-neutral collisions. Physical Review Special Topics-Accelerators and
Beams, 2001. 410(10): p. art. no.-101302.

160

98.

99.

100.

101.

103.

104.

105.

106.

107.

108.

109.

110.

111.

Bruhwiler, D.L., et al., Particle-in—cell simulations of tunneling ionization effects
in plasma-based accelerators. Physics of Plasmas, 2003. 10(5): p. 2022-2030.

Huang, C ., et al., QUICKPIC: A highly eﬂicient particle-in-cell code for modeling
wakeﬁeld acceleration in plasmas. Joumal of Computational Physics, 2006.
217(2): p. 658-679.

Lefebvre, E., et al., Electron and photon production from relativistic laser-plasma
interactions. Nuclear Fusion, 2003. 43(7): p. 629—633.

d'Humieres, E., et al., Proton acceleration mechanisms in high-intensitv laser
interaction with thinfoils. Physics of Plasmas, 2005. 12(6): p. 062704.

Baton, S.D., et al., Inhibition of fast electron energy deposition due to preplasma
filling ofcone-attached targets. Physics of Plasmas, 2008. 15(4): p. 042706.

Morse, R.L. and CW. Nielson, Numerical Simulation of Weibel Instability in One
and 2 Dimensions. Physics of Fluids, 1971. 14(4): p. 830-840.

Decker, C.D., W.B. Mori, and T. Katsouleas, Particle-in-Cell Simulations of
Raman Forward Scattering from Short-Pulse High-Intensity Lasers. Physical
Review E, 1994. 50(5): p. R3338-R3341.

Fonscca, RA, et al., OSIRIS: A three-dimensional, fully relativistic particle in
cell code for modeling plasma based accelerators. Computational Science-Ices
2002, Pt lii, Proceedings, 2002. 2331: p. 342-351.

Le febvre, E., et al., Numerical modeling and applications of laser-accelerated ion
beams. Computer Physics Communications, 2007. 177(1-2): p. 60-63.

Petrov, G.M., et al., Modeling of clusters in a strong 248-nm laser field by a
three-dimensional relativistic molecular dynamic model. Physical Review E, 2005.
71(3): p. 036411.

David, N. and SM. Hooker, Effects of polarization on inverse Bremsstrahlung
heating ofa plasma. Physical Review E, 2005. 72(3): p. 036402.

David, N., D.J. Spence, and SM. Hooker, Molecular-dynamic calculation of the
inverse-bremsstrahlung heating of non-weakly-coupled plasmas. Physical Review
E, 2004. 70(5): p. 056411.

Gibbon, P., et al., Tree-code simulations of proton acceleration from laser-
irradiated wire targets. Physics of Plasmas, 2004. 11(8): p. 4032-4040.

Matte, JP. and J. Virmont, Electron Heat-Transport down Steep T emperature-
Gradients. Physical Review Letters, 1982. 49(26): p. 1936-1939.

161

114.

115.

116.

117.

118.

122.

123.

124.

Matte, J .P., et al., Electron Heat-Flow with Inverse Bremsstrahlung and Ion
Motion. Physical Review Letters, 1984. 53(15): p. 1461-1464.

Liu, J .M., Measurements of inverse bremsstrahlung absorption and non-
Maxwellian electron velocity distributions. Physical Review Letters, 1994. 72(17):
p. 2717-2720.

Matte, J .P., et al., Spectroscopic Signature of Non-Maxwellian and Nonstationarv
Effects in Plasmas Heated by Intense, Ultras/tort Laser-Pulses. Physical Review
Letters, 1994. 72(8): p. 1208-1211.

Ethier, S. and JP. Matte, Electron kinetic sinndations of solid density Al plasmas
produced by intense subpicosecond laser pulses. I. Ionization dynamics in 30
femtosecondpulses. Physics of Plasmas, 2001. 8(5): p. 1650—1658.

Marchand, R., S. Lafortune, and X. Bonnin, Average Ion Approximation for
Modeling Impurity Transport in Tokamaks. Computer Physics Communications,
1993. 76(2): p. 203-214.

Town, R.P.J., A.R. Bell, and SJ. Rose, Fokker-Planck Simulations of Short-
Pulse-Laser-Solid Experiments. Physical Review E, 1994. 50(2): p. 1413-1421.

Town, R.P.J., A.R. Bell, and SJ. Rose, Fokker-Planck Calculations with
Ionization Dynamics of Short-Pulse Laser-Solid Interactions. Physical Review
Letters, 1995. 74(6): p. 924-927.

Batani, D., et al., A cellular automaton model of laser-plasma interactions. Laser
and Particle Beams, 2001. 19(4): p. 631-642.

Gusarov, A.V., A.G. Gnedovets, and I. Smurov, Gas dynamics of laser ablation:
Influence of ambient atmosphere. Journal of Applied Physics, 2000. 88(7): p.
4352-4364.

ltina, T.E., et al., Monte Carlo simulation study of the effects of nonequilibrium

chemical reactions during pulsed laser desorption. Journal of Chemical Physics,
1997. 106(21): p. 8905-8912.

Ruhl, H., et a1. Kinetic approach to saperintense laser-solid interaction. in In
Superstrong Fields in Plasmas. 1998. Woodbury, CT: American Institute of
Physics.

Chen, S. and GD. Doolen, Lattice Boltzmann method for fluid flows. Annual
Review of Fluid Mechanics, 1998. 30: p. 329-364.

He, X.Y. and LS. Luo, A priori derivation of the lattice Boltzmann equation.
Physical Review E, 1997. 55(6): p. R6333-R6336.

162

 

126.

127.

134.

135.

136.

137.

138.

He, X.Y. and LS. Luo, Theory of the lattice Boltzmann method: From the
Boltzmann equation to the lattice Boltzmann equation. Physical Review E, 1997.

56(6): p. 6811-6817.

Gombosi, T.I., Gaskinetic theory. 1994, Cambridge [England] ; New York:
Cambridge University Press. xiv, 297.

Versteeg, H.K. and W. Malalasekera, An introduction to computational ﬂuid
dynamics .' the finite volume method. 1995, Harlow, Essex, England Longman
Scientific & Technical ;: New York : Wiley. x, 257 p.

Bhatnagar, P.L., E.P. Gross, and M. Krook, A Model for Collision Processes in
Gases .1. Small Amplitude Processes in Charged and Neutral One-Component
Systems. Physical Review, 1954. 94(3): p. 511-525.

Bittencourt, J.A., Fundamentals of plasma physics. 3rd ed. 2004, New York:
Springer. xxiii, 678 p.

Oxenius, J ., Kinetic theory of particles and photons .' theoretical foundations of
non-LTE plasma spectroscopy. 1986, Berlin ; New York: Springer-Verlag. xii,
353 p.

Epperlein, E.M., Implicit and Conservative Difference Scheme for the Fokker-
Planck Equation. Journal of Computational Physics, 1994. 112(2): p. 291-297.

Mousseau, VA. and DA. Knoll, Fully implicit kinetic solution of collisional
plasmas. Journal of Computational Physics, 1997. 136(2): p. 308-323.

Bobylev, A.V. and K. Nanbu, Theory of collision algorithms for gases and
plasmas based on the Boltzmann equation and the Landau-Fokker-Planck
equation. Physical Review E, 2000. 61(4): p. 4576-4586.

He, X.Y. and LS. Luo, Lattice Boltzmann model for the incompressible Navier-
Stokes equation. Journal of Statistical Physics, 1997. 88(3-4): p. 927-944.

Hou, S.L., et al., Simulation of Cavity Flow by the Lattice Boltzmann Method.
Journal of Computational Physics, 1995. 118(2): p. 329-347.

Chen, R.H., et al., Simulation of viscous ﬂow in a 2D-cavity by the lattice
Boltzmann method. Acta Physica Sinica, 2000. 49(4): p. 631-635.

Yu, D.Z., et al., Viscous ﬂow computations with the method of lattice Boltzmann
equation. Progress in Aerospace Sciences, 2003. 39(5): p. 329-367.

Guo, Z.L., C.G. Zheng, and BC. Shi, Discrete lattice effects on the forcing term
in the lattice Boltzmann method. Physical Review E, 2002. 65(4): p. 046308.

163

143.

144.

145.

146.

147.

148.

149.

150.

He, X.Y., et al., Analytic solutions of simple ﬂows and analysis of nonslip
boundary conditions for the lattice Boltzmann BGK model. Journal of Statistical
Physics, 1997. 87(1-2): p. 115-136.

Martys, N.S., X.W. Shan, and H.D. Chen, E valuation of the external force term in
the discrete Boltzmann equation. Physical Review E, 1998. 58(5): p. 6855-6857.

Luo, L.S., Unified theory of lattice boltzmann models for nonideal gases. Physical
Review Letters, 1998. 81(8): p. 1618-1621.

Luo, L.S., Theory of the lattice Boltzmann method: Lattice Boltzmann models for
nonideal gases. Physical Review E, 2000. 62(4): p. 4982-4996.

Ladd, A.J.C. and R. Verberg, Lattice—Boltzmann simulations of particle-ﬂuid
suspensions. Journal of Statistical Physics, 2001. 104(5-6): p. 1191-1251.

Breyiannis, G. and D. Valougeorgis, Lattice kinetic simulations in three-
dimensional magnetohydrodynamics. Physical Review E, 2004. 69(6): p.
065702(R).

Buick, J .M. and CA. Greated, Gravity in a lattice Boltzmann model. Physical
Review E, 2000. 61(5): p. 5307-5320.

He, X.Y., X.W. Shan, and GD. Doolen, Discrete Boltzmann equation model for
nonideal gases. Physical Review E, 1998. 57(1): p. R13-R16.

Nourgaliev, R.R., et al., The lattice Boltzmann equation method: theoretical
interpretation, numerics and implications. lntemational Journal of Multiphase
Flow, 2003. 29(1): p. 117-169.

Li, B.M. and D.Y. Kwok, Discrete Boltzmann equation for microﬂuidics.
Physical Review Letters, 2003. 90(12): p. 124502.

Li, B.M. and D.Y. Kwok, Lattice Boltzmann model of microﬂuidics with high
Reynolds numbers in the presence of external forces. Langmuir, 2003. 19(7): p.
3041-3048.

Gunstensen, A.K., et al., Lattice Boltzmann Model of Innniscible Fluids. Physical
Review A, 1991. 43(8): p. 4320-4327.

Swift, M.R., W.R. Osborn, and J .M. Yeomans, Lattice Boltzmann Simulation of
Nonideal Fluids. Physical Review Letters, 1995. 75(5): p. 830-833.

Grunau, D., S.Y. Chen, and K. Eggert, A Lattice Boltzmann Model for Multiphase
Fluid-Flows. Physics of Fluids a-Fluid Dynamics, 1993. 5(10): p. 2557-2562.

164

153.

154.

155.

156.

157.

158.

159.

160.

161.

162.

163.

164.

165.

166.

167.

Succi, S., The lattice Boltzmann equation for ﬂuid dynamics and beyond.
Numerical mathematics and scientiﬁc computation. 2001, Oxford, New York:
Clarendon Press ; Oxford University Press. xvi, 288 p.

Makabe, T., Plasma Electronics .' Applications in Microelectronic Device
Fabrication. 2006, Hoboken: CRC Press. 355.

Chen, S.Y., et al., Lattice Boltzmann Model for Simulation of
Magneto/it'ilt'0ilimtniics. Physical Review Letters, 1991. 67(27): p. 3776-3779.

Martinez, D.O., S.Y. Chen, and W.H. Matthaeus, Lattice Boltzmann
Magnetohydrodviuzmics. Physics of Plasmas, 1994. 1(6): p. 1850-1867.

Schaffenberger, W. and A. Hanslmeier, T wo-dimensional lattice Boltzmann model
for magnetohydrodvnamics. Physical Review E, 2002. 66(4): p. 046702.

Dellar, P.J., Lattice kinetic schemes for magnetohydrodynamics. Journal of
Computational Physics, 2002. 179(1): p. 95-126.

Orszag, SA. and CM. Tang, Small-Scale Structure of 2-Dimensional
Magnetohydrodynamic Turbulence. Journal of Fluid Mechanics, 1979. 90(JAN):
p. 129-143.

Dahlburg, RB. and J .M. Picone, Evolution of the Orszag-Tang Vortex System in a
Compressible Medium .1. Initial Average Subsonic Flow. Physics of Fluids B-
Plasma Physics, 1989. 1(11): p. 2153-2171.

Mei, R., et al., Lattice Boltzmann method for 3-D ﬂows with curved boundary.
Journal ofComputational Physics, 2000. 161(2): p. 680-699.

Grauer, R. and C. Marliani, C urrent-sheet formation in 3D ideal incompressible
magnetohydrodynamics. Physical Review Letters, 2000. 84(21): p. 4850-4853.

Finn, J .M. and PK. Kaw, C oalescence Instability of Magnetic Islands. Physics of
Fluids, 1977. 20(1): p. 72-78.

Longcope, D.W. and HR. Strauss, The Coalescence Instability and the
Development of Current Sheets in 2-Dimensional Magnetohydrodynamics.
Physics of Fluids B-Plasma Physics, 1993. 5(8): p. 2858-2869.

Marliani, C. and HR. Strauss, Reconnection of coalescing magnetic islands.
Physics of Plasmas, 1999. 6(2): p. 495.

Luo, LS. and SS. Girimaji, Lattice Boltzmann model for binary mixtures.
Physical Review E, 2002. 66(3): p. 035301.

Sofonea, V. and RF Sekerka, BGK models for diffusion in isothermal binary
fluid systems. Physica A, 2001. 299(3-4): p. 494-520.

165

168.

169.

170.

171.

173.

174.

176.

177.

178.

179.

180.

Guo, Z.L. and TS. Zhao, Finite-difference-based lattice Boltzmann model for
dense binary mixtures. Physical Review E, 2005. 71(2): p. 026701.

Xu, A.G., Finite-difference lattice-Boltznutnn methods for binary ﬂuids. Physical
Review E, 2005. 71(6): p. 066706.

Asinari, P., Viscous coupling based lattice Boltzmann model for binary mixtures.
Physics of Fluids, 2005. 17(6): p. 067102.

Sofonea, V. and RF. Sekerka, Diffuse-reﬂection boundary conditions for a
thermal lattice Boltzmann model in two dimensions: Evidence of temperature
jump and slip velocity in microchannels. Physical Review E, 2005. 71(6): p.
066709.

He, X.Y., L.S. Luo, and M. Dembo, Some progress in lattice Boltzmann
method .I. Nonunif'orm mesh grids. Journal of Computational Physics, 1996.
129(2): p. 357-363.

Guo, Z.L., C.G. Zheng, and TS. Zhao, A lattice BGK scheme with general
propagation. Journal of Scientiﬁc Computing, 2001. 16(4): p. 569.

He, X.Y. and GD. Doolen, Lattice Boltzmann method on a curvilinear coordinate
system: Vortex shedding behind a circular cylinder. Physical Review E, 1997.
56(1): p. 434-440.

He, X.Y., L.S. Luo, and M. Dembo, Some progress in the lattice Boltzmann
method: Reynolds number enhancement in simulations. Physica A, 1997. 239(1-3):
p. 276-285.

Dunweg, B., U.D. Schiller, and A.J.C. Ladd, Statistical mechanics of the
ﬂuctuating lattice Boltzmann equation. Physical Review E, 2007. 76(3): p.
036704.

Boghosian, B.A., et al., Galilean-invariant multi-speed entropic lattice Boltznumn
models. Physica D-Nonlinear Phenomena, 2004. 193(1-4): p. 169-181.

Chikatamarla, S.S., S. Ansumali, and IV. Karlin, Entropic lattice Boltzmann
models for hydrod}'namics in three dimensions. Physical Review Letters, 2006.
97(1): p. 010201.

Ghia, U., K.N. Ghia, and C .T. Shin, High-Re Solutions for Incompressible-Flow
Using the Navier Stokes Equations and a Multigrid Method. Journal of
Computational Physics, 1982. 48(3): p. 387-411.

He, X.Y. and G. Doolen, Lattice Boltzmann method on curvilinear coordinates
system: Flow around a circular cylinder. Journal of Computational Physics, 1997.
134(2): p. 306-315.

166

184.

186.

187.

188.

189.

Li. H.Y. and H. Ki, Lattice Boltzmann method for weakly ionized isothermal
plasmas. Physical Review E, 2007. 76(6): p. 066707.

Kim, Y.K. and ME. Rudd, Bimtrv—Encozinter-Dipole Model for Electron-lmpact
Ionization. Physical Review A, 1994. 50(5): p. 3954-3967.

Shah, M.B., et al., Single and Double Ionization of Helium by Electron-Impact.
Journal of Physics B-Atomic Molecular and Optical Physics, 1988. 21(15): p.
2751-2761.

Montague, R.G., M.F.A. Harrison, and A.C.H. Smith, A Measurement of the
C ross-Section for Ionization of Helium by Electron-Impact Using a Fast Crossed
Beam Technique. Joumal of Physics B-Atomic Molecular and Optical Physics,
1984. 17(16): p. 3295-3310.

Kemp, A.J., R.E.W. Pfund, and J. Meyer-ter-Vehn, Modeling ultrafast laser-
driven ionization dynamics with Monte Carlo collisional particle—in-cell
simulations. Physics ofPlasmas, 2004. 11(12): p. 5648-5657.

Taflove, A. and SC. Hagness, Computational electrodvnamics .' the finite-
difference time-domain method. 2nd ed. Artech House antennas and propagation
library. 2000, Boston: Artech House. xxiii, 852 p.

Berenger, J.P., A Perfectly Matched Layerfor the Absorption of Electromagnetic-
ll’aves. Joumal of Computational Physics, 1994. 114(2): p. 185-200.

Sacks, Z.S., et al., A perfectly matched anisotropic absorber for use as an
absorbing boundary condition. leee Transactions on Antennas and Propagation,
1995. 43(12): p. 1460-1463.

Kaganov, M.I., I.M. Lifshitz, and L.V. Tanatarov, Relaxation between Electrons
and the C rystalline Lattice. Soviet Physics Jetp-Ussr, 1957. 4(2): p. 173-178.

167