A Collective Collision
Operator
for
Michael A. Gallis and John R. Torczynski
Engineering Sciences Center
Sandia National Laboratories”
Albuquerque New Mexico 87111-0826
USA
Abstract. A new schemeto simulateelasticcollisionsin particlesimulationcodesis presented.The newschemeaimsat
simulatingthe collisions in the highly collisional regime, in which particle simulationtechniquestypically become
computationrdlyexpensive.The new schemeis based on the conceptof a grid-basedcollision field.Accordingto this
scheme, the particles perform a single coltision with the backgroundgrid during a time step. The properties of the
backgroundfield are calculatedfrom the momentsof the distributionfimctionaccumulatedon the grid. The collision
operator is based on the Langevinequation.Based on comparisonswith other methods, it is found that the Langevin
methodoverestimatesthe collisiontkquency for dilutegases.
INTRODUCTION
Particle simulation methods and in particular the Direct Simulation Monte Carlo method (DSLMC)have proved
very successful in simulating complicated physical phenomena in a raretied environment. However. many problems
of practical interest involve dense flows that make DSMC simulation computationally very: demanding, if not
impossible.
Even in the most favorable problems, particle simulation methods have the disadvantage of being
computationally expensive. The result is that the enormous potential of these methods in physical flow modeling has
not been exploited by industry. Although massively parallel computers have contributed significantly in speeding up
particle simulation codes, three-dimensional applications remain out of the scope of most modem computer
platforms at present.
The most important limitation in modeling the highdensity regime comes from the large number of collisions
that particles undergo, which adds a very significant overhead to the calculation. In the case of non-reacting flows.
these collisions keep reproducing the equilibrium distributions without adding any new information to the flow.
In a previous paper by Gallis and Torczynski (1), the BGK equation (2) and Cercignani’s extension (3) w’ere
examined as possible collision operators. In this paper, a method based on the work of Langevin (-!) to model the
collisions in a collective manner is examined and compared with the previous methods.
MODELS FOR THE COLLISION TERM OF THE BOLTZMANN EQUATION
The complexity of the nonlinear structure of the collision integml of the Boltzmarm equation has
encouraged the development of simpltiled models that retain most of the properties of the Boltzmann equation but
are not derived directly from it. The BGK method uses a simplified collision operator that takes the form:
H
~(nf )
=
n v(fO – f)
(1)
colliszon
*Sandia is a multiprogram laboratory operated by the Ssmdia Coloration, a Lockheed Martin Company, for the United States
Department of Energy under Contract DE-AC04-94AL85000.
DISCLAIMER
This report was prepared as an account of work sponsored
by an agency of the United States Government. Neither
the United States Government nor any agency thereof, nor
any of their employees, make any warranty, express or
implied, or assumes any legal liability or responsibility for
the accuracy, completeness, or usefulness of any
information, apparatus, product, or process disclosed, or
represents that its use would not infringe privately owned
rights.
Reference herein to any specific commercial
product, process, or service by trade name, trademark,
manufacturer, or otherwise does not necessarily constitute
or imply its endorsement, recommendation, or favoring by
the United States Government or any agency thereof. The
views and opinions of authors expressed herein do not
necessarily state or reflect those of the United States
Government or any agency thereof.
DISCLAIMER
Portions of this document may be illegible
in electronic image products. Images are
produced from the best available original
document.
. .
where f. is the Maxwellian distribution corresponding to the local temperature and average velocity, y is the initial
distribution of the particles, n the number density, and
v is the collision frequency. Examination of the BGK model
indicates
thatit reproduces the equilibrium solution f =f. atequilibrium. The BGK equation also reproduces the
correct moments imd satisfies the H-theorem. For the highly collisional regime, and when used with the Chapman-
Enskog expansion, it leads to the Navier-Stokes equations. In this equation, the term rrv~ represents the collisions
replenis.tig the equilibrium distribution
f., and the term mfrepresents the collisions depleting the distribution~
The physical interpretation of the assumption that the replenishing number of collisions is given by
nvfo is
that post-collision molecules are selected from a -elli~ distribution at the local average velocity and
temperature. It is evident that this behavior can readily be incorporated in a particle code. By doing so, the burden of
calculating collisions is reduced.
Despite the shortcomings of the BGK equation, it is able to reproduce the main features of the Boltzmarm
equation. Several modifications to the BGK equation have been proposed at the expense of the simplicity of the
model. Cercignani (3) proposed a modification to the BGK equation that would allow it to reproduce the same
viscosity and thermal conductivity as the fill Boltzm equatio~ here referred to as the BGKC method. To do this
the Prandtl number is introduced as a parameter in the equation. The generalization of the BGK equation is obtained
by the replacement of the Maxwellian distribution with a local anisotropic threedimensional Gaussianj referred to
as the ellipsoidal statistical (ES) model:
r. 1
[
f. = m“3/2(detA)”2 e.xp - ~ Aqu,uj
1,]=1
1
I
where the matrix A is:
Av = [(2RZ’/Pr)d,j - 2(1- Pr)p,, /(P Pr)]-l
(2)
(3)
I
where R is the universal gas constant, Z7is the velocity in the center-of-mass frame of reference. p is the local density
of the flow, Pr is the Prandtl number, and
PUis the pressure tensor. Recently this method ~vasintroduced into particle
simulations with reasomble success (1).
I
LANGEVIN’S EQUATION
Langevin proposed a method for Brownian motion to model collisions between particles (4). According to
this method, the influence of the surroundings on a particle that moves with a velocity Z7can be split up into two
parts: a deterministic part representing the dynamical fiction experienced by a particle (–A-), and a probabilistic
part characteristic of the Brownian motion (~ ). The Langevin equation takes the form:
~ = ~Oe
-@t ●fi
(4)
in
which Al is a time interval over which motion occurs and 10 is the \elocity of the particle at the beginning of the
interval. The dynamical friction coefficient P and the probabilistic part ~ are assumed to be independent of the
velocity of the particle. The probabilistic part ~ is chosen from the modified Mamvellian distribution (4):
P(Z7,At,Zio) =
~2~’’fle-2@’~2e’Tl-~f~~~:rJ “)
in which
T is local temperature, k the Boltzrnann constant, and m the mass of the particle. Effectively this
distribution is Maxwelliam at temperature T(l – e-2@t). In this study the dynamic friction coefficient is interpreted
as the equilibrium collision frequency.
It can be shown that the Langevin equation forms the basis for deriving all the physically significant
relations concerning the motion of Brownian particles. The Langevin method has been applied successfully in the
past for the modeling of charged flows with PIC by Jones et al. (5) and for DSMC by Gallis et al. (6). It should be
noted, however, that the Brownian-motion paradigm corresponds more closely to the motion of charge particles in a
plasma than to that of neutrals in rarefied
flow, More specifically, Brownian motion describes the motion of a
particle in very dense fluid (liquid density) in which multiple collisional interactions are ve~ frequent. whereas the
Boltzmann equation describes the motion of particles when only binary collisions are important. Thus, charged
flows due to the long range interaction of the Coulomb forces fit the description of the Brownian motion more
closely than neutral particles.
Another signillcant difference exists between the conventional DSMC collision operator and the Laugevin
model. In the DSMC method a fraction of the particles modi@ their properties signtilcantly in every time step while
the rest-remain unchanged, but in the Langevin model every particle in a cell receives a small modification of its
velocity in every time step. The modification constitutes a small change in their properties. It should be noted that
the random velocities are not selected from a M=wellian with the temperature
T of the background gas, but with a
‘2@ that depends on the time interval and the collision frequency. In cases where the
temperature of I’(1– e )
collision frequency is very low, the probabilistic contribution O adds a vev small vector to the particle velocity. the
deterministic contribution similarly applies a small modification to the particle velocity for evexytime step.
The Lrmgevin method assumes that the thermal energy in a cell is always distributed according to a
Maxwellian distribution. This assumption mturally restricts the method to equilibrium or near-equilibrium
situations. The worst-case scenario for the method would be the collision of two monochromatic particle beams.
While the DSMC method predicts a gradual depletion of particles from the initial distributions to an equilibrium
distribution, the Langevin model evolves all particles gmdually towards an equilibrium distribution at the average
velocity and tempemture.
APPLICATION TEST CASE
The first test case examined is that of the flow over a 10° wedge (see Figures 1 and 2). The flow travels
from Iefl to right and the sides of the wedge are assumed to be specularly reflecting. The gas is pure argon at an
ambient temperature of 180 K and a number density of 1022m-3with an upstream velocity of 1365 mk (see Table 1).
The domain is discretized with 25,000 cells, and a time step of 10-8s is used. The downstream boundary condition is
assumed to be non-reentrant (vacuum), which is acceptable for the supersonic flow downstream of an oblique shock.
This test case is simulated with the DSMC, BGKC and Langevin methods. Steady-state conditions are observed
after 5,000 time steps while the code executes a total of 65,000 time steps.
Molet
m
Number 1
Tempera~
Speed of ~
I Upstream Conditions
I
cular Mass 39.99 mu
ific Heat Ratio 5/3
)ensitv
1022m-3
ure 180 K
Sound
249.75 m/s
Free Stream Velocity 1365
dS
Mach Number 5.465
I
Downstream Conditions
I
~
Table 1. Shock wave conditions
Figure l(a-c) presents the number-density contour lines as calculated by the DSMC. BGKC and Langevin
methods. Figure 2(a-c) in the same fashion presents the corresponding translational temperature contour lines. At .x
= Y = O,the location of the leading edge of the wedge, the shock layer forms,
A quantitative comparison indicates that all three simulations are in agreement about the downstream
conditions and the shock angle. For example, with aIl three methods, the flow downstream eventually reaches a
maximum number density that is about double the upweam nuber density, and the shock angle is 19-20°, as