Transcript A low Mach compressible two-fluid model with phase

A COMPRESSIBLE MODEL FOR
LOW MACH TWO-PHASE FLOW WITH
HEAT AND MASS EXCHANGES
N. GRENIER, J.P. VILA & Ph. VILLEDIEU
CONTEXT AND MOTIVATION
Context : COMPERE program from CNES &
DLR
External Heat flux
Gas phase
Evaporation
Research partners : ONERA, ZARM, CNRS,
Erlangen university, Air LIquide (Grenoble),
Astrium ST (Bremen)
Capillary raise
Liquid phase
Objectives : development of numerical tools for
simulating complex fluid behavior inside space
Separated two-phase launcher tanks:
flow with a free
• dynamical behavior  sloshing
moving interface
• thermal effects  heat and mass exchanges
• low gravity effects  capillary effects,
Marangoni convection …
2
MULTIMAT 2011 - 5-9 septembre 2011
OUTLINE OF THE PRESENTATION
• 1. Presentation of the model
• 2. Numerical method
• 3. Numerical test cases
• 4. Conclusion
3
MULTIMAT 2011 - 5-9 septembre 2011
OUTLINE OF THE PRESENTATION
• 1. Presentation of the model
• 2. Numerical method
• 3. Numerical test cases
• 4. Conclusion
4
MULTIMAT 2011 - 5-9 septembre 2011
1. Presentation of the model
Modeling choices
• Two fluid model  diffuse interface model
• Advantage : not necessary to localize (level set method) or reconstruct (VOF
method) the interface between the two fluids  easy to implement
• Drawback : interface diffusion  necessary to define a “mixture” physical model
and to use low diffusive numerical scheme
• Compressible model
• Advantage : more general, easier to implement into a gas dynamics code (ONERA
context)
• Drawback : ill conditioned for low Mach number flows  low Mach Scheme
• Same velocity field for both fluids
• Advantage : hyperbolic model, no closure assumption needed
• Drawback : impossible to deal with subscale phenomena (subgrid bubbles or
droplets …)
5
MULTIMAT 2011 - 5-9 septembre 2011
1. Presentation of the model
Inviscid two-fluid Model
  g
    g v   0


t

 
Liquid bulk density
 t      v   0
(1) 
  v      v  v  pI    g
 t
  E

     Ev  pv    g.v
 t
Gas bulk density
with  E   e 
  g  l
1 2
 v being the mixture total energy per unit volume and
2
the mixture bulk density.
To get a close model, it is now necessary to give a relation between the
“mixture” pressure p, the bulk densities  g ,l and the mixture specific
internal energy e.
6
MULTIMAT 2011 - 5-9 septembre 2011
1. Presentation of the model
Extension to non isothermal flows
Let T denote the mixture temperature and  the gas volume fraction. 
and T are assumed to be the unique solution of the following system :
(2)
  g 
 

p
,
T

p
,
T
 g



1


 




g
l


e


e
(
,
T
)


e
(
,T )
g g
l l


1
Local mechanical equilibrium
Local thermal equilibrium
where p = pg(g,T), p = pl(l,T) denote the gas and liquid EOS and
e=eg(g,T), e=el(l,T) denote the gas and liquid colorific laws.
The mixture EOS is then (implicitly) defined by :
 ρg 
 ρ

(3) p(ρ g ,ρ ,  e)  pg  , T   p 
,T 
 1 


7
MULTIMAT 2011 - 5-9 septembre 2011
1. Presentation of the model
Other interpretation of closure equations (5)-(6)
Closure relations (2)-(3) can also be interpreted as a direct consequence of the
following modeling assumption for the mixture Gibbs potential :
g mix ( p, T )  yg g g ( p, T )  yl gl ( p, T )
with
yg 
g

, yl  l


Ideal mixture assumption
Important consequence : System (1) with pressure law given by (2)-(3) is
thermodynamically consistent in the sense that it has a convex entropy in the
sense of Lax defined as :

 (  g , l ,  E ,  v)    g s g (  g ,  g eg )  l sl ( l , l el )


Φ(  g , l ,  E ,  v)   (  g , l ,  E )V
where eg, el are the gas and liquid specific internal energies (implicitly
defined by the solution of (2)), sg and sl are specific entropies, and


are the real densities.
 g = g , l = l

1
8
MULTIMAT 2011 - 5-9 septembre 2011
1. Presentation of the model
Inclusion of diffusion and capillary effects.
  g
    g v   0


t

 
 t      v   0
(4) 
  v      v  v  pI    g  div( τ  τ )
v
c
 t
  E

     Ev  pv    g.v  div(φc )  div( τ v : v) + div( τ c : v)
 t
with :
τv   div(v) I + (v  t v)

   
τ c   (T )   I 

2





φc   T
Viscous stress tensor
Capillary stress tensor
(body force formulation)
Heat flux
9
MULTIMAT 2011 - 5-9 septembre 2011
1. Presentation of the model
Approximate enthalpy equation for low Mach flows
Neglecting viscous and capillary effects, the energy equation is equivalent
to :
 h
p
     hv      T  +
 v.p
t
t
Heat flux
presssure contribution
which is the Eulerian formulation of the well-known thermodynamic relation :
dh   q 
1

dp
For low Mach number flow, with imposed pressure on one of the boundaries,
one generally has : 1/ dp << q, and therefore the energy equation can be
replaced by the heat equation :
 h
     hv      T   0
t
10
MULTIMAT 2011 - 5-9 septembre 2011
1. Presentation of the model
Phase change modeling
Phase change phenomena can be included in model (7) by just adding
a relaxation source term in the r.h.s. :
U
 (U )  U
   F(U , U ) 
 S (U )
t

where (U) is the thermodynamic equilibrium state corresponding to U,
defined as the state which maximizes the mixture entropy under the
constraints of imposed total volume, total mass and total energy :
Max yg sg (eg , vg )  yl sl (el , vl )
yg  yl 1
yg eg  yl el  e
yg vg  yl vl  v
1
, e  e(U )
where v 
 (U )
In practice, the thermodynamic equilibrium time scale is assumed to be
infinitely small compared to the macroscopic time scale  local thermodynamic
equilibrium assumption.
This idea was first proposed in : HELLUY P., SEGUIN N., “Relaxation model of phase transition flows”, M2AN, Math. Model. Numer. Anal., vol. 40,
num. 2, 2006, p. 331–352.
11
MULTIMAT 2011 - 5-9 septembre 2011
OUTLINE OF THE PRESENTATION
• 1. Presentation of the model
• 2. Numerical method
• 3. Numerical test cases
• 4. Conclusion
12
MULTIMAT 2011 - 5-9 septembre 2011
2. Numerical scheme
A finite volume relaxation scheme
Each time step is divided in two stages :
Transport step  Eulerian finite volume scheme
U U
*
K
n
K
t

mK

eK
Gen1/ 2 me  t S Kn1/ 2
K
e
Numerical flux on edge e
ne,K
Ke
Relaxation step  local thermodynamic equilibrium
U Kn1  (U K* )
Note that, by construction, the second step is entropy diminishing.
13
MULTIMAT 2011 - 5-9 septembre 2011
2. Numerical scheme
Expression of the hyperbolic numerical flux
(isothermal case  no energy equation)
0



G(U L ,U R , ne )  ve+U L  ve-U R   0

pn 
advection term
 e e
Low Mach Scheme
Which expression choosing for pe and ve ?
pe =
1
2
1
ve =
2
 pL  pR 
 v L  v R  .n e
e
ne
UL
UR
pressure term
Centered scheme for pressure (see Dellacherie
(2011) recent work on low Mach number schemes)
-  e  pR  pL 
Centered expression + stabilizing
pressure term. Expression of the
positive parameter e will be
given later.
Remark : a similar idea has been proposed by Liou (AUSM+up scheme, JCP 2006) and by Li & Gu (all Mach Roe type
scheme, JCP 2008) for the compressible gas dynamics system.
14
MULTIMAT 2011 - 5-9 septembre 2011
2. Numerical scheme
Semi-implicit version of the scheme
(isothermal case)
To avoid a restrictive stability condition based on the sound celerity, mass
conservation equations are solved with a implicit scheme.
  g* , K 
 * 
 
 l ,K 
  gn, K  t
 n 
   mK
 l ,K 

+

v

e

eK

  g* , K  -   g* , K 
 *   ve  * 
 
 
 l ,K 
 l ,K 

 me

Newton
algorithm
An explicit scheme is used to compute the new velocity from the
momentum equation :
 K* v*K  Kn v nK 
with
ve =
1
2
t
mK
 v
eK

e
*
K* v nK + ve Ke
v nKe  pene,K  me
 v nK  v nKe  .ne -  e  pK*  pKe* 
and
pe =
1
2
15
MULTIMAT 2011 - 5-9 septembre 2011
p
*
K
*
 pKe

2. Numerical scheme
Formal justification of the stabilizing role of “-  (pR - pL) “
Let us consider the following modified system for isothermal flows
Modified convective velocity
 ρ
 t     ( v   v)ρ   0
 g 
(1') 
where ρ   
 l 
  v     ( v   v)  (  v)  pI   0
 t
Proposition : the term  v   h p has a stabilizing effect in the sense for that any
smooth solution of (1’) one has the following free energy balance equation :

1

1


2
2

F


v

di
v

F


v

p
(
v


v
)



  p. v

t 
2
2




Dissipative source term if v
is proportional to – grad(p)




g
l
 + l fl 
 denotes the free energy
where  F (  g , l )   g f g 

(

,

)
1


(

,

)
g
l 
g
l 


of the mixture.
Remark : the same property holds for the non isothermal case but with the entropy
instead of the free energy.
16
MULTIMAT 2011 - 5-9 septembre 2011
2. Numerical scheme
How to choose the value of e ?
Stability theoretical result : Under the two following conditions
(i)
 t mK
t mK 
 e  Max 
,
*
* 
2
m

2
m

K
K
K
K 


e
e
(ii)
2t v K mK
1
Sup m
K
K
e
with
vK 
CFL like condition
1
 K* mK

*
Ke
ve me
e
mean cell velocity
the semi-implicit scheme is entropic (in the sense of Lax).
In practice, we take :
 t mK t mKe
1
 e  Max 
,
n
n
2
m

 K K mKe  Ke
 t ( vnK  cKn )mK 

  cfl
 and Sup 
m



K
with cfl much larger than 1 for low mach number flows.
17
MULTIMAT 2011 - 5-9 septembre 2011
2. Numerical scheme
Discretization of the enthalpy equation
To respect the maximum principle on the temperature, we use the following
upwind scheme based on the sign of the mass fluxes :
 K* hK*   Kn hKn 
t
mK
 (h
*
g ,K
max( g* ,e , 0)  hg* , Ke min( g* ,e , 0))me 
e
t
mK
 (h
*
l ,K
max(l*,e , 0)  hl*, Ke min(l*,e , 0))me
e
*
*


Where g ,K , respectively l ,K , denote the gas, respectively the
liquid, mass numerical flux.
In practice, two variants of the scheme can be used :
• an explicit scheme with respect to the fluid temperature,
• a fully implicit scheme with respect to all thermodynamic
variables  g , l , T
18
MULTIMAT 2011 - 5-9 septembre 2011
2. Numerical scheme
Relaxation step U*  Un+1
, v and h are left unchanged during this step. We thus have :
n 1
n 1
*
*
*










l
g
l
 g
 * *
n 1
n 1
n 1

h


h
(
p
,
T
)



g
g
l hl ( p, T )

If both phases can coexist (gas – liquid thermodynamic equilibrium)
p(ln1,gn1,T n1 )=psat (T n1 )
System of 3 equations
and 3 unknowns
else only one phase can be present in the cell at the end of the time step
 gn1  * , ln1  0 or ln1  * , gn1  0
Remark : in practice, for numerical purpose, a minimal lower value is imposed for
gas and liquid mass fractions.
19
MULTIMAT 2011 - 5-9 septembre 2011
OUTLINE OF THE PRESENTATION
• 1. Presentation of the model
• 2. Numerical method
• 3. Numerical test cases
• 4. Conclusion
20
MULTIMAT 2011 - 5-9 septembre 2011
3. Numerical test cases
Linear oscillations in a 2D rectangular tank
•ρ1=1 kg.m-3 ; c1=300 m.s-1
•ρ2=1000 kg.m-3 ; c2=1200 m.s-1
•Transverse acceleration : a0 = 0.01 g
•Coarse cartesian grid : 40 X 20
•Ma = 2 10-5
2
Possibility to compute an analytical solution
as a série expansion by potential flow theory.
(see for example Landau & Lifschitz T6, fluid
Mechanics)
21
MULTIMAT 2011 - 5-9 septembre 2011
3. Numerical test cases
Linear oscillations in a 2D rectangular tank
Exact solution
Numerical
Scheme
Godunov
scheme
Low Mach
scheme
Time step
0.0005
0.05
Total CPU
time
Second order low Mach scheme
Second order Godunov type scheme
22
MULTIMAT 2011 - 5-9 septembre 2011
20
1
3. Numerical test cases
Dynamical test case :
bubble rise inside a liquid : Sussman et al test case
(Sussman, M. and Smereka, P. and Osher, S., A Level Set approach for computing solutions to
incompressible two-phase flow, Journal of Computational Physics, 114, 146-159, 1994)
Explicit
Godunov type
scheme with
real EOS
Explicit
Godunov type
scheme with
modified EOS
Semi-implicit
low Mach
scheme with
real EOS
Cartesian mesh
Cartesian mesh
140 X 233
Cartesian mesh
140 X 233
140 X 233
23
MULTIMAT 2011 - 5-9 septembre 2011
3. Numerical test cases
Bubble rise inside a liquid : Sussman et al test case
(Sussman, M. and Smereka, P. and Osher, S., A Level Set approach for computing solutions to
incompressible two-phase flow, Journal of Computational Physics, 114, 146-159, 1994)
Semi-implicit
low Mach scheme
with real EOS
Sussman et al
Solution with
Level Set
method and
incompressible
model
Usual Godunov
type scheme with
real EOS
24
MULTIMAT 2011 - 5-9 septembre 2011
3. Numerical test cases
Rayleigh-Bénard instability
Wall with
imposed
temperature
Periodic
boundary
conditions
Gas phase
g
Liquid phase
Ra  g 
TH 3

with

1   
  T  p
Critical Rayleigh number for instability : Rac  1708
25
MULTIMAT 2011 - 5-9 septembre 2011
3. Numerical test cases
Rayleigh-Bénard instability
Stable
Unstable
Unstable
Stable
Stable
Unstable
26
MULTIMAT 2011 - 5-9 septembre 2011
3. Numerical test cases
Marangoni convection test case
Adiabatic Wall
Gas phase
Wall with
imposed
temperature
T = T0
Liquid phase
Wall with
imposed
temperature
T = T1<T0
Adiabatic Wall
No gravity. Static contact angle : q = 90°
27
MULTIMAT 2011 - 5-9 septembre 2011
3. Numerical test cases
Marangoni convection test case
Coarse grid
Medium grid
Fine grid
Temperature
field
Volume fraction
field
28
MULTIMAT 2011 - 5-9 septembre 2011
3. Numerical test cases
1D Evaporation test case
Evaporation front
Outlet with
imposed
Liquid phase
pressure :
p = p0
=l , p = p0 , T = Tsat(p0),
u = uI
Gas phase
=v , p= p0 , u= 0,
T = f(x)
Approximate theoretical solution
qw
m
Lv
uI  
m
v
 v 
ul  u I  1  
 l 
29
MULTIMAT 2011 - 5-9 septembre 2011
Wall with
imposed
heat flux
qw
3. Numerical test cases
1D Evaporation test case
l  1000kg m3
v  1 kg m3
Interface position vs time
for several values of Lv and qw.
30
MULTIMAT 2011 - 5-9 septembre 2011
3. Numerical test cases
1D Evaporation test case
31
MULTIMAT 2011 - 5-9 septembre 2011
OUTLINE OF THE PRESENTATION
• 1. Presentation of the model
• 2. Numerical method
• 3. Numerical test cases
• 4. Conclusion
32
MULTIMAT 2011 - 5-9 septembre 2011
CONCLUSIONS AND FUTURE PROSPECTS
• An Eulerian two-fluid model with diffuse interface has been applied to the
simulation of low Mach separated two-phase flows with heat and mass
transfers.
•Using formal arguments, a simple semi-implicit low Mach scheme has been
proposed for this model. For isothermal flows, this scheme has been proved to
be entropy diminishing under a CFL condition which do not depend on the
sound celerity.
• This methodology can be very easily implemented in existing industrial
compressible CFD codes for multi-physics applications (work in progress at
ONERA). It is a very interesting alternative to classical approaches based on
one-fluid incompressible model with VOF or Level Set methods.
33
MULTIMAT 2011 - 5-9 septembre 2011
CONCLUSIONS AND FUTURE PROSPECTS
•This two-fluid approach has been successfully applied to several academic
problems for low Mach two-phase flows.
• Future works will be devoted to the
• assessment of the method for more complex phase change problems.
• extension of the model to more complex physical problems : multicomponent gas phase with an incondensable specie, 3 phases
problems …
• parallelization of the code for 3D applications
34
MULTIMAT 2011 - 5-9 septembre 2011
Thank you for your attention …..
35
MULTIMAT 2011 - 5-9 septembre 2011
Back – up
36
MULTIMAT 2011 - 5-9 septembre 2011
1. Presentation of the model
Purely Dynamical model
(inviscid Isothermal flow)
  g
    g v   0

 t
 
(1) 
    v  0
 t
  v
 t      v  v  pI    g

with
  g  l the mixture bulk density.
To get a close model, it is necessary to give a relation between the
“mixture” pressure p and the bulk gas density  g and the bulk liquid
density l .
Remark: The gas volume fraction  is not explicitly transported in this model.
37
MULTIMAT 2011 - 5-9 septembre 2011
1. Presentation of the model
Purely dynamical model
• Let  denote the gas volume fraction : g  g ; l  (1   )l
 is defined as the solution of
(2)
 g
pg 


  

p



1





Local pressure equilibrium between
the two non miscible fluids
where p = pg(g) and p = pl(l) denote the gas and liquid equation of state.
Mixture EOS
 ρg
(3) p(ρ g , ρ )  pg 


 ρ 

p



1





Remark : if the expressions of pg and pl are complex, p is only implicitly defined in
function of the bulk densities.
38
MULTIMAT 2011 - 5-9 septembre 2011
1. Presentation of the model
Example : stiffened gas model for both fluids.
Expression of the Gibbs potential for each fluid :
gi ( p, T )   i cvi T 1  ln(T )  ( i 1)cvi T ln( p  i )  ui0  Tsi0
Fluid i Equation of state
1
T
 vi ( p, T )  ( i  1)cvi
 i ( p, T )
p  i
Fluid i calorific law
p   i i
ei ( p, T )  ei0  cvi T
p  i
With these notations, system (5)-(6) is equivalent to :
Mixture specific
volume
g

l
1
v


v
(
p
,
T
)

vl ( p, T )
g

 g  l  g  l
 g  l

(7) 
 e   g e ( p , T )   l e ( p, T )

 g  l g
 g  l l

39
MULTIMAT 2011 - 5-9 septembre 2011
1. Presentation of the model
Other interpretation of closure equations (5)-(6)
Closure relations (5)-(6) can also be interpreted as a direct consequence of the
Ideal mixture assumption
following modeling assumption for the mixture Gibbs potential :
g mix ( p, T )  yg g g ( p, T )  yl gl ( p, T ) - T ( yg , yl )
with
yg 
g

, yl  l


Important property : System (4) with pressure law given by (5)-(6) is
thermodynamically consistent in the sense that it has a infinite set of convex
entropies in the sense of Lax defined as :
 g l


(

,

,

E
,

v
)



s
(

,
e
)


s
(

,
e
)


(
, )

g
l
g
g
g
g
l l
l
l



Φ(  ,  ,  E ,  v)   (  ,  ,  E )V
g
l
g
l

where  is an arbitrary concave function, eg, el are the specific internal
energies, implicitly defined by the solution of (5), sg and sl are the specific


entropies, and
are the real fluid densities.
 g = g , l = l

1
40
MULTIMAT 2011 - 5-9 septembre 2011
1. Presentation of the model
Purely dynamical model (3/3)
Proposition : If pg and pl are strictly non decreasing functions, model
(1)-(2)-(3) is hyperbolic and has a convex entropy in the sense of Lax
defined as :

 g 
1
 l 
2
 (  g , l ,  V)   V   g f g   + l fl 

2

1





 
Φ(  ,  ,  V )   (  ,  ,  V )V  pV
g
l
g
l

Lax entropy (convex function
of the conservative variables)
Entropy flux
where fg and fl are the free energy of the gas and liquid phases and are
defined as :
pi (  )
pi
fi ( ) 
d
dfi   pi d i  2 d i
2
i


41
MULTIMAT 2011 - 5-9 septembre 2011
4. APPLICATIONS
Linear oscillations in a 2D rectangular tank
•ρ1=1 kg.m-3 ; c1=300 m.s-1
•ρ2=1000 kg.m-3 ; c2=1200 m.s-1
•Transverse acceleration : a0 = 0.01 g
•Coarse cartesian grid : 40 X 20
•Ma = 2 10-5
2
Possibility to compute an analytical solution
as a série expansion by potential flow theory.
(see for example Landau & Lifschitz T6, fluid
Mechanics)
42
MULTIMAT 2011 - 5-9 septembre 2011
4. APPLICATIONS
Linear oscillations in a 2D rectangular tank
Exact solution
Numerical
Scheme
Godunov
scheme
Low Mach
scheme
Time step
0.0005
0.05
Total CPU
time
Second order low Mach scheme
Second order Godunov type scheme
43
MULTIMAT 2011 - 5-9 septembre 2011
20
1
1. Presentation of the model
References
• R. ABGRALL, R. SAUREL. A simple method for compressible multifuid flows, SIAM J. Sci.
Comput. 21 (3) : 1115-1145, (1999). 66
• G. ALLAIRE, G. FACCANONI et S. KOKH, A strictly hyperbolic equilibrium phase transition model.
C. R. Acad. Sci. Paris Sér. I, 344 pp. 135–140, 2007.
• CARO F., COQUEL F., JAMET D., KOKH S., “A Simple Finite-Volume Method for
Compressible Isothermal Two-Phase Flows Simulation”, Int. J. on Finite Volumes, vol. 3,
num. 1, 2006, p. 1–37.
• HELLUY P., SEGUIN N., “Relaxation model of phase transition flows”, M2AN, Math. Model. Numer.
Anal., vol. 40, num. 2, 2006, p. 331–352.
• LE METAYER O., MASSONI J., SAUREL R., “Elaborating equations of state of a liquid
and its vapor for two-phase flow models”, Int. J. of Th. Sci., vol. 43, num. 3, 2004, p. 265–276.
• G. CHANTEPERDRIX, JP VILA, P. VILLEDIEU, A compressible model for separated two-phase
flow computations, FEDSM02, 14-18 July, Montreal, Quebec, Canada, 2002
44
MULTIMAT 2011 - 5-9 septembre 2011