Transcript Multiphase and Reactive Flow Modelling

Multiphase and Reactive Flow
Modelling
BMEGEÁTMW07
K. G. Szabó
Dept. of Hydraulic and Water Management
Engineering,
Faculty of Civil Engineering
Spring semester, 2012
Basic notions and terminology
Ordinary phases:
preserve volume
Condensed phases
expands
– Solid
– Liquid
– Gaseous
There also exist extraordinary
phases, like plastics and
similarly complex materials
preserves shape
Fluid phases
deform
The property of fluidity serves in
the definition of fluids
Properties and models of solids
Properties of solids:
• Mass (inertia),
position, translation
• Extension (density, volume),
rotation, inertial momentum
• Elastic deformations (small, reversible
and linear), deformation and stress fields
• Inelastic deformations (large, irreversible
and nonlinear), dislocations, failure etc.
Modelled features:
1. Mechanics
•
•
Statics: mechanical equilibrium is necessary
Dynamics: governed by deviation from
mechanical equilibrium
2. Thermodynamics of solids
Mass point model
Rigid body model
The simplest
continuum model
Even more
complex models
Models and properties of fluids
Key properties of fluids:
• Large, irreversible deformations
• Density, pressure, viscosity, thermal conductivity, etc.
(are these properties or states?)
Features to be modelled:
1. Statics
•
•
Hydrostatics: definition of fluid (inhomogeneous [pressure and
density])
Thermostatics: thermal equilibrium (homogenous state)
2. Dynamics
1. Mechanical dynamics: motion governed by deviation from
equilibrium of forces
2. Thermodynamics of fluids:
•
•
Deviation from global thermodynamic equilibrium often governs
processes multiphase, multi-component systems
Local thermodynamic equilibrium is (almost always) maintained
Only continuum models are appropriate!
Modelling Simple Fluids
• Inside the fluid:
– Transport equations
Mass, momentum and energy balances
  

5 PDE’s for p(t, r ), u(t, r ) and T (t, r )
– Constitutive equations
Primary (direct)
field variables
Algebraic equations for  ( p, T ),  ( p, T ), k ( p, T ), 
• Boundary conditions
Secondary (indirect) field variables
On explicitly or implicitly specified surfaces
• Initial conditions
Thermodynamical representations
Representation (independent variables)
TD potential
enthropy and volume (s,1/ρ)
internal energy
temperature and volume (T,1/ρ)
free energy
enthropy and pressure (s,p)
enthalpy
temperature and pressure (T,p)
free enthalpy
• All of these are equivalent:
can be transformed to each other by appropriate formulæ
• Use the one which is most practicable:
e.g., (s,p) in acoustics: s = const  ρ(s,p)  ρ(p).
We prefer (T,p)
Note
Some models of fluids
•
  const ,   const
•
 ( p ),   const
Stoksean fluid
In both of these, the heat transport problem can
be solved separately (one-way coupling):
fluid dynamical
equations
•
compressible
(or barotropic) fluid
heat transport
equation (1 PDE)
 ( p, T ),  ( p, T ), k ( p, T ), 
general simple fluid
Mutually coupled thermo-hydraulic equations:
fluid dynamical
equations
heat transport
equation
• Non-Newtonian behaviour etc.
models for complex fluids
Phase transitions
• Evaporation, incl.
– Boiling
– Cavitation
• Condensation
• Freezing
• Melting
• Solidification
• Sublimation
All phase transitions
involve latent heat
deposition or release
Several solid phases
(crystal structures) may exist
Typical phase diagrams of a pure material:
 1

In equilibrium 1, 2 or 3 phases can exist together
Complete mechanical and thermal equilibrium
Conditions of local phase equilibrium
in a contact point
in case of a pure material
• 2 phases:
T(1)=T(2)=:T
p(1)=p(2)=:p
μ(1)(T,p)= μ(2)(T,p)
Locus of solution:
a line Ts(p) or ps(T),
the saturation
temperature or
pressure (e.g.
‘boiling point´).
• 3 phases:
T(1)=T(2) =T(3)=:T
p(1)=p(2)=p(3) =:p
μ(1)(T,p)= μ(2)(T,p) = μ(3)(T,p)
Locus of solution:
a point (Tt,pt), the triple
point.
Multiple components
Almost all systems have more than 1 chemical components
Phases are typically multi-component mixtures
• Concentration(s): measure(s) of composition
There are lot of practical concentrations in use, e.g.
– Mass fraction (we prefer this!)
c1  m1 m , c2  m2 m ,  ck  mk m , 
 ck   mk
k
m 1
k
– Volume fraction (good only if volume is conserved upon mixing!)
1  V1 V ,  2  V2 V ,   k  Vk V ,   k  Vk V  1
k
k
Concentration fields appear
 as new primary field variables
in the equation: ck ( t ,r ) for k  2,, N p
One of them (usually that of the solvent) is redundant, not used.
Notations to be used
(or at least attempted)
• Phase index (upper):
– (p) or
– (s), (l), (g), (v), (f) for solid, liquid, gas, fluid, vapour
• Component index (lower): k
• Coordinate index (lower): i, j or t
Examples:  ( s) , c( p ) , u( p)
k
i
• Partial differentiation:
 t , i (1   x ,  2   y ,  3   z )
Note
Material properties in
multicomponent mixtures
• One needs constitutional equations for
each phase
• These algebraic equations depend also on
the concentrations
For each phase (p) one needs to know:
– the equation of state
– the viscosity
– the thermal conductivity
– the diffusion coefficients


   p ,T ,c   ,c   ,
k   p ,T ,c   ,c   ,
D   p ,T ,c   ,c   ,
  p  p ,T ,c1 p  ,c2 p  ,
p
k
p
p
1
2
p
1
p
k ,
p
p
2
p
1
p
2
Conditions of local phase equilibrium
in a contact point
in case of multiple components
• Suppose N phases and K components:
• Thermal and mechanical equilibrium on the interfaces:
T(1)=T(2) =T(3)=:T
p(1)=p(2)=p(3)=:p
• Mass balance for each component among all phases:
 
 
















 T , p ,c ,c , ,c   T , p ,c ,c , ,c     T , p ,c 

  , c   , , c       T , p ,c   ,c   , ,c         T , p ,c 

T
,
p
,
c



1
1
1
2
1
1
1
2
1
K
1
1
1
2
1
K
2
1
2
2
2
1
2
2
2
K
2
1
2
2
2
K
N
1
N
1
2
1


N
,c2 N  , ,c K N 
N
,c2 N  , ,c K N
N
,c2 N  , ,c K N 

  , c   , , c       T , p ,c   ,c   , ,c         T , p ,c 

T
,
p
,
c



1
K
1
1
1
2
1
K
2
K
2
1
2
2
2
K
(N-1)K equations for 2+N(K-1) unknowns
N
K
1

Phase equilibrium
in a multi-component mixture
Gibbs’ Rule of Phases, in equilibrium:
# phases N  # components 2  K  2
TD limit on the # of phases
If there is no (global) TD equilibrium:
additional phases may also exist
– in transient metastable state or
– spatially separated, in distant points
Miscibility
The number of phases in a given system is also
influenced by the miscibility of the components:
• Gases always mix →
Typically there is at most 1 contiguous gas phase
•
Liquids maybe miscible or immiscible →
Liquids may separate into more than 1 phases
(e.g. polar water + apolar oil)
1. Surface tension (gas-liquid interface)
2. Interfacial tension (liquid-liquid interface)
(In general: Interfacial tension on fluid-liquid interfaces)
•
Solids typically remain granular
Topology of phases and interfaces
A phase may be
• Contiguous
(more than 1 contiguous
phases can coexist)
Interfaces are
• 2D interface surfaces
separating 2 phases
– gas-liquid: surface
– liquid-liquid: interface
– solid-fluid: wall
• Dispersed:
– solid particles, droplets
• 1D contact lines separating
or bubbles
3 phases and 3 interfaces
– of small size
• 0D contact points with 4
– usually surrounded by
phases, 6 interfaces and 4
a contiguous phase
contact lines
• Compound
Topological limit on the # of phases
(always local)
Special Features to Be Modelled
• Multiple components →
– chemical reactions
– molecular diffusion of constituents
• Multiple phases → inter-phase processes
– momentum transport,
– mass transport and
– energy (heat) transfer
across interfaces.
(Local deviation from total TD equilibrium is typical)
Class 3 outline
•
•
•
•
•
•
Balance equations
Mass balance — equation of continuity
Component balance
Advection
Molecular diffusion
Chemical reactions
•
•
•
•
•
extensive quantity: F
density: φ=F/V=ρ∙f
specific value f=F/m
molar value f=F/n
molecular value F*=F/N
Differential forms of balance
equations
Conservation of F:
• equations for the
density
– general
– only convective flux
• equation for the
specific value
 t  jF   0
 if jF  v  
 t  v     0
 if m is conserved
Dt f   t f  v     0
These forms describe
passive advection of F
Class 4
• Diffusion — continued
!
Incomplete
without class
notes
– further diffusion models
– the advection—diffusion equations
• Chemical reactions
–
–
–
–
–
–
the advection—diffusion—reaction equations
stochiometric equation
reaction heat
chemical equilibrium
reaction kinetics
frozen and fast reactions
Further diffusion models
Thermodiffusion and/or barodiffusion
Occur(s) at
• high concentrations
• high T and/or p gradients
For a binary mixture:
Analogous cross effects
jdiff   D c  kT T T  k p p p 
appear in the heat
D  kT :coefficient of thermodiffusion
D  k p :coefficient of barodiffusion
conduction equation
Further diffusion models
~
~
M k K k  K kk
Nonlinear diffusion model
jk    

 x k
det K 
k M
Cross effect among species’ ~
K  adjK 
diffusion
x
M
x
Valid at
K k  k     s if k  
Dk M k s k Dks
• high concentrations
• more than 2 components K kk  0
• low T and/or p gradients
Dk : binary diffusion coefficient
(For a binary mixture it falls M   xk M k : mean molar mass
k
back to Fick’s law.)
M
xk 
 ck : mole fraction
Mk
The advection–diffusion equations
local rate
of change
advective
flux
diffusive
flux
 t   ck   v    ck   jk 
 since m is conserved
1
Dt ck   t ck  v   ck   jk 

e.g. D  2ck
The concentrations are conserved but not passive quantities
The advection–diffusion–reaction
equations
local rate
of change
advective
flux
diffusive
flux
 t   ck   v    ck   jk 
reactive source terms
 mass production rate density
 since m is conserved
1
Dt ck   t ck  v   ck   jk   local specific production rate

The concentrations are not conserved quantities
Class 5
• Mathematical description of interfaces
–
–
–
–
implicit description
parametric description (homework)
normal, tangent, curvature
interface motion
• Transport through interfaces
–
–
–
–
continuity and jump conditions
mass balance
heat balance
force balance
!
Incomplete
without class
notes
Interfaces and their motion
• Description of interface surfaces:
New
– parametrically
primary(?)
field
– by implicit function
variables
– (the explicit description is the common case of
the two)
• Moving phase interface:
(only!) the normal velocity component
makes sense
!
Incomplete
without class
notes
Description of an interface by an
implicit function
F (t , x, y , z )  0
n  F F (unit normal)
1 1
1 
         n (mean curvature)
2  R1 R2 
Type I surfaceintegrals :
 f (t, x, y, z )  dA   f (t, x, y, z )   F (t, x, y, z )  F (t, x, y, z )  dV
Type II surfaceintegrals :
 v(t, x, y, z )  dA   v(t, x, y, z )   F (t, x, y, z )  F (t, x, y, z )  dV
Equation of motion of an interface
given by implicit function
• Equation of interface
• Path of the point that
remains on the interface
(but not necessarily a
fluid particle)
• Differentiate
• For any such point the
normal velocity
component must be the
same
• Propagation speed and
velocity of the interface
F (t , r )  0
r (t )
F (t , r (t ))  0
d
F (t , r (t ))   t F  r (t )  F  0
dt

 t F  u  n  F  0
u  n  r (t )


u  nu


Only the normal
component makes sense
Parametric description of an
interface and its motion
Homework:
Try to set it up analogously
Mass balance through an interface
Steps of the derivation:
• describe in a reference frame that moves
with the interface (e.g. keep the position of
the origin on the interface)
• velocities inside the phases in the moving
frame
• mass fluxes in the moving frame
• flatten the control volume onto the
!
interface
Incomplete
without class
notes
The kinematical boundary conditions
The net mass flux through the interface:
!

jmass
  1 u1  u  n   2  u2   u  n
def
Incomplete
without class
notes
 u  u  n  0


For the tangential components :
This condition does not follow
u  n  0
t F  u  F  0
from mass conservation
Without (net) mass transfer:

jmass
 0  u 1  n  u 2   n  u  n  u
u  n  0
tangential components :
u  n  0
 u  0
continuity of
velocity
t F  u  F  0
conservation
of interface
Diffusion through an interface
Mass flux of component k in the co-moving reference
frame:
k uk  u   ck  u  u  uk  u  ck  u  u   jk
Case of conservation of component mass:
• on a pure interface



u

u
k
k
  n  0
(no surface phase, no

surfactants)
• without surface
reactions




ck  jmass
 jk  n  0
(not a reaction front)
The component flux
j   k1 uk1  u  n   k2  uk2   u  n
through the interface:


 ck1  jmass
 jk1  n  ck2   jmass
 jk2   n

k
def
Momentum balance through an
interface
Effects due to
• surface tension
• surface viscosity
• surface
compressibility
• mass transfer
Surface tension
• The origin and interpretation of surface
tension
!
Incomplete
without class
notes
Dynamical boundary conditions
with surface/interfacial tension
• Fluids in rest
– normal component:
• Moving fluids without
interfacial mass
transfer
 p  2 S 
Modifies the
thermodynamic phase
equilibrium conditions
 p  n  τ n   2 S 
– normal component:
– tangential components:  t  τ n   t  S t  n 
The viscous stress tensor:  ij     i u j   j ui 
!
Incomplete
without class
notes
The heat conduction equation
The equation
• Fourier’s formula
– (thermodiffusion not included!)
• Volumetric heat sources:
– viscous dissipation
– direct heating
– chemical reaction heat
Boundary conditions
• Thermal equilibrium
• Heat flux:
– continuity (simplest)
– latent heat (phase transition of
pure substance)
Even more complex cases:
– chemical component diffusion
– chemical reactions on surface
– direct heating of surface
c p tT  u  T   jheat  qheat
jheat    T
T   0
n  jheat   0

n  jheat   L  jmass
n  jheat   
Jump
conditions
Approaches of fine models
Phase-by-phase
• Separate sets of
governing equations for
each phase
• Each phase is treated as
a simple fluid
• Describing/capturing
moving interfaces
• Prescribing jump
conditions at the
interfaces
One-fluid
• A single set of governing
equation for all phases
• Complicated
constitutional equations
• Describing/capturing
moving interfaces
• Jumps on the interfaces
are described as singular
source terms in the
governing equations
Phase-by-phase
mathematical models
1.
2.
3.
4.
A separate phase domain for each phase
p    p 
p t , r , u t , r ,
A separate set of balance equations for each
phase domain, for the primary field variables
 p 
T t , r ,
introduced for the single phase problems,
supplemented by the constitutional relations
describing the material properties of the given
 ( p ) T , p,,
phase
The sub-model for the motion of phase domains
 ( p ) T , p,,
and phase boundaries
(further primary model variables)
k ( p ) T , p,,
Prescribing the moving boundary conditions:

coupling among the field variables of the
e.g. F t, r  0
neighbouring phase domains and the interface
F t, r  0
variables
F t, r   0
The one-fluid mathematical model

  p  ( t , r )  1 or 0
A single fluid domain
Characteristic function for each phase   p  ( t, r )  1
p
Material properties expressed by the
properties of individual phases and the
characteristic functions
    ( p)   ( p)
4. A single set of balance equations for the
p
primary field variables introduced for the
single phase problems, supplemented
    ( p)   ( p)
by discrete source terms describing
p
interface processes
( p)
( p)
k



k

5. The sub-model for the motion of phase
p
domains and phase boundaries
(further primary model variables)
  
pt, r , ut, r ,

p 

 ( t, r ) or something
equivalent
T t, r ,
1.
2.
3.
Numerical implementations
of interface sub-models
Main categories
• Grid manipulation
• Front capturing:
implicit interface
representation
• Front-tracking:
parametric interface
representation
• Full Lagrangian
Specific methods
• MAC: (Marker-And-Cell)
• VOF: (Volume-of-Fluid)
• level-set
• phase-field
• CIP
E.g. SPH
Front tracking methods
on a fixed grid
by connected marker points
(Suits the parametric
mathematical description)
• In 3D: triangulated unstructured
grid represents the surface
Tasks to solve:
• Advecting the front
• Interaction with the grid
(efficient data structures are
needed!)
• Merging and splitting (hard!)
MAC
(Marker-And-Cell method)
• An interface reconstruction — front capturing —
model (the primary variable is the characteristic
function of the phase domain, the interface is
reconstructed from this information)
• The naive numerical implementation of the
mathematical transport equation  t     u  0 :
– 1st (later 2nd) order upwind differential scheme
• Errors (characteristic to other methods as well!):
!
– numerical diffusion in the 1st order
– numerical oscillation in higher orders
Incomplete
without class
notes
Due to the
discontinuities
of the function
MAC
!
Incomplete
without class
notes
C nj 1  C nj 
u  t
 C nj  C nj1 
h
VOF
(Volume-Of-Fluid method)
1D version (1st order explicit in time):
• Gives a sharp interface, conserves mass
• Requires special algorithmic handling
The scheme of evolution:
VOF in 2D and 3D
SLIC:
Simple Line
Interface Construction
Hirt & Nichols
PLIC:
Piecewise Linear
Interface Construction
Numerical steps of VOF
1.
Interface reconstruction
within the cell
1.
determine n
•
2.
2.
several schemes
position straight interface
Interface advection
•
several schemes exist,
goals:
•
•
•
conserve mass exactly
avoid diffusion
avoid oscillations
3. Compute the surface
tension force in the
Navier–Stokes eqs.
•
several schemes
Implementation of VOF in
• Any number of phases can be
present
• The transport equation for is
adapted to allow
– variable density of phases
– mass transport between
phases
• Contact angle model at solid
walls is coupled
• Special (`open channel´)
boundary conditions for VOF
• Surface tension is
implemented as a continuous
surface force in the momentum
equation
• For the flux calculations
ANSYS FLUENT can use one
of the following schemes:
– Geometric Reconstruction:
PLIC, adapted to nonstructured grids
– Donor-Acceptor:
Hirt & Nichols, for quadrilateral
or hexahedral grid only
– Compressive Interface
Capturing Scheme for
Arbitrary Meshes (CICSAM):
a general purpose sheme for
sharp jumps (e.g. high ratios
of viscosities) for arbitrary
meshes
– Any of its standard schemes
(probably diffuse and oscillate)
The level set method
[hu: nívófelület-módszer]
F (t , x, y , z )  0
n  F F
 t F  u  n  F  0
1 1
1 
         n
2  R1 R2 
• the interface is implicit
• F is continuous
– standard advection
schemes work fine
• the curvature can be
obtained easily
 f (t, x, y, z )  dA   f (t, x, y, z )   F (t, x, y, z )  F (t, x, y, z )  dV
 v(t, x, y, z )  dA   v(t, x, y, z )   F (t, x, y, z )  F (t, x, y, z )  dV
• the effect of surface
 2 S    F ()  F  dV
tension within a cell
can be computed
The level set method
F (t , x, y , z )  0
n  F
 t F  u  0
F ( t , x , y , z )  1
• If
then the computational
demand can be
substantially decreased
1 1
1 
         n  2 F
2  R1 R2 
 f (t, x, y, z )  dA   f (t, x, y, z )   F (t, x, y, z )  dV
 v(t, x, y, z )  dA   v(t, x, y, z )   F (t, x, y, z )   F (t, x, y, z )  dV
 2 S    F ()  F  dV
Signed distance function
as an implicit level-set function
• What kind of function is it?
F (t , x, y , z )  0 , F (t , x, y , z )  1
Signed distance from the
interface!
• Alas, F  1 is not
 t F  u  F  0
conserved.
• Generating F: τ is pseudo F  S F    F  1  0
time (t is not changed)
• Apply alternatively!
S F   sgn F 
• Unfortunately, mass is not
F
conserved in the numeric


S F 
2
implementation.
F 2  h 2 F
• A better numeric scheme
Numerical implementation of the interfacial
source terms in the transport equations

0
if
 1 F
1
  F 
H  F    

sin
 if
  
 2 2 2

1
if

0

 1
1
  F 
  F     cos

  
 2 2

1
1
t, r 
H  F t , r  



0
  F  
  F 
0
F  
F 
F  
if F  
if F  
if F  
With ε = 1.5h, the interface forces are
smeared out to a three-cell thick band
Only first order accurate in h
• For example, the normal jump condition due to
surface tension can be expressed as an
embedded singular source term in the Navier–
Stokes equation:
 Dt v   g  p    τ  2 S  F n
– contribution to a single cell
in a finite volume model:
C.f. VOF
 2  S  F F dV
cell
• Other source terms (latent heat, mass flux) in
the transport equations can be treated
analogously.
Level set demo simulations
Evaluation criteria for comparison
Not only for VOF and Level Set
• Ability to
– conserve mass/volume
exactly
– numerical stability
– keep interfaces sharp
(avoid numerical diffusion
and oscillation)
• Ability and complexity to
model
– more than 2 phases
– phase transitions
– compressible fluid phases
• Demands on resources
–
–
–
–
–
number of equations
grid spacing
grid structure
time stepping
differentiation schemes
• Limitations of applicability
– grid types
– differential schemes
– accuracy
Recommended books
• Stanley Osher, Ronald Fedkiw:
Level Set Methods and Dynamic Implicit
Surfaces
Applied Mathematical Sciences, Vol. 153
(Springer, 2003). ISBN 978-0-387-95482-0
– Details on the level set method
• Grétar Tryggvason, Ruben Scardovelli,
Stéphane Zaleski:
Direct Numerical Simulations of Gas–
Liquid Multiphase Flows
(Cambridge, 2011). ISBN 9780521782401
– Modern solutions in VOF and front tracking
SPH
Smoothed Particle Hydrodynamics
• The other extreme — a meshless method:
The fluid is entirely modelled by moving
representative fluid particles — fully Lagrangian
• There are no
–
–
–
–
mesh cells
interfaces
PDE
field variables
• Everything is described via ODE’s
SPH simulation of hydraulic jump
Fr1 = 1.15
Fr1 = 1.37
Fr1 = 1.88
SPH simulation of dam-break
Void bubble
Vacuum
Entrapped air
Air
Liquid vs. liquid-gas simulation
Evaluation of SPH
Advantages
• Conceptually easy
• Best suits problems
– in which inertia dominates
(violent motion, transients,
impacts)
• FSI modelling
– with free surface or liquid–
gas interface
• Interface develops
naturally
• Computationally fast
– Easy to parallelise
– Can be adapted to GPU’s
Disadvantages
• High number of particles
• Hard to achieve
incompressibility
• Some important boundary
conditions are not
realised so far