Study of instability driven mixing via improved tracking and transport control

Study of instability driven mixing via improved tracking and transport control

Xiaolin Li, Wurigen Bo, James Glimm and the FronTier Team

Department of Applied Math and Statistics SUNY at Stony Brook

2

Outline of the talk

1. PDE, conservation law, and discontinuity 2. Improving the front tracking method

3. Comparison and benchmarks

4. Comparison of Rayleigh-Taylor instability 5. Transport control with tracking

6. Conclusion

7. Application to other scientific and engineering problems

PDE

• Hyperbolic equation (wave equation)

• Parabolic equation (diffusion equation)

• Elliptic equation (steady state equation)

Parabolic equations

In 1-D: ₂

2

x D u t

u

∂

= ∂

∂

∂ or

u

= Du

In multi-dimension:

u

= D ∆ u

∑

= ∂

∆ ^d

i ₁ xi² 2

where is the Laplace operator

Parabolic equation flattens all variation (variation deminishing.

Physically, it is the diffusion equation originated from the heat transfer equation.

Solution in infinite space

Initial condition:

1. Singularity:

2. Discontinuity:

) ( )

( x δ x

ϕ =

) ( )

( x = h x

ϕ

Solution in infinite domain

( −∞ , ∞ )

dy y

Dt e t

x

u

y x

) 4 (

) 1 ,

(

)

( ²

π ∫

ϕ

∞

−

− −

=

) ( )

0 ,

( x x

u = ϕ

Singular initial condition

 



≠

=

= ∞

= 0 0

) 0 ( )

( x

x x x δ ϕ

Dt x

Dt e t

x

u

2

4 ) 1

,

( =

π

Discontinuous initial condition





>

= ≤

= 0 0

0 ) 1

( )

( x

x x h

ϕ

dy e

t x u

Dt x

∫

∞

−

−

= 1

1 )

,

( π

Hyperbolic equation

1-D wave equation:

0 0

2

= 0 ⇒ + = − =

−

tt

a u u au u au

u

) (

) , ( )

( )

,

( x t x at u x t x at

u = ϕ − = ψ +

Traveling wave solution:

0

, =

= dt

a du dt

dx

Characteristics, along the curve

Linear and nonlinear equations

Linear equation:

a = a ( t x , )

Typical equation:

u

+ au

= 0 a = const

Nonlinear equation:

a = a ( x , t , u )

Typical equation (Burgers equation)

u a

uu

u

+

= 0 , =

Conservation law:

) ( ' ,

0 )

( u a f u

f

u

+

= =

Shock Wave

1. Shock is a result of the intersection of characteristics.

2. Shock is a discontinuity across which physics change sharply.

3. Shock speed is derived from conservation—Rankine- Hugoniot condition

s is the shock speed.

) ( )

( )]

( [ ]

[

)]

( [ ] [

L R

L

R u f u f u f u

u u

u f u

s

−

=

−

=

=

Examples of conservation law

1. Traffic Flow in a highway 2. Flood wave

3. Glaciers motion

4. Chemical exchange process 5. Oil reservoir

6. Gas dynamics

Traffic flow

0 )

( =

+

t

Q ρ

ρ

Greenberg (1959) studied the traffic of Lincoln Tunnel and found:

ρ ρ ρ

ρ a

Q ( ) = log

) (

228 ),

( 2 .

17 mph vpm

a = ρ

=

Traffic relaxation:

t 1

0

∝

− ρ

ρ

Flood wave

= 0

∂ + ∂

∂

∂

t Q t

A

A: The cross sectional area of the river bed Q: Water flux in volume

Kleitz (1858) and Seddon (1900) used balance

Between gravitational force and friction force derived:

2 / 3

3 sin

PC A g vA A

Q

f

∝

=

= α

Equation for gas dynamics

0 ))

( (

0 )

( )

(

0 )

(

2

= +

+

= +

+

= +

x t

x t

x t

P e

u e

P u

u

u ρ ρ

ρ ρ

) ,

( e

P

P = ρ

Mass, momemtum and energy conservation:

Equation of state (EOS):

Riemann problem

















 =















=





>

= <

R R

R R

L L

L L

R L

p u U

p u U

x U

x x U

U

ρ ρ

0 ) 0

0 , (



 



= 

t V x t

x U( , ) 1. Initial Condition:

2. Invariance of solution:

R R

L

L

U U U

U ,

,

,

3. Four states:

4. Three waves: left wave, contact, right wave.

Glimm Scheme

n

U j

Given states at the nth time step:

Glimm’s scheme advances the state via:

h t jh

V U

h h

t j V U

n n

j

n n

j

ϑ ς ς

ϑ ξ ξ

+

 =



 



= 

+ +

 =



 



= 

+ +

+ ++

,

) 2 / 1 (

,

1 1

2 / 1 2

/ 1

2 / 1

is a random variable in

ϑ  

 



 −

2 1 , 2

1

The convergence of Glimm scheme is through large

Number theorem and is the first significant convergence Proof for the gas dynamics equations.

Godunov scheme

2

0

/ 1 2

/ 1 1

∆ = + −

∆

−

+

x f f

t U

U

)) /

(

(

2 /

1 + +

+

=

j

f V x t

f

20

The Discrete Representation of The Front Tracking

Volume filling rectangular mesh (Eulerian Coord.)

(N-1) dimensional Lagrangian mesh (interface)

A 3D Interface A 2D Representation

Y

X (i,j)

21

The 3D interface

“Three Dimensional Front Tracking”,J. Glimm, J. Grove, X. L. Li, K. Shyue, Y. Zeng and Q. Zhang, SIAM J. Sci. Comp., 19, 1998.

Front Tracking Method

 Front tracking method is implemented in code FronTier.

Major components:

1. A moving mesh to represent interface 2. Navier-Stokes equations

3. Dynamic subgrid scale models Procedure to solve:

1. Propagate points on interface 2. Redistribute surface mesh

3. Reconstruct the tangled part in surface mesh

4. Solve equations for liquid and gas separately with ghost fluid method Interface

Fluid 2 Fluid 1

Numerical methods related to front tracking:

1. Coupling fluid solver with interface propagation 2. Handling topological changes

Ghost Fluid Method

 The ghost states on the other side of the interface is constructed by a ghost fluid method (B.Khoo et.al. 2005).

 Stencil across the interface

 Solving a Riemann problem

 Using the middle states from the Riemann problem to construct the ghost states

 2D and 3D

 Project interface normal vectors onto cell centers.

 Construct ghost states along normal directions.

 Surface tension force is modeled in the Riemann problem by

Interface Point Propagation

Fluid 1

Fluid 2

Ghost fluid 2

Fluid 2

Fluid 1

Ghost fluid 1 Real fluid states Reconstruct the left

interface state

Reconstruct the right interface state

 Interface states are reconstructed from the interpolation of real and ghost states

Propagate interface point:

Advantage: No need to solve Navier-Stokes equations on the interface.

More robust and efficient than our previous front tracking method.

A Riemann problem with s_l, ^s_r ^{as its}

data is solved to determine the interface point speed v_n

25

• The challenge to Lagrangian method

• Eulerian level set method

• The Marching cube method

• Reverse engineering: grid-based tracking

• Combining the best of Lagrange and Euler

• The locally grid-based tracking method

Geometry and Topology

Level Set Methods

• A popular and powerful scheme for interface tracking is to compute the interface position by propagating a function whose level set

corresponds to the interface position

– The interface location at time t is given by : where

– Commonly used for material interfaces

– See the books of Sethian: “Level Set Methods and Fast Marching Methods”

or Osher and Fedkiw: “Level set methods and dynamic implicit surfaces”

• Designed to handle interface topology changes automatically

– However interface are limited to shapes that can be represented by level sets

• Couples to a numerical scheme for updating flows states on a volume mesh via a ghost fluid (extrapolation across interfaces) method

• When fully developed has similar features to explicit interface methods in many aspects

Slide 26

t F 0

φ + ∇ =φ

( )^{,t 0}

φ x =

The idea of grid-based

untangling

Grid-based Front Tracking

1. The common agreement: interface is greatly simplified in Eulerian grid.

2. Marching Cubes, Lorenson and Cline, 1987, (Static, Computer Graphics).

3. Level set method, Osher and Sethian, 1988, (Implicit).

4. Grid-based front tracking, SJSC, 21, 6 2000, (Explicit and Dynamic).

The 14 isomorphically distinct cases

Grid-based topological correction

Grid-based Tracking

Interface bifurcation under grid-based front tracking method

“Robust Computational Algorithm for Dynamic Interface Tracking in Three dimensions”, J. Glimm, J. Grove, X. L. Li and D. C. Tan, SIAM J. Sci. Comp., 21, 2000.

Basic FronTier Test Simulations

Interface Topological Changes

 Grid based tracking is robust but too diffusive.

 Challenge: Robustness of the algorithm is crucial for large scale computing.

Grid based tracking Grid free tracking

Interface Topological Changes

 Algorithms to handle topological changes

 Grid free tracking (GF)

 Grid based tracking (GB)

 Locally grid based tracking (LGB)

Tangled

interface GF

LGB GB

Robust Locally Grid Based (LGB) Untangle

 Advantage

 Local, it is suitable for large scale computing.

 Robust, It generates topologically valid surface mesh.

 A robust algorithm to reconnect a grid based surface mesh with a grid free surface mesh

“A Simple Package for Front Tracking”, J. Du, B. Fix, J. Glimm. X. L. Li, Y. Li, L. Wu, JCP, 213, 2006.

37

Benchmark Plus

38

3D rotation of slotted sphere

39

Fifth order level set (WENO) vs. fourth order front tracking (Runge-Kutta)

40

Front tracking reversal test of interface in the deformation velocity field

41

Resolution Test

42

Front tracking reversal test of interface in 3D deformation velocity field

64 64

64× ×

128 128

128× ×

0 .

= 0

t t =1.5 t = 3.0

43

Topological bifurcation: it’s done!

44

模拟维模拟模拟晶体结晶过程

Topological merging of 3D surface mesh

Examples

 Interface bifurcation and merging are commonly observed in multiphase flow

mesh bifurcation in a curvature dependent surface propagation

mesh merging in a droplet collision simulation

46

Conservative Front

Tracking

47

The extended stencil method

For ghost-cell scheme:

(

)

n j n

j F F

x u t

u ₊⁺₁¹ _{+1 +3/2} − _+1/2

∆

− ∆

=

(

)

1 + −

+ −

∆

− ∆

= ⁿ_{j j}^L _j

n

j F F

x u t

u

) ,

, , (

) ,

, , (

2 1

1 2

/ 1

2 1

1 2

/ 1

n j n

j n j n

j R

j

n j n

j n

j n

j L

j

u u

u u

F F

u u

u u

F F

+ +

− +

+ +

− +

=

=

R j L

j

F

F

≠

Previous works

1. Chern and Colella, LLNL Report,1987 2. D. K. Mao, JCP, 1991, 1993

3. Pember, JCP, 1995

Neglect higher order term and note that t

dS uv udV

SM

n V

∆

=

∫

∆

∫

We have the integral form of conservation 0 )

( =

+

∂ −

∂_t

∫

∫

∫

S n s

n

V _M

This can also be written as

0 )

) ( ( )

( + − =

∂ +

∂_t

∫

∫

∫

S

n S

n

V _{F M}

Or simply

0 )

) (

( − =

∂ +

∂_t

∫

∫

S n V

Conservative Interface-Interior Coupling

0 )

( =

+

t

f u

u

The conservation law:

The Rankine-Hugoniot condition:

R R

L

L

su f u su

u

f ( ) − = ( ) −

( )

(

)

n j n

j n

j n

j n

j

n j L

n I n

j n

j n

j n

j

F F

t U

x U

x

F F

t U

x U

x

, 2 / 1 2

/ 1

2 / 3 1

1 1

1 1

1

2 / 1

2 / 1 ,

2 / 1 1

1

+ +

+ +

+ +

+ +

+

−+ +

+ +

−

∆

−

∆

=

∆

−

∆

−

∆

=

∆

In 1-D, this lead to:

where

R n

I L

n I

n R n

n R R

n I

n L n

n L L

n I

F F

u s

u f F

u s

u f F

, 2 / 1 ,

2 / 1

2 / 1 2

/ 1 2

/ 1 ,

2 / 1

2 / 1 2

/ 1 2

/ 1 ,

2 / 1

) (

) (

+ +

+ +

+ +

+ +

+ +

=

−

=

−

=

Due to Rankine-Hugoniot condition

1D Conservative Front Tracking Geometry

Two cases

• Fronts do not cross the cell center in one time step.

• Fronts do cross the cell center in one time step.

New cell average and : v_iⁿ v_iⁿ₊₁

∫

∫

+

−

= −

= −

+ +

−

2 / 3 2 / 1

) 2 (

/ 3 1

) (

2 / 1

) , )) (

( (

1

) , ) (

) ( (

1

i

n n

i

x

t n

n i

n i

t

x n

i n

n i

dx t

x t U

v x

dx t x x U

v t

σ σ

σ σ

n (mesh)

shock position error

order

error

L1 order

50 4.8e-2

100 1.3e-4 8.5

3.8 200 4.3e-5 1.6

1.4 400 1.6e-6 1.4

1.6

Convergence test of conservative tracking

In multi-dimensional case, we consider the time space equation:

ˆ 0 ) ˆ (

) ˆ (

)

ˆ + ( + + ∇⋅ =

= un f u n g u n h u n F

F _{t x y z} 

On a time-space cell

ˆ = 0

∫

TS

tsdS n

F

Cell-merge is needed if the volume of the time-space cell is less than half of the regular cell, in 2-D the time space cell is constructed the same way as the 3-D grid based interface.

The time-space interface between n and n+1 time steps

Before cell merger

After cell merger

Conserva

Tracking

Nonconsertive Tracking Mass

Error 0.0 0.21%

X-Mom

Error 0.0 0.21%

Energy

Error 0.0 0.21%

62

Rayleigh-Taylor Instability

Inertial Confinement Fusion (ICF)

64

FronTier application: chaotic mixing

65

FronTier application: chaotic mixing

Chaotic mixing is not only important to ICF, but also a test of large scale FronTier

application to petascale computing. We have implemented a load balanced parallel algorithm and ran up to 1024 processors on New York Blue. Collaboration with B. Cheng, John Grove, and D. Sharp at LANL.

66

3D Turbulent Mixing

67

Acceleration Driven Mixing

• Rayleigh-Taylor (RT), steady acceleration:

2 2 1

2 1

;

h α Agt A ρ ρ ρ ρ

= = −

+

The Paradox α

Agt

h

= α

David Youngs and K.Read

(1984)

Read’s Experiment (1984)

3D alcohol/air

3D NaI soln./Pentane 3D NaI soln./Hexane

Exp # 29 39 58 Exp # 33 35 Exp # 62 60

Alpha = 0.073 0.076 0.077 Alpha = 0.066 0.071 Alpha = 0.063 0.073

Summary of Experiments

The Alpha of Bubbles

FronTier:

Alpha = 0.08

TVD:

Alpha = 0.025-0.045

72

FronTier TVD

Agt = 25.3 h = 4.16 Density plot ²

73

74

Goal of mixing study

• Predict large scale features. Size of mixing zone

• Predict statistics (means, variances) of fluid quantities

• For use in combustion

– Predict full probability distribution (PDF) of species concentration and temperatures

• Accurate models down to atomic level of mix are needed

• These must be sensitive to transport,

Reynolds number, Schmidt number

75

Real vs. Ideal Mixing Physical vs.

Numerical Scale Breaking

• Numerical nonideal effects

– Numerical surface tension – Numerical mass diffusion

• Physically nonideal effects

– Surface tension

• Surfactants, variable surface tension, Marangoli force

– Mass diffusion

• Initial or for all times

– Viscosity

– Compressibility

– Long wave length initial perturbations

76

Main New Results

• Systematic agreement of simulation with experiment and theory

• Alpha, bubble width, bubble height fluctuations – Most relevant experiments included in

agreement

• Reed-Youngs, Smeeton-Youngs, Andrews experiments

• Omitted:

– Immiscible with surfactant (Dimonte and Smeeton- Youngs)

– Initial diffusion layer (in progress) – Subgrid models

– Nonideal initial conditions

77

Scale Breaking: Experiments and Simulations

Scale breaking physics

Alpha-

experimental

Alpha-

simulation

Experiment s

Fluids

Surface tension 0.050-0.077 0.067 RY, SY Liquid/liquid;

liquid/gas Surface tension

with surfactant

0.050-0.061 ???? SY,DS Liquid/liquid

Mass diffusion 0.070 0.069 BA Gas/gas

Initial mass diffusion

0.062 ???? SY Liquid/liquid

Viscosity 0.070 ???? SA Liquid/liquid

Compressibility Up to 0.2 Lasers plasmas

78

Comparison of Mixing Rates:

Comparison, Simulation and Theory

Theory Experiment Simulation

height 0.06 0.067 0.062

radius 0.01 0.01 0.01

fluctuations in height

0.028 0.034

79

Turbulent Mixing

• Acceleration driven mixing

– Steady acceleration – Rayleigh-Taylor mixing

– Impulsive acceleration – Richtmyer-Meshkov mixing

• Most RT computations underpredict mixing rates relative to experiments

– Simulation analysis using time dependent densities (Atwood number) makes this point

• Cause appears to be numerical mass diffusion, which reduces the local density contrast and thus the large scale mixing rates

– Numerical surface tension also significant

• Questions raised about the role of initial noise in the experiments

80

Numerical Non-Ideal Effects

• Numerical mass diffusion

– Removed by tracking

– Errors modify density contrast by a factor of 2 for typical grids

• Numerical surface tension

– Reduced by local grid based (LGB) tracking – Errors proportional to curvature x Delta x

– Arises from approximation of interface by a line segment within each mesh block

– Arises from grid level description of interface and thus occurs for all untracked methods

81

Time Dependent Atwood Number

• Atwood number

• For each z

– Compute the maximum and minimum density

– Form a space (height) and time dependent A(z,t) from min/max

• Average A(z,t) over bubble region to get A(t)

• Untracked A(t) is about ½ nominal value due to mass

diffusion; (incompressible) tracked A(t) is virtually constant

• If A(t) is used in definition of alpha, all low compressible simulations agree (with each other, with experiment, with theory)

• If A(t) is used in compressible simulations, all simulations are self similar, but self similar growth rate depends on compressibility

2 1

2 1

A ρ ρ

ρ ρ

= −

+

82

Physical Non-Ideal Effects

• Viscosity, mass diffusion, surface tension

• Compressibility

– Solution depends on initial temperature stratification; assume isothermal. Initial density depends on height z, so that Atwood number is z dependent.

– Data interpretation using a time dependent Atwood number restores self similarity, but the mixing rate alpha increases with

compressibility.

83

Turbulent Mixing: Predictions of Gross Features (Mixing Rate alpha, etc.)

• Systematic agreement of theory, simulation and experiment for RT turbulent mixing

– Scale breaking physics important to this agreement

• Tracking is the best of practical interface methods

– Control over numerical mass diffusion and numerical surface tension needed for RT agreement

• Validation studies still in progress (viscosity)

84

Other Applications of

Front Tracking

85

模拟三维内燃机喷嘴

Ask not what the earth can do for us, ask what we can do for the earth

American consumes about 200 billion gallons per year, a 10% saving will be 20 billion gallon amounts to more than 40 billion dollars, not to mention the benefit to the environment.

86

3D Simulation of a Real Fuel Injection

 All parameters are from an experiment performed by Parker*

nozzle radius (R) 0.1mm

grid 20/R

fuel density 0.66 g/cm³ gas density 0.0165 g/cm³ fluid viscosity 0.013 Poise surface tension 24 mN/m²

Reynolds number 20,300 Weber number 2.2×10⁶ Ohnesorge number 0.073

Density ratio 40

* P. Parker, Atomization and Sprays, 8, (1998)

87

Verification: Rayleigh Instability

Number of cells on radius

FT/GF M (2D)

FT/GFM (3D) 5

10 20

0.1396 0.0607 0.0321

0.2853 0.1702 0.0672 The relative errors of the growth rate

Comparison with the dispersion relation dispersion relation

Incompressible fluid equation

0

=

⋅

∇

∇ +

−∇

=

⋅

∇ +

u

u u

u

u _t ( ) p ν ²

This is a mixed hyperbolic and elliptic equation U is the velocity of fluid and p is the pressure.

89

Incompressible Rayleigh-Taylor instability on Reynold number

(from left: 14,140,1400)

90

Incompressible code in 3D

91

模拟维模拟模拟晶体结晶过程

Two Dimensional Solute Precipitation

92

模拟维模拟模拟晶体结晶过程

Three Dimensional Solute Precipitation

93

Dissolution is the opposite process of

deposition

94

The Spring Model for two dimensional surface

:

The X Parachute

96

风力发电机的数值模拟

97

Fluid-Rigid body interaction

98

Simulation of Cell Migration

99

American and other exotic options

( )

1 , )

0 , (

, 0 ) 0 , (

, 0 )

2 ( 1

2 2 2 2

∂ =

∂

−

=

>

−

=

≤

=

−

=

=

−

∂ +

− ∂

∂

− ∂

∂

∂

S C

E S C

E S E

S S

C

E S S

C

t T

C S D

S C S

S C C

f f

τ

γ σ

τ γ

The Black-Schole Equation:

The interface

Condition at all time:

Initial Condition:

American and many other exotic options are PDE free boundary problems. Front

tracking provides an accurate tool to solve the basket hedging problems. We have already established

1-D and 2-D computational platform for such problems.

100

One Dimensional American options

Front tracking on 1-D American call option

Front tracking on 1-D American put option

101

World’s fastest computers (top 500)

From mega scale, peta scale, to exa sacle

102

模拟维模拟模拟晶体结晶过程

Parallelization of Front Tracking

103

Parallel load balancing

Like AMR, FronTier has encountered great obstacle in load balancing and parallel scaling. One important development is adaptive partition load balancing.Up to 8196 processors have been tested. No better for number larger than that.

Parallel Performance of FT

Grid Partition nCores Time to

solution(s) Ideal(s) 256×256×128 16×16×8 2048 157.1 157.1 256×256×256 16×16×16 4096 157.5 157.1 256×256×512 16×16×32 8192 158.2 157.1 256×256×1024 16×16×64 16384 159.8 157.1 Performance of LGB

 Jet simulation

 300-3million Triangles

 Bluegene/L 4096 cores

Weak scaling

 Rayleigh Instability

 Bluegene/L

105

A quotation from Albert Einstein

1. Stony Brook, AMS Department, galaxy cluster (over 500 processors)

2. Stony Brook, CEAS, Seawulf cluster 3. New York Blue:

103.22 teraflops

Major Computing Resources:

106

Thank you for your attention