VOLUME 74, NUMBER 26 PHYSICAL REVIEW LETTERS 26 JvxE 1995
A Consistent Boltzmann Algorithm
Francis J. Alexander, Alejandro L. Garcia, * and Hemi J. AlderLawrence Livermore National Laboratory, Livermore, California 94550
(Received 14 February 1995)
The direct simulation Monte Carlo method for the Boltzmann equation is modified by an additionaldisplacement in the advection process and an enhanced collision rate in order to obtain the exact hardsphere equation of state at all densities. This leads to consistent thermodynamic and transport propertiesin the low density (Boltzmann) regime. At higher densities transport properties are comparable to thepredictions of the Enskog model. The algorithm is faster than molecular dynamics at low and moderatedensities and readily run on a parallel architecture.
PACS numbers: 05.70.Ce, 02.50.Ng, 51.10.+y
The direct simulation Monte Carlo (DSMC) method is aparticle-based, numerical scheme for solving the nonlinearBoltzmann equation for hard spheres (HS) . Ratherthan exactly calculating successive HS collisions, as inmolecular dynamics (MD) , DSMC generates collisionsstochastically with scattering rates and postco11ision ve-locity distributions determined from the kinetic theory ofa dilute HS gas. DSMC encounters the usual inconsis-tency of the Boltzmann equation; namely, it yields thetransport properties for a dilute HS gas with diameter cr,yet has an ideal gas equation of state (implying o. = 0). In this Letter a modification to DSMC is introducedwhich removes this inconsistency and recovers the exactHS equation of state at all densities with virtually no ad-ditional computational cost. This consistent Boltzmannalgorithm (CBA) has transport properties that are in simi-lar (in some cases better) agreement with HS MD thanEnskog theory .
In the standard DSMC method the positions andvelocities (r, , v;) of the particles (mass m) are evolved intime by two steps: advection and collisions. During theadvection step all particles are simultaneously propagateda distance v;6t, where the time step 6t is typically onthe order of the mean collision time. The particles aresorted into (fixed) spatial cells of dimension 6x, whichis typically on the order of A, the mean free path.Within each cell pairs of particles are then randomlyselected as possible collision partners with a HS collisionprobability that is dependent on their relative velocities.Once a pair is selected, the postcollision relative velocitiesare also stochastically determined, consistent with theconservation of momentum and energy. The collision isexecuted with the particles remaining in place.
Since in the Boltzmann equation the advection processcorresponds to that of point particles, the virial 0(Av; r;, ) is zero, giving an ideal gas equation of state(Av; is the change in velocity of particle i, and r;, isthe line connecting the centers of the colliding particles).To obtain the correct HS virial, the CBA includes theextra displacement in the advection step that the particleswould have experienced if they had collided as hard
FIG. 1. Schematic illustration of the displacement occurringafter a collision.
d =l l
o (I)V Vq
where v, = v] v2 and v,' = v& v2 are the precolli-sion and postcollision relative velocities, respectively .Particle 1 is displaced by the vector distance d and par-ticle 2 by d, as shown in Fig. 1. Equation (1) impliesthat in a one-dimensional system, when two hard rods oflength 0 collide, that after the collision the distance be-tween centers will be larger than the separation betweensimilarly colliding point particles by a distance 2cr .For hard spheres in three dimensions this effect general-izes to the displacement d above. In the low density limitthe displacement yields the HS second virial coefficientb2 = (2'/3) o.3 and hence consistency.
At higher densities the collision rate, I, in a HS gas isenhanced due to the volume occupied by the spheres ;I = AY, where A is the Boltzmann (i.e., low density)collision rate. These are, of course, functions of n, thenumber density; the Y factor is known from Monte Carloand MD simulations. When the collision rate I is used inthe CBA, the correct virial and hence, equation of state, isobtained at all densities.
Kinetic theory calculations and a series of computersimulations were carried out to obtain quantitative re-sults from this model. Pressure measured by the nor-mal momentum transfer across a plane confirmed that the
5212 0031-9007/95/74(26)/5212(4) $06.00 1995 The American Physical Society
VOLUME 74, NUMBER 26 PH YS ICAL REVIEW LETTERS 26 JUNE 1995
simulation reproduced the HS equation of state. From thehydrodynamic expression for the direct scattering func-tion, S(k, cu), the sound speed obtained from the locationof the Brillouin peak is in agreement with HS MD at lowdensities. At higher densities the Rayleigh and Brillouinpeaks are not well separated. Thus, accurate measure-ments of the sound speed cannot be made in this way.The radial distribution (pair correlation) function is thatof a perfect gas, so the sound speed determined from theequal-time density fluctuations (via the compressibility) isnot in agreement with the correct value obtained directlyfrom the equation of state.
The transport coefficients, namely, the shear viscosity(q), thermal conductivity (~), and self-diffusion (D),have been measured numerically as well as determinedanalytically from their kinetic theory expressions [10,11].They are of the Enskog form; that is, there are threeseparate contributions to the total transport coefficient akinetic (K), a potential (P), and a kinetic-potential cross(C) term. For example, for the viscosity,
rt/rto = + rt b2n + g (b2n) Y, (2)
q' =A+ B exp( st) dt, (3)where t is measured in units of the mean collision time.The term A represents the delta function contribution fromthe initial displacement (t = 0) and is proportional tothe Boltzmann average of the displacement squared fordiffusion, of the momentum Aux squared for viscosity,and of the energy Aux squared for thermal conductivity.The coefficient B is proportional to the initial decay inthe autocorrelation function determined from the nextcollision in which a common particle participates. Theintegral has the usual representation of a Markov process[10,11], where the exponential s represents the decay ofcorrelations with further stochastic collisions.
CBA calculations for the shear viscosity yield A =144/257r, or 3 times its Enskog value, and for the ther-mal conductivity, A = 64/257r, or twice its Enskog value.For shear viscosity, B = 32/25(3~3 + ~) = 0.1535,or 1.28 times its Enskog value, and for thermal con-ductivity (numerically) B = 0.104 ~ 0.004 or 0.559 ~0.022 times its Enskog value. Previous work  hasshown that the decay constant s, for the potential termis the same as the decay constant for the kinetic and crossterms, namely, s = 4/5 for shear viscosity, and s = 8/15for thermal conductivity. This leads to a potential contri-bution to g~ of 1.64 or 2.15 times the Enskog result and apotential contribution to ~ of 1.01 ~ 0.01 or 1.34 ~ 0.01
where go is the shear viscosity in the low density (Boltz-mann) approximation. The kinetic and cross contributionsare identical to those given by Enskog theory . Specifi-cally, rt =1, ~ =1, D =1, tt =4/5, ~ =6/5,and Dc = 0. The potential contributions must also be ofthe Enskog form; for example, for the viscosity,
0.4 0.6 0.8
FIG. 2. HS shear velocity as a function of number density(with kT = m = o. = 1, n* = no') . Th. e solid curve repre-sents CBA, the dashed curve Enskog theory, and the solidcircles HS MD (from ).
times the Enskog result. In these calculations the artifi-cial contribution to momentum and energy transfer due tothe finite distance separation between colliding particleswithin the same cell has been eliminated.
The shear viscosity measured in nonequilibrium Aows(Poiseuille and relaxing velocity sine wave) was inagreement with the kinetic theory results. The shearviscosity is in good agreement with both Enskog theoryand HS MD at lower densities; see Fig. 2. At the highestdensities the shear viscosity of the CBA shows betteragreement with HS MD than does Enskog theory. Forthe thermal conductivity good agreement with HS MD isfound at all densities; see Fig. 3.
The self-diffusion coefficient can be represented byD = DF + aDO, where DE is the self-diffusion coefficientfor HS in the Enskog approximation  and Dp =o.21 /6 is the self-diffusion coefficient for a randomwalker in three dimensions with jump rate I and steplength a-. The value of aD0, nonexisting in the Enskog
00 0.2 0.4 0,6
FIG. 3. HS thermal conductivity as a function of numberdensity (with kT = m = cr = 1, n* = ncr3). The solid curverepresents CBA, the dashed curve Enskog theory, and the solidcircles HS MD (from ).
VOLUME 74, NUMBER 26 PH YS ICAL REVIEW LETTERS 26 JUNE l995
0.60 0.1 0.2 0.3
0.4 0.5 0.6
0I I I I I
-8 -4 0x/1
FIG. 4. HS self-diffusion (normalized by MD self-diffusion)as a function of number density. The solid curve representsCBA and the dashed curve Enskog theory.
theory, can be derived from A = 1, B = 47r/3(3~3 +vr), and the numerically determined value of s = 0.628 ~0.001. The numerical value of a = 0.200 + 0.001 is lessthan 1 because of negative correlation between successivedisplacements (i.e. , B ( 0). This self-diffusion coefficientis in better agreement with HS MD than Enskog theoryfor number densities up to no. 3 = n* = 0.5; see Fig. 4.At higher densities the agreement fails because thedisplacement becomes of greater magnitude than the meanfree path. The self-diffusion is also too large at higherdensities because backscattering events connected withstructural effects are absent in this model (i.e., there isno "caging").
The CBA runs with nearly the same efficiency as stan-dard DSMC at low densities, since the calculation of dis-placements and the use of the Y factor only increase thecomputational cost by a few percent. At low densities HSMD is inefficient because of the large number of possiblecollision partners within a neighborhood of a mean freepath . Thus, the number of operations per collision perparticle with HS MD goes as n 2 at low densities, while itis independent of density for CBA. In comparison with ascalar HS MD code the CBA runs 2 orders of magnitudefaster for n* = 0.01414. This advantage can be further en-hanced by running on a parallel architecture .
At high densities the CBA becomes inefficient comparedwith HS MD. The reason is that a cell the size of a meanfree path, for example, one which is roughly 1/10 of a HSdiameter, represents only a small fraction (1/1000) of asingle HS particle. Thus 20 X 10 particles are requiredto represent 1000 HS particles, assuming 20 particles percell. On a single processor computer HS MD and CBA areof comparable efficiency at n* = 0.3, while on a massivelyparallel machine (with 1000 processors) this "break-even"density increases to n* = 0.7.
In conclusion, DSMC has been a popular method forthe simulation of aerodynamic Bows where conventionalNavier-Stokes solvers are inaccurate. The CBA willextend its applicability to a variety of new problems that
FIG. 5. Density (solid curves) and temperature (dashedcurves) versus normalized position for a Mach 2 shock wave.The total shock tube length is 70 mean free paths (70k), theupwind number density is n* = 0.1, and the downwind density,determined from Hugoniot conditions, is n = 0.187 for CBAand n* = 0.229 for DSMC. CBA is represented by filledsymbols, standard DSMC (with 1' factor enhanced collisionrate) open symbol.
involve moderate density gases. These include the studyof cold boundary layers in high altitude Aows and denseshocks [15,16). As an example, the normalized densityand temperature profiles for a Mach 2 shock wave arecompared in Fig. 5 to the profiles obtained from standardDSMC. The methodology described in this Letter canbe extended to more realistic intermolecular potentials byvarying the displacement as a function of density andtemperature. The implementation of such an extensionis in progress.
We thank Wm. Crutchfield, G. L. Eyink, A. J.C. Ladd,M. Malek Mansour, M. Mareschal, and T E. Wain-wright for a number of very helpful discussions and J.B.Bell for support. This work was camed out under theauspices of the Department of Energy at Lawrence Liv-ermore National Laboratory under Contract No. W-7405-EN 6-48.
*Permanent address: Department of Physics, San JoseState University, San Jose, CA 95192-0106
 G. A. Bird, Molecular Gas Dynamics and the DirectSimulation of Gas F/ows (Clarendon, Oxford, 1994).
 A. L. Garcia, Numerical Methods for Physics (PrenticeHall, Englewood Cliffs, NJ, 1994), Chap. 10.
 W. Wagner, J. Stat. Phys. 66, 1011 (1992). M. P. Allen and D. J. Tildesley, Computer Simulation of
Liquids (Clarendon, Oxford, 1987). Modified collision rates are commonly used in DSMC
to reproduce the transport properties of nonhard spheregases, e.g. , of Maxwell molecules and of the variablehard sphere model . However, these modificationsretain the ideal gas equation of state. See also F. Baras,M. M. Mansour, and A. L. Garcia, Phys. Rev. E 49, 3512(1994).
VOLUME 74, NUMBER 26 PH YS ICAL REVIEW LETTERS 26 JUNE 1995
 P. Resibois and M. De Leener, Classical Kinetic Theoryof Fluids (John Wiley and Sons, New York, 1977).
 Note that the Landau-Boltzmann equation also includesparticle interactions by modifying the advection. See L. D.Landau, Sov. Phys. JETP 3, 920 (1957); S. Grossmann, IlNuovo Cimento 37, 698 (1965).
 F.J. Alexander, A. L. Garcia, and B.J. Alder, in 25 Yearsof Non Equilibrium Statistical Mechanics, edited byJ. Brey, J. Marro, M. Rubi, and M. San Miguel, LectureNotes in Physics (Springer-verlag, Berlin, 1995).
 Moving to contact, each point particle travels an extradistance o./2 (as compared with hard rods). Moving apartafter the collision, each point particle must also travel anadditional distance o./2.
 T. E. Wainwright, J. Chem. Phys. 40, 2932 (1964).
 B.J. Alder, D. M. Gass, and T. E. Wainwright, J. Chem.Phys. 53, 3813 (1970).
 J.J. Erpenbeck and W. W. Wood, Phys. Rev. A 43, 4254(1991).
 M. Reed and K. Flurchick, Comput. Phys. Commun. 81,56 (1994).
 M. A. Fallavollita, J.D. McDonald, and D. Baganoff,Comp. Sci. Eng. 69, 269 (1992).
 E. Salomons and M. Mareschal, Phys. Rev. Lett. 69, 269(1992).
 A. Frezzotti and C. Sgarra, J. Stat. Phys. 73, 193 (1993). The Mach number for the shock wave is defined as the
upwind fluid speed divided by the sound speed in theupwind gas.