Characteristics of reaction zone in a dual-mode scramjet combustor during mode transitions

Aerospace Science and Technology 99 (2020) 105779

Contents lists available at ScienceDirect

Aerospace Science and Technology www.elsevier.com/locate/aescte

Characteristics of reaction zone in a dual-mode scramjet combustor during mode transitions Wubingyi Shen, Yue Huang ∗ , Yancheng You, Lizhe Yi School of Aerospace Engineering, Xiamen University, Xiamen, 361102, China

a r t i c l e

i n f o

Article history: Received 16 September 2019 Received in revised form 15 January 2020 Accepted 11 February 2020 Available online 14 February 2020 Communicated by Kai Liu Keywords: Dual-mode scramjet Combustion modes transitions Combustor shock-train Reaction zones

a b s t r a c t The mechanism of controllable combustion modes (scramjet and ramjet) is one of the most significant challenges for dual-mode ramjet/scramjet engines. The detailed characteristics of the reaction zone during the transitions of combustion modes within the Hyshot II combustor are numerically investigated with 7 species-7 steps hydrogen/air reaction mechanism and the thickened flame model (TFM). In the case of continuous increasing or decreasing equivalence ratios (ERs), the critical ERs for the mode transitions are 0.474 or 0.539, respectively. As the ER increases or decreases, the location of the leading shock of the combustor shock-train moves upstream or downstream, respectively. When the ER continuously increases, the pressure ratio at the leading shock to the outlet of the isolator is greater than when the ER continuously decreases. Inversely, the mass-weighted average Mach number inside the combustor for a continuously increasing ER is less than that for the continuously decreasing ER. In the scramjet mode, the primary combustion reaction zone is located at the rear of the combustor. In the ramjet mode, combustion occurs in the separation zone at the middle of the combustion chamber, while the reaction zone is closer to equilibrium at the critical ER of 0.474 as compared to the scramjet mode. © 2020 Elsevier Masson SAS. All rights reserved.

1. Introduction The dual-mode scramjet engine has the potential to provide airbreathing propulsion for a combined cycle engine over a wide range of flight Mach numbers [1,2]. The dual-mode scramjet can operate in both scramjet and ramjet combustion modes inside its combustor. The scramjet to ramjet mode transition occurs in a dual-mode combustor caused by a pressure rise from the combustion heat addition, which leads to thermal choking of the combustor. The inverse ramjet to scramjet mode transition occurs as the flight Mach number increases, which is characterized by the absence of a pre-combustion shock train. An isolator is placed between the fuel injector and the inlet, which allows the dual-mode concept to be conceived with a fixed-geometry combustor. This serves to house the pre-combustion shock train when the engine is operating away from the purely scramjet mode. The approach to control the operability of a fixed-geometry dual-mode scramjet is to modify the choking throat and pre-combustion shock train by adjusting the fuel injection law. Therefore, understanding the continuous changes in the equivalence ratios on the transient com-

*

Corresponding author. E-mail address: [email protected] (Y. Huang).

https://doi.org/10.1016/j.ast.2020.105779 1270-9638/© 2020 Elsevier Masson SAS. All rights reserved.

bustion is crucial to enable robust control over the operability of the dual-mode combustor. Combustion mode transitions in a dual-mode combustor is a complex, transient, fluid-combustion process that involves thermal choking [3–7], shock-train formation [5,6,8], boundary layer separation [9], flame dynamics [10–12], and turbulence–chemistry interactions [11,12]. Previous experimental investigations have examined the effects of the flow conditions coming from the inlet, the fuel type and mixing characteristics, the range of fuel-air ratios, and the geometry of the combustion chamber on the transition characteristics. Li et al. [3,4] studied the combustion mode transitions for different inlet total pressure conditions and fuel equivalence ratios, as well as changes to the hydrogen and kerosene flow rates on the pressure distribution under the mode transitions. Laurence et al. [5,6] systematically investigated the response of the HyShot II combustor to different ERs. The critical ER where the onset of thermal choking occurs and the high ER conditions leading up to the inlet unstart were identified. Zhang et al. [9,10] found that a high ER slope controlled by strut injectors and the incoming flow Mach number slope as controlled using a Laval nozzle of the alterable throat area would not lead to slope variations for higher pressures. Turner et al. [11] investigated variations in thrust and combustion efficiency of two combustor configurations of different lengths and constant area sections before divergent operations in different combustion modes. Other studies analyzed the flame

2

W. Shen et al. / Aerospace Science and Technology 99 (2020) 105779

dynamics, mixing and reaction zone features during the transition process, and heating addition-triggered transitions in dual-mode combustors using non-intrusive measurement methods [12–15]. Due to the limited number of parameters that can be measured experimentally, high-resolution numerical simulations are a feasible way to study combustion mode transitions. Several numerical studies have been performed to investigate the influences of the injection strategy, fuel equivalence ratio, inlet boundary conditions at the entrance of the isolator, ignition capability, and air throttling on the mode transitions. Huang et al. [16,17] explored the mixing and combustion in a dual-mode combustor with backwards-facing steps under different inlet boundary conditions at the entrance of the isolator, injection strategies, and fuel equivalence ratios on the mode transition. In addition, nonlinear characteristics of the mode transition processes in both a strut-based combustor and cavity-strut-based combustor were analyzed in detail [18,19]. Larsson et al. [20] captured the combustor shock-train for the HyShot II combustor for ER between 0.40 and 0.52 and found it allowed for shorter combustor lengths with a minimal compromise to the combustion efficiency. Cao et al. [21] obtained different combustion modes in a scramjet engine with a two-stage injection by changing the total amount of fuel added or adjusting the fuel distribution between the two injectors. Kouchi et al. [22] demonstrated that the presence of additional igniters on the sidewalls led to an intensive mode over all ERs above 0.4. Tian et al. [23] studied the combustion mode formation and transitions for a dual-mode scramjet with air throttling for different mass fluxes of the throttling air and throttling off time. Other studies analyzed the flame structure, heat release distribution, and shock reflection during mode transitions using one-dimensional and quasione-dimensional models [24–28]. The transitions among different combustion modes as observed through experiments and simulations are often accompanied by hysteresis phenomena that are critical for flight safety and are of great importance. Yu et al. [29] described a class of nonlinear phenomena and proposed the topological invariance for the stability boundaries of mode transitions. Chang et al. [30] and Bao et al. [31] experimentally investigated the hysteresis phenomenon caused by excessive heat release and analyzed its effect on restarting the control. Zhang et al. [32] explained the mechanism of hysteresis as the transition of the pre-combustion shock train in the isolator. Cui et al. [33] developed a mechanical model based on bistability by assuming the unstable positive feedback mechanism of the shock motion in the converging channel of the Laval nozzle can lead to the hysteresis phenomenon. However, many problems concerning mode transitions still remain to be resolved. The Reynolds average Navier–Stokes (RANS) has been used in the majority of numerical studies regarding the parametric effects under the combustion mode transition processes. A few studies have used the large eddy simulation (LES) to better understand the detailed transient fluid-combustion process with very few variable parameters. The application of RANS for combustion in dual-mode scramjets makes it difficult to capture complex flow field characteristics. In addition, it is also prohibitively expensive for LES to resolve wall boundary layers at high Reynolds numbers or to determine the impact of multi-parameter variables. An alternative approach is the detached eddy simulation (DES), which is validated as a viable compromise between accuracy and cost. This paper is to explore the detailed characteristics of the reaction zone during combustion mode transitions in a Hyshot II combustor through numerical investigations with two- and three-dimensional models. The two-dimensional model was used to determine the critical ER of the mode transition, and the three-dimensional calculations analyzed the characteristics of the reaction zone during mode transitions using the improved delayed detached eddy simulation (IDDES).

2. Numerical methods and setup 2.1. Computational approach 2.1.1. Governing equations Since the application in question is a high Mach number flow with heat release, the compressible form of the equations is used. The continuity equation handles the conservation of mass and is defined as:

∂ ρ ∂ ρ ui + =0 ∂t ∂xj

(1)

where ρ is the fluid density, t is time, u i is the velocity in direction j and x j is the spatial vector. The momentum equation is defined as:

∂ ρ ui u j ∂ τi j ∂p ∂ ρ ui + =− + ∂t ∂xj ∂ xi ∂xj where p is the static pressure and given by:



τi j = 2μl S i j −

1 ∂ uk 3 xk

(2)

τi j is the viscous stress tensor,



(3)

δi j

where μl is the molecular viscosity and δi j is the Kronecker delta. The strain rate is defined as:

Sij =



∂u j ∂ ui + 2 ∂xj ∂ xi 1



(4)

The energy equation is given by:

∂q j ∂ τi j u i ∂ ρ e0 ∂ ρ u j h0 + =− + ∂t ∂xj ∂xj ∂xj

(5)

where q j is the heat flux vector. The specific stagnation internal energy, e 0 , and specific stagnation enthalpy, h0 , are defined as:

e0 = h −

p

ρ

h0 = e0 +

1

+ ui ui

(6)

2

p

(7)

ρ

with the specific static enthalpy, h, given by:

h=

Ns 

Y α hα

(8)

α =1

where Y α and hα are the mass fraction and specific enthalpy of species α , respectively, and N s is the number of chemical species. The species transport equation is defined using the species mass fraction:

∂ρu j Yα ∂ ρ V α, j Y α ∂ρYα + = ω˙ α − , ∂t ∂xj ∂xj

α = 1, . . . , N s − 1

(9)

˙ α is the reaction rate for species α . V α , j is the j compowhere ω nent of the diffusion velocity for species α , which can be described using Fick’s law: V α, j Y α = − D α

∂ Yα ∂xj

(10)

where the species diffusion coefficient, D α , is defined as:

Dα =

λ

ρ c p Leα

(11)

with the Lewis number, Le α , assumed to be constant for each species. λ and c p are the thermal conductivity and specific heat capacity at constant pressure of the mixture, respectively.

W. Shen et al. / Aerospace Science and Technology 99 (2020) 105779

3

with k being the mean turbulence kinetic energy per unit mass, defined as:

Table 1 Mechanism with 7 species and 7 reactions [34]. Reaction

A cm3 /mole·sec

b

T

H2 + O2 ↔ OH + OH H + O2 ↔ O + OH OH + H2 ↔ H2 O + H O + H2 ↔ H + OH OH + OH ↔ H2 O + O H + OH + M ↔ H2 O + M H + H + M ↔ H2 + M

1.7 × 1013 1.2 × 1017 2.2 × 1013 5.06 × 104 6.30 × 1012 2.21 × 1022 7.3 × 1017

0

24157.0 8310.5 2591.8 3165.6 548.6 0 0

−0.91 0 2.67 0 −2.00 −1.00

Third body efficiencies for reactions: 2.5 for M = H2 , 16.0 for M = H2 O, 1.0 for all other M.

The chemical kinetic model used in this paper is the H2 /air with 7 species (H2 , O2 , H2 O, O, OH, H, N2 ) and 7 elementary reactions [34]. The Arrhenius coefficients for this kinetic mechanism are given in Table 1. The chemical kinetic model is demonstrated to be sufficient to consider the turbulence/chemistry interactions in supersonic combustion simulations [35].

k=

1 u u 2 i i

(20)

The mean stagnation energy is defined as:

e˜ 0 = h˜ −

p˜

ρ¯

1

+ u˜ i u˜ i + k

(21)

2

The turbulence kinetic energy gradient can be used to approximate the molecular diffusion term in the time averaged energy equation:

  ∂ τi j u i ∂ ∂ k˜ ≈− μl ∂xj ∂xj ∂xj

(22)

and the mean heat flux vector is defined as: s  ∂T + ρ V α , j Y α hα ∂xj

N

2.1.2. Reynolds Averaged Navier Stokes (RANS) The instantaneous governing equations can be averaged in time to produce the Reynolds Averaged Navies Stokes (RANS) equations. The time averaging of a quantity φ is defined as:

φ¯ =

1

t +T φ(τ )dτ

T

(12)

t

where T is a time scale much greater than any physical time scales in the flow. The Reynolds Averaged Navier Stokes equations are obtained when this averaging is applied to the instantaneous system of equations and they describe the evolution of mean quantities. In order to avoid the introduction of extra terms into the continuity equation the standard Reynolds averaging of Eq. (12) can be replaced by Favre averaging, defined as:

φ˜ =

ρφ ρ¯

(13)

where the overbar and tilde indicate Reynolds averaged and Favre averaged quantities, respectively. An instantaneous quantity can be decomposed into a mean and fluctuating contribution:

φ = φ˜ + φ 

(14)

with the fluctuations denoted by  . Applying Favre averaging to the instantaneous system of Eq. (1), (2), (5) and (9), the compressible RANS equations can be obtained:

∂ ρ¯ u˜ j ∂ ρ¯ + =0 ∂t ∂xj u i u j ∂ ρ¯  ∂ ρ¯ u˜ i u˜ j ∂ τ¯i j ∂ p¯ ∂ ρ¯ u˜ i + =− + − ∂t ∂xj ∂ xi ∂xj ∂xj  ∂  ∂ ρ¯ e˜ 0 ∂ ρ¯ u˜ j h˜ o + = τ¯i j u˜ i + τi j u  i − q¯ j − ρ u j h0 ∂t ∂xj ∂xj   ∂ ρ¯ u˜ j Y˜ α ∂ ρ¯ Y˜ α  Y ¯˙ α − ∂ ρ V α , j Y α + ρ¯ u + =ω j α ∂t ∂xj ∂xj

(15)

(16) (17) (18)

The turbulent Reynolds stress term, ρ¯  u i u j , in the momentum equation is often modeled using the concept proposed by Boussinesq, where the turbulent stresses are proportional to gradients in mean velocity:

  1 ∂ u˜ k 2 −ρ¯  u i u j = 2μt S˜ i j − δi j − ρ¯ kδi j 3 ∂ xk 3

(19)

q¯ j = −λ

(23)

α =1

where the heat conduction and energy flux due to gradients in species mass fraction can respectively be approximated as:

λ

∂T ∂ T˜ ≈λ ∂xj ∂xj

Ns 

(24)

ρ V α , j Y α hα ≈ −

α =1

Ns 

ρ¯ D α h˜ α

α =1

∂ Yα ∂xj

(25)

The turbulent flux of stagnation enthalpy,

ρ u j h 0 , is defined as:

 h  + ρ¯ u  k ˜ i ρ u j h 0 = ρ¯ u u i u j + ρ¯ u j j

(26)

The turbulent enthalpy flux and turbulent kinetic energy flux, along with the turbulent diffusion term from the time-average species transport equation, can all be modeled using the gradient diffusion approximation:  h = − ρ¯ u j  k = − ρ¯ u j

μt c p ∂ T˜

(27)

Prt ∂ x j

μt ∂ k˜ σk ∂ x j

 Y = − ρ¯ u j α

(28)

μt ∂ Y˜ α

(29)

Sct ∂ x j

where σk is a model constant. Prt and Sct are the turbulent Prandtl and Schmidt numbers, respectively. The turbulent viscosity μt , is computed using the SpalartAllmaras one equation model [36], which solves a transport equation for the kinematic eddy viscosity v˜ :

∂ ρ¯ u˜ j v˜ ∂ ∂ ρ¯ v˜ + = ∂t ∂xj ∂xj Q = cb1 ρ¯ S v˜ + ρ¯





μl + μt ∂ v˜ + Q, σ ∂xj

cb2 ∂ v˜ ∂ v˜

σ ∂xj ∂xj

− c w1 f w ρ¯

 2 v˜

(30)

d

with d being the wall distance. The turbulent viscosity is defined as:

μt = ρ¯ f v1 v˜

(31)

4

W. Shen et al. / Aerospace Science and Technology 99 (2020) 105779

L IDDES = ˜f d (1 + ˜f e ) L RANS + (1 − ˜f d ) L LES

The auxiliary relations used are:

f v1 =

χ3 χ

3

S = +

+ c 3v1 v˜

k2 d2

f v2 = 1 −

,

c w1 =

f v2 ,

χ , 1 + χ f v1

cb1 k2

+

χ=

cb2 + 1

ρ¯ v˜ , μl

(32)

σ

where is the magnitude of the vorticity of the velocity vector. The wall blockage function for the formulation of the destruction term is:



fw = g r=

1 + c 3w3

1/6

g 6 + c 6w3





g = r + c w2 r 6 − r ,

,

(33)

v˜ max( S , 10−16 )k2 d2

The model constants are:

k cb1 cb2 c w2 0.41 0.1355 0.622 0.3

c w3 2.0

σ

c v1 0.6667 7.1

2.1.3. Large eddy simulation (LES) The equations for LES can be obtained by applying a filter to the instantaneous governing equations. A filtered variable is defined as:

¯ x) = φ(



    φ x∗ G x − x∗ dx∗

(34)

D

where G (x − x∗ ) is the filter function. Finite volume LES can be thought of as using a top-hat filter in physical space. The top-hat filter is a sharp cutoff filter, where the turbulent structures are either big enough to be captured by the mesh, or they are not [37]:



G (x) =

1/ |x| ≤ /2 0 otherwise



(35)

where is the filter width. Similarly to Favre averaging, Favre filtering is used to obtain the compressible form of the governing equations, in order to avoid the introduction of additional terms into the continuity equation. ˜ , where Eq. (12) is also used to define a Favre filtered variable (φ) the overbar now corresponds to the filtering defined in Eq. (34). Although Favre filtering is conceptually different to the Favre averaging process and the resulting system of equations has a different physical meaning, they are visually identical to the RANS equations. The only modification required is the way in which the turbulent viscosity is calculated, with this parameter now referred to as the subgrid-scale (SGS) viscosity (μ S G S ). The Spalart-Allmaras (SA) original model with detached eddy simulation (DES) modification [38] is a hybrid RANS-LES method which makes use of a single turbulence model for both the near wall RANS regions and the freestream LES calculations. The equation of the Spalart-Allmaras original model is still used for the DES modified Spalart-Allmaras model. However, the length scale d defined in Eq. (30) has different definition in the DES model. In the SA-model, d is the distance to the nearest wall. In the DES model, ˜ which is defined as: d is replaced with d,

d˜ = min(d, C DES ),

=

 3

δx δ y δ z

(36)

where δx , δ y and δz are the cell dimensions in the x, y and z directions, respectively. Calibration of the LES-mode for decaying homogeneous isotropic turbulence suggests the value C DES = 0.65. Improved Delayed Detached-Eddy Simulation (IDDES) [39]is developed from the Delayed Detached-Eddy Simulation (DDES) and its hybrid turbulent length scale is given as:

(37)

L LES is written as the grid scale . ˜f d and ˜f e are both blending functions. SA model based IDDES is used to simulate the three-dimensional unsteady reaction flow during the transitions of combustion modes within the Hyshot II combustor. The steady RANS solutions obtained with the SA model are used as initial condition for the IDDES. 2.1.4. Combustion modeling The α species transport equation in RANS or LES can be written as:

  ∂ ρ¯ u˜ j Y˜ α ∂ ∂ ρ¯ Y˜ α ∂ Y˜ α ¯ ˜ ˙ α(Q ) + + =ω ρ¯ D ∂t ∂xj ∂xj ∂xj

(38)

where Q denotes any quantity entering the computation of the reaction rate. The thickened flame model is an extension of the initial model proposed by Butler & O’ Rourke [40] for premixed flames. These authors showed that a correct laminar premixed flame propagation on a coarse grid may be achieved by artificially thickening the flame front. The α species transport equation (38) is rewritten as:

 ˜  ∂ ρ¯ u˜ j Y˜ α ∂ ρ¯ Y˜ α 1 ¯˙ α ( Q˜ ) + ∂ F ρ¯ D ∂ Y α + = ω ∂t ∂xj F ∂xj ∂xj

(39)

As the diffusion and chemical source terms have been multiplied and divided by the thickening factor F , respectively, the thickened flame front propagates at the same laminar flame speed S 0L as the original flame of thickness δ L0 but its thickness becomes F δ L0 . For sufficiently large values of F , the flame structure can be resolved on a coarse grid. 2.1.5. Numerical methods The flow solver is based on a finite-volume cell center scheme. The Navier-Stokes equations is solved with an implicit time stepping technique and the convective fluxes is computed with a HLLC approximate Riemann solver. The inviscid flux function is a secondorder upwind scheme and the diffusive fluxes are approximated using a second-order central scheme. The two-layer wall function is used to solve for the viscous flux at the wall and it is reliably on fine grids as well as on coarse ones (300 > Y+ > 0.1) [41]. The Minmod limiter is utilized to construct a second-order method in space and all final solutions are achieved using second order discretization. A second-order implicit scheme relying on a dual time stepping methodology is applied and the local time-stepping technique defined through a unique CFL number, whereas the global time step must be determined to capture accurately the transient combustion flow. 2.2. Calculation setup A three-dimensional physical model of the HyShot II combustor is composed of a wedge-shaped intake ramp, a rectilinear combustor, and a single expansion exhaust nozzle. The computational domain, boundary conditions, and mesh details are shown in Fig. 1. The isothermal no-slip wall boundary condition is used at the lower and upper walls of the combustor, and the wall temperature is set to 300 K. The injector boundary condition is specified as the total pressure and the total temperature. The exhaust nozzle is the supersonic outlet, and the parameters are interpolated from the internal values of the combustor. Dirichlet boundary conditions are used for all variables at the isolator inlet, with the data derived from the two-dimensional RANS [42].

W. Shen et al. / Aerospace Science and Technology 99 (2020) 105779

5

Fig. 1. HyShot II combustor computational domain, boundary conditions, and mesh details in the region around the fuel injector. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)

Fig. 2. Comparison of numerical results to the experiment data at ER = 0.295 under different spatial resolutions corresponding to 2.8, 12, and 30 million total cell-counts.

The detailed parameters of the isolator inlet are shown in Table 2. The entire computational domain is initialized from the isolator inlet parameters. The turbulent Schmidt number and laminar/turbulent Prandtl number ratio are fixed at 0.7 and 0.89, respectively, and the time step for the calculations is 5 × 10−7 s. The steady-state convergence criterion is that numerical residues drop at least 5 order of magnitude and the difference between the inlet and outlet mass flux is required to drop below 1.5%. The maximum number of inner iterations (local) is set to 20∼30 to guarantee the reduction of at least 2∼3 order of magnitude on the numerical residues. The boundary and initialization conditions of the twodimensional numerical simulation are the same as those in three dimensions. The ER is defined as

Fig. 3. Comparison of numerical density gradient contours in the injector symmetry plane (bottom) with an experimental schlieren image (top).

6

W. Shen et al. / Aerospace Science and Technology 99 (2020) 105779

Table 2 Detailed parameters of the HyShot II combustor for the isolator inlet. Static pressure

Density

Flow velocity

Temperature

Mach number

130.9 kPa

0.328 kg/m3

1801 m/s

1377 K

2.49

N2 mass fraction

O2 mass fraction

NO mass fraction

N mass fraction

O mass fraction

0.752

0.216

0.032

0

0

Fig. 4. Comparison of different combustion model numerical results to experiment data at ER = 0.295.

 ˙ H2 = m ER =

2

γ +1



 2γγ+−11 Pt A

˙ H2 8m ˙ O2 m

γ R Tt

(40)

Table 3 Minimum mesh size and corresponding wall y+ for each computational grid. Grid

Total

(41)

2.3. Validation of calculations The HyShot II combustor experimental data from the high enthalpy tunnel in Gottingen (HEG) in German Aerospace Center (DLR) are used to verify the calculations. Three different grids are used to estimate the influence of the grid resolution on the results, with total cell-counts of 2.8, 12, and 30 million. Fig. 2 shows the pressure and heat flux distributions of the upper and lower walls for different grids, and the numerical results agree well with the experiment data [42]. The minimum mesh size and the corresponding wall Y+ in each computational grid are shown in Table 3. Using 12 million cells with Y+ = 3.5 could accurately simulate the flow field and shock structure while reducing the computational burden. Based on the computational density gradient and experimental schlieren image [6], the shock-train in Fig. 3 is observed to have an acceptable agreement. These comparisons provide a high confidence that the numerical methods can be used to perform the subsequent investigations.

Grid_1 Grid_2 Grid_3 Grid_4 Grid_5

2.8M 12M 12M 12M 30M

Y+

Minimum mesh size (mm)

x

y

z

0.02 0.02 0.02 0.02 0.02

0.01 0.01 0.0075 0.005 0.005

0.03 0.03 0.03 0.03 0.03

5 5 3.5 2 2

The idea in using Thickened Flame Closure is due to its ability to handle chemical effects and dynamics of the all speed combustion processes through the Arrhenius law [43], and TFM is more amenable to the use of RANS, DES, LES. The comparison between the finite-rate model and the thickened flame model is shown in Fig. 4 and higher accuracy is obtained with the thickened flame model compared to the finite rate chemistry. The thickened flame closure approach is implemented in the present IDDES of the transient combustion mode transitions. 3. Results and discussion In the different combustion modes, the Mach numbers for different slices of the combustion chamber and the wall pressure

W. Shen et al. / Aerospace Science and Technology 99 (2020) 105779

7

Fig. 5. Upper (left) and lower (right) wall pressure distributions as the ER increases.

Fig. 6. Upper (left) and lower (right) wall pressure distributions as the ER decreases.

Fig. 7. Mach distributions as the ER increases (left) and decreases (right).

distribution change significantly. Therefore, it is possible to use the pressure change in the central slice of the combustor to distinguish between the scramjet or ramjet modes [44,45]. The computational costs are unacceptable for direct use of three-dimensional numerical simulations to determine the potential critical equivalence ratio of the modal transitions. To simplify the calculations, two dimensions could be used to determine the essential characteristics of the combustion chamber and the critical ER of the mode transition during the direct-reverse continuous changes of the ER from 0.122 to 0.903. Based on the identified critical ER of the mode transition, the characteristics for the combustion reaction zone of the scramjet combustor near the critical equivalence ratio (scramjet mode, ramjet mode, and critical transition process) are studied using three-dimensional high-resolution simulations with IDDES. 3.1. Two-dimensional simulations The pressure distributions of the upper and lower walls as the ER continuously increases are shown in Fig. 5. As the ER increases,

Fig. 8. Pressure ratio (P2 /P1 ) and mass average Mach number inside the combustor for increasing/decreasing ER.

8

W. Shen et al. / Aerospace Science and Technology 99 (2020) 105779

Fig. 9. Comparison of the hydrogen mass fraction and density gradient distributions between ER = 0.295 and 0.539.

the pressure of the combustion chamber rises due to the released heat. Increases in the combustion-induced pressure at the isolator can propagate along the boundary-layer separation as driven by the shock-train structure. It is seen in Fig. 5 that the pressure quickly rises on the lower and upper walls for ER = 0.474. When the equivalence ratio changes from 0.474 to 0.539, the average Mach number over the entire combustion chamber drops from 1.09 to 0.98. As the equivalence ratio continues to increase from 0.539 to 0.903, it is observed that the pressure distributions on the upper and lower walls remain unchanged, and the combustion mode transitions from scramjet into ramjet. The pressure distributions on the upper and lower walls as the ER continuously decreases are shown in Fig. 6. As the ER decreases, the pressure of the entire combustor drops, while the boundary layer separation phenomenon within the isolator and the shock wave generated by the action of the upstream isolator gradually disappears. Fig. 6 reveals that the pressures on the upper and lower walls decrease sharply and fluctuate at ER = 0.539. When the equivalence ratio changes from 0.659 to 0.539, the average Mach number over the entire combustion chamber rises from 0.98 to 1.1. The shock-train in the isolator disappears when ER = 0.34. As the ER changes from 0.474 to 0.122, it is seen that the pressure distributions on the upper and lower walls of the HyShot II combustor remain unchanged. As the ER increases or decreases, the position of the leading shock moves upstream or downstream, respectively. The pressure of the leading shock of the combustor shock-train at x = −0.0136 m on the lower wall is denoted as P1 , and the pressure at the isolator outlet (x = −0.001 m) is P2 . The positions of P1 and P2 are labeled in Fig. 3. It is seen from Fig. 7 and 8 that the pressure ratio (Mach number) for an increasing ER is greater (less than) than values for the decreasing ER. For a continuously increasing or decreasing ER, the critical ER of the mode transition is 0.474 or 0.539, respectively. Correspondingly, with an increased ER, the pressure ratios at the leading shock

and isolator outlet are greater than for a decreased ER. Compared with the variations of the pressure ratio, the variations in the mass-weighted average Mach number for the combustor are exactly opposite. 3.2. Three-dimensional simulations The critical parameters used to describe the transverse jet are the velocity ratio u j /u 0 and the momentum flux ratio ρ j u 2j /ρ0 u 20 . With a small velocity ratio, the transverse jet transforms into a mixed jet and wake stream. Based on the jet velocity ratio, a general jet trajectory formula [46] is proposed as:

     ρ j u 2j   ρ j u 2j β  y/ D ) = α x/( D 

ρ0 u 20

ρ0 u 20

(42)

For ER = 0.295, α is 1.6, β is 1/3, and D is the injector diameter. The hydrogen trajectory of the transverse jet and density gradient distributions are shown in Fig. 9. It is seen from Fig. 8 that the hydrogen trajectory obtained from the empirical formula is consistent with the RANS results. The shock waves within the combustion zone are derived from the shock-trains of the isolator and the bow shock wave caused by the transverse jet. The shock waves are reflected between the wall and the jet in the combustion chamber. When u j /u 0 < 3, the jet shear layer appears to be wrinkled, and the flow is globally unstable [47,48]. The hydrogen trajectory differs from the theoretical value determined from Eq. (42), which causes the shock reflection to become unstable. By comparing the density gradient distributions, it is seen that the shock reflection is affected primarily by the hydrogen trajectory. The flow is divided into the recirculation zone (R-Z), injector mixing zone (I-M), self-ignition zone (S-I), combustion zone (C-Z), and shock reflection zone (S-R). The choking of the transverse jet causes a bow shock wave (B-W) in the jet upstream. Under the

W. Shen et al. / Aerospace Science and Technology 99 (2020) 105779

9

Fig. 10. Contours in the injector plane, from x = −0.025 m to 0.1 m at ER = 0.355 and 0.539.

bow shock wave compression, the high-temperature zone forms behind the shock, which ignites and stables the flame. The interactions of the bow shock wave with the upstream boundary layer form a reverse pressure gradient that results in a separation of the boundary layer and produces a Lambda shock wave (L-W) and a low-velocity recirculation zone. The recirculation zone forms an upstream separation vortex, which helps enhance fuel mixing in the recirculation zone upstream of the jet. Subsequently, the mixture gas self-ignites under the influence of the shock-train (S-W) from the isolator. The shock wave is constantly reflected in the shock reflection zone, which is between the walls of the combustion chamber and the trajectory of the hydrogen gas. The first collision point of the isolator shock-train with the hydrogen trajectory is called the CP1 point. The second (CP2 ) and first collision points of the bow shock wave and the hydrogen trajectory is called the CP3 point. Fig. 10 shows the temperature, Takeno Flame Index (TFI) [49], and density gradient distribution as functions of time for ER =

0.355 and 0.539. A positive TFI indicates the flame is premixed, and a negative TFI means the flame is non-premixed. The TFI is negative in the injector plane, which suggests the flame is dominantly non-premixed [49]. With an increased ER, the backpressure in the combustion chamber increases, the separation of the boundary layer upstream of the transverse jet intensifies, and a new Lambda wave appears on the upper wall. The angle of the bow shock wave increases from 29.18◦ to 37.34◦ . The combustor shock train merges with the bow shock wave, and the reaction zone of the combustion chamber advances. Thus, the shock reflection area is reduced, and the area of the non-premixed flame in the combustion chamber increases significantly. The subsonic region first appears near the lower wall downstream of the combustion chamber. As the ER increases, the subsonic region appears near the upper wall and expands forwards, which eventually results in a thermal chock and mode transition. The development of the subsonic region is observed in Fig. 11.

10

W. Shen et al. / Aerospace Science and Technology 99 (2020) 105779

Fig. 11. Contours of the instantaneous Mach number in the injector plane (sonic line is colored white).

Fig. 12. Scatter plots of temperature and YH2O at x = 117 mm. Points within 0.5 mm from the lower wall are colored red, points within 0.5 mm from the upper wall are colored blue, the left points are colored black, and the equilibrium flame line is colored orange.

The distribution of the mixture fraction Z, temperature, and

affected by the boundary layer, and the cooling effects of the upper

mass fraction of H2 O are shown in the slice at x = 117 mm in

wall are not apparent at ER = 0.295. There is a small backpressure

Fig. 12. Fig. 12(a) shows that the region near the walls is primarily

increase in the combustion chamber, and the combustion zone oc-

W. Shen et al. / Aerospace Science and Technology 99 (2020) 105779

curs mainly in the middle of the combustor. With an increased ER, the combustor backpressure due to additional heat leads to the boundary layer separation, while the combustor shock-train moves forward as shown in Figs. 9 and 10. The subsonic region appears at the rear of the combustor and moves forwards. This is also seen in Fig. 11, which shows how the subsonic region moves for a changing ER. The fuel burns primarily in the boundary layer separation region and the middle part of the combustor, causing the mixture fraction of the part near the upper wall to increase. At the same time, the wall cooling affects the burnt gases so that the mass fraction of H2 O and the temperature distributions in Fig. 12 present a certain gradient. The temperature of the combustion chamber gradually approaches the equilibrium flame temperature, and the H2 O distribution is consistent with the equilibrium flame value in Fig. 12(b). The subsonic region in the combustor and combustor shock-train move forwards until they merge with the fuel injector bow shock wave at ER = 0.539. Due to the presence of thermal choking (shown in Fig. 6), the combustion is completed at the front of the combustion chamber. This causes a large amount of H2 O to accumulate at the rear of the combustion chamber. In Fig. 12(c), the temperature of the combustion chamber is greater than the equilibrium flame temperature because of the complete reaction products. Therefore, the mixture fraction decreases, which indicates the fuel mixing efficiency controls the combustion. 4. Conclusion The numerical results presented in this paper confirm the primary qualitative findings for the transitions of combustion modes and wall pressure oscillations from previous experiments. For a continuously increasing or decreasing ER, the critical ER for the mode transition is 0.474 or 0.539, respectively. In the scramjet mode, the main combustion reaction zone is at the rear of the combustor. In the ramjet mode, combustion occurs in the separation zone and in the middle of the combustion chamber. Therefore, the reaction zone is closer to equilibrium at the critical ER of 0.474 compared to the scramjet mode. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements This research is sponsored by National Natural Science Foundation of China [grant numbers 51876182 and 51406171], the Foundation of Science and Technology on Scramjet Laboratory [grant number 614270302060417], the Fundamental Research Funds for the Central Universities of China [grant number 20720180058], and the Aeronautics Power Foundation [grant number 6141B090325]. References [1] W. Huang, Z. Du, L. Yan, R. Moradi, Flame propagation and stabilization in dualmode scramjet combustors: a survey, Prog. Aerosp. Sci. 101 (2018) 13–30. [2] W. Huang, L. Yan, J. Tan, Survey on the mode transition technique in combined cycle propulsion systems, Aerosp. Sci. Technol. 39 (2014) 685–691. [3] J. Li, G. Jiao, J. Luo, W. Song, Effects of total pressure on mode transition in a dual-mode combustor, Acta Astronaut. 155 (2019) 55–62. [4] J. Li, D. Shen, Q. Fu, Y. Wang, W. Song, Mode transition of fuel control test in a dual-mode combustor, Appl. Therm. Eng. 111 (2017) 1312–1319. [5] S.J. Laurence, D. Liebera, J.M. Schramm, K. Hannemann, J. Larsson, Incipient thermal choking and stable shock-train formation in the heat-release region of a scramjet combustor. Part I. Shock tunnel experiments, Combust. Flame 162 (2015) 921–931. [6] S.J. Laurence, S. Karl, J.M. Schramm, K. Hannemann, Transient fluid-combustion phenomena in a model scramjet, J. Fluid Mech. 722 (2013) 85–120.

11

[7] C. Zhang, J. Chang, S. Feng, W. Bao, D. Yu, Investigation of performance and mode transition in a variable divergence ratio dual-mode combustor, Aerosp. Sci. Technol. 80 (2018) 496–507. [8] B. Xiao, C. He, J. Xing, S. Zhang, Experimental and numerical investigations of freejet and direct-connect dual-mode scramjet, Aerosp. Sci. Technol. 75 (2018) 297–303. [9] C. Zhang, J. Chang, J. Ma, W. Bao, D. Yu, J. Tang, Effect of Mach number and equivalence ratio on the pressure rising variation during combustion mode transition in a dual-mode combustor, Aerosp. Sci. Technol. 72 (2018) 516–524. [10] C. Zhang, J. Chang, J. Zhang, W. Bao, D. Yu, Effect of continuous Mach number variation of incoming flow on ram-scram transition in a dual-mode combustor, Aerosp. Sci. Technol. 76 (2018) 433–441. [11] J.C. Turner, M.K. Smart, Mode change characteristics of a three-dimensional scramjet at Mach 8, J. Propuls. Power 29 (2013) 982–990. [12] C. Aguilera, K.H. Yu, Scramjet to ramjet transition in a dual-mode combustor with fin-guided injection, Proc. Combust. Inst. 36 (2017) 2911–2918. [13] R.D. Rockwell, C.P. Goyne, H. Chelliah, J.C. McDaniel, B.E. Rice, J.R. Edwards, L.M.L. Cantu, E.C.A. Gallo, A.D. Cutler, P.M. Danehy, Development of a premixed combustion capability for dual-mode scramjet experiments, J. Propuls. Power 34 (2018) 438–448. [14] J. Ye, D. Shi, W. Song, G. Li, Z. Zhang, Z. Hu, Investigation of turbulence flow characteristics in a dual-mode scramjet combustor using hydroxyl tagging velocimetry, Acta Astronaut. 157 (2019) 276–281. [15] Y. Tian, X. Zeng, S. Yang, F. Zhong, J. Le, Experimental study on the effect of equivalence ratio and injector position on flow structure and flame development in the scramjet combustor, Aerosp. Sci. Technol. 82–83 (2018) 9–19. [16] W. Huang, Q. Li, L. Yan, L. Liao, Numerical exploration of mixing and combustion in a dual-mode combustor with backward-facing steps, Acta Astronaut. 127 (2016) 572–578. [17] W. Huang, L. Yan, Numerical investigation on the ram–scram transition mechanism in a strut-based dual-mode scramjet combustor, Int. J. Hydrog. Energy 41 (2016) 4799–4807. [18] L. Yan, L. Liao, W. Huang, L. Li, Nonlinear process in the mode transition in typical strut-based and cavity-strut based scramjet combustors, Acta Astronaut. 145 (2018) 250–262. [19] L. Liao, L. Yan, W. Huang, L. Li, Mode transition process in a typical strut-based scramjet combustor based on a parametric study, J. Zhejiang Univ. Sci. A 19 (2018) 431–451. [20] J. Larsson, S. Laurence, I. Bermejo-Moreno, J. Bodart, S. Karl, R. Vicquelin, Incipient thermal choking and stable shock-train formation in the heat-release region of a scramjet combustor. Part II. Large eddy simulations, Combust. Flame 162 (2015) 907–920. [21] R. Cao, J. Chang, W. Bao, M. Guo, J. Qin, D. Yu, Z. Wang, Analysis of combustion mode and operating route for hydrogen fueled scramjet engine, Int. J. Hydrog. Energy 38 (2013) 5928–5935. [22] T. Kouchi, G. Masuya, T. Mitani, S. Tomioka, Mechanism and control of combustion-mode transition in a scramjet engine, J. Propuls. Power 28 (2012) 106–112. [23] Y. Tian, S. Yang, J. Le, Numerical study on effect of air throttling on combustion mode formation and transition in a dual-mode scramjet combustor, Aerosp. Sci. Technol. 52 (2016) 173–180. [24] K. Nordin-Bates, C. Fureby, S. Karl, K. Hannemann, Understanding scramjet combustion using LES of the HyShot II combustor, Proc. Combust. Inst. 36 (2017) 2893–2900. [25] C. Fureby, On the supersonic flame structure in the hyshot II scramjet combustor, in: 26th ICDERS Conference, 2017. [26] D.J. Micka, S.M. Torrez, J.F. Driscoll, Heat release distribution in a dual-mode scramjet, in: 16th AIAA/DLR/DGLR International Space Planes and Hypersonic Systems and Technologies Conference, 2009. [27] L. Tian, L. Chen, Q. Chen, F. Li, X. Chang, Quasi-one-dimensional multimodes analysis for dual-mode scramjet, J. Propuls. Power 30 (2014) 1559–1567. [28] J. Ma, J. Chang, J. Zhang, W. Bao, D. Yu, Control-oriented unsteady onedimensional model for a hydrocarbon regeneratively-cooled scramjet engine, Aerosp. Sci. Technol. 85 (2019) 158–170. [29] D. Yu, T. Cui, W. Bao, Catastrophe, hysteresis and bifurcation of mode transition in scramjet engines and its model, Sci. China, Technol. Sci. 52 (2009) 1543–1550. [30] J. Chang, L. Wang, W. Bao, Q. Yang, J. Qin, Experimental investigation of hysteresis phenomenon for scramjet engine, AIAA J. 52 (2014) 447–451. [31] W. Bao, Q. Yang, J. Chang, Y. Zong, J. Hu, Dynamic characteristics of combustion mode transitions in a strut-based scramjet combustor model, J. Propuls. Power 29 (2013) 1244–1248. [32] Y. Zhang, S. Zhu, B. Chen, X. Xu, Hysteresis of mode transition in a dual-struts based scramjet, Acta Astronaut. 128 (2016) 147–159. [33] T. Cui, S. Tang, C. Zhang, D. Yu, Hysteresis phenomenon of mode transition in ramjet engines and its topological rules, J. Propuls. Power 28 (2012) 1277–1284. [34] R.A. Barurle, S.S. Girimaji, Assumed PDF turbulence-chemistry closure with temperature-composition correlations, Combust. Flame 134 (2003) 131–148. [35] C.J. Jachimowski, An analytical study of the hydrogen-air reaction mechanism with application to scramjet combustion, Tech. Rep., NASA, 1988, p. 2791.

12

W. Shen et al. / Aerospace Science and Technology 99 (2020) 105779

[36] P.R. Spalart, S.R. Allmaras, A one-equation turbulence model for aerodynamic flows, Rech. Aérosp. 1 (1994) 5–21. [37] U. Piomelli, J.R. Chasnov, Large-eddy simulation of turbulent flows, in: Large Eddy Simulation and Related Techniques: Theory and Applications, Von Karman Institute Lecture Series, 2008. [38] P.R. Spalart, Young-Person’s Guide to Detached-Eddy Simulation Grids, Tech. Rep., NASA, 2001, CR-2001-211032. [39] M.L. Shur, P.R. Spalart, M.K. Strelets, A.K. Travin, A hybrid RANS-LES approach with delayed-DES and wall-modeled LES capabilities, Int. J. Heat Fluid Flow 29 (2008) 1638–1649. [40] T.D. Butler, P.J. O’ Rourke, A numerical method for two-dimensional unsteady reacting flows, Proc. Combust. Inst. 16 (1977) 1503–1515. [41] B.E. Launder, On the computation of convective heat transfer in complex turbulent flows, J. Heat Transf. 110 (1988) 1112–1128. [42] S. Karl, Numerical Investigation of a generic scramjet configuration, Ph.D. Thesis, Technical University of Dresden, 2011. [43] O. Colin, F. Ducros, D. Veynante, T. Poinsot, A thickened flame model for large eddy simulations of turbulent premixed combustion, Phys. Fluids 12 (2000) 1843–1863.

[44] Y. Zhang, S. Zhu, G. Liu, X. Li, X. Xu, An overview on mode transition in dual mode ramjet, J. Propuls. Technol. 34 (2013) 1719–1728. [45] B. Xiao, J. Xing, Y. Tian, X. Wang, Experimental and numerical investigations of combustion mode transition in a direct-connect scramjet combustor, Aerosp. Sci. Technol. 46 (2015) 331–338. [46] E.F. Hasselbrink, M.G. Mungal, Transverse jets and jet flames. Part 2. Velocity and OH field imaging, J. Fluid Mech. 443 (2001) 27–68. [47] S. Narayanan, P. Barooah, J.M. Cohen, Dynamics and control of an isolated jet in crossflow, AIAA J. 41 (2003) 2316–2330. [48] A.R. Karagozian, Transverse jets and their control, Prog. Energy Combust. Sci. 36 (2010) 531–553. [49] H. Yamashita, M. Shimada, T. Takeno, A numerical study on flame stability at the transition point of jet diffusion flames, Proc. Combust. Inst. 26 (1996) 27–34.