modeling by nuhman10


									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
               Department 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

                             Principal results of the authors and coworkers for
                         the modeling and simulation of turbulent mixing are
                         presented. We present simulations and theories which
                         agree with experiment in several macroscopic flow
                         variables for 3D Rayleigh-Taylor unstable flows. The
                         simulations depend on advanced numerical
                         techniques and on accurate physical modeling.
                         Sensitivity to the choice of numerical method is
                         demonstrated. Richtmyer-Meshkov unstable flows (in
                         2D) show dependence on physical modeling and on
                         the numerical method for temperature and
                         concentration probability density functions.
2                                                                                  author name et al.

Correspondence/Reprint request: Prof. James Glimm, Department of Applied Mathematics and Statistics,
Stony Brook University, Stony Brook NY 11794-3600, USA. E-mail :

     Rayleigh-Taylor (RT) and Richtmyer-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 simulations 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 interfacial area for the
surface separating the two fluids. Such 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 tracking: control of numerical diffusion
     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 years, and is now highly robust, and openly
available for downloading. 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 differential
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:LiuLiGli05], 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
     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 stimulations
     We demonstrate 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:LiuGeoBo06] . With      A   2  1  /  2  1  and g denoting gravity,
the self similar bubble side penetration distance h has the form h   Agt ,

where       b    is the dimensionless bubble side growth rate. See Table 1.
Moreover, the simulations agree with other experimentally observed
quantities: the bubble radius and the bubble height fluctuations. In the self
similar flow regime, these quantities also have a growth rate proportional to
t 2 , leading to the definition of an  r and an  h to characterize the growth of
the bubble radii and the bubble height fluctuations.
                   Table 1. Comparison of RT growth rates.
              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

    Simulation results of others are summarized in [7:alpha, 8:c&c]. Briefly,
these investigators find a growth rate  b half or less of the experimental value.
Uncertainties 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 simulations 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 factor of 3 with change of numerical method
(tracked vs. untracked), and over a factor of 50% within the untracked
4                                                                 author name et al.

   Figure 1. Bubble growth rate vs. dimensionless surface tension.
Comparison among experiment and several simulations.

    Comparison 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 tracked and untracked simulations. To
identify numerical mass diffusion as the leading cause of this discrepancy, we
determine a time dependent Atwood number A(t ) from the simulation (tracked
and untracked) itself [10: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 simulation. 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 time 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  .
    Presently 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

4. Theory to predict 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
neighboring 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 ensemble 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
                                                bubble height to diameter ratio
                                                is given by  b / 2 r  3.0 in
                                                predictions of the model and in
                                                agreement with experiments.

                                                Figure    2.     Above:     Two
                                                simulations 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.
6                                                                  author name et al.

     The two key inputs to the model are a formula for the velocity of 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 candidates 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 velocity 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,
experimentally, 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 from 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 single 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 phenomenological
aspect of this model, has been studied separately through experimental data
analysis and simulations of bubble interactions [14:GliLi88,15:GliLiMen90].

    The above conceptual definition of the bubble merger model allows the
derivation of the formula
                              1            1 1
                        b  cb r 2  
                                                   h
                              2            2k 2 
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 (the 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

      Here      cb  vterminal /  Agr   0.532 for hexagonal geometry bubbles
[16:Abarazi]. Note is taken of uncertainties in the “post terminal”
oscillations observed in the terminal velocity [17,18:Lin] and extrapolation
of the A  1 formulas for bubble velocity to values of A less than one.
Also k  0.42  0.43 is a parameter relating the new to the old bubble
radii after an area preserving 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, k  2  1  0.414 .
    A related bubble merger model [19;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
presumably as in [20:alon] by a potential theory analysis of bubble merger.
This model does not agree with experimental values for the bubble width to
bubble height ratio, a fact that has been noted [21: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 the extreme bubble, not the average bubble
velocity, in contrast to [19]. 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 models
     Buoyancy 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 Z i (t ) , i  1, 2 or i  b, s to denote the position of the
edge of the mixing zone at which fluid i vanishes, In this notation, bubbles
( b ) correspond to i  1 . In a self similar scaling regime, we expect
              (1)i Z i (t )   i Agt 2 (RT); ( 1)i Z i (t )   i Agt i (RM) .
Here sign conventions are fixed by an assumption that the light fluid is above
the heavy and g  0 . The buoyancy equations are ordinary differential
equations for Z i (t ) and the above are constraints on the choice of parameters
in these equations.
     We have favored an approach which minimizes the number of free
parameters in the buoyancy drag equations. For Vi  Z i  dVi / dt , we write
8                                                                 author name et al.

the equation derived from force balance among inertial, buoyancy and drag
forces, with the inertial term corrected for added mass effects.
                                 dVi        C  V2
                         (1)i        Agt  i i ' i
                                  dt        1  2
Here i  i  3 is the index complementary to i .

    Even this minimal equation has four undetermined parameters, the four
values of the drag coefficient C i for i  1, 2 and for RT and RM cases. Our
main result is to reduce this number to one, for example the RT drag
coefficient C1 , 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 achieve
full agreement with experiment and with the theoretically required A = 1 spike
values [27:Cheng, 28: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 asymptotes for the
buoyancy drag equation link the RT drag coefficients and the RT growth rates.
                         i 
                                 2 1  Ci (1  (1)i A) 
The first step is a postulate that the RT values of the drag coefficients C i 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             Ci (1  (1)i A)
                              1
                        i          1  A  (1  (1)i 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 C i , previously specified by the RT grow rates.
    The second main step is to relate the bubble and spike data, thereby
eliminating one of the two remaining parameters. Here we introduce another
hypothesis: that the center of mass Z COM (t ) of the mixing zone is stationary.
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 makes
the same frame assumption. This hypothesis is satisfactory for A not too close
to 1 , but it is not exactly correct, especially for A  1 . Thus, we introduce a
refined      assumption,      that    the    COM      satisfies    self    similar
scaling, Z COM (t )   COM Agt .

      A further analysis shows that     s /  b satisfies a quadratic equation whose
coefficients depend on         A and  COM /  s . Taking the limit A  1, we can
evaluate       s ( A  1)  0.5 (free fall) in the solution of this quadratic equation,
and obtain
                                                 1      b          7
                 COM ( A  1) /  s ( A  1)  
                                                 6  s ( A  1)     60
   b  0.05 . With the values at A  1 fixed, we further postulate
                              COM  A
where a fit to the LEM data [27:DS] requires   10  7

Figure 3.      s / b   as a function of   A . Comparison of experiments [25:DS]
and theory [27:cheng-a].
10                                                                     author name et al.

         We compare the experimental data for        s / b   and    as a
function of A in Figures 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  COM  0
for A  0.8 .
short title                                                                  11

5. Averaged equations and mixture models
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 proposed 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 each
phase separately, are retained in the averaged flow equations. Each 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 hyperbolically 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 model 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 conducted
[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 model [31], which is hyperbolically stable and has been
successfully applied to the modeling of ICF experimental data.
5.2 Formulation
    Let X i denote the characteristic function of phase i , so that with
averages given by angle brackets,   i  X i    is the volume average of phase

i . Density k  X k  / k and pressure pk are defined as volume
averages while velocity vk    X k vz / k , energy, internal energy and
entropy, Ek , ek , S k are defined as mass weighted quantities. Multiplying the
primitive (Euler equations of compressible fluids) equations by X i leads to
the equations:
                               k       
                                    v* k  0
                                t        z
                      (  k k ) (  k k vk )
                                                0
                          t           z
12                                                                                  author name et al.

(  k k vk ) ( k k vk vk )    ( k pk )        
                                            p * k   k k g
     t              z               z              z
 (  k k Ek ) (  k k vk Ek )    (  k pk vk )         
                                                  ( pv)* k   k k vk g
       t              z                 z                 z
Here we have introduced expressions v * , p * and ( pv)* still 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
mixing 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
                               q*  1q q2  2
                                                            q  v, p
( pv)*  p *(  v   v )  v *(  p2  2pv p1 )  ( 1pv p2v2  2pv p2v1 )
                                 2 1                 1

                         kq                        ,     q  v, p, pv
                                   k  d kq  k '
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. Conservation 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 mixing. But there is an entropy inequality for flows
which are isentropic before averaging [32:Jin].
     The identity d k d k '  1 is required for the identity  k   k '  1 , and
                     q     q                                           q       q

                                                                           p   pv
so the closure has three free parameters. Of these, two ( d k , d k ) play no
role for RT and RM data, and so we set these parameters to unity. There
remains one free parameter.
     The equations are missing one condition at each edge of the mixing
zone. A count of the number of characteristics 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 d k to the edge motions, with
short title                                                                   13

                                          Vk '
                                 d kv 
in the incompressible limit. For RM flows, d k 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 [37: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 derived 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 comparison, 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),
while for the Abrall-Saurel [41:Saurel] these same relative errors 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 * (left) 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.
14                                                                 author name et al.

Figure 7. Density plot for circular RM simulation (FronTier), shown at late
Table 2. L1 norm relative errors for closure terms as compared to simulation
                                    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 model problem consists of an
incoming 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 described 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 problem is defined in detail in
[42:Yu,43:masser] Convergence 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 averages. For most variables studied, a small interval of
angular averaging 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
remove these discrepancies. It 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 differences remain even after some level
of physical mass diffusion is allowed in the FronTier simulation.
     To illustrate the problems created by mixed cell thermodynamics, we
present in Table 3 the width, in mesh units, for the mass and temperature
discontinuities 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 diffusion layer thickness in mesh units
after reshock for a 1D radially symmetric RAGE simulation.
                    mass diffusion layer    5r
                    thermal diffusion layer 63r
    In Figure 8, we show a side by side comparison of the RAGE (left) and
FronTier (right) simulations at a time shortly after reshock, contrasting the
temperature fields at late time.

                                              Figure 8. RAGE (left) and
                                              FronTier (right). Temperature
                                             fields compared at a time
                                             following reshock for a tin-
                                             lucite simulation.

We show in Fig. 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 simulation. More comparable are the FronTier simulations of Sect. 7 with
finite Reynolds and Schmidt numbers. The RAGE PDF is roughly intermediate
between our Sc  1and Sc  0.1tracked simulations, a fact consistent with
16                                                                 author name et al.

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 x after reshock. We show in Figure 10 the ratio of the
interface length to mixing zone area in mesh units, i.e. with this divergent
factor of x scaled out. The flow is highly chaotic. On this basis, we
postulate a simple model for simulation error,
               Error  C1  x  Interface Length  C1C2
for the unregularized simulation, i.e. for the simulation with zero transport

Figure 10. The
interface, after
reshock, occupies
about 20-30% of
the available grid
blocks, uniformly
under mesh
refinement, for
short title                                                                        17

The C i 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-zero transport
coefficients. To 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
     For the separation of numerical from physical Schmidt 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 scale
K mesh  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         Remesh  1 or equivalently Kmesh  K / x  1 ,
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
18                                                                author name et al.

trends, we plot the mesh length to area ratio, at a fixed time t  90 , 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), depends 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 ( Sc  0.1 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 (slower)

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 wish 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 length to
                                      area ratio (see Fig. 11) at t  90 , i.e. at
                                      the beginning 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

                                             The first is a correlation length,
                                        for the bubble 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 point, by a
Poisson distribution. Thus the PDF of distance to exit is        1exp( x /  )   and
here    d(Probability to exit)/dx . Also   corr is the correlation

length for this exit distance.  is determined by a fit to simulation data. In Fig.
14, we show corr as a function of Re mesh , to indicate the degree of mesh
convergence achieved for corr and thus for the ratio diff / corr .

Figure 14. Plot of     corr vs. Remesh for Sc  10,1.0,0.1 .
      The diffusion distance is conveniently defined as diff           (t  t0 )
where         is coefficient of mass diffusion and t0 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 t0 , but in view of the strong increase in
interface length shortly after reshock, we take t0 to be the time of reshock. A
ratio diff / corr  1 , which is satisfied for all LES cases ( Sc  0.1 is the
smallest value we consider). For DNS cases, with Sc  10 , are not diffusive,
while the opposite case, DNS with Sc  1,0.1 , is diffusive. The diffusion
length scale     diff has converged on all meshes considered here.

7.5 Concentration PDF’s
    A principal finding is that concentration and temperature PDFs show
dependence on the details of physical values for transport coefficients. The
20                                                                author name et al.

goal of a predictive simulation, of course, is to compute independently of the
mesh, but without the prohibitive expense 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 determined from
experimental data. For the present purposes, this prescription, especially the
dependence on experimental data, is outside of the objectives of the

                                          Figure15. Light fluid concentration
                                          PDF for DNS simulations and
                                          Sc  10,1.0,0.1 .

According 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. We
hope to obtain concentration 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 physical 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 level 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.
Re mesh in Figure 16. The relative error for each principal wave is defined
as   computed - exact / exact . Here              is a time L1 norm; and exact is
interpreted to be the result of the finest grid computed with zero transport
terms. The total wave error results from averaging over the shock, mixing zone
center and mixing zone width.
Figure 16. Plot of relative wave error vs. Re mesh
      It is evident 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 achieved both in the same calculation. We can predict the level of mesh
needed to follow this simulation strategy and achieve convergence for both in
the same simulation. It appears that this can be accomplished with hardware
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 simulation, 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 converged 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 numerical mass diffusion, and so the algorithm would
not be feasible, or would require additional levels of mesh refinement.
8. Discussion, comments and conclusions
The good news is that many macroscopic flow observables for RM mixing
appear to be insensitive to details of numerical and physical modeling. The bad
news is that some important RM observables, including temperature and
22                                                                 author name et al.

concentration PDFs do appear to be sensitive to numerical and physical
modeling details. For RT mixing, it appears that many flow 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 mixed 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 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, there 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 original
     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 levels, 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 both 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 mesh 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
quantities 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 be 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 proposal, 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 fluctuations) does
not close the RT simulation issues. 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.

   This work is supported in part by the U.S. Department of Energy and the
Army Research Organization.
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. Phys. 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,
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, Institute 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 Applicatae 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.
24                                                                    author name et al.

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 10th International Workshop on the
      Physics of Compressible Turbulent Mixing, Ed M. Legrand and M.
      Vandenboomgaerde. Commisariat a l’Energie Atomique, Bruyerer-le-Chatel,
30.   Keyfitz, B. 1991, SIAM J. Math. Anal. 22, 1284.
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 Differential
      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
      Grove, 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. 2007, 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.

To top