Transcript Document

A high resolutive scheme for
reactive flows
Guoxi Ni
IAPCM,Beijing,2014-5-22
With Xihua Xu,Song Jiang
1
Content
1. Model for reactive fluids
2. ALE type scheme on moving
mesh
3. GRP solvers for Fluxes
4. Numerical examples
2
1 Model for reactive fluids
Compressible model for reactive fluids
 
 t    (  u )  0,

  u
   (  uu )       F ,

 t
  E
 t    (  Eu )    (u )   u  F    (k )   Q,

where , ,Q are stress,temprature and heat 。
3
For Newton flow，
 ij  ij p  ij
 ij   (
ui u j
u
2

)  ( '  ) ij i
x j xi
3
x j

x
p  p(  , e) and
q  k
where
for fluids with reactive ，
p  p(  , e,  )
4
For reactive fluids，a reactive rate should
be given,usually,
d
 r ( p, e,  )
dt
Such as Arrhenius function
E
(
)
d
RT
 k0  e .
dt
where  is mass fraction，E is active energy。
5
In practice,some other complex reactive
rates are to be used ,for examples,
Lee-Tarver reactive rate function
d
 I (1   ) x   G (1   ) x  y p z ,
dt
Or modified reactive rate function
d
 I (1   ) x   G1 (1   ) x  y p z  G2 (1  F ) x  y p z ,
dt
6
This reactive rate cause many difficulties,
in numerical simulation,which includes
the following:
1) More complex structure in the
reactive region.
2) Different fluids in the region,
especially fluids with comples
equation of states
7
Typical landscape for ZND.
8
Standard EOS’ are likes the following
1.ideal gas
p  (  1)  e,
2.stiffen
p  ( 1) e   p ,
p  ( 1)e  c02 (  0 ),
3.Mie-Gruneisen
4.JWL
c02 0
 
p  (  1)  e 
{(  1)(  1) 
(1   )  },
 1

p  A exp(R1v)  B exp(R2v)  wcT / v,
9
Some numerical methods
W.C.Davis, foundation work,1979
D.L.Youngs, VOF
W.Bao,S.Jin,projection method
M.BenArtzi,GRP scheme,
A.J.Majda, R.J.Leveque,fraction steps
R.Fedkiw, J.J. Quirk,Ghost fluid method
R. Abgrall,Saurel,7 equation model,1999
B.N. Azarnok,T.Tang,moving meshes,2005
10
2 ALE type scheme on moving mesh
The target of the current method is to
present a discrete scheme for the reactive
fluids on moving mesh method. Different
from previous approach, the method is
constructed directly from the integral
equations of the conservation laws.
11
For example,a triangle cell move with this
velocity U g
12
Integral form governing equations for reactive fluids
d
 dt  (t ) dv   S (t )  (U  U g )  nds,

d
Udv    U (U  U g )  nds    nds,
S (t )
S (t )
 dt  (t )

d
 Edv     E (U  U g )  nds   U nds   q nds,
S (t )
S (t )
S (t )
 dt  (t )
d
   udv     u (U  U g )  nds    K (T )dv.
S (t )
 (t )
 dt  (t )
13
If U g
0
then we get the syetem for Euler method
 
 t     u  0,

  u     uu    ,
 t

  E     Eu    ( u )    ( K  )   Q,
 t
  u

    uu   K (T ),
 t
14
If U g  U then we get the syetem for (modified)
Lagrange method
dJ
 0,
 dt

 d  Ju  J    uu  0,
 dt

 d  JE  J   ( P  u )  0,

dt
 d  u

  K (T ),
dt

15
Denote U  U U g , then for mass equation
d
dt

 (t )
 dv   
S (t )
 
 (U  U g )  nds
4
S (t )
U  nds   F ,e (t ),
where
i 1
i
F ,ei (t )   U  nds,
ei
16
Momentum equation

d
dt  ( t )
Udv   U (U  U g )  nds   nds
S (t )
S (t )
   UU  nds   nds   U gU  nds
S (t )
S (t )
S (t )
4
  ( FU ,ei (t )  U g F ,ei (t ))
i 1
where
FU ,ei (t )   UU  nds   nds,
ei
ei
17
Energy equation

d
dt  ( t )
 Edv     E (U  U g )  nds     U  nds   q  nds
S (t )
 
S (t )
S (t )
 ( 12 U 2  e)U  nds     U  nds   q  nds
S (t )
    U g  nds  
S (t )
S (t )
S (t )
S (t )
U 2U g  nds  12  U g 2U  nds
S (t )
4
  ( F E ,ei (t )  U gi FU ,ei (t )  12 U g 2 F ,ei (t )),
where
i 1
F E ,ei (t )    ( 12 U 2  e)U  nds   U  nds   q  nds.
S (t )
S (t )
S (t )
18
In summary
4
d
 dt  (t )  dv   F ,ei (t ),
i 1

4
d
 dt  (t ) Udv   ( FU ,ei (t )  U g F ,ei (t )),

i 1

4
i
2
1
d
dt  ( t )  Edv   ( F E , ei (t )  U g FU ,ei (t )  2 U g F , ei (t )),

i 1

4
d
 u dv   F ,ei (t )  x   K (T ) |i .
 dt  (t )
i 1
19
3 DRP solvers for fluxes
In the following , we are going to
determine fluxes for the above discretize
system.
The strategy is to use GRP scheme for
the mass,momentum and energy
equations,and to update the reactive rate
seperately.
20
The GRP scheme is a numerical method based
on the analytic resolution of the associated
generalized Riemann problem.
1. It is a high order analytic extension of the
Godunov scheme
2. Main gradients are the Lax-Wendroff
approach plus singularity tracking
21
Ben-Artzi, Matania; Falcovitz, Joseph’ J.Comput. Phys. 55
(1984), no. 1, 1–32.
Ben-Artzi, Matania; Falcovitz, Joseph,Cambridge monographs
on Applied and Computational Mathematics, 11, 2003.
LeFloch, Ph.; Raviart, P.-A.: , Non LinWaire 5 (1988),no. 2,
179–207. ; Non LinWaire 6,(1989), no. 6, 437–480.
LeFloch, Philippe; Li, Tatsien, Asymptotic Anal. 3 (1991), no.
4, 321–340..
Li and Z. Sun, J. Comput. Phys. 222 (2007), no. 2,796–808.
M. Ben-Artzi, J. Li and G. Warnecke, J. Comput. Phys.
218(2006), no. 1, 19–43.
22
.
23
The target is to get the following values
u
p

( ) , ( )  , ( ) 
t
t
t
then we can get the numerical flux for
mass,momentum and energy which based
on the Taylor expansion of time.
24
For the reactive rate equation,consider
the contact velocity,

2,i 1/ 2 ,u  0.
i  

1,i 1/ 2 ,u  0.
then we can get the numerical flux for
reactive rate.
25
.
26
what remained is to determine the mesh
velocity in the above discrete equation, a
direct implementation of numerical
techniques on the discretized level to
prevent mesh distortion is preferred.
27
Moving mesh method (H.Z.Tang,T.Tang 2003)
E ( , ) 
1
2


(t xG1x  t yG2y )d d ,
 (G1x)  0,  (G2y)  0,
G1  G2  diag (w1 , w2 ),
Where
w1  (11 | |2 1 |  |2 )1 /2 , w2  (12 | |2 2 |  |2 ) 2 /2.
28
4 Numerical tests
Case 1 steady state problem
pressure
29
density
30
Case 2 strong detonation
pressure
31
density
32
Case 3 unsteady problem
Fixed meshes
33
Results by moving meshes
34
Case 4 two dimension unsteady problem
35
pressure
36
.
37
meshes
38
Thanks
谢谢！
39