Research Signpost
37/661 (2), Fort P.O., Trivandrum

695 023, Kerala, India
Recent Res. Devel.
Immunology, 3(2001): 145

159 ISBN: 81

7736

073

6
Modeling and Simulation of
Turbulent Mixing
Hyunkyung Lim,
Yan Yu
, James Glimm* and Xiaolin Li
Departme
nt of
Applied Mathematics and Statistics, Stony Brook
University, Stony Brook NY 11794

3600 USA
Thomas Masser,
Baolian
Cheng, John Grove
,
and
David H. Sharp
Los Alamos National Laboratory, Los Alamos NM
, 87545
, USA
Abstract
Principal results of
the authors and coworkers for
the modeling and simulation of turbulent mi
xing are
presented. We present
simulations and theories which
agree with experiment in several macroscopic flow
variables for 3D Rayleigh

Taylor unstable flows. The
simulations depen
d on advanced n
umerical
techniques and on
accurate
physical
modeling
.
Sensitivity
to the
choice of
numerical method
is
demonstrated
.
Richtmyer

M
eshkov unstable flows
(in
2D)
show dependence on
ph
ysical modeling and on
the n
umerical method
for temperature
and
concentration probability
density functions.
author name
et al.
2
Co
rrespondence/Reprint request: Prof
.
James Glimm
, Department of
Applied Mathematics and Statistics,
Stony Brook University, Stony Brook NY 11794

3600, USA.
E

mail
: glimm@ams.sunysb.edu.
Introduction
Rayleigh

Taylor (RT) and Ric
h
tmyer

Meshkov (RM) unstable flows
have
been a source of continued interest [
1:
Sh84]. Our main results include theories
and simulations for 3D RT unstable flows in agreement with experiment.
We
identify sensitivities in the sim
ulations to numerical modeling (numerical mass
diffusion and numerical surface tension) and to physical modeling (transport
and surface tension effects
and mixed phase thermal equilibrium
). The
sensitivities arise from a chaotic and nonconvergent interfaci
al area for the
surface separating the two fluids
. S
uch macro observables as the overall
growth of the
RT
mixing zone
show sensitive behavior
. This same issue is
explored more extensively in 2D RM simulations, where sensitivity to physical
modeling and to
numerical modeling show up in the temperature and species
concentration probability density functions
(PDFs)
.
As observed in [
2:
Pitsch],
these PDFs are a necessary input for the modeling of combustion.
2.
Front t
racking:
control of n
umerical
d
iffusion
Front tracking is a numerical algorithm in which a dynamically moving
discontinuity surface is represented numerically as a lower dimensional grid,
moving freely through a volume grid. Software to support this concept has
been improved over a number of yea
rs, and is now highly robust, and openly
available for downloadin
g. With the use of this algorithm
, finite difference
expressions are not computed across discontinuity interfaces. Rather, the
physically based jump conditions (the weak solution to the diffe
rential
equations) provide the connection between the regions on each side of the
interface. Not only does this algorithm eliminate numerical mass (or thermal)
diffusion, but it can be modified to allow arbitrary amounts of (physical) mass
diffusion
[
3:
Liu
LiGli05], a capability that
is useful when the physical mass
diffusion is small but nonzero.
Rigorous testing and comparison
to
alternate methods for interface
handling have
demonstrated the quality of the front tracking
method
[
4:
DuFixGli05].
From the
point of view of numerical algorithms, we use a higher order
Godunov scheme (MUSCL) and from the point of view of turbulence
simulations, the method is implicit LES, or ILES, i.e. no explicit subgrid
models have been employed.
3.
3D RT
s
t
imulations
We dem
onstrate sensitivity of the RT growth rate to transport (miscible
fluids) and surface tension (immiscible fluids) and to the numerical analogues
short title
3
of these effects. With setting of the physical modeling parameters to their
experimental values and with Front
Tracking to control the numerical issues,
we obtain agreement with experiment for the overall RT growth rate
[5:GeoGli
06,
6:LiuGeo
Bo06]
.
With
2 1 2 1
/
A
and
g
denoting gravity,
the self similar bubble side penetration distance
h
has the fo
rm
2
h Agt
,
where
b
is the dimensionless bubble side growth rate.
See Table 1.
Moreover, the simulations agree with other experimentally o
bserved
quantities:
the
bubble radius
and
the
bubble height fluctuations
.
In the self
similar flow regime, these quantities also have a growth rate proportional to
2
t
, leading to the definition of an
r
and an
h
to characterize
the growth of
the bubble
radii
and the bubble h
eight fluctuations.
Table 1. Comparison of RT
growth rates
.
Simulation results of others are summarized in
[7
:alpha,
8
:c&c
]. Briefly,
these inves
tigators find a growth rate
b
half or l
ess of the
experimental value.
U
ncertainties in the initial conditions
has been mentioned as a possible cause
of the discrepancy
. Here we duplicate the discrepancy in our own untracked
simulation
[9:Li93]
, and in this simulation we show that numerical mass
diffusion is the primary cause of the discrepancy, in the sense that if its effects
are compensated for, even in the untracked simulations, experimentally
consistent results are obtained
[10:GeoGli05]
. Numerical surface smoothing in
the untracked simulati
ons also plays a role in the discrepancy.
Sensitivity to both physical modeling and numerical effects is displayed in
Fig. 1. The simulated value of
b
varies by a factor of 25%
with change of
dimensionless
surface tension,
a fac
tor of
3
with change of numerical method
(tracked vs. untracked)
, and over a factor of 50%
within the
un
tracked
simulations
.
Fluid type
Experiment
Simulation
Theory
Miscible:
b
0.070
0.070
Immiscible:
b
0.050

0.077
0.067
0.060
Immiscible
:
r
0.01
0.01
0.01
Immiscible:
h
0.028
0.034
0.025
author name
et al.
4
Figure 1. Bubble growth rate vs. dimensionles
s surface tension.
Comparison among
experiment and several simulations.
Comparis
on of tracked to untracked simulations (which disagree with
experiments and with tracked simulations in the overall growth rate) clearly
indicate that numerical mass diffusion in untracked simulations is a major
contributor to the discrepancy between track
ed and untracked simulations.
To
identify numerical mass diffu
sion as the leading cause of this discrepancy
, we
determine a time dependent Atwood number
( )
A t
from the simulation (tracked
and untracked) itself
[1
0
:George]
.
Using
( )
A t
, we redetermine
b
and find
consistency among all simulations and between simulation and experiment.
See Figure 2.
The use of a time dependent Atwood number is a post processing step that
can be added to any si
mulation. For this purpose, we outline the main steps in
its definition. In each horizontal layer and for a fixed time step, we determine
extreme densities or densities of extreme values of the two fluid
concentrations. These densities serve to define a ti
me and layer dependent
Atwood number
(,)
A z t
. For analysis of the bubble growth rate, we average
(,)
A z t
over horizontal layers in the upper half of the bubble region, which in
most cases is about the top quarter of
the mixing layer. This average defines
( )
A t
. Because
( )
A t
is defined in terms of extreme, not average values, it
gives mixing information different from the local mixing rate
.
Prese
ntly missing from our simulations are two possibly canceling effects:
the role of unobserved noise in the initial conditions and the role of subgrid
short title
5
models, to introduce additional transport from unresolved scales in an ILES
simulation.
4. Theory to predi
ct growth rates
4.1
Bubble merger models
Bubble merger models were first introduced by Sharp and Wheeler
[11:SW]. They were given quantitative validity [12:GS] with the assumption
of a bubble interaction velocity derived from the relative positions
of
n
eighboring bubbles. These models are used to predict the RT bubble side
mixing zone growth rate as well as the bubble width, both in agreement with
experiment. As we formulate
them [13:chaos], they have no
free parameter. These models
postulate an ensembl
e of
bubbles of different heights and
(due to bubble merger) an
evolving mean radius. The
excellent agreement of the
model with experiment and
simulation is displayed in Table
1. Note that agreement
comprises three statistical
measures of the mixing rate
:
growth rates
for bubble height,
bubble radii
and bubble height
fluctuations. Specifically, the
b
ubble height to diameter
ratio
is given by
/2 3.0
b r
in
predictions of the model and in
agreement with experiments.
Figure 2. Above: Two
s
imulations with ideal physics,
tracked and untracked,
compared to one (tracked) with
physical mass diffusion. Below:
The same simulations, but
compared using a time
dependent Atwood number.
author name
et al.
6
The two key inputs to the model are a formula for the velocity o
f each
bubble and the criteria and method of bubble merger. The method of merger is
the simplest. It is by preservation of cross sectional area among the pair of
merging bubbles. The criteria for merger is a function of the velocities of the
two bubble can
didates for merger. When one of them acquires a negative
velocity, it is merged into its (generally larger) neighbor.
This brings us to the formula for the bubble velocity. It is composed of two
parts. The first is the single bubble velocity, i.e. the ve
locity of the bubble tip
in a periodic array of bubbles of the same radius. This quantity has been
extensively researched, and in dimensionless terms is known as the Froude
number. It is independent of the Atwood number, and is well characterized,
experime
ntally, theoretically and through simulation. The second component
to the bubble velocity, the bubble interaction velocity, has to do with non

periodic perturbations of the ensemble of bubbles occupying the light fluid
edge of the growing RT mixing layer.
We assume a hexagonal array for the
bubbles, in which each bubble is surrounded by six neighbor bubbles. The
influence of each of these neighboring bubbles on the velocity of the central
bubble is treated additively. For this purpose, we pick randomly fro
m the
ensemble a candidate bubble (i.e. a bubble height) to interact with a given
bubble. For
a
chosen pair of bubbles with
a
given relative height displacement,
we compute an interaction velocity (which can be positive or negative, relative
to the basic s
ingle bubble motion). It is positive if the central bubble is more
advanced (with a larger height), and otherwise negative. To compute this
interaction velocity, we regard the interacting bubble pair as a bubble

spike
combination at a larger length scale.
This bubble

spike combination has a
velocity determined by the previously mentioned single bubble theory, but
now applied at intermediate times and not at the time of terminal velocity. The
bubble interaction velocity, which is perhaps the most phenomenol
ogical
aspect of this model, has been studied separately through experimental data
analysis and simulations of bubble interactions [14:GliLi88,15:GliLiMen90].
T
he
above
conceptual
definition
of th
e bubble merger model
allows
the
derivation of
the
f
or
m
ula
1/2
1 1 1
2 2 2
b b r h
c
k
for the bubble growth rate
b
in terms of the radius
r
and height
fluctuation
h
growth rates. Through solution of the self

similarity
conditions (th
e RNG fixed point), we obtain the values for these growth
rates as in Table 1, and complete the definition of the model.
short title
7
H
ere
terminal
/0.532
b
c v Agr
fo
r hexagonal geometry bubbles
[1
6
:Abarazi]. N
ote is taken of uncertainties in the “post termin
al”
oscillations obser
ved in the terminal velocity [1
7
,1
8
:Lin] and extrapolation
of the
1
A
formulas for bubble
velocity to values of
A
less than one.
Also
0.42 0.43
k
is a parameter relating
the new to the old bubble
radii after an area pr
eserving bubble merger. The minor degree of
uncertainty for
this ratio
arises from its
weak dependence on fluctuations of
the bubble radius, a phenomena not included in the model; for uniform
bubbles,
2 1 0.414
k
.
A
related bubble merger model [
1
9
;
svartz]
is
based on the conservation of
probability for a distribution of bubbles of different sizes. The key input, the
bubble merger rate, as a function of the merging bubble radii,
is supplied
pr
esumably as in [20
:alon] by a
potential
theory analysis of bubble merger.
T
his
model
does not agree with
experimental values for
the bubble width
to
bubble height ratio
,
a fact that has been noted
[2
1
:Dimonte]
.
Both models
predict values for
b
close to experiment.
Our model has a richer structure, in
that the merger rate depends on the height separation as well as the bubble
radii. We also use the knowledge of height fluctuations to advance the mixing
zone edge by the velocity of th
e extreme
bubble
, not the average bub
ble
velocity, in contrast to [1
9
]. These differences presumably account for the
difference in
the conclusions of the models
and the superior agreement our
model achieves in comparison to experiment
.
4.2
Buoyancy drag mo
dels
Bu
oyancy drag models
predict bubble and spike side mixing zone growth
rates for both RT and RM mixing, both asymptotic and transient.
They have
been considered by a number of authors [22,23:Alon 24:DS,25:D,26:DS].
We
introduce the notation
( )
i
Z t
,
1,2
i
or
,
i b s
to denote the position
of the
edge of the mixing zone at which fluid
i
vanishes, I
n this notation, bubbles
(
b
) correspond to
1
i
. In a self similar scaling regime, we expect
2
( 1) ( ) (RT); ( 1) ( ) (RM)
i
i i
i i i i
Z t Agt Z t Agt
.
Here sign conventions are fixed by an assumption that the light fluid is above
the heavy and
0
g
.
The buoyancy equations are
ordinary
differential
equation
s
for
( )
i
Z t
and the
above are
constraint
s
on the choic
e of parameters
in these equations.
We have favored an approach which minimizes the number of free
parameters in the buoyancy drag equations. For
.
/
i
i i
V Z dV dt
, we write
author name
et al.
8
the equation derived from force balance among inertial, buoyancy and drag
forces
, with the inertial term corrected for added mass effects
.
2
'
1 2
( 1)
i
i i i i
dV C V
Agt
dt
Here
'
3
i i
is the index complementary t
o
i
.
Even this minimal equation has four
undetermined
parameters, the four
values of the drag coefficient
i
C
for
1,2
i
and for RT and RM cases. Our
main result is to reduce this num
ber to one, for example the RT drag
coefficient
1
C
, or equivalently the RT bubble growth rate coefficient
b
determined above from the bubble merger model. This goal is achieved for
most values of
A
, and a uniform theory valid for all values of
A
requires one
additional parameter. With these choices, we predict all available RT and RM
bubble and spike data
for the full range of allowed values of
A
,
and
achiev
e
full agreement with experiment and with the theoretically required
A
= 1 spike
values
[2
7
:Cheng
,
2
8
:Cheng

a
]
.
There are two steps to our reduction in the number of parameters.
At the
outset, we note that
an
elementary
analysis of the
large time asympt
otes for the
buoyancy drag equation link the RT drag coefficients and the RT growth rates.
1
2 1 (1 ( 1) )
i
i
i
C A
The first
step
is a postulate that the RT values of the d
rag coefficients
i
C
and
the growth rate coefficients
i
apply to the RM
data also. This postulate is
justified by data analysis, and appears to have excellent agreement with the
data.
An
asymptotic analysis of the RM solutions of the buoyancy drag
equation leads to the formula
(1 ( 1) )
1
1
1 (1 ( 1) )
i
i
i
i
C A
A A
This formula is a prediction
of RM mixing rates in terms of RT mixing
rates
, as it gives the value of the observed quantity
i
in terms of the
parameters
i
C
, previously specified by t
he RT grow rates
.
The second main step is to relate the bubble and spike data, the
reby
eliminating one
of the two remaining
parameter
s
. H
ere we introduce another
hypothesis: that the center of mass
( )
COM
Z t
of the mixing zone is station
ary.
short title
9
This statement is frame dependent and applies in the frame of the experimental
container. The buoyancy equation itself is also frame dependent, and ma
kes
the same frame assump
tion. T
his hypothesis is satisfactory for
A
not to
o close
to
1
, but it
is not exactly correct,
especially for
1
A
.
Thus,
we introduce a
refined assumption, that the COM satisfies self similar
scaling,
2
( )
COM COM
Z t Agt
.
A further analysis shows
that
/
s b
satisfies a quadratic equatio
n whose
coefficients depend on
A
and
COM
/
s
. Taking the limit
1
A
, we can
evaluate
( 1) 0.5
s
A
(
free fall) in t
he solution of this quadratic equation,
and obtain
COM
1 7
( 1)/( 1)
6 ( 1) 60
b
s
s
A A
A
for
0.05
b
. With the values at
1
A
fixed, we further postulate
COM
7
60
A
w
here a fit to the LEM data
[27:DS] re
quires
10 7
Figure 3.
/
s b
as a function of
A
. Comparison of experiments [25:DS]
and theory [27:cheng

a
].
author name
et al.
10
We compare the experimental data for
/
s b
and
as a
function of
A
in Figure
s
3
and 4
.
In Figure 5
, we compare our theory for
COM
to inferred data from the LES experiments [25:DS].
Figure 4.
vs.
A
; theory and experimental data [25:DS].
Figure 5
.
COM
vs.
A
; theory and experimental data. Note
0
COM
for
0.8
A
.
short title
11
5.
Averaged equa
tions and mixture
m
odels
5.1
Introduction and discussion of mix models
Averaged equations suppress many solution details and allow simplified
solutions to problems through averaging over the details of the fine scale
structure. A discussion of many propos
ed models and closures, with a listing of
the number of adjustable parameters and other key features
was given in
[29
:Cheng]. We propose mix models based on a
complete first order
multiphase closure. This means that all primitive variables, averaged over e
ach
ph
a
se separately, are retained
in the averaged flow equations
.
E
ach phase has
an independent pressure and complete thermodynamical variables. Mixed
phase thermodynamics is not part of this model. The model has the pleasant
feature of being hyperbolical
ly stable, meaning that independently of the
viscosity assumed in the equation, the time propagation has real characteristics
and is stable.
In contrast, a closure with fewer independent variables is called a reduced
model. A commonly studied reduced mod
el with equal phase pressures has
complex characteristics and is stable only with assumption of a sufficiently
large viscosity. The instability has been known for many years. Detailed
studies, mathematical and numerical, of the instability have been conduc
ted
[30
:Keyfitz], and show phase separation (of mixed phase
flow regimes
into
separate
unmixed
regions of pure phase) as the instability mode. However, not
all single pressure reduced models are hyperbolically unstable. We mention the
Scannapieco

Cheng mod
el [31
],
which is hyperbolically stable and
has been
successfully applied
to the modelin
g of ICF experimental data
.
5.2
Formulation
Let
i
X
denote the characteristic function of phase
i
, so that with
averages
given by angle brackets,
i i
X
is the volume average of phase
i
. Density
/
k k k
X
and pressure
k
p
are defined as volume
averages while
velocity
/
k k z k
v X v
, energy, inte
r
nal energy and
entropy,
,,
k k k
E e S
are defined as mass weighted quantities. Multiplying the
primitive (Euler equations of compressible fluids) equations by
i
X
leads to
the equations:
* 0
k k
v
t z
( ) ( )
0
k k k k k
v
t z
author name
et al.
12
( ) ( ) ( )
*
k k k k k k k k k k
k k
v v v p
p g
t z z z
( ) ( ) ( )
( )*
k k k k k k k k k k k
k k k
E v E p v
pv v g
t z z z
Here
we have
introduced expressions
*
v
,
*
p
and
( )*
pv
s
till to be
defined, and we have omitted the within

phase Reynolds and related heat
flux tensors, as these appear to be small in the data we study. Extensions to
higher dimensions (assuming a preferred normal direction (
z
) to the
m
ixing layer) have been studied.
Mathematical analysis
[31,32:Jin]
of the unclosed equations suggests
a functional form for the closure terms as convex combinations of the
corresponding
pure phase quantities
1 2 2
* ,
q q
q q q v p
1 2 2 1 1 2 2 1 1 2 2 2 2 1
( )* *( ) *( ) ( )
pv pv pv pv pv pv
pv p v v v p p p v p v
where
'
, ,,
q
k
k
q
k k k
q v p pv
d
These equations satisfy boundary conditions at the edge of the mixing
zone, and are hyperbolic. Modifications for inclusion of transport in the
primitive (unaveraged) equations have been derived.
Conserv
ation of mass,
momentum and total energy follows from the form of the equations.
Entropy is not conserved, even for flows which are isentropic before
averaging. The reason for this is that averaging itself is not isentropic, as it
modifies the entropy of m
ixing. But there is an entropy inequality for flows
which are isentropic before averaging [32:Jin].
The identity
'
1
q q
k k
d d
is required for the identity
'
1
q q
k k
, and
so the closure has three free parameters. Of thes
e, two
(,)
p pv
k k
d d
play no
role for RT and RM data, and so
we set these parameters to unity. T
he
re
remains
one free parameter.
The equations are missing one condition at each edge of the mixing
zone. A count of the number of character
istics shows that at those
locations, there is a missing characteristic, associated with the incoming
sound wave of the fluid missing in the single phase region.
Thus we couple
the model to the buoyancy drag equations for the mixing zone edges, to
give an
extra condition there. A mathematical analysis
for the RT case
shows the connection of
v
k
d
to the edge motions, with
short title
13
'
v
k
k
k
V
d
V
in the incompressible limit.
For RM flows,
v
k
d
is not important
in the
model and is conveniently set to one.
The equations, in
the 1D
RT
incompressible case, admit a closed form
solution
[
34,35,36:saltz] The closed form solution was extended to weakly
compressible flows in a singular perturbation expansion [3
7:Jin
]
.
The
extension of these ideas to multiple phases was given in [38:Cheng].
5.3
Numerical studies of mix models
We solve
numerically the
equations
for
the
complete first order closure
[39
:Jin]. In this study, we compare the solutions to the theory de
rived above,
and we solve a sample higher dimensional flow problem.
Using the RT and RM two fluid simulation data sets discussed here, we
compare
[40
:Bo]
closure model residual errors between our closure and that of
Abgrall and Saurel. In this compar
ison,
our errors are about two
times smaller
than theirs.
The average relative error for the three closure terms for our model
is 12% (closure residuals compared to RT simulation data) and 9% (RM),
w
hile for the Abrall

Saurel [41
:Saurel] these same relative er
rors are 30%
(RT)
and 13%
(RM)
.
See Table 2.
The Abrall

Saurel closure model errors are
computed assuming no relaxation terms. If their relaxation terms are included
their errors are larger.
See Figure 6.
Figure 6. Plot of
*
v
(l
eft) and
( )*
pv
(right)
times
/
d dz
as an exact
expression computed from the 3D RT FronTier simulation data and closed
expressions to define
*
v
and
( )*
pv
approximately.
Our own closure is
compared to two interpretations of the Abgrall

Saurel closure, with and
without their relaxation terms.
author name
et al.
14
Figure 7. Density plot for circular RM simulation
(FronTier)
, shown at late
time.
Table 2.
1
L
norm relative
errors for closure terms as compared to simulation
data.
v*
p*
(pv)*
Average
RT data comparison
18%
0%
18%
12%
RM data comparison
7%
0%
20%
9%
6. RM code comparisons
Code comparison of RAGE to FronTier
was conducted for a circular
Richtmyer

Meshkov
problem in 2D. This mo
del problem consists of an
incom
ing circular shock wave crossing a
perturbed circular interface, reflecting
from the origin and reshocking the now strongly perturbed interface and
mixing zone.
The inner
and outer materials are describ
ed
by a stiffened gamma
law gas with parameters chosen to
model
either lucite or air as the
interior
fluid
and tin
as the
exterior
fluid
.
For the purpose of this comparison, all transport
coefficients have been set to zero.
The
probl
em is defi
ned in detail
in
[42
:Yu
,43
:masser] C
onvergence under mesh refinement for the FronTier code
was analyzed
in [42
:Yu]
. The solution develops new structures at each new
length scale under mesh refinement, and so convergence was considered in the
sense of local
spatial aver
ages. For most variables studied
, a small interval of
angular ave
raging was sufficient to allow
convergence.
Density contours
showing the highly chaotic mixing zone at
late time are presented
in Figure
7
.
The comparison study shows approximate agreement
for
such macro
solution features as
the locations of the lead shock waves and mixing zone
edges
. But
micro solution features such as the
temperature and concentration
PDFs are not comparable,
and systematic mesh refinement does not appear to
re
move these d
iscrepancies. I
t appears that ILES simulations (e.g those with
short title
15
zero or under resolved regularizing transport
and no explicit subgrid model
)
are sensitive to numerical transport coefficients and models for mixed cell
thermodynamics [43
:Masser07]. The differ
ences remain even after some level
of physical mass diffusion is allowed in the FronTier simulation.
To illustrate the problems created by mixed cell therm
odynamics, we
present in Table 3
the width, in mesh units, for the mass and temperature
discontinu
ities after reshock, for an unperturbed (radially
symmetric
1D)
simulation of the same physics.
Use of an EOS to model an air

tin rather than
a lucite

tin problem leads to significantly larger diffusion widths [43].
Table 3
. Numerical mass and thermal di
ffusion layer thickness in mesh units
after reshock
for a 1D radially symmetric RAGE simulation
.
Lucite

Tin
mass diffusion layer
5
r
thermal diffusion layer
63
r
In Figure 8
, we show a side by side compa
rison of the RAGE (left) and
FronTier (right) simulations a
t a time shortly after reshock, contrasting
the
temperature fields at late time.
Figure 8. RAGE (left) and
FronTier (r
ight). Temperature
fields compared at a time
following reshock for a tin

luc
ite simulation.
W
e show
in F
ig. 9
the RAGE light fluid concentration PDF for the ideal
tin

lucite
simulation. The comparable FronTier simulation has a PDF delta
function showing zero fluid mixture, in accordance to the physical modeling of
the sim
ulation.
More comparable are the FronTier simulations of Sect. 7 with
finite Reynolds and Schmidt
numbers. The RAGE PDF
is roughly intermediate
between our
1
Sc
and
0.1
Sc
tracked simulations, a fact consistent with
author name
et al.
16
the mass diffusion width of the RAGE calculation, suggesting a
numerical
Schmidt number of perhaps 0.3.
Figure 9
. RAGE PDF for
light fluid mass
concentration, taken along a
mid line of the mixing zone.
7. RM sensitivity to physical modeling
7.1
Introduction
Under mesh refinememt, the interface (50% concentration isosurface)
between the two fluids fails to converge, and its length is approximately
proportional to
1
x
after reshock. We show
in Figure 1
0
the
ratio of the
inter
face length to mixin
g zone area
in mesh units, i.e. with this divergent
factor of
1
x
scaled out
. The flow is
highly chaotic. On this basis, we
postulate a simple model for simulation error,
1 1 2
Error x Interface Length
C CC
for the
unre
gularized simulation, i.e. for the simulation with zero transport
coefficients
.
Figure 1
0
. The
interface, after
reshock, occupies
about 20

30% of
the available grid
blocks, uniformly
under mesh
refinement, for
th
unregularized
simulation.
short title
17
The
i
C
are order 1, apparently mesh independent quantities.
This
formula
raises questions for the meaning of the simulation,
especially for interface sensitive quantities, such as the degree of
atomic
scale
fluid mixing.
It also raises the
possibility that for quantities depending
on physical transport, such as Reynolds and Schmidt numbers, there is a
possibility of apparent convergence to a nonunique solution selected in a
code dependent manner by numerical Reynolds and Schmidt numbers.
7
.2
The simulation study
To investigate grid convergence and sensitivity to physical transport
properties
, we studied the same
2D RM
problem with non

ze
ro transport
coefficients. T
o determine the dependence on these coefficients, we
allowed them to
vary over wide ranges, as needed to illustrate a range of
LES and DNS simulations
, and to separate numerical from physical
transport
[44
:Lim]
.
See Figure 11
for
parameters describing the
simulation
s.
For the separation of numerical from physical Sc
hmidt number, the
FronTier code is especially convenient, as it can simulate any Schmidt
number through the controlled interface diffusion algorithm, without the
expense of additional mesh refinement.
Figure 11
.
The
Reynolds number and Kolmogorov
s
ca
le
mesh
/
K K
x
vs. mesh Reynolds number.
7.3
The DNS

LES transition
We divide the simulations into those in an LES regime and those in a DNS
regime. The division between the two is not sharp, but is characterized roughly
by the condition
mesh
Re 1
or equivalently
mesh
/1
K K
x
,
where
K
is the Kolmogorov length scale. We find signs of interface
convergence in the DNS regime but not in the LES regime. To show these
author name
et al.
18
trends, we plot the mesh l
ength to area ratio, at a fixed time
90
t
, after
reshock and at the beginning of
the chaotic regime. In Figure 12
, we see that
convergence of the interface (as measured by a decrease in the mesh interface
length to area ratio), depen
ds primarily on the mesh Reynolds number or
equivalently on the mesh Kolmogorov length, and is otherwise largely
independent of mesh and Schmidt number.
To understand numerical convergence for the DNS regime, we replot the
same data (
0.1
Sc
only) in Figure 13
, interface fraction (length/area) in
physical units vs.
Re
. By this measure,
the interface has converged
in the
DNS regime. Higher
Sc
numbers give similar behavior,
but
delayed (sl
ower)
convergence.
7.4
Physical and numerical mass diffusion
Starting with a mesh Reynolds number of about 1, it appears that the interface
stops growing as the mesh is further refined. Imagining that we are dealing
with a mesh converged
interface, we wis
h to understand the degree of mass
diffusion to expect in the simulation.
Intuitively, this is measured by the ratio
of two length scales.
Figure 12
. The mesh interface l
ength to
area ratio (see Fig. 11
) at
90
t
, i.e. at
the begin
ning of the chaotic regime, for
simulations conducted with a variety of
grid levels, mesh Reynolds numbers and
Schmidt numbers.
short title
19
Figure 13. Replot of data from
Figure 11
, using physical, not mesh
units.
The first is a correlation length,
for the bubb
le size distribution, and
the second is a diffusion length, for
the diffusion distance achieved in
the time which has elapsed.
To define a correlation length,
we model the
probability
a
of
minimum distance to exit
a given phase, given a random starting po
int, by a
Poisson distribution.
Thus the PDF of distance to exit is
1
exp(/)
x
and
here
1
d(Probability to exit)/dx
. Also
corr
is the correlation
length for this exit distance.
is determin
ed by a fit to simulation data.
In Fig.
14
, we show
corr
as a function of
mesh
Re
, to indicate the degree of mesh
convergence achieved for
corr
and thus for the ratio
diff corr
/
.
Figure 14
. Plot of
corr
vs.
mesh
Re
for
10,1.0,0.1
Sc
.
The diffusion distance is conveniently defined as
diff 0
( )
t t
where
is coefficient of mass diffu
sion and
0
t
is the time for creation of an
element of interface length. Since the interface length grows continuously with
time, there is no precise way to set
0
t
, but in view of the strong increase in
inter
face leng
th shortly after reshock, we take
0
t
to be the time of reshock.
A
ratio
diff corr
/1
, which is satisfied for all LES cases (
0.1
Sc
is the
sm
allest value we consider). F
or DNS cases
,
with
10
Sc
, are
not diffusive,
wh
ile the opposite case, DNS
with
1,0.1
Sc
,
is
diffusive. The diffusion
length scale
diff
has converged on all meshes considered here.
7.
5
Concentration PDF’s
A
pr
incipal finding is that c
oncentration and temperature PDFs show
dependence on the details of physical values for transport coefficients.
The
author name
et al.
20
goal of a predictive simulation, of course, is to compute
independently of the
mesh, but without the prohibitive e
xpense of DNS. According to [2:Pitsch],
subgrid models will compensate for the missing scales, as far as viscosity is
concerned, while the atomic level mass diffusion is fit to a convenient
distribution (the beta

distribution) with coefficients determine
d from
experimental data. For the present purposes, this prescription, especially the
dependence o
n experimental data, is outside
of the objectives of the
computation.
F
igure15
. Light fluid concentration
PDF for DNS simulations and
10,1.0,0.1
Sc
.
A
ccording to [2], the subgrid model is defined in terms of
filtered variables
with the filter length
larger than the grid spacing. We approach this latter
objective indirectly, replacing the filter and its length scale larger than the
mesh
with large (artificial)
coefficient in a physical
viscosity
term
, to turn an
LES into a (numerical) DNS. According to Figs.
12
and 13
, this will cut off
the
formation of small interface structures under continued mesh refinement
.
W
e
hope to obtain
concent
ration PDFs by using the correct Schmidt number for
the simulation.
The general shape of the concentration PDFs appear to be
short title
21
largely determined by the ratio
diff corr
/
. To force this ratio to its converged
value, we want to preserve the p
hysical Schmidt number for the simulation,
rather than the physical value mass diffusion directly.
7.6
Convergence of macro solution features
The convergence of principal solution waves (shock and mixing zone
edges) is adversely affected by the high leve
l of viscous damping in the DNS
simulations performed at moderate grid resolution. To understand this
degradation of solution accuracy, we plot the relative position error vs.
mesh
Re
in Figure 16
. The relative error for each principal
wave is defined
as
computed  exact/exact
. Here
is a time
1
L
norm
; and exact is
interpreted to be the result of the finest grid computed with zero transport
terms. The total wave error results from a
veraging over the s
hock, mixing zone
center and mix
ing zone width.
Figure 16
. Plot of relative wave error vs.
mesh
Re
It is evide
nt from
comparison of Figures 13 and 16
that we can achieve
convergence with either the micro or the macro
observables, but
we have
not
yet ac
hieve
d
both in the same
calculation. We can predict the level of mesh
needed to follow this simulation strategy and achieve converge
nce for both in
the same simula
tion
.
I
t appears that this can be accomplished with hardw
are
currently existing
. But we observe that this plan involves subgrid modification
of the momentum equation only with a numerically chosen viscosity, while the
mass diffusivity is chosen to give a physical Schmidt number.
For a more practical LES simula
tion, we regard both the viscosity and
mass diffusion coefficients as defining a subgrid term, for the momentum and
concentration equations respectively. For this strategy, we determine the ratio
diff corr
/
from a simulation that has conve
rged in its
micro variables, as in
Figure 13
,
and force the ratio to have the same value in a macro

converged
LES simulation by adjustment (as a mass diffusion subgrid model) of the
coefficient of mass diffusivity. We note that this proposed LES algorithm
depends on a unique capability of front tracking: to simulate arbitrary (even
very small) levels of mass diffusion. For many codes, the combined physical
and subgrid levels of mass diffusion (especially for large Schmidt numer) are
less than the inherent n
umerical mass diffusion, and so the algorithm would
not be feasible
, or would require additional levels of mesh refinement
.
8.
Discussion, comments and c
onclusions
The good news is that many macroscopic flow observables for RM mixing
appear to be insensiti
ve to details of numerical and physical modeling. The bad
news is that some important RM observables, including temperature and
author name
et al.
22
concentration PDFs do appear to be sensitive to numerical and physical
modeling details. For RT mixing, it appears that many flo
w observables,
including the much studied overall growth rate, are sensitive to numerical and
physical modeling details. Sensitivity is generally traced to numerical mass
diffusion, exacerbated in the RM case by a temperature equilibrium hypothesis
for mix
ed states. Both numerical mass diffusion and mixed state thermal
equilibrium are common features of untracked simulations. Because these
mixing flows are chaotic, and have a divergent interfacial area, apparently
proportional to
1
x
in the absence of physical regularization, errors
associated with the interface, of necessity unavoidable especially for untracked
simulations, run the risk of nonconvergence. Where these simulations are
regularized by grids and numerical code features, t
here is a risk of an apparent
convergence to code dependent values.
In the language of turbulent modeling,
this feature of LES has been systematized through the use of filtered equations,
and LES that converge to the solution of the filtered, but not the o
riginal
equations.
Solutions are assessed for convergence in both macro
and micro variables.
While we have not achieved acceptable convergence for both types of
variables simultaneously, we have developed a methodology which allows
prediction of mesh leve
ls, numerical algorithms and physical or subgrid
modeling procedures which in combination will achieve this goal.
With
subgrid modification of both the momentum and concentration equations, these
algorithms are expected to yield convergence results for bot
h macro and micro
variables with feasible grids.
To obtain practical LES simulations for quantities sensitive to numerical
and physical modeling, the classical ideas of turbulence modeling propose
filtered equations with a filter length larger than the me
sh length. These
equations allow convergence under mesh refinement, but to a solution of the
filtered, not the original equations.
In the spirit of LES, convergence of the
filtered equations to the true ones is achieved (for statistically observable
quanti
ties only) by decreasing the filter length (and the mesh in proportion
)
.
We propose here to deviate from this program only in the use of an artificial
value for physical viscosity in place of the filter
. The purpose of this step is to
control the growth of
interface length under (further) mesh refinement. An
alternate route would be the use of surface tension, with a numerically chosen
(artificial) non

physical coefficient. With the interface thus stabilized,
we then
propose that
the mass diffusion should b
e set
to preserve physical values for
corr diff
/
. It appears that the
use of physical values for the Schmidt number
will accomplish this goal. In any case, it is necessary to replace the step of
calibration from experimental data. This pro
posal, of course, requires
exploration beyond the information presented here.
short title
23
The agreement of front tracking RT simulations with experiment for
several growth rate measures (bubble height, width, height fl
uctuations) does
not close the RT simulation
issu
es. There remain initial condition modeling
uncertainties and possible effects of missing subgrid models. Since these
effects will move the growth rate in opposite directions, the possibility of
canceling errors is not yet excluded.
Acknowledgment
This
work is
supported
in part by
t
he U.S. Department of Energy and the
Army Research Organization
.
References
1.
Sharp, D.H. 1984, Physica D 12, 3.
2.
Pitsch, H.
2006,
Ann. Rev.
Fluid Mech. 38, 453.
3.
Liu, X. F., Li, Y. H., Glimm, J. and Li, X. L. 2007, J. Comp. Phy
s. In Press.
4.
Du, J., Fix, B., Glimm, J., Jia, X., Li, X., Li, Y. and Wu, L. 2006, J. Comp. Phys
.
213, 613.
5.
George, E., Glimm, J, Li,
X. L., Li, Y. H. and Liu, X. F.
2006, Phys. Rev. E. 73,
016304.
6.
Liu, X. F., George, E. Bo, W. and Glimm,
J.
2006,
Phys.
Rev. E 73, 056301.
7.
Dimonte, G,, Youngs, D. L., Dimits, A., Weber, S., Marinak, M., Wunsch, S.,
Garsi, C., Robinson, A., Nadrews, M., Ramaprabhu, P., Calder, A., Fryxell, B.,
Bielle, J., Dursi, L., MacNiece, P., Olson, K., Ricker, P., Rosner, R., Timmes, F
.,
Tubo, H., Young, Y.

N. and Zingale, M. 2004, Phys. Fluids 16, 1668.
8.
Cabot, W. and Cook, A.
2006, Nature Physics 2, 562.
9.
Li, X. L. 1993, Phys. Fluids 5, 1904.
10.
George, E. and Glimm, J. 2005, Phys. Fluids 17, 054101.
11.
Sharp, D.H. and Wheeler, J 1961, Instit
ute of Defense Analysis technical report,
AIP Document No. PAPSFLDA

31447

50.
12.
Glimm, J. and Sharp, D.H. 1990, Phys. Rev. Lett. 64, 2137.
13.
Cheng,
B., Glimm, J. and Sharp, D. H.
2002, Chaos, 12, 267.
14.
Glimm, J and Li, X.L.
1998, Phys. Fluids 31, 2077.
15.
Glimm,
J., Li, X.L., Menikoff, R., Sharp, D. H.
and Zhang, Q.
1990, Phys. Fluids
A, 2
, 2046.
16.
Abarzhi, S.I.
1998, Phys. Rev.Lett. 81, 337.
17.
Lin, A. 2000, State University of New York at Stony Brook Thesis.
18.
Glimm, J. Li, X.L., Lin, A. D. 2002, Acta Mathematicae Appl
icatae Sinica 18, 1.
19.
Oron, D., Sadot, O., Srebo, Y., Rikanati, A., Yedvab, Y., Alon, U. Erez, L. Erez,
G. Ben

Dor, G., Levin, L.A., Ofer, D. and Shvarts, D.
1999, Laser Part. Beams
17, 465.
20.
Alon, U, Hecht, J. Ofer, D and Shvarts, D. 1995, Phys. Rev. Lett.
74, 535.
.
21.
Dimonte, G. 2000, Phys. Plasmas, 7, 2255.
22.
Alon, U., Hecht, J., Mukamel, D. and D. Shvarts 1994, Phys. Rev. Lett. 72, 2867.
23.
Alon, U., Hecht, J., Ofer, D. and D. Shvarts 1995, Phys. Rev. Lett. 74, 534.
24.
Dimonte, G. and Schneider, M. 1996, Phys. Rev.
E 54, 3740.
25.
Dimonte, G. 1999, Phys. Plasmas, 6, 2009.
author name
et al.
24
26.
Dimonte, G. and Schneider, M. 2000, Phys. Fluids 12, 304.
27.
Cheng,
B., Glimm, J. and Sharp, D. H.
2000, Phys. Lett. A 268, 366.
28.
Cheng, B., Glimm, J. and Sharp,
D. H.
2002, Phys. Rev. E, 66, 036312.
29.
Cheng
, B. 2007, In: Proceedings of the 10
th
International Workshop on the
Physics of Compressible Turbulent Mixing, Ed M. Legrand and M.
Vandenboomgaerde. Commisariat a l’Energie Atomique, Bruyerer

le

Chatel,
France.
30.
Keyfitz, B. 1991, SIAM J. Math. Anal. 22, 12
84.
31.
Scannapieco, A. and Cheng, B. 2002, Phys Lett. A 299, 49.
32.
Jin, H. Glimm, J. and Sharp, H. 2006, Phys. Lett. A 353, 469.
33.
Jin, H. Glimm, J. and Sharp, H. 2006, Phys. Lett. A 360, 114.
34.
Glimm, J., Saltz, D. and Sharp, D. H. 1998, In: Nonlinear Partial Diff
erential
Equations, Ed. Chen, G. Q., Li, Y. and Zhu, X. World Scientific, Singapore, 124.
35.
Glimm, J.
,
Saltz, D. and Sharp, D. H.
1998, Phys. Rev. Lett. 80, 712,
36.
Glimm, J., Saltz, D. and Sharp, D. H. 1999, J. Fluid Mech. 378, 119.
37.
Glimm, J. and Jin, H. 2001,
Bol. Soc. Bras. Mat. 32, 213.
38.
Cheng, B., Glimm, J. and Sharp, D. H. 2005, Phys. Fluids 17, 087102.
39.
Glimm, J. Jin, H. Laforest, M., Tangerman, F. and Zhang, Y. 2003, Multiscale
Model. Simul. 1 458.
40.
Bo, W., Jin, H., Kim, D., Liu, X., Lee, H., Pestieau, N.,
Yu, Y. Glimm, J. and
G
rove, J. 2007, Stony Brook University Technical Report
.
41.
Saurel, A. and Abgrall, R. 1999, J. Comput. Phys. 150, 425.
42.
Yu, Y, Zhao, M., Lee, T., Pestieau, N., Bo, W., Glimm, J. and Grove, J.W. 2006,
J. Comp. Phys. 217, 200.
43.
Masser, T. 2
007, Stony Brook University Ph.D. Thesis.
44.
Lim, H., Yu, Y., Jin, H., Kim, D., Lee, H. Glimm, J
. and D. H. Sharp 2007, Stony
Brook University Technical Report
.
Comments 0
Log in to post a comment