Large eddy simulation of hydrogen combustion in supersonic flows using an Eulerian stochastic fields method

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y x x x ( 2 0 1 6 ) 1 e1 2

Available online at www.sciencedirect.com

ScienceDirect journal homepage: www.elsevier.com/locate/he

Large eddy simulation of hydrogen combustion in supersonic flows using an Eulerian stochastic fields method Cheng Gong a,*, Mehdi Jangi b, Xue-Song Bai a, Jian-Han Liang c, Ming-Bo Sun c a

Division of Fluid Mechanics, Dep. of Energy Sciences, Lund University, Box 118, SE 221 00 Lund, Sweden School of Engineering and Information Technology, Murdoch University, WA 6150, Australia c Science and Technology on Scramjet Laboratory, National University of Defense Technology, Changsha 410073, China b

article info

abstract

Article history:

An Eulerian Monte-Carlo approach, the so-called Eulerian stochastic fields (ESF) method is

Received 5 July 2016

implemented and evaluated for simulation of non-premixed hydrogen/air combustion in

Received in revised form

supersonic flows. The ESF method is integrated into a compressible flow large eddy

22 August 2016

simulation (LES) solver, and validated on a supersonic combustor with a strut as flame-

Accepted 4 September 2016

holder. Comparison with experimental data and with results from a well-stirred reactor

Available online xxx

(WSR) model demonstrates the advantage of the LES-ESF method for simulation of localextinction and re-ignition phenomena. The hydrogen/air flame structure and the stabili-

Keywords:

zation of the combustion process in the supersonic combustor are analysed based on the

Transported probability density

present LES-ESF method. Oscillation of the recirculation zones is found to be the dominant

function

mechanism for the local-extinction/re-ignition and the flame stabilization under the pre-

Eulerian stochastic fields

sent condition.

Large eddy simulation

© 2016 Hydrogen Energy Publications LLC. Published by Elsevier Ltd. All rights reserved.

Supersonic Local-extinction Re-ignition

Introduction Hydrogen combustion in supersonic flows is seen in highspeed propulsion systems, for example, scramjet engines [1,2]. The process involves mixing of the fuel and the oxidizer, ignition and stabilization of the flames, and shock/flame interaction. Due to the short residence time and high scalar dissipation rate/strain rate in the flow field the flame can be

locally quenched (which may be followed by re-ignition) especially in the leading front of the flame. It is therefore a challenge to stabilize the flame. Simulation of supersonic combustion is important for the development of high-speed propulsion engines, which requires the development of robust and predictive models for the flow transport and flow/ flame interaction. Predictive modelling of turbulent combustion provides a powerful tool in the design of many combustion facilities, e.g., internal combustion engines and gas

* Corresponding author. Fax: þ46 462224717. E-mail address: [email protected] (C. Gong). http://dx.doi.org/10.1016/j.ijhydene.2016.09.017 0360-3199/© 2016 Hydrogen Energy Publications LLC. Published by Elsevier Ltd. All rights reserved. Please cite this article in press as: Gong C, et al., Large eddy simulation of hydrogen combustion in supersonic flows using an Eulerian stochastic fields method, International Journal of Hydrogen Energy (2016), http://dx.doi.org/10.1016/j.ijhydene.2016.09.017

2

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y x x x ( 2 0 1 6 ) 1 e1 2

turbine engines. A time-consuming but “model-free” method for the combustion modelling is the direct numerical simulation (DNS) approach, in which all flow and reaction scales in turbulent combustion are resolved directly. Due to the limitation of computation resources, however, it is impossible to carry out DNS on practical combustion engines. In practice, Reynolds-Averaged NaviereStokes (RANS) or large eddy simulation (LES) are employed in engineering simulations to reduce the computational cost. LES has been a promising approach in the modelling of supersonic combustion due to its compromise between computational cost and accuracy [3,4]. When the RANS or LES approach is employed, a turbulent combustion model is needed to take into account the turbulence-combustion interaction and the turbulent transport process. Probability density function (PDF) method [5] has been shown as one of the predictive and robust approaches for accommodating the effect of turbulence on the reaction rates in RANS and LES-based modelling of turbulent combustion. Recent progress in transported PDF method has been reviewed by Haworth [6]. In the PDF method, transport equations for the one-point, one-time Eulerian joint PDF of a set of variables that describe the reactive flow field are employed. The PDF transport equations are high-dimensional, upto Ns þ 5 dimensions, where Ns is the number of species in the chemistry mechanism considered in the simulations. Ns is typically much larger than the number of flow variables. Consequently, the conventional finite-volume or finitedifference schemes cannot be used to solve the PDF equations directly. Multiple strategies have been developed to solve the transported PDF equations. The most common one is the Lagrangian particle Monte-Carlo method [5]. In this method, the joint PDF is represented by a large number of notional particles that evolve according to the prescribed stochastic equations. Since the reactive variables are solved in Lagrangian frame work, the method is mesh-independent, and does not have errors arising from spatial discretization. However, the implementation of the Lagragian particle method into the conventional Eulerian CFD code is not straightforward. Special attentions are required to address issues of numerical consistency, statistical error control, and stability [6], in particular in the simulation of compressible flows. Motivated by this, the Eulerian fields based approaches have been developed for the transported PDF method. An alternative Eulerian method is the multi-environment PDF method that proposed by Fox [7]. In the multi-environment PDF method, the joint PDF is solved by a set of deterministic Eulerian transported equations, in contrast to the stochastic particle method. Recently, Koo et al. implemented the multienvironment PDF method to couple with the compressible LES simulation [8], and the method was used in simulation of supersonic combustions [9]. Another Eulerian fields method is ~o the Eulerian stochastic fields (ESF) method proposed by Valin [10] and Sabel'nikov et al. [11]. In the ESF approach, the joint PDF is represented by a number of stochastic equations, which are solved on the same mesh as the flow solver. The ESF method has been successfully used in a wide range of combustion problems, including the gas jet flames [12e16] and spray combustions [17e19]. The stochastic fields are

continuous and differentiable in space and continuous though not differentiable in time. These features of ESF method make it easy to be implemented in conventional CFD codes, with good numerical stability when coupling with the flow solver. Jaishree and Haworth [20,21] have compared the three approaches on the simulations of non-premixed flames. The main advantage of multi-environment PDF method is the high computational efficiency, while its model accuracy is lower than the other two stochastic methods. The ESF method was found to obtain the similar agreement with experiment as the Lagrangian particle method. The ESF method is usually coupled with the solver for low-Mach number flow, in which a simplified version of the enthalpy equation is applied, and solved in stochastic fields. In the simplified enthalpy equation, the exchange between kinetic energy of the flow (macro-scale motion of the molecules) and the internal energy of the mixture is neglected. This simplification is not suitable for compressible flows (with Ma > 0.3). In this work, an ESF method is presented for LES simulations in fully coupled compressible formulation. In the present work, the LES-ESF method is applied to modelling of the combustion process in the German Aerospace Center (DLR) model scramjet combustor [22]. In this supersonic combustor, hydrogen is injected into the supersonic flow from a strut. The flame is stabilized downstream the strut. Local-extinction/re-ignition was observed near the flame base [23], which results in a moderate temperature (~1000 K) at the shear-layer near the strut base. This combustor has been simulated using different turbulent combustion models, e.g., the partially-premixed reactor (PaSR) model [24], the flamelet-type models [25e30], and the well-stirred reactor (WSR) model [23]. Most of these works focused on the model validation based on the prediction of the mean reactive fields, e.g., the mean profile of temperature downstream the strut. Though reasonable agreement with the experiment could be obtained for the mean temperature in the downstream reaction region, the mean temperature near the strut base was usually over-predicted [25,27] or underpredicted [23,28,29]. The related analysis on the localextinction/re-ignition process and its role in the flame structure is relative rare. Motivated by the above observations, the present study is aimed at the following two goals. The first is to evaluate the capability of the ESF method in the modelling of turbulent combustion in supersonic flows. The second goal is to investigate the local-extinction/re-ignition process and its role in the flame structure and stabilization in the supersonic hydrogen combustion process.

Numerical methodology LES equations for compressible flow The formulation for compressible flow LES is based on the Favre-filtered conservation equations of mass, momentum and energy. These equations are obtained by separating the small-scale dynamics from the resolved scale motion using a box-filter. The equations can be expressed in the following form

Please cite this article in press as: Gong C, et al., Large eddy simulation of hydrogen combustion in supersonic flows using an Eulerian stochastic fields method, International Journal of Hydrogen Energy (2016), http://dx.doi.org/10.1016/j.ijhydene.2016.09.017

3

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y x x x ( 2 0 1 6 ) 1 e1 2

~j vr vru þ ¼0 vt vxj

(1)

 
 sgs v ti;j  tij ~j þ pdij ~i u ~i v ru vru þ ¼ vxj vxj vt ~iu ~ i vrY ~j vrY v þ ¼ vxj vt vxj

~ vY rD vxj

!

v  vxj

(2)

g ~ ~j rY i uj  rY i u

!

e_ i þ ru


  ~ v rE~u ~j þ pu ~j vrE v vksgs sgs sgs ~j ti;j þ rn þ ¼ þ qj  Hj þ sij u vt vxj vxj vxj

(3)

(4)

where overbars and tildes denote resolved-scale and Favre~ i and u e_ i are filtered resolved-scale variables, respectively. Y the resolved mass fraction and chemical reaction rate for the ~i þ rksgs is the total energy; e~ is the ~i u e þ 12 ru ith species. rE~ ¼ r~ resolved internal energy; ksgs is the sub-grid scale (SGS) turbulent kinetic energy. The Lewis numbers are assumed unity here, which results in the heat flux in Eq. (4) as qj ¼ r

n vh~ Pr vxj

~j ~i u tij ¼ rui uj  ru     sgs ~u ~j þ puj  pu ~j Hj ¼ rEuj  rE sgs

(6)

~j ti;j ¼ uj ti;j  u

The SGS stresses are closed using a one-equation eddy model: vrk vt

sgs

þ



sgs ~

vrk uj v rnt vk ¼ vxj Prt vxj vxj

sgs



~j sgs vu

 tij

vxj

  1 2 sgs tij ¼ 2rnt S~ij  S~kk dij þ rksgs dij 3 3

r 3=2  cd ðksgs Þ D

(7)

(8)

pffiffiffiffiffiffiffiffi (9) nt ¼ cm D ksgs ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p sgs 3 where D ¼ DxDyDz is the filter scale. The SGS energy flux Hj sgs and SGS viscosity dissipation work term sij are modelled as sgs

Hj

sgs

þ sij ¼ r

(12)

This model is often referred to as the well-stirred reactor (WSR) model. This model is sometimes used in numerical simulation of volumetric ignition process, while it is shown not accurate in prediction of many turbulent flames.

Eulerian stochastic fields method The composition PDF method [5] is employed to model turbulent combustion. For the one-point, one-time joint PDF P~sgs ðj; x; tÞ of the composition vector j, the modelled transport equation is given by Ref. [5].

r

  vP~sgs j

  vP~sgs j

Ns X v h   ~  i ru_ j Psgs j vja vt vxj a¼1  3 2   ~ Ns  i  v 4 m msgs vPsgs j 5 Cf X v h
~ a P~sgs j þ  r ja  f ¼ vxi s ssgs tsgs a¼1 vja vxi

~j þ ru

þ

(13) (5)

where h~ ¼ e~ þ p=r is the resolved enthalpy of the mixture. The unclosed SGS terms in Eqs. (1)e(3), including the SGS stresses sgs sgs term tij , the SGS energy fluxes Hj , and the SGS viscous sgs dissipation work sij , are given as

sgs sij

  ~ e_ i ¼ ru_ i j ru

nt Prt

~ vksgs vh~ vu ~i i þ þu vxj vxj vxj

! (10)

The SGS flux in the species mass fraction equation could be modelled in the similar way   ~ g ~ ~j ¼ r nt vY i  rY i uj  rY i u Sct vxj

(11)

The main challenge of the turbulent combustion modelling is how to closure the filtered chemical reaction rate in Eq. (3). In the simplest treatment, the SGS fluctuation effects on the filtered reaction rate may be neglected. The filtered reaction rate is calculated based on the resolved reactive scalars, i.e.,

Here, the sub-grid eddy-viscosity model and the interaction by exchanging with the mean (IEM) model are adapted to model the sub-grid transport and the micro-mixing terms, respectively. Equation (13) is solved using Eulerian stochastic fields (ESF) method proposed by Valino [10]. In the ESF method, the reactive fields are represented by N stochastic fields for each of the Ns scalars, namely xna ðx; tÞ for 1  n  N; 1  a  Ns . In this way, the joint PDF, P~sgs ðj; x; tÞ, is represented by an ensemble of the stochastic fields, N Y Ns   1 X 
P~sgs j; x; t ¼ d ja  xna N n¼1 a¼1

(14)

where d is the Dirac delta function. The mean and the moments of each variable for the LES field can be approximated from the ensemble of N stochastic fields. For example, the mean is, N X ~a ¼ 1 xn f N n¼1 a

(15)

The governing equation for the nth stochastic field is:     vxna v vxn 1 r
n dt þ Gt a dt þ ru_ a jn dt  Cf x vxi 2 tsgs a vxj vxi sffiffiffiffiffiffiffiffi n  ~ a dt þ r 2 Gt vxa dWn ; f i r vxi

~j rdxna ¼ ru

(16)

n where jn ¼ ½Y1n ; Y2n ; /; YNsp ; hn  is the composition vector for the nth field. Gt represents the total diffusion coefficient. dWn represents a vector Wiener process that is spatially uniform but different for each field. The term with Cf describes the micro-mixing process using the IEM model. The micro-mixing terms decrease the SGS variance of the reactive scalars. In the limit of infinitely fast micro-mixing the SGS PDF relaxes towards a delta Dirac, and the ESF method is equivalent to the WSR model. Cf is the model constant for the degree of micromixing with the standard value of 2.0. However, varying

Please cite this article in press as: Gong C, et al., Large eddy simulation of hydrogen combustion in supersonic flows using an Eulerian stochastic fields method, International Journal of Hydrogen Energy (2016), http://dx.doi.org/10.1016/j.ijhydene.2016.09.017

4

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y x x x ( 2 0 1 6 ) 1 e1 2

values for Cf could be found in literatures to model the combustion process under the different conditions [15,31]. The sub-grid mixing time scale tsgs is calculated as [12,14]. 2

tsgs ¼

rD m þ msgs

(17)

Numerical implementation The open source CFD library, OpenFOAM [32], is adopted for the numerical integration of the governing equations presented in Section Eulerian stochastic fields method. The compressible flow NaviereStokes equations (Eqs. (1)e(3)) are solved using a density-based finite volume scheme [33]. The convective terms are discretized using the Kurganov and Tadmor scheme [34] with the Van Leer limiter [35]. The second-order Gauss filtered linear scheme and Euler backward scheme are used for the diffusion term and temporal integration, respectively. The stochastic field equation can be written in the following general form: r

 


 vxna ~ a þ ru_ a jn ¼ [ad xna þ [wien xna þ [mix xna ; f vt

Fig. 1 e Schematic of the combustor geometry (a) and mesh combined with the iso-contour of stoichiometric mixture (b); the stoichiometric mixture iso-contour is coloured by the temperature.

(18)

with  
 vxn v vxn ~j a þ Gt a [ad xna ¼ ru vxj vxi vxi sffiffiffiffiffiffiffiffi
 Gt vxna dWin [wien xna ¼ r 2 r vxi Dt

Step 4: Chemical reaction

(19)

xna ðtnþ1 Þ ¼ xnð3Þ þ a

ZDt   u_ j t;

t ¼ 0 : j ¼ jnð3Þ

(23)

0



 ~ a ¼ 1Cf r xna  f ~a [mix xna ; f 2 tsgs Equation (18) is solved using an approximate factorization scheme, in a similar way as in Ref. [15]. Since the stochastic field equation does not consider the energy source term due to the compressibility, an energy correction step is added before the integration of chemical reaction during each time step. The stochastic fields is integrated from time tn to tnþ1 ¼ tn þ Dt via the following four steps:

In step 3, hn and h~ are the enthalpy from the nth stochastic field and the LES field, respectively. In step 4, the chemical reaction rate is integrated using Euler-implicit scheme, with constant internal energy. The Wiener process in step 1 is represented by time-step increment Dt1=2 hni , where hni ¼ f  1; 1g is a dichotomise random vector. The random vectors are generated using a randomly rearrange operation

Step 1: Advection-diffusion-Wiener 

rxnð1Þ ¼ rxna ðtn Þ þ [ad xnð1Þ Dt þ [wien xna ðtn Þ Dt a a |fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflffl} implicit

(20)

Eulerexp licit

Step 2: Micro-mixing

¼ xnð1Þ  0:5Cf xnð2Þ a a

 Dt
nð2Þ  fð1Þ x a tsgs a

(21)

Step 3: Energy correction

hnð3Þ ¼ hnð2Þ þ

N 1 X h~  hmð2Þ N m¼1

! (22)

Fig. 2 e Evaluation of the chemical mechanism on the laminar flame speed of hydrogen/air mixture at 1 atm and 300 K. The experimental data is from Aung et al. [40].

Please cite this article in press as: Gong C, et al., Large eddy simulation of hydrogen combustion in supersonic flows using an Eulerian stochastic fields method, International Journal of Hydrogen Energy (2016), http://dx.doi.org/10.1016/j.ijhydene.2016.09.017

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y x x x ( 2 0 1 6 ) 1 e1 2

Fig. 3 e Instantaneous density gradient (down) compared experimental shadow image (up) in the non-reactive flow.

on N seeds (the fields number N is always set as an even number) with half seeds as “1”, and other half as “1”. This method can provide the mean and the variance of the random vectors equal zero and unity, respectively.

Simulation case and conditions The LES-ESF method described above is applied to simulate reactive and non-reactive flows in a hydrogen fuelled supersonic combustor that has been studied experimentally by Waidmann et al. [22]. A schematic of the combustor geometry is given in Fig. 1a. The combustor consists of a one-side divergent channel and a wedge strut. Vitiated air (23.2% O2, 73.6% N2, 3.2% H2O in mass) enters from the constant area section with the cross section of 50 mm  45 mm. The base of the strut is placed at x ¼ 0. Hydrogen is injected at the centre of the strut base through 15 holes with diameter of 1 mm placed with a uniform distance of 24 mm in the spanwise direction. Since no specific interaction between the neighbouring fuel injection was noticed in a previous numerical simulation [23], only one injector is included in the present simulation, with the periodic boundary condition set on the spanwise direction. A mesh containing around one million hexahedral cells

5

is used to discrete the domain. The cells are clustered towards the wedge wall as well as the wake region downstream of the strut base. The grid resolution for the wake region downstream the strut base is about 0.1 mm  0.1 mm  0.1 mm. In the downstream region, the resolution for the wake region is about 0.5 mm  0.2 mm  0.1 mm. This grid resolution is similar to that used in Ref. [24]. As shown in Fig. 1b, the reaction zone is solved inside the refined region in the present mesh. Adiabatic condition is set on all walls. Since the combustion only takes place around the wake region under the present configuration, no-slip condition is set on the wedge wall, and the slip condition is set on the upper and lower walls. Following the work of Genin and Menon [23] as well as Oevermann [26], the free stream velocity is fixed to 732 m/s, with a static pressure of 100 kPa and a static temperature of 340 K. The ports of fuel injection system are choked. The fuel is assumed to have an injection speed of 1200 m/s, with static pressure and temperature of 100 kPa and 250 K, respectively. A kinetic mechanism with 7 species and 7 reactions [36] is adopted for the chemical reactions. This mechanism was developed for supersonic combustion, and has been successfully used in a number of supersonic combustion studies [23,37e39]. The mechanism is evaluated for the laminar flame speed at 1 atm and 300 K by comparing with the experimental data from Aung et al. [40]. As shown in Fig. 2, the flame speed predicted using the present mechanism agree well with the experiment, except a slight over-prediction in the fuel-rich region with f > 2:5. Following Mustata et al. [12], 8 stochastic fields for the LES simulation are adopted in the present simulation. In their work, small difference was observed between the results with 8 fields and 16 fields, however the computational cost was doubled with 16 fields. The Bilger mixture fraction is used to quantify the mixing in the reactive flow field, which is defined as





YH  YHo 2WH  YO  YOo WO . .  ; Z¼ f f YH  YHo 2WH  YO  YOo WO

(24)

where YH and YO are the element mass fractions; WH and WO are the atomic weights for the elements; the superscripts “o” and “f” indicate the oxidizer stream and the fuel stream,

Fig. 4 e Profiles of mean axial velocity in the non-reactive flow compared with available experimental data and the numerical result from Genin and Menon [23]. Please cite this article in press as: Gong C, et al., Large eddy simulation of hydrogen combustion in supersonic flows using an Eulerian stochastic fields method, International Journal of Hydrogen Energy (2016), http://dx.doi.org/10.1016/j.ijhydene.2016.09.017

6

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y x x x ( 2 0 1 6 ) 1 e1 2

Fig. 5 e Instantaneous density gradient (down) compared with experimental shadow image (up) in the reactive flow.

respectively. The stoichiometric mixture fraction (Zst ) in the present condition is 0.0285. The simulations were carried out in a Linux cluster comprising of 16-core 2.3 GHz XEON processors using 128 cores. Typical CPU times for the simulation with 8 fields were approximately 5000 core-hours for 1.0 ms in physical time, which is approximately 8 times of the CPU time consumption for the simulation with the well-stirred reactor model, WSR.

Results and discussion The LES solver was first validated under a non-reactive flow condition. Fig. 3 shows the density gradient from the LES at an

instance of time, and a Schlieren image from the experiment [22], which also indicates the density gradient in the flow field. The shock waves from the numerical simulations (indicated by the high density gradient) show similarity to that in the Schlieren shadowgraphy. The supersonic incoming flow is deflected at the leading front of the strut, and two shocks are generated symmetrically on each side. These shock waves reflect at the upper/lower wall, respectively. The reflected shock waves interact with the expansion waves generated at the trailing edge of the strut, which generates a series of Mach waves at the downstream field. It is noticeable that due to the boundary layers of the upper and lower walls are not resolved in LES, some differences between the numerical simulation and experiment exist, especially at the reflecting point of the shock waves. This feature has some effect on the position of the downstream pressure waves. However, its influence on the shock wake region of the strut is rather small. To further evaluate the performance of the present LES solver, the predicted mean axial velocity at different axial positions is compared with the measurement using the laser Doppler Velocimetry (LDV) technique. For comparison the numerical results from Genin and Menon [23] (cf. Fig. 4) are also shown. The predicted velocity profile at the first axial location differs from the experiment in the wake region. The present result is similar to that from Genin and Menon [23]. In this region the velocity gradient is high. It is likely that the experiment failed to capture the peak velocity. At the two downstream locations the predicted velocity profiles are shown to be in closer agreement with the experiment as well

Fig. 6 e Time averaged temperature distributions on the central plane from the simulations with different combustion models. Please cite this article in press as: Gong C, et al., Large eddy simulation of hydrogen combustion in supersonic flows using an Eulerian stochastic fields method, International Journal of Hydrogen Energy (2016), http://dx.doi.org/10.1016/j.ijhydene.2016.09.017

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y x x x ( 2 0 1 6 ) 1 e1 2

as the numerical results of Genin and Menon. The shock waves reported earlier have a clear impact on the velocity in the freestream. This feature has been captured in the numerical simulation (cf. Fig. 4aeb). The reactive flow field was then simulated using both the WSR and the ESF model. To investigate the sensitivity of micro-mixing effect in the ESF model, the model constant Cf ¼ 2:0 (the standard value) and 3.0 was used in the simulation, respectively. Fig. 5 shows the density gradient in the reactive flow, and compared with the Schlieren image from the experiment [22]. Due to the heat-release in the wake region, the wake shear layer in the reactive flow expands with a small angle, which is in contrast to the wake shrinking towards the centerline in the non-reactive flow. As a consequence of gas expansion, the wake flow in the reacting case is wider and the deflected shock waves become weaker in the downstream field, with quasi uniform pressure along the vertical direction. The time-averaged temperature distribution on the central plane from LES-WSR and LES-ESF simulations are shown in Fig. 6. It is shown that the fuel injection has very short penetration (~10 mm), and surrounded by mixture with moderately high temperature (~1000 K). Some difference in the mean temperature predicted by the WSR and the ESF models can be

7

identified around the main flame base (15 mm < x < 50 mm). In this region, the LES-ESF model with Cf ¼ 2:0 predicts a lower temperature than that from LES-WSR. The LES-ESF method with Cf ¼ 3:0 predicts the temperature in between of those predicted from LES-WSR and LES-ESF with Cf ¼ 2:0. The time-averaged mean temperature profiles at different axial positions from the LES simulations are shown in Fig. 7 and compared with the experimental data. It is shown that at two downstream position (x ¼ 58 mm and 166 mm), where the turbulent flame is fully developed, all simulations using different models are in reasonable agreement with the experiment. The width of the turbulent flame brush is found to be over-predicted in all simulations, and the ESF method predicted slightly better result of the flame brush width (x ¼ 58 mm). At the first experimental location (x ¼ 11 mm), where local-extinction/re-ignition occurs, the WSR model under-predicted the temperature, while the results from the ESF model with Cf ¼ 2:0 are in reasonable agreement with the experiment. At x ¼ 30 mm, where there is no available experimental data, the ESF method with Cf ¼ 2:0 predicted a lower temperature in the reaction zone than that from the WSR simulation. Since the effect of sub-grid scale turbulence on the combustion is neglected in the WSR model, the difference of temperature predicted by ESF and WSR models

Fig. 7 e Profiles of mean temperature in the reactive flow compared with available experimental data. Please cite this article in press as: Gong C, et al., Large eddy simulation of hydrogen combustion in supersonic flows using an Eulerian stochastic fields method, International Journal of Hydrogen Energy (2016), http://dx.doi.org/10.1016/j.ijhydene.2016.09.017

8

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y x x x ( 2 0 1 6 ) 1 e1 2

Fig. 8 e The distribution of mean heat-release rate (HRR) and streamlines around the flame base from the LES-ESF simulation with Cf ¼ 2:0; the white lines illustrate the shock waves reflecting on the combustion region.

comes from the sub-grid scale fluctuation modelled in the ESF method. At all positions, the result predicted by ESF method with Cf ¼ 3:0 falls in between of the results from the WSR and the ESF model with Cf ¼ 2:0, which indicates that the predicted sub-grid scale effect decreases with the increase of Cf . A comparison with experiment shows that the ESF method with Cf ¼ 2:0 provides the best agreement with the experiment under the present condition. Therefore, the flame structure and stabilization are investigated based on the results from the ESF with the standard model constant Cf ¼ 2:0. To understand the stabilization mechanism of the flame, the streamlines and chemical heat-release rate (HRR) of the time averaged mean flow field around the flame base are plotted in Fig. 8. It is shown that there are two pairs of recirculation zones in the wake. The upstream one is small in size, and it is generated owing to the fuel injection. The fuel mixes with air in the shear layer of both the upstream and the downstream recirculation zones. The downstream recirculation zone is large in size, and it is generated owing to the air flow sudden expansion when passing through the strut. The strut-generated recirculation zone is closed to the base of the strut in the non-reactive flow (cf. Fig. 3), but is pushed downstream in the reactive flow. The main mass flow direction in the recirculation zones is denoted by the thick array in the figure. The mass flow in the strut-generated recirculation zone exchanges with that in the injection-generated recirculation zone. The fuel-rich hot product from the strutgenerated recirculation zone helps igniting the combustible mixture from the injection-generated recirculation zone. Although moderately high mean temperature (~1000 K) is observed in the shear layer of the upstream recirculation zone, the instantaneous temperature distribution in Fig. 9 shows that the flame is not stabilized continuously in this region. Local-extinction can be observed at each side of the shear layer of the upstream recirculation zones (cf. Fig. 9). It appears that the turbulent flame in the present case is stabilized owing to the downstream recirculation zones. Major HRR is in the shear layer of the downstream recirculation zone around the stoichiometric mixture. In particular, the peak HRR is found at the position where the shock wave is reflected in the shear layer of the downstream recirculation zone (cf. Fig. 8). The shock-wave appears to affect the flow field and the flame stabilization. The local-extinction/re-ignition behaviour around the upstream recirculation zones was also observed by Genin and

Menon [23]. In their work, the re-ignition event was thought to be generated by the large scale flow motion from the downstream recirculation zone into the upstream recirculation zone. To understand the re-ignition process, the instantaneous distributions of velocity and mixture fraction around the strut wake region are plotted in Fig. 10. The main recirculation zones can be identified from the velocity field. The distribution of mixture fraction shows that the mixture fraction in the fuel-injection generated recirculation zone is much higher than the stoichiometric value. In the ignition process, the leading edge of the strut-generated recirculation zone moves towards the strut base, which brings the product to the shear-layer near the strut base. The temperature distribution in Fig. 9c shows that the ignition event is first initialized at the fuel-rich mixture inside the recirculation zone. This is consistent with the flow direction of the recirculation zones shown in Figs. 8 and 10. To understand the local-extinction/ re-ignition process, the temperature of stoichiometric mixture as a function of the local scalar dissipation rate at the instants of time of that shown in Fig. 9 is plotted in Fig. 11. Here the scalar dissipation rate is made up of the resolved part and the sub-grid scale part of the mixing field, i.e., c ¼ cres þ csgs ¼ 2:0

Gt  ~2 Cf f VZ þ Z00 2 r tsgs

(25)

00 2 ~ and Zf where Z are the resolved mixture fraction and its variance, respectively. Both steady burning state and transient behaviours can be identified from the scatter plot shown in Fig. 11. The inert

Fig. 9 e Temperature distributions at two couples of inherent instantaneous times from the LES-ESF simulation; the white lines denote the iso-contour of stoichiometric mixture.

Please cite this article in press as: Gong C, et al., Large eddy simulation of hydrogen combustion in supersonic flows using an Eulerian stochastic fields method, International Journal of Hydrogen Energy (2016), http://dx.doi.org/10.1016/j.ijhydene.2016.09.017

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y x x x ( 2 0 1 6 ) 1 e1 2

9

Fig. 10 e The mixture fraction and velocity vectors around the wake region. The vectors are plotted at the mixture with equivalence ratio higher than 0.5. The pink arrays marks the main recirculation zones. The solid black lines denote the isocontour of T ¼ 1000 K. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

mixing state is shown as the low temperature cloud (around 500 K). The steady burning state is shown as the high temperature cloud (above 1500 K). The transient parts in Fig. 11a show two branches of igniting clouds: one at lower scalar dissipation rate, around 1e5 1/s, and one at higher scalar dissipation rate, around 10e40 1/s. As marked by different colours in Fig. 11a, the igniting branch with lower scalar dissipation rate corresponds the flame base on the upper shear-layer at this instant, where the stoichiometric mixture is ignited by the fuel-rich hot product from the strut-generated recirculation zone at x z 50 mm. The branch with higher scalar dissipation rate corresponds to the re-ignition event on the lower side shearlayer. Both of the igniting branches fall to the regimes with lower scalar dissipation rate than the critical quenching scalar dissipation rate. The critical scalar dissipation rate for quenching (ccrt ¼ 227 s1 ) is calculated based on 1D laminar counter-flow burner with the same mechanism as the 3D simulation. Comparison between Fig. 11a and b shows that the ignition branch with higher scalar dissipation rate is still developing in this instants of time. The temperature distributions in Fig. 9c and d illustrate the early evolution of a re-ignition event on the upper side shearlayer. In this time, the igniting branch with high scalar

dissipation rate on the stoichiometric mixture did not appear (cf. Fig. 11c,d). The entrained fuel-rich hot package can ignite the local combustible mixture; and a flame front is formed in the inhomogeneous mixture, which can be illustrated by the scatter plots in Fig. 12. The scatter plot in Fig. 12a shows that the flame front near the strut base is propagating from the fuel-rich mixture towards the stoichiometric mixture. Since the fuel-injection generated recirculation zone cannot provide hot product due to the high equivalence ratio, the flame front near the strut base will be blown off by the high convectional speed in the shear-layer once the leading edge of strutgenerated recirculation zone moves downstream to its mean position. Based on the above discussion, one can conclude that the local-extinction process under the present condition is dictated by the recirculation zones rather than the high scalar dissipation rate. The reaction front propagation in the non-premixed mixture is strongly dependent on the scalar dissipation rate. Fig. 12b shows the resolved and the SGS scalar dissipation rate in the reaction front (defined as T > 700 K). It is shown that at the reaction front the SGS scalar dissipation rate is much higher than the resolved one. As discussed by Mastorakos [41], the propagation speed of flame front increases with scalar dissipation rate in non-premixed mixture. This reveals that the SGS

Please cite this article in press as: Gong C, et al., Large eddy simulation of hydrogen combustion in supersonic flows using an Eulerian stochastic fields method, International Journal of Hydrogen Energy (2016), http://dx.doi.org/10.1016/j.ijhydene.2016.09.017

10

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y x x x ( 2 0 1 6 ) 1 e1 2

Fig. 11 e Scatter plot of temperature vs. scalar dissipation rate at stoichiometric mixture; the data is sampled at the same instantaneous times as Fig. 9; the clouds from upper and lower shear layer at t1 instant is plotted separately using different colours.

Fig. 12 e Scatter plot of the reactive scalars sample in 0 < x < 7 mm at t ¼ t2 in Fig. 9; (a) temperature vs. mixture fraction; (b) scalar dissipation rate vs. mixture fraction with T > 700 K at t ¼ t2. Please cite this article in press as: Gong C, et al., Large eddy simulation of hydrogen combustion in supersonic flows using an Eulerian stochastic fields method, International Journal of Hydrogen Energy (2016), http://dx.doi.org/10.1016/j.ijhydene.2016.09.017

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y x x x ( 2 0 1 6 ) 1 e1 2

scalar dissipation rate modelled by the present ESF method can promote the early development of the flame front shown in Fig. 12a. This feature could explain the under-predicted temperature near the strut base in the WSR simulation, in which the SGS scalar dissipation rate is neglected (cf. Fig. 7a).

[5] [6]

[7]

Conclusions An Eulerian stochastic fields (ESF) method has been developed for large eddy simulation of turbulent combustion in supersonic flow, and validated based on a hydrogen flame combustor with strut as flame-holder. The developed ESF method is shown to be able to capture the shock-wave/ combustion interaction, and the model gives results in closer agreement with the experiment than that from the well-stirred reactor (WSR) model. The flame structure and stabilization mechanism were investigated using the ESF method. Two pairs of recirculation zones are shown in the reaction zone, namely, the fuelinjection generated recirculation zone in the upstream and the strut generated recirculation zone downstream. The fuelinjection generated recirculation zone could not stabilize the flame, owing to the local high equivalence ratio. The turbulent flame is stabilized by the strut generated recirculation zone where the shock-waves from the free-stream region are reflected by the flame and interacting with the flame. Localextinction and re-ignition processes were observed in the shear-layer near the strut base. The oscillations of the recirculation zones were found to be the control mechanism for the local-extinction and re-ignition under the present condition.

Acknowledgements

[8]

[9]

[10]

[11]

[12]

[13]

[14]

[15]

[16]

This work was supported by the Swedish Research Funding Organizations (CeCOST, KCFP and VR). Cheng Gong was sponsored by China Scholarship Council (CSC) with the grant number of No. 201206110003. The computations were performed using the computer facilities provided by High Performance Computing Center North (HPC2N), and Center for Parallel Computers (PDC).

[17]

references

[20]

[1] Tian Y, Yang S, Le J, Su T, Yue M, Zhong F, et al. Investigation of combustion and flame stabilization modes in a hydrogen fueled scramjet combustor. Int J Hydrogen Energy 2016. http://dx.doi.org/10.1016/j.ijhydene.2016.07.219 [in press]. [2] Sun MB, Gong C, Zhang SP, Liang JH, Liu WD, Wang ZG. Spark ignition process in a scramjet combustor fueled by hydrogen and equipped with multi-cavities at Mach 4 flight condition. Exp Therm Fluid Sci 2012;43:90e6. [3] Wang H, Wang Z, Sun M, Qin N. Large-Eddy/Reynoldsaveraged NaviereStokes simulation of combustion oscillations in a cavity-based supersonic combustor. Int J Hydrogen Energy 2013;38:5918e27. [4] Saghafian A, Shunn L, Philips DA, Ham F. Large eddy simulations of the HIFiRE scramjet using a compressible

[18]

[19]

[21]

[22]

[23] [24]

11

flamelet/progress variable approach. Proc Combust Inst 2015;35:2163e72. Pope SB. PDF methods for turbulent reactive flows. Prog Energy Combust Sci 1985;11:119e92. Haworth DC. Progress in probability density function methods for turbulent reacting flows. Prog Energy Combust Sci 2010;36:168e259. Fox RO, Stiles HL. Computational models for turbulent reacting flows, vol. 419. Cambridge: Cambridge University Press; 2003. Koo H, Donde P, Raman V. A quadrature-based LES/ transported probability density function approach for modeling supersonic combustion. Proc Combust Inst 2011;33:2203e10. Koo H, Donde P, Raman V. LES-based Eulerian PDF approach for the simulation of scramjet combustors. Proc Combust Inst 2013;34:2093e100. ~ o L. A field Monte Carlo formulation for calculating the Valin probability density function of a single scalar in a turbulent flow. Flow Turbul Combust 1998;60:157e72. Sabel'nikov V, Soulard O. Rapidly decorrelating velocity-field model as a tool for solving one-point Fokker-Planck equations for probability density functions of turbulent reactive scalars. Phys Rev E 2005;72:016301. ~ o L, Jime nez C, Jones WP, Bondi S. A Mustata R, Valin probability density function Eulerian Monte Carlo field method for large eddy simulations: application to a turbulent piloted methane/air diffusion flame (Sandia D). Combust Flame 2006;145:88e104. Jones WP, Marquis AJ, Vogiatzaki K. Large-eddy simulation of spray combustion in a gas turbine combustor. Combust Flame 2014;161:222e39. Jones WP, Prasad VN. Large Eddy simulation of the sandia flame series (DeF) using the Eulerian stochastic field method. Combust Flame 2010;157:1621e36. Jones WP, Marquis AJ, Prasad VN. LES of a turbulent premixed swirl burner using the Eulerian stochastic field method. Combust Flame 2012;159:3079e95. Dodoulas I, Navarro-Martinez S. Large eddy simulation of premixed turbulent flames using the probability density function approach, Flow. Turbul Combust 2013;90:645e78. Jangi M, Lucchini T, Gong C, Bai X-S. Effects of fuel cetane number on the structure of diesel spray combustion: an accelerated Eulerian stochastic fields method. Combust Theory Model 2015:1e19. Gong C, Jangi M, Bai X-S. Diesel flame lift-off stabilization in the presence of laser-ignition: a numerical study. Combust Theory Model 2015;19:696e713. Jones WP, Marquis AJ, Noh D. LES of a methanol spray flame with a stochastic sub-grid model. Proc Combust Inst 2015;35:1685e91. Jaishree J, Haworth D. Comparisons of Lagrangian and Eulerian PDF methods in simulations of non-premixed turbulent jet flames with moderate-to-strong turbulencechemistry interactions. Combust Theory Model 2012;16:435e63. Jaishree J. Lagrangian and Eulerian probability density function methods for turbulent reacting flows [Ph.D. thesis]. The Pennsylvania State University; 2011. € hm M, Brummund U, Clauss W, Waidmann W, Alff F, Bo Oschwald M. Supersonic combustion of hydrogen/air in a scramjet combustion chamber. Space Technol 1995;6:421e9. nin F, Menon S. Simulation of turbulent mixing behind a Ge strut injector in supersonic flow. AIAA J 2010;48:526e39. Huang Z-w, He G-q, Qin F, Wei X-g. Large eddy simulation of flame structure and combustion mode in a hydrogen fueled supersonic combustor. Int J Hydrogen Energy 2015;40:9815e24.

Please cite this article in press as: Gong C, et al., Large eddy simulation of hydrogen combustion in supersonic flows using an Eulerian stochastic fields method, International Journal of Hydrogen Energy (2016), http://dx.doi.org/10.1016/j.ijhydene.2016.09.017

12

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y x x x ( 2 0 1 6 ) 1 e1 2

[25] Hou L, Niu D, Ren Z. Partially premixed flamelet modeling in a hydrogen-fueled supersonic combustor. Int J Hydrogen Energy 2014;39:9497e504. [26] Oevermann M. Numerical investigation of turbulent hydrogen combustion in a SCRAMJET using flamelet modeling. Aerosp Sci Technol 2000;4:463e80. [27] Wu J, Wang Z, Bai X, Sun M, Wang H. The hybrid RANS/LES of partially premixed supersonic combustion using G/Z flamelet model. Acta Astronaut 2016;127:375e83. [28] Cao C, Ye T, Zhao M. Large eddy simulation of hydrogen/air scramjet combustion using tabulated thermo-chemistry approach. Chin J Aeronautics 2015;28:1316e27. [29] Berglund M, Fureby C. LES of supersonic combustion in a scramjet engine model. Proc Combust Inst 2007;31:2497e504. [30] Wang H, Shan F, Piao Y, Hou L, Niu J. IDDES simulation of hydrogen-fueled supersonic combustion using flamelet modeling. Int J Hydrogen Energy 2015;40:683e91. [31] Pei Y, Hawkes ER, Kook S. Transported probability density function modelling of the vapour phase of an n-heptane jet at diesel engine conditions. Proc Combust Inst 2013;34:3039e47. [32] The open source CFD toolbox available at. 2016. http://www. openfoam.com. [33] Greenshields CJ, Weller HG, Gasparini L, Reese JM. Implementation of semi-discrete, non-staggered central

[34]

[35]

[36]

[37]

[38]

[39] [40]

[41]

schemes in a colocated, polyhedral, finite volume framework, for high-speed viscous flows. Int J Numer Methods Fluids 2010;63:1e21. Kurganov A, Tadmor E. New high-resolution central schemes for nonlinear conservation laws and convectionediffusion equations. J Comput Phys 2000;160:241e82. Van Leer B. Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method. J Comput Phys 1979;32:101e36. Jachimowski CJ. An analytical study of the hydrogen-air reaction mechanism with application to scramjet combustion. 1988. NASA Technica Paper 2791. Baurle RA, Girimaji SS. Assumed PDF turbulence-chemistry closure with temperature-composition correlations. Combust Flame 2003;134:131e48. € bus H, Gerlinger P, Bru¨ggemann D. Scalar and joint scalarMo velocity-frequency Monte Carlo PDF simulation of supersonic combustion. Combust Flame 2003;132:3e24. Potturi A, Edwards JR. LES/RANS simulation of a supersonic combustion experiment. AIAA paper 2012-0611. 2012. p. 2012. Aung K, Hassan M, Faeth G. Flame stretch interactions of laminar premixed hydrogen/air flames at normal temperature and pressure. Combust Flame 1997;109:1e24. Mastorakos E. Ignition of turbulent non-premixed flames. Prog Energy Combust Sci 2009;35:57e97.

Please cite this article in press as: Gong C, et al., Large eddy simulation of hydrogen combustion in supersonic flows using an Eulerian stochastic fields method, International Journal of Hydrogen Energy (2016), http://dx.doi.org/10.1016/j.ijhydene.2016.09.017