Transcript Direct Simulation Monte Carlo: theory, methods, & open challenges Alejandro L. Garcia San Jose State University Lawrence Berkeley Nat.

Direct Simulation Monte Carlo:
theory, methods, & open challenges
Alejandro L. Garcia
San Jose State University
Lawrence Berkeley Nat. Lab.
RTO-AVT-VKI Lecture Series, January 25th 2011
von Karman Institute for Fluid Dynamics
Physical Scales for Dilute Gases
DSMC is the dominant numerical
algorithm at the kinetic scale
T/T
Collision
Gradient Scale
Collision
Molecular
Diameter
System Size
Quantum scale
Kinetic scale
Hydrodynamic scale
DSMC applications are expanding to multi-scale problems
Molecular Dynamics for Dilute Gases
Molecular dynamics
inefficient for simulating
the kinetic scale.
Collision
Relevant time scale is
mean free time but MD
computational time step
limited by time of collision.
DSMC time step is large
because collisions are
evaluated stochastically.
Volume of potential
collision partners
Direct Simulation Monte Carlo
Development of DSMC
•
•
•
•
•
•
•
•
DSMC developed by Graeme Bird (late 60’s)
Popular in aerospace engineering (70’s)
Variants & improvements (early 80’s)
Applications in physics & chemistry (late 80’s)
Used for micro/nano-scale flows (early 90’s)
Extended to dense gases & liquids (late 90’s)
Used for granular gas simulations (early 00’s)
Multi-scale modeling of complex fluids (late 00’s)
DSMC is the dominant numerical method
for molecular simulations of dilute gases
Particle vs. Continuum
When is the continuum
description of a gas not
accurate?
Knudsen
Number
Rarefied/
Aerospace
We
are
here
Mean Free Path
=
System Length
When Kn > 0.1, continuum
description is not accurate
Microscale
Flows
Equilibrium fluctuations
are noticeable when the
number of particles < 106 .
From G.A. Bird
 – inter-atomic spacing
d – atomic diameter
DSMC Algorithm
• Initialize system with particles
• Loop over time steps
– Create particles at open
boundaries
– Move all the particles
– Process any interactions of
particle & boundaries
– Sort particles into cells
– Sample statistical values
– Select and execute random
collisions
Example: Flow past a sphere
G.A. Bird, Molecular Gas Dynamics and Direct Simulation of Gas Flows, Clarendon, Oxford (1994)
F. Alexander and A. Garcia, Computers in Physics, 11 588 (1997)
Random Numbers
Need a high-quality random number generator for the uniform
(0,1) distribution, such as Mersenne Twister.
Many distributions (e.g., Gaussian, exponential) may be
generated by the inversion method:
Generate random value x with distribution P(x) as
x = f (R) where R is uniformly distributed in (0,1).
Most other distributions are generated by the accept-reject
method:
Draw xtry uniformly in the range of x;
Accept it if P(xtry) > max{P(x)}R else draw again.
Be careful to use high-quality algorithms and be sure that you
verify your implementation with independent testing.
Initialization
Divide the system into cells and generate particles in each cell
according to desired density, fluid velocity, and temperature.
From density, determine number of particles in cell volume, N,
either rounding to nearest integer or from Poisson distribution.
Assign each particle a position in the cell, either uniformly or
from the linear distribution using the density gradient.
From fluid velocity and temperature, assign each particle a
velocity from Maxwell-Boltzmann distribution P(v; {u,T}) or
from the Chapman-Enskog distribution P( v;{u, T , u, T }).
Be careful initializing particles for initial value problems.
Ballistic Motion
Particles motion is ballistic;
during a time step, t, particle
positions are updated as,
r(t + Dt) = r(t) + v(t) t
The particles move without
interaction and can even overlap.
For transient flows, on the first time step use ½ t (Strang
splitting) to maintain accuracy. If measuring non-conserved
variables (e.g., fluxes) then time-center the sampling
(half move, sample, half move, collisions) for all steps.
Simple Boundaries
With periodic boundaries particles
are re-introduced on the opposite
side of the system.
Specular surfaces modeled by
ballistic (mirror) reflection of
point particles.
Be careful with corners.
Be careful with body forces.
Inject
Remove
Thermal Walls
A more realistic treatment of a
material surface is a thermal wall,
which resets the velocity of a
particle as a biased-Maxwellian
distribution,
uw
y
x
z
These distributions
(exponential and
Gaussian) are simple to
generate by inversion.
Walls can also be part-thermal, part-specular (accommodation).
Reservoir Boundary
Inflow/outflow boundary conditions commonly treated as a
reservoir with given density, fluid velocity, temperature.
Particles in the main system are removed if they cross the
boundary into the reservoir.
Particles injected from reservoir to main system by either:
• Surface generator: From number flux determines number
to be injected; generate particle velocities from surface
distribution (e.g., inflow Maxwellian).
• Volume generator: Initialize a “ghost cell” with particles
before the ballistic move; discard any that do not cross the
boundary into the main system during the move phase.
Traditional DSMC Collisions
• Sort particles into spatial
collision cells
• Loop over collision cells
Selected
collision
partners
– Compute collision
frequency in a cell
– Select random collision
partners within cell
– Process each collision
Collision pair with large relative velocity are more likely
to collide but they do not have to be on a collision trajectory.
Collision Rate
Number of collisions that should occur during time step is
M coll  t RN ,V , d , vr ,...
where the functional form of the collision rate, R,
depends on the intermolecular model.
For hard spheres of diameter d, the collision rate is,
RN ,V , d , vr  
where is
vr
N ( N  1) d 2 vr
2V
the average relative velocity.
Early DSMC implementations used N Nav instead of
N(N 1), where Nav is an estimated average value of N.
Collision Selections
To avoid having to compute the average relative velocity
for all particle pairs in a cell, a larger number of
attempted collisions are selected and some are rejected.
Number of collisions attempted during a time step is,
M cand  t RN ,V , d , vmax 
where vmax ≥ max{vr} is estimated maximum relative velocity.
An attempted collision is accepted with probability,
P( v i , v j ) 

R N , V , d , vr  v i  v j
R N , V , d , vmax 

Early DSMC implementations rounded Mcand down and
carried fraction to next time step. Modern approach is to
randomly round to nearest integer (or use Poisson dist.).
Post-Collision Velocities
Post-collision velocities
(6 variables) given by:
• Conservation of
momentum (3 constraints)
• Conservation of energy
(1 constraint)
• Random collision solid
angle (2 choices)
v1’
Selection of the post-collision velocities
must satisfy detailed balance.
v1
Direction of vr’
is uniformly
distributed in
the unit sphere
Vcm
vr ’
v2
vr
v2’
Random Solid Angle
Post-collision relative velocity is,
z

The azimuthal angle is just,
Polar angle distribution is,
y

But with change of variable,
x
So,
Generated by inversion method
Molecules & “Simulators”
In DSMC the number of simulation particles (“simulators”) is
typically a small fraction of the number physical molecules.
Each simulator represents Nef physical molecules.
Nef = 2
Physical Molecules
DSMC Simulators
Accuracy of DSMC goes as 1/N; for traditional DSMC
about 20 particles per collision cell is the rule-of-thumb.
DSMC “Parliament”
DSMC dynamics is correct if:
• The DSMC simulators are an
unbiased sample of the physical
population (unbiased parliament).
• Collision rate is increased by Nef so the number of collisions
per unit time for a simulator is same as for a physical molecule.
• In sampling, each simulator counts as Nef physical molecules.
Early DSMC implementations used a different representation,
rescaling the simulator diameter and mass to maintain the
same physical mean free path and mass density.
Ballistic & Collisional Transport
By their ballistic motion particles
carry mass, momentum and energy.
In a dilute gas, this is the only
source of transport.
In DSMC, momentum and energy are
also transported by the collisions.
The larger the collision cell, the more
collisional transport (greater average
separation between particle pairs).
v
vt
Cell Size and Time Step
Can calculate collisional transport by Green-Kubo theory;
error is quadratic in cell size and time step.
Collisional transport is incorrect so to minimize it the cell
size in DSMC is limited to a fraction of a mean free path.
For similar reasons, the time step is limited to a fraction of
a mean collision time.
Due to symmetry, the collisional transport does not affect
the pressure. However, if we restrict collisions to only
particles moving towards each other then this symmetry is
broken and DSMC has a non-ideal gas equation of state.
Nearest Particle Collision Partner
New DSMC implementations minimize collisional transport
error by choosing collision partner as closest particle in cell.
To avoid bias due to re-collisions, a particle pair is not
allowed to collide twice (choose next-nearest particle).
Two common implementations are:
• Transient Adaptive Sub-cells (introduced by Bird)
• Virtual Sub-cells (introduced by LeBeau, et al)
M.A. Gallis, J.R. Torczynski, D.J. Rader, G.A. Bird, J. Comp. Phys., 228 4532-48 (2009)
Fluctuations in DSMC
• Hydrodynamic fluctuations (density, temperature,
etc.) have nothing to do with Monte Carlo aspect
of DSMC.
• Variance of fluctuations in DSMC is exact at
equilibrium (due to uniform distribution for
position and Maxwell-Boltzmann for velocity).
• Time-correlations correct (at hydrodynamic scale)
• Non-equilibrium fluctuations correct (at
hydrodynamic scale)
Sampling and Fluctuations
Measurements in DSMC are done by statistical sampling.
For volume measurements
the particles are sorted into
sampling cells and polled.
vi
For surface measurements,
particles crossing a surface
during a time interval are
counted.
Many sample measurements are required due to fluctuations.
Error Bars
Error bars may be estimated by equilibrium variances
since non-equilibrium corrections are small.
In general, the standard deviation for sampled valued
goes as 1/N where N is the number of simulators.
If each simulator represents one molecule then the
fluctuations are the same as in the physical system.
If Nef >> 1 then the fluctuations are much greater than
in the physical system, with correspondingly larger
error bars.
Statistical Error (Fluid Velocity)
Fractional error in fluid velocity
Eu 
u
ux

u x2 / S
ux

1
1
SN Ma
where S is number of samples, Ma is Mach number.
For desired accuracy of Eu = 1% with N = 100 simulators/cell
1
1
S

2 2
NMa Eu Ma2
S  102 samples for Ma =1.0
S  108 samples for Ma = 0.001
(Aerospace flow)
(Microscale flow)
N. Hadjiconstantinou, A. Garcia, M. Bazant, and G. He, J. Comp. Phys. 187 274-297 (2003).
Statistical Error (Other Variables)
Fractional error in density, temperature, pressure
E 
1
SN
ET 
C
SN
C
EP 
SN
C , C   O(1)
Fractional error in temperature difference
EDT
Br

Ma 2
1
SN
where Brinkman number; Br  1, if DT due to viscous heating
For given EDT , number of samples S  1/Ma4
Statistical Error (Fluxes)
Fractional error in stress and heat flux
1
Et 
KnMa
1
SN
Br
Eq 
KnMa 2
1
SN
Comparing state versus flux variables
Eu
Et 
Kn
EDT
Eq 
Kn
Typically Kn < 0.1 so error bars for fluxes significantly
greater; measurements such as drag force are difficult for
low Mach number flows.
Variance Reduction in DSMC
Variance reduction in DSMC has been difficult to achieve,
in part because DSMC is already an importance
sampling algorithm for the Boltzmann equation.
Attempts to mollify the fluctuations in DSMC, while
preserving accuracy at kinetic scales, have been mostly
unsuccessful.
A promising approach is Hadjiconstantinou’s low-variance
algorithm, which is loosely based on DSMC.
Fluctuations and Statistical Bias
Given the presence of fluctuations in DSMC, we need
to be careful to avoid all sources of statistical bias.
Suppose we dynamically vary the cell sizes so that each cell has the
same number of particles.
This is helpful in that
DSMC is not accurate
when N is too small.
Yet there could be unintended consequences when we replace
fluctuations in N with fluctuations in V.
DSMC Collision Rate
The average number of collisions in a cell is
M  M TRY PACCEPT
Use this
N ( N  1) vr MAX t
vr


vr MAX
2V
Want this
N  vr t
N ( N  1) vrt


2V
2V
?
In general 1/ V  1/ V
so equality does not hold.
2
For Poisson
N ( N  1)  N
2
Brownian Systems
Fluctuations in DSMC are not always a nuisance; there
are interesting phenomena that rely on fluctuations.
Pressure
Hot
P
i
s
t
o
n
Pressure
Cold
“Adiabatic”
Piston
Problem
Chambers have gases at different temperatures, equal pressures.
Walls are perfectly elastic yet gases come to common temperature.
How?
Heat is conducted between the chambers by the
non-equilibrium Brownian motion of the piston.
Adiabatic Piston by DSMC
Initial State: X = L/4, M = 64 m
NL = NR = 320, TR = 3 TL
1.8
0.55
Temperatures
1.6
0.5
0.5
TR(t)
1.4
0.45
(x-piston)/L
energy per particle
Left
Right
1.2
1
0.35
0.8
TL(t)
0.6
0.4
Piston Position
0.4
0
1
2
3
time
4
Single run
5
time
6
7
8
9
20 run ensemble
0.3
10
5
x 10
0.25
0
1
2
3
time
4
5
time
6
7
8
9
10
5
x 10
Feynman’s Ratchet & Pawl
Carnot*
engine
driven by
fluctuations
Brownian
motors and
nanoscale
machines
* almost
COLD
HOT
Violate
NO. Fluctuations
also lift the pawl,
dropping the
weight back down.
pawl
nd
2
Law of Thermo?
WARM
WARM
Triangula Brownian Motor
Feynman’s complicated mechanical geometry not needed.
An asymmetrically shaped Brownian object in a non-equilibrium
system (e.g., dual-temperature distribution) is enough.
Hot gas
Cold gas
P. Meurs, C. Van den Broeck, and A. Garcia, Physical Review E 70 051109 (2004).
Algorithm Refinement
Algorithm Refinement is a multi-algorithm hybrid
methodology based on Adaptive Mesh Refinement.
At the finest level or
resolution, instead
of refining the mesh
you “refine” the
algorithm (change
to a model with
more physics).
This refinement
may be adaptive.
Particle/PDE AR Hybrid
http://cims.nyu.edu/~donev/FluctHydro/Hybrid.2D.sphere.plug.inst.mov
Particle/PDE
Algorithm
Refinement for
flow past a sphere.
Stationary particle
Molecular simulation
of solvent fluid
Interface
Continuum simulation
of solvent fluid
Note: Continuum calculation done everywhere
Advances in Algorithms Refinement
Stochastic Particle Algorithms – Our original AR method
was limited to dilute gases using the Direct Simulation
Monte Carlo scheme. Have developed more advanced
stochastic particle schemes for non-ideal fluids.
Stochastic Continuum Algorithms – Our original AR method
used a deterministic, explicit scheme for the full
Navier-Stokes equations. Have developed stochastic
PDE schemes to capture hydrodynamic fluctuations.
Coupling Issues – Perfecting the coupling of particle and
PDE schemes is challenging due to fluctuations.
Stochastic Navier-Stokes PDEs
Landau introduced fluctuations into the Navier-Stokes
equations by adding white noise fluxes of stress and heat.
Hyperbolic Fluxes
Parabolic Fluxes
Stochastic Fluxes
Fluctuating Hydrodynamic Solvers
We now have simple, accurate, and efficient finite
volume schemes for solving the stochastic NavierStokes PDEs of fluctuating hydrodynamics.
where
J.B. Bell, ALG, and S. Williams, Physical Review E 76 016708 (2007)
A. Donev, E. Vanden-Eijnden, ALG, and J. B. Bell, CAMCOS, 5(2):149-197, (2010)
Comparison with DSMC
DSMC used in testing the numerical
schemes for stochastic Navier-Stokes
Spatial
correlation
of density
Time
correlation
of density
Stochastic vs. Deterministic PDEs?
Question:
Is it necessary to
use stochastic
PDEs in the
continuum region
given that the
particle region
has fluctuations?
Answer: YES!
Simple Brownian Motion
http://cims.nyu.edu/~donev/FluctHydro/Hybrid.2D.sphere.diffusion.mov
First test is the
calculation of the
Brownian motion
of a spherical
particle.
Measure velocity
auto-correlation
function.
AR with Stochastic PDEs
Excellent agreement between a hybrid using stochastic PDE
solver and an (expensive) pure particle calculation.
A. Donev, J.B. Bell, ALG, and B. Alder, SIAM Multiscale Mod. Sim. 8 871-911 (2010).
AR with Deterministic PDEs
A hybrid using a deterministic PDE solver significantly
under-predicts the velocity auto-correlation function.
Adiabatic Piston
Cold, dense gas
PDE region
Piston
Particle region
Hot, dilute gas
PDE region
Initially the gas pressure is equal on both sides of the piston.
Sample Run of Adiabatic Piston
Cold, dense gas
Piston
Hot, dilute gas
http://cims.nyu.edu/~donev/FluctHydro/Piston.2D.m=100.rigid.hybrid.mov
 = 4/3 eq
T = ¾ Teq
 = ¾ eq
T = 4/3 Teq
Note: Adiabatic Piston is a simple Brownian heat engine
Time Relaxation of the Piston
Hybrid using deterministic PDEs underpredicts the time relaxation of the piston.
Hybrid using stochastic PDEs in excellent
agreement with pure particle calculations.
A. Donev, J.B. Bell, ALG, and B. Alder, SIAM Multiscale Mod. Sim. 8 871-911 (2010).
Relaxation from Mechanical
Non-Equilibrium
DSMC Variants for Dense Gases
DSMC variants have been developed for dense gases
of hard spheres,
* Consistent Boltzmann Algorithm (CBA)
* Enskog-DSMC
and for general potentials,
* Consistent Universal Boltzmann Algorithm (CUBA)
Basic idea is to modify the collision process so that
the collisional transport produces the desired non-ideal
equation of state.
Consistent Boltzmann Algorithm (CBA)
In CBA, a particle’s position as well as its velocity
changes upon collision.
The displacement of
position has magnitude
equal to the diameter.
Before
After
Direction is along the apse line (line between
their centers) for a hard sphere collision with
the same change in the relative velocity.
After
Before
Very easy to implement in DSMC (really!).
F. Alexander, A. Garcia and B. Alder, Physical Review Letters 74 5212 (1995)
Consistent Universal Boltzmann
Algorithm (CUBA)
Making CBA displacement a
function of density and
temperature allows you to
choose the equation of state.
Vapor condensation into droplet
Density contours at t =2,100,300,5000
Using van der Waals EOS we
can even form condensation
into a liquid with DSMC.
Not computationally efficient
versus Molecular Dynamics.
F. Alexander, A. Garcia and B. Alder, Physica A 240 196 (1997).
", N. Hadjiconstantinou, A. Garcia, and B. Alder, Physica A 281 337-47 (2000).
Stochastic Hard-Sphere Dynamics
vij
vn
vj
vi
dS
d
When two particles, i and j, are less than a diameter apart (|rij| < d)
there is a probability rate (3χ/d) vn for them to collide, where vn =
–vij . rij/|rij|; particles only collide if they are approaching (i.e., if
vn > 0). Accepted collisions are evaluated deterministically as if the
particles had a hard-sphere diameter of ds = |rij|.
A. Donev, B.J. Alder, and ALG, Physical Review Letters 101, 075902 (2008).
Properties of SHSD
Stochastic Hard Sphere Dynamics (SHSD) is
equivalent to a fluid with a linear core pair potential.
Pair Correlation
Function ( = 1)
Fluctuations of
density are
consistent with
the equation of
state (i.e.,
compressibility)
A. Donev, B.J. Alder, and ALG, J. Statistical Mechanics P11008 (2009).
References and Spam
Reprints, pre-prints and
slides available:
www.algarcia.org
DSMC tutorial &
programs in my
textbook.
Japanese
RGD 2012 in Zaragoza, Spain
Summer 2012
Hosted by ZCAM, the Spanish node of the European Centers for
Atomic and Molecular Calculations (CECAM)
DSMC 2011 Workshop
Late September 2011
Santa Fe, New Mexico
Hosted by Sandia Nat. Lab.
Von Neumann Symposium on
Multi-scale Algorithms
July 4-8, 2011
Snowmass, Utah
Sponsored by
American
Mathematical
Society