DTIC ADA444012: Evaluation of a Particle Method for the Ellipsoidal Statistical Bhatnagar-Gross-Krook Equation (PREPRINT)

0.46 MB·English

by  Defense Technical Information Center

Form Approved
REPORT DOCUMENTATION PAGE
OMB No. 0704-0188
Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and
maintaining the data needed, and completing and reviewing this collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information,
including suggestions for reducing this burden to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson Davis
Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a
collection of information if it does not display a currently valid OMB control number. PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS.
1. REPORT DATE (DD-MM-YYYY) 2. REPORT TYPE 3. DATES COVERED (From - To)
27-12-2005 Technical Paper
4. TITLE AND SUBTITLE 5a. CONTRACT NUMBER

Evaluation of a Particle Method for the Ellipsoidal Statistical Bhatnagar-Gross-Krook 5b. GRANT NUMBER

Equation (PREPRINT)
5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S) 5d. PROJECT NUMBER
Jonathan M. Burt and Iain D. Boyd (Dept. Aerospace Engineering, Univ. of MI) 5503
5e. TASK NUMBER
000P
5f. WORK UNIT NUMBER

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION
REPORT NUMBER
Air Force Research Laboratory (AFMC)
AFRL/PRSA AFRL-PR-ED-TP-2006-001
10 E. Saturn Blvd.
Edwards AFB CA 93524-7680
9. SPONSORING / MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR’S
ACRONYM(S)

Air Force Research Laboratory (AFMC)
AFRL/PRS 11. SPONSOR/MONITOR’S
5 Pollux Drive NUMBER(S)
Edwards AFB CA 93524-7048 AFRL-PR-ED-TP-2006-001
12. DISTRIBUTION / AVAILABILITY STATEMENT

Approved for public release; distribution unlimited

13. SUPPLEMENTARY NOTES
Presented at the 44th AIAA Aerospace Sciences Meeting and Exhibit, 9-12 Jan 2006, Reno, NV.
14. ABSTRACT
A particle method is presented for the numerical simulation of rarefied gas flows, based on the ellipsoidal statistical Bhatnagar-
Gross-Krook (ES-BGK) model of the Boltzmann equation. The method includes consideration of rotational nonequilibrium, and
enforces exact momentum and energy conservation for a mixture involving monatomic and diatomic species. This method is
applied to the simulation of a nozzle flow of the type associated with cold-gas spacecraft thrusters, and flowfield characteristics
are compared with experimental data as well as results from direct simulation Monte Carlo (DSMC) and Navier-Stokes
simulations of the same flow. The ES-BGK method is shown to allow for a relatively high degree of accuracy in transitional flow
regimes, while avoiding the intermolecular collision calculations which typically make the DSMC simulation of low Knudsen
number flows prohibitively expensive.
15. SUBJECT TERMS

16. SECURITY CLASSIFICATION OF: 17. LIMITATION 18. NUMBER 19a. NAME OF RESPONSIBLE
OF ABSTRACT OF PAGES PERSON
Thomas A. Smith
a. REPORT b. ABSTRACT c. THIS PAGE 19b. TELEPHONE NUMBER
A 16 (include area code)
Unclassified Unclassified Unclassified [phone]
Standard Form 298 (Rev. 8-98)

Prescribed by ANSI Std. 239.18
Evaluation of a Particle Method for the Ellipsoidal
Statistical Bhatnagar-Gross-Krook Equation (Preprint)

Jonathan M. Burt 1 and Iain D. Boyd 2
Department of Aerospace Engineering
University of Michigan, Ann Arbor, MI 48109

A particle method is presented for the numerical simulation of rarefied gas flows, based
on the ellipsoidal statistical Bhatnagar-Gross-Krook (ES-BGK) model of the Boltzmann
equation. The simulation procedure includes consideration of rotational nonequilibrium,
and enforces exact momentum and energy conservation for a mixture involving monatomic
and diatomic species. This method is applied to the simulation of a nozzle flow of the type
associated with cold-gas spacecraft thrusters, and flowfield characteristics are compared
with experimental data as well as results from direct simulation Monte Carlo (DSMC) and
Navier-Stokes simulations of the same flow. The ES-BGK method is shown to allow for a
relatively high degree of accuracy in transitional flow regimes, while avoiding the
intermolecular collision calculations which typically make the DSMC simulation of low
Knudsen number flows prohibitively expensive.

I. Introduction
I
n the design and performance analysis of low-thrust spacecraft propulsion systems, various numerical simulation
techniques may be employed to determine efficiency, thrust characteristics, or the potential for plume
impingement and contamination on spacecraft surfaces. A particular challenge in the simulation of small thrusters
involving a chemically inert neutral gas, such as electro-thermal or cold-gas thrusters, is the accurate consideration
of a wide range of Knudsen number regimes. In a typical thruster of this type, gas is expelled through a convergent-
divergent (Laval) nozzle into a near vacuum, with subsonic near-equilibrium flow in the convergent section of the
nozzle. As the gas accelerates through the divergent section beyond the nozzle throat, a subsonic boundary layer
grows along the wall, while a supersonic core-flow region exists around the nozzle centerline. The gas density
continues to decrease with downstream distance through the nozzle, and rarefaction effects become more prominent
within both the supersonic and subsonic regions. Here the gas velocity distribution begins to diverge significantly
from the equilibrium limit, and thermal energy is increasingly distributed non-uniformly among the translational and
internal degrees of freedom. Beyond the nozzle exit plane, these nonequilibrium effects continue to increase as the
gas rapidly expands and thermal energy is converted into energy associated with bulk motion of the exhaust flow.
Rotational freezing occurs due to the large gradients and low collision frequency, and the flow approaches the free
molecular limit within a short distance of the nozzle exit, particularly at points far from the nozzle centerline.
The simulation of highly rarefied flows, as described above for the divergent nozzle region and plume, is
typically performed using the direct simulation Monte Carlo (DSMC) method of Bird.1 This method approximates a
numerical solution to the Boltzmann equation – the governing equation for dilute gas flows based on a statistical
representation of molecular velocities – by decoupling in time the collision and advection terms in the equation. A
large number of particles, each representing a large number of atoms or molecules, are tracked through a
computational grid, and are sorted into cells according to their location. During each time step, some fraction of the
particles in a cell collide with each other, and probabilistic techniques are used for calculations of individual
collisions. All particles are then moved through the grid according to assigned velocities, and particles are created or
removed at inflow and outflow boundaries. Finally, macroscopic quantities are sampled by averaging various
particle properties in each cell, and the process is then repeated at the next time step.
The DSMC method has been shown to provide accurate solutions for highly rarefied nozzle and plume flows,2
and retains validity through the full range of possible Knudsen numbers for a dilute gas. In theory then, DSMC may
be applied to the entire flowfield for a small spacecraft thruster as described above. The main obstacle, however, to

1 Graduate student, AIAA student member.
2 Professor, AIAA associate fellow.
1
American Institute of Aeronautics and Astronautics
Approved for public release; distribution unlimited.
the DSMC simulation of the higher density flow near and upstream of the nozzle throat is the enormous
computational cost. There are three major factors which tend to make DSMC prohibitively expensive, even on large
parallel machines, for the simulation of relatively high-density flows: First, for accurate simulation the cell size
should ideally be much smaller than any characteristic length scale for gradients in macroscopic flow quantities, and
nonphysical diffusion effects become increasingly significant as the cell size increases. This generally limits cell
dimensions to roughly one mean free path, so that a very large number of cells may be required for the simulation of
even small-scale flows when the gas density is high. Second, the temporal decoupling of collision and advection
processes is in general a valid approximation only if the time step is no larger than the local mean collision time.
This means that for a high density or low Knudsen number flow, a very large number of time steps may be required
to reach steady state conditions. These two requirements may be relaxed somewhat when macroscopic gradient
length scales are very large compared to the mean free path, so that the flow approaches a near-homogeneous state.
However, an additional limiting factor for high density flows cannot be relaxed under near-homogeneous conditions:
The collisional process driving a flow toward equilibrium is modeled in DSMC through the probabilistic calculation
of individual collisions between pairs of representative particles. The collision frequency increases with the gas
density, so that for a fixed set of simulation parameters, a higher density flow corresponds to an increase in the
average number of collisions experienced per particle during each time step. Under conditions typical in the
convergent nozzle region and around the nozzle throat in a small electro-thermal or cold-gas spacecraft thruster, this
last requirement can make DSMC simulation of such flows prohibitively expensive on all but the largest parallel
machines.
In practice, these small-scale exhaust flow simulations may be performed using a finite volume representation of
the Navier-Stokes (NS) equations in the high density regions where the continuum and quasi-equilibrium
assumptions underlying these equations are considered valid. The DSMC method is then used only in the lower
density regions where these assumptions break down, and where a DSMC flowfield calculation is feasible using
available computational resources. The interface between NS and DSMC regions is modeled either through an
uncoupled matching of boundary conditions for two independent simulations, or using a single hybrid code where
the interface location may be varied dynamically during the simulation through the evaluation of one or more
continuum breakdown parameters. For either approach, the accurate determination of the boundary beyond which
the NS solution is invalid may be difficult and expensive. As a result, there is generally some significant tradeoff
between overall simulation accuracy and efficiency: Increased accuracy results from a more precise or conservative
determination of the proper boundary for the NS domain. This in turn requires a greater computational load, either
through an iterative procedure to evaluate the boundary then update the overall flowfield solution, or through an
extension of the DSMC domain into regions which may be accurately handled by the NS solver.
This tradeoff between accuracy and efficiency may to some extent be avoided if an alternate simulation approach
is used in flowfield regions where the validity of the Navier-Stokes equations is questionable, but the application of
DSMC is particularly expensive. This alternate method must be more efficient than DSMC in modeling relatively
high density flows, able to account for nonequilibrium phenomena which characterize a transition Knudsen number
regime out of range of NS solutions, and in good quantitative agreement with both NS and DSMC for near-
equilibrium flows. One such method is presented here, using a particle-based probabilistic solution to the ellipsoidal
statistical Bhatnagar-Gross-Krook (BGK) model of the Boltzmann equation.
In recent years, several authors have developed particle methods in which standard DSMC collision operations
are replaced by alternate procedures that account for the collisional drive toward equilibrium without modeling
individual collisions.3-9 In place of collision calculations on pairs of representative particles, some fraction of
particles are randomly selected during each time step, and are assigned new velocities according to set distribution
functions. In several of these methods, the collision process in each cell within the computational grid is based on
the homogeneous form of the BGK equation:

d(nf)/dt =νn(f − f) (1)
e

Here f is the velocity distribution function for a simple monatomic gas, f is a Maxwellian distribution with the
e
same average and variance as f, n is the number density, and ν is a characteristic relaxation frequency. This is a
simplified linear model of the highly nonlinear homogeneous Boltzmann equation, where the detailed physics of the
collision term is ignored in favor of a Maxwell-type molecular model in which the entire distribution function f
decays toward equilibrium at a velocity-independent rate. As proposed by Macrossan6 and independently
implemented by Gallis and Torczynski8 and Andries et al.,9 the method may be modified for better quantitative
agreement with the Boltzmann equation as a solution to
2
American Institute of Aeronautics and Astronautics
Approved for public release; distribution unlimited.
d(nf)/dt =νn(f − f) (2)
G

where f is an anisotropic Gaussian distribution function. This modification to the BGK equation was proposed by
G
Holway10 and Cercignani11 to correct for the nonphysical constraint in the BGK model of unity Prandtl number, and
has recently been shown to satisfy Boltzmann’s H-theorem for entropy production.12 Following Cercignani, we refer
to this as the ellipsoidal statistical BGK (ES-BGK) model.
In the first section of this paper, we extend the methods of Macrossan6 and Gallis and Torczynski8 for a gas
mixture involving monatomic and diatomic species, and describe a procedure to enforce exact momentum and
energy conservation while allowing for effects of rotational nonequilibrium. A simple homogeneous test case is then
used to verify the accuracy of the proposed rotational relaxation model. Next, the method is applied to a nozzle flow
test case involving a small cold-gas thruster with N propellant exhausting into a vacuum, and results are compared
2
with DSMC, NS, and experimental data for the same flowfield. We demonstrate relatively good agreement between
ES-BGK, DSMC and experimental results, particularly in higher density regions near the nozzle throat, and show
that the ES-BGK flowfield solution incorporates some nonequilibrium effects which are also captured by DSMC but
neglected in the NS calculation.

II. Numerical Method
As discussed above, all procedures unrelated to intermolecular collisions in the ES-BGK particle method are
identical to those in DSMC. However, the process of calculating individual representative collisions is replaced by a
procedure which involves the resampling of particle properties from predetermined distributions. Consider the ES-
BGK collision procedure in a cell containing N particles, where each particle may be assigned to one of several
different monatomic or diatomic gas species. We assume vibrational and electronic energy modes are unexcited, and
the variable hard sphere (VHS) collision model1 is used to approximate the variation of dynamic viscosity µ with
translational temperature T. As a first step in the method, T is calculated as a function of the average momentum
t t
and translational energy of all particles in the cell:

1 N
T = ( m c2 − m c 2/ m ) (3)
t 3k(N−1) a a a

Here m and c are the mass and velocity respectively of an individual particle, k is Boltzmann’s constant, the operator
designates an average over all particles in the cell, and the notation x2 for any vector x is used to represent x⋅x.
a
The coefficient N/(N–1) in Eq. (3) corrects for a statistical error which appears when the traditional formula for T
t
based on the velocity distribution is applied to a discrete set of particle velocities.13
Using a first order Chapman-Enskog expansion of the ES-BGK equation, it can be shown that the correct
relations for heat flux and shear stress are recovered in a near-equilibrium flow if the relaxation frequency is given
by

p
ν = Pr⋅ (4)
µ

where Pr is the Prandtl number and p is the pressure. From the ideal gas equation of state and the VHS
approximation that µ scales with temperature to some constant power ω, Eq. (4) may be rewritten as

⎛Tω ⎞
ν = Pr⋅nk⎜⎜⎝⎜µref ⎠⎟⎟⎟Tt1−ω (5)
ref

where T is a reference temperature and µ is the dynamic viscosity at T . For a simple VHS gas, this viscosity
ref ref ref
can be calculated as a function of ω and the molecular collision diameter d at the reference temperature:
ref

30(mkT /π)1/2
µ = ref (6)
ref 4(5−2ω)(7−2ω)d2
ref

3
American Institute of Aeronautics and Astronautics
Approved for public release; distribution unlimited.
Note that in most cases the value of d is itself determined from experimental measurements of µ . Eq. (6) is
ref ref
useful however in a hybrid code (discussed below) involving both DSMC and the method proposed here, as
reference collision diameters are among the required input parameters for DSMC collision calculations. A simple
approximate extension of Eq. (6) to a gas mixture may be obtained by summing the product of µ and the time-
ref
averaged species mass fraction in the cell over all gas species.
Once the relaxation frequency ν is known, Eq. (2) may be integrated to give the distribution function f after some
finite time step ∆t as a linear combination of the initial distribution f(0) and a corresponding anisotropic Gaussian
distribution f .
G

f(∆t)= f(0)exp(−ν∆t)+ f (1−exp(−ν∆t)) (7)
G

From Eq. (7), we determine the number of particles N selected for reassignment of velocity components as
v

N = floor[N(1−exp(−ν∆t))+R] (8)
v

where ν is evaluated using Eq. (5), R is a random number in [0,1] and the operator “floor” rounds to the nearest
smaller integer. Note that a time-averaged value of the number density n should be used in Eq. (5) to avoid an error
associated with statistical fluxuations in the calculated relaxation frequency. Once N is known, this number of
v
unique particles is then randomly selected and assigned to a list for velocity resampling. (In order to reduce calls to
the random number generator, if N > ½N then N – N randomly chosen particles are removed from a list which
v v
initially includes all particles in the cell, and the remaining particles are used for resampling.) In the special case
where N = 1, momentum conservation requires that the velocity of the selected particle remains constant, so that N
v v
values of zero and one are equivalent during the resampling procedure. To correct for the resulting error, if N is
v
calculated as one through Eq. (8) then an additional random number is generated and N is set with equal probability
v
to either zero or two.
After N particles have been selected, all velocity components of these particles are resampled from a normalized
v
Maxwellian distribution. Each velocity component is set to either

c* = cos(2πR ) −ln(R )/m or c* = sin(2πR ) −ln(R )/m (9)
i 1 2 i 1 2

where the subscript i denotes the direction, m is the particle mass, and R and R are random numbers between 0 and
1 2
1. Two statistically independent values of c* may be obtained if Eqs. (9) are evaluated for the same values of R and
i 1
R , so that only a single random number must be generated per velocity component. Following Gallis and
2
Torczynski,8 we then modify all resampled velocities to new values c(1) which conform to an anisotropic Gaussian
i
distribution:

c(1) = S c* (10)
i ij j
where
1−Pr⎡ 1 N ⎤
Sij = δij − 2Pr ⎢⎢⎣kTt N−1( mcicj a − mci a mcj a / m a)−δij⎥⎥⎦ (11)

Here δ is the Kronecker delta, c and c are particle velocity components before resampling, T is evaluated using
ij i j t
Eq. (3), and the superscript (1) indicates a newly reassigned value.
To account for the effect of collisions on the rotational energy distribution among diatomic gas species, we
define a rotational relaxation frequency ν such that, in a homogeneous gas, the instantaneous cell-averaged
r
rotational temperature T will vary as
r

dT /dt =−ν(T −T ) (12)
r r r eq

Here T is the equilibrium temperature to which both T and T converge after a long period of time, and T is given
eq r t r
by

4
American Institute of Aeronautics and Astronautics
Approved for public release; distribution unlimited.
2
T = E / ζ (13)
r k r a r a

The above symbols E and ζ denote rotational energy and rotational degrees of freedom, respectively, for an
r r
individual particle. The equilibrium temperature can be expressed as a weighted sum of the rotational and
translational temperatures, such that

ζ T + ζ T =(ζ + ζ )T (14)
t,eff t r a r t,eff r a eq

where ζ is the effective number of translational degrees of freedom for all particles in the cell. If all N particles
t,eff
are together considered as an isolated system, then it follows from Eqs. (3) and (13) and energy conservation
arguments that

d ⎡ (N-1) ⎤
dt ⎢⎣ mc 2a/ m a +3 N kTt + ζr a kTr⎦⎥ = 0 (15)

while momentum and mass conservation require that the first term within the brackets be constant in time as well.
By comparing Eqs. (14) and (15) and noting that T is a time invariant quantity, we can show that there are
eq
effectively 3(N-1)/N translational degrees of freedom in the cell, so that the equilibrium temperature may be given
from Eq. (14) as

3((N-1)/N)T + ζ T
T = t r a r (16)
eq 3(N-1)/N+ ζ
r a

Note that the arguments used to derive Eq. (16) make no assumptions regarding how the translational and
rotational temperatures relax toward the equilibrium state, so that this equation should also be valid for DSMC and
other methods for which a discrete set of particle velocities is used to describe the velocity distribution of a gas in
rotational nonequilibrium. As there are exactly three translational degrees of freedom for any real gas under
consideration here, this derivation reveals the existence of a deterministic error which scales as (N-1)/N and may
significantly influence calculated bulk flow properties when the average number of particles per cell is very small.
To the authors’ knowledge no previous references to this error source exist in the literature, although the upper
bound of the resulting error should be relatively small (less than 5%) when particle numerical weights are set such
that over 20 particles are assigned to each cell, as is standard practice in DSMC.
The rotational relaxation frequency ν shown in Eq. (12) is defined here as the frequency of inelastic collisions,
r
given as the ratio of the collision frequency F to the rotational collision number Z. To determine F , consider a
coll r coll
cell containing particles of N different species, each with a mass m, VHS reference collision diameter d and
spec i ref,i
number density n averaged over a large number of time steps. Assuming a near-equilibrium flow, we can efficiently
i
calculate the collision frequency as

4 kT ⎛ T ⎞1−ω Nspec i
Fcoll = n 2πref ⎜⎜⎜⎝Tt ⎠⎟⎟⎟ ∑ ∑ninjAij(2−δij) (17)
ref i=1 j=1
where
A = π(d +d )2 mi +mj and n = N∑specn .
ij 4 ref,i ref,j mm i
i j i=1

Note that no implied summations are used in the index notation of Eq. (17).
The rotational collision number Z is calculated using a modified form of a formula given by Parker.14
r

3/5 Z∞
Z = r (18)
r 1+(π1/2/2)(T*/T )1/2 +(π +π2/4)(T*/T )
eq eq

5
American Institute of Aeronautics and Astronautics
Approved for public release; distribution unlimited.
The above constants Z∞ and T* are determined for N as Z∞=23.0 and T*=91.5 K.1 As a modification to Parker’s
r 2 r
formula, the equilibrium temperature T is used here in place of T following observations of Boyd15 that, for a
eq t
single intermolecular collision, Z tends to vary as a function of both the translational and rotational collision
r
energies. The 3/5 factor in the numerator, also not in Parker’s original formula, accounts for the fact that Z is
r
defined here to characterize the rate at which T relaxes toward T , not toward T as used by Parker and others.1
r eq t
By analogy with Eq. (8), Eq. (12) can be integrated in time to give an initial estimate of the number of particles
N for which to resample rotational energy as
r

⎛ ⎛ F ⎞⎞
Nr ≈ N⎜⎜⎜⎝1−exp⎜⎜⎜⎝− Zcoll ∆t⎠⎟⎟⎟⎠⎟⎟⎟ (19)
r

where F and Z are evaluated using Eqs. (17) and (18) respectively. (The exact determination of N requires a
coll r r
more involved procedure, and is discussed in detail below.) In our implementation, N particles are randomly
r
selected from the list of N particles which have already been tagged for velocity resampling, and the rotational
v
energy E of each selected diatomic particle is reassigned to a normalized value
r

E(1) = −ln(R) (20)
r

where a unique random number R∈[0,1] is used for each particle.
To avoid random walk errors and permit overall accuracy for a near-equilibrium flow comparable to that of
either DSMC or NS simulations, we desire for both momentum and energy to be exactly conserved during the
resampling procedures. It follows that all reassigned velocities c(1) and rotational energies E(1) must be modified to
r
final values c(2) and E(2) such that
r

mc(2) = mc(0) (21)
v v
and
1 1
N m(c(2))2 +N E(2) = N m(c(0))2 +N E(0) (22)
2 v v r r r 2 v v r r r

where the operators and denote averages over all particles used in the resampling of velocity and
v r
rotational energy respectively. Further, to satisfy energy equipartition at the equilibrium limit, we require that
translational and rotational temperatures based on reassigned values be equal. This condition can be written as

1 Nv ( m(c(2))2 − mc(2) 2 / m )= 2 E(2) / ζ (23)
3k(N −1) v v v k r r r r
v

where ζ is the number of rotational degrees of freedom for an individual particle. For each of the N particles
r v
selected for velocity resampling, the final velocity is calculated as

c(2) = B +C(c(1) − mc(1) / m ) (24)
v
v

while each of the N particles with resampled rotational energies is given a final value of
r

E(2) = DE(1). (25)
r r

By algebraic manipulation of Eqs. (21) through (25), we can show that momentum conservation, energy
conservation and energy equipartition are all satisfied if the vector B and scalar constants C and D are chosen as
follows:

B = mc(0) / m (26)
v
v

6
American Institute of Aeronautics and Astronautics
Approved for public release; distribution unlimited.
C = ⎢⎢⎣⎢⎢⎡⎝⎜⎜⎜⎛3(Nv3−(N1)v+−N1)r ζr r⎠⎞⎟⎟⎟⎟⎜⎜⎝⎜⎜⎜⎛⎜ m(c(0))2 v −m(cm(1)c)(20)v2v−/ mmcv(1+) 2v2/(Nmr/vNv) Er(0) r⎠⎞⎟⎟⎟⎟⎟⎟⎟⎤⎥⎥⎦⎥⎥1/2 (27)

D =⎛⎜⎜⎜⎝3(Nv −N1v)+ζ Nr r ζr r⎞⎠⎟⎟⎟⎟⎛⎜⎜⎜⎜⎜⎜⎝ m(c(0))2 v − mc(0) 22v/Er(m1) vr+2(Nr/Nv) Er(0) r⎞⎠⎟⎟⎟⎟⎟⎟⎟ (28)

From the definition of rotational temperature indicated by Eq. (13), the right hand side of Eq. (28) may be
expressed as a ratio of temperatures: The denominator gives the rotational temperature based on values of E(1),
r
while the numerator gives the temperature corresponding to the conditions of Eq. (23). Thus, Eq. (28) may be used
to show that the temperature T(2) based on resampled and renormalized values of either particle velocity c(2) or
rotational energy E(2) can be given as
r

T(2) =⎛⎜⎜⎜⎝3(N −N1v)/+kN ζ ⎞⎠⎟⎟⎟⎟⎛⎜⎜⎜⎝ m(c(0))2 v − mc(0) 2v/ m v +2⎛⎜⎜⎜⎝NNr ⎞⎠⎟⎟⎟ Er(0) r⎞⎠⎟⎟⎟ (29)
v r r r v

As particles are randomly chosen for velocity or rotational energy resampling, the averaged quantities in Eq. (29)
will, on average, equal the corresponding quantities averaged over all N particles in the cell. If we denote the
equivalent temperature based on cell-averaged values as T(2), then it follows from substitution of Eqs. (3) and (13)
a
into Eq. (29) that

3((N -1)/N )T + ζ T
T(2) = v r t r a r (30)
a 3(N -1)/N + ζ
v r r a

As the temperature T(2) is generally not equal to the equilibrium temperature T determined from Eq. (16), the
a eq
procedure used to calculate N must be based on a relaxation frequency other than the true characteristic frequency
r
of rotational relaxation ν = F / Z. A correction factor β is therefore introduced as the ratio of this required
r coll r
relaxation frequency to ν, so that in a homogeneous adiabatic flow,
r

dT /dt =−βν(T −T(2)) (31)
r r r a

From a comparison of Eqs. (12) and (31), it follows that

β=(T −T ) (T −T(2)) (32)
r eq r a

Substituting Eqs. (16) and (30) into Eq. (32), we find

3+ ζ N */(N -1)
β= r a r v (33)
3+ ζ N/(N-1)
r a

where N* is an approximation to the integer value N based on an exact solution to Eq. (31). If Eq. (31) is integrated
r r
over the time step interval ∆t, then the same reasoning used to find N by Eq. (8) may be used to calculate N* as
v r

⎛ ⎛ F ⎞⎞
Nr* = N⎜⎜⎝⎜1−exp⎜⎜⎝⎜−β Zcoll ∆t⎠⎟⎟⎟⎠⎟⎟⎟ (34)
r

Note from Eqs. (33) and (34) that β and N* cannot be independently determined, but must instead be calculated
r
simultaneously using an iterative method. In the implementation proposed here, β is initially set to one and N* is
r
found by Eq. (34), after which β and N* are updated using Eqs. (33) and (34) respectively. This procedure is
r
repeated four additional times, after which the number of particles to select for rotational energy resampling is
7
American Institute of Aeronautics and Astronautics
Approved for public release; distribution unlimited.
computed as N = floor[N* +R] where R is a random number in [0,1]. This simple “brute force” approach is found to
r r
give a reasonable balance of accuracy and efficiency over a wide range of conditions, and more complicated
iterative schemes with faster convergence should not generally be required due to the insignificant contribution of
this procedure to the total computational expense.

III. Homogeneous Rotational Relaxation
To verify the accuracy of the procedure described above for rotational relaxation, this procedure is implemented
in a modified version of the DSMC code MONACO16 and applied to a zero-dimensional test case based on
simulations of Bird.1 An isolated control volume of N gas is given initial conditions involving translational
2
equilibrium at 500 K and a rotational temperature of 0 K. The simulation domain consists of a single cell with
specularly reflecting walls and 10,000 particles, and the time step is set to 1/2 of the initial mean collision time. The
rotational collision number is found from Eqs. (16) and (18) to hold a constant value of Z = 4.31, so that the
r
rotational and translational temperatures will vary in time according to analytical expressions similar to those
derived by Bird.1

⎛ ⎛ ⎞⎞ ⎛ ⎞
Tr = 300⎜⎜⎜⎝⎜⎜1−exp⎜⎜⎜⎝⎜⎜−4.31tτcoll⎠⎟⎟⎟⎟⎟⎠⎟⎟⎟⎟⎟ , Tt = 300+200exp⎝⎜⎜⎜⎜⎜−4.31tτcoll⎠⎟⎟⎟⎟⎟ (35)

Here both temperatures are given in K, t is the total elapsed time and τ is a characteristic mean collision time. The
coll
variable τ may be defined in terms of the instantaneous mean collision frequency F as
coll coll

⎡ t ⎤ -1
τcoll = t⎢⎣∫0 Fcoll(t')dt'⎥⎦ (36)

where F is calculated for a VHS gas using Eq. (17).
coll
Figure (1) shows the variation in both T and T as functions of the normalized time t/τ . Results from the ES-
r t coll
BGK simulation are presented along with results from a DSMC simulation under identical conditions, and
theoretical curves based on Eqs. (35) are also shown for comparison. For consistency, a constant rotational collision
number of 4.31 is used in the DSMC simulation. The ES-BGK data points are found to be in close agreement with
both theory and DSMC, while a slight underprediction is observed in the relaxation rate from the ES-BGK
simulation. This discrepancy can be traced to the fact that only a fraction of particles are chosen for velocity
resampling during each time step. If the artificial condition N =N is imposed, then nearly perfect agreement is found
v
between the ES-BGK and theoretical curves of Fig. (1). As discussed by Bird,1 even DSMC results for
homogeneous rotational relaxation with constant Z can only be approximated by an analytical solution, and the
r
scale of the disagreement between DSMC and theoretical curves is comparable to that between theory and ES-BGK.
Both DSMC and ES-BGK results have some dependence on the VHS temperature exponent ω, so by comparison
with DSMC we can conclude that rotational relaxation is modeled with acceptable accuracy by the method proposed
here.

IV. Nozzle Flow Simulation
Few experimental studies are available in the literature involving the measurement of flowfield characteristics
within a nozzle of the type employed in cold-gas spacecraft thrusters. To the authors’ knowledge, the most detailed
is that of Rothe,17 where an electron beam technique was used to measure density and rotational temperature for the
flow of N through a very small nozzle into a vacuum chamber. Both NS and DSMC simulations of the internal
2
nozzle flow have been performed by Chung et al.18 In the DSMC simulation by these authors, the computational
domain is restricted to the lower density regions beyond the nozzle throat, with inflow conditions at the throat based
on a NS simulation of the flow within the convergent section of the nozzle. Generally good agreement is found
between DSMC results and experimental data, while a NS simulation of the entire nozzle flow gives relatively poor
agreement with the experiment in areas downstream of the throat where nonequilibrium effects become significant.
To evaluate the accuracy of the ES-BGK approach presented here, we perform ES-BGK and DSMC simulations
for the divergent nozzle and nearfield plume flow investigated by Rothe.17 Identical grids, boundary conditions and
numerical weights are used for the two simulations, with the same grid geometry and inflow conditions at the nozzle
throat employed in the DSMC simulation of Chung et al.18 As the ES-BGK particle method is based on assumptions
8
American Institute of Aeronautics and Astronautics
Approved for public release; distribution unlimited.
of near-equilibrium flow, we can investigate the applicability of the method to higher Knudsen number regimes by
comparing ES-BGK results with both DSMC and experimental data.
Following Chung et al. we use a VHS temperature exponent (ω) of 0.24 and assume fully diffuse reflection for
particle collisions with the nozzle wall. These authors investigated a number of different models for translational
energy transfer during particle-wall collisions, and found best agreement with experimental data when all collisions
are assumed diffuse, with a 10% probability of thermal accommodation to the wall temperature of 300 K. Energy
exchange between translational and rotational modes was calculated using a detailed model which allows the
exchange probability to vary as a function of the total incident energy of the colliding particle. As noted by these
authors, a lack of experimental measurements on the interaction of N gas with solid graphite (the nozzle wall
2
material) leaves some doubt as to the validity of this energy exchange model. A more commonly accepted procedure
for modeling rotational energy transfer during particle-wall collisions in DSMC involves an assumption that particle
rotational energy is resampled from an equilibrium distribution at the wall temperature with some uniform
probability. Therefore, as a simplification to the procedure described by Chung et al., we assume that 10% of
particle-wall collisions involve full accommodation of both translational and rotational energy to the wall
temperature, and the remaining 90% of these collisions involve adiabatic diffuse reflection. In the latter type of
collision, the post-collision direction of the particle trajectory is randomly assigned according to the equilibrium
velocity distribution for particles which pass through a plane, while the incident values of both translational and
rotational energy are retained.
Axisymmetric ES-BGK and DSMC simulations are performed on a 1.4 GHz AMD Athlon cluster, with static
domain decomposition among eight separate tasks. Each calculation is run for 24 hours, with roughly 15 million
particles at steady state and an unstructured grid of about 41,000 triangular cells scaled according to the local mean
free path. The grid geometry is shown in Fig. (2) following that given by Chung et al. Not visible in the figure is the
wall curvature at the throat, where the radius of curvature is one half the throat radius of 2.55 mm. The nozzle
divergence angle is 20°, the exit-to-throat area ratio is 66, and the simulation domain extends 14 mm downstream of
the exit plane. To avoid unphysical flow characteristics beyond the exit plane due to the small domain size, we
neglect the effect of ambient pressure in the vacuum chamber and use outflow boundaries within the plume region.
Results for both ES-BGK and DSMC simulations are shown in Figs. (3) through (10).

V. Simulation Results
Mach number contours for both simulations are shown in Fig. (3), where the top half and bottom half of the
figure give ES-BGK and DSMC results, respectively. Excellent agreement between the two methods is found in the
area around the nozzle throat. In both simulations, a significant reduction in Mach number is found toward the
centerline in this region, as results from a weak compression wave due to the curvature of streamlines at the throat.18
Mach number contours along the nozzle wall are also nearly identical, and both results show a termination of the
sonic line at the nozzle lip due to the rapid acceleration expected at the lip in any highly underexpanded nozzle flow.
However, prominent differences between the two results appear along the centerline beyond a point about 1 cm
downstream of the throat. As the degree of rarefaction increases with downstream distance, the centerline Mach
number gradient in the DSMC simulation becomes significantly greater than in the ES-BGK simulation, so that the
intersection of the Mach 4 contour line with the axis occurs roughly 2.6 cm further upstream in DSMC than in ES-
BGK results. This difference in Mach number is likely a result of the unphysical lack of velocity dependence in the
collisional relaxation rate of Eq. (2), which should correspond to some reduction in the streamwise temperature
gradient within areas of low collision frequency. The difference between the results is however largely restricted to a
narrow region along the axis, and very good agreement between the two sets of contours is still found well
downstream of the throat at points further from the centerline. Note that no contour lines are shown within the plume
backflow region beyond the nozzle lip, due to the small sample size in this very low density region.
Contours of density are shown in Fig. (4), again with ES-BGK results in the top half and DSMC results in the
bottom half of the figure. All values are normalized by the stagnation density ρ = 5.323×10-3 kg/m3. As in Fig. (3)
o
we find excellent agreement around the nozzle throat and further downstream in regions far from the centerline.
Both results show an off-axis local maximum in the density profile just downstream of the throat, which is
consistent with a weak compression wave described by Chung et al.18 and mentioned above in the discussion of
Mach number contours. As above, the greatest difference between the two results is found along the axis at points
far from the throat. Here there is a noticeable increase in the ES-BGK expansion rate relative to that of DSMC, so
that at a given point on the axis the density is lower in the ES-BGK results.
This trend is also found in Fig. (5), which shows the normalized density variation along the centerline. The axial
coordinate x in the plot is normalized by the throat radius R = 2.55 mm, and a dotted vertical line gives the location
t
of the nozzle exit plane. Here results for the ES-BGK and DSMC simulations are plotted against those for DSMC
9
American Institute of Aeronautics and Astronautics
Approved for public release; distribution unlimited.