gas-pilot flame scramjet

CJA 873 21 June 2017 Chinese Journal of Aeronautics, (2017), xxx(xx): xxx–xxx

No. of Pages 18

1

Chinese Society of Aeronautics and Astronautics & Beihang University

Chinese Journal of Aeronautics [email protected] www.sciencedirect.com

4

Numerical simulations of turbulent flows in aeroramp injector/gas-pilot flame scramjet

5

Bing Chen a,*, Xu Xu a, Baoxi Wei a,b, Yan Zhang a

3

6 7

8

a b

School of Astronautics, Beihang University, Beijing 100083, China Science and Technology on Scramjet Laboratory, Beijing Power Machinery Institute, Beijing 100074, China

Received 8 April 2016; revised 24 January 2017; accepted 16 March 2017

9

11 12

KEYWORDS

13

Combustion mode; Laminar flamelet model; Menter’s SST k-x model; Scramjet; Turbulent combustion

14 15 16 17

Abstract To uncover the internal flow characteristics in an ethylene-fueled aeroramp injector/gaspilot (ARI/G-P) flame scramjet, a Reynolds-averaged Navier-Stokes (RANS) solver is constructed under a hybrid polyhedral cell finite volume frame. The shear stress transport (SST) k-x model is used to predict the turbulence, while the Overmann’s compressibility corrected laminar flamelet model is adopted to simulate the turbulent combustion. Nonreactive computations for Case 1 (G-P jet on), Case 2 (ARI jets on), and Case 3 (both ARI and G-P jets on) were conducted to analyze the mixing mechanism, while reactive Cases 4–7 at equivalent ratios of 0.380, 0.278, 0.199 and 0.167 respectively were calculated to investigate the flame structure and combustion modes. The numerical results are compared well to those of the experiments. It is shown that the G-P jet plays significant role in both the fuel/air mixing and flame holding processes; the combustion for the four reactive cases takes place intensively in the regions downstream of the ARI/G-P unit; Cases 4 and 5 are under subsonic combustion mode, whereas Cases 6 and 7 are mode transition critical and supersonic combustion cases, respectively; the mode transition equivalent ratio is approximately 0.20. Ó 2017 Production and hosting by Elsevier Ltd. on behalf of Chinese Society of Aeronautics and Astronautics. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/ licenses/by-nc-nd/4.0/).

18

19

1. Introduction

20

Supersonic combustion ramjet (i.e. scramjet) has long been recognized as one of the most suitable systems for the hypersonic propulsion. Of the special consideration in the develop-

21 22

* Corresponding author. E-mail address: [email protected] (B. Chen). Peer review under responsibility of Editorial Committee of CJA.

Production and hosting by Elsevier

ment of a scramjet are the efficient mixing and combustion processes, given our limited experience with sustained hypersonic propulsion. Hydrocarbon-fueled scramjet has received considerable attention in recent years due to the high volumetric energy density, low cost, and relative simplicity of operation over hydrogen-fueled one. However, on account of the longer residence time required for mixing and completion of chemical reactions for hydrocarbon fuels, there are still several challenges in the development of a high-performance scramjet. Air and fuel must mix at a molecular level before combustion, so turbulent mixing and combustion are at the heart of scramjet operation. Recently, numerical tools are playing a more and more important role in the predictions of this kind of combus-

http://dx.doi.org/10.1016/j.cja.2017.06.010 1000-9361 Ó 2017 Production and hosting by Elsevier Ltd. on behalf of Chinese Society of Aeronautics and Astronautics. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Please cite this article in press as: Chen B et al. Numerical simulations of turbulent flows in aeroramp injector/gas-pilot flame scramjet, Chin J Aeronaut (2017), http:// dx.doi.org/10.1016/j.cja.2017.06.010

23 24 25 26 27 28 29 30 31 32 33 34 35

CJA 873 21 June 2017

No. of Pages 18

2 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97

tion flows.1–5 One of the main focuses in these research activities for modeling the turbulent combustion has been the complex interaction between turbulence and chemical reaction. Due to a very short fuel residence time in the combustor, the flame stabilization mechanisms are usually governed by autoignition. Thus, detailed chemical kinetics is usually required to accurately model the ignition and extinction phenomena.3–5 In this case, the traditional species transported finite-rate kinetics method becomes useless for several reasons. First, when a detailed chemical kinetics mechanism is adopted, there will be a large amount of species involved in the species transport equations, which will evoke a vast amount of computational work even beyond the state-of-the-art computer hardware capabilities. Second, to get the chemical source terms in the species transport equations, the Arrhenius law is often used with the turbulence-chemistry interaction ignored, which may cause severe errors. To account for this interaction, some researchers used the eddy dissipation concept (EDC) model6, linear eddy mixing (LEM) model7, transported probability density function (PDF) model8, etc. Third, the coupled Reynolds-averaged Navier-Stokes (RANS) equations and species transport equations are very stiff and difficult to solve accurately due to the strong nonlinearity of the source terms as well as the wide range of time scales associated with both chemistry and turbulence. One alternative approach to calculate the vast number of species involved hydrocarbon/air turbulent combustion flow is the laminar flamelet model, and there has been extensive research regarding this aspect.1–5,9–15 This model successfully separates the chemical time scale and turbulence time scale, allowing the chemistry to be solved independently before the combustion flow computation. The chemistry results can be stored in a tabulated form as a function of a limited number of indexing scalars. During the real-time combustion simulation, just these scalars are needed to be determined in addition to the RANS calculation, while the species transport equations are not needed to be solved anymore. Because the number of the indexing scalars is independent of chemical mechanism, the computational effort is not proportional to the species number. This is very attractive in the detailed chemical kinetics involved combustion flow calculation. The laminar flamelet model was originally established for the low Mach number turbulent combustion flow, and has been successfully applied to the simulation of turbulent diffusion flames in the subsonic flow. In the supersonic turbulent combustion flow, due to complex flow patterns such as shock waves, contact discontinuities and flame fronts, the model’s applicability is questioned by many experts. The key doubtful point is whether the thickness of the flamelet is really smaller than that of the Kolmogorov vortices or not. But according to the investigations conducted by Bray et al.16, Waidmann et al.17, and Williams18, the combustion in a typical scramjet was approximately in the flamelet regime. More recently, Fan et al.19 conducted a careful theoretical and quantitative comparison of these scales in a scramjet combustor and argued that: (A) the flamelet model assumption is valid for the premixed combustion in the recirculation zones and/or the shear layers for all flight Mach numbers; (B) the assumption is also accurate for the non-premixed combustion for most of the flight Mach numbers except for very high Mach number under which the slow reaction mode dominates the combustion. Moreover, several authors have reported successful computa-

B. Chen et al. tions for scramjet turbulent combustions with this model in recent years.1–5,14,15,17,19 All these studies confirmed that the laminar flamelet model could model the supersonic turbulent combustion flows very well. The aim of this paper is to develop a parallel finite volume RANS/flamelet computational fluid dynamics (CFD) code for supersonic turbulent combustion flow and validate the code against the experimental measurements on an ethylene (C2H4)-fueled aeroramp injector/gas-pilot (denoted as ARI/ G-P) flame dual-mode scramjet combustor developed early this decade at Beihang University (BUAA).20–22 Wei et al.20,21 and Zhang et al.22 experimentally studied the operation properties and mode transition influencing factors. It is worth doing careful three-dimensional (3D) CFD simulation to further clarify the flow structure, fuel/air mixing mechanism, and turbulent flame structure, especially to uncover the flow details in the vicinity of the ARI/G-P injection unit and the intensive combustion region. Therefore, the other more important objective of the present study is to numerically investigate the internal flow characteristics in the ARI/G-P combustor. In the RANS/flamelet simulation, the ethylene/air chemical model23 at the University of California at San Diego (UCSD) was used to generate the flamelet table with the FlameMaster24 code. The ratio of chemical reaction characteristic time scale to turbulence characteristic time scale was evaluated. Numerical results were analyzed with the main emphasis focusing on the mixing mechanism, flame structure, and combustion mode judgment. The present paper is organized in the following. Descriptions on the flamelet combustion model, flow governing equations, and numerical algorithm are given in tion 2. In Section 3, the ARI/G-P scramjet combustor configuration and the computation setups (such as computational grids, modeling of G-P jet, chemical reaction model, and boundary conditions) are presented. In Section 4, the numerical results of three nonreactive cases, i.e. Case 1 (only G-P jet on), Case 2 (only ARI jets on at equivalent ratio of u = 0.380), and Case 3 (both ARI and G-P jets on at u = 0.380), are given to analyze the mixing mechanism and the role of the G-P jet. In Section 5, the numerical results of four reactive cases (Cases 4–7) at u = 0.380, 0.278, 0.199 and 0.167 respectively are presented in detail. Finally, some conclusions are drawn in Section 6.

98

2. Computational methodology

140

2.1. Laminar flamelet model

141

The main assumptions of the laminar flamelet model are as follows1–5,9–15:

142

(1) The turbulent flame can be regarded as a statistical ensemble of steady laminar diffusion flame structures (named as flamelets) embedded in the turbulent flow field. (2) Each flamelet can be approximated as a 1D structure (Fig. 1) with respect to the mixture fraction according to the laminar diffusion flamelet concept presented by Peters et al.9–11 (3) All chemical time scales are short compared to the smallest turbulent time scale, and the combustion is in equilibrium relative to the turbulence; the thickness scale

144

Please cite this article in press as: Chen B et al. Numerical simulations of turbulent flows in aeroramp injector/gas-pilot flame scramjet, Chin J Aeronaut (2017), http:// dx.doi.org/10.1016/j.cja.2017.06.010

99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139

143

145 146 147 148 149 150 151 152 153 154

CJA 873 21 June 2017

No. of Pages 18

Numerical simulations of turbulent flows in aeroramp injector/gas-pilot flame scramjet

3

where the subscripts F and O represent the fuel and oxidizer stream, respectively. The FlameMaster24 code is used to solve Eqs. (1)–(4) to get the two-parameter indexed instantaneous laminar flamelet tables, i.e. Yi(Z, v) and T(Z, v). To account for the turbulence effects, the assumed b- and d-PDF models9–11,25 are used to get the Faver-averaged quantities. Subsequently, the tables are e the variance parameterized by the mean mixture fraction Z, 002 e of the mixture fraction Z , and the mean scalar dissipation rate ~ v as ei ¼ Y ei ð Z; e Z e002 ; ~ Y vÞ; Fig. 1

155 156 157 158 159 160 161 162

One-dimensional structure of laminar flamelet.

165 166 167 168 169 170 171 172 173 174 175

of a flamelet is short compared to the scale of the smallest turbulent eddies (i.e., the Kolmogorov vortices), and the inner structure of the flamelet cannot be destroyed by turbulence. (4) The turbulence influences the flamelets by transporting and twisting them in space; the interaction between turbulence and chemistry can be accounted for with the assumed PDF models.

177 178

180 181

183

The 1D structure of a laminar flamelet is illustrated in Fig. 1 where Z is the mixture fraction, a general mixing variable representing the relative amount that each inflow stream contributes to the local mixture. Z is also a conserved scalar representing other scalars such as the total enthalpy or a particular element mass fraction. Zst is the stoichiometric mixture fraction. n is the unit normal direction vector to the flamelet surface. The governing equations for the flamelet structure can be obtained from the full 1D species transported NavierStokes equations using the low-Mach number asymptotic arguments as9–13 1 @ 2 Yi þ x_ i ¼ 0; qv 2 @Z2

Ns v @2T X hi x_ i   ¼0 2 2 @Z qcp i¼1 Ns X h¼ Yi hi ðTÞ

187 188 189 190 191 192 193 194 195 196 197

199

ð1Þ ð2Þ

ð3Þ

i¼1

184

186

i ¼ 1; 2; . . . ; Ns

p ¼ qRu T

Ns X Yi i¼1

Wi

ð6Þ

202 203 204 205 206 207 208 209 210 212

11

163 164

e Z; e Z e002 ; ~ Te ¼ Tð vÞ

200 201

ð4Þ

where p is the pressure, T is the temperature, q is the density, cp is the specific heat capacity at constant pressure, h is the enthalpy of the mixture, Yi, hi and Wi represent the mass fraction, enthalpy and molecular weight of species i respectively, Ru denotes the universal gas constant, x_ i is the chemical reaction source of species i which can be determined by the chemical kinetics, v is the scalar dissipation rate whose reciprocal represents the time scale of the molecular diffusion, and Ns is the total number of species. The boundary condition for Eqs. (1)–(4) defined in the mixture fraction space is  Z ¼ 0; Yi ¼ Yi;O ; h ¼ hO ð5Þ Z ¼ 1; Yi ¼ Yi;F ; h ¼ hF

where the three indexing parameters can be determined by e and Z e002 , and the tursolving the transporting equations of Z bulence model. To account for the compressibility, only the species mass fractions in the flamelet table are used and the local temperature is calculated implicitly from the total energy1–5 in the present context. In other word, the flamelet table is used as a complicated equation of state describing the relation of the local temperature to the compositions, density, pressure, etc. for the compressible CFD solver.1–5,14,17 Assumption (3) is the key point behind the flamelet model approach. In the present research, the ARI/G-P combustor has a relatively low simulated flight Mach number of about 5.0, under which the assumption is believed to be valid.19 For the non-premixed flame case, the chemistry time scale sc can be defined as the time scale of the main exothermic reactions in the chemical kinetics system, and can be calculated with10 sc ¼

Z2st ð1

 Zst Þ vq

213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231

2

ð7Þ

233

where vq is the corresponding extinction scalar dissipation rate. The turbulence time scale can be stood for by the Kolmogorov vortices time scale26, sg, m1=2 sg ¼ ð8Þ e

234

where m is the kinematic viscosity and e the turbulent kinetic dissipation. The ratio of the chemistry time scale to the Kolmogorov time scale is defined as the Karlovitz number Ka, sc Ka ¼ ð9Þ sg

240

If Ka is less than 1.0, the key flamelet assumption is valid.

246

2.2. Governing equations

247

2.2.1. Mean flow equations

248

The mean flow governing equations consist of the conservation equations of mass, momentum and energy, and the transporting equations for the mixture fraction and its variance. In the Faver-averaged form, the system of equations are expressed as

249

@U þ r  ½FðUÞ  GðUÞ ¼ S @t

ð10Þ

where t is the time, and U, FðUÞ, GðUÞ and S are the mean conservation variables, inviscid and viscous fluxes, and source terms, respectively. These vectors are given by

Please cite this article in press as: Chen B et al. Numerical simulations of turbulent flows in aeroramp injector/gas-pilot flame scramjet, Chin J Aeronaut (2017), http:// dx.doi.org/10.1016/j.cja.2017.06.010

235 236 237 239

241 242 243 245

250 251 252 253 255 256 257 258 259

CJA 873 21 June 2017

No. of Pages 18

4

B. Chen et al. h

iT

261

e q e002 ; q u~1 ; q u~2 ; q u~3 ; q e~t ; q  Z; Z U¼ q

262 264

FðUÞ ¼ F1 e1 þ F2 e2 þ F3 e3 ¼ Fl el

ð12Þ

265 267

GðUÞ ¼ G1 e1 þ G2 e2 þ G3 e3 ¼ Gl el

ð13Þ

268 269 270 271 272 273

275 276 277 278 279 280 282 283

285 286 288 289

291 292

294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312

314

ð11Þ

 is the Reynolds-averaged density, el (l = 1, 2, 3) the where q unit vector of the Cartesian coordinate xl (l = 1, 2, 3), ~ u ¼ u~l el the Faver-averaged velocity vector, and e~t the Faveraveraged specific total energy (including the chemical energy). The inviscid flux Fl and viscous flux Gl can be expressed as 2

3 u~l q 6q 7 6 u~l u~1 þ pdl1 7 6 7 6q 7 6 u~l u~2 þ pdl2 7 6 7 u~l u~3 þ pdl3 7; Fl ¼ 6 q 6 7 6 7 u~l h~t q 6 7 6 7 e u~l Z q 4 5 e002  u~l Z q

2

3 0 6 7 sl1 6 7 6 7 6 7 s l2 6 7 6 7 sl3 Gl ¼ 6 7 6 7 6 slj u~j  ql þ -l ðkÞ 7 6 7 6 7 e -l ð ZÞ 4 5

ð14Þ

e002 Þ -l ð Z

where dij = ei s_ ej (i,j = 1, 2, 3) is the Kronecker operator, and p and h~t are the Reynolds-averaged pressure and Faveraveraged total enthalpy per unit mass, respectively. The viscous stress tensor sij, heat flux ql and operator -l are given by   @ u~i @ u~j 2 @ u~l 2 kdij  q sij ¼ ðl þ lt Þ þ  dij ð15Þ @xj @xi 3 @xl 3 

 ~ X  Ns  ei j l @h j ~ @Y  Di  hi ql ¼  þ t þ q @xl cp Prt @xl i¼1 cp;l

ð16Þ

  lt @k -l ðkÞ ¼ l þ rk @xl

ð17Þ

  e l @Z e ¼ q D þ t -l ð ZÞ Sct @xl

ð18Þ

  e002 l @Z e002 Þ ¼ q D þ t -l ð Z Sct @xl

ð19Þ

where Prt and Sct are the turbulent Prandtl number and Schmidt number, respectively, rk is the Prandtl number for k which is the turbulence kinetic energy, and lt is the turbulence viscosity. Turbulent parameters can be determined by the (shear stress transport) SST k-x turbulence model.27–32 Note that the second term on the right-hand side in the ql expression vanishes under the unity Lewis number (Lel = 1) assumption, so the species mass fractions do not need to be explicitly computed. The conductivity ji and viscosity li of species i are determined with the Eucken’s formula24 and Lennard-Jones formula24, respectively. The mixture properties such as the viscosity l, heat conduction coefficient j, and diffusion coefficient D are computed from the Wilke’s mixing rules.24 The specific heat capacity cp,l and enthalpy hl are expressed by the polynomial functions of the temperature as follows.24 e002 transport equation has a source In Eq. (10), only the Z term component as e @Z e l @Z ~v s7 ¼ 2 t   q Sct2 @xj @xj

ð20Þ

where Sct2 is the turbulent Prandtl number for Z. The scalar dissipation rate, ~ v, measures the decay of the mixture fraction fluctuations, and plays for Z the equivalent role of the turbulent kinetic energy dissipation rate e for the velocity. The assoe002 =~ ciated time scale Z v is often set to be proportional to the turbulence time scale k/e, which leads to the following expressing:10 e e002 ~ v ¼ Cv Z ð21Þ k

315 316 317 318 319 320 321 322 324

where Cv is a constant of order unity. The general equation of state is adopted to close Eq. (10). All the model constants are given in Table 1.

325

2.2.2. Turbulence model

328

The turbulence model is very crucial for the numerical simulation in the scramjet engine, especially when transverse injections are involved.27,28 The SST k-x turbulence model29 is one of the most widely used two-equation Boussinesq concept based turbulence models in internal flow simulations. Moreover, the success of this turbulence model used with flamelet model had been proved in Refs. 2–5. Therefore, the SST k-x model was used in this paper. This model came from the combination of the k-x model and k-e model using a switching function designed by Menter. Consequently, the model is identical to the k-x model inside the boundary layer and the k-e model elsewhere. Moreover, the Bradshaw’s assumption, i.e. the principal shear stress is proportional to the turbulent kinetic energy, was used in the deriving of the model. Hence, the model can account for the transport of the turbulent shear stress in turbulent boundary layer, and is very applicable to separated flows in the adverse pressure gradient boundary-layer. Full descriptions of the original SST k-x model can be found in Ref. 29. Some modifications reported in the literatures are adopted in the present research. The shear strain rate magnitude |S| is used in the turbulence eddy viscosity expression instead of that of the vorticity, |X|. Consequently, the new expression is given by30

329

lt ¼

a1 qk maxða1 x; jSjF2 Þ

327

330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353

ð22Þ 355

where x is the specific dissipation rate for k, |S| the magnitude of velocity strain rate, a1 = 0.31 a model constant, and F2 a blending function given in Ref. 29. This is because there are full of strong shock/boundary-layer interactions (SBLIs), and subsequent separation phenomena in a scramjet combustor, while the original SST k-x model has sometimes yielded over-prediction of the separated flow region where strong SBLIs occur.30 The modification makes the model be more consistent with the Boussinesq eddy-viscosity assumption, and achieve rotational invariance, which seems to be more rational for SBLIs evoked separation flows. The compressibility correction26,31 and low Reynolds number correction26 are also used presently. When the SST k-x turbulence model2–5,27–32 is used, the mean scalar dissipation rate can be finally modeled as e002 ~ v ¼ Cv b  x Z

326

ð23Þ

where the model constant is b* = 0.09. The Kolmogorov time scale can be further expressed as19,26

Please cite this article in press as: Chen B et al. Numerical simulations of turbulent flows in aeroramp injector/gas-pilot flame scramjet, Chin J Aeronaut (2017), http:// dx.doi.org/10.1016/j.cja.2017.06.010

356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 373 374 375

CJA 873 21 June 2017

No. of Pages 18

Numerical simulations of turbulent flows in aeroramp injector/gas-pilot flame scramjet Table 1 Prt

Sct

Sct2

rk

Cv

0.90

0.5

0.50

1.0

2.0

376

378 379 380 381 382

384 385 386 387 388 389 390 391

392

393 394 395

Constants used in model.

 sg ¼

l b qkx

1=2 ð24Þ

Furthermore, the turbulent Reynolds number Rel is used to check whether the flow is a fully developed turbulence or not. It is defined as19,26 pffiffiffiffiffi uðlÞl q 2kl ¼ ð25Þ Rel ¼ m l pffiffiffi where the integration length scale l ¼ k=x, u(l) is the corresponding velocity scale of the energy-containing eddies. It is obvious that Rel is the ratio of the inertial force to the viscous force. Therefore, those vortices with Rel greater than 1.0 are unstable and will gradually break down into smaller scale eddies until the balance between these two forces is established when the eddy scale reaches the Kolmogorov scale. 2.3. Numerical algorithms The computational domain is meshed with unstructured polyhedral cells, i.e., a cell element can be a tetrahedron, wedge, pyramid or hexahedron in the grid system for 3D case. The

Fig. 2

5

lower-upper symmetric Gauss–Seidel (LU-SGS) implicit algorithm33 is employed as the time integration method, and the Harten-Lax-van Leer Contact (HLLC) scheme34 as the spatial discretization method. The state parameters on either side of a cell face are reconstructed with a second-order interpolation in conjunction with the Venkatakrishnam slope limiters35–37 to maintain monotonicity. The cell center gradients are computed by the Green-Gauss theorem, whereas the normal gradients at a cell face are approximated as in Refs. 5,35,36. The SST k-x turbulence model27–32 is solved separately after each timestepping of the RANS equations. In the implicit integration process for both the turbulence and fluid motion equations, the local time step method is used to determine the time step Dt for each cell.

396

3. ARI/G-P model and computation setup

410

3.1. ARI/G-P combustor model description

411

The ARI/G-P combustor sketched in Fig. 220–22 consists of 6 segments with a rectangular cross-section. The top and bottom walls are parallel to each other with a height of 32 mm, while the other two side walls are divergent at angles of 1.14°, 6°, 1.4°, 3.2°, 1.2° and 12° for each segment. The entrance and exit dimensions are 32 mm  54 mm and 32 mm  145 mm, respectively. Segment 1 is used as an isolator. The integrated ARI/G-P injection unit, consisting of an aeroramp injector (ARI) and a gas-pilot (G-P) flame, is located on the bottom wall, and plays as the injection, ignition, and flame holding

Schematic of ARI/G-P combustor.20–22

Please cite this article in press as: Chen B et al. Numerical simulations of turbulent flows in aeroramp injector/gas-pilot flame scramjet, Chin J Aeronaut (2017), http:// dx.doi.org/10.1016/j.cja.2017.06.010

397 398 399 400 401 402 403 404 405 406 407 408 409

412 413 414 415 416 417 418 419 420 421

CJA 873 21 June 2017

No. of Pages 18

6

B. Chen et al. Table 2 air.

Free stream conditions of ethylene jets and vitiated

Parameter

Ethylene jet

Air stream

Mach number Total temperature (K) Total pressure (MPa) C2H4 mass fraction H2O mass fraction O2 mass fraction N2 mass fraction Turbulence intensity Turbulence viscosity ratio Equivalent ratio Momentum flux ratio

1.0 300

2.0 1200 0.87 0 0.080 0.232 0.688 0.04 0.1

1.0 0 0 0 0.05 0.02

450

451

3.2. Computation setup

452

3.2.1. Computational grids

453

An unstructured storage multi-block hexahedron grid system was used for accuracy consideration, although the present CFD code is based on the unstructured hybrid grid. Only half of the model was taken into consideration. Three computational grids, consisting of 2.5 million (coarse), 4.78 million (medium), and 7.66 million (fine) cells, respectively, were used to verify grid convergence. The medium grid is presented in Fig. 3. The coordinate origin was set at the centroid point of the first row of aeroramp injectors, the x-, y- and z-coordinate were along the combustor main flow direction, the spanwise direction, and the direction normal to the ARI/G-P unit wall, respectively. Each ARI injector and G-P

423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449

454 455 456 457 458 459 460 461 462 463 464

Parameters of G-P flame.

Parameter

Value

Diameter (mm) Momentum flux ratio Total temperature (K) Total pressure (MPa) Fuel mass flow rate (g/s) Oxidizer mass flow rate (g/s)

5 1.3 3000 0.9 5 12

0.1–0.5 0.7–3.5

device. The ARI is composed of four flush-wall sonic fuel jets arranged in a two-by-two matrix pattern. The jets are set at prescribed transverse (pitch) angles of 20° and 40°, and toein (yaw) angles of 15° and 30°. The ARI unit has an equivalent diameter deq = 2.8 mm, which is defined as deq = djm0.5, where dj = 1.4 mm is the diameter of each injector and m is the total number of the injectors. The space between injectors is 4deq in the streamwise direction and 2deq in the spanwise direction. The G-P flame, a kerosene/oxygen combustion jet with a diameter of dg = 5 mm, is placed 8deq downstream of the ARI jet center. Several reactive and nonreactive experiments on this scramjet had been performed in the direct-connect supersonic combustion facility at BUAA.20–22 In the reactive tests, the high enthalpy vitiated air was heated by a hydrogen-oxygen heater, entering the isolator with a Mach number of Ma = 2.0, and a static pressure of 0.115 MPa which is used as the reference pressure as p0. The ethylene (C2H4) was injected into the combustor at a Mach number of 1.0 and total temperature Tt of 300 K. The G-P chamber pressure was about 0.9 MPa. The side wall centerline values of pressure were measured. Detailed conditions of the ethylene jets, vitiated air, and G-P flame are summarized in Tables 2 and 3. In the nonreactive experiments, both the ARI ethylene jets and G-P flame jet were replaced with nitrogen at the same momentum flux ratios as those in the reactive cases. Wall oil flow visualization and Schlieren images in the ARI/G-P adjacent region were obtained. More information regarding these ARI/G-P combustor ground experiments can be found in Refs. 20–22.

422

Table 3

Fig. 3

Computational grid near ARI/G-P jets, medium grid.

injector were modeled as cylinders with lengths of Li = 3.57deq and Lg = 2.86deq, respectively. A five-block butterfly-topology was adopted for each fuel injector, and a four-block one for the half G-P injector. Grid cells were refined in the wallnormal direction. Three grids were partitioned into subdomains for the message-passing-interface (MPI)-based parallel computations by MeTiS38 which is an open source software package for partitioning unstructured grids based on a multilevel graph partitioning algorithm.

465

3.2.2. Modeling G-P flame jet

474

Previous research had uncovered that the G-P jet plays several key roles for the ARI/G-P combustor20–22:

475

(1) It enhances the fuel/air mixing remarkably. The G-P jet lifts up the ARI plume greatly, which leads to an obviously large fuel penetration height, big plume area, and subsequently high mixing efficiency. (2) It acts as an igniter. The high enthalpy G-P gas, as the kerosene/oxygen combustion product with a total temperature of 3000 K, is a continuous pilot torch flame for the C2H4/air combustion. (3) It also operates as a flame holder. There are several recirculation zones generated by the interactions among the ARI and G-P jets. The one pilot torch plus several recirculation zones flow pattern is the root of flame holding capability of the ARI/G-P injection unit.

477

466 467 468 469 470 471 472 473

476

478 479 480 481 482 483 484 485 486 487 488 489 490

Theoretically, all properties of the G-P gas should be simulated accurately. Unfortunately, the laminar flamelet model can only cope with a two-stream combustion system.1–5,9–15

Please cite this article in press as: Chen B et al. Numerical simulations of turbulent flows in aeroramp injector/gas-pilot flame scramjet, Chin J Aeronaut (2017), http:// dx.doi.org/10.1016/j.cja.2017.06.010

491 492 493

CJA 873 21 June 2017

No. of Pages 18

Numerical simulations of turbulent flows in aeroramp injector/gas-pilot flame scramjet 494 495 496 497 498 499 500 501 502 503 504 505

506 507 508 509 510 511 512 513 514 515 516 517 518 519 520

To overcome this obstacle, the G-P flame was simply modeled as a hot stream with the same species composition of the vitiated air. Therefore, the chemical composition and high enthalpy effects of the G-P gas on the chemical mechanism were neglected. To the best of our knowledge, the overall simulation accuracy does not deteriorate. There are two known supporting factors. First, the G-P flame mass flow rate is relatively small compared to that of either the fuel or the main air stream. Second, previous experiments have demonstrated that the overall combustion efficiency is not sensitive to the G-P momentum flux ratio21, illustrating that the G-P stream does not influence the chemical process significantly. 3.2.3. Chemistry mechanism and flamelet table The UCSD ethylene/air chemistry model23, consisting of 34 reactions with 19 species, was used in the flamelet table generation. Here, the Lel = 1 condition for all species was assumed. The background pressure was chosen at approximately 1.5  105 Pa. The temperature was set at 730 K for the oxidizer stream and 300 K for the fuel stream. The table included 56 laminar flamelets with the scalar dissipation rate varying from v1 = 0.1 s1 to vq = 1061.0 s1, which led to an upper bound of 30 K for the discrepancy between the peak temperatures from flamelet to flamelet. The laminar flamelet table was further convoluted with the assumed b-PDF and d-PDF25 applied to Z and v, respectively. Then the final chemical species were subsequently tabulated with respect to three mean variables, e Z e002 and ~v. i.e. Z,

521

3.2.4. Boundary conditions

522

The flow at the ARI/G-P combustor entrance is supersonic, so all flow variables were fixed under the free stream condition of the vitiated air. Mass flow inlet conditions were applied for both the ARI and G-P jets. The entrance boundary-layer had been ignored. Supersonic extrapolation condition was adopted at the combustor outlet boundary. Mirror reflection condition was implemented on the symmetry boundary. On the wall, no-slip velocity, adiabatic temperature, and zero-gradient pressure conditions were imposed. The density on the wall was obtained from the equation of state. The mixture fraction Z was fixed at 0 at the entrances of both the air stream and G-P jet, and at 1.0 at the ARI injector entrances. Extrapolation condition for the mixture fraction variance was applied on open flow boundary. On the wall, e002 ¼ 0 were the zero-gradient extrapolation Z condition and Z imposed. The turbulent kinetic energy was fixed to a turbulence intensity of Tu = 4% for the air stream, and Tu = 5% for both the ARI and G-P jets. The turbulence kinetic energy dissipation rate was based on the integration length scale l0 = k1/2/(b*x) using the corresponding hydraulic diameter. Then, the specific dissipation rate can be calculated through x = e/(b*k). On the wall, kw = 0, and

523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545

7

4. Nonreactive results

551

Three typical nonreactive cases listed in Table 4 were calculated to investigate the mixing mechanism. Only G-P valve was turned on for Case 1, and ARI valve for Case 2, whereas both valves were open for Case 3. The equivalent ratio was u = 0.380 in Cases 2 and 3. Grid convergence study was performed for Case 3.

552

4.1. Grid convergence

558

The pressure profiles obtained on three grids are presented in Fig. 4. The pressure p increases slightly in the isolator due to the viscous effect in the boundary-layers, and decreases abruptly at around x/deq = 71.4 for the wall expansion from Segment 1 to Segment 2. Another obvious pressure decrease occurs at approximately x/deq = 28.0, where the expansion waves originating from the opposite wall impinge on the side wall. Then, the pressure abruptly increases at about x/deq = 13.5 as a result of the impingement of the merged bow shock. Downstream of the strong shock/boundary interaction region, the pressure turns back rapidly at about x/deq = 41.0, followed downstream by a geometry divergence resulted overall gradually decrease with several wiggles. It is obvious that the pressure discrepancies from the coarse grid through the fine grid are relatively small, indicating that good grid convergence has been achieved. Furthermore, the grid convergence index (GCI) devised by Roache39 and also cited by Wilcox32 was used to evaluate the accuracy regarding the grid resolution. For a given flow property w, the GCI is defined as39

559

GCI ¼ 1:25

jeh j ra  1

ð27Þ

where eh is the relative difference between two grids for w, r the fine-to-coarse grid point ratio, and a the spatial algorithm order.

Table 4

Three typical nonreactive cases.

Case

G-P

ARI

Equivalent ratio

1 2 3

On Off On

Off On On

0.380 0.380

546

548 549 550

xw ¼ 10

6l b qd21

ð26Þ

where xw is the turbulent dissipation rate for k at the wall, and d1 the normal distance of the cell center away from the wall.

Fig. 4

Pressure profiles along divergent side wall centerline, Case 3.

Please cite this article in press as: Chen B et al. Numerical simulations of turbulent flows in aeroramp injector/gas-pilot flame scramjet, Chin J Aeronaut (2017), http:// dx.doi.org/10.1016/j.cja.2017.06.010

553 554 555 556 557

560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 581 582 583 584

CJA 873 21 June 2017

No. of Pages 18

8

B. Chen et al. Table 5

Grid-resolution results at entrance centers of Segments 3 and 5 (r = 1.602, a = 2).

Variable

Point A 3

Density (kg/m ) Pressure (Pa) Temperature (K) x-velocity (m/s) Mach number C2H4 mass fraction H2O mass fraction O2 mass fraction N2 mass fraction

585 586 587 588 589 590 591 592 593 594

Point B

Medium grid

Fine grid

GCI (%)

Medium grid

Fine grid

GCI (%)

0.2784 83479.7 999.4 784.3 1.310 0.2408 0.0607 0.1761 0.5224

0.2794 82847.9 988.4 797.6 1.343 0.2467 0.0603 0.1748 0.5182

0.29 0.61 0.89 1.33 1.96 1.91 0.53 0.60 0.65

0.2556 61302.2 796.7 1017.5 1.827 0.0773 0.0738 0.2142 0.6347

0.2621 61977.5 785.3 1020.2 1.842 0.0758 0.0741 0.2151 0.6350

1.98 0.87 1.15 0.21 0.67 1.57 0.32 0.33 0.03

The GCI values of 9 flow variables at points A (xA/deq = 53.6) and B (xB/deq = 253.6) from the medium grid to fine grid are presented in Table 5. Points A (entrance center of Segment 3) and B (entrance center of Segment 5) locate in the ARI/G-P adjacent and far downstream regions, respectively. The grid point ratio is r = 7.66/4.78  1.602. The numerical method order is approximately a = 2. Table 5 reveals that the maximums of the GCIs are less than 2% at both points. Therefore, it can be claimed that the GCIs for all flow properties remain the same low level throughout the

Fig. 5

flow field, which illustrates that the cold flow computations on the medium grid are of good accuracy.

595

4.2. Results and discussion

597

4.2.1. Flow structure comparison

598

Fig. 5 presents the numerical and experimental Schlieren images in the ARI/G-P adjacent region. The separation streamlines and sonic lines are also plotted. It should be noted

599

Experimental and numerical schlierens, Cases 1–3.

Please cite this article in press as: Chen B et al. Numerical simulations of turbulent flows in aeroramp injector/gas-pilot flame scramjet, Chin J Aeronaut (2017), http:// dx.doi.org/10.1016/j.cja.2017.06.010

596

600 601

CJA 873 21 June 2017

No. of Pages 18

Numerical simulations of turbulent flows in aeroramp injector/gas-pilot flame scramjet

658

that though the ARI jets for Case 1 and G-P jet for Case 2 were closed, the same computational grid as that for Case 3 was used except that the entrances of the ARI jets and G-P jet were treated as no-slip, no-penetration walls for Cases 1 and 2, respectively. Therefore, the four blinded ARI injector tubes act as low length/depth ratio cavities for Case 1, while the blinded G-P tube for Case 2. In Case 1, the under-expansion G-P jet enters the main stream normally, accelerated out of the G-P tube, and generates a barrel shock and Mach disk. Due to the barrel shock blockage, upstream of the G-P jet originates a strong detached bow shock, which is so strong that the downstream flow is subsonic near the bottom wall. The bow shock interacts with the boundary-layer, generating a strong adverse pressure gradient which propagates upstream through the boundary-layer. Consequently, a relatively large scale boundary-layer separation bubble is generated with a separation shock emanating from its leading-edge. The separation shock impinges and merges with the bow shock. When the bow shock hits the top wall, another boundary-layer separation bubble with a leadingedge separation shock attached is produced. A shear layer between the jet plume and surrounding stream is also generated. Consequently, several subsonic bubbles are generated. In Case 2, two bow shock waves are formed due to the blockages of the ARI fuel jets. The strengths are relatively weak due to the streamwise-biasing angled injections. Therefore, the boundary-layer separation bubbles are small, and the attached separation shocks are weak. Some noise can be seen apparently in the G-P exit adjacent region due to the blinded G-P hole effect. Case 3 is a combination of Cases 1 and 2. The bow shocks and separation bubbles evoked by the ARI jets are the same as those in Case 2. Because the air flow has been altered by the ARI jets, the flow pattern is a little different from Case 1 downstream of the ARI jets: (A) the separation bubble upstream of the G-P barrel shock is not clear; (B) the edge between the barrel shock and shear layer is not clear either; (C) the shapes of the subsonic bubbles upstream of the G-P jet, downstream of the barrel shock and Mach disk, are also altered apparently. A reasonable level of agreement has been achieved between the numerical and experimental schlierens generally. However, due to the influence of the turbulence model, two obvious discrepancies exist: (A) the numerical bow shock wave for Case 1 is displaced slightly in the x-direction; and (B) the numerical boundary separation bubble upstream of each jet is slightly larger. The flow structure in the ARI/G-P adjacent region is further illustrated by the wall friction streamlines in Fig. 6. The numerical flow topology is quite similar to that of the experiment. The incoming flow is blocked by the ARI/G-P jets successively, making the streamlines deviate transversely. The numerical reattachment line angle is a little larger. Two pairs of counter-rotating zones are formed just upstream of the G-P jet. These recirculation zones and their successive downstream structures play a positive role not only to the fuel/air mixing, but also to the flame holding of the following reactive cases.

659

4.2.2. ARI/G-P jets mixing mechanism

660

The ethylene mass fraction (YC2H4) and streamlines in 8 x-planes for Case 3 are presented in Fig. 7. In the x/deq = 0

602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657

661

9

Fig. 6 Experimental and numerical wall friction streamlines for Case 3.

plane, the fuel penetration height is relatively small, i.e. z/deq < 1.0. Both jets have a 20° injection angle of pitch and a 15° angle of yaw in the opposite direction, which results in two relatively small counter-rotating vortices, uplifting the fuel plume and enhancing the fuel/air mixing. In the x/deq = 4 plane, another two fuel jets with a 40° angle of pitch and a 30° angle of yaw inject into the flow stream, interact with the first two and intensify the fuel/air mixing further. Fig. 7(d) reveals that two groups of jet plumes mix perfectly in the transverse direction, and merge into a relatively large-scale fuel plume at the mid-plane at x/deq = 7. The scale of the counter-rotating vortices is also enlarged apparently. At the x/deq = 10 plane, the G-P jet normally injects into the combustor at a mass flow rate of 17 g/s, resulting in a relatively large penetration height of z/deq > 8. When the fuel plume interacts with the G-P jet, the fuel plume is highly uplifted and the counter-rotating vortices are pushed aside to the two opposite sides of the G-P jet. Downstream of the G-P jet the two vortices rapidly move close, and the scale is enlarged significantly, during which the fuel plume mixes intensively with the air stream due to the entrainment of the vortices. Fig. 8 presents the ethylene mass fraction for Case 3 in three z-planes at z/deq = 0, 2.857 and 5.714, respectively. It is shown that in the z/deq = 0 plane the fuel jets are blocked and split by the G-P jet, making the fuel plume scale increase rapidly in the y-direction. Downstream of the G-P jet, ethylene quickly fills the wake region. In Figs. 7 and 8, the C2H4 mass fraction iso-lines at YC2 H4 = 1.0  105 and 0.001 are used to show the edge of the fuel plume, whereas the one at YC2 H4 = 0.068 is used to illustrate the ethylene/air stoichiometric location. It is apparent that due to the interactions among the ARI and G-P jets, (A) the fuel plume scale increases continuously in the transverse direction, and (B) the peak fuel concentration decays rapidly along the streamwise direction. To further clarify the G-P jet role, numerical results in four x-planes for both Cases 1 and 2 are given in Fig. 9. The flood contour variable for Case 1 is chosen as the static temperature

Please cite this article in press as: Chen B et al. Numerical simulations of turbulent flows in aeroramp injector/gas-pilot flame scramjet, Chin J Aeronaut (2017), http:// dx.doi.org/10.1016/j.cja.2017.06.010

662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698

CJA 873 21 June 2017

No. of Pages 18

10

B. Chen et al.

Fig. 7

Fig. 8 grid.

699 700 701

C2H4 mass fraction and streamlines in eight x-planes for Case 3, medium grid.

C2H4 mass fraction in three z-planes for Case 3, medium

to identify the hot G-P jet plume approximately for there was not any fuel injected, while the C2H4 mass fraction is used as flood contour variable for Case 2 for convenient comparison

to Fig. 7. It is shown that the G-P jet in Case 1 and ARI jets in Case 2 separately generate a pair of counter-rotating vortices whose rotating direction is consistent to each other. When these two pairs of vortices interact and overlap, the vortices’ structures for Case 3 are generated, and the fuel/air mixing is enhanced. Consequently, the fuel plume area of Case 3 is much larger than that of Case 2, indicating that the mixing efficiency of the ARI/G-P unit is much better than that of the pure ARI unit.

702

5. Reactive results

711

Then, four typical reactive cases in Table 6, i.e. Cases 4–7 at u = 0.380, 0.278, 0.199 and 0.167, respectively, were calculated to demonstrate the flow structure, flame structure, and combustion mode. The corresponding experimental results had been published in Refs. 20–22.

712

Please cite this article in press as: Chen B et al. Numerical simulations of turbulent flows in aeroramp injector/gas-pilot flame scramjet, Chin J Aeronaut (2017), http:// dx.doi.org/10.1016/j.cja.2017.06.010

703 704 705 706 707 708 709 710

713 714 715 716

CJA 873 21 June 2017

No. of Pages 18

Numerical simulations of turbulent flows in aeroramp injector/gas-pilot flame scramjet 717

5.1. Reactive cases at u = 0.380

718

5.1.1. Flow characteristics

719

Fig. 10 presents the static pressure, Mach number, and the sonic lines in the combustor symmetric planes. Due to the heat addition, a relatively large adverse pressure gradient is generated and propagates upstream through the wall boundarylayer, pushing the main air stream aside from the wall in the transverse direction. Consequently, the effective supersonic core flow passage in Segment 2 looks like a convergentdivergent- convergent-divergent (C-D-C-D) tube, and there is a large subsonic flow region in Segments 2 and 3. This subsonic region is much larger than that in Case 3. Downstream of the G-P jet, the Mach number increases due to the combustion heat releases. The flow gets to the thermal choking state at around the Segment 3/4 conjunction station. The sonic lines reveal that the thermal choking line looks like a ‘‘W” aligning

720 721 722 723 724 725 726 727 728 729 730 731 732

11

in the streamwise direction in the z/deq = 5.714 plane. For the geometry compression, weak oblique shock waves originate at the entrance of Segment 5, and interact intensively with the side wall boundary-layers, leading to a complex oblique shock train structure. Fig. 11 further plots the centerline and averaged total temperature Tt,c,Tt,av and Mach number Mac,Maav distributions. The Mach number in the isolator is around 2.0, and the centerline Mach number increases to 2.16 at about x/deq = 55.0 for the geometry expansion. When the air flows through the C-D-C parts (37.2  x/deq  0.2) of the boundary-layer altered effective passage, the Mach number profile exhibits an increase-decrease-increase variation characteristics. Then, the air stream interacts with the ARI and G-P jets successively in the last divergent part where the combustion originates. Due to the combined effects of the geometry divergent, air/jet interactions, and heat release, the flow rapidly

Fig. 9 Flows in four x-planes for Case 1 and Case 2, medium grid. Note: Contour variables: static temperature (Case 1); C2H4 mass fraction (Case 2).

Please cite this article in press as: Chen B et al. Numerical simulations of turbulent flows in aeroramp injector/gas-pilot flame scramjet, Chin J Aeronaut (2017), http:// dx.doi.org/10.1016/j.cja.2017.06.010

733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749

CJA 873 21 June 2017

No. of Pages 18

12 Table 6

B. Chen et al. Four typical reactive cases.

Case

G-P

ARI

Equivalent ratio

4 5 6 7

On On On On

On On On On

0.380 0.278 0.199 0.167

Fig. 11 Centerline and mass-flow-weighted average total temperature and Mach number distributions for Case 4, medium grid.

Fig. 10 Static pressure, Mach number and sonic lines (dash lines) in two symmetric planes for Case 4, medium grid.

777

turns to subsonic, and the centerline Mach number reaches Ma = 0.23 at about x/deq = 24.9. Two wiggles on the Mach number profile exist near the air/jet interaction positions. Downstream of the ARI/G-P unit, the Mach number keeps increasing for heat addition and/or geometry expansion except a slight decrease in Segment 5 for viscous effect. The total temperature distributions indicate that the heat releases intensively in the subsonic region. That is to say the chemical reaction mainly takes place in this subsonic region, and hence subsonic combustion is the leading combustion mode. This coincides with the experimental conclusion drawn in Refs. 20,21. The pressure profiles along the side wall centerline for Case 4 are presented in Fig. 12. In Fig. 12(a), the pressure distributions obtained on three grids agree well with each other, although there are some small discrepancies. Numerical results on the coarse grid exhibit an apparent unreasonable local peak pressure in Segment 5. This is because the wall-normal grid resolution of the coarse grid is too poor to accurately capture the complex SBLIs. Due to space limitation, further GCI analysis on the reactive case is not presented. Generally speaking, a good grid convergence has also been accomplished for the reactive case. Fig. 12(b) reveals that the present numerical results agree well with those of the experiment, and the peak pressure is captured accurately, whereas the original point of numerical pressure rise and the overall peak pressure presented by Wei et al.20 are a little more upstream than those of the experiment.

778

5.1.2. Flame structures

779

Because turbulence stochastically twists and transports each individual flamelet, the instantaneous turbulent flame often

750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776

780

exhibits strong unsteadiness, and is quite complicated in practice. Fortunately, the mean flame as a whole usually remains stationary at a certain fixed position.39 The high hydroxyl (OH) radical concentration regions can be used to stand for the most intensive combustion regions approximately. Fig. 13 gives the numerical OH concentration YOH and total temperature distributions in the symmetric planes. Three iso-lines are e = 0.048, 0.068 and 0.088, respectively, where also plotted at Z e Z st ¼ 0:068 is the stoichiometric mixture fraction. It is apparent that (A) the high OH concentration region is consistent with the high total temperature region (i.e. the flame core); (B) the flame core falls into the relatively narrow range of e 6 0:088 with Z est being the centroid point. 0:048 6 Z The mixedness and flame index40 are used to further investigate the flame structure in the ARI/G-P combustor. The

Fig. 12

Side wall pressure profiles for Case 4.

Please cite this article in press as: Chen B et al. Numerical simulations of turbulent flows in aeroramp injector/gas-pilot flame scramjet, Chin J Aeronaut (2017), http:// dx.doi.org/10.1016/j.cja.2017.06.010

781 782 783 784 785 786 787 788 789 790 791 792 793 794 795

CJA 873 21 June 2017

No. of Pages 18

Numerical simulations of turbulent flows in aeroramp injector/gas-pilot flame scramjet 796 797 798 799

801 802 803 804 805 806 807 808 809 810 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831

mixedness ZFO is a measure of the mixing extent of the unburned fuel with oxygen coming from the surrounding air. It is defined as40 8Y Y O2 < O2 6 YC2 H4 fst fst ZFO ¼ ð28Þ Y O : YC H 2 > YC2 H4 2 4 fst where fst is the stoichiometric oxygen-to-fuel mass ratio. Obviously, ZFO is positive in the fuel-rich region and negative in the fuel-lean region. The absolute value of ZFO increases as the mixedness increases. The flame index (denoted as GFO) is designed to distinguish premixed flame from non-premixed diffusion flame.40 It is defined as the dot-product of the mass fraction gradients of both fuel and oxygen, GFO ¼ rYC2 H4  rYO2

ð29Þ

The sign of GFO relies on the relative supplying directions of fuel and oxygen. The fuel and air come from the opposite direction for a typical non-premixed diffusion flame, whereas both of them come from the same direction for a typical premixed flame. Therefore, the flame index GFO is positive for a premixed flame and negative for a diffusion flame. A high value of |GFO| means a rapid supplying rate of fuel and oxygen by molecular diffusion. The mixedness and flame index distributions in the z/deq = 5.714 plane are presented in Fig. 14. ZFO evenly distributed in [0.03, 0.03]; GFO exponentially distributed in [2000.0, 0), while evenly in (0, 70.0]. Fig. 14(a) reveals that there is a large amount of fuel left unburned in the flow region downstream of the ARI/G-P jets. In Fig. 14(b), the negative flame index is plotted in colored floods with its magnitude being exponentially distributed, while the evenly distributed positive flame index is shown in gray floods. Three iso-lines at GFO = 1, 10 and 20 respectively are also included. It is shown that the flame index is negative in almost all the

Fig. 13 OH mass fraction, total temperature and mixture e = 0.048, 0.068 and 0.088 for Case 4. fraction iso-lines at Z

13

intensive combustion region. Therefore, the flames in the combustor are mainly diffusion flames. Near the ARI/G-P joint jet stream, there are some small scattered positive GFO regions where the flames are premixed flames. These well premixed flame regions are generated by the entraining mechanism39 of the large-scale counter-rotating vortices. Generally speaking, the complex flame structure is illustrated very well by the GFO distributions and non-premixed flame is the pronounced flame type in the combustor. Furthermore, Fig. 15 gives the mixedness and flame index plots in 8 x-planes to clarify the streamwise evolution characteristics. ZFO evenly distributed in [0.03, 0.03]; GFO exponentially distributed in [2000.0, 0), while evenly in (0, 70.0]. Upstream of the G-P jet, a large negative GFO zone with | GFO| > 2000.0 is generated since the fuel injects into the main air stream. This means a rapid supplying rate of fuel and oxygen by molecular diffusion due to their large gradients from the opposite direction. Consequently, a typical non-premixed flame structure is formed. Then, the thin fuel/air shear layer flame interacts with the G-P jet introduced at x/deq = 10.0. Downstream of the G-P jet, due to both the intensified fuel/air mixing and fuel/oxygen combustion consumption, the absolute value of the negative GOF rapidly decreases while the zone length scale gradually increases in the transverse direction. As a result, the fuel/oxygen molecular diffusion rates decrease rapidly. A rapid decrease rate of the ZFO absolute value along the streamwise direction is also clearly exhibited. There are also some small regions with positive GFO in the core flow for the entraining role of the vortex motion. Some small well-mixed regions are also generated in the boundary-layer where the vortices hit the bottom wall; hence, several small premixed flame regions are formed (see Fig. 15(e)–(h)).

Fig. 14 Mixedness, flame index distributions in symmetric planes for Case 4.

Please cite this article in press as: Chen B et al. Numerical simulations of turbulent flows in aeroramp injector/gas-pilot flame scramjet, Chin J Aeronaut (2017), http:// dx.doi.org/10.1016/j.cja.2017.06.010

832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863

CJA 873 21 June 2017

No. of Pages 18

14 864

B. Chen et al.

5.1.3. Characteristic scales

886

est ¼ 0:068 iso-surface are Only the characteristic scales on the Z taken into consideration because the most intensive chemical reaction mainly takes place near the stoichiometric regions. The turbulence Reynolds number Rel and Karlovitz number Ka distributions on the stoichiometric surface are presented in Fig. 16. It is shown that the turbulence Reynolds number Rel is much larger than 1.0 on almost all the surface except a relatively small part near the bottom and top walls where Rel is near 1.0. This implies that the combustion takes place in the fully developed turbulence flow regions for Case 4. As also shown in Fig. 16, the Karlovitz number Ka on the stoichiometric surface is much less than 1.0, which means the chemical time scale is much shorter than the Kolmogorov vortices time scale. This generally further indicates that the combustion for Case 4 is in equilibrium relative to the turbulent time scale; hence the key assumption on flamelet is held. Fig. 16 reveals that the stoichiometric mixture faction surface mainly stays in the large-scale subsonic region shown in Fig. 10, but a relatively small rear part extends into the supersonic region. Therefore, though Case 4 is a subsonic combustion dominating case, the intensive combustion takes place not only in the subsonic region but also in the supersonic flow region.

887

5.2. Other equivalent ratio cases

865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885

888 889

Cases 5–7 were further simulated on the medium grid to study the fuel mass flow rate influence. Fig. 17 presents the OH mass

Fig. 15

fraction floods, sonic lines, and stoichiometric iso-line at e = 0.068 for Cases 4–7. Due to different equivalent ratios Z and hence different amounts of heat release, the flow structures are different from case to case. In all the four reactive cases, the subsonic regions downstream of the G-P jet are enlarged for the heat addition. As a result, the static pressure rapidly increases, and produces an adverse pressure gradient which in turn prompts the generation of subsonic bubbles near the combustor wall, or even evokes the boundary-layer separations. In Case 7, the heat release amount is relatively small, so the resulted adverse pressure gradient is not large enough to generate subsonic flow bubbles near two divergent side walls. However, in Case 6, due to the increase in heat release compared to Case 7, the pressure increment is high enough to generate a small subsonic flow region near each side wall. The flows are still supersonic between the side wall subsonic bubbles and the core subsonic region downstream of the G-P jet. From Case 6 to Case 5, both the side wall subsonic bubbles and the core subsonic flow region are enlarged for more heat release. When fuel mass flux and heat release increases to the level of Case 4, the side wall subsonic bubbles and the core subsonic flow region are enlarged to merge into one large subsonic region. It is shown that the critical equivalent ratio at which the side wall subsonic bubble originates is approximately 0.20, which is consistent with the conclusion drawn in Ref. 21. All the stoichiometric lines for the three low equivalent ratio cases are in the G-P jet downstream subsonic flow regions where the chemical reaction takes place intensively.

Mixedness (left) and flame index (right) distributions in 8 x-planes for Case 4.

Please cite this article in press as: Chen B et al. Numerical simulations of turbulent flows in aeroramp injector/gas-pilot flame scramjet, Chin J Aeronaut (2017), http:// dx.doi.org/10.1016/j.cja.2017.06.010

890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917

CJA 873 21 June 2017

No. of Pages 18

Numerical simulations of turbulent flows in aeroramp injector/gas-pilot flame scramjet

15

Fig. 16 Rel and Ka floods on stoichiometric iso-surface at e = 0.068 for Case 4. Z

Fig. 18

Fig. 17 OH mass fraction floods, sonic lines, and stoichiometric e = 0.068, Cases 4–7. mixture fraction iso-lines at Z

918 919 920 921 922

The flame structures for Cases 5–7 are illustrated by the flame index distributions in Fig. 18. GFO exponentially distributed in [2000.0, 0), while evenly in (0, 70.0]. It is shown that all the three are non-premixed flame dominating cases because the negative GFO flood almost fulfills the combustion

Flame index GFO flood in symmetric planes, Cases 5–7.

regions. There are also some scattered premixed flame regions (denoted by positive GFO). The flame structures are similar to that for Case 4 (Fig. 14), but the combustion zone length scale in the streamwise direction becomes shorter and shorter from Case 4 to Case 7. In Fig. 18(b), the premixed flame zone in the z/deq = 5.714 plane is split into several scattered spots due to the occurrence of the edge-wiggled side wall subsonic bubble (Fig. 17(c)). Fig. 19 further plots the pressure distributions along the side wall centerline for the three low equivalent ratio cases. In general, the numerical results agree well with those of the experiments, although there are some differences. The numerical original points of pressure rise for the three cases are slightly upstream than those of the experiments, but the numerical result of Case 4 agrees well with that of experiment (Fig. 12(b)). It is probably because there might be some kinds of flame lifts which cannot be accurately captured by the steady flamelet model. It also may be the results of the turbulence model influence.

923

5.3. Combustion mode

942

Although the combustion mode has been studied for a long time22,41–43, there is not any generally accepted criterion to classify the combustion mode yet. Different conclusion would probably be drawn according to different classification criteria for a given case. One widely used method is based on the

943

Please cite this article in press as: Chen B et al. Numerical simulations of turbulent flows in aeroramp injector/gas-pilot flame scramjet, Chin J Aeronaut (2017), http:// dx.doi.org/10.1016/j.cja.2017.06.010

924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941

944 945 946 947

CJA 873 21 June 2017

No. of Pages 18

16

B. Chen et al.

Fig. 19 Pressure distributions along divergent side wall centerlines for Cases 5–7.

Fig. 20

948 949 950 951

Mach number profiles along passage centerline.

minimum of the 1D Mach number. Wei et al.20,21 applied this method to the ARI/G-P combustor, and concluded that one case at an equivalent ratio u  0.26 was under supersonic combustion mode (e.g. Cases 6 and 7). Zhang et al.22 pointed out

Fig. 21

Mass-flow-weighted average profiles.

that the conclusion was strongly influenced by the accuracy of the 1D result. Here, to further study the ARI/G-P combustion mode, not only the 1D Mach number but also the static pressure rise position relative to injector43 are taken into consideration. Fig. 20 presents the passage centerline Mach number for Cases 3–7. Only the first three segments of Case 3 are plotted. It is apparent that the reactive Mach number profiles are identical to the nonreactive upstream of the original point to decrease. But the original points for Cases 4–7 are different: the original points for Cases 6 and 7 are both at approximately x/deq = 6.2 even though different amounts of heat are added, whereas, for Cases 4 and 5, the original points are at about x/deq = 39.1 and 14.9, respectively. This indicates that there is a threshold value for the heat amount to let its influence propagate to the upstream area of the ARI injectors, and the critical case is Case 6. Fig. 21 further presents the mass-flow-weighted total temperature, Mach number and static pressure profiles for

Please cite this article in press as: Chen B et al. Numerical simulations of turbulent flows in aeroramp injector/gas-pilot flame scramjet, Chin J Aeronaut (2017), http:// dx.doi.org/10.1016/j.cja.2017.06.010

952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970

CJA 873 21 June 2017

No. of Pages 18

Numerical simulations of turbulent flows in aeroramp injector/gas-pilot flame scramjet

992

the five cases. It is shown that the original points of total temperature rise are almost at the same x-station of about x/deq = 6.2, although the peak total temperature decreases from Case 4 to Case 7. This phenomenon indicates that the combustion (and hence heat release) originates from the same x-position for the four reactive cases. However, the original points of pressure rise are very different: (A) the original points of Cases 4 and 5 are at x/deq = 50.7 and 30.1, respectively, both of which are upstream of the fuel injectors; (B) those of Cases 6 and 7 are almost at the same position of about x/deq = 2.6 which is slightly upstream of the total temperature rise point but still downstream of the fuel injectors. The original points of Mach number decrease are the same as those of the pressure rise for each reactive case. It is obvious that the minimums of the averaged 1D Mach number for Cases 4 and 5 are less than 1.0, while that of Case 7 is larger than 1.0, and Case 6 is a critical case. Both Figs. 20 and 21 indicate that (A) Cases 4 and 5 are under subsonic combustion mode, while Case 7 is a typical supersonic combustion case, and Case 6 is approximately a critical transition mode case; (B) the mode transition equivalent ratio is about 0.20, which is lower than that in Ref. 20.

993

6. Conclusions

994

The flow in the ethylene(C2H4)-fueled ARI/G-P scramjet combustor at BUAA has been calculated with a parallel RANS/Flamelet Code, in which RANS equations, SST k-x turbulence model, and the two scalar transport equations are solved with the LU-SGS time algorithm and HLLC spatial scheme under a hybrid polyhedral cell finite volume frame. The compressibility corrected steady flamelet model with the assumed PDF models to account for the interaction between turbulence and chemistry is used to simulate the turbulent combustion. Three nonreactive and four reactive cases are simulated. Pressure distributions and GCIs calculation on three grids with different resolution have been performed. Numerical results show that good grid convergences have been achieved for both nonreactive and reactive cases. Numerical results (such as wall pressure distributions, Schlieren images and viscous streamlines) are compared well with those of the experiments published before. The mixedness and flame index are calculated to analyze the turbulent flame structure in the ARI/G-P combustor. It is proved that the turbulent flame of four reactive cases is dominated by diffusion flame, although there are some small scattered premixed flame spots embedded in the core flow and bottom wall boundary-layers downstream of the G-P jet. The intensive combustion region is proved to be near the stoichiometric mixture fraction, and falls into the narrow range e = 0.048 to 0.088. The turbulence Reynolds number from Z and Karlovitz number distributions on the stoichiometric isosurfaces indicate that flows in the intensive combustion region are in fully developed turbulent state and the chemical time scale is short compared to the smallest turbulent time scale. It is also shown that the combustion in the ARI/G-P downstream flow regions is intensive for all the four reactive cases, and the G-P jet plays significant role in both the fuel/air mixing and flame holding processes. The passage centerline Mach number and mass-flow -weighted average profiles indicate that Cases 4 and 5 are

971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991

995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029

17

under subsonic combustion mode, whereas Cases 6 and 7 are mode transition critical and supersonic combustion cases, respectively. It is proved that the mode transition equivalent ratio is about 0.20, which is the same as the critical equivalent ratio for the side wall subsonic bubble occurrence. This implies that the mode transition of the combustor takes place when the relative large-scale side wall subsonic bubble occurs.

1030

Acknowledgments

1037

This study was co-supported by the National Natural Science Foundation of China (Nos. 51176003 and 51276007), and the Fundamental Research Funds for the Central Universities of China (No. YWF-15-GFY).

1038

References

1042

1. Oevermann M. Numerical investigation of turbulent combustion in a scramjet using flamelet modeling. Aerospace Sci Technol 2000;4(7):463–80. 2. Gao ZX, Lee CH. A flamelet model for turbulent diffusion combustion in supersonic flow. Sci China: Technol Sci 2010;53 (12):3379–88. 3. Pecˇnik R, Terrapon VE, Ham F, Iaccarino G, Pitsch H. Reynoldsaveraged Navier-Stokes simulations of the HyShot II scramjet. AIAA J 2012;50(8):1717–32. 4. Pecˇnik R, Terrapon VE, Ham F, Iaccarino G. Full system scramjet simulation. CRT Ann Res Brief 2009;33–46. 5. Terrapon VE, Han F, Pecˇnik R, Pitsch H. A flamelet-based model for supersonic combustion. CRT Ann Res Brief 2009;47–58. 6. Chakraborty D, Paul P, Mukunda H. Evaluation of combustion models for high speed H2/air confined mixing layer using DNS data. Combust Flame 2000;121(1–2):195–209. 7. Ge´nin F, Chernyavsky B, Menson S. Large eddy simulation of scramjet combustion using a subgrid mixing/combustion model. In: 12th AIAA international space planes and hypersonic systems and technologies; 2003 December 15–19; Norfolk, Virginian, USA. Reston: AIAA; 2003 [Report no AIAA-2003-7035]. 8. Baurle R, Hsu A, Hassan H. Assumed and evolution probability density functions in supersonic turbulent combustion calculations. J Propul Power 1995;11(6):1132–8. 9. Peters N. Laminar diffusion flamelet models in non-premixed turbulent combustion. Prog Energy Combust Sci 1984;10 (3):319–39. 10. Peters N. Turbulent combustion. Cambridge: Cambridge University Press; 2000. 11. Peters N. Laminar flamelet concepts in turbulent combustion. Twenty-first symposium (international) on combustion, vol. 21(1). p. 1231–50. 12. Pierce CD, Moin P. Progress-variable approach for large-eddy simulation of non-premixed turbulent combustion. J Fluid Mech 2004;504:73–9. 13. Ihme M, See YC. Prediction of autoignition in a lifted methane/air flame using an unsteady flamelet/progress variable model. Combust Flame 2010;157(10):1850–62. 14. Saghafian A, Shunn L, Philips DA, Ham F. Large eddy simulations of the HIFiRE scramjet using a compressible flamelet/ progress variable approach. Proc Combust Inst 2015;35 (2):2163–72. 15. Saghafian A, Terrapon VE, Pitsch H. An efficient flamelet-based combustion model for compressible flows. Combust Flame 2015;162(3):652–67. 16. Bray KNC, Libby PA, Williams FA. High speed turbulent combustion. In: Libby PA, Williams FA, editors. Turbulent reacting flow. New York: Academic Press; 1994.

1043

Please cite this article in press as: Chen B et al. Numerical simulations of turbulent flows in aeroramp injector/gas-pilot flame scramjet, Chin J Aeronaut (2017), http:// dx.doi.org/10.1016/j.cja.2017.06.010

1031 1032 1033 1034 1035 1036

1039 1040 1041

1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090

CJA 873 21 June 2017

18 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130

17. Waidmann W, Alff F, Bo¨hm M, Brummund U, Claub W, Oschwald M. Supersonic combustion of hydrogen/air in a scramjet combustion chamber. Space Technol 1996;15(6):421–9. 18. Williams FA. Progress in knowledge of flamelet structure and extinction. Prog Energy Combust Sci 2000;26(4–6):657–82. 19. Fan ZQ, Liu WD, Sun MB, Wang ZG, Luo WL. Theoretical analysis of flamelet model for supersonic turbulent combustion. Sci China: Technol Sci 2012;55(1):193–205. 20. Wei BX, Xu X, Yan ML, Shi XX, Yang Y. Study on aeroramp injector/gas-pilot flame in a supersonic combustor. J Propul Power 2012;28(3):486–95. 21. Wei BX, Chen B, Yan ML, Shi XX, Zhang Y, Xu X. Operational sensitivities of an integrated aerodynamic-ramp-injector/gas portfire flame holder in a supersonic combustor. Acta Astronaut 2012;81(1):102–10. 22. Zhang Y, Chen B, Liu G, Wei BX, Xu X. Influencing factors on the mode transition in a dual-mode scramjet. Acta Astronaut 2014;103:1–15. 23. UCSD. A short chemical-kinetic mechanism for C2H4 ignition and detonation, mechanical and aerospace engineering. [Internet]. [cited 2014 Jan 1]. 24. Pitsch H, Seiser R, Varatharajan B, A guide to FlameMaster 3.3. [Internet]. [cited 2013 June 13]. 25. Lien FS, Liu H, Chui E, McCartney CJ. Development of an analytical b-function PDF integration algorithm of simulation of non-premixed turbulent combustion. Flow Turb Combust 2009;83 (2):205–26. 26. Wilcox DC. Turbulence modeling for CFD. California: DCW Industries Inc; 2006. 27. Huang W, Liu WD, Li SB, Xia ZX, Liu J, Wang ZG. Influences of the turbulence model and the slot width on the transverse slot injection flow field in supersonic flows. Acta Astronaut 2012;73:1–9. 28. Huang W, Yan L. Progress in research on mixing techniques for transverse injection flow fields in supersonic crossflows. J Zhejiang Univ – Sci A 2013;14(8):554–64. 29. Menter FR. Two-equation eddy-viscosity turbulence models for engineering applications. AIAA J 1994;32(8):1598–605.

No. of Pages 18

B. Chen et al. 30. Kim SD, Song DJ. Modified shear-stress transport turbulence model for supersonic flows. J Aircraft 2005;42(5):1118–25. 31. Rumsey CL. Compressibility considerations for k-x turbulence models in hypersonic boundary-layer applications. J Spacecraft Rock 2010;47(1):11–20. 32. Wilcox DC. Formulation of the k-x turbulence model revisited. AIAA J 2008;46(11):2823–38. 33. Shuen JS. Upwind differencing and LU factorization for chemical non-equilibrium Navier-Stokes equations. J Comput Phys 1992;99 (2):233–50. 34. Toro EF. Riemann solvers and numerical methods for fluid dynamics. Berlin: Springer Press; 2009. 35. Kitamura K, Shima E. Simple and parameter-free second slope limiter for unstructured grid aerodynamic simulations. AIAA J 2012;50(6):1415–26. 36. Mavriplis DJ. Unstructured-mesh discretizations and solvers for computational aerodynamics. AIAA J 2008;46(6):1281–98. 37. Venkatakrishnan V. Convergence to steady state solutions of the Euler equations on unstructured grids with limiters. J Comput Phys 1995;118(1):120–30. 38. George K, Kirk S, Vipin K. PARMETIS: parallel graph partitioning and sparse matrix ordering library, Version 3.1. Twin Cities: University of Minnesota; 2003. 39. Roache PJ. Verification and validation in computational science and engineering. Albuquerque: Hermosa Pub; 1998. 40. Yamashita H, Shimada M, Takeno T. A numerical study on flame stability at the transition point of jet diffusion flames. Twenty-sixth symposium (international) on combustion, vol. 26(1). p. 27–34. 41. Kouchi T, Masuya G, Mitani T, Tomioka S. Mechanism and control of combustion-mode transition in a scramjet engine. J Propul Power 2012;28(1):106–12. 42. Bonanos AM, Schetz JA, O’Brien WF, Goyne CP. Dual-mode combustion experiments with an integrated aeroramp-injector/plasma-torch igniter. J Propul Power 2008;24(2):267–73. 43. Cabell K, Hass N, Storch A, Gruber M. HIFiRE Direct-Connect Rig (HDCR) phase I scramjet test results from the NASA Langley arc-heated scramjet test facility. Reston: AIAA; 2011 [Report no AIAA-2011-2248].

Please cite this article in press as: Chen B et al. Numerical simulations of turbulent flows in aeroramp injector/gas-pilot flame scramjet, Chin J Aeronaut (2017), http:// dx.doi.org/10.1016/j.cja.2017.06.010

1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164 1165 1166 1167 1168 1169