Acceleration methods for multi-physics compressible ...

Journal of Computational Physics 358 (2018) 201–234

Contents lists available at ScienceDirect

Journal of Computational Physics www.elsevier.com/locate/jcp

Acceleration methods for multi-physics compressible flow Oren Peles ∗ , Eli Turkel Department of Applied Mathematics, Tel Aviv University, Israel

a r t i c l e

i n f o

Article history: Received 16 February 2017 Received in revised form 8 October 2017 Accepted 9 October 2017 Available online 20 December 2017 Keywords: Multi-physics flows Turbulent flows Reactive flows Convergence acceleration

a b s t r a c t In this work we investigate the Runge–Kutta (RK)/Implicit smoother scheme as a convergence accelerator for complex multi-physics flow problems including turbulent, reactive and also two-phase flows. The flows considered are subsonic, transonic and supersonic flows in complex geometries, and also can be either steady or unsteady flows. All of these problems are considered to be a very stiff. We then introduce an acceleration method for the compressible Navier–Stokes equations. We start with the multigrid method for pure subsonic flow, including reactive flows. We then add the Rossow–Swanson–Turkel RK/Implicit smoother that enables performing all these complex flow simulations with a reasonable CFL number. We next discuss the RK/Implicit smoother for time dependent problem and also for low Mach numbers. The preconditioner includes an intrinsic low Mach number treatment inside the smoother operator. We also develop a modified Roe scheme with a corresponding flux Jacobian matrix. We then give the extension of the method for real gas and reactive flow. Reactive flows are governed by a system of inhomogeneous Navier–Stokes equations with very stiff source terms. The extension of the RK/Implicit smoother requires an approximation of the source term Jacobian. The properties of the Jacobian are very important for the stability of the method. We discuss what the chemical physics theory of chemical kinetics tells about the mathematical properties of the Jacobian matrix. We focus on the implication of the Le-Chatelier’s principle on the sign of the diagonal entries of the Jacobian. We present the implementation of the method for turbulent flow. We use a two RANS turbulent model – one equation model – Spalart–Allmaras and a two-equation model – k–ω SST model. The last extension is for two-phase flows with a gas as a main phase and Eulerian representation of a dispersed particles phase (EDP). We present some examples for such flow computations inside a ballistic evaluation rocket motor. The numerical examples in this work include transonic flow about a RAE2822 airfoil, about a M6 Onera wing, NACA0012 airfoil at very low Mach number, two-phase flow inside a Ballistic evaluation motor (BEM), a turbulent reactive shear layer and a time dependent Sod’s tube problem. © 2017 Elsevier Inc. All rights reserved.

1. Introduction Every steady CFD computation has two aspects. One is the accuracy of the final numerical solution which is determined by the details of the spatial (and perhaps physical time) algorithm. The second aspect is the scheme used to reach the steady state (or dual time step). The Navier–Stokes equations form a nonlinear multidimensional system of equations. Fur-

*

Corresponding author. E-mail address: [email protected] (O. Peles).

https://doi.org/10.1016/j.jcp.2017.10.011 0021-9991/© 2017 Elsevier Inc. All rights reserved.

202

O. Peles, E. Turkel / Journal of Computational Physics 358 (2018) 201–234

thermore, depending on the physical situations the system can behave as a mixed type system including elliptic, parabolic and hyperbolic portions. In some cases such as transonic flows, supersonic flow with boundary layers and choked combustion chambers, the mixed types appears in one physical domain. Therefore, it is very difficult to solve the equations directly. Thus, a common way to solve the steady state equations is to add a pseudo-time variable and to march the system in this pseudo-time until the steady state is (approximately) reached. Given a complex three dimensional configuration the time to march the equations in pseudo-time can be quite large. This is especially true if multiple solutions are needed for example, to solve an optimization problem. A time dependent problem using a dual time step procedure also needs many steady state solves. Hence, this paper will concentrate on the efficiency of the scheme to reach the steady state and we later consider the application to the dual time step procedure. The above description implies that the spatial algorithm is given. There exist many spatial schemes ranging from second order finite difference and finite volume schemes to high order schemes frequently based on various Galerkin finite element approaches. In this research we use the second order finite volume scheme popularized by Jameson–Schmidt–Turkel (JST) [1] in 1981. In the original paper the equations were advanced in pseudo-time by a four stage Runge–Kutta algorithm using a variable local time step to accelerate the approach to the steady state. A two dimensional flow around a NACA0012 profile required about 150 iterations to reduce the residual about 4 orders of magnitude when enthalpy damping (which is not appropriate for viscous flows) was added and only a two order reduction without enthalpy reduction. A typical application [2] for an algebraic turbulence model requires about 200 iterations to reduce the residual by three orders of magnitude. In 1984 Jameson and Baker [3] introduced the use of multigrid into CFD computations. For an Euler calculation about an ONERA M6 wing 100 iterations were required to reduce the residual to machine zero and only about 20 cycles for the lift and drag to reach their asymptotic level. In 2001 Jameson and Caughey [4] achieved multigrid efficiency for inviscid flow by combining LU-SGS with multigrid. However, even laminar flow proved difficult for their algorithm. Jameson and Baker [3] also added a residual smoothing technique to allow a larger time step. This can be considered a preconditioning of the fluid equations. One way of accelerating the convergence to the steady state is to introduce a preconditioning [5,6]. A local time step can already be considered a scalar preconditioning. Later Jameson introduced residual smoothing, which is a preconditioning based on the heat equation, to enable a larger time step. Various other algebraic preconditioners have been introduced most with limited improvement. In 2006 Rossow [7] introduced a preconditioner based on an upwind first order approximation to the Navier–Stokes equations which demonstrated a significant improvement. For viscous flow about a RAE2822 airfoil, using an algebraic turbulence model, machine zero was reached after about 150 cycles. Swanson, Turkel and Rossow [8] extended this to three dimensions with further improvements. The physical model we consider is the Reynolds–Averaged Navier–Stokes (RANS) equations for steady state solutions and Unsteady RANS (URANS) for time dependent problems. Unlike Direct Numerical Simulation (DNS), Large-Eddy Simulation (LES) or an algebraic turbulent model such as Baldwin–Lomax model we add to the Navier–Stokes equations, turbulence model equations. We will study the numerical solution of a one equation model (Spalart–Allmaras) and a two-equation model (k–ω SST model). The turbulence phenomena are modeled in the entire spectrum and produce an additional viscosity (called turbulent or eddy viscosity). The turbulence equations are inhomogeneous advection equations with diffusion terms where the advection velocity is the fluid velocity. The system will be solved uncoupled, i.e. we alternately solve the fluid and turbulence equations with the other one frozen. We extend the previous preconditioner [7,8] to both steady state solutions and time dependent flows for complex flows, including turbulent two-phase, multi-species, reacting and real gas flows. We include both external flow, for example around airfoils and wings, and internal flow in a complex geometry including internal combustion chambers. Instead of one continuity equations we then need to solve individual continuity equations for each chemical specie including diffusion terms that result from the unequal spatial concentrations of the species. The internal energy is no longer linear with the temperature but is represented by a polynomial relation with empirical polynomial coefficients taken from thermodynamic databases. The Arrhenius chemical kinetics model is used for the chemical reactions. The chemical reactions add source terms to the continuity equations of the species and to the energy equation. The last extension of the physical system that we consider is the coupling to an Eulerian representation of dispersed phase flow (EDP – Eulerian Dispersed Phase). This model represents the spreading of solid particles or liquid droplets in the computational domain. The particles have their own spatial density and advection velocity. The governing equations for the motion of the dispersed phase are a continuity equation together with momentum equations. The momentum equation is controlled by the drag force between the particles and the gas. The drag force acts on both the dispersed phase and the main gas phase. Therefore, the same force with negative sign is added to the fluid momentum equations. Some models include mutual interactions between the particles and the turbulence. The particles produce turbulence and their motion is influenced by the eddy viscosity. Stiffness is a wide range of space and time scales in a physical model and mathematical formulation and it causes difficulties and challenges for the numerical modelling. The stiffness in the system is composed of physical stiffness can be due to wide range of velocities (from Mach number zero to supersonic flow). A recent study of the accuracy of a similar low Mach code is given by Langer [9]. Other causes of stiffness include large scales of reaction rates of the species, and strong source terms relative to the advection and diffusion. In many problems, including boundary layers, the numerical grid is characterized by a large grid aspect ratio that causes numerical stiffness. As previously mentioned the emphasis of this study is the iterative method used to march the equations that involve the various physical models just described. The iterative numerical method should yield an accurate solution in both time and


203

space and do it fast in terms of iterations and CPU time while being robust. The numerical method we develop and study is based on the explicit Runge–Kutta (RK) scheme with an implicit smoothing, marching with a pseudo time stepping for a steady state and using a dual time stepping approach for unsteady problems. Implicit smoothers are usually more stable (in many cases unconditionally stable) but too large time steps can lead to non-physical results. Explicit solvers are only conditionally stable. We study the interaction of the RK/Implicit smoothing method with all of the above mentioned physical phenomena and also with a range of numerical methods including spatial difference schemes – scalar and matrix dissipation model, Roe scheme and flux splitting method (we use AUSM family scheme). We also study the interaction between the smoother and the MUSCL approach for second and third order accurate solutions and the interaction with other acceleration techniques such as low Mach preconditioning and multigrid methods. In summary, in this study we extend the previous work of Rossow [7] and Swanson, Turkel and Rossow [8] for acceleration the convergence of the Navier–Stokes equations for complex flows to a steady state and use with a dual time step. The extensions include one and two equation turbulence models, real gas and reactive materials and multi-phase problems. We also include improvements for the dual time step time dependent problem and also the application to low Mach flows. In the CFD Vision study [10] they note “Increasingly, multicomponent and multiphysics simulations are performed during the design cycle.” They add “Furthermore, research should continue into methods to accelerate the calculation of chemical kinetics so that the CFD solution progression is not limited by these ordinary differential equations.” Similarly in [11] the authors write “CFD still faces several challenges that need to be addressed. The turnaround time associated with CFD is one of the factors limiting the use of CFD in design and creation of databases and also in multi-disciplinary applications.” In this paper we take a step towards increasing the efficiency for such problems using fast numerical solutions of multiphysics problems. 2. Flow equations 2.1. Navier–Stokes equations We considered the 3D compressible, non-homogeneous Navier–Stokes equations in conservative form for n chemical species

∂Q + ∇( F − F V ) = S CR − S f ∂t S CR is the vector of source terms due to chemical reaction and S f is the vector of source terms from the interaction with the disperse phase. The conservative variables are,

Q = {ρ1

ρ2 . . . ρn ρ u ρ v ρ w E }T .

ρi are the partial densities of the species, the total density ρ is the sum over all partial densities, ρ = gas the energy is related to the pressure and to the velocity by

E=

P

(γ − 1)

n

i =1

ρi . For an ideal

1 + ρ u2 + v 2 + w 2 2

while for a real gas we have

E=

N

1

ρk ek ( T ) + ρ u 2 + v 2 + w 2 .

k =1

2

ek ( T ) are polynomials in temperature with empirical coefficients from chemical databases. F is the inviscid fluxes vector, and F v is the viscous flux, given by

⎛

ρ1 V ⎜ .. ⎜ . ⎜ ⎜ ρ ⎜ NV F =⎜ ⎜ ρ u V + xˆ P ⎜ ⎜ ρ v V + yˆ P ⎜ ⎝ ρ w V + zˆ P ( E + P ) V

⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠

⎛

⎞

⎜ ⎜ ⎜ ⎜ ⎜ ⎜ Fv = ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝

−ρ1 υ1 .. . −ρ N υ N τx τy τz −q + τx u + τy v + τz w −

N k =1

ρk υk hk

⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟. ⎟ ⎟ ⎟

⎟ ⎟ ⎠

204


υ j are the diffusion velocities of the species j, given by μ ρ νj = − ∇ i .

ρ

sc

μ is the viscosity and Sc is the Schmidt number (the ratio of momentum diffusivity and mass diffusivity). In the case of μ μ μ turbulence, Sc is replaced by Scl + Sct where the subscripts l and t refer to the laminar and turbulent viscosities and Schmidt t l numbers respectively. h j are the specific enthalpy per unit mass of the species j, q = −k∇ T is the heat flux and τx , τy and τz are the stress vectors. In addition, the equation of state is P=

ρ RT W

= RT

n

ωi

i =1

ωi are the molar concentrations with the mean molecular weight W −1 =

n

i =1

ρi 1 ρ wi .

2.2. Turbulent models 2.2.1. k–ω SST model The k–ω SST turbulence model is a 2-equation model in which two variables, k and ω , need to be solved. k is the turbulent kinetic energy and ω is the dissipation of the turbulence. The additional turbulence equations in non-conservative form are [12]:

∂ ∂k ∂k ∗ (ν + σk νT ) + u · ∇ k = P k − β kω + ∂t ∂xj ∂xj ∂ ∂ω ∂ω 2 2 (ν + σω νT ) + C D /ρ + u · ∇ ω = α S − β ω + ∂t ∂xj ∂xj S is the vorticity vector magnitude and the eddy viscosity is computed from

(2.1)

νT =

a1 k max(a1 ω, S F 2 )

and C D is the cross-diffusion 1 ∂k ∂ ω term: C D = 2(1 − F 1 )σω2 ρ ω . ∂x j ∂x j The boundary condition for k and ω are extrapolation for outflow, k = k∞ and ω = ω∞ at inflow condition. The wall boundary condition is k = 0 (no kinetic energy) but ω also needs to be given. In order to match the behavior of the turbulence close to the wall, the value of ω on the wall is dependent on the distance of the first cell from the wall, the standard model gives:

ωwall =

60ν

β1 y 21

(2.2)

y 1 is the distance at the first grid cell above the surface. 2.2.2. Spalart–Allmaras (SA) model The SA model is a one equation model, in which one variable – ν˜ is computed –

2 ∂ ν˜ 1 C b1 ν˜ + u · ∇ ν˜ = C b1 [1 − f t2 ] S˜ ν˜ + {Viscous terms} − C ω1 f ω − 2 f t2 ∂t σ d κ ˜ ν S˜ = S + 2 2 f v2 .

κ d

(2.3)

S is the magnitude of the rotation tensor, d is the wall distance. The turbulence viscosity is given by μt = μ f v1 , f v1 = f v1 (ν˜ /ν ). All the other quantities in the model can be found in the literature [13]. The boundary conditions are ν˜ = 0 on viscous walls and ν˜ = ν˜ 0 for the inflow boundary condition and ∂n ν˜ = 0 for the outflow condition. The results from the SA turbulence model seem to be relatively sensitive to the free stream and initial values. For fully turbulent flow, Allmaras et al. [14] suggested using 3ν < ν˜ < 5ν which is a very narrow range of possible values. 2.3. Chemistry source term The chemical reaction source terms of the inhomogeneous Navier–Stokes equations are:

S = (ρ˙1 , ρ˙2 , ρ˙3 , . . . , 0, 0) T .


205

We introduce the primitive variables q = (ωi , u , T ) where ωi are the molar concentrations. We now define a different set of primitive variables U = (ρi , u , P ). The Jacobian matrix from q-primitive variables to the conservative variables is

⎛

⎜ ⎜ ⎜ dQ = ⎜ ⎜ ⎜ ⎝

w1 0 0 u w1

1

w 1 e1 (T ) +

2

u2

0 w2 0 u w2

w 2 e2 (T ) +

1 2

u2

0 0 w3 u w3

w 3 e3 (T ) +

0 0 0 1 2

u2

ρ

⎞

0 0 0 0

ρu ρcv

⎟ ⎟ ⎟ ⎟ dq = J dq ⎟ ⎟ ⎠

So, in conservative variables the source terms satisfy:

⎛

⎛

⎞

ρ1 ⎜ ⎜ ρ ⎜ d ⎜ 2⎟ ⎟ ⎜ ρ = ⎜ 3⎟ ⎜ ⎜ dt ⎝ ρu ⎠ ⎜ ⎝ E

⎞

⎛ ⎞ ⎟ ω˙ 1 ⎟ ⎜ ω˙ 2 ⎟ ⎟⎜ ⎟ ⎟ ⎜ ω˙ 3 ⎟ ⎟⎜ ⎟ ρ ⎟⎝ 0 ⎠ ⎠ 1 1 1 T˙ w 2 e2 (T ) + u 2 w 3 e3 (T ) + u 2 ρu ρcv w 1 e1 (T ) + u 2 w1 0 0 u w1

0 w2 0 u w2

2

⎛

⎞ ˙1 w1ω ⎜ w 2 ω˙ 2 ⎟ ⎜ ⎟ = ⎜ w 3 ω˙ 3 ⎟ ⎝ 0 ⎠

0 0 w3 u w3

2

0 0 0

0 0 0 0

2

0 Each reaction of the N r reactions in the chemical model is described by the reaction equation:

υik χk

⇔

υik χk

υ and υ are the forward and backward stochiometric coefficients respectively. The Arrhenius model provides the chemical ˙ 1 , ω˙ 2 , ω˙ 3 , 0, T˙ ). These are given by: source terms for (ω

ω˙ k =

Nr

υik qi υik ≡ υik − υik

(2.4)

i =1

T˙ = −

1

ρcv

˙k ek ( T ) w k ω

(2.5)

k

The progress variables, q i are given by

qi = αmi ωim K f ,i ( T )

υ ik

ωk − K r ,i ( T )

k

υ ik

ωk

(2.6)

k

where

k f ,i = A i T βi e − E i /RT

(2.7a)

kr ,i = k f ,i /kc ,i .

(2.7b)

A i , βi and E i are the Arrhenius constants where A i is the rate constant, βi is the temperature exponent and E i is the activation energy. ωim is the molar concentration of a third body that participates in the reaction in some dissociation reactions and αim is the third-body efficiency (αim = 1 for reactions without a third-body). In addition –

k c ,i = K p

K p = exp

S =

k∈ K υik i

P atm RT

S R

−

υik S k

k∈ K i

H =

k∈ K i

υik H k

H RT

206


2.4. Eulerian Dispersed Phase (EDP) In many engineering applications two phase flow occurs. There are two approaches to solve two phase flow involving a gas phase and a dispersed phase. The first one is the Eulerian–Lagrangian approach [15]. In this approach, the particles are treated as a point mass and their trajectories are calculated from the Newton equation of motion. The second approach is the Eulerian–Eulerian [16]. In this approach the dispersed phase is assumed to be continuous. The main advantage of the Eulerian approach is that it is simplifies the coupling between the dispersed and the gas phase and post-processing calculation as mass fluxes and accumulation rate on the model boundary such as walls. Unfortunately, there are also some disadvantages of the model. The main one is that the continuous assumption is not always true. The second one is that some physical phenomena cannot easily be simulated. For example, cross-flow (trajectories can intersect each other) and non-elastic collision with walls. Numerically, the problem of two phase flow is a stiff problem. The reason is the wide range of scales involved in the physical phenomena. In the current work we developed a numerical algorithm, based on the RK/implicit smoother method in order to simulate EDP systems. The physical model consists of a two-phase flow in which the main phase is a gas phase, governing by the Navier–Stokes equations. The second phase is a dispersed phase, including a mixture of liquid and solid particles. The solid particles can change their state, melt and freeze. This phase is governing by a “fluidized” set of equations for motion and heat transfer. We solve for the variables of the dispersed phase {ρ p , v p , T p , χ p } T . ρ p , v p , T p are the particles density, velocity vector and temperature respectively and χ p is the solid fraction, 0 ≤ χ p ≤ 1. When χ p = 1the particles are totally solid and when χ p = 0 they are totally liquid. The governing equations for the dispersed phase density and velocity are: Conservation of mass:

∂ ρp + ∇(ρ p v p ) = 0 ∂t Conservation of momentum:

∂ v p ,i fD + v p · ∇ v p ,i = · ( v g ,i − v p ,i ); ∂t τu

i = (x, y , z)

For the temperature and solid fraction, define the heat flux between the particles and the gas phase:

Q˙ = m p C p , p

fN

τT

(T f − T p )

The thermal relaxation time is given by

3

C p, p

2

C p, f

τT = pr

τu

p r is the Prandtl number, C p , p and C p , f are the heat capacities of the dispersed and gas phases respectively. The Stokes relaxation time is

τu =

ρ˜ p d2p 18μ

d p is the particle’s diameter, ρ˜ p is the particle’s bulk (material) density and for the temperature (conservation of energy) is given by

∂Tp + v p · ∇ T p = ∂t

Q˙ m p C p, p

T = T melt

0

else

μ is the laminar viscosity. Finally, the equation

T melt is the melting temperature. The equation for the solid fraction is given by

∂ χp + v p · ∇ χ p = ∂t

˙

− m pQLm

T = T melt

0

else

L m is the latent heat for melting. 3. Numerical scheme The basic scheme we used is a second order finite volume based on the work of Jameson–Schmidt–Turkel (JST) [1]. In order to stabilize the central differences and allow convergence to a steady state together with reducing oscillations near shocks an artificial dissipation term is added that consists of a combination of second and fourth differences in each direction. The original paper used scalar coefficients to control the amount of dissipation. Later this was improved [2] to


207

matrix valued coefficients. We shall also consider the AUSM schemes [17–19]. The advantage of the use of AUSM schemes (as well as other flux splitting method such as CUSP) is that it is not required the transformation to primitive variables. To march the scheme to a steady state an explicit Runge–Kutta method is used. To accelerate the convergence to a steady state a local time step is used based on the local stability condition which is mesh dependent. In this work the residual smoothing is replaced by an improved preconditioning given by the implicit first order upwind smoother described by Rossow [7] and later Swanson, Turkel and Rossow [8]. In addition multigrid is used to further accelerate the convergence rate. 4. Acceleration methods 4.1. Multigrid The multigrid method is a multi-level numerical method that can accelerate the convergence of partial differential equations; especially elliptic problems such as the Laplace/Poisson equations (see [20,21,3,22]). For M < 1, the steady state Euler equation is an elliptic system and so the multigrid method is very helpful. Multigrid will be a basic part of all the schemes to be discussed. 4.2. Rossow–Swanson–Turkel RK/implicit smoother method We review the concepts of the Rossow–Swanson–Turkel RK/implicit smoother. JST [1] introduced the basics of the finitevolume together with the Runge–Kutta time marching scheme. The implicit smoother was introduced for inviscid, viscous laminar and an algebraic turbulence model flow by Rossow [7] and Swanson et al. [8]. It was extended to turbulent flow with a k–ω SST model by [23], turbulent and reactive flow in [24] and for the Spalart–Allmaras turbulence model in [25]. The method enables not only faster convergence but also a more stable numerical scheme. In fact, in some cases the explicit scheme didn’t converge and we needed the RK/Implicit scheme. Consider a general system in conservation form

∂Q + ∇ F = S ( Q ). ∂t Applying the Gauss theorem for an arbitrary control volume (computational grid cell) and discretizing in time yields the semi-discrete scheme

Q 1 =S− t V

ds = S − F ·n

1 V

S

ds F ·n

(4.1)

all faces

We use a low-storage RK scheme for time marching

Q k+1 = Q 0 − αk t R k

k = 0... p

αk is the RK coefficient of the k-th step and the residual of the step is given by Rk = Sk −

1 V

ds. Fk · n

all faces

For accelerating the calculations using large CFL numbers and in order to make the computations of the stiff problem more robust, the residual is replaced by a spatially smoothed residual R˜ k . Following Rossow [7], Swanson et al. [8] we start with the spatially discretized equation (4.1) and linearize F and S in time to obtain

I+

t V

ds − t A ·n

all faces

∂S ∂Q

Q˜ = R k

(4.2)

A ≡ ∂∂FQ·n is the flux Jacobian and ∂∂QS is the source term Jacobian. Transforming the system to a set of primitive variables,

we can write the flux Jacobian matrix A as a sum of positive and negative matrixes A = A + + A − where A ± = 1/2( A ± | A |). Since A − Q represent the flux associated with waves that enter the cell from the outside and A + Q represent the flux associated with waves that leave the cell (see [7]), we move the A − matrix to the RHS of equation (4.2) multiply by Q˜ of the neighbor cells (NB) and leave the A + matrix on the LHS multiplied by Q˜ of the local cell. We then get the implicit smoothing operator

I +ε

t V

all faces

dS − t A +n

∂S ∂Q

Q˜ local = R k − ε

t V

all faces

dS Q˜ N B A −n

(4.3)

208


In (4.3) ε is a relaxation parameter and might be scheme dependent. For small ε , the convergence is better but the algorithm can become unstable. ε > 1 is an over-relaxation and one can use it for very stiff problems for example reactive flow and two-phase flows. This is an implicit linear equation for the smoothed residual Q˜ . It is solved iteratively by a symmetric Gauss–Seidel method or a Red–Black (the red black method is easy to parallelize for shared memory architectures). One performs only a few GS iterations for smoothing, typically two iterations for regular (N–S) problems or 3–4 iterations for stiffer problems, reactive flows for example. 4.2.1. Entropy fix In the smoothing operator, for very low speed flow near stagnation points, the Mach number must be bounded away from zero to prevent zero dissipation. The velocity at the cell interface is chosen as qd = M d · c · dfak where n M d = max( M , dcut ); dcut = max(c 1 min( f AR , 1), c 2 ) and f AR depends on the ratio of the eigenvalues: f AR = ( λλnr )1/2 . λn is the eigenvalue in the cell face direction and λr is the eigenvalue in the tangent direction for two dimensions. For the three dimensional case we use the harmonic average of the two tangent directions. The constants c 1 and c 2 may be scheme dependent.

4.2.2. Multi-block communication strategy The described methods can be parallelized using a multi-block approach. The method assumes that there are sufficient points in each subdomain to describe the physics within that subdomain. For multi-block parallelism, especially for a large number of small blocks, we can exchange information across blocks interfaces not at every RK stage but in the end of the last RK iteration. This saves lot of CPU time with almost no effect on the convergence rate and results. 4.2.3. Coupling types For two-phase turbulent flow, in which the main flow is a real gas composed of n chemical species there are a total of n + 4 + 2 + 6 variables. This implies that the flux and source terms Jacobian matrixes are (n + 12) × (n + 12) matrices. For the source term Jacobian matrix, it is difficult to derive the relations between the variables of the three main systems (main phase, dispersed phase and turbulence model). Computing this matrix can also be very inefficient in computational terms. This Jacobian matrix takes the form:

⎛ ⎜ ⎜ ⎜ ∂S ⎜ =⎜ ⎜ ∂Q ⎜ ⎝

∂ S gas.phase ∂ Q gas.phase ∂ S gas.phase ∂ Q gas.phase ∂ S turb ∂ Q gas.phase

∂ S gas.phase ∂ Q dis.phase ∂ S dis.phase ∂ Q dis.phase 0

⎞ ∂ S gas.phase ∂ Q turb ⎟ ⎟ ⎟ ⎟ 0 ⎟ ⎟ ⎟ ∂ S turb ⎠ ∂ Q turb

For the implicit smoothing of the problems considered in this work, it was found to be reasonable to approximate the Jacobian by a simple matrix. There are several approaches to these approximations. In the current work we assume that the system is loosely coupled. Thus, only the main block matrices are taken into account. Thus, ∂∂QS simplifies to:

⎞

⎛

∂ S gas.phase ⎜ ∂Q gas.phase ⎜ ⎜ ∂S ⎜ ≈⎜ 0 ⎜ ∂Q ⎜ ⎝ 0

0

0

∂ S dis.phase ∂ Q dis.phase

0

0

∂ S turb ∂ Q turb

⎟ ⎟ ⎟ ⎟ ⎟. ⎟ ⎟ ⎠

For two phase flows, the derivative of the gas-phase source terms with respect to the dispersed phase variables and the derivatives of the dispersed phase source terms with respect to the gas variables shouldn’t be ignored. A more reasonable choice is to approximate the Jacobian as

⎛ ⎜ ⎜ ⎜ ∂S ⎜ ≈⎜ ⎜ ∂Q ⎜ ⎝

∂ S gas.phase ∂ Q gas.phase ∂ S gas.phase ∂ Q gas.phase

∂ S gas.phase ∂ Q dis.phase ∂ S dis.phase ∂ Q dis.phase

0

0

⎞ 0 0

∂ S turb ∂ Q turb

⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠

We do not currently have a conclusion about the advantage or disadvantage of the two options.


209

4.3. Jacobian approximation for the RK/implicit smoother In primitive variables we use only some of the Jacobian entries:

⎛ ∂ ω˙ 1 ⎜ ∂ ω1 ⎜ ⎜ ⎜ 0 ⎜ ⎜ R J {ω , T , u } = ⎜ ⎜ 0 ⎜ ⎜ ⎜ ⎜ 0 ⎝ 0

0

0

0

∂ ω˙ 2 ∂ ω2

0

0

0

∂ ω˙ 3 ∂ ω3

0

0

0

0

0

0

0

∂ ω˙ 1 ∂T ∂ ω˙ 2 ∂T ∂ ω˙ 3 ∂T

⎞

⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ 0 ⎟ ⎠ ∂ T˙ ∂T

As previously discussed, most of the entries are expected to be non-positive. Examples will be described in the coming sections. 4.4. RK/implicit smoother for real gas We will now define the extension for the RK/Implicit smoother for multi-components, real gas problems. We use the , E ), and the set of primitive variables U = (ρi , u , P ). The Jacobian matrix transforming conservative variables Q = (ρi , ρ u the primitive variables δ U to the conservative variables δ Q is given by

⎛

1 /ρ

0

0

⎜ 0 1 /ρ 0 0 ⎜ ⎜ ∂U 0 0 1 /ρ 0 ⎜ =⎜ ⎜ ∂Q ⎜ u (1 − γ ) v (1 − γ ) w (1 − γ ) (γ − 1) ⎝ 0

0

0

⎞

− u /ρ

0

⎟ ⎟ ⎟ − w /ρ ⎟ ⎟. ⎟ U2 (γ − 1) ⎟ ⎠ 2 − v /ρ

0

1

The inverse transformation Jacobian is ∂∂ Q . The flux Jacobian in primitive variables is given by U

⎛

ρ uk + nx ρ u nyρu nz ρ u ⎜ nx ρ v ρ uk + n y ρ v nz ρ v ⎜ ⎜ ⎜ nx ρ w n ρ w ρ u y k + nz ρ w ⎜

∂ F = ∂U ⎜ ⎜ ⎜ nx h + ρ uuk ⎝

n y h + ρ vuk

n z h + ρ wuk

n y ρ yi

nz ρ y i

nx ρ y i | u|

nx

⎞

uuk

ny

⎟ ⎟ ⎟ ⎟ wuk 2 ⎟ ⎟ | u| ⎟ uk − Πi =1.. N ⎟ ⎠ 2 vuk

nz uk γ

(γ − 1) 0

uk I NxN

2

+ ρp + e ) and Πi = (γ −RT1) w i − e i (= 0 for perfect gas). U F The flux Jacobian in conservative variables is given by A c = ∂∂W = ∂∂ UF ∂∂W . The flux Jacobian in the primitive variables is

The total energy is h = ρ (

2

⎛

uk ⎜ ∂U ∂ U ∂ F ⎜ 0 ∂W Pn = = =⎜ 0 Ac ∂W ∂U ∂ W ∂ U ⎝ n pγ x nx ρ y i

0 uk 0 n y pγ n y ρ yi

0 0 uk nz pγ nz ρ y i

n x /ρ n y /ρ n z /ρ uk 0

0 0 0 0 uk I NxN

⎞ ⎟ ⎟ ⎟. ⎠

The absolute value of P n is defined by | P n | ≡ T |Λ| T −1 where T is the matrix of eigenvectors of P n and |Λ| is its absolute values of the eigenvalues along a diagonal matrix, hence,

⎛ ⎜ ⎜ |Pn| = ⎜ ⎝

|uk | + a1 cnxnx a1 cnxn y a1 cnxn z M 0nx p γ M 0nx ρ y i

a1 cnxn y |uk | + a1 cn y n y a1 cn z n y M 0n y p γ M 0n y ρ y i

a1 cnxn z a1 cn y n z |uk | + a1 cn z n z M 0nz p γ M 0nz ρ y i

M 0 n x /ρ M 0 n y /ρ M 0 n z /ρ |u k | + a 1 c a1 y i /c

0 0 0 0 |uk | I NxN

⎞ ⎟ ⎟ ⎟ ⎠

210


where

M 0 = min | M k |, 1 sgn( M k ) M k = uk /c

· n uk = u a1 = 1 − | M 0 |. The matrices P n± = 12 ( P n ± | P n |) are given by

⎛

a1 cnxnx a ⎜ 1 cn x n y 1⎜ 1 a1 cnxn z P n± = u k ± |u k | + ⎜ 2 2⎝ (1 ± M 0 )nx p γ (1 ± M 0 )nx ρ y i

a1 cnxn y a1 cn y n y a1 cn z n y (1 ± M 0 )n y p γ (1 ± M 0 )n y ρ y i

a1 cnxn z a1 cn y n z a1 cn z n z (1 ± M 0 )n z p γ (1 ± M 0 )n z ρ y i

⎞ (1 ± M 0 )nx /ρ 0 (1 ± M 0 )n y /ρ 0 ⎟ ⎟ (1 ± M 0 )n z /ρ 0 ⎟ . ⎠ a c 0 1

a1 y i /c

0

4.5. General systems and unstructured grids The above subsections give specific examples of the generalization of the NS equations. If one introduces another model one would need to modify the smoother using the following considerations: The matrixes P n and | P n | should be modified according to the corrected eigenvalues of the system. If the generalization of the system contains a physical source term, a low-cost approximation of the source term Jacobian should be used in the right-hand side of the smoother operator (see equation (4.3)). An appropriate approximation will depend on physical considerations for each case separately. It is recommended that the diagonal entries of the Jacobian contain the negative parts of the source terms (sink terms). For an unstructured grid, equation (4.2) defines a linear system for the smoothed residuals. We can consider this system as a general system of the form Ax = b. The order of the rows of A need not have a geometric significance. We then use any iteration method to solve this system (Jacobi, Gauss–Seidel, Red–Black GS). For some of these iterative methods the convergence depends of the ordering of the rows. For a general system, there is no clear choice of ordering the rows (hence the internal ordering of the cells) to increase the rate of convergence. [26] discusses the implementation of the RK/Implicit smoother method for unstructured grids and suggests that there is an advantage to ordering the system in the direction of strong coupling (for example in the direction perpendicular to the airfoils because of the grid aspect ratio). The difficulty with this suggestion is that for general cases, such as internal flows, there is no “natural” direction to the flow. 5. Low Mach number and time dependent flows 5.1. Low Mach properties of Euler/N–S equations The Navier–Stokes equation for a steady-state inviscid flow change type, at a Mach number equal to one, from elliptic for subsonic to parabolic for supersonic flows. The pseudo-time stepping approach is used to overcome this problem. As the Mach number becomes small, one expects, under suitable conditions, the compressible flow to behave like an incompressible flow. A main difficulty in terms of convergence rate is the incompressible limit where the Mach number approaches to zero. The numerical issue is to solve almost incompressible flow with a compressible solver. One solution to the low Mach difficulty is the preconditioning method [5,6,27–29]. We multiply the time derivative with respect to the pseudo-time by any non-singular matrix P without affecting the solution in the steady state. The difficulty in the convergence is the spread in the eigenvalues of A. We wish to choose P so that the eigenvalues of P A have less of a spread. The preconditioning matrix is not unique. There are several matrices in the literature including [5,30–32]. Since the pseudo-time serves only to march the system toward the steady-state solution, the accuracy in pseudo time is not important. Therefore, every convergence acceleration technique which does not affect the final spatial accuracy can be used. An example of such an acceleration technique is the use of a local time step, where the time step is chosen separately for each cell, and not globally. Until now we have discussed acceleration to a steady state. However, a second difficulty is caused by the lack of accuracy in the steady state [6,5,27]. This is caused by the artificial dissipation. The Jacobian flux matrix in the numerical scheme (such as Roe scheme, Matrix dissipation scheme etc.) should be | P A | instead of | A |. Numerically, the modified equation is then

P −1 w t + A w x = − P −1 | P A | w x x x

(5.1)

Turkel et al. [27] and Darmofal and van Leer [32] show that if the primitive variables dw = (dp /ρ a, du , dv , dS ) are chosen, then in the low Mach limit, where dS = dp − a2 dρ is the entropy

T

A w x = O (1/ M ), O (1), O (1), O (1)


211

but

T | A | w xx x = O (1), O (1/ M ), O (1/ M ), O (1) . Without preconditioning ( P = I ), the left and the right sides of equation (5.1) do not agree at the limit of M → 0. In order to match the left and right sides of the equation P −1 and P −1 | P A | should behave as

⎛

O (1/ M 2 ) ⎜ O (1/ M ) ⎜ ⎝ O (1/ M ) O (1/ M )

O (1/ M ) O (1) O (1) O (1)

O (1/ M ) O (1) O (1) O (1)

⎞

O (1/ M ) O (1) ⎟ ⎟ O (1) ⎠ O (1)

We see that P −1 | P A | is responsible for the accuracy maintenance on the incompressible regime. We now can use R −1 to multiply the time derivative w t . This operator will serve as the convergence accelerator (such as RK implicit smoother)

R −1 w t + A w x = − P −1 | P A | w x x x

(5.2)

5.2. Roe-c scheme w/MUSCL Rossow [7], suggested modifying the A + and A − so that it will behave properly for low Mach numbers. He replaced the speed of sound c by a modified c :

c =

α=

α 2 q2 + M r2 c 2 1 2

1 − M r2

2 ,1 M r2 = min max M 2 , kM ∞

In this work the same idea has been extended. We use first order “upwind” scheme with | A |(c ) instead of P −1 | P A | where A ± = 12 ( A ± | A |) and second order is achieving using the MUSCL approach. A similar approach is used for the Jacobian in the RK/implicit smoother (4.3). We have found that the Roe-c scheme is more robust than the standard preconditioning (5.2). In particular, we could achieve lower Mach numbers using c , see also [9]. 5.3. Time dependent problems When the physical time scale is large relative to the CFL time then the time dependent problem can be solved using the dual time steps technique (see in Jameson [33]). The physical system we solve, in conservation form, is:

∂Q + ∇ F = S ( Q ). ∂t We add an artificial time so the equation becomes

∂Q ∂Q + ∇ F = S(Q ) − . ∂τ ∂t Here

τ is a pseudo-time and t is the physical time. The discrete version is Q Q 1 1 Q ds = S − ds − =S− − F ·n F ·n τ t V V t S

(5.3)

all faces

Q /t is a forward discrete approximation to the physical time: ct Q n+1 − 2Q n + 12 Q n−1 ∂Q Q ≈ = ; ∂t τ t

ct = 3/2 for a second order method.

We approximate (4.3) using preconditioning and a dual time steps algorithm:

I + ct αk

τ ct ωck − F (ωn , . . .) τ + ct αk P ωck+1 = ωc0 − αk τ P R kQ + P ωck t t t

(5.4)

P is the low Mach preconditioning matrix and R Q is the residual in the conservative variables. For low Mach flow we can still choose P = I since the low Mach is treated with by using c . Then, equation (5.4) can be then be rearranged to

212


Fig. 5.1. Convergence history for series of CFL numbers for three physical time steps.

obtain k +1 c

ω

0 c

= ω − αk τ

{ R kQ +

ct ωc0 − F (ωn ,...) } t

τ (1 + ct αk t )

.

Define

R k = R kQ +

ct ωc0 − F (ωn , . . .)

(5.5a)

t

For time dependent problems the residual of the k-th step is:

R kQ = S k −

1 V

ds − Fk · n

all faces

Q t

(5.5b)

The final RK/Implicit smoother is then

I +ε

t V

all faces

∂S τ A ndS − t + ct αk I ∂Q t +

t − dS Q˜ N B Q˜ local = R k − ε A n V

(5.6)

all faces

We note that this formulation differs from the standard dual time step, e.g. Turkel and Vatsa [34], in two different but connected ways. First the lower order term on the left-hand side of eq. (5.6) is multiplied by αk . Second in (5.5a) ωc0 appears and not ωck i.e.:

R k = R kQ + and

I +ε

ct ωck − F (ωn , . . .)

t

t V

all faces

dS − t A +n

∂S τ + ct I ∂Q t

Q˜ local = R k − ε

t V

dS Q˜ N B . A −n

all faces

To validate the dual time stepping scheme, included within the RK/implicit smoother, we consider the solution of the Riemann problem (Sod’s shock tube [35]). The initial conditions are p = 1 PSI and T = 231.11 K in the right side of the tube and p = 10 PSI and T = 288.89 K in the left side with Mach number zero at both sides. For three physical time steps, the convergence for a series of CFL numbers is presented in Fig. 5.1. Improved convergence is obtained as the CFL number increases until reaching an asymptotically convergence. Comparison of the results to the analytic solution is presented in Fig. 5.2. Computations demonstrate that the new approach is more robust especially at very low Mach numbers. Fig. 5.3 shows a comparison between the conversion histories of the two methods for the Shu–Osher problem (see details in [36]). For ten pseudo-time iterations with CFL = 64, the results are practically the same but the new method converged much better. Furthermore, it reaches the same error after about five iterations so it can save almost half of the CPU time.


213

Fig. 5.2. Sod’s tube results – Analytical solution vs. numerical simulation. Pressure (top) and density (bottom).

Fig. 5.3. Convergence history comparison between the standard dual time-stepping method (A) and the new dual-time stepping method (B) for the onedimension Shu–Osher problem. The CFL used was 64.

214


6. Real gas and reactive flow 6.1. Numerical consideration for stiffness Reactive flow is multi-components (chemical species) flow, in which chemical reactions occurs within the flow. Physically, the reactions are modeled with source terms that represent the change in time of the several species. Numerically, these source terms make the problem a very stiff one because of the high diversity of time scales between individual reactions and between the chemical time scale to the advection and diffusion time scales. Another type of chemical reactions that will not be discussed here is surface reactions. In [10] they state “In the area of turbulent reactive flows, investment needs to continue toward the development of a validated predictive multiscale combustion modeling capability to optimize design and operation of evolving fuels for advanced engines.” The common index to the evaluate how fast the chemistry process is relative to the flow is the Damköhler number which is defined as

Da ≡

τadvection τreactions τ

The second Damköhler number is DaII ≡ τ diffusion . Da 1 implies that the problem is flow controlled and the chemistry reactions tends to equilibrium chemistry. Non-premixed combustion is considered to be equilibrium chemistry. Da 1 implies that the reactions control the system; this is finite rate chemistry. Premixed combustion is considered to be an example of finite rate chemistry. Equations (2.4) and (2.5) are an example of a stiff system of ordinary differential equations for chemical kinetics. The system (chemical mechanism) can be very large. The number of species defines the size of the system. However, the complexity, in term of CPU time, is also determined by the number of reactions. For most reaction mechanisms, the number of reactions is greater than the number of species – for example the O–O2 –O3 system contains three species and four reactions. The H2 O–O2 –H2 process can be defined with six species and eight reactions or together with H2 O2 there are nine species and 19 reactions. The NOx system (used for hypersonic boundary layer calculations) involves six species and five reactions. Hydrogen-Air combustion contains 13 species and 33 reactions. Methane (CH4 ) combustions include 15 species and 58 reactions. The GRI-Mech 3.0 model which is used to model natural gas combustion contains 53 species and 325 reactions. Furthermore, concerning finite rate chemistry, Yee [37] states that: “for stiff reactions containing shock waves, it is possible to obtain stable solutions that look reasonable and yet are completely wrong”. Unfortunately, many of the benchmarks for reactive flow such as the blunt-projectile example and the tiled wedge are shock-induced combustion problems. 6.2. Approaches to reactive flow solvers 6.2.1. Operator splitting Given the equation ut + A (u ) = 0 where the operator A is a sum of operators with different properties – A (u ) = j A j (u ), the idea of operator splitting is to decompose the problem into simpler problems and to solve each one separately with a scheme that ensures stability. For example:

ut = Au + Bu ⇒ u (t ) = et ( A + B ) u 0 Lie splitting: v (t ) = et A et B u 0 or equivalently v (t ) = et B et A u 0 . Strang splitting: v (t ) = et A /2 et B et A /2 u 0 or v (t ) = et B /2 et A et B /2 u 0 . Unless the operators A and B commute, the results depend on the order of the splitting and the error is a function of AB–BA. In general, u (t ) − v (t ) = O (t p ). We would like p to be at least as large as the accuracy of the basic schemes. The advantages of using splitting methods are that it is easy to implement and that the flow and chemistry are solved with solvers that are known to be stable and accurate for each “stage”. The disadvantage of the splitting methods is the intrinsic splitting errors. 6.2.2. Coupled solver For a coupled problem one can consider three methods: fully implicit, including multigrid, point or (semi) implicit (treating only the source term implicitly) or explicit with an implicit smoother (current work). All of the implicit methods, including schemes for stiff ODE systems use the source term Jacobian. For some of the methods, such as point implicit and for the implicit smoother, the Jacobian can be approximated. Most researchers take a diagonal form for the approximation of the Jacobian. For example, Bussing and Murman [38] used only the diagonal elements and approximated the Jacobian as R j j = τ1 . τ j is the characteristic time of the creation or destruction of specie j. Fig. 6.1 presents the sparsity of the Jacobian j

as a function of the number of species (from Perini et al. [39]). For all sizes of chemical mechanisms all of the entries on the diagonal are included. Other properties of the matrix are rarely discussed. According to Lian et al. [40] one can use the knowledge that some of the eigenvalues of the matrix are negative


215

Fig. 6.1. Sparsity of the Jacobian matrix vs. reaction mechanism dimension. From Perini et al. (2004).

6.3. The sources term Jacobian In primitive variables {ω, T } the full Jacobian matrix has the following structure:

⎛

∂ ω˙ j ... ⎜ ∂ ωk ⎜ ⎜ NxN ⎜ ⎜ ∂ T˙ ⎜ ⎜ ... ⎝ ∂ω j

...

∂ ω˙ j ∂T

⎞

1xN

⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ˙ . . . ∂∂ TT 1x1 ⎟ ⎠

Nx1

We now show that under reasonable physical assumptions, the diagonal entries are non-positive for most situations. This is also expected from Le Chatelier’s principle which states that when a chemical system in equilibrium is subject to a change in conditions such as temperature, pressure, concentrations etc., it will reach a new equilibrium that minimizes the changes that caused the equilibrium violation. In particular, we prove that individual reactions that behave according to the Le ∂ ω˙

˙

Chatelier’s principle contribute negative parts to ∂ ωk and to ∂∂ TT . Reactions which do not obey the principle, may contribute k a positive part. ∂ ω˙

6.4. The sign of ∂ ωk k Define the forward and backward progress variables q f and qr :

q f i ≡ K f ,i ( T )

υ ik

ωk

qri ≡ K r ,i ( T )

k

υ ik

ωk

k

so the rate-of-progress variable defined in equation (2.6) can be written as q i = αmi ωim {q f i − qri }. From equation (2.4) r ∂ ω˙ j ∂q = υi j i ∂ω j ∂ω j

N

(6.1)

i =1

With either

υij = 0 or υij = 0 but not both so if υij = 0 then υik = υik and if υij = 0 then υik = −υik . ∂q

We compute ∂ ωi for three types of reactions: j A. Reactions where species j does not participate as a product or reactant ∂q

For these reactions ∂ ωi = 0 and j

υi j = 0. These reactions do not contribute to the sum in (6.1).

216


B. Reactions where species j participates in the reaction and

ω im = ω j

∂ qi ωmi = αmi υi j q f ,i − υij qr ,i ∂ω j ωj ∂ ω˙

The contribution to ∂ ωk is k

υik αmi

So,

⎧ ωmi 2 ⎪ ⎪ ⎨ −υik αmi ω q f ,i υik = 0

ωmi j {υ q − υij qr ,i } = ωmi ⎪ 2 ω j i j f ,i ⎪ − υ α qr ,i υik = 0 ⎩ mi ik ωj

υik αmi ωωmij {υij q f ,i − υij qr ,i } < 0. Hence, reactions in group B contribute negative parts to

C. Reactions where species j participates in the reaction and

∂ ω˙ k ∂ ωk .

ω im = ω j

∂ qi = αmi υij q f ,i − υij qr ,i + {q f ,i − qr ,i } ∂ω j

(6.2) ∂ ω˙

The first term of equation (6.2) contributes a negative part to ∂ ωk similar to reactions of group B. The second term q f − qr is k positive if the reaction is forward going and negative if the reaction is backward going. The only possible positive contribu∂ ω˙

tions to ∂ ω j is the outcome of backward going reactions with a third body that is also reactant or forward going reactions j with a third body that is also a product in the reaction. From the physical point of view, the third body is a chemical catalyst so this is the situation in which the species accelerates its own production. These reactions contradict the Le Chatelier’s principle. ˙ 6.5. The sign of ∂∂ TT

For a constant internal energy process,

T˙ = −

W

ρCv

˙k = − Ukω

k

W

ρCv

Uk

Nr

υik qi = −

i =1

k

Nr W

ρCv

q i U i .

i =1

Assuming that C v = ∂∂ UT is not a function of the temperature (as for perfect gas) some useful relations are:

∂k f E kf = β+ ∂T RT T ∂q ∂ kr E U kr ,i E q U k r , i ν3 ν4 = β+ − = = β+ + ω ω ∂T RT RT T ∂T RT T RT T C D

Differentiating T˙ with respect to the temperature yields

ρ C v ∂ T˙ W R ∂T

=−

Nr U i RT

i =1

1 + βi +

Ei

RT

qi +

U i RT

qir .

The Arrhenius coefficients A, E and β are given from empirical data, but they don’t have significant physical meaning. The more physical activation energy is defined by

E a (T ) = −

d ln k d(1/RT )

(6.3)

.

Substituting k f from equation (2.7) into (6.3) yields the connection between E, E a and β

E a ( T ) = β RT + E In terms of the activation energy E a , we have


=−

Nr U i i =1

RT

1+

E ai RT

qi +

U i RT

qir .


217

From chemical kinetics theory,

k( T ) =

1

3/2 #∞

2

σ (ε; E 0 )ε exp(−ε/RT )dε

(πμ)1/2 k B T

E0

E 0 is the positive threshold energy for the reaction and σ (ε ; E 0 ) is the collision cross-section between the molecules. The cross section is given from theoretical chemistry. For example, the transitional energy along the center line is:

σ (ε) =

1 − E 0 /ε 0

ε ≥ E0 ε < E0

gives k( T ) ∝ T 1/2 e − E 0 /RT so that β = 1/2 and E = E 0 > 0. More complicated models for the collision cross section predict values that can be both positive and negative for the temperature exponent β but gives a strong restriction on the relation of β and E to have negative Ea. Negative Activation energy is the limit of threshold energy approaching zero (barrierless reactions) (Jaffe [41]). Define the backward activation energy as E bi = E ai − U i . E bi is also assumed to be positive. Reactions with negative U i are called exothermic and reactions with positive U i are called endothermic reactions. In terms of E a and E b we have


=−

Nr E ai − E bi

E ai

RT

RT

i =1

1+

qfi − 1+

E bi

qri .

RT

˙

We now analyze the contribution of each individual reaction to the sign of ∂∂ TT . The following reactions yield a non-positive contribution: 1. Reactions in equilibrium – q = q f − qr → 0 2. Energy neutral reaction – U = E ai − E bi → 0 3. Forward going (q > 0) endothermic (E a > E b ) reactions and backward going (q < 0) exothermic (E a < E b ) reactions. 6.5.1. The limit of temperature approach to zero ˙

We can rewrite T˙ and ∂∂ TT as


ρ C v T˙ W R T

=−

=−

Nr E ai − E bi

RT

i =1 Nr

E ai − E bi RT

i =1

Sˆ k = S k0 − R ln

A i T βi e βi

P atm

ωk

k

A i T βi e βi

υ ik

ωk

1+

E ai RT

e − E ai /RT − 1 +

e − E ai /RT − e − E bi /RT exp −

k

P

υ ik

1 R

E bi

RT

Sˆ +

e − E bi /RT exp −

1 R

Sˆ +

υik

k∈ K i

υik

k∈ K i

− R ln Xk . ˙

∂ ω˙

˙

As T approaches zero, both T˙ and ∂∂ TT are positive, similarly for ∂ ωk . A positive contribution to ∂∂ TT can be due to backward k going endothermic reactions or else to forward going exothermic reactions. These reactions also contradict the Le Chatelier’s principle. The Le Chatelier’s principle itself does not have an exact mathematical formulation [42] and many works discuss the validity and possible violations of the principle [43–46]. 6.5.2. Numerical examples

˙ Although the sign of ∂∂ TT is determined from empirical data, as mentioned above, and we can’t prove it for the general ∂ p˙

case, nevertheless, in most cases it is negative. The value of ∂ p is given by

∂ T˙ ∂ P˙ T ∂ P˙ RT Ea Eb = = σi 1 + q f ,i − 1 + qr + . ∂P P ∂T P RT RT ∂T i

∂ p˙

In the following two examples ∂ p is presented. The first example is a flow of stoichiometric mixture of hydrogen and oxygen in a rapid expanded diffuser [47]. Typically CFL number used for this example is about 80. The second example is a flow of stoichiometric mixture around blunt projectile. Typically CFL number is about 10. This is an example for shock induced ∂ p˙ combustion. As can be seen Figs. 6.2 and 6.3, indeed ∂ p is negative in most of the computational domain.

218


Fig. 6.2. ∂ p˙ /∂ p for rapid diffuser problem [47].

Fig. 6.3. ∂ p˙ /∂ p for the shock induces blunt projectile problem [108].

6.6. Determination of the temperature from the internal energy For a given internal energy e 0 , we want to determine the temperature (and/or the pressure). Since for non-ideal gases the internal energy is a non-linear function of the temperature, we solve it numerically using the Newton–Rapson method. We solve the equation f ( T ) = e ( T ) − e 0 = 0 iteratively:

T n +1 = T n −

f (T ) f (T )

= Tn −

f (T n ) f (T n )

.

c ( T ) T −e ( T )+e

de ( T )

Since, dT = c v we get T n+1 = v n cn ( T ) n 0 . v n When the algorithm converges, it gives also the temperature and the heat capacity ratio

γ=

cp cv

=

c v + R /W cv

γ


219

7. Turbulence simulations 7.1. Considerations for turbulence modeling Turbulent flow is a flow that is characterized by a random velocity perturbation in all spatial and time scales (eddies). Those perturbations are measured in terms of the turbulent kinetic energy (TKE) which is defined as:

k=

3 2 u + v 2 + w 2 2

u , v and w are the components of the perturbation vector to the velocity. The transition from laminar to turbulence depends strongly on the Reynolds number. For large enough Re, the flow is considered to be turbulent. The randomness of the turbulent flow implies that it is a time dependent phenomena. However, turbulence can be modeled as a steady state flow or as a time dependent flow. The randomness of the turbulent flow implies that it is a time dependent phenomena. However, turbulence can be modeled as a steady state flow or as a time dependent flow. There are three basic alternatives to treat turbulent flow: a. Reynolds Averaged Navier–Stokes (RANS) – Modeling all turbulence scales – spatial and temporal. In this case there is an additive viscosity called the turbulent viscosity. The turbulent viscosity is calculated via an algebraic turbulent model such as Baldwin–Lomax model [48], or by a differential model, with additional turbulent flow variables. Time dependent RANS is called URANS (Unsteady RANS). b. Large Eddies Simulations (LES) – Modeling small scale turbulence (sub-grid) and resolving large scales. For LES, a model to the sub-grid eddy viscosity needs to be supplied. LES is a time dependent solution. c. Direct Numerical Simulation (DNS) – In the DNS, all turbulent scales are resolved. There is no turbulence model. LES requires much more computer resources than URANS. A typical time dependent LES simulation can take weeks. DNS take even more of computational efforts. In this work we are focus on a RANS model and work with a one equation model – Spalart–Allmaras [13] model and with a 2-equation model – SST k–ω model [49]. 7.2. Grid requirements for RANS modeling RANS simulations require that the non-dimensional wall distance y + will be of the order or less than one where y + is defined as:

y+ =

u∗ y ∗ u

ν

is the friction velocity.

For transonic flow, this leads to a mesh y normal to the wall of about 10−6 . If the length of an airfoil is about 1 and we use about 100 points along the airfoil, the grid aspect ratio is in the order of 1000. ρ ux ρ ux The thickness of a turbulence boundary layer can be estimated from δ = 0.382 x0.2 Rex = μ . δ = 0.382 x0.2 , Rex = μ . Rex

Rex

If we let the grid spacing increases geometrically, then yn = y 1 qn−1 where q is the grid growth rate. We then obtain: δ y1

=

q N −1 . q−1

We estimate the number of points required in the boundary layer by

N = ln

δ y1

0.382 L (q − 1) + 1 / ln q = ln (q − 1) + 1 / ln q 0 .2 Re L

y1

For q = 1.1, L = 1, Re = 1.5 × 106 and y 1 = 10−6 this gives about 80 points inside the boundary layer. In practical grids for computations for airfoils and wings there are about 20 points inside the boundary layer. Large aspect ratios (AR) are a numerical source of stiffness because of the different scales in the normal and tangential directions. Multi-block structuredgrid based solvers are very sensitive to a high aspect ratio since it is not easy to avoid the prolongation of the AR to irrelevant parts of the grid [50]. For example, Multi-block, body-fitted, structured grid for NLR-7301 airfoil with flap, can be seen in Fig. 7.1. High aspect ratio grid cells from the wake of the main airfoil are prolonged above the flap. 7.3. k–ω SST model 7.3.1. Boundary conditions and multigrid connection The dependency of the wall boundary condition for the k–ω SST turbulent model on the first grid cell distance y 1 , is a difficulty for multigrid where various size grids are used. There are several strategies for computing turbulent flow with multigrid: (1) Use multigrid only for the mean flow variables. Then one computes k and ω only on the finest grid and restricts the eddy viscosity or k and ω and compute the turbulent viscosity on the coarse grids. (2) Use multigrid only for the mean flow variables and compute k and ω on each

220


Fig. 7.1. Multi-block structured grid for NLR 7301 airfoil with flap. In top-left window is a zoom on the flap.

grid level without coarse grid corrections. (3) Use multigrid method also for the turbulent variables. Options 2 and 3 require the specification of ω at the boundary also on the coarse grids. However, the dependency of the specified ω on the grid size is not well matched with multigrid theory. Moreover, the coarse grids are often not suitable for turbulent simulations in the sense of y + 1) is needed. In Fig. 8.15, a zoom on the motor’s head is presented for a sequence of three particle diameters (1, 5 and 20 microns). The colors represent the disperse phase density and the black lines are streamlines. Close to the injection corner there is a zone with a relatively low density of disperse a phase. We call this area a void zone. The void zone becomes “deeper” as the particle diameter increases. In Fig. 8.16, we present a zoom on the motor aft section. There are now three void zones

226


Fig. 8.10. Typical temperature (top figure) and Mach number (bottom) contour maps.

Fig. 8.11. Convergence history for BEM simulation.

Table 8.1 Inflow condition for the reactive shear layer problem.

T (K) U (m/sec) P (pa) O2 mass fraction N2 mass fraction H2 mass fraction

Oxidizer

Fuel

Torch

811 390 1e+5 0.21 0.79 0

367 140 1e+5 0 0.9606 0.0394

1250 140 1e+5 0.21 0.777 0.013


227

Fig. 8.12. Typical dispersed phase density (5 microns particles).

Fig. 8.13. Convergence history – gas phase density residual.

– close to the injection corner, near the wall before the throat and after the throat close to the wall of the nozzle. The void zone becomes deeper as the particle diameter increases. As the particles get bigger, they tend to concentrate along the motor axis and not to fill the whole nozzle. This is due to their increased inertia. In Fig. 8.17, a zoom on the motor head is presented for the one micron particle diameter. In this figure, the red lines are the Lagrangian trajectories. The blue lines are the Eulerian streamlines. This way we can compare the results of two very different calculation methods in order to verify the results. In Fig. 8.18, the Lagrangian trajectories vs. the Eulerian streamlines in the motor aft section are presented. A closer view near the wall is also shown. It can be seen that the Lagrangian trajectories intersect. This is a physical phenomenon that cannot be simulated in the EDP model. 8.5. Reactive shear layer This example demonstrates the improvement using multigrid for a mid-subsonic reactive flow case. Detailed discussion about reactive flows appears in section 4. The literature is quite limited in such examples. Most of the test cases are either for premixed combustions where the ignition is due to shock waves (shock induced ignition) or else for very low Mach number (as flames). We consider an experiment with data from Chang et al. [62]. The results, are compared to a numerical simulation of Edwards [63] and qualitatively to the UV image obtained by Chang et al. The problem consists of a two layer fluid. The top one is a slow speed and cold flow of fuel (mostly molecular hydrogen). The bottom layer is a fast and hot oxidizer (air) flow. In between, there exists a third, thin layer of a mixed and very hot

228


Fig. 8.14. Convergence history – dispersed phase density residual for various particle diameters. The convergence is getting better for small particles.

Fig. 8.15. Dispersed phase density distribution in a zoom on the motor head for three particle diameters. The black lines are streamlines.


229

Fig. 8.16. Dispersed phase density distribution in a zoom on the motor aft section for three particle diameters. The black lines are streamlines.

Fig. 8.17. Lagrangian trajectories vs. Eulerian streamlines in the motor head. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

layer. Numerically, this middle layer represents the stream that was produced by a torch in the experiment. Chang et al. reported that no simultaneous ignition occurred without the torch. The boundary conditions are described in Fig. 8.19. The physical data that is given, for this case, is for the primitive set of p, T and u. For subsonic flow, we have to define two variables and extrapolate the third one from the exterior domain. A possible way to define these three variables ( P , T , u) γ p is to specify the entropy S ( p , T ) together with the total enthalpy – h = γ −1 ρ + 12 u 2 , which is a conserved variable and extrapolate the pressure from the interior. The detailed boundary conditions are presented in Table 8.1. For this case we use a 128 × 128 grid with an equal distribution in the x direction. In the y direction the mesh is clustered toward the centerline.

230


Fig. 8.18. Lagrangian trajectories vs. Eulerian streamlines in the motor aft section.

Fig. 8.19. Boundary condition description for the reactive shear layer problem.

Fig. 8.20. Temperature contour map for the reactive shear layer. The shiny spot at x = 0 in the middle is the numerical “torch”.

The spatial scheme is AUSMD/V with a v-cycle multigrid. The problem is turbulent and the Spalart–Allmaras one-equation model is applied. The CFL number is 64. The resulting temperature distribution is presented in Fig. 8.20. Comparison to Edwards of the temperature at certain cut-line transverse to the flow direction (constant x) is shown in Fig. 8.21. In Fig. 8.22, the OH mass fraction is presented and we also present a UV OH emission photograph (from Chang et al.). A quantitative comparison to Edwards for the H2 and H2 O mass fraction at x = 0.44 is presented in Fig. 8.23. The convergence history for a single grid run, two and three grids multigrid is shown in Fig. 8.24. There are some differences between the results of this work and the results from Edward, especially in the temperature. In many reactive flow simulations, the result for the chemical composition is more


231

Fig. 8.21. Cut-line temperature distribution at constant x = 0.44.

Fig. 8.22. A) Simulated OH mass fraction distribution. (B) Experimental OH UV emission from Chang et al. In the simulation results, the shiny spot is the numerical “torch” which is not the same as in the experiment.

accurate than for the temperature. One possible reason for the difference is that Edwards take the width of the torch to be one millimeter and in this work the width of the torch is 4 millimeters. We didn’t investigate the sensitivity to the torch’s width. 8.6. UV radiation from rocket motor plume RANS simulation results are a time averaged procedure. We present the computational results of a BEM motor plume. In Fig. 8.25, the top right image is a single frame UV image. In the top-left, a 0.4 second averaged frame is presented. The bottom figures present the simulated OH mole fraction and the temperature distribution in the plume. We note the similarity to the averaged frame.

232


Fig. 8.23. Cut-line H2 O and H2 distribution at constant x = 0.44.

Fig. 8.24. Convergence history for reactive shear layer problem.

9. Conclusions We solved complex physical flows including flow speeds that vary from almost zero to hypersonic flow. Reactive flow, both premixed and non-premixed was successfully solved including the acceleration of the multigrid method for pure subsonic reactive flow. We did not use multigrid to accelerate the convergence of supersonic flow. The computation was successfully done for the turbulent variables using the RK implicit preconditioning (smoother) for the k and ω variables of the k–ω SST model and for ν˜ of the Spalart–Allmaras model. The main difficulty in applying the multigrid method for the k–ω model is the boundary conditions which depend on the grid size. For low Mach numbers flow, the Roe scheme using c coupled with the RK/Implicit smoother operator that we developed converged rapidly and accurately. The reactive flow problem is the most complex problem that was solved in the sense of stiffness and numerical issues of numerical instability. The solution of reactive flow using the RK/Implicit smoother requires the evaluation of the source term Jacobian. It was found that the construction of the Jacobian matrix consumes so much of the computational time that it eliminates the algorithmic gain in terms of CPU time that is achieved using the smoother. Hence, the Jacobian must be approximated. The approximated Jacobian can affect the convergence but it also can affect the solution so it must be done carefully. We studied analytically its properties based on physical/chemical considerations. The main obstacle in this


233

Fig. 8.25. UV images and simulation results of RANS calculation if BEM plume.

physical study of the source term Jacobian is the empirical nature of the Arrhenius kinetic model which produces entries that are not necessarily physical. The dispersed phase flow gives another kind of coupling between the mean flow system and added model variables (symmetric coupling via source term). One advantage of the RK/Implicit smoother algorithm for multi-physics problems is in its modularity. In order to add new physical considerations to the problem, which can be represented either by new equations to be solved for the added variables or new source terms one should formulate the new equations and the proper Jacobian and add it to the code. The implicit smoother enabled us to solve the complex physical problems with relatively low CPU time cost and a numerical robustness. References [1] A. Jameson, W. Schmidt, E. Turkel, Numerical solutions of the Euler equations by finite volume methods using Runge–Kutta time-stepping schemes, in: AIAA Paper, 1981, p. 1259. [2] R. Swanson, E. Turkel, On central-difference and upwind schemes, J. Comput. Phys. 101 (2) (1992) 292–306. [3] A. Jameson, T. Baker, Multigrid solution of the Euler equations for aircraft configurations, in: AIAA Paper, 1984, p. 93. [4] D. Caughey, A. Jameson, How many steps are required to solve the Euler equations of steady, compressible flow-in search of a fast solution algorithm, in: 15th AIAA Computational Fluid Dynamics Conference, 2001, p. 2673. [5] E. Turkel, Preconditioned methods for solving the incompressible and low speed compressible equations, J. Comput. Phys. 72 (2) (1987) 277–298. [6] E. Turkel, Review of preconditioning methods for fluid dynamics, Appl. Numer. Math. 12 (1993) 257–284. [7] C.C. Rossow, Efficient computation of compressible and incompressible flows, J. Comput. Phys. 220 (2) (2007) 879–899. [8] R. Swanson, E. Turkel, C.C. Rossow, Convergence acceleration of Runge–Kutta schemes for solving the Navier–Stokes equations, J. Comput. Phys. 224 (1) (2007) 365–388. [9] S. Langer, Accuracy investigations of a compressible second order finite volume code towards the compressible limit, Comput. Fluids 149 (2017) 119–137. [10] J. Slotnick, A. Khodadoust, J. Alonso, D. Darmofal, W. Gropp, E. Lurie, D. Mavriplis, CFD Vision 2030 Study: A Path to Revolutionary Computational Aerosciences, 2014. [11] P.R. Spalart, V. Venkatakrishnan, On the role and challenges of CFD in the aerospace industry, Aeronaut. J. 120 (1223) (2016) 209–232. [12] F. Menter, C. Rumsey, Assessment of two-equation turbulence models for transonic flows, in: Fluid Dynamics Conference, 1994. [13] P. Spalart, S. Allmaras, A one equation turbulence model for aerodynamic flows, AIAA J. 94 (1992). [14] S.R. Allmaras, F. Johnson, Modifications and clarifications for the implementation of the Spalart–Allmaras turbulence model, in: ICCFD7-Seventh International Conference on Computational Fluid Dynamics, Big Island, Hawaii, 2012. [15] M. Golafshani, H. Loh, Computation of two-phase viscous flow in solid rocket motors using a flux-split Eulerian–Lagrangian technique, in: AIAA Paper, 1989, p. 2785. [16] S. Dash, D.E. Wolf, R.A. Beddini, H.S. Pergament, Analysis of two-phase flow processes in rocket exhaust plumes, J. Spacecr. Rockets 22 (3) (1985) 367–380. [17] M. Liou, C.J. Steffen Jr, A new flux splitting scheme, J. Comput. Phys. 107 (1) (1993) 23–39. [18] M. Liou, A sequel to AUSM: AUSM+, J. Comput. Phys. 129 (2) (1996) 364–382. [19] M. Liou, A sequel to AUSM, Part II: AUSM+-up for all speeds, J. Comput. Phys. 214 (1) (2006) 137–170. [20] A. Brandt, Multi-level adaptive solutions to boundary-value problems, Math. Comput. 31 (138) (1977) 333–390. [21] U. Trottenberg, C.W. Oosterlee, A. Schuller, Multigrid, Academic Press, 2000. [22] M. Wasserman, Y. Mor-Yossef, J. Greenberg, A robust implicit multigrid method for RANS equations with two-equation turbulence models, J. Comput. Phys. 229 (6) (2010) 5820–5842. [23] O. Peles, S. Yaniv, E. Turkel, Convergence acceleration of Runge–Kutta schemes using RK/implicit smoother for Navier–Stokes equations with SST turbulence, in: Proceedings of 52nd Israel Annual Conference on Aerospace Sciences, 2012.

234


[24] O. Peles, E. Turkel, S. Yaniv, Fast iterative methods for Navier–Stokes equations with SST turbulence model and chemistry, in: Seventh International Conference on Computational Fluid Dynamics (ICCFD7), 2012. [25] R. Swanson, C. Rossow, An Initial Investigation of the Effects of Turbulence Models on the Convergence of the RK/Implicit Scheme, NASA/TM-2008-215342, 2008. [26] S. Langer, Agglomeration multigrid methods with implicit Runge–Kutta smoothers applied to aerodynamic simulations on unstructured grids, Comput. Phys. 277 (2014) 72–100. [27] E. Turkel, A. Fiterman, B. Van Leer, Preconditioning and the Limit to the Incompressible Flow Equations (No. ICASE-93-42), Institute for Computer Applications in Science and Engineering Hampton VA, 1993. [28] E. Turkel, V. Vatsa, Local preconditioners for steady state and dual time stepping, Math. Model. Numer. Anal., ESAIM: M2AN 39 (3) (2005) 515–536. [29] E. Turkel, R. Radespiel, N. Kroll, Assessment of preconditioning methods for multidimensional aerodynamics, Comput. Fluids 26 (2) (1997) 613–634. [30] B. Van Leer, W. Lee, P. Roe, Characteristic time-stepping or local preconditioning of the Euler equations, in: 10th Computational Fluid Dynamics Conference, 1991. [31] Y. Choi, C. Merkle, The application of preconditioning in viscous flows, J. Comput. Phys. 105 (2) (1993) 207–223. [32] D. Darmofal, B. Van Leer, Local preconditioning: manipulating mother nature to fool father time, in: Computing the Future II: Advances and Prospects in Computational Aerodynamics, 1998. [33] A. Jameson, Time dependent calculations using multigrid, with applications to unsteady flows past airfoils and wings, in: AIAA Paper, 1991, p. 1596. [34] E. Turkel, V. Vatsa, Local preconditioners for steady and unsteady flow applications, ESAIM: Math. Model. Numer. Anal. 39 (3) (2005) 515–535. [35] G. Sod, A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws, J. Comput. Phys. 27 (1) (1978) 1–31. [36] C. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes, II, J. Comput. Phys. 83 (1989) 32–78. [37] H. Yee, Building blocks for reliable complex nonlinear numerical simulations, Turbulent Flow Comput. (2004) 199–236. [38] T. Bussing, E. Murman, Finite-volume method for the calculation of compressible chemically reacting flows, AIAA J. 26 (9) (1988) 1070–1078. [39] F. Perini, E. Galligani, R.D. Reitz, A study of direct and Krylov iterative sparse solver techniques to approach linear scaling of the integration of chemical kinetics with detailed combustion mechanisms, Combust. Flame 161 (5) (2014) 1180–1195, https://doi.org/10.1016/j.combustflame.2013.11.017. [40] C. Lian, G. Xia, C. Merkle, Impact of source terms on reliability of CFD algorithms, Comput. Fluids 39 (10) (2010) 1909–1922. [41] R. Jaffe, On the Possibility of Negative Activation Energies in Bimolecular Reactions, NASA-TM-78509, 1978. [42] C. Olivera-Fuentes, C. Colina, Stability, Displacement and Moderation of Chemical Equilibrium: Rediscovering Le Chatelier’s Principle, 2007. [43] D.S. Corti, E.I. Franses, Exceptions to the Le Chatelier principle, Chem. Eng. Educat. 37 (4) (2003) 290–299. [44] Y. Liu, Y. Liu, M.G. Drew, A mathematical approach to chemical equilibrium theory for gaseous systems IV: a mathematical clarification of Le Chatelier’s principle, J. Math. Chem. 53 (8) (2015) 1835–1870. [45] Z.K. Liu, J. Ågren, M. Hillert, Application of the Le Chatelier principle on gas reactions, Fluid Phase Equilib. 121 (1) (1996) 167–177. [46] M. Hillert, Le Chatelier’s principle—restated and illustrated with phase diagrams, J. Phase Equilib. 16 (5) (1995) 403–410. [47] J.P. Drummond, M.Y. Hussaini, T.A. Zang, Spectral methods for modeling supersonic chemically reacting flowfields, AIAA J. 24 (9) (1986) 1461–1467. [48] H. Lomax, B. Baldwin, Thin layer approximation and algebraic model for separated turbulent flows, AIAA J. (1978) 78–257. [49] F. Menter, Improved Two-Equation k-Turbulence Models for Aerodynamic Flows, NASA Technical Memorandum, 103975, 1992. [50] J. Vierendeels, K. Riemslagh, E. Dick, A multigrid semi-implicit line-method for viscous incompressible and low-Mach-number flows on high aspect ratio grids, J. Comput. Phys. 154 (2) (1999) 310–341. [51] F. Menter, Improved Two-Equation k-Omega Turbulence Models for Aerodynamic Flows, NASA STI/Recon Technical Report N, 93, 22809, 1992. [52] J. Kok, S. Spekreijse, Efficient and Accurate Implementation of the k–ω Turbulence Model in the NLR Multi-Block Navier–Stokes System, 2000. [53] K. Yonezawa, Y. Yamashita, Y. Tsujimoto, Y. Watanabe, K. Yokota, Effect of nozzle contour on flow separation in overexpanded rocket nozzles, J. Fluid Sci. Technol. 2 (2007) 97–108. [54] H. Lüdeke, J.B. Calvo, A. Filimon, Fluid structure interaction at the Ariane-5 nozzle section by advanced turbulence models, in: H. Lüdeke, J.B. Calvo, A. Filimon (Eds.), European Conference on Computational Fluid Dynamics, Egmond aan Zee, the Netherlands, 2006. [55] K. Squires, J. Forsythe, S. Morton, D.S.M. Blake, K. Wurtzler, B. Airplanes, Analysis of full aircraft with massive separation using detached-eddy simulation, in: Proceedings of the High Performance Computing Modern, 2002. [56] E. Marineau, J. Schetz, R. Neel, Turbulent Navier–Stokes simulations of heat transfer with complex wall temperature variations, J. Thermophys. Heat Transf. 21 (3) (2007) 525–535. [57] O. Peles, Y. Sudman, S. Yaniv, Conjugate heat transfer and turbulent Navier–Stokes simulations in a convergent divergent nozzle, in: Proceedings of the 54nd Israel Annual Conference on Aerospace Sciences, 2014. [58] S. Deck, P. Duveau, P. d’Espiney, P. Guillen, Development and application of Spalart–Allmaras one equation turbulence model to three-dimensional supersonic complex configurations, Aerosp. Sci. Technol. 6 (3) (2002) 171–183. [59] S. Catris, B. Aupoix, Density corrections for turbulence models, Aerosp. Sci. Technol. 4 (1) (2000) 1–11. [60] V. Schmitt, C. F., Pressure distributions on the ONERA-M6-Wing at transonic Mach numbers, in: Experimental Data Base for Computer Program Assessment, vol. 4, 1979. [61] N. Gregory, C. O’reilly, Low-Speed Aerodynamic Characteristics of NACA 0012 Aerofoil Section, Including the Effects of Upper-Surface Roughness Simulating Hoar Frost, HM Stationery Office, 1973. [62] C. Chang, C. Marek, C. Wey, R. Jones, M. Smith, Turbulence Measurement in a Reacting and Non-Reacting Shear Layer at a High Subsonic Mach Number, 1993. [63] J. Edwards, A low-diffusion flux-splitting scheme for Navier–Stokes calculations, Comput. Fluids 26 (2) (1997) 635–659.