SVG-Viewer needed.

SVG-Viewer needed.

SVG-Viewer needed.

Equations & Algorithms


Contents

Contents
1 Incompressible Navier–Stokes Equations
 1.1 Single Phase Conservation Equations
 1.2 The Scalar Equation
 1.3 Pressure Solver Implicit Formulation
  1.3.1 Divergence-Free Flows
  1.3.2 Discretisation
 1.4 Axisymmetric Flows
  1.4.1 Axisymmetry
2 Compressible Flow Modelling
 2.1 Single Phase Conservation Equations
 2.2 N-Phase Conservation Equations
 2.3 Level-Set Conservation Equations
 2.4 Equations of state
  2.4.1 Ideal Gas equation of state
  2.4.2 Stiffened Gas Equation of State
  2.4.3 Tait Equation of State for Water
  2.4.4 Tait Equation of State for Ammonia
  2.4.5 Incompressible
  2.4.6 Peng-Robinson Equation of State for Water (Liquid)
  2.4.7 Peng-Robinson Equation of State for Water (Vapour)
  2.4.8 Sutherland’s law
 2.5 Pressure Solver Formulation
  2.5.1 Mass conservation
  2.5.2 Pressure equation
  2.5.3 Limit of incompressible multiphase flow
  2.5.4 Limit of compressible single phase flow
  2.5.5 Pressure Solver Formulation for reactive flows
 2.6 Numerical Schemes
  2.6.1 Interpolation schemes for compressible density
   Upwind scheme
   Linear scheme
  2.6.2 Unsteady simulation
3 Numerical Methods
 3.1 Governing Equations for Steady Flows
 3.2 Discretisation of Governing Equations
  3.2.1 Flux balance equations
  3.2.2 Approximation of convection terms
  3.2.3 Approximation of diffusion terms
  3.2.4 Approximation of source terms
  3.2.5 Final form of discretised equations
 3.3 Boundary Conditions
  3.3.1 Inflow planes
  3.3.2 Outflow planes
  3.3.3 Symmetry planes
  3.3.4 Rigid wall
  3.3.5 Total pressure boundary condition
   Ideal gas
   Incompressible or weakly compressible fluids
   Mixture of compressible fluids
   Linearization for pressure coefficient and flux update
 3.4 Solution Procedure
  3.4.1 Convergence criterion
  3.4.2 Calculation sequence
4 Reactive Flow Modelling
 4.1 Conservation equations
  4.1.1 Species equation
  4.1.2 Temperature equation
 4.2 Turbulent reacting flow
  4.2.1 RANS
  4.2.2 LES
 4.3 Non–Premixed Reactive Flows
  4.3.1 Slow Chemistry
  4.3.2 Infinitely Fast Chemistry
5 Boussinesq Approximation
6 Non-Newtonian Models
7 Viscoelasticity
 7.1 Olroyd-B. model
 7.2 Exponential Phan-Thien-Tanner model (PTT)
 7.3 Linear PTT model
 7.4 Temperature correction using Arrhenius law
 7.5 Multiphase viscoelasticity
 7.6 Logarithm of Conformation Tensor (LCT)
8 Radiation
9 Immersed Surfaces Technique (IST)
 9.1 Introduction
 9.2 Equations and Implementation
  9.2.1 Immersed Level Set and solid Heaviside functions
  9.2.2 Solid temperature equation
10 Smart Adaptive Run Parametrisation (SArP) algorithm
 10.1 Theory
 10.2 Genetic Algorithm
  10.2.1 Structure
  10.2.2 Triggers and stopping criteria
   Triggers
   Stopping criteria
  10.2.3 Genetic Operations
   Mutation Operation
   Crossover Operation
  10.2.4 Selection Process
   Fitness function
   Selection Operators
  10.2.5 Prediction and Fitness approximation
   Strategy
   Learning Methods
  10.2.6 Speeding up the EGA
References

Chapter 1
Incompressible Navier–Stokes Equations

1.1 Single Phase Conservation Equations

Considering the motion of a Newtonian fluid with constant properties (density and viscosity) velocity ui = (u,v,w), pressure p and temperature T are obtained by solving the incompressible Navier-Stokes equations which, in Cartesian–coordinates, are respectively:

                 ∂uj
                 ---- =   0                                        (1.1)
                 ∂xj
∂-(ρu ) + -∂--(ρuu )  =   - ∂p-+  ∂(σij)-+ ρg ;  i = 1,2,3          (1.2)
∂t   i    ∂xj    i j        ∂xi    ∂xj      i
      ∂T         ∂T        ∂  (  ∂T )    ′′′
  ρCp ---+ ρCpuj ---- =   ---- λ ---- + ˙Q                          (1.3)
      ∂t         ∂xj      ∂xj ◟--∂◝x◜j-◞
                                - q′j′
where ρ represents the density, gi is the force of gravity, μ is the molecular viscosity, σij is the viscous stress tensor and Q˙′′′ is a volumetric heat source. ρ and μ are considered as constant for these incompressible equations. Note that the Einstein summation convention applies to repeated indices.

There is a need for closure laws for the viscous stress tensor σij. This will in turn express the way a fluid environment is deformed (or strained) under a given stress and thus these constitutive equations are also known as the stress–strain relations. Making use of the Stokes’ hypothesis yields:

       (          )
         ∂ui-  ∂uj-    2-   ∂uk-
σij = μ  ∂xj + ∂xi   - 3μ δij∂xk                          (1.4)

Fourier law is also used, which assumes that the thermal conductivity λ is a positive constant, which means that isotropic heat conductivity is assumed.

In some cases, solids or fluids do not have isotropic properties for heat conduction. In this case, thermal conductivity λ can differ according to the conduction direction. The energy equation 1.3 becomes:

   ∂T          ∂T      ∂qj    ′′′
ρCp-∂t + ρCpuj ∂x--= - ∂x--+ Q˙                          (1.5)
                 j       j

where

qj = - λij ∂T
         ∂xi
(1.6)

takes into account the anisotropy of heat conduction. λij is a tensor, based on the thermal conductivity along the principle axes.

1.2 The Scalar Equation

The evolution equation for the scalar C has the form:

                        (      )
∂C     ∂             ∂      ∂C
∂t-+  ∂x-(ujC ) =   ∂x--  D∂x--  + FC                     (1.7)
        j              j      j
where D = ν∕Sc is the molecular diffusivity of the scalar, and ρ is the density of the fluid. ν is the kinematic viscosity, and Sc is the Schmidt number. The source term FC represents any source term applied on the scalar, depending on what C is actually representing.

1.3 Pressure Solver Implicit Formulation

We start with the discretised version of the momentum equations given by,

           ∑                      k+1
AP,ijuk+1 =    Anb,ijuk+1+ Su - ∂p-P--ΩP
     P,j    nb       nb,j         ∂xi
(1.8)

where Aij is the coefficient (diagonal) tensor with the diagonals being, Au, Av, and Aw, subscript P denotes the point being discretised, and nb denotes its neighbours. ΩP is the cell volume. We split the unknown pressure (pk+1) into two parts, pk+1 = pk + p and introduce an intermediate velocity uik* to obtain:

                         ∑                      k        ′
AP,ij(ukP*,j + uk+1 - ukP*,j) =   Anb,ijuk+1+ Su -  ∂pPΩP  - ∂pP-ΩP
            P,j           nb       nb,j        ∂xi      ∂xi
(1.9)

Applying operator splitting we obtain two equations,

              k*     ∑         k*         ∂pkP-
        AP,iju P,j  =      Anb,ijunb,j + Su - ∂xi ΩP                (1.10)
                      nb
      k+1    k*      ∑          k+1    k*    ∂p-′P
AP,ij(uP,j - uP,j)  =      Anb,ij(unb,j - unb,j)- ∂xi ΩP             (1.11)
                      nb
If the first term in the RHS is approximated as,
∑    Anb,ij(uk+1- uk*  )
---nb--∑----nb,j----nb,j = (ukP+,1j - ukP*,j)
         nbAnb,ij
(1.12)

based on the idea that the cell–centre velocity correction is a coefficient weighted average of the velocity corrections at the neighbouring cells, we obtain the SIMPLEC method. If this term is neglected then we obtain the SIMPLE method.

       ∑                         ∂p ′
(AP,ij -    Anb,ij)(uk+P1,j - ukP*,j) = ----PΩP
        nb                       ∂xi
(1.13)

Denoting AP,ij = (AP,ij - nbAnb,ij) and multiplying Eq. (1.13) by AP,ij′-1 and taking the divergence, we obtain the pressure–correction equation as,

∂uk+1   ∂uk*       ∂  (     ∂p′   )
--P,j-- ---P,j = - ---- A ′-P1,ij--P-ΩP
 ∂xj     ∂xj      ∂xj       ∂xi
(1.14)

1.3.1 Divergence-Free Flows

For incompressible single–phase flows, two–phase flows using interface tracking, or for homogeneous flow models, we set the expectation on the divergence of uk+1 to be zero, thus obtaining the pressure correction equation as,

∂ukP*,j    ∂  (     ∂p′   )
-----=  ---- A ′-P1,ij--P-ΩP
∂xj     ∂xj       ∂xi
(1.15)

1.3.2 Discretisation

The integral form of Eq. (1.15) is solved. Integrating using the divergence theorem we get the following discretised equation,

                          ′
∑  ˜uk*n   = ∑   n  A -1 ∂pPΩ
    f,j f,j       f,j  f,ij∂xi  f
 f           f
(1.16)

where, the velocity at the faces ũ is calculated with the Rhie–Chow method, nf,j is the area vector of the cell face f, Af,ij-1 is the interpolated coefficient tensor at the face and Ωf is the volume of dual control volume centered at the face.

1.4 Axisymmetric Flows

1.4.1 Axisymmetry

For axisymmetric geometries, the continuity equation is given by

∂ρ-+ -∂-(ρu)+  ∂-(ρv)+  ρv-= 0
∂t   ∂x        ∂r       r
(1.17)

or

∂ρ-   ∂--      1-∂-
∂t +  ∂x(ρu) + r∂r (rρv) = 0
(1.18)

where x is the axial coordinate, r is the radial coordinate, u is the axial velocity, and v is the radial velocity.

The axial and radial momentum conservation equations are given by

                                          [  (               )]
∂-(ρu) + ∂--(ρuu )+ 1-∂-(rρvu ) = - ∂p-+-∂- μ   2∂u-- 2-(∇ ⋅⃗v)
∂t       ∂x        r ∂r           ∂x   ∂x       ∂x   3
                                         1 ∂ [   ( ∂u   ∂v )]
                                       + r∂r- rμ   ∂r-+ ∂x-   +         (1.19)

                                                                        (1.20)
and
 ∂        ∂        1 ∂            ∂p    ∂ [  ( ∂v   ∂u )]
--(ρv) + ---(ρuv )+ ----(rρvv) = - ---+ --- μ   ---+ ---
∂t       ∂x     [  r(∂r           ∂r)] ∂x      ∂x   ∂r
          + 1-∂--rμ   2∂v-- 2-(∇ ⋅⃗v)   - 2μ-v + 2-μ-(∇ ⋅⃗v)           (1.21)
            r ∂r       ∂r   3              r2   3 r
where
∇ ⋅⃗v = ∂u-+  ∂v-+ v-
       ∂x    ∂r   r
(1.22)

or

       ∂u-  1-∂--
∇ ⋅⃗v = ∂x + r ∂r(rv)
(1.23)

Chapter 2
Compressible Flow Modelling

This section describes the compressible flow modelling capabilities of TransAT.

2.1 Single Phase Conservation Equations

Equations that are solved for compressible flows are the following (Poinsot & Veynante2001):

     -∂ ρ+  -∂-(ρuj)  =  0                                             (2.1)
     ∂t     ∂xj
∂-        -∂--             -∂p-  ∂-σij
∂t (ρui)+ ∂xj (ρuiuj) =  - ∂xi +  ∂xj + ρgi                            (2.2)
                              (     )   (            )
  ρCp∂T- + ρCpuj ∂T-- =   -∂-- λ ∂T-- +   ∂p-+ uj-∂p-  + Φ + Q˙′′′       (2.3)
      ∂t         ∂xj      ∂xj    ∂xj      ∂t     ∂xj
where
       (          )
σ  = μ   ∂ui+  ∂uj-  - 2μδ  ∂uk-
 ij      ∂xj   ∂xi     3  ij∂xk
(2.4)

and Φ is the viscous dissipation of kinetic energy. Thermodynamical properties ρ, p and T are linked using a specified equation of state

ρ = ρ(p,T )
(2.5)

2.2 N-Phase Conservation Equations

Equations that are solved for compressible n-phase flows are given belows. In the current version, compressible flow is only available for the homogeneous model (no drift velocity), without phase change.

        ∂  m    ∂  ( m  m)
        ∂tρ  + ∂x-- ρ  uj   =   0                                    (2.6)
                  j
     ∂ρkαk-+  -∂-(ρkαkum )  =   0                                    (2.7)
       ∂t     ∂xj       j
 ∂            ∂ (        )        ∂p    ∂σmij
-- (ρmumi )+  ----ρmumi umj   =   - ---+  ----+  ρgi                   (2.8)
∂t          (∂xj         )        ∂xi(   ∂xj)
     ρmC  m   ∂T-+ u -∂T-   =   -∂-- λ ∂T-- +  Dp-+ Φ + Q˙′′′         (2.9)
         p    ∂t    j∂xj        ∂xj    ∂xj     Dt
                     N∑
                        α   =   1                                   (2.10)
                     k=1 k
with the subscript k denoting the phase index and the superscript m denoting a mixture variable, and α being the volume fraction. Mixture variables are defined as follows
                                      ∑N
             Mixture density      ρm =    αk ρk                  (2.11)
                                      k=1
                                      ∑N αkρkuk
            Mixture velocity     um =  k=1ρm----                 (2.12)
                                       N
       Mixture heat capacity  ρmCm   = ∑  αkρkCp,k               (2.13)
                                 p    k=1
                                      ∑N
            Mixture viscosity      μm =    αk μk                  (2.14)
                                      k=1
                                  m   ∑N
Mixture thermal conductivity      λ  =    αk λk                  (2.15)
                                      k=1

2.3 Level-Set Conservation Equations

The equations used to model compressible flows with Level-Set method are the following. In the current version, phase change is not available for compressible flows.

       ∂-     I-∂--
       ∂tϕ + uj∂xj ϕ  =  0                                       (2.16)

     -∂ρ + -∂--(ρuj)  =  0                                       (2.17)
     ∂t    ∂xj
∂-       -∂--              -∂p-  ∂σij           I
∂t (ρui)+ ∂xj (ρuiuj) =  - ∂xi + ∂xj  + ρgi + γκδ (ϕ)ni          (2.18)
     (            )          (      )
 ρCp   ∂T-+  uj ∂T-   =  --∂-  λ-∂T-  + Dp- + Φ + ˙Q′′′            (2.19)
       ∂t      ∂xj       ∂xj    ∂xj     Dt
where γ is the surface tension coefficient, κ is the interface curvature, n is the normal vector to the interface, δI is a smoothed Dirac delta function centered at the interface, T is the temperature, Cp is the heat capacity of the fluids, λ is the heat conductivity, and Q˙ is the volumetric heat source. Material properties are updated using χ:

            ∑    k k
 ρ(χ,t) =      χ  ρ                              (2.20)
 μ(χ,t) =   ∑  χk μk                             (2.21)
            ∑
               χkρkCkp
Cp(χ,t) =   --∑-χkρk--                           (2.22)
            ∑
 λ(χ,t) =      χk λk                             (2.23)

2.4 Equations of state

2.4.1 Ideal Gas equation of state

The simplest way to link density, pressure and temperature is to use the ideal gas equation of state. This is applicable only for gases. It reads

    MP
ρ = RT---
(2.24)

where M is the molecular weight and R the universal gas constant. The different parameters must respect the equality

γ = ---Cp---
         R--
    Cp - M
(2.25)

where R = 8.314472J∕Kmol is the universal gas constant. The speed of sound is then written

    ∘ -----
c =   γRT--
       M
(2.26)

2.4.2 Stiffened Gas Equation of State

The Stiffened Gas Equation is an equation of state used for water under very high pressures. This equation can be applied for certain other liquids, but specific heat capacity is likely to be miscalculated. It reads

ρ = -γ(p-+-p∞-)-
    Cp (γ - 1) T
(2.27)

The fluid will behave as an ideal gas at pressure p = p.

The different parameters must respect the equality

γ = ---Cp---
         R--
    Cp - M
(2.28)

where R = 8.314472J∕Kmol is the universal gas constant. The speed of sound can be written

    ∘ ------------
c =   Cp (γ - 1) T
(2.29)

2.4.3 Tait Equation of State for Water

The Tait equation is an equation of state derived for water and sea-water Cooke & Chen (1992). Density of water as a function of pressure can then be expressed as

      ( p + B ) 1N-
ρ = ρ0  -------
        p0 + B
(2.30)

where ρ0 = 1000.0kgm3, p0 = 0.1MPa, B = 295.9MPa, and N = 7.415. For this equation of state, speed of sound can be written

    ┌│ -----------1--
    │∘ N- (p0 +-B-)N-
c =   ρ0 (p + B )1-NN-
(2.31)

2.4.4 Tait Equation of State for Ammonia

The Tait equation of state is also available in TransAT for ammonia. The density is expressed as:

ρ = --------ρ0(------)
    1 - C.log  p+B--
               p0+B0
(2.32)

where:
ρ0 = 683.44kgm3,
p0 = 257.79kPa,
C(T) = 1.45(          - 13.0628+4.4657log10(T))
 0.1180 - 10,
B(T) = max(1,8241.8 - 38.241T + 0.043507T2)atm,
B0 = B(T0),
T0 = 238.7K

2.4.5 Incompressible

TransAT also offers the possibility to use an incompressible equation of state. In that case the density and other properties do not vary with pressure and the constant values specified by the user will be used.

2.4.6 Peng-Robinson Equation of State for Water (Liquid)

Another possibility when working with water is to use the Peng-Robinson equation of state.
In general the Peng-Robinson equation of state writes:

p = -RT----- ------aα-------
    Vm - b   Vm2+ 2bVm  - b2
(2.33)

Started from this equation a cubic equation for the molar specific volume is obtained and solved. The density is then obtained using following relation:

      M
ρ = ----w---
    Vm - Vt
(2.34)

with V t = 2.58775.10-6m3mol the volume translation and Mw = 18gmol the molecular weight

2.4.7 Peng-Robinson Equation of State for Water (Vapour)

The Peng-Robinson equation of state is also available for water vapour. The implementation is similar to the liquid state implementation.

2.4.8 Sutherland’s law

The viscosity μ and the thermal conductivity λ can be evaluated using Sutherland’s law. They are then defined by the equations:

       (    )3 ∕2
μ = μ0  -T--     Tμ,0-+-Sμ-                           (2.35)
        T μ,0      T + Sμ
       (  T )3 ∕2 T   + S
λ = λ0   ----    -λ,0---λ-                           (2.36)
         Tλ,0      T + Sλ
where T is the local temperature. These expressions can also be used for incompressible flows.

2.5 Pressure Solver Formulation

2.5.1 Mass conservation

The phasic continuity equation is written below in the volume fraction and mass fraction forms as,

∂           ∂         m
∂t (αkρk) + ∂x-(αkρkuj ) = 0
             j
(2.37)

 ∂           ∂
-- (ρmYk) + ----(ρmYkujm ) = 0
∂t          ∂xj
(2.38)

where superscript m stands for mixture. It can be verified for both Eqs. (2.37) and (2.38) that summation over the phase index k gives the mixture continuity equation,

∂ ρm    ∂
---- + ----(ρmumj ) = 0
 ∂t    ∂xj
(2.39)

Equation (2.38) can also be written in a material derivative form as follows, which will be useful in the derivation of the pressure equation,

   [             ]     [   m      m  m]
ρm  ∂Yk- + um ∂Yk- + Y   ∂ρ--+  ∂ρ-u--- = 0
     ∂t     j ∂xj     k   ∂t     ∂xj
(2.40)

and plugging Eq. 2.39 gives

   [             ]
 m  ∂Yk-   m ∂Yk-
ρ    ∂t + uj ∂xj  =  0
(2.41)

We introduce now the notation of material derivative advected by the mixture velocity as,

Dm ϕ    ∂ϕ    m ∂ ϕ
-Dt--=  ∂t-+ uj ∂x--
                   j
(2.42)

where ϕ is any advected quantity and the superscript m denotes that the advection velocity is the mixture velocity. Using this notation we can write the material derivative of Y k as,

    m
ρmD---Yk = 0
    Dt
(2.43)

2.5.2 Pressure equation

An equation for the evolution of the pressure of a multiphase system is derived here. It is based on a single pressure, single temperature assumption where the pressure p and temperature T are the mixture pressure and temperature, respectively. The starting point to derive a pressure equation for compressible multiphase flow is the mixture continuity equation (Eq. (2.39)), but written in a non-conservative form as,

∂-ρm    m ∂ρm-    m∂umj
 ∂t  + uj ∂xj +  ρ  ∂xj = 0
(2.44)

This can be rewritten in the material derivative notation as,

Dm ρm       m∂umj
--Dt-- = - ρ -∂x-
                j
(2.45)

Taking ρm = ρm(p,T,Y k) we can expand the above equation as a linearization based on the independent unknowns as,

                                                          m
∂ρm-    Dmp--   ∂ρm-    DmT--  ∑   ∂ρm-   DmYk--      m∂u-j
 ∂p |T,Yk Dt  +  ∂T  |p,Yk Dt  +     ∂Yk |p,T  Dt   = - ρ  ∂xj
                                k
(2.46)

To evaluate the constrained derivatives, we start with the definition of density with respect to phasic densities and mass fractions.

 1    ∑  Yk
-m-=     ---
ρ      k ρk
(2.47)

where ρk = ρk(p,T). Taking the derivative of ρm with respect to p at constant (T,Y k) gives,

∂ρm-         m 2∑   Yk-∂ρk-     m ∑  αk-∂ρk-
 ∂p |T,Yk = (ρ )     ρ2k ∂p |T = ρ      ρk ∂p |T
                  k                k
(2.48)

Similarly the derivative at constant (p,Y k) can be written as,

∂ρm             ∑  Y  ∂ρ         ∑  α  ∂ρ
----|p,Y = (ρm )2    -2k--k|p = ρm    --k--k-|p
∂T    k          k ρk ∂T          k ρk ∂T
(2.49)

The derivative of density with respect to Y k at constant (p,T) is obtained as,

   m          m 2
∂ρ--|   = - (ρ-)--
∂Yk  p,T       ρk
(2.50)

Additionally we write the mixture temperature equation 2.9 using material derivatives

       DmT     Dmp     ∂  (  ∂T )         ′′′
ρmCpm  -----=  -----+ ---- λ ---- + Φ + Q
        Dt      Dt    ∂xj    ∂xj
(2.51)

where Φ is the viscous dissipation of kinetic energy. The material derivative of T advected by the mixture velocity can then be written as,

DmT        1    (Dmp      )
-Dt-- = ρmC--m-  -Dt--+ Ξ
            p
(2.52)

where Ξ represents the entropic part of heat transfer (the last three terms in the RHS of Eq. (2.51)).

Plugging Eqs. (2.48)-(2.50), Eq. (2.52) and Eq. (2.43) into Eq. (2.46), we obtain

[         (                     ) ]       [               ]                 m
 ∑   m αk-  ∂ρk-    ---1---∂ρk-    Dmp--    ∑   m αk-∂ρk-  ---Ξ---      m∂u-j
    ρ  ρk   ∂p |T + ρmCpm  ∂T  |p    Dt  +      ρ  ρk ∂T |p ρmCpm   = - ρ  ∂xj
  k                                          k
(2.53)

where the identity ρmY k = ρkαk is used. This is the final form of the pressure equation for n-phase compressible flow.

2.5.3 Limit of incompressible multiphase flow

In the limit of incompressible flow, the LHS of Eq. (2.53) can be set to zero, giving,

  m
∂uj- = 0
∂xj
(2.54)

2.5.4 Limit of compressible single phase flow

For the limit of compressible single-phase flow, we set ρk = ρ, αk = 1 in Eq. (2.53) to obtain,

(                 )
  ∂ρ-    -1--∂ρ-    Dmp--  ∂ρ-  -Ξ--     ∂umj-
  ∂p|T + ρCp ∂T |p   Dt  + ∂T |pρCp = - ρ∂xj
(2.55)

Note that for an ideal gas p = ρRT, it can be shown that,

∂ρ      1  ∂ρ     ∂ρ
--|T + -------|p = ---|s
∂p     ρCp ∂T     ∂p
(2.56)

where s denotes entropy. Substituting the partial derivatives, we get

∂ρ      1     1  ( - p )    1  (     R )     1
--|s = ---+  ----  ---2- = ----  1- ---  =  -----
∂p     RT    ρCp   RT      RT       Cp      γRT
(2.57)

where we have used the relationships Cp - Cv = R and Cp∕Cv = γ.

2.5.5 Pressure Solver Formulation for reactive flows

A conservative formulation is used for the pressure solver in the case of reactive, compressible flows. We start with Eq. (1.13) and multiply it through by ρPk+1 to obtain:

 k+1  k+1    k*      k+1  ′- 1∂p′
ρP  (uP,j - u P,j) = - ρP A P,ij∂xP-ΩP
                                i
(2.58)

Taking the divergence of the above equation we obtain:

                                     (                )
-∂-- k+1 k+1    -∂-- k+1  k*     --∂-   k+1  ′-1 ∂p′P-
∂xj(ρP  uP,j ) - ∂xj (ρP  uP,j) = -∂xj   ρP  A P,ij∂xi ΩP
(2.59)

Splitting ρPk+1 = ρPk* + ρP and setting,

 ∂                ∂ρk+1
---(ρkP+1uk+P1,j ) = ---P---
∂xj                 ∂t
(2.60)

we obtain,

                                                 (                      )
  ∂ρkP*-  ∂ρ′P-   -∂-- k* k*    -∂-- ′  k*      -∂--   k*   ′   ′- 1∂p′P-
-  ∂t  -  ∂t -  ∂xj(ρP uP,j) - ∂xj (ρP uP,j) = - ∂xj  (ρP + ρP )A P,ij∂xi ΩP
(2.61)

Neglecting the ρp non–linear terms in the RHS and substituting ρP = pP∕cP2, where cP is the sound speed, we obtain,

    k*                     ′        ( uk*   )        (          ′   )
- ∂-ρP-- -∂-(ρk*ukP*,j) = ∂-(pP-)+ --∂-  -P,jp′  -  -∂-- ρk*A ′-1 ∂pPΩP
   ∂t    ∂xj  P        ∂t c2P    ∂xj    c2P  P     ∂xj   P   P,ij∂xi
(2.62)

2.6 Numerical Schemes

The solver for compressible flows is based on the work of Karki & Patankar (1989) andMichelassi et al. (1994). The starting point is the incompressible solver presented in Chapter Numerical Methods, developed by Zhu (1992). The mass conservation equation for steady-state conditions reads

    (       )(       )
∑    ρ*+ ρ ′  c* + c′  .n  = 0
 n    n    n   n    n   n
(2.63)

where ρ is the density and cn.nn is the volume flux through the face n. Superscript * denotes guessed properties while denotes corrections to ensure continuity. This equation can be rewritten as follows

∑   * *     ∑    * ′     ∑   ′ *      ∑   ′ ′
   ρncn.nn+     ρncn.nn +    ρncn.nn +    ρncn.nn = 0
◟n--◝◜---◞  ◟n---◝◜----◞  ◟n--◝◜---◞   ◟n--◝◜---◞
    T0           T1          T2           T3
(2.64)

Where T0 and T1 are the incompressible contribution, T2 is the compressible contribution and T3 is neglected. ρ is linked to p using the equation of state relation

 ′  ( ∂ρ )   ′
ρ =   ∂p-   p
          Φ
(2.65)

where Φ is the temperature or entropy, whether we use the hypothesis of isothermal or entropic transformation.

Equation 3.13 is then rewritten for compressible flows

  C ′   ∑    C ′   ∑   * *
A p pp -   A npn =    ρncn.nn
         n          n
(2.66)

where

      ∑
ACp =    ACn - SCp                                (2.67)
       n
   C     I    corr
  An =  An + An                                  (2.68)
where AnI are the incompressible coefficients given in section 3.2, and Ancorr and SpC are correction terms due to compressibility effects. They depend on the way of evaluating ρat the cell face.

2.6.1 Interpolation schemes for compressible density

In this section, we note P the current cell, W and E the neighbouring cells and w and e the west and east faces of cell p.

Upwind scheme

For high velocities, an upwind scheme can be used. We can thus write

ρ′w = f+w ρ′W + (1 - f+w)ρ′P                             (2.69)
 ρ′ = f+ρ′ + (1 - f+)ρ′                             (2.70)
  e    e P         e  E
with
        +   1-+-sgn(cn.nn)
      f n =       2                                  (2.71)
         {  1   if  ϕ > 0
sgn(ϕ) =   - 1  if  ϕ < 0                            (2.72)
With this scheme, we obtain
                        [(∂ ρ)  ]
Acowrr  =  max (cw.nw,0)   ---                                                (2.73)
                        [(∂p ) Φ]W
  corr                      ∂ρ-
A e    =  - min (ce.ne,0)   ∂p                                                (2.74)
                               Φ(E                                          )
   C      ∑    corr   [( ∂ρ)  ]     ∑                      ∑
  SP   =     A n   +    ---     (        min (cn.nn,0) -        max (cn.nn,0))(2.75)
           n            ∂p  Φ P   n=W,S,B                n=E,N,T
Linear scheme

This scheme can be used for compressible flows at low velocities. We have

ρ′w =  fwρ′W + (1 - fw)ρ′P                             (2.76)
  ρ′e = feρ′P + (1- fe)ρ′E                             (2.77)
with fw and fe the linear interpolation factors for the W - P and P - E intervals, respectively. With this scheme, we obtain
                     (   [( ∂ρ )  ]            [( ∂ρ)  ]  )
Acowrr  =  fw (cw.nw,)  fw   ---      + (1- fw )   ---                       (2.78)
                           (∂p[ Φ(  W)  ]          ∂p[( Φ P)  ]  )
  corr                            ∂ρ-                  ∂ρ-
A e    =  - (1- fe)(ce.ne, ) fe   ∂p  Φ  P + (1 - fe)   ∂p  Φ E               (2.79)
          ∑                          (   [(    ) ]             [(   )  ] )
  SCP   =      Acnorr + (1- fw )(cw.nw, )  fw   ∂ρ-      + (1- fw )   ∂ρ-       (2.80)
           n                                ∂p  Φ W               ∂p  Φ P
                           (   [( ∂ρ)  ]            [(∂ ρ)  ] )
              - (fe)(ce.ne,)  fe   ---     + (1-  fe)   ---                   (2.81)
                                  ∂p  Φ P             ∂p   Φ E

2.6.2 Unsteady simulation

When performing an unsteady simulation, an extra-term must be added to the compressible contribution SP C. The latter becomes:

                         [(    ) ]
SCP,unsteady = SCP,steady - -Ω-   ∂ρ-
                      Δt    ∂p  Φ  P
(2.82)

where Ω is the cell volume and Δt the time step.

Chapter 3
Numerical Methods

3.1 Governing Equations for Steady Flows

The partial differential equations governing steady incompressible flows in nonorthogonal coordinates may be written in the following general form (Patankar1980):

-∂--(Ciϕ- Di ϕ) = S ϕ,   i = 1,2,3
∂xi
(3.1)

where the coefficients Ci and D are related to convection and diffusion, respectively and Sϕ the source terms for any primary unknown ϕ. Eq. 3.1 uses the velocity components V i, which are along the coordinates yi instead of along the grid-aligned directions xi.

Eq. 3.1 has two important features which are worthy of mention: (i) it is in conservation form which implies all terms arising from the divergence operator are under differential operators, and (ii) it does not contain second-derivatives of coordinates (curvature terms) which are very sensitive to grid smoothness.

3.2 Discretisation of Governing Equations

3.2.1 Flux balance equations

Integrating Eq. 3.1 over a typical control volume centered at P (Figure 3.1) leads to a flux balance equation

                          ∫
Ie - Iw + In - Is + It - Ib =  SϕdV
                           ΔV
(3.2)

where If represents the total flux of ϕ across the cell-face f (= e, w, n, s, t or b). Each of the surface fluxes If contains a convective contribution IfC and a diffusive contribution IfD, that is

If = ICf + IDf
(3.3)

Eq. 3.2 involves no approximation and represents a finite-volume analogue of the differential equation ( 3.1).

3.2.2 Approximation of convection terms

The convective contribution in Eq. 3.3 can be written as:

 C
If = Cfϕf
(3.4)

where Cf is the mass flux across the cell-face f, and can be calculated, for w-, s- and b-faces, as:

C  = (ρV )
 w       jw

Cs = (ρVj)s

Cb = (ρVj)b
(3.5)

The determination of ϕf is a key element for both accuracy and stability of numerical solutions. The more accurate schemes tend to be less stable, and vice versa. Four schemes are incorporated in the program, leaving a choice for the user. They are the standard hybrid differencing (Spalding1972) which is basically of first-order accuracy for convection-dominant flows, third-order unbounded QUICK (Leonard1979), and second-order bounded HLPA (Zhu1991). Consider the w-face of the control volume shown in Figure 3.1 and assume, without loss of generality, that V w 0 (V w is the normal velocity at the w-face).


PIC

Figure 3.1: Nodes required by convection schemes in x1–direction


The HLPA scheme can be written as

                       ϕ   - ϕ
ϕw = ϕW  + γw(ϕP - ϕW )-W-----W-W-
                       ϕP - ϕW W
(3.6)

where

     {
       1          if|ϕˆW  - 0.5| < 0.5
γw =   0          otherwise
(3.7)

All the schemes but Hybrid are implemented in the program in a deferred way proposed by Khosla & Rubin (1974)

ϕl+f1 = ϕuf,l+1 + λ(ϕhf,l- ϕuf,l)
(3.8)

where u and h indicate the upwind and higher-order schemes, and l represents the iteration level; the parameter λ blends the two schemes with limiting values λ = 0 for the upwind and λ = 1 for the higher-order scheme. The differed correction leads to an always diagonally dominant coefficient matrix, thus lending the necessary stability to the numerical process while restoring higher-order solution at convergence.

All the schemes but Hybrid require two upstream nodes for each cell-face, which will involve a value outside the solution domain for a near-boundary control volume. Therefore in the program the Hybrid scheme is used for all the control volumes adjacent to boundaries.

3.2.3 Approximation of diffusion terms

The diffusive flux IfD in Eq. 3.3 can be divided into two parts

ID = IDN  + IDC
 f    f      f
(3.9)

The first part IfDN contains only one term which has the first derivative of ϕ in the direction ”normal” to the cell-face f. It can be written, for the w-face for example, as:

IDwN =  Dw (ϕP -  ϕW );withDw  = [Γ ϕ∕ΔV ]w
(3.10)

The second part IfDC contains all the other terms. Only the normal derivative diffusion flux, IfDN, is coupled with the convective flux, IfC, to calculate the main coefficients of the difference equations, while the cross-derivative diffusion flux, IfDC, is treated explicitly as a pseudo-source term to avoid the possibility of producing negative coefficients in an implicit treatment.

3.2.4 Approximation of source terms

The source terms Sϕ are linearised as

      U     P
Sϕ = Sϕ + S ϕ ϕP
(3.11)

The coefficient SϕP is defined so that it is always not more than zero for all the conservation equations, and SϕU always assumes non-negative values for both k- and ϵ-equations. This enhances the stability of the numerical process and prevents the calculated values of k and ϵ from becoming negative in low turbulence regions. The volume integral of the source term can therefore be approximated as

∫
    SϕdV  ≃ SUϕ ΔV  + SPϕ ΔV ϕP
 ΔV
(3.12)

3.2.5 Final form of discretised equations

After replacing all the terms in Eq. 3.2 by their discretised analogies, the following difference equation results in:

        ∑
AP ϕP =     Anbϕnb + SU ,   nb = W, E, S,N,B, T
         nb
(3.13)

where

      ∑
AP =     Anb - SP
      nb
(3.14)

The main coefficients Anb that relate the principal unknown ϕP to its neighbours ϕnb contain the combined contribution from convection and normal diffusion. The physical source term (Eq. 3.12) and the cross-derivative diffusion flux IDC (Eq. 3.9) are included in the coefficients SU and SP .

In order to stabilise the iterative solution process, it is often necessary to under-relax the current solution by

ϕP := ϕoP + α ϕ(ϕP - ϕoP)
(3.15)

where αϕ([0,1]) is an under-relaxation factor and the superscript ”o” refers to the old value at the previous iterative level. Introducing Eq. 3.15 into Eq. 3.13 leads to an under-relaxed equation which is of the same form as Eq. 3.13 except that the coefficients SU and AP are replaced by

            1- αϕ
SU := SU +  ------AP ϕoP
             α ϕ
(3.16)

AP  :=  AP ∕αϕ
(3.17)

3.3 Boundary Conditions

The most frequently encountered boundary conditions, i.e., inflow, outflow, symmetry and rigid wall conditions, are implemented in the program. The pressures at all the boundary nodes are evaluated by linear extrapolation from values at interior nodes. The program also provides a provision allowing the user to specify boundary conditions other than these four types.

3.3.1 Inflow planes

At the inflow planes, the value of each variable ϕ is prescribed. Furthermore, the program defines the following ’n’ factors to normalise the residual:

     ∑               ∑
F1 =    m˙ii,    Fn =     ˙miiϕn,ii,    n = 2,3,⋅⋅⋅,n
      ii               ii
(3.18)

with

ϕ = ϕin
(3.19)

where the summation index ii runs over all the inlet computational nodes, ṁ is the mass flux and the subscript ”in” refers to the inlet values. Fn (n=1, 2, ... or n) is set to 1 if the corresponding value calculated by Eq. 3.18 is 0.

3.3.2 Outflow planes

At the outflow planes, the stream-wise gradients of all variables are set to zero, implying a fully developed flow condition. In addition, the velocities V j(j = 1,2,3) are corrected in the following way to ensure the overall mass conservation at the outflow planes:

Vj,ii := Vj,iiFa
(3.20)

         ∑
Fa = F1∕    ˙mout,ii
          ii
(3.21)

where ii runs over all the outlet computational nodes and ṁout is the outflow mass flux.

For pressure outflows with a pressure loss coefficient ζ, the pressure at the boundary is defined as

pBC = p∞  + 1ζρu¯2
            2
(3.22)

where p is the set pressure value, ζ is the value set for the pressure loss coefficient, ρ is the fluid density and u is the bulk outflow velocity, i.e. the average outflow velocity accounting for all the cells that are part of the same outflow.

3.3.3 Symmetry planes

At the symmetry planes, the normal velocity component and the normal gradients of other variables are set to zero. The convective and normal diffusive fluxes IfC (Eq. 3.4) and IfDN (Eq. 3.9) are also explicitly set to zero, where ”f” refers to the boundary surface.

3.3.4 Rigid wall

At the rigid walls, the no-slip condition is used for laminar flows, i.e., the velocity components of the flow are set to those of the wall. The standard wall-function approach (Launder & Spalding1974) is used for turbulent flows. In this, the resultant wall shear stress ⃗τw is related to the flow velocity vector ⃗V by

⃗τw = - λw ⃗VP
(3.23)

where

      {                         +
        μ ∕yP                if yP < 11.6
λw =    ρC1μ∕4k1∕P2κ∕ ln(Ey+P)  otherwise
(3.24)

 +      1∕4 1∕2
yP = ρC μ kP  yP∕μ,     κ = 0.41,  E  = 9.793
(3.25)

the subscript P refers to the first control volume centre from the wall, and yP is the normal distance from the wall. Further, the diffusive flux of k is set to zero at the wall, and the near-wall values of the production rate G and the dissipation rate ϵ are determined from:

      --τw2-
GP  = κμy+
          P
(3.26)

       3∕4 3∕2
ϵP = C-μ-k-P-
       κyP
(3.27)

3.3.5 Total pressure boundary condition

For inflow conditions one can define the total pressure P0 or total temperature T0. The mass flow rate is not known a priori. The static pressure at the inlet pi is extrapolated from inside as long as subsonic conditions prevail at the inlet. In case of a supesionic flow, pi must be defined by the user.

Ideal gas

For an ideal gas, isentropic expansion (pvγ = constant) from P0 to pi is assumed and knowing this relationship, the Mach number can be solved for.

      (              )-γ-
P0-        (γ --1)  2 γ-1
 pi =  1 +    2   M
(3.28)

Velocity can then be obtained by substituting ρ = ρ0(p∕p0)1∕γ

  2   (     ) (  ) (   (   ) 1    )
u--=   --γ--    1-   p   p-  γ - p
 2     γ - 1    ρ     0  p0
(3.29)

Incompressible or weakly compressible fluids

For a liquid phase, either with compressible and incompressible formulation, the density is a very weak function of pressure, or in other words, density remains almost constant even for large variations in pressure. The relationship between P0 and pi is given by the standard form of the Bernoulli equation,

P     p   1
-0-=  -i+ --u2
ρ0    ρi  2
(3.30)

Mixture of compressible fluids

The P0 to pi relationships for the two equations of state depends on the constant entropy expansion process, however they are both based on the general form of Bernoulli’s equation.

The general form of the Bernoulli equation can be derived as follows from the conservation of mass and conservation of total energy equations for steady one-dimensional compressible flow where entropy producing terms such as heat conduction, heat sources, viscous dissipation etc. have been set to zero.

d-(ρu) = 0
dx
(3.31)

and

   (                )       (               )
-d-           1- 2       d--      p-  1-2
dx  (ρe + p+  2ρu )u  =  dx  (e+  ρ + 2u )ρu  =  0
(3.32)

Using the product rule and using the conservation of mass equation, we obtain,

 d     p   1  2
dx-(e+ ρ-+ 2-u ) = 0
(3.33)

Since we are interested in the process from state 1 to 2 we can write the above equation as,

      ( p)    ( u2)
de+ d  ρ-  + d  -2-  = 0
(3.34)

Temperature:  To derive the relationship between the temperatures at the two states, we rewrite the equation in terms of enthalpy h = e + p∕ρ and integrate from state 1 to 2,

∫ 2     ∫ 2  (  2)
   dh +    d  u--  = 0
 1       1     2
           u22   u21
 h2 - h1 + 2--- 2--= 0                             (3.35)
For constant Cp and assuming state 2 to be stagnation condition (u2 = 0), we obtain
            2
T0 = T1 + u-1-
          2Cp
(3.36)
If Cp is not constant we can write the enthalpy difference as,
         ∫  T2       ∫ T1        ∫ T2        --
h2 - h1 =     CpdT -      CpdT =      CpdT = Cp [T1 : T2 ](T2 - T1)
           T0         T0           T1
(3.37)
Using Eq. (3.37) in Eq. (3.35) and as setting as before state 2 to be the stagnation state (u2 = 0) we get,
          u21
T0 = T1 + ----
          2Cp
(3.38)
Pressure:  To derive the relationship between the pressures at the two states, one has to incorporate the relationship between pressure and density for an isentropic expansion. We start with the thermodynamic relationship,
T dS = de + pdv
(3.39)
For an isentropic process (dS = 0), we obtain,
                 (  )
de = - pdv = - pd  1- =  p-dρ
                   ρ     ρ2
(3.40)
Using the above relationship and substituting into Eq. (3.34) and further simplifying we obtain,
dp    ( u2)
---+ d  ---  = 0
ρ        2
(3.41)
In the above equation if ρ = constant, integrating from state 1 to 2 gives the standard Bernoulli equation (Eq. (3.30)) for incompressible flow and assuming p ργ gives Eq. (3.28). For a homogeneous multiphase mixture one can use Eq. (3.41) to calculate the inflow velocity as follows. We write the mixture density as,
-1-   ∑  Yk-
ρm =     ρk
       k
(3.42)
where k is the phase index, Y k is the mass fraction of phase k, and ρk is the density of phase k. We assume that the mixture has a single pressure, temperature and velocity. The mass fractions of each phase is assumed to be specified. Integrating Eq. (3.41) between states s1 and s2 for a two phase mixture, we obtain,
∫ s2Y      ∫ s2Y       u2    u2
    -1dp +     --2dp = -s1-  -s2
 s1 ρ1      s1 ρ2       2     2
(3.43)
Assuming the phase 1 to be an ideal gas and phase 2 to be a constant density liquid, and further noting that the mass fractions do not change during the expansion, we obtain,
                 (      )(            )     2     2
Y2(ps2 - ps1)+ Y1  -γ---   -ps2 - -ps1  =  us1-  us2
ρ2                 γ - 1   ρ1,s2   ρ1,s1      2     2
(3.44)
Setting state 2 to be stagnation condition (us2 = 0) and state 1 the inflow conditions (u1 = uinflow, we obtain,
 2                        (      ) (            )
uinflow-   Y2-               --γ--    -ps2-   -ps1-
   2   =  ρ2(ps2 - ps1)+ Y1 γ - 1    ρ1,s2 - ρ1,s1
(3.45)
Linearization for pressure coefficient and flux update

To find the pressure coefficient arising out of specifying the inflow velocity using the isentropic relation for the ideal gas, we take the derivative of Eq. (3.29) with respect to pressure to obtain,

          (   ) (   ) 1
ui∂ui-= -   -1    p0  γ
   ∂p       ρ0     p
(3.46)

Multiplying both sides by the square of the inflow cell area AiAi and defining the flux as ũ = uiAi, we get,

       (     ) (   ) 1γ
∂˜u-= -   AiAi-   p0
∂p       ρ0˜u     p
(3.47)

or rewriting in terms of ρ we get,

∂˜u     AiAi
∂p-= - -ρ˜u--
(3.48)

3.4 Solution Procedure

3.4.1 Convergence criterion

The residuals are defined for the continuity equation as:

      [                                 ]
       ∑                                  1∕2
RL2 =     (Ce - Cw + Cn - Cs + Ct - Cb)2ii   ∕F1
        ii
(3.49)

RL ∞ = maxii (|Ce - Cw  + Cn - Cs + Ct - Cb |ii)∕F1
(3.50)

For all the other equations the residuals are defined as:

          [                             ]1∕2
           ∑   ∑                       2
Rm (L2 ) =    (    Anbϕnb + SU - AP ϕP)ii    ∕Fm     (nb = E, W,N, S,T, B)
            ii  nb
(3.51)

                 ∑
Rm (L ∞) = max (|   Anbϕnb + SU - AP ϕP |ii)∕Fm     (nb = E, W,N, S,T, B)
            ii   nb
(3.52)

where ii runs over all the computational nodes and Fm (m=1 to 7) are the normalisation factors defined in Eq. 3.18. The calculation is declared converged when,

Rmax =  MAX  (R1, ...,R7 ) ≤ EPS
(3.53)

where EPS is the convergence criterion prescribed by the user. TransAT makes use of the L norm, which is a stricter convergence condition than the L2 norm.

3.4.2 Calculation sequence

The sequence in which the calculation is carried out is as follows:

  a.
Initialise all field values by guess.
  b.
Solve the V 1,V 2 and V 3-momentum equations using the guessed pressure field.
  c.
Solve the pressure-correction equation to obtain the pressure-correction at the cell-centers; correct the convective fluxes at the cell-faces, the velocities and pressure at the cell-centers.
  d.
Solve the k- and ϵ-equations and update μt, if the flow is turbulent.
  e.
Solve the scalar S-equation, if required.
  f.
Return to step b with updated field values.

Sequence of steps b to f is repeated until the convergence criterion (Eq. 3.53) is satisfied.

A complete sequence of steps b to f is defined as an outer loop, and a solution of the system of Eqs. 3.13 for each variable over the entire solution domain as an inner loop. The system of Eqs. 3.13 is solved by using the strongly implicit procedure of Stone (1968). For all variables but the pressure-correction, one inner loop is normally enough to reduce the corresponding residue to a reasonable low level, especially in case of highly convective flows. However, sufficient inner loops are required for the pressure-correction p-equation, mainly due to the fact that the p-equation is of diffusion type in which the coefficient matrix is symmetric but without any convective fluxes, thus making convergence relatively slow. The inner loops of solving the p-equation are stopped when either of the following two criteria is satisfied (van Doormaal & Raithby1984):

  m      i
R   ≤ λR
(3.54)

m > NSWC1
(3.55)

where Ri and Rm are the residuals calculated in the same way as in Eq. 3.52 at the initial stage and after the m-th loop, λ = 0.1 and NSWC1=20 are used in the program.

Chapter 4
Reactive Flow Modelling

4.1 Conservation equations

4.1.1 Species equation

General formulation: Mass fraction

The species equation for every species k can be written as (Poinsot & Veynante2001)

∂ρYk-+  -∂-(ρ(um + V  )Y ) = ˙ωm
 ∂t     ∂xi    i    k,i  k     k
(4.1)

Where V k,im is the diffusion velocity and the source term ω˙km has units of [kgsm3]. The flow velocity uim is defined in mass weighted fashion as

 m    N∑     k
ui =     Ykui
      k=1
(4.2)

where uik is the velocity of species k. The sum of the mass weighted diffusion velocities is then zero by definition:

∑N   m
    Vk,iYk = 0
k=1
(4.3)

The diffusion velocities are determined by the multicomponent diffusion equation which can be derived for a low density ideal gas (Williams1958R.J. Kee2003):

          N∑  XkXj                     ∇p
∇Xk   =      ------(vj - vk + (Yk - Xk)--
          j=1  Dkj                      p
           ρ ∑N                N∑  XkXj   DT,j   DT,k  ∇T
          +--    YkYj(fk - fj) +    [-----(------ -----)]----
            pj=1               j=1 ρDkj   Yj     Yk    T
(4.4)

where DT is the thermal diffusion coefficient, Dkj is the binary diffusion coefficient of k in j and fk is the volume force. The diffusion velocity can also be approximated with Fick’s law where the binary diffusion coefficients is replaced by a diffusion coefficient of species k in the mixture. Furthermore pressure and temperature gradients are neglected. It is important to realise that Fick’s law is an approximation for the multicomponent diffusion equation and the constraint for diffusion velocities. Eq.(4.3) is no more fulfilled in general (Poinsot & Veynante2001).

 m          ∂Yk
Vi,kYk = - Dk ∂x-
               i
(4.5)

The species equation is then:

∂ρY    ∂ρum Y      ∂     ∂Y
---k-+ ----i-k-=  ---ρDk --k-+ ω˙mk
 ∂t      ∂xi      ∂xi    ∂xi
(4.6)

General formulation: Mole fraction

The species equation in terms of mole fractions reads:

∂cXk     ∂
----- + ----(cXk (uci + V ck,i)) = ˙ωck
  ∂t    ∂xi
(4.7)

Where c is the total concentration, Xk the mole fraction of species k and ω˙k has units of [molsm3]. The concentration weighted velocity is then:

 c   ∑N     k
ui =    Xku i
     k=1
(4.8)

For the diffusion velocities the resulting constraint is:

∑N   c
    Vi,kXk = 0
k=1
(4.9)

In general V i,kc = V i,km is only true if the molar weights of all species are equal. Again the diffusion velocities can be approximated with Ficks law where the same assumptions as in the mass fraction formulation are made:

             ∂Xk
Vic,kXk  = - Dk----
              ∂xi
(4.10)

This finally results in:

           c
∂cXk-+  ∂cuiXk-=  -∂-cDk ∂Xk-+  ˙ωck
 ∂t       ∂xi     ∂xi    ∂xi
(4.11)

As this equation was derived under the assumption of ideal gas, constant temperature and constant pressure, it can only be applied in the case of constant concentration.

General formulation: Concentration

One can write the transport equation directly for concentration of species k (ck)

∂c     ∂
--k-+ ----(ck(uci + Vck,i)) = ω˙ck
 ∂t   ∂xi
(4.12)

The concentration weighted velocity is then

    ∑N
uc=     ckuk
 i  k=1    i
(4.13)

with the constraint:

N∑    c
   Vi,kck = 0
k=1
(4.14)

The diffusion velocity is then written according to Ficks law as:

Vic,kck = - Dk ∂ck
            ∂xi
(4.15)

Which results in the concentration equation:

∂ck    ∂ucck    ∂    ∂ck    c
-∂t-+  -∂ix--=  ∂x-Dk ∂x--+ ˙ωk
          i      i      i
(4.16)

For constant concentration, Eq.(4.16) and Eq.(4.11) are the same but the latter can be used to derive an transport equation for a passive scalar.

Passive scalars: Mass fraction

In a passive scalar formulation the flow velocity is not determined by mass or concentration averaged velocities. The velocity field ub is assumed to be known and independent of the species. Therefore the diffusion velocities are simply added to this base flow field. In this case Dk is the diffusion coefficient of species k in the base flow. For passive mass fractions the following equation is solved:

∂ρYk-   ∂ρubiYk-  -∂--    ∂Yk-   m
 ∂t  +   ∂xi   = ∂xi ρDk ∂xi + ˙ωk
(4.17)

One component can be constrained to ensure mass conservation.

Passive scalars: Concentration

As mentioned above the mole fraction formulation of the species equation can not be used for a general passive scalar as the total concentration can vary in space and time. However by replacing the concentration averaged velocity by the base velocity uib in Eq. (4.16) one obtains a transport equation at the passive concentration (Fox2003).

∂ck-+  ∂ubick=  -∂-D  ∂ck-+ ˙ωc
 ∂t     ∂xi    ∂xi  k∂xi    k
(4.18)

Passive scalars: Mixture fraction and progress variable

In order to apply the mixture fraction approach the passive scalar must obey an equation of the form (Fox2003)

∂ Φ    b∂Φ        ∂2Φ
-∂t + ui∂x--= Dk ∂x-x- + ˙ωΦ
          i        i i
(4.19)

Assuming constant D and using the continuity equation this can be rewritten in a form suitable for TransAT

∂ Φ   ∂Φubi    ∂   ∂ Φ          ∂ubi
-∂t + -∂x-- = ∂x-D ∂x--+ ω˙Φ + Φ ∂x--
         i      i     i           i
(4.20)

This equation is solved for the mixture fraction and progress variable where in the case of the mixture fraction the source term will vanish.

Chemical source term

In a laminar flow the chemical source term is determined by Arrhenius rate constants:

          r
˙ωm = W   ∑  ν ˙ω
 k     k l=1  il l
(4.21)

where r is the number of reactions, Wk is the molecular weights of component k, νkl = νkl′′- νkl is the difference between forward and backward stoichiometric coefficients and ˙ωl is given by

        ∏n ( ρYj) ν′jl     ∏n ( ρYj) ν′′jl
ω˙l = kfl     ----   -  kbl     ----
        j=1  Wj          l=1  Wj
(4.22)

and kfl,kbl are the rate coefficients of the forward and backward reaction, respectively. It results that the sum i=1nω˙km = 0. The rate coefficients are usually represented in the Arrhenius equation

      --Ea
k = Ae RT
(4.23)

where A is the pre-exponential factor, T is the temperature, Ea is the activation energy and R is the gas constant.

Definition of mixture fraction and progress variable based on concentration

Besides the specific derivation for a combustion system, mixture fraction and progress variable for a system can be derived in a general way. One tries to make a coordinate transformation from the physical space (in terms of concentration or mass fractions) to a space where as few transported scalars as possible have a chemical source term. Consider a system with three species: reactant A, reactant B, and the product P which take part in the reaction:

νAA +  νBB = νP P
(4.24)

All three species fulfill a transport equation of the form

∂ Φ   ∂Φub     ∂   ∂ Φ          ∂ub
--- + ----i = ---D ----+ ω˙cΦ + Φ --i-
 ∂t    ∂xi    ∂xi  ∂xi          ∂xi
(4.25)

where Φ can be replaced by the concentration of the species (cA,cB or cP ). The inlets of the system contain either species A with concentration cA0 or species B with cB0. A conserved scalar, (mixture fraction Z) can then be defined as:

Z  = νAcB---νBcA-+-νBcA0-
        νAcB0 + νBcA0
(4.26)

To obtain a transport equation for the mixture fraction, Φ in Eq.(4.25) can be replaced by Z where ω˙Z = 0. Without a reaction the pure mixing solution is then cA = cA0(1 -Z) and cB = cB0Z (Poinsot & Veynante2001).
As there is only one reaction in this system the source term of all species is directly linked to the reaction rate of this reaction. R is defined as the molar reaction rate [molsm3] with which species the unity stoichiometric coefficient would be produced or destroyed. If the reaction rates are determined with an Arrhenius formulation we have R = kcAcB in the case of a one step reaction of two reactants. One can then define the vector ϒ as ω˙ic = ϒiR. For the one step reaction this yields:

    (       )
    (  - νA )
ϒ =    - νB
        νP
(4.27)

The progress variable (ξ) describes the extend of reaction. It is defined such that it is zero in the unreacted state and one when the reaction has finished. The source term of the progress variable is defined to be R∕γ where γ is an appropriate scaling factor to obtain the correct scaling (0 < ξ < 1). The concentration vector c can be written as a linear combination of mixture fraction and progress variable (Fox2003):

c = c0 + MZ Z + M ξξ
(4.28)

c0 and MZ can be determined from the pure mixing solution as

     (      )
c0 =    cA0
         0
(4.29)

and

      (  - cA0 )
MZ  =    c
          B0
(4.30)

If we write the the transport equation (4.25) for c the source term which remains is Mξ∕γR as Z and c0 are conserved scalars. As the source term for the equations was defined as ϒR, Mξ is equal to γϒ. For Eq.(4.24) the appropriate scaling factor is (Fox2003)

γ = ----cA0cB0----
    νAcB0 + νAcB0
(4.31)

This finally yields the concentration of all three species expressed with mixture fraction and progress variable:

cA  =   cA0(1- Z - (1 - Zst)ξ)                        (4.32)
cB  =   cB0(Z - Zstξ)                                (4.33)
c   =   ν c  Z  ξ                                    (4.34)
 P       P B0  st
where Zst is the stoichiometric mixture fraction defined as
Z  =  ----νBcA0-----
 st   νAcB0 + νBcA0
(4.35)

The source term for the progress variable in Arrhenius formulation is then

˙ωξ  =   R-                                              (4.36)
        γ
        kcAcB-
    =   cB0Zst                                          (4.37)
                1 - Z        Z
˙ωξ  =   kZstcB0(-------- ξ)(----  ξ)                    (4.38)
                1- Zst      Zst
Definition of mixture fraction and progress variable based on mass fraction

If we define the mixture fraction for mass fractions the mixture fraction definition is slightly different

    νAMAYB   - νBMBYA   + νBMBYA0
Z = ------ν-M--Y---+-ν--M--Y--------
           A  A B0    B   B A0
(4.39)

As the mixture definition is now in terms of masses the mass stoichiometric coefficients νiMi have to be used where Mi is the molar weight of species i. Again we have the species A,B and P which take part in the reaction defined in Eq.(4.24). The inlets of the system contain either species A with mass fraction Y A0 or species B with Y B0. Again a linear combination of mixture fraction and progress variable can be defined for the mass fractions:

Y = Y   + M  Z +  M  ξ
      0     Z       ξ
(4.40)

From the pure mixing solution we can obtain:

     (      )
Y0 =    YA0
         0
(4.41)

and

      (       )
        - YA0
MZ  =    YB0
(4.42)

Using the same definition as above reaction vector is then

     (    M ν  )
        - -Aρ-A-
ϒ  = |(  - MBρνB |)
         MPνP-
           ρ
(4.43)

As the product mass fraction has already the desired scaling (0 < Y P < 1) it can directly be used as reaction progress variable. The scaling factor γ can then be obtained from:

˙ωξ = ˙ωP = R ∕γ
(4.44)

Therefore

1∕γ = MP νP ∕ρ
(4.45)

and

     (    MA-νA )
     |  - MMP ννP |
M ξ = ( - MPBνBP )
           1
(4.46)

From this we get:

YA  =   YA0(1-  Z)-  MA-νA-ξ                        (4.47)
                     MP νP
                MB-νB-
YB  =   YB0Z -  MP νP ξ                             (4.48)
YP  =   ξ                                           (4.49)
The source term for the progress variable in an Arrhenius form reaction rate expression is then
        R
˙ωξ  =   γ-                                                       (4.50)

    =   νPMP--kcAcB                                              (4.51)
          ρ
        νPMP-k-ρYA-ρYB-
    =     ρ    MA  MB                                            (4.52)
        νPMP ρk               MA νA           MB νB
˙ωξ  =   --------(YA0(1- Z )- ------ξ)(YB0Z  - ------ξ)           (4.53)
        MAMB                 MP  νP           MP νP
For the eddy dissipation the source term is always divided by one of the three mass stoichiometric coefficients. Therefore the mass stoichiometric coefficients can be normalised by the the mass stoichiometric coefficient of component A:
r = νBMB--                                   (4.54)
    νAMA
    νPMP
s = ν-M---                                   (4.55)
     A  A
According to the eddy dissipation theory, the mass reaction rate of A is
           ϵ    (    Y   BY  )
ω˙mA = - Aρ -min  YA, -B-,---P-
           k          r    s
(4.56)

From this we conclude that

       ω˙m
R = - ---A--
      νAMA
(4.57)

Therefore the source term for the progress variable is:

        R
ω˙ξ =   γ-                                            (4.58)
                m
    =   νP-MP-ω˙A                                     (4.59)
        νAMA   ρ(             )
           ϵ-        YB- BYP--
ω˙ξ =   sA kmin  YA,  r ,  s                          (4.60)

4.1.2 Temperature equation

A balance equation can be also written for the temperature

   ∂T               ∂p               ∑n
ρcp∂t-+ ρcpv ⋅∇T  = ∂t-+ ∇  ⋅(λ ∇T )-     ρcp,iD ∇Yi ⋅∇T  + ωT
                                     i=1
(4.61)

where λ is the thermal conductivity and cp is the heat capacity at constant pressure. An equation can be written for temperature source term

     1-∑ r
ωT = cp    Qkwk
        k=1
(4.62)

where Qk = - i=1nνikWihi is the heat of combustion of reaction k.

The vector (Y 1,Y 2,,Y n,T) can be called the “reactive scalars” vector ψ.

If different diffusivity must considered, the convection velocity of the species must be corrected by adding a correction velocity V ic = k=1NDk∂Yk
∂xi to ensure global mass conservation.

4.2 Turbulent reacting flow

4.2.1 RANS

When the species have different densities, the usual Reynolds average gives new unclosed terms ρuj in the mass balance equation corresponding to the correlation between density and velocity fluctuations. To avoid this difficulty, mass-weighted averages (called Favre averages) are usually preferred,

     ---
f^=  ρf-
     ρ
(4.63)

where f is a generic quantity. Splitting Y i into a Favre mean ^Y
 i and a fluctuation Y i′′ Eq.(4.17) becomes

-∂ ^Yi  --            ---------     ( --    )  --
ρ----+ ρ^v ⋅∇ ^Yi = ∇ ⋅(ρDi∇Yi )- ∇ ⋅  ρv′′Yi′′ + ρS^i
  ∂t
(4.64)

In this equation the r.h.s. terms must be modelled.

The species turbulent fluxes  ′′ ′′
v Yi are generally closed using a classical gradient assumption

v′′Y ′′i = - Dt ∇ ˜Yi
(4.65)

where Dt is a turbulent diffusivity which is modelled by analogy to the eddy viscosity as Dt = νStct where Sc is a turbulent Schmidt number. For a reactive scalar this closure is not always acceptable and can be be improved upon (Peters2000).

The laminar diffusive fluxes are generally neglected or retained by adding a laminar diffusivity to the turbulent viscosity in Eq.(4.65).

--
ρ^v ⋅∇ ^Yi = ρ¯D¯i ∇ ˜Yi
(4.66)

where Di is a “mean” species molecular diffusion coefficient.

4.2.2 LES

The same equation obtained for RANS approach can be obtained with LES, substituting the filter operation to the Favre average and the same terms needs a closure that is achieved using the sub-grid scale viscosity and additional sub-grid scale model for source terms.

4.3 Non–Premixed Reactive Flows

In combustion this case is also called Diffusion Flames since the reactants are initially separated and the reaction takes place in a narrow region or surface.

In non-premixed flows the reactions are strictly connected with the mixing and by comparing these two phenomena we can distinguish different reaction regimes. With reference to Figure 4.1, we define “slow” chemistry for the reactions with time scales larger than the micromixing time scale, “fast” chemistry when the chemical time scales are all smaller than the Kolmogorov time scale and “finite-rate” chemistry in the middle. In the following we will consider first the two limiting cases, “frozen” (very-slow) chemistry and “infinitely-fast” chemistry.


PIC


Figure 4.1: Comparison example between flow and chemical time scales. The flow time scales range from the Kolmogorov time scale τν, through the turbulence time scale τu, up to the mean residence time τres. The micromixing time τm will usually lie between τν and τu. Figure taken by Fox (2003)


For non–premixed reacting flows, it is often possible to define a mixture-fraction vector Z that can be employed to develop a simpler set of equations with chemical source term. In Fox (2003) a general method for finding the mixture-fraction vector (when it exists) for a given set of initial and inlet conditions. When a mixture-fraction vector exists, the system is fully described for infinitely-fast chemistry. For finite-rate chemistry a reaction-progress vector ξ must be taken into account. This vector is null for all initial and inlet conditions.

The importance of the mixture fraction is that it is a conserved quantity, a scalar quantity that is neither created nor destroyed so the equations for Zi reduces to

 ∂Zi
ρ-∂t-+ ρv ⋅∇Z  = - ∇ ⋅(- ρD∇ZYi  )
(4.67)

Many computational models for combustion are developed with a single mixture fraction which can be used only for two-feed systems. If more than two feeds enter into a combustion chamber, the concept of single mixture fraction can no longer be used. In fact, even if it can formally be extended to multiple mixture fraction variables, it becomes less attractive for modelling because the resulting flamelet equations are defined in a multi-dimensional space.

The usual definition of mixture fraction for non-premixed combustion is best derived for an homogeneous system in the absence of diffusion. In this case the reaction equation for complete combustion of hydrocarbon fuel CmHn can be written as

 dY       dY
--′-O- = --′-F-
νOWO     νF WF
(4.68)

where F stands for fuel and O for oxidiser and νare the stoichiometric coefficients. For a homogeneous system this equation may be integrated to obtain

νYF - YO  = νYF,u - YO,u
(4.69)

where ν = ν′OW′O--
νF WF is the stoichiometric oxygen-to-fuel mass ratio and the subscript u denotes the initial conditions in the unburnt mixture. If the diffusivities of fuel and oxidiser are equal, 4.69 is valid also for spatially inhomogeneous systems such as a diffusion flame. So the quantity νY F - Y O is a conserved variable and the mixture fraction can be expressed as

      β - βO,u
Z =  β-----β----
      F,u    O,u
(4.70)

where β is a conserved variable. In our case it becomes

Z = νYF----YO-+-YO,2
      νYF,1 + YO,2
(4.71)

where subscript 1 denotes the the fuel stream and 2 the oxidiser stream.

The stoichiometric mixture fraction Zst can be defined as the value of Z for a stoichiometric mixture. In this case the quantity νY F - Y O vanishes and we can obtain

       YO,2
Zst = νY---
         F,1
(4.72)

The mixture fraction can be related to the commonly used equivalence ratio r, which is defined as the fuel-to-air ratio in the unburnt mixture normalised by that of the stoichiometric mixture

    νYF,u
r =  YO,u
(4.73)

The mixture fraction can be also expressed starting from a quantity related to chemical elements. In fact, while the mass of chemical species may change due to chemical reactions, the mass of elements is conserved. (Poinsot & Veynante2001Peters2000).

Combustion regime indicators

4.3.1 Slow Chemistry

Arrhenius approach

In our discussion, the most important closure will be the source term ^
Si. The simpler model is to neglect the fluctuations and write the Arrhenius equation with the mean quantity Y˜i. This approximation is called the “Arrhenius approach” and the closure relation is simply

------     --
ωi(Y ) = ωi(y)
(4.74)

This model is relevant only for laminar flow or when chemical time scale are larger than turbulent time scale (τϕ > τu, low Damkohler number limit)

4.3.2 Infinitely Fast Chemistry

In this representation one or more chemical reactions are specified but they are assumed to be infinitely fast. In modelling of combustion processes, this approach is known in terms of the “mixed is burned” assumption. The reactions are completely controlled by the mixture fraction since reactants are converted directly to final products with a one-step description. So it is enough to solve the balance equation for the mixture fraction variable and then calculate mass fractions and temperature. This approach is valid only for low Mach number flow with constant thermodynamic pressure, same species heat capacities, Lewis number equal to unity and molecular diffusivities all equal to D.

Laminar flames

The equation for the mixture fraction is:

 ∂Z
ρ--- + ρv ⋅∇Z =  ∇ ⋅(ρD ∇Z )
 ∂t
(4.75)

Once the solution is known in the entire flow field, the flame surface is defined as the surface of stoichiometric mixture, which is obtained by setting Z(x,t) = Zst. In the vicinity of that surface the reactive-diffusive structure can be described by the so-called “flamelet equations”

 ∂Yi-   -ρ-χ-∂2Yi
ρ ∂t =  Lei2 ∂Z2 +  ωi
(4.76)

where χ is the scalar instantaneous dissipation rate defined as χ = 2D|∇Z|2.

Turbulent flames

The complexity added by turbulence comes from the averaging (or filtering) procedure. To determine the average mass fractions and temperature, the mean value of the mixture fraction is not sufficient: higher moments are needed. The two possible approaches are

Eddy Dissipation Concept

The eddy dissipation concept, devised by Magnussen and Mjertager (1997), directly extends the Eddy Break-Up model to non-premixed combustion. The main idea is to replace the chemical time scale of an assumed one-step reaction by the turbulent time scale k
ϵ. Thereby the model eliminates the influence of chemical kinetics, representing the fast chemistry limit only. The fuel mean burning rate is estimated from fuel (Y F ), oxidiser (Y O) and products (Y P ) mean mass fractions as

              (     ¯     ¯ )
ω¯F =  A¯ρ-ϵmin  Y¯F ,YO-, B-YP
        k           ν  1+ ν
(4.77)

where A and B are modelling constants. For a highly exothermic reaction (combustion), BY¯P-
1+ν is included because it is proportional to the heat release and therefore models the temperature dependence of the reaction. For isothermal reactions this term can be dropped.
The eddy dissipation concept is based on the assumption that scalar mixing lengthscales and turbulent lengthscales are the same, i.e. the concentration is uniform within a Kolmogorov length-scale. This is only valid for unity Schmidt number (observed in gas phase reactions) but not for Sc>>1 which is the case in liquid fluids. In this case k and ϵ have to be replaced by the scalar variance (< )2 >) and by the scalar dissipation rate (ϵϕ). In order to model (--ϵϕ′2--
<(Φ )>) one can derive an algebraic equation depending on the Schmidt number. There also exist more complex models which involve the solution of transport equations for the variance in different parts of the Kolmogorov energy spectrum.

Chapter 5
Boussinesq Approximation

The Boussinesq approximation can be used to solve natural convection problems without having to solve the full compressible Navier-Stokes equations. This is achieved by selectively considering local variations in density, due to local variations of temperature, only for the gravity term in the Navier-Stokes equations. All remaining properties and terms remain unaffected.

The variation of density with temperature is linearized around a set reference density ρ0 and a corresponding reference Temperature T0. The linear constant of proportionality is the thermal expansion coefficient, denoted by β0.
Variations of density are then approximated through the linear relation

^ρ(T) = ρ0 ⋅(1 - β0(T - T0))
(5.1)

Please note: The linear approximation of density variations with temperature is generally only valid for a small temperature range. The expected temperature range should also be considered for the choice of reference properties and thermal expansion coefficient. For large temperature variations it is recommended to solve for compressible flows.

Chapter 6
Non-Newtonian Models

For Newtonian fluids, the viscosity is a constant that defines the ratio of stress-to-strain-rate is constant. This property is false for so-called Non-Newtonian fluids. This section describes the capabilities of TransAT for handling such Non-Newtonian fluids.

Bingham model

The Bingham model may be used to describe plastic materials, which feature a yield stress. In this model, the viscosity η may be written

    {
       ∞       if  |τ| ≤ τb
η =    ηp + τb˙γ- if  |τ| > τb
(6.1)

where γ˙ is the shear rate, τ is the shear stress. τb and ηp are two model constants standing for the yield stress and the plastic viscosity, respectively.

Casson model

The Casson model has been developed to simulate plastic materials. In this model, the viscosity is defined by

        1∕2
η1∕2 = τb--+ η1∕2
       ˙γ1∕2    p
(6.2)

where γ˙ is the shear rate, τ is the shear stress. τb and ηp are two model constants standing for the yield stress and the plastic viscosity, respectively.

Power-law model

The power-law model is a simple model for the simulation of pseudo-plastic or dilatant fluids. Viscosity in the power-law model reads

η = ηcγ˙n
(6.3)

where ˙γ is the shear rate, and ηc and the power law exponent n are constants of the model.

Carreau model

The Carreau model may be used for polymeric liquids. it reads

                  (        )n-1
η = η∞ + (η0 - η∞) 1 + λγ˙)2  2
(6.4)

where γ˙ is the shear rate. The zero-shear viscosity η0 , the infinite shear viscosity η ,the relaxation time of the fluid λ and the exponent n are model constants.

Cross model

It reads

η = -----η0----
    (1 + (λ˙γ)m)
(6.5)

where γ˙ is the shear rate. The zero-shear viscosity η0 , the relaxation time of the fluid λ and the exponent n are model constants.

Carreau-Yasuda model

The Carreau-Yasuda model may be used for the simulation of polymers. The viscosity is then written

                           a 1-n-
η = η∞ + (η0 - η∞ )(1 + (λ ˙γ) ) a
(6.6)

where γ˙ is the shear rate. The zero-shear viscosity η0, the infinite-shear viscosityη, the relaxation time λ and the exponents a and n are model constants.

Thixotropic model

The Thixotropic model may be used for the simulation of crude oil. The viscosity is is a function of a scalar quantity λ

η = η∞ + αλ
(6.7)

where the infinite-shear viscosity η and α are fluid properties and are considered to be model-constants. The scalar λ evolves according to the equation

D λ
Dt--= - aλ ˙γ + (1 - λ)b˙γ
(6.8)

where a is the break-down rate of λ, b is the build-up rate and Dλ-
Dt is the substantial derivative of λ.The value of λ ranges from 0 to 1. At 0 shear-rates, the value of λ is 1 and gradually as shear increases λ decays to 0.

Hershel-Bulkley model

The Hershel-Bulkley model may be used to describe In this model, the viscosity η may be written

    {
      ∞             if   |τ| ≤ τb
η =   k(γ˙)n-1 + τγb˙  if   |τ| > τb
(6.9)

where γ˙ is the shear rate, τ is the shear stress. τb ,k and n are three model constants. taub refers to the yield-stress of the fluid, n is an exponent and depending on whether it is greater than or less than 1 determines whether the fluid is shear-thickening or shear-thinning and k(γ˙)n-1 is the plastic viscosity.

Chapter 7
Viscoelasticity

For Newtonian fluids, the stress is simply the sum of the pressure gradient and the newtonian stress. A more general description is:

∂ρui-  ∂ρuiuj-  ∂-σij
 ∂t  +   ∂xj  =  ∂xj + ρgi
(7.1)

where σ is the total stress tensor, which can be decomposed:

σij = - pδij + τij
(7.2)

where p is the pressure and τ is the viscoelastic stress tensor.

For a newtonian or viscoplastic fluid (no elastic behaviour), the viscoelastic stress becomes the viscous stress:

              2
τisj = 2μsSij --μsSkk δij
              3
(7.3)

For viscoelastic fluids the viscoelastic stress becomes:

      s    e
τij = τij + τij
(7.4)

where τe is the polymeric or viscoelastic part, determined via the stress governing-equation depending on the viscoelastic model used.

In general, the response of a viscoelastic material cannot be described with a single characteristic time scale. This is why the more accurate viscoelastic models are the so-called multimode models which can be described as:

  e  ∑    e
τ  =     τm
      m
(7.5)

where each τme (mode) correspond to a specific elastic relaxation time scale and obeys its own constitutive equation.

The momentum equation then reads:

∂ρui   ∂ρuiuj     ∂p    ∂τsij   M∑  ∂ τeij,m
-----+ -------= - ----+ ----+     ------
 ∂t     ∂xj       ∂xi   ∂xj   m=1  ∂xj
(7.6)

The role of the viscoelastic solver in TransAT is to solve for the stress equation and couple it with the general momentum equation written above. The stress equation depends on the model used. The following models are available:

7.1 Olroyd-B. model

The stress equation (for one mode) for the Olroyd-B model is:

     τe          2
ˇτe + -- = 2GS -  -Gtr(S )I
     λ           3
(7.7)

The only model parameters that are required (for each individual mode) are the elastic time scale λ and the polymeric viscosity μp (since G = μp
 λ).

7.2 Exponential Phan-Thien-Tanner model (PTT)

The exponential PTT model writes:

                       (        )
 e       e    e          ϵtr(τ e)    e         2
ˇτ + ξ(S.τ + τ .S) + exp  --μ---- τ  = 2GS  - 3Gtr (S)I
                            p
(7.8)

7.3 Linear PTT model

The linear PTT is the linearized version of the exponential PTT and writes:

                    τe   ϵtr(τ e)            2
ˇτe + ξ(S.τe + τ e.S) +- + ------τ e = 2GS - -Gtr (S)I
                    λ      μp              3
(7.9)

For both PTT models the required parameters are (for each mode): ϵ, ξ, λ and μP .

7.4 Temperature correction using Arrhenius law

Usually viscosity correction using the Arrhenius law writes:

μ = μ0e- α(T -T0)
(7.10)

where α and T0 are given parameters.

7.5 Multiphase viscoelasticity

The viscoelastic solver is currently not available for homogeneous and algebraic slip multiphase models, only with level-set. The properties of each phase are weighted using the level-set function and the procedure is the same than for single phase viscoelastic flows.

7.6 Logarithm of Conformation Tensor (LCT)

The solution of the viscoelastic constitutive equation is made difficult by the high Weissenberg number problem (HWNP) which occurs when the Weissenberg number reaches moderately high values. The origin of this numerical is not well explained but seems to correspond to loss of positivity of conformation tensor A.

    1 - ξ
A = --G--τ + I

The formulation of the constitutive equation in terms of the matrix logarithm of the conformation tensor has the advantage of not exhibiting such unstable behavior. Transforming the stress equation to obtain the LCT equation results in an equivalent for the stress equation, that needs to be solved and coupled to the momentum equation.

Chapter 8
Radiation

The Rosseland Model (also known as diffusion approximation model) for radiation has been implemented in TransAT. This model can be quickly described as follow:

The Rosseland model assumes that the radiation can be taken into account as a temperature-dependent conductivity. This means that no additional equation is solved, but an effective thermal conductivity is calculated to account for radiative heat transfer.

The radiative heat transfer writes:

      16σT 3
qrad = ------∇T
        3β

where σ is the Stefan-Boltzmann constant and β is the Rosseland extinction coefficient. The diffusion approximation states that the radiative contribution to the heat transfer can be treated similarly to the conductive part provided that a radiative conductivity can be calculated. This conductivity depends on temperature (T3 contribution) and on the extinction coefficient β. In general β itself is temperature-dependent and must be computed. The Rosseland approximation states that β can be expressed as a polynomial of temperature:

    ∑N
β =    aiT i
    i=0

The ai are phase-dependent coefficients called the Rosseland coefficients, given as user-input. If the extinction coefficient is to be temperature-independent it is equal to a0 with all others Rosseland coefficients set to 0.

This is often a rough approximation because the radiative properties of media are also spectrum-dependent. This means that the above approximation for β for a set of Rosseland coefficients is only valid in a certain range of wavelength. To remedy that the full spectrum can be cut into several sub-ranges (called radiation bands) for each of which a different set of Rosseland coefficients is provided.

The radiative conductivity can then be calculated as a sum of contributions over the radiation bands.

Chapter 9
Immersed Surfaces Technique (IST)

9.1 Introduction

Immersed boundary (IB) techniques were developed in the late nineties for the simulation of flow interacting with solid boundary (see for review Mittal & Iaccarino (2005)), under various formulations. The most widely recognised formulation for instance employs a mixture of Eulerian and Lagrangian variables, where the solid boundary is represented by discrete Lagrangian markers embedding in and exerting forces to the Eulerian fluid domain. The interactions between the Lagrangian markers and the fluid variables are linked by a simple discretised delta function. The IB methods are all based on the direct momentum forcing (penalty approach) on the Eulerian grids, with various forcing formulations on the Lagrangian marker. The forcing should be performed because it ensures the satisfaction of the no-slip boundary condition on the immersed boundary in the intermediate time step. This forcing procedure involves solving a banded linear system of equations whose unknowns consist of the boundary forces on the Lagrangian markers, thus, the order of the unknowns is one-dimensional lower than the fluid variables. Numerical experiments show that the stability limit can in general be altered by the proposed force formulation, in that the second-order accuracy of the adopted numerical scheme is degraded to first order, and in the best conditions to 1.5.

The Immersed Surfaces Technique (IST) has been developed as part of TransAT, although other similar approaches have been developed in parallel. It is implemented in the CMFD code TransAT. The underpinning idea is inspired from the Level Set interface tracking technique for two-phase flows, where free surfaces are described by a hyperbolic convection equation. In the IST, the solid is described as the second phase, with its own thermo-mechanical properties. The technique differs substantially from the IB method discussed above, in that the jump condition at the solid surface is implicitly accounted for, not via direct momentum forcing (using the penalty approach) on the Eulerian grids. It has the major advantage of being able to solve conjugate heat transfer problems, in that conduction inside the body is directly linked to external fluid convection.

A simple example is reported in Figure 9.1, consisting of the flow past a conical shape or a round-shaped structure. The solid is first immersed into a cubical grid covered by a Cartesian mesh. The solid is defined by its external boundaries, which will be described by a sort of level set function denoted by ϕ. The level set is a signed distance to the surface. It is zero at the surface, negative in the fluid and positive in the solid. The treatment of viscous shear at the solid surfaces is handled very much the same way as in all CFD codes, where wall indices are known from a boundary condition file. In the IST these wall cells are identified as those in which ϕ = 0.


PIC

Figure 9.1: Representation of the solid surface within a cubical grid in the Immersed Surfaces Technique. ϕ = 0 at the solid boundaries.


9.2 Equations and Implementation

9.2.1 Immersed Level Set and solid Heaviside functions

The immersed surface is represented on the fluid grid by a Level Set function (ϕs), where ϕs = 0 represents the fluid–solid interface. ϕs is a signed distance function which is positive in the solid phase and negative in the fluid phase. The solid part of the domain is excluded during the solution of the fluid conservation equations through the use of a smooth Heaviside function H(ϕs) which has value 1 in the fluid phase and 0 in the solid phase.

H(ϕs) = 1(1 - tanh (2ϕs))
        2          δsf
(9.1)

where, δsf is the finite solid–fluid interface thickness. This function is referred to as Hf since it indicates the fluid region.

We present the fluid conservation equations for the case where the solid is stationary and the solid velocity is set to zero (uis = 0). The standard Navier–Stokes equations are used for the fluid phase:

 ∂
---(Hf ufj) = 0
∂xj
(9.2)

   f  f f       (         )          f       (          )
∂H--ρ--ui+  -∂-- Hf ρfufuf  = - Hf ∂p--+  -∂-- 2μfHf Sf  + Hf ρfgi - 2μfSf njδ(ϕs)
    ∂t      ∂xj        i j         ∂xi    ∂xj         ij                 ij
(9.3)

the last term in the RHS is a viscous shear at the wall, where nj is the normal to the fluid–solid interface and δ(ϕs) is the Dirac delta function representing the location of the interface. The wall shear itself is modelled as (Beckermann et al.1999),

             (    )
  f  f      f  2ui
2μ Sijnj = μ   δsf
(9.4)

When used in combination with RANS turbulence modelling and wall–functions, the wall shear is calculated using the logarithmic law of the wall.

9.2.2 Solid temperature equation

For an immersed surface TransAT can solve for the temperature distribution inside the solid. This separate equation for solid temperature is termed TSolid for short in the user manual and the graphical user interface.
For TSolid the steady state heat conduction equations are solved inside all solids selected

0 = ∇ ⋅(λ ∇T )+ ωT
(9.5)

where λ is the thermal conductivity and ωT is a heat source term.
Boundary conditions can be of the Dirichlet (constant temperature) or the Neumann (constant heat flux) type. Spatial variation in the boundary condition value is allowed.

Please note: Only steady state temperature distributions are solved inside the solid.

Please note: For non-zero heat flux through the immersed surface, an appropriate heat source term has to be applied inside the solid, in order to balance the heat equation.

Chapter 10
Smart Adaptive Run Parametrisation (SArP) algorithm

10.1 Theory

Evolutionary Genetic Algorithms (EGAs) are a family of computational models inspired by Darwin’s evolution theory. These algorithms liken potential solutions to specific problems chromosome-like structures. Recombination operators are applied to these structures such that critical information is preserved. Despite being often viewed as function optimizers, genetic algorithms are applied to solve a wide range of problems.

A genetic algorithm implementation begins with a population of chromosomes (i.e. sets of parameters) which is expanded from generation to generation. The fitness of each chromosome is evaluated at each generation. The fitness of chromosomes is the determining factor for the generation of the following generations. Chromosomes of the next generation are produced by a reproduction operator.

From one generation to the other, the reproduction operator produces new chromosomes by means of mutation and crossover. The reproduction operator is implemented in such a way that it biased towards the fittest solutions of the population.

As the EGA progresses, the fittest solutions should prevail. The EGA is stopped if the convergence criteria of the EGA or the maximum number of generation are reached. The EGA can also be stopped earlier if the solution obtained from a chromosome meets the targets set for a satisfactory solution. In the case of the EGA implemented in TransAT, satisfactory solutions can be obtained for steady simulation if residuals are below convergence threshold and for unsteady simulations if the physical time to be reached has been overcome.


PIC

Figure 10.1: Overview of the genetic algorithm process


10.2 Genetic Algorithm

The EGA of the SArP module is a controller built around the numerical solver of TransAT with the aim of automatically selecting the numerical parameters to either optimise the convergence or stabilise a simulation. It selects the fittest solution in a finite population of solutions which evolves over different generations accordingly to a natural genetic process-like algorithm.

The only information required to make the EGA progress over generations is the knowledge of the fitness of the individuals to sort them. The fitness of the individuals is obtained by evaluating the so-called fitness function. The fitness function is based on objectives (e.g. target residuals, physical time constraint, etc.).

The problem statement can be expressed in terms of fitness function as follows:

argxm∈aΩx f (x )
(10.1)

where Ω is the set of chromosomes (i.e. the set of numerical parameters) and f is the fitness function.

10.2.1 Structure

The flowcharts in Figure 10.2 and Figure 10.3 show the different stages. Figure 10.2 gives a general overview of the process while Figure 10.3 goes into more details about the selection process.

As shown in Figure 10.3, the initial population of the EGA is generated by altering the original set of parameters i.e. chromosome. New chromosomes are produced at every generation by means of mutating the fittest individual with the aim of enhancing the fitness. In a nutshell, a population made of N individuals, comprises one original chromosome and N-1 mutations of that original chromosome. The mutation rate, probability to carry out a mutation on a chromosome, is set at a high value to guarantee a variety of chromosomes in the population. With high mutation rates, the search space is wider and it reduces the risk of encountering a local optimisation problem.


PIC

Figure 10.2: Genetic Algorithm Process



PIC

Figure 10.3: Genetic algorithm flowchart


Chromosomes are composed of genes. The genes are the individual parameters used to set numerical solvers. These parameters can be of one of the three types:

The ranges of the different genes were set accordingly to the options available in TransAT and from experience.

10.2.2 Triggers and stopping criteria

Triggers

To avoid unnecessarily increasing the computation time, the EGA should only be used when a simulation diverges, crashes or does not converge fast enough. The activation conditions and the steps performed by the EGA vary depending on the type of simulations, stationary or transient.

The EGA has two modes which are depicted in the flowcharts of Figure ?? and Figure 10.5:

In stability mode, the EGA is launched upon occurrence of one of the following conditions depending on the type of simulation:

Similarly to residual convergence in steady cases, the time-step duration decrease is measured by fitting a linear function to the point cloud of time-step duration.


PIC

Figure 10.4: SArP Stability Mode Flowchart



PIC

Figure 10.5: SArP Optimisation Mode Flowchart


Stopping criteria

The EGA stops when one of the stopping criteria is met. The stopping criteria vary accordingly to the nature of the simulation, steady or unsteady. For steady cases, the EGA stops if:

A typical steady simulation fitness evaluation is shown in Figure 10.6. In that example the residuals exit condition is met by the variable represented with red disks and green crosses since their residuals drop below γResmax where γ is the user-defined residual drop factor. The exit condition is however not met by the variable represented with purple diamonds.


PIC

Figure 10.6: Steady Simulation Evaluation Example


For unsteady cases, the EGA stops if:

10.2.3 Genetic Operations

Genetic operations are used to create new chromosomes based on an existing population. The genetic operations consist in applying chromosome crossover followed by mutation in each generation except the first one where the mutation operation is applied directly on the original set of parameters.

Mutation Operation

Mutation is only applied on a random number of new chromosomes. Mutation is then applied on a random set of genes, which size (ratio of number of genes in set to total number of genes in chromosome) is equal to the mutation rate. Thanks to mutation, the parameters space is more thoroughly explored which in turn prevents the algorithm from getting stuck in an area around a local minimum.

There is a mutation operator for each type of parameter. They operate as follows:

Crossover Operation

The implemented crossover operation is based on a one-point crossover. According to the encoding method and the predefined physical rules, the overall consistency of the underlying solution and the validity of the individuals in the population must be maintained. To make sure that genes and combination of genes are within the possible ranges, an extra correction step is added in the breeding process of the EGA. These corrections allow avoiding crashes or unphysical solutions in TransAT. In short, the corrections are akin to eliminating unfit chromosomes a priori. In cases where two individuals have the exact same chromosome, mutation is applied to the latter individual.

10.2.4 Selection Process

The basic part of the selection process is to stochastically select from one generation to create the basis of the next generation. The requirement is that the selection process should make the chances of survival of the fittest individuals greater than weaker ones. This replicates natural selection process in that fitter individuals will tend to have a better probability of survival and will go forward to form the mating pool for the next generation. Weaker individuals are not without a chance. In nature such individuals may have genetic coding that may prove to be useful for future generations.

Fitness function

In each generation, chromosomes are accordingly evaluated with TransAT to measure their fitness. Simulations are run for each chromosome over a window which size corresponds to a certain number of iterations or time-step depending on the simulation type. Residual values and convergence rates are analysed over the window to compute the fitness of the chromosome.

The fitness function indicates the fitness of a chromosome and serves as a metaheuristic to initiate the genetic algorithm. The multi-objective fitness function used to deduce the fitness of chromosomes is defined as a linear combination of simulation runtime efficiency and convergence quality. The distance to the target, quality of convergence and stability are the input variables to the fitness function.

Different indicators determining how good a simulation converges are isolated and weighted in the fitness function. In order to take the different weights of the indicators into account, the indicators are multiplied by weighted coefficients in the fitness function definition.

The stability and the quality of convergence at a point i (i.e. at a specific time-step or iteration depending on the simulation type), is measured with the least squares approximation method which is applied on the window [i - w,i], with w the size of the window.

The fitness function reads:

∀i ∈ [w, N ] fi(di,ci,ei) = - αi|log(di)|- βici - γ|log(ei)|
(10.2)

Where:

For stationary simulations, the fitness function is exactly as defined in Eq. 10.2. For transient simulations, two extra terms are appended to the fitness function definition in Eq. 10.2. One is related to the convergence of the last time step and the other to the duration of time steps.

For steady cases, the following sequence of filtering is performed to select the fittest chromosome:

  1. Select the chromosomes that reached the target residuals if any
  2. Select the chromosomes that do not trigger SArP if the set created in 1 is empty
  3. Select the chromosomes with the lowest runtime out of the set created in 2
  4. Select the chromosome with the highest fitness at the last iteration out of the set created in 3

For unsteady test cases, the selection process is different than with steady test cases. The selection is done sequentially:

  1. Select the chromosomes with converged time steps
  2. Select the chromosomes that do not trigger SArP out of the set created in 1
  3. Select the chromosomes with the highest time step to runtime ratio out of the set created in 2
  4. Select the chromosome with the highest fitness at the last time step out of the set created in 3

Selection Operators

To extract the fittest individuals, different selection strategies from the DEAP module (Fortin et al. (2012)) are implemented:

The elitism strategy always enhances the fitness of the individuals and guarantees that the fittest solution is constantly selected for the next generation. In the selection process, the elitism strategy is applied alongside a selection operator set by the user among the ones above in each generation to add randomness to the algorithm. This hybrid solution combines the advantages of the two employed strategies: favouring the fittest solution and deep exploration of the search space.

10.2.5 Prediction and Fitness approximation

In theory, the fitness of a chromosome is obtained by evaluating the fitness function which means running a simulation. Bearing in mind the number of simulations required to run to evaluate the fitness function (one for each chromosome produced during the EGA process), running the EGA can be computationally costly.

To offset the potential excessive cost of the EGA, estimation models predicting the fitness based on past fitness evaluations implemented in the SArP module can be used. Instead of letting the EGA progress naturally – i.e. by evaluating the fitness with a TransAT run for every chromosome – the rules of evolution are redefined by integrating an associative memory that evaluates the fitness of a chromosome by predicting it based on past evaluations. The full flow chart of the EGA process with prediction enabled is depicted in Figure 10.7.


PIC

Figure 10.7: Genetic algorithm flowchart


The aim is to select more fitting individuals to be eventually evaluated with TransAT. Such a bias in the evolution process gives rise to a self-adapting strategy and decreases the number of calls to the original fitness function which is expensive. Such approaches have been the focus of several academic studies which have shown that they can improve the computational performances.

The proposed solution in the SArP module is based on machine learning techniques such as non-linear regression models and decision trees.

Strategy

The idea, as detailed in the previous diagram, is to launch the EGA in order to feed the machine learning algorithm with its output data. Indeed, the data used to feed the machine learning algorithms comes from a local database that is expanded by runs of the SArP module with fitness prediction enabled. When fitness prediction is enabled, fitness evaluations of chromosomes evaluated during the EGA are added to that local database.

If the database is large enough, it can be used to refine the machine-learning models. When prediction is enabled, the fitness of the chromosomes of a generation is first estimated using the prediction algorithms of machine-learning models. A set of half the population size consisting of the fittest elements is then extracted from the population to be evaluated using TransAT. This way, the number of calls to the TransAT solver is significantly reduced.


PIC

Figure 10.8: Predictive evaluation flowchart


As shown in Figure 10.8, the SArP module executes the following steps when prediction is enabled:

  1. Verify whether data history exists for this particular simulation or if only the chromosomes created during the current EGA run should be relied on
  2. Once the appropriate sample of data is located, proceed to selecting chromosomes that are similar to the current chromosome in terms of state and progress level. The selection is solely applied on data from chromosomes that are in the same phase, progress or state as the current chromosome for which the fitness is to be predicted
  3. Check that the training samples size complies with the learning threshold. Learning thresholds were defined empirically from the study of the different models performances in varying the size of the training sample and analysing the results. As a result, minimum thresholds were fixed for each model in order to get satisfactory prediction results. If the training sample selected after a cross-validation process is not large enough to get satisfactory model performances, the fitness prediction model is not applied and the EGA proceeds to evaluate all the chromosomes
  4. Train the machine learning algorithm on the selected data to generate the prediction model
  5. The fitness of the chromosomes of the current generation are predicted based on the applied model
  6. Return the chromosomes with the best predicted fitness to be passed on to the next generation

Learning Methods

Two types of machine-learning models were implemented: classifiers and regressors.

The former type of model has a binary behaviour. Chromosomes can be either fit or unfit depending on their fitness. These chromosome classes are uniformly distributed within the data from the database. This classification method is an efficient way to eliminate cases that are quickly diverging. Although this model cannot order two solutions in the same class, it remains faster than regressors.

With regressors, regression is used to predict the actual value of the fitness value. It gives more accurate results than with classifiers. The estimation of the fitness function as a continuous function allows the differentiation and ordering of chromosomes. Non-linear regression models are used to handle the complexity of the problems at hand.

The list of implemented algorithms, based on the scikit-learn API (Buitinck et al. (2013)), in the SarP module is as follows:

10.2.6 Speeding up the EGA

The main strategy to speed up the EGA is to reduce the number of evaluations performed during the correction and optimisation process. There are three main options to do this:

Note that all three options can be coupled as thy are independent from one another.

Alternatively, the EGA can be stopped in the middle of the process if a satisfactory solution has been found. This can be done by clicking PIC in the Execute tab of TransATUI or creating the file transat_mb.rti and setting its content to stop_sarp (only one word in the file).

Bibliography

   Beckermann, C., Diepers, H.-J., Steinbach, I., Karma, A. & Tong, X. 1999 Modeling melt convection in phase-field simulations of solidification. J. Comp. Phys. 154, 468–496.

   Buitinck, L., Louppe, G., Blondel, M., Pedregosa, F., Mueller, A., Grisel, O., Niculae, V., Prettenhofer, P., Gramfort, A., Grobler, J., Layton, R., VanderPlas, J., Joly, A., Holt, B. & Varoquaux, G. 2013 API design for machine learning software: experiences from the scikit-learn project. In ECML PKDD Workshop: Languages for Data Mining and Machine Learning, pp. 108–122.

   Cooke, C. H. & Chen, T. 1992 On shock capturing for pure water with general equation of state. Comm. Appl. Num. Meth. 8, 219–233.

   Deb, K., Agrawal, S., Pratap, A. & Meyarivan, T. 2000 A Fast Elitist Non-dominated Sorting Genetic Algorithm for Multi-objective Optimization: NSGA-II, pp. 849–858. Berlin, Heidelberg: Springer Berlin Heidelberg.

   van Doormaal, J. & Raithby, G. D. 1984 Enhancements of the simple method for predicting incompressible fluid flows. Num. Heat Transfer 7, 147–163.

   Fortin, F.-A., De Rainville, F.-M., Gardner, M.-A., Parizeau, M. & Gagné, C. 2012 DEAP: Evolutionary algorithms made easy. Journal of Machine Learning Research 13, 2171–2175.

   Fox, R. 2003 Computational Models for Turbulent Reacting Flows. Cambridge University Press.

   Karki, K. & Patankar, S. 1989 Pressure based calculation procedure for viscous flows at all speeds in arbitrary configurations. AIAA J. 27, 1167–1178.

   Khosla, P. K. & Rubin, S. G. 1974 A diagonally dominant second-order accurate implicit scheme. Computers & Fluids 2, 207–209.

   Launder, B. & Spalding, D. 1974 The numerical computation of turbulent flows. Comp. Meth. Appl. Mech. Eng. 3, 269–289.

   Leonard, B. P. 1979 A stable and accurate convective modelling procedure based on quadratic upstream interpolation. Comp. Meth. Appl. Mech. Eng. 19, 59–98.

   Michelassi, V., Theodoridis, G. & Papanicoulaou, E. 1994 Comparison of time–marching and pressure–correction solvers for compressible flows – part I: low-speed formulation. In ASME 1995 International Mechanical Engineering Congress and Exposition. San Francisco, CA, USA.

   Mittal, R. & Iaccarino, G. 2005 Immersed boundary methods. Ann. Review Fluid Mech. 37, 239–261.

   Patankar, S. V. 1980 Numerical heat transfer and fluid flow. New York, USA: Hemisphere.

   Peters, N. 2000 Turbulent Combustion. Cambridge University Press.

   Poinsot, T. & Veynante, D. 2001 Theoretical and Numerical Combustion. R. T. Edwards Inc.

   R.J. Kee, M.E. Coltrin, P. G. 2003 Chemically reacting flow. John Wiley and Sons.

   Spalding, D. B. 1972 A novel finite difference formulation for differential expressions involving both first and second derivatives. Int. J. Num. Meth. Eng. 4, 551–559.

   Stone, H. L. 1968 Iterative solution of implicit approximations of multidimensional partial differential equations. SIAM J. Num. Anal. 5, 530–558.

   Williams, F. A. 1958 Spray combustion and atomization. Phys. Fluids 1, 541.

   Zhu, J. 1991 A low diffusive and oscillation-free convection scheme. Comm. Appl. Num. Meth. 7, 225–232.

   Zhu, J. 1992 An introduction and guide to the computer program fast-3d. Tech. Rep. 691. Institut for Hydromechanics, University of Karlsruhe, Germany.

   Zitzler, E., Laumanns, M. & Thiele, L. 2001 Spea2: Improving the strength pareto evolutionary algorithm. Tech. Rep..