The Mathematics and Numerical Principles for Turbulent Mixing James Glimm1,3 With thanks to: Wurigen Bo1, Baolian Cheng2, Jian Du7, Bryan Fix1, Erwin George4, John Grove2, Xicheng Jia1, Hyeonseong Jin5, T.

Download Report

Transcript The Mathematics and Numerical Principles for Turbulent Mixing James Glimm1,3 With thanks to: Wurigen Bo1, Baolian Cheng2, Jian Du7, Bryan Fix1, Erwin George4, John Grove2, Xicheng Jia1, Hyeonseong Jin5, T.

The Mathematics and Numerical
Principles for Turbulent Mixing
James Glimm1,3
With thanks to:
Wurigen
Bo1,
Baolian
Cheng2,
Jian
Du7,
Bryan
Fix1,
Erwin George4, John Grove2, Xicheng Jia1, Hyeonseong Jin5, T. Kaman1, Dongyung
Kim1, Hyun-Kyung Lim1, Xaolin Li1, Yuanhua Li1, Xinfeng Liu6, Xingtao Liu1, Thomas
Masser1,2, Roman Samulyak3, David H. Sharp2, Justin Iwerks1,Yan Yu1
1.
2.
3.
4.
5.
6.
7.
SUNY at Stony Brook
Los Alamos National Laboratory
Brookhaven National Laboratory
Warwick University, UK
Cheju University, Korea
University of Southern California
University of South Carolina
1
Hyperbolic Conservation Laws
U t   F (U )  0
Status of existence theory:
1D and small data:
Perturbation expansion in wave interactions
existence, uniqueness,
smooth dependence on data
(G., Liu, Brezan, others)
1D and 2x2 system:
Existence (Diperna, Ding, Chen, others)
Measure valued solutions, then shown to be
classical weak solutions
2D: special solutions only
2
Nature of Solutions (1D)




Solutions are typically discontinuous
Shock waves form
Solution space is functions of bounded
variation and L_infty
Solutions are interpreted in the weak sense,
as distributions
3
Compensated Compactness

Proofs (for 1D, large data, 2x2 systems)


First establish existence of solutions in a very weak space of measure
valued distributions
Then (more difficult) show that such a weak solution is a classical
weak solution, as a distribution:

U



F
(
U
)

0
all
test
functions



t


Theorem (Gangbo and Westdickenberg):



For 2D, 3D: only first step (measured valued solution)
Isentropic equations Comm. PDEs (In press)
Limit may not proved to satisfy the original equation
4
Obstacles to extensions

Some 3x3 and larger systems are unstable in BV
norm



(but perhaps these are not physically motivated)
In 2D, even for gas dynamics, special solutions can
be unbounded. Also BV is unstable.
Nonuniqueness, 2D:




Scheffer (1989), Snirelman
DeLellis and Szelkelyhidi (Archives)
Volker Elling (numerical)
Lopez, Lopez, Lowengrub and Zheng (numerical)
5
Numerical implications: Standard view




Typically numerical solutions appear to be
convergent
Problems with existence seems to be a
strictly mathematical concern
This point of view is incorrect as
mathematics and as physics
This lecture: to identify and cure problems as
mathematics and as physics
6
Implications, continued

For many problems, physical regularization (viscosity, mass diffusion,
etc.) is very small.


Even if included in the simulation, the effects are under resolved. Thus
effects of physical regularization are dominated by numerical effects such as
numerical mass diffusion.
Large effort in V&V = Verification and Validation

Verification is proof that the numerical solutions are bona fida solutions of
the mathematical equations.



This step fails if mathematical solutions are nonunique or discontinuous in
dependence on initial conditions (loss of well posedness) or fail to exist or if the
nonuniqueness etc. is resolved numerically by artifacts of algorithm
Validation is agreement with physical experiments
In practice, validation fails dramatically for turbulent mixing (Rayleigh-Taylor
instability).
7
Chaotic Mixing:
A challenge to the standard view


Solutions are unstable on all length scales
Under mesh refinement, new structures
emerge



In this sense there is no convergence
Optimistically, we hope that the large scale
structures converge and the new small scale
ones that emerge under mesh refinement will not
influence the large scale ones
This hope is partially correct and partially not
8
Failure of the standard view

Turbulent mixing and turbulent combustion



Atomic or molecular level mixing requires a new
length scale (the atoms or molecules)
And a change in the laws of physics at these
length scales or above.
New terms in the conservation laws, to express
mass diffusion, viscosity, heat conduction:
U t   F (U )   D  U
9
Parabolic Conservation Laws
U t   F (U )   D  U
Even in 1D, and in the limit
D
distinct solutions occur (Smoller, Conley)
Ratios of diffusion terms at the root of issue



 0
For combustion, dimensionless ratios in D influence mixing properties
and combustion rates, thus the global dynamics (Schmidt and Prandtl
numbers). In this sense the small scales can easily affect the large ones.
For turbulent mixing, parabolic transport terms can affect solution,
especially the ratios of transport terms

Density differences drive instabilities, and so mass diffusion diminishes the
instability. Viscosity limits the complexity of the flow and hence the amount
of the interfacial mass diffusion. In the limit D  0, the Schmidt number
(viscosity/mass diffusion) determines the net amount of mass diffusion and
hence the RT instability growth rate.
10
Requirements for simulation of
turbulent combustion



Averaged flow velocities and pressures as a function
of x,y,z,t
Joint probability distributions of the concentrations
of species (fuel, oxidizer) and temperature, as a
function of x,y,z,t
This depends on the diffusion matrix D and if the
calculation is under resolved, it depends on the
numerical code and the numerical D, not the
physical one
11
Mathematical implications for
hyperbolic conservation laws




The ultra weak solutions as measure valued
distributions (PDFs) depending on x,y,z,t is exactly
the framework of compensated compactness
Existence proofs in this framework should be easier
than for classical (weak) solutions
Uniqueness, which typically fails for such weak
solution methods, should fail, and is not correct (in
the hyperbolic setting), neither as mathematics nor
as physics.
On numerical grounds, appears to correctly express
the physics of unregularized solutions.
12
Numerical example: Richtmyer-Meshkov
(shock accelerated) turbulent mixing





Circular shock starts at outer edge of (1/2)
circle.
Moves inward, crosses a perturbed circular
density difference
Reaches origin, reflects, recrosses the density
difference
Highly unstable flow
Illustrates all previous points
13
Numerical example:
Primary breakup of jet (diesel fuel);
Jet in cross flow (jet engine)





Liquid leaves narrow nozzle at speed near Mach 1. Breaks
up rapidly into small droplets.
Breakup is essential feature of flow to predict combustion
rates. Burning is mainly film burning, and so is rate limited
by surface area of droplets.
For scram jet, fuel leaves combustion chamber in 1 millisec,
and must complete combustion before that.
Flows are only stabilized by transport (viscosity, mass
diffusion or surface tension and heat conduction.
Illustrates main points for chaotic mixing.
14
Numerical example:
Chemical fuel separation
in a mixer/contactor






High shear rate in rotating device (Couette flow) produces
high levels of Kelvin-Helmholtz instability
Interface between immiscible fluids (oil and water based)
produces fine scale droplets.
Chemical reactions occur on the interface and are thus
interface area limited.
Prediction of total droplet area and distribution of droplet
sizes and shapes is needed to predict the chemical reaction
rates.
Flow is stabilized by surface tension (Weber number).
Illustrates main points for chaotic flow.
15
Numerical Example:
Rayleigh-Taylor Turbulent Mixing


Steady acceleration of a density difference
Unstable if acceleration points from heavy to
light


Water over air
Penetration distance h of light fluid into heavy
(bubbles) is a simple measure of mixing rate.
  h / Agt
2
16
2D Chaotic Solutions: Shock implosion of
perturbed interface—4 grid levels
17
Circular
RM instability
Initial (left)
and after
reshock
(right) density
plots. Upper and
lower inserts
show enlarged
details of flow.
18
Definitions for Turbulent Mixing



An observable of the flow is SENSITIVE if it shows
dependence on numerical algorithms and/or on
physical modeling (transport).
MACRO observables: edges of mixing zones,
trajectories of leading shock waves, mean
densities, velocities of each fluid species
MICRO observables: probability density functions
for temperature and for species concentration;
properties dependent on atomic mixing
19
Explanation


Macro observables describe the average
mixing properties of the large scale flow
Micro observables describe the atomic level
mixing properties
20
Sensitive and Insensitive Observables
Rayleigh-Taylor
Richtmyer-Meshkov
Macro
Micro
Sensitive
Insensitive
Sensitive
Sensitive
RM has a highly unstable interface, in the
absence of regularization, diverging as 1/delta x.
This divergent interface length or area causes sensitivity
to atomic scale observables.
21
Simplistic Error Analysis for Micro
Observables in Chaotic Flow
E R R O R  C 1   x  Interface length (area)
Interface length (area) = C 2   x
E R R O R  C1   x  C 2   x
1
1
 C 1C 2  O  1 
Error as results from numerical or physical modeling, e.g.
numerical mass diffusion or ideal vs. physical transport coefficients
22
Chaotic dependence of interface length on
mesh: Unregularized simulation
S urface Length/A rea
vs. tim e
Surface Length
Area
 x
23
Verification (Macro Observables)



A mesh convergence study
Analyze heavy fluid density
Errors are of two types





Position errors due to errors in locations of shock waves
Density errors in the smooth regions between the shock waves
Large reduction in error variance results from separating
these types of errors
Some statistical averages performed to achieve convergence
Convergence achieved in the sense of mean quantities
24
Convergence of mean density of heavy fluid
Heavy fluid relative density errors, at t = 60
25
Convergence of wave trajectories
  K olm ogorov length =
= size of sm allest eddy
 K m esh  K olm ogorov length /  x
 K m esh  1 : D N S
<- change in physics ->
26
Validation: Agreement between
simulation and experiment

Rayleigh-Taylor instability


Classical hydro instability, acceleration driven
Three simulation studies based on Front Tracking
show agreement with experiment




Immiscible fluids (surface tension is regularizer)
Miscible fluids high Schmidt number (viscosity and initially
mass diffused layer are regularizers)
Miscible fluids (helium-air) (mass diffusion is regularizer)
All three cases give essentially perfect agreement to
experiment on overall mixing rate.
27
Simulation compared to experiment:
Rayleigh-Taylor mixing with small levels of
mass diffusion (large Schmidt number).
Plot of h = penetration
distance of light fluid
into heavy.
Physical values of
viscosity and initial
mass diffusion are used.
The two simulations are
quarter and half of
physical problem. 28
Comparison of simulations and experiment
(Validation): For Rayleigh-Taylor unstable mixing
Immiscible fluids (surface tension)
  grow th rate
for m ixing zone
  dim ensionless
surface tension
29
Other RT Simulation Studies


50 year old problem
Most simulations disagree with experiment by
factor of 2 on overall mixing rate
  h / Agt

2
A few agree with experiment



Front Tracking (control numerical mass diffusion)
Particle methods (control numerical mass diffusion)
Mueschke-Andrews-Schilling (control numerical mass
diffusion and use experimental initial conditions)
30
Comparison of FronTier and RAGE
for 2D RM instability

T. Masser:



Macro variables not sensitive
Micro: Temperature is sensitive
Stony Brook:



Macro variables: not sensitive
Interface length divergent
Micro: Mixture concentrations sensitive
31
Major Temperature Differences at Reshock



At reshock the fingers of tin are heated to a
much higher temperature in the FronTier
simulation than the corresponding fingers in the
RAGE simulation.
There are at least three possible mechanisms
that might be responsible.

Velocity shear in FronTier missing in RAGE

Thermal and Mass diffusion at the interface
in RAGE

Differences in the hyperbolic solver
After reshock FronTier continues to have a
significantly higher maximum temperature.
32
DNS (regularized simulation) convergence of interface
lambda_C = avg. distance to exit a phase
lambda_Cmesh = lambda_C/mesh size
lambda_K = Kolmogorov length (smallest eddies)
LHS of left frame shows uniform behavior in mesh units,
independent of mesh and physics, for LES regime.
RHS of right frame shows uniform behavior relative to mesh
for fixed physics, ie mesh convergence for DNS regime.
33
Atomic Scale Mixing Properties



Atomic scale mixing properties are sensitive to
physical modeling and to numerical methods unless
fully resolved (direct numerical simulation = DNS)
Correct simulation for both micro and macro
variables: use sub grid models (large eddy
simulation = LES)
LES modify equations to compensate for physics
occurring on small scales (below the grid size) but
not present in the computation.
34
Combine strengths from two
numerical traditions

Capturing likes steep gradients, rapid time scales



Turbulence models often use smooth solutions, slow
time scales with significant levels of physical mass
diffusion


Tracking is an extreme version of this idea
Often, no subgrid model and so not physically accurate
for under resolved (LES) simulations
Often, too many zones to transit through a concentration
gradient
Best part of two ideas combined in present study
35
Subgrid models for turbulence, etc.


Typical equations have the form U t   F (U )    U
Averaged equations:
U t   F (U )    U
F (U )  F (U )
F (U )  F (U )  FSG S (U )
FSGS is the subgrid scale model and corrects for grid errors
36
Subgrid models for turbulent flow
   F U   F U 
  F U     U   F
U 
 U (key m odeling step)
U  
  F U      
 U
FSG S U
Ut
FSG S
Ut
SG S
turbulent
turbulent
37
Subgrid Scale Models (Moin et al.)



No free (adjustable) parameters in the SGS terms
Parameters are found dynamically from the simulation itself
After computing at level Delta x, average solution onto
coarser mesh. On coarse mesh, the SGS terms are
computed two ways:




Directly as on the fine mesh with a formula
Indirectly, by averaging the closure terms onto the coarse grid.
Identity of two determinations for SGS terms becomes an equation
for the coefficient, otherwise missing.
Assume: coefficient has a known relation to Delta x and otherwise is
determined by an asymptotic coefficient. Thus on a fine LES grid the
coefficient is known by above algorithm.
38
Typical atomic level observable

Mean chemical reaction rate (no subgrid
model needed for computation of w)
w  f1 f 2 e
T / TAC
T A C  A ctivation T em perature
f i  m ass fraction of species i
T 
 f1 f 2 
 f1   f 2 
;
 defined at fixed T
d (T )  probability distribution for T

w

 f1   f 2   T e
T / TAC
d (T )  m ean reaction rate
TAC
39
Re = 300. Theta(T) vs. T (left);
Pdf for T (right)    f f 
1
2
 f 1  f 2 
40
Re = 3000. Theta(T) vs. T (left);
Pdf for T (right)
41
Re = 300k. Theta(T) vs. T (left);
Pdf for T (right)
42
Kolmogorov-Smirnov Metric for
comparison of PDFs


Sup norm of integral of PDF differences
PDF data is noisy and the smoothing in
the K-S metric is needed to obtain
coherent results
x
|| p1  p 2 || K  S  ||

p1 ( y )  p 2 ( y ) dy || 

43
Convergence properties for reaction
rate pdf for w  const. f1 f 2 exp(T / T AC )
Errors for pdf for reaction rate w,
compare coarse to fine, medium to
fine and relative fluctuations in coarse
grid
Re
c to f m to f fluct. c
300
0.04
0.03
0.24
3K
0.49
0.04
0.49
300K
0.09
0.03
0.25
Conclusion:
numerical
convergence of
chemical
reaction rates,
using LES SGS
models for high
Schmidt
number flows
44
Convergence of w pdf


Not just mean converges
Moments to all order, and full distribution
converges
45
Microphysics to study combustion:
Primary breakup of fuel jet injection to engine
Parameters from diesel engine
Above: no cavitation bubbles. Below: with cavitation
46
Comparison to experiment:
simulation with and without cavitation bubbles
Velocity of jet tip
Mass flux through
Reynolds number
observation window as
vs. time
a measure of jet spreading
47
Base case: specification of numerical parameters for insertion of cavitation bubbles
3D Simulation
48
Comparison of simulation and
experiment (Argon National Laboratory)
49
Simplified engine geometry
50
Jet in Cross Flow
Mach 1 cross flow with a
2D (planar) jet of diesel
fuel.
51
Conclusions







Chaotic flow simulations are sensitive to numerical and physical
modeling
Several examples of natural physical problems
RM mixing: Temperature, concentration and chemical reaction rate PDFs
are sensitive to transport, to numerical algorithms (under resolved) and
are convergent with use of SGS model (no adjustable parameters).
RM mixing: mesh convergence including microphysical variables
(verification)
RT mixing: agreement of simulation with experiment (validation)
Solutions as measure valued distributions from the compensated
compactness theory provides a useful framework for interpretation of
simulations.
Simulations suggest that this framework might be a basis for 2D/3D
existence proofs.
52
Thank You
Smiling Face: FronTier art simulation
Courtesy of Y. H. Zhao
53