- Email: [email protected]
Z flamelet model
Acta Astronautica 127 (2016) 375–383
Contents lists available at ScienceDirect
Acta Astronautica journal homepage: www.elsevier.com/locate/aa
The hybrid RANS/LES of partially premixed supersonic combustion using G/Z flamelet model Jinshui Wu a, Zhenguo Wang a,n, Xuesong Bai b, Mingbo Sun a,n, Hongbo Wang a a b
Science and Technology on Scramjet Laboratory, National University of Defense Technology, Changsha, 410073, China Division of Fluid Mechanics, Lund University, P. O. Box 118, S-221 00 Lund, Sweden
art ic l e i nf o
a b s t r a c t
Article history: Received 25 February 2016 Received in revised form 1 June 2016 Accepted 10 June 2016 Available online 11 June 2016
In order to describe partially premixed supersonic combustion numerically, G/Z flamelet model is developed and compared with finite rate model in hybrid RANS/LES simulation to study the strut-injection supersonic combustion flow field designed by the German Aerospace Center. A new temperature calculation method based on time-splitting method of total energy is introduced in G/Z flamelet model. Simulation results show that temperature predictions in partially premixed zone by G/Z flamelet model are more consistent with experiment than finite rate model. It is worth mentioning that low temperature reaction zone behind the strut is well reproduced. Other quantities such as average velocity and average velocity fluctuation obtained by developed G/Z flamelet model are also in good agreement with experiment. Besides, simulation results by G/Z flamelet also reveal the mechanism of partially premixed supersonic combustion by the analyses of the interaction between turbulent burning velocity and flow field. & 2016 IAA. Published by Elsevier Ltd. All rights reserved.
Keywords: Supersonic Combustion Partially premixed G/Z flamelet model
1. Introduction In scramjet engines, in order to achieve flameholding, different kinds of injection methods and flameholding devices are introduced. The fully mixing of fuel and air will take a long distance under supersonic condition, however, during the mixing process combustion will take place. This kind of combustion is partially premixed. As a combination of the premixed and the non-premixed mode of combustion, partially premixed combustion shows features of both. One of the most popular flame stabilization mechanisms in partially premixed combustion under low speed condition, i.e. the triple flame, considers that there exist a fuel lean premixed branch, a fuel rich premixed branch and a diffusion branch. In order to describe to this kind of combustion under low speed condition, many work have been done. Müller et al. [1] have developed a model for partial premixed turbulent combustion that they validate with lift-off heights in turbulent jet diffusion flames. Domingo et al. [2] have proposed a flame index to distinguish combustion mode in the entire flow field. Pierce et al. [3] developed a progress variable approach to account for flame quenching in non-premixed combustion. Pitsch and Ihme [4] have proposed an unsteady/ flamelet progress variable method in non-premixed turbulent combustion. However, all these methods are developed under low speed condition. For high speed condition, especially n
Corresponding authors. E-mail addresses: [email protected] (Z. Wang), [email protected] (M. Sun). http://dx.doi.org/10.1016/j.actaastro.2016.06.021 0094-5765/& 2016 IAA. Published by Elsevier Ltd. All rights reserved.
supersonic flow, whether these methods are valid is still unknown. There are two main issues considering the numerical simulation and combustion modeling of partially premixed supersonic combustion. Firstly, the combustion model chosen to describe the chemical reaction should be carefully considered. In partially premixed zones under supersonic condition, the range of the temperature distribution changes quite a lot. The maximum temperature sometimes could even reach 3000 K while the minimum temperature may be in the same magnitude with the temperature of cold flow. Under these conditions, the applicability of the traditional finite rate chemical reaction model is restricted. It may fail to effectively predict chemical reaction with broad temperature changes in partially premixed zones under supersonic condition by using finite rate chemical reaction model. Secondly, the flame propagation must be taken into account when there exists the partially premixed zones. Under supersonic condition, the flame propagation speed is much larger than low speed flow condition. The thickness of the premixed flame front could be thinner than a grid size of the LES. The grid resolution of LES is not capable of capturing the flame front without modeling. In order to investigate partially premixed supersonic combustion in scramjet engines, the framework of this essay is based on the model proposed by Müller et al. [1] which combines a level set method of describing premixed flame front and flamelet model. In this paper we have also drawn on many previous experiences of the LES of supersonic flow and combustion in our group [5–11]. A strutinjection supersonic combustion which has been experimentally investigated in German Aerospace Center was numerically simulated
376
J. Wu et al. / Acta Astronautica 127 (2016) 375–383
by G/Z flamelet model and finite rate chemistry model. In finite rate chemistry model, the reacting species transport equations together with their chemical source terms are solved, for simplicity the description of species transport equations are omitted and no turbulent combustion modeling is considered this essay. By our previous study [12], the maximum value of the Karlovitz number in main regions of this combustor is smaller than 1.0, indicating that the flame scales are smaller than the Kolmogorov scales so that the assumption of the flamelet model is satisfied under the given condition. Many other authors have also investigated this combustor numerically by flamelet model or finite rate chemistry model [13–17], however, they did not focus on the partially premixed supersonic combustion. Some simulation results could not well reproduce the characteristic of the partially premixed combustion under supersonic condition. In this essay, the practicability and accuracy of applying G/Z flamelet model into partially premixed supersonic combustion will be analyzed in the following part.
According to the discussion of the Kim et al. [21], modeling constant α ¼ 0.2 and β ¼ 6.67. u′ is local velocity fluctuation, in this paper the ksgs distribution is used to provide an estimate of the local velocity fluctuation u′ = 2k sgs/3 . In this essay, the SGS model chosen is the one-equation model of Yoshizawa and Horiuti [22]. Menter BSL (Baseline) model [23] is considered for the RANS portion. A similar blending function formulation [7] based on the hybrid RANS/LES methodology proposed by Baurle et al. [24] is used in this work. For brevity the governing equations for SGS model and RANS model and hybrid RANS/LES formula are omitted. 2.2. Re-initialization and a modified sub-cell-fix approach
∼ To obtain G equation in the whole field the so-called re-in∼ itialization equation ∇G = 1 is always used. A re-distancing technique introduced by Sussman et al. [25] is used to iterate on ∼ ∼ ∼ the equation Gτ = sign (G ) 1 − ∇G until a steady solution is
(
2. Model descriptions and numerical methods 2.1. Governing equations Flamelet model for premixed combustion in this paper is based on the level set approach using the scalar G. The distance function G(x, t) represents the distance of a certain point to the flame front. Although G equation model could well capture the interaction between large eddy structures and the premixed flame fronts, however, it could not be used to describe flame differences due to fuel concentration. Another scalar equation which considers diffusion effect, i.e. mixture fraction equation Z, is combined with G equation to account for partially premixed combustion zones. The Favre filtered compressible Navier-Stokes equations and flamelet equations [18] can be written as below:
∂ρ¯ ∂(ρ¯ u˜ i ) + =0 ∂t ∂xi ∂ ⎡⎣ ρ¯ u˜ i u˜ j + p¯ δij − τ¯ij − τijHybrid ⎤⎦ ∂(ρ¯ u˜ i ) + =0 ∂t ∂xj ∂ ⎡⎣ ρ¯ E˜ + p¯ u˜ i + q¯i − u˜ j τ¯ji + HiHybrid + σiHybrid ⎤⎦ ∂ρ¯ E˜ =0 + ∂xi ∂t Hybrid ∂ρ¯ Z˜ ∂(ρ¯ u˜ i Z˜ + Z˜ ) ∂ ⎛ ¯ ∂Z˜ ⎞ + = ⎜ ρ¯ D ⎟ ∂t ∂xi ∂xi ⎝ ∂xi ⎠ ∼ ∼ ∼ ∼ ∼ ∂ρ¯ G ∂(ρ¯ u˜ i G ) + = ρ¯ st ∇G − ρ¯ DtG κ c ∇G ∂t ∂xi
(
)
(1)
where is ρ¯ the filtered density, u˜ i is the filtered velocity, p¯ is the filtered pressure, E˜ is the filtered total energy, D Is the laminar thermal diffusion coefficient, τ¯ij is the filtered laminar viscous stress tensor, τijHybrid is the unclosed stress tensor, q¯i is the energy flux caused by heat conduction and mass diffusion, HiHybrid is the unclosed enthalpy flux vector, σiHybrid is the unclosed heat flux vector, is the unclosed mixture fraction flux vector which could be ∼ ∼ Hybrid expressed as Z˜ ≈ − ρ¯ DtZ ∇Z using gradient diffusion assump∼ ∼ tion and κc is filtered flame curvature ( κc = ∇⋅n = ∇⋅ −∇G / ∇G , n
(
)
is the normal direction vector). st is the turbulent burning velocity. The other parameters and unclosed items mentioned above can be found in reference [19]. In this essay, the turbulent burning velocity model proposed by Pocheau [20] is applied.
2.3. Partially premixed combustion modeling for supersonic flow Particularly, the iso-surface of G ¼ 0 has defined the position of the instantaneous premixed flame front, and divided the whole flow field into two region, where G 4 0 stands for the burnt region and G o 0 stands for the unburnt region. However, it has to be pointed out that for partially premixed flame, G equation has no meaning in regions exceeding lean or rich limit. However, for simplicity, G(x, t) ¼ 0 surfaces are still calculated in whole flow field, only in G 4 0 regions the flamelet libraries are applied. The flamelet model adopts a series of flamelet parameters to represent the chemical thermodynamic state, for example, the scalar dissipation rate is used to characterize the non-equilibrium effects, and the joint probability density function to characterize the interaction between the turbulence and the combustion. Considering the shock wave capturing algorithm in the supersonic flow, the temperature is calculated from the energy equation. The procedure is described as follows. The total energy could be written as NS
E=
∑ ρi (hi0 + hiT ) − p + i=1
1 ρ (u2 + v2) 2
(2)
(3)
hi0
where is non-standard enthalpy of formation which is independent of temperature, hiT is heat enthalpy which relates with temperature. A new total energy format is E′ introduced by neglecting, hi0 , NS
E′ =
∑ ρi hiT i=1
−p+
1 ρ (u2 + v2) 2 NS − ∑i = 1 ωi hi0
1
⎛ u′ ⎞α ⎤ α st ⎡ = ⎢1 + β⎜ ⎟ ⎥ ⎝ sL ⎠ ⎦ sL ⎣
)
reached. The converged steady solution gives the signed distance function. Efforts to improve the accuracy of distance function in level-set methods have led to a variety of re-initialization methods different from the original re-distancing algorithm. In this essay, a modified sub-cell fix method proposed by Sun et al. [26] based on Russo [27] and Hartmann [28] is applied. This modified sub-cellfix scheme independent of local curvature makes use of a combination of the points adjacent to zero level-set surfaces and preserves the interface in a second-order accuracy. The scheme is capable of handling large local curvature, and as a result it demonstrates satisfactory performance on several challenging test cases [26]. For brevity, the modified sub-cell-fix process will not be repeated here.
(4)
Then, a source term will appear in the new transport equation for E′. A time-splitting method is used to uncouple the flow and reaction as below.
J. Wu et al. / Acta Astronautica 127 (2016) 375–383
⎛ Δt ⎞ ⎛ Δt ⎞ Qn + 1 = L c ⎜ ⎟ L f (Δt ) L c ⎜ ⎟ Qn ⎝ 2⎠ ⎝ 2⎠
(5)
where Lf and L c denote the solver for flow and reaction, respecNS k tively. The filtered source term − ∑i = 1 ωihi0 appears in the reaction solver. During the reaction step, the new total energy will be updated due to reaction heat release. And the temperature will be calculated implicitly by newton iteration. The mass fractions of the species are interpolated from the flamelet libraries. The filtered ∼ mass fraction of the ith species Yi in the flamelet libraries is expressed by an ensemble of the laminar diffusion flamelets:
∼∼ Yi (Z , k Z″ 2, χ˜st ) =
˜ χ˜ ∬ Yi (Z˜ , χ˜st ) P (Z˜ , χ˜st ) dZd st
Cχ μ C μt 2 ∇Z˜ + 2Z t2 ε k Z″ 2 Δ Sct Cμ Sct
C χ Z = 2.0
Cε = 0.5 Cμ
(7)
k Z ″2 is calculated by scale similarity model [30]:
k Z ″2 = CZ Δ2 ∇Z˜ ⋅∇Z˜
Fig. 1. Sketch of the DLR scramjet combustor. Table 1 Flow conditions of the airstream and the hydrogen jet.
(6)
where k Z″ 2 is the filtered mixture fraction variance, χ˜st is the fil˜ χ˜ ) is the species mass fraction of tered scalar dissipation rate, Yi (Z, st ˜ χ˜ ) is the joint probability the laminar flamelet libraries and P (Z, st ∼ ∼ ˜ χ˜ ), Z and χ˜st density function of Z and χ˜st . When calculating P (Z, st are assumed to be statistically independent of each other, and they are subject to the Beta distribution and the Dirac function, respectively. In hybrid LES/RANS simulation, a direct calculation of scalar dissipation rate is not possible. Instead a filtered value must be calculated which contains both resolved and SGS-modeled parts [29]:
χ˜ = 2
377
(8)
where CZ is the model parameter and equals 0.2. Δ is the filtered scale. The flamelet libraries are generated by FlameMaster V3.9 software [31] using the mechanism containing 9 species and 19 step reactions [32]. Lewis numbers of these species are assumed to be unity. Radiation heat loss and pressure changes with time are not considered in the flamelet libraries. According to former study [33], the pressure boundary condition for generating flamelet libraries is 1.5 bar and the temperature of the oxidizer and fuel are 1000 K and 250 K, respectively. The size of the flamelet libraries is ∼ 1000*100*40 in Z , k Z″ 2 and χ˜ space. Besides, for some quenching regions where burning velocity is zero, when recirculation zone exists, G˜ (x, t) ¼ 0 surfaces will move with flow recirculation to form uncertain flame front. Therefore, ∼ on one side, G˜ equation should combine with Z to determine flame front; on the other side, a method is introduced into G˜ ∼ ∼ ∂ρ¯ G equation in quenching regions [34], where + ∇ρ¯ u˜ *G = 0, ∂t ˜1 , u˜2 , u˜3 ), u˜1 is the velocity point to the flow u˜ * = ( u˜1*, u˜2*, u˜* 3) = ( u recirculation direction. 2.4. Computational model and numerical method The sketch of the DLR scramjet combustor [35,36] is given in Fig. 1. In the numerical simulations of Ref. [35], in order to consider the shock capture scheme in high speed flow with the use of flamelet model, only the species mass fractions are interpolated from the libraries, the local temperature is implicitly calculated from the solution of the energy using the caloric equation of state. This method is also used in the modeling of partially premixed supersonic combustion in above section. In this combustor, vitiated air at a nominal temperature of 340 K flows over a strut. Hydrogen at sonic conditions exits through fifteen holes uniformly spaced
Air
H2
Ma T (K) P (kPa) YO2
2.0 340 100 0.232
1.0 250 100 0
YN2
0.736
0
YH2 O
0.032
0
YH2
0
1.0
along the back plane of the fuel injector, where it mixes with freestream air. The combustor was designed and experimentally investigated by German Aerospace Center. Experimental pressure distributions, axial velocity profiles, velocity fluctuation profiles, and static temperature profiles are available for both non-reacting (mixing) and reacting (hydrogen oxidation) cases. Detail descriptions of the combustor are omitted here. The original 15 holes are substituted by three holes to save the computational resources, and the computational domain has a length of 22 mm in the spanwise direction and periodic boundary condition. The bottom wall of the combustor has the coordinate of y ¼0 mm and the tip of the strut has the coordinates x ¼ 35 mm and y ¼ 25 mm. The structured mesh has approximately 6.7 million cells according to research earlier [13,14]. The mesh at the vicinity of the inlet, the strut, the wake flow region and the shear layer region is refined. The whole configuration is divided into 24 sections which are calculated in parallel. The flow conditions of the incoming airstream and the hydrogen jet are given in Table 1. In supersonic flow, the fluctuation in the downstream cannot propagate upstream. Therefore, the inflow boundary condition the boundary conditions for the hydrogen jet are fixed with the given parameters. The outflow is assumed to be supersonic such that its boundary condition is treated using extrapolation from the interior of the domain. According to the discussion of Génin et al. [14], the top and bottom walls of the chamber are assumed to be adiabatic and slip wall so as to avoid the interaction between the shock wave and the boundary layer. The surfaces of the strut are treated as no-slip and adiabatic walls. The fifth-order WENO (weighed essentially non-oscillation) scheme is used for the inviscid fluxes in order to improve the resolution of the flow field. The time integration is performed by means of a second-order accurate total-variation-diminishing (TVD) Runge-Kutta method. The Courant-Friedrichs-Lwey (CFL) number is fixed to be 0.5.
3. Results and discussions In our group, we have obtained many results about this combustor without and with combustion, and compared these results with experiment to test the validity and efficiency of our programming code for supersonic combustion [12]. In this essay, we
378
J. Wu et al. / Acta Astronautica 127 (2016) 375–383
Fig. 2. Perspective view in terms of contours of the axial momentum on the back plane ∼ ∼ ∼ and, an iso-surface of G = 0 colored by temperature and an iso-surface of Z = Zst .
only cared about the partially premixed combustion model we used and partially combustion phenomenon in this combustor, results without combustion will no longer be presented. Fig. 2 presents the perspective view of contours of the axial ∼ momentum, an iso-surface of G = 0 colored by temperature and an ∼ ∼ iso-surface of Z = Zst . It could be seen from Fig. 2 that the flow ∼ field has a great impact on the evolution of G equation, the interaction between jet and airstream in shear layer has wrinkled ∼ the iso-surface of G = 0 . The heat release in the recirculation zone behind the strut has increased the local pressure and expanded ∼ the recirculation zone, which also leads to the expansion of the G surface towards the airstream. Fig. 2 could also verify that the reinitialization and modified sub-cell-fix approach mentioned above ∼ for G equation are capable of capturing the surfaces of premixed front under complex computational model and complex flow condition. From the contour of the axial momentum on the back plane, the shock waves in the entire flow are well captured from ∼ ∼ calculation. It is shown in Fig. 3(a) that iso-surface of Z = Zst only covers a small region of the entire flow field, however, the iso∼ surface of G = 0 showed in Fig. 2 could extend to the entire region, most of the reactions occur on fuel lean condition. It could be seen ∼ in the Fig. 3(b) that the iso-surface of YOH = 0.1 indicates that there ∼ ∼ exists strong reaction zone in the vinicity of iso-surface of Z = Zst showed in Fig. 3(a). Besides, from the contours of the temperature on six cross planes showed in Fig. 3(c), we could see that the temperature in the second slice reaches maximum corresponding to the most intensive reaction zone compared with the other slices. The maximum temperature value is about 2500 K evaluated from Fig. 3(c). Fig. 4 further visualizes the instantaneous contour of tem∼ ∼ ∼ perature distribution together with Z = Zst line and G = 0 line at different spanwise locations. All high temperature zones at dif∼ ferent locations are surrounded G = 0 line. In the vicinity of the ∼ injector orifice, G = 0 lines drawn in the enlarged graph of Figs. 4 and 5 have engulfed the jet fluid and could extend to the thin mixing layer between the jet and the airstream. In the thin mixing layer after the strut, the mixture reaches the stoichiometric condition so that the turbulent velocity is large enough to overcome the flow speed in the shear layer and pushes the iso-surface of ∼ G = 0 toward to the strut edge. Due to the existence of fuel jet, iso∼ ∼ surface of G = 0 around the jet is broken into small G = 0 islands which will be transported to downstream as non-reacting circles as shown in both Figs. 4 and 5. As hydrogen jet develops downstream, mixing of the fuel and the air is promoted due to the largeeddy structures, mixture comes close to the stoichiometric state ∼ ∼ ( Z = Zst ) and mixing layers are broadened and the reaction zone moves from the upper and lower layers shown in Fig. 5 to the center of the flow field. It can be seen that the temperature in the center of the flow field achieves its peak value at the vicinity of
Fig. 3. Iso-surfaces and perspective view of computational results.(a) Iso-surface of ∼ ∼ ∼ Z = Zst (b) Iso-surface of YOH = 0.1 (c) Perspective view in terms of temperature on the six cross planes.
∼ ∼ line Z = Zst . In further downstream, combustion intensity decreases, reaction zone begins to shrink and breaks into small pieces and the temperature in the flow field decreases gradually due to continuous combustion consumption of fuel in the upstream. It also could be seen in the global graphs at different locations in Fig. 4 that the reaction wake flows are quite not the ∼ same due to the unsteadiness of the G = 0 surface, at Z/ZD ¼ 0.25, the reaction zone is close to the upper wall while at Z/ZD ¼ 0.75 the reaction zone is close to the lower wall. ∼ The turbulent burning velocity affected by Z is calculated using ∼ Eq. (2) and it combines with velocity to determine the evolution G and affects the combustion process indirectly. Fig. 6 gives the instantaneous contours of turbulent burning velocity and Mach number at different spanwise locations. It could be seen in Fig. 6 that the subsonic region downstream the strut mainly locates in the fuel rich region. Compared to white lines shown in Fig. 5, in the subsonic region, fuel is injected continuously to sustain the combustion and subsonic local speed condition is in favor of flame stabilization. It is illustrated in Fig. 6 that turbulent burning velocity differs from each other at different locations. Summarized from Fig. 6, the schematic of supersonic partially premixed combustion is also depicted in Fig. 7. We could see that turbulent burning velocity st1 (in shear layer near strut) are generally larger
J. Wu et al. / Acta Astronautica 127 (2016) 375–383
379
Fig. 4. Calculated front view of spanwise slice at Z/ZD ¼ 0.25, 0.5 and 0.75 (from up to down, ZD ¼22mm), instantaneous contour of temperature distribution together with ∼ ∼ ∼ Z = Zst line (white) and G = 0 line (black).
∼ Fig. 5. Snapshots of temperature contour evolvement together with G = 0 line (black) on central plane behind the strut at sixe different times ( Δt = 1.769 × 10−5s ).
than st2 (in subsonic zone) and st3 (close to the main flow). Under the interaction of turbulent burning velocity st1 and the recirculation zones generated near jet injection in the shear layer ∼ close to the strut, the iso-surface of G = 0 can propagate to the region very close to the strut edge. In subsonic zone, turbulent burning velocity st2 is small due to fuel rich condition, and the fuel ∼ injection velocity will push the iso-surface of G = 0 away from the strut. Together with the interaction of the recirculation zones ∼ generated by injection, the iso-surface of G = 0 is broken into very small iso-surfacs as shown in Figs. 4–6. In the region closes to the main flow, the turbulent burning velocity st3 is in competition with ∼ local velocity's projection on the iso-surface of G = 0 . Hot reaction zone almost exists in the region fuel/air mixture reaches
stoichiometry condition. Further downstream, the turbulent burning velocity decreases due to the consumption of the fuel and the large eddy structures generated have divided the iso-surface of ∼ G = 0 into small iso-surface islands. In these further regions, the flow field characteristics plays a more important role on the evo∼ lution of G equation. Fig. 8 presents time average temperature profiles using two different models compared with experiment at different locations. At x ¼78 mm, the peak value of temperature in the shear layer predicted by both models are higher than experimental data. This is also obtained by Berglund using two different flamelet models [13] and by Fureby et al. using different PaSR models [16]. There are two main factors that account for the over prediction of the
380
J. Wu et al. / Acta Astronautica 127 (2016) 375–383
∼ ∼ ∼ Fig. 6. Calculated instantaneous front view of spanwise slices at Z/ZD ¼ 0.25, 0.5 and 0.75 (from up to down, ZD ¼ 22 mm) together with Z = Zst line (black) and G = 0 line (purple).
temperature peak in the shear layer of the G/Z model. First, in the shear layer, the local equivalent ratio is close to stoichiometry equivalent ratio, so that the turbulent burning velocity is large ∼ enough to push the propagation of G = 0 surface forward to the edge of the strut as shown in Fig. 4. Then a more fierce combustion than experiment in the shear layer will be predicted. Second,
according to flamelet theory, if the local scalar dissipation rate is large enough, the flamelet solution will move close to the extinction. It is reasonable that the filtered scalar dissipation rate is under estimated so that the temperature peak obtained by flamelet model is larger than experiment. However, at this same location, the temperature profiles in the center zone after the strut predicted by the G/Z model are in better agreement with
J. Wu et al. / Acta Astronautica 127 (2016) 375–383
experiment than finite rate model. The finite rate model we used has over predicted the temperature which is similar with the results of Dharavath et al. [17]. In this zone, the combustion is partially premixed, the application of partially premixed model in this article could well reproduce this low temperature reaction zone by considering the interaction between G equation and Z equation. The simulation results obtained Berglund using two different
Fig. 7. Schematic of partially premixed supersonic combustion.
381
flamelet models have under predicted the temperature at this location. Further investigations by other authors show no better improvement [14–16], they have also under estimated the temperature in this partially premixed zone. At x¼ 125 mm, the width of the reaction zone predicted by G/Z model is better than finite rate model compared with experimental data. The results by finite rate model have broaden the reaction zone, this is also obtained by the work of Génin and Menon [14] and Dharavath et al. [17]. The peak value of the temperature predicted by both models shown in Fig. 8 are in good agreement with experimental data. At x¼ 233 mm, the width of the reaction zone and the peak value of the temperature predicted are smaller and lower than the experimental data compared with finite rate model. This will be explained in the following analysis. In order to analyze the influence of flow-field on the devel∼ opment of G and combustion, the average velocity and velocity fluctuation obtained by G/Z flamelet model and finite rate model are compared with experiment shown in Figs. 9–11. The average
Fig. 8. Average temperature compared with experiment using different combustion models at different axial locations.
Fig. 9. Average axial velocity compared with experiment using different combustion models at different axial locations.
Fig. 10. Average axial velocity fluctuation compared with experiment using different combustion models at different axial locations.
382
J. Wu et al. / Acta Astronautica 127 (2016) 375–383
Fig. 11. Average normal velocity fluctuation compared with experiment using different combustion models at different axial locations.
axial velocity predictions shown in Fig. 9 at x ¼ 78 mm in the shear layer are much smaller than the experiment by both models, even negative velocity values which could cause the generation of the recirculation are predicted. According to the ∼ ∼ properties of G equation, the surface of G = 0 could reach further upstream which causes the inaccurate prediction of reaction. This is consistent with the static temperature contours shown in Fig. 4. The average velocity fluctuation predictions in the shear layer are slightly over estimated by both models as shown in Figs. 10 and 11. The turbulent buring velocity are proportional to the velocity fluctuation. The over prediction of the velocity fluctuation could also contribute to the propagation of the surface of ∼ G = 0. At x ¼ 125 mm in Figs. 9–11, it is shown that deviations of average velocity and average velocity fluctuation compared with experiment by G/Z flamelet model are in more acceptable level than finite rate model so that G/Z flamelet model could well reproduce the overall combustion properties. Further downstream at x ¼233 mm, the average velocity fluctuation predicted by finite rate model is larger G/Z flamelet model, which may contribute to more intense reaction and higher temperature. The peak value of temperature and the width of the reaction zone predicted by G/Z flamelet model shown in Fig. 8 are smaller and lower than the experimental data. Besides the underestimation of velocity fluc∼ tuation, this may also be explained by the inaccuracy of G ∼ equation in the downstream for the surfaces of G = 0 are broken into pieces shown in Fig. 4 which could not sustain continous combustion.
4. Conclusions The partially premixed supersonic combustion in the DLR scramjet combustor is investigated numerically by hybrid RAN/LES using two different combustion models. Simulation results show that the combination of G equation method and flamelet libraries could capture the characteristics of partially premixed supersonic combustion. Further analyses of the comparsion between two models and experimental data show that the G/Z flamelet model used in this combustor could reproduce the low temperature reaction zone in the region after the strut. In this region, the predictions of temperature profiles by G/Z model are in better agreement with experiment than finite rate model and other simulation results, while the temperature in the shear layer are over predicted by the cause that filtered scalar dissipation rate might be under estimated. Although there exist some differences compared with experment in predicting the width of the reaction zone and temperature at further downstream, the deviations of average velocity and average velocity fluctuations compared with experiment are in acceptable level so that G/Z flamelet model could reproduce the overall combustion properties.
This paper validates the rationality and feasibility of the G/Z flamelet model in partially premixed supersonic combustion. On the one hand, it implies that the partially premixed combustion should be carefully dealt with. In order to reproduce the overall properties of the partially premixed combustion under supersonic condition, the modeling for numerical simulation is necessary. The research provides a guidance for the application of the G/Z flamelet model in partially premixed supersonic combustion. On the other hand, it points out there are some problems when applying the G/Z flamelet model in the shear layer near the strut and further downstream under supersonic condition. Some other modeling aspects should be taken into account for the quenching effect near the strut and elimating the strong dependence on G equation in further downstream.
Acknowledgments The research work was funded by the National Natural Science Foundation of China (Nos. 11472305, 11522222 and 91541101).
References [1] C.M. Muller, H. Breitbach, N. Peters, Partially premixed turbulent flame propagation in jet flames, Symposium (International) on Combustion, 25 (1994) 1099-1106. [2] P. Domingo, L. Vervisch, K. Bray, Partially premixed flamelets in LES of nonpremixed turbulent combustion, Combust. Theory Model. 6 (2002) 529–551. [3] C.D. Pierce, P. Moin, Progress-variable approach for large eddy simulation of non-premixed turbulent combustion, J. Fluid Mech. 504 (2004) 73–97. [4] H. Pitsch, M. Ihme, An unsteady/flamelet progress variable method for LES of nonpremixed turbulent combustion, in: Proceedings of the AIAA Paper 20052557, 2005. [5] H.B. Wang, Z.G. Wang, M.B. Sun, N. Qing, Combustion characteristics in a supersonic combustor with hydrogen injection upstream of cavity flameholder, Proc. Combust. Inst. (2013) 2073–2082. [6] M.B. Sun, H. Geng, J.H. Liang, Z.G. Wang, Mixing characteristics in a supersonic combustor with gaseous fuel injection upstream of a cavity flameholder, Flow., Turbul. Combust. 82 (2009) 271–286. [7] M.B. Sun, Z.G. Wang, J.H. Liang, H. Geng, Flame characteristics in a supersonic combustor with hydrogen injection upstream of a cavity flameholder, J. Propuls. Power 24 (2008) 688–696. [8] H.B. Wang, Z.G. Wang, M.B. Sun, N. Qin, Large eddy simulation of a hydrogenfueled scramjet combustor with dual cavity, Acta Astronaut. 108 (2012) 119–128. [9] H.B. Wang, Z.G. Wang, M.B. Sun, N. Qin, Large eddy simulation based studies of jet-cavity interactions in a supersonic flow, Acta Astronaut. 93 (2014) 182–192. [10] C.Y. Liu, Z.G. Wang, H.B. Wang, M.B. Sun, Numerical investigation on mixing and combustion of transverse hydrogen jet in a high-enthalpy supersonic crossflow, Acta Astronaut. 116 (2015) 93–105. [11] Y.X. Yang, Z.G. Wang, M.B. Sun, H.B. Wang, L. Li, Numerical and experimental study on flame structure characteristics in a supersonic combustor with dualcavity, Acta Astronaut. 117 (2015) 376–389. [12] Z.Q. Fan, W.D. Liu, M.B. Sun, Z.G. Wang, F.C. Zhuang, W.L. Luo, Theoretical analysis of flamelet model for supersonic turbulent combustion, Sci. China Tech. Sci. 55 (2012) 193–205. [13] M. Berglund, C. Fureby, LES of supersonic combustion in a scramjet engine
J. Wu et al. / Acta Astronautica 127 (2016) 375–383
model, Proc. Combust. Inst. 31 (2007) 2497–2504. [14] F. Génin, S. Menon, Simulation of turbulent mixing behind a strut injector in supersonic flow, AIAA J. 48 (2010) 526–539. [15] A.S. Potturi, J.R. Edwards, LES/RANS simulation of a supersonic combustion experiment, in: 50th AIAA Aerospace Sciences Meeting, AIAA Paper 2012-0611, 2012. [16] C. Fureby, E. Fedina, J. Tegner, A computational study of supersonic combustion behind a wedge-shaped flameholder, Shock. Waves 24 (2014) 41–50. [17] M. Dharavth, P. Manna, D. Chakraborty, Thermochemical exploration of hydrogen combustion in generic scramjet combustor, Aerosp. Sci. Technol. 24 (2013) 264–274. [18] N. Peters, Turbulent Combustion, Cambridge University Press, 2000. [19] M.B. Sun, Flow and flameholding mechanisms of cavity flameholders in supersonic flows (PhD thesis), National University of Defense Technology, 2008. [20] A. Pocheau, Scale invariance in turbulent front propagation, Phys. Rev. E 49 (1994) 1109–1122. [21] W.W. Kim, S. Menon, H. Mongia, Large eddy simulations of a gas turbine combustor flow, Combust. Sci. Technol. 143 (1999) 25–62. [22] A. Yoshizawa, K. Horiuti, A statistically-derived subgrid scale kinetic energy model for large-eddy simualtion of turbulent flows, J. Phys. Soc. Jpn. 54 (1985) 2834–2839. [23] F.R. Menter, Two equation eddy viscosity turbulence models for engineering applications, AIAA J. 32 (1994) 1589–1605. [24] R.A. Baurle, C.J. Tam, J.R. Edwards, H. Hassan, Hybrid simulation approach for cavity flows: blending, algorithm, and boundary treatment issues, AIAA J. 41 (2003) 1463–1484. [25] M. Sussman, P. Smereka, S. Osher, A level set approach for computing solutions to incompressible two-phase flow, J. Comp. Phys. 114 (1994) 146–159.
383
[26] M.B. Sun, Z.G. Wang, X.S. Bai, Assessment and modification of sub-cell-fix method for re-initialization of level-set distance function, Int. J. Numer. Methods Fluids 62 (2010) 211–236. [27] G. Russo, P. Smereka, A remark on computing distance functions, J. Comput. Phys. 163 (2000) 51–67. [28] D. Hartmann, M. Meinke, W. Schröder, Differential equation based constrained reinitialization for level set methods, J. Comput. Phys. 227 (2008) 6821–6846. [29] M. Ihme, Y.C. See, Large-eddy simulation of a turbulent lifted flame in a vitiated co-flow, in: 47th AIAA Aerospace Sciences Meeting 2009-2239, 2009. [30] F. di Mare, W.P. Jones, K.R. Menzies, Large eddy simulation of a model gas turbine combustor, Combust. Flame 137 (2004) 278–294. [31] H. Pitsch, 〈http://www.stanford.edu/group/pitsch/FlameMaster.htm〉, 2008. [32] U. Maas, J. Warnatz, Ignition processes in hydrogen-oxygen mixtures, Combust. Flame 74 (1988) 53–69. [33] V.E. Terrapon, F. Ham, R. Pecnik, H. Pitsch, A flamelet-based model for supersonic combustion, in, Center for Turbulence Research, Annu. Res. Briefs (2009). [34] K.J. Nogenmyr, P. Petersson, X.S. Bai, A. Nauert, J. Olofsson, C. Brackman, H. Seyfried, J. Zetterberg, Z.S. Li, M. Richter, A. Dreizler, M. Linne, M. Aldén, Large eddy simulation and experiments of stratified lean premixed methane/ air turbulent flames, Proc. Combust. Inst. 31 (2006) 1467–1475. [35] M. Oevermann, Numerical investigation of turbulent hydrogen combustion in a SCRAMJET using flamelet modeling, Aerosp. Sci. Technol. (2000) 463–480. [36] R. Guerra, W. Waidmann, C. Laible, An experimental investigation of the combustion of a hydrogen jet injected parallel in a supersonic air stream, in: 3rd International areospace planes conference, AIAA Paper 1991-5102, 1991.
Contents lists available at ScienceDirect
Acta Astronautica journal homepage: www.elsevier.com/locate/aa
The hybrid RANS/LES of partially premixed supersonic combustion using G/Z flamelet model Jinshui Wu a, Zhenguo Wang a,n, Xuesong Bai b, Mingbo Sun a,n, Hongbo Wang a a b
Science and Technology on Scramjet Laboratory, National University of Defense Technology, Changsha, 410073, China Division of Fluid Mechanics, Lund University, P. O. Box 118, S-221 00 Lund, Sweden
art ic l e i nf o
a b s t r a c t
Article history: Received 25 February 2016 Received in revised form 1 June 2016 Accepted 10 June 2016 Available online 11 June 2016
In order to describe partially premixed supersonic combustion numerically, G/Z flamelet model is developed and compared with finite rate model in hybrid RANS/LES simulation to study the strut-injection supersonic combustion flow field designed by the German Aerospace Center. A new temperature calculation method based on time-splitting method of total energy is introduced in G/Z flamelet model. Simulation results show that temperature predictions in partially premixed zone by G/Z flamelet model are more consistent with experiment than finite rate model. It is worth mentioning that low temperature reaction zone behind the strut is well reproduced. Other quantities such as average velocity and average velocity fluctuation obtained by developed G/Z flamelet model are also in good agreement with experiment. Besides, simulation results by G/Z flamelet also reveal the mechanism of partially premixed supersonic combustion by the analyses of the interaction between turbulent burning velocity and flow field. & 2016 IAA. Published by Elsevier Ltd. All rights reserved.
Keywords: Supersonic Combustion Partially premixed G/Z flamelet model
1. Introduction In scramjet engines, in order to achieve flameholding, different kinds of injection methods and flameholding devices are introduced. The fully mixing of fuel and air will take a long distance under supersonic condition, however, during the mixing process combustion will take place. This kind of combustion is partially premixed. As a combination of the premixed and the non-premixed mode of combustion, partially premixed combustion shows features of both. One of the most popular flame stabilization mechanisms in partially premixed combustion under low speed condition, i.e. the triple flame, considers that there exist a fuel lean premixed branch, a fuel rich premixed branch and a diffusion branch. In order to describe to this kind of combustion under low speed condition, many work have been done. Müller et al. [1] have developed a model for partial premixed turbulent combustion that they validate with lift-off heights in turbulent jet diffusion flames. Domingo et al. [2] have proposed a flame index to distinguish combustion mode in the entire flow field. Pierce et al. [3] developed a progress variable approach to account for flame quenching in non-premixed combustion. Pitsch and Ihme [4] have proposed an unsteady/ flamelet progress variable method in non-premixed turbulent combustion. However, all these methods are developed under low speed condition. For high speed condition, especially n
Corresponding authors. E-mail addresses: [email protected] (Z. Wang), [email protected] (M. Sun). http://dx.doi.org/10.1016/j.actaastro.2016.06.021 0094-5765/& 2016 IAA. Published by Elsevier Ltd. All rights reserved.
supersonic flow, whether these methods are valid is still unknown. There are two main issues considering the numerical simulation and combustion modeling of partially premixed supersonic combustion. Firstly, the combustion model chosen to describe the chemical reaction should be carefully considered. In partially premixed zones under supersonic condition, the range of the temperature distribution changes quite a lot. The maximum temperature sometimes could even reach 3000 K while the minimum temperature may be in the same magnitude with the temperature of cold flow. Under these conditions, the applicability of the traditional finite rate chemical reaction model is restricted. It may fail to effectively predict chemical reaction with broad temperature changes in partially premixed zones under supersonic condition by using finite rate chemical reaction model. Secondly, the flame propagation must be taken into account when there exists the partially premixed zones. Under supersonic condition, the flame propagation speed is much larger than low speed flow condition. The thickness of the premixed flame front could be thinner than a grid size of the LES. The grid resolution of LES is not capable of capturing the flame front without modeling. In order to investigate partially premixed supersonic combustion in scramjet engines, the framework of this essay is based on the model proposed by Müller et al. [1] which combines a level set method of describing premixed flame front and flamelet model. In this paper we have also drawn on many previous experiences of the LES of supersonic flow and combustion in our group [5–11]. A strutinjection supersonic combustion which has been experimentally investigated in German Aerospace Center was numerically simulated
376
J. Wu et al. / Acta Astronautica 127 (2016) 375–383
by G/Z flamelet model and finite rate chemistry model. In finite rate chemistry model, the reacting species transport equations together with their chemical source terms are solved, for simplicity the description of species transport equations are omitted and no turbulent combustion modeling is considered this essay. By our previous study [12], the maximum value of the Karlovitz number in main regions of this combustor is smaller than 1.0, indicating that the flame scales are smaller than the Kolmogorov scales so that the assumption of the flamelet model is satisfied under the given condition. Many other authors have also investigated this combustor numerically by flamelet model or finite rate chemistry model [13–17], however, they did not focus on the partially premixed supersonic combustion. Some simulation results could not well reproduce the characteristic of the partially premixed combustion under supersonic condition. In this essay, the practicability and accuracy of applying G/Z flamelet model into partially premixed supersonic combustion will be analyzed in the following part.
According to the discussion of the Kim et al. [21], modeling constant α ¼ 0.2 and β ¼ 6.67. u′ is local velocity fluctuation, in this paper the ksgs distribution is used to provide an estimate of the local velocity fluctuation u′ = 2k sgs/3 . In this essay, the SGS model chosen is the one-equation model of Yoshizawa and Horiuti [22]. Menter BSL (Baseline) model [23] is considered for the RANS portion. A similar blending function formulation [7] based on the hybrid RANS/LES methodology proposed by Baurle et al. [24] is used in this work. For brevity the governing equations for SGS model and RANS model and hybrid RANS/LES formula are omitted. 2.2. Re-initialization and a modified sub-cell-fix approach
∼ To obtain G equation in the whole field the so-called re-in∼ itialization equation ∇G = 1 is always used. A re-distancing technique introduced by Sussman et al. [25] is used to iterate on ∼ ∼ ∼ the equation Gτ = sign (G ) 1 − ∇G until a steady solution is
(
2. Model descriptions and numerical methods 2.1. Governing equations Flamelet model for premixed combustion in this paper is based on the level set approach using the scalar G. The distance function G(x, t) represents the distance of a certain point to the flame front. Although G equation model could well capture the interaction between large eddy structures and the premixed flame fronts, however, it could not be used to describe flame differences due to fuel concentration. Another scalar equation which considers diffusion effect, i.e. mixture fraction equation Z, is combined with G equation to account for partially premixed combustion zones. The Favre filtered compressible Navier-Stokes equations and flamelet equations [18] can be written as below:
∂ρ¯ ∂(ρ¯ u˜ i ) + =0 ∂t ∂xi ∂ ⎡⎣ ρ¯ u˜ i u˜ j + p¯ δij − τ¯ij − τijHybrid ⎤⎦ ∂(ρ¯ u˜ i ) + =0 ∂t ∂xj ∂ ⎡⎣ ρ¯ E˜ + p¯ u˜ i + q¯i − u˜ j τ¯ji + HiHybrid + σiHybrid ⎤⎦ ∂ρ¯ E˜ =0 + ∂xi ∂t Hybrid ∂ρ¯ Z˜ ∂(ρ¯ u˜ i Z˜ + Z˜ ) ∂ ⎛ ¯ ∂Z˜ ⎞ + = ⎜ ρ¯ D ⎟ ∂t ∂xi ∂xi ⎝ ∂xi ⎠ ∼ ∼ ∼ ∼ ∼ ∂ρ¯ G ∂(ρ¯ u˜ i G ) + = ρ¯ st ∇G − ρ¯ DtG κ c ∇G ∂t ∂xi
(
)
(1)
where is ρ¯ the filtered density, u˜ i is the filtered velocity, p¯ is the filtered pressure, E˜ is the filtered total energy, D Is the laminar thermal diffusion coefficient, τ¯ij is the filtered laminar viscous stress tensor, τijHybrid is the unclosed stress tensor, q¯i is the energy flux caused by heat conduction and mass diffusion, HiHybrid is the unclosed enthalpy flux vector, σiHybrid is the unclosed heat flux vector, is the unclosed mixture fraction flux vector which could be ∼ ∼ Hybrid expressed as Z˜ ≈ − ρ¯ DtZ ∇Z using gradient diffusion assump∼ ∼ tion and κc is filtered flame curvature ( κc = ∇⋅n = ∇⋅ −∇G / ∇G , n
(
)
is the normal direction vector). st is the turbulent burning velocity. The other parameters and unclosed items mentioned above can be found in reference [19]. In this essay, the turbulent burning velocity model proposed by Pocheau [20] is applied.
2.3. Partially premixed combustion modeling for supersonic flow Particularly, the iso-surface of G ¼ 0 has defined the position of the instantaneous premixed flame front, and divided the whole flow field into two region, where G 4 0 stands for the burnt region and G o 0 stands for the unburnt region. However, it has to be pointed out that for partially premixed flame, G equation has no meaning in regions exceeding lean or rich limit. However, for simplicity, G(x, t) ¼ 0 surfaces are still calculated in whole flow field, only in G 4 0 regions the flamelet libraries are applied. The flamelet model adopts a series of flamelet parameters to represent the chemical thermodynamic state, for example, the scalar dissipation rate is used to characterize the non-equilibrium effects, and the joint probability density function to characterize the interaction between the turbulence and the combustion. Considering the shock wave capturing algorithm in the supersonic flow, the temperature is calculated from the energy equation. The procedure is described as follows. The total energy could be written as NS
E=
∑ ρi (hi0 + hiT ) − p + i=1
1 ρ (u2 + v2) 2
(2)
(3)
hi0
where is non-standard enthalpy of formation which is independent of temperature, hiT is heat enthalpy which relates with temperature. A new total energy format is E′ introduced by neglecting, hi0 , NS
E′ =
∑ ρi hiT i=1
−p+
1 ρ (u2 + v2) 2 NS − ∑i = 1 ωi hi0
1
⎛ u′ ⎞α ⎤ α st ⎡ = ⎢1 + β⎜ ⎟ ⎥ ⎝ sL ⎠ ⎦ sL ⎣
)
reached. The converged steady solution gives the signed distance function. Efforts to improve the accuracy of distance function in level-set methods have led to a variety of re-initialization methods different from the original re-distancing algorithm. In this essay, a modified sub-cell fix method proposed by Sun et al. [26] based on Russo [27] and Hartmann [28] is applied. This modified sub-cellfix scheme independent of local curvature makes use of a combination of the points adjacent to zero level-set surfaces and preserves the interface in a second-order accuracy. The scheme is capable of handling large local curvature, and as a result it demonstrates satisfactory performance on several challenging test cases [26]. For brevity, the modified sub-cell-fix process will not be repeated here.
(4)
Then, a source term will appear in the new transport equation for E′. A time-splitting method is used to uncouple the flow and reaction as below.
J. Wu et al. / Acta Astronautica 127 (2016) 375–383
⎛ Δt ⎞ ⎛ Δt ⎞ Qn + 1 = L c ⎜ ⎟ L f (Δt ) L c ⎜ ⎟ Qn ⎝ 2⎠ ⎝ 2⎠
(5)
where Lf and L c denote the solver for flow and reaction, respecNS k tively. The filtered source term − ∑i = 1 ωihi0 appears in the reaction solver. During the reaction step, the new total energy will be updated due to reaction heat release. And the temperature will be calculated implicitly by newton iteration. The mass fractions of the species are interpolated from the flamelet libraries. The filtered ∼ mass fraction of the ith species Yi in the flamelet libraries is expressed by an ensemble of the laminar diffusion flamelets:
∼∼ Yi (Z , k Z″ 2, χ˜st ) =
˜ χ˜ ∬ Yi (Z˜ , χ˜st ) P (Z˜ , χ˜st ) dZd st
Cχ μ C μt 2 ∇Z˜ + 2Z t2 ε k Z″ 2 Δ Sct Cμ Sct
C χ Z = 2.0
Cε = 0.5 Cμ
(7)
k Z ″2 is calculated by scale similarity model [30]:
k Z ″2 = CZ Δ2 ∇Z˜ ⋅∇Z˜
Fig. 1. Sketch of the DLR scramjet combustor. Table 1 Flow conditions of the airstream and the hydrogen jet.
(6)
where k Z″ 2 is the filtered mixture fraction variance, χ˜st is the fil˜ χ˜ ) is the species mass fraction of tered scalar dissipation rate, Yi (Z, st ˜ χ˜ ) is the joint probability the laminar flamelet libraries and P (Z, st ∼ ∼ ˜ χ˜ ), Z and χ˜st density function of Z and χ˜st . When calculating P (Z, st are assumed to be statistically independent of each other, and they are subject to the Beta distribution and the Dirac function, respectively. In hybrid LES/RANS simulation, a direct calculation of scalar dissipation rate is not possible. Instead a filtered value must be calculated which contains both resolved and SGS-modeled parts [29]:
χ˜ = 2
377
(8)
where CZ is the model parameter and equals 0.2. Δ is the filtered scale. The flamelet libraries are generated by FlameMaster V3.9 software [31] using the mechanism containing 9 species and 19 step reactions [32]. Lewis numbers of these species are assumed to be unity. Radiation heat loss and pressure changes with time are not considered in the flamelet libraries. According to former study [33], the pressure boundary condition for generating flamelet libraries is 1.5 bar and the temperature of the oxidizer and fuel are 1000 K and 250 K, respectively. The size of the flamelet libraries is ∼ 1000*100*40 in Z , k Z″ 2 and χ˜ space. Besides, for some quenching regions where burning velocity is zero, when recirculation zone exists, G˜ (x, t) ¼ 0 surfaces will move with flow recirculation to form uncertain flame front. Therefore, ∼ on one side, G˜ equation should combine with Z to determine flame front; on the other side, a method is introduced into G˜ ∼ ∼ ∂ρ¯ G equation in quenching regions [34], where + ∇ρ¯ u˜ *G = 0, ∂t ˜1 , u˜2 , u˜3 ), u˜1 is the velocity point to the flow u˜ * = ( u˜1*, u˜2*, u˜* 3) = ( u recirculation direction. 2.4. Computational model and numerical method The sketch of the DLR scramjet combustor [35,36] is given in Fig. 1. In the numerical simulations of Ref. [35], in order to consider the shock capture scheme in high speed flow with the use of flamelet model, only the species mass fractions are interpolated from the libraries, the local temperature is implicitly calculated from the solution of the energy using the caloric equation of state. This method is also used in the modeling of partially premixed supersonic combustion in above section. In this combustor, vitiated air at a nominal temperature of 340 K flows over a strut. Hydrogen at sonic conditions exits through fifteen holes uniformly spaced
Air
H2
Ma T (K) P (kPa) YO2
2.0 340 100 0.232
1.0 250 100 0
YN2
0.736
0
YH2 O
0.032
0
YH2
0
1.0
along the back plane of the fuel injector, where it mixes with freestream air. The combustor was designed and experimentally investigated by German Aerospace Center. Experimental pressure distributions, axial velocity profiles, velocity fluctuation profiles, and static temperature profiles are available for both non-reacting (mixing) and reacting (hydrogen oxidation) cases. Detail descriptions of the combustor are omitted here. The original 15 holes are substituted by three holes to save the computational resources, and the computational domain has a length of 22 mm in the spanwise direction and periodic boundary condition. The bottom wall of the combustor has the coordinate of y ¼0 mm and the tip of the strut has the coordinates x ¼ 35 mm and y ¼ 25 mm. The structured mesh has approximately 6.7 million cells according to research earlier [13,14]. The mesh at the vicinity of the inlet, the strut, the wake flow region and the shear layer region is refined. The whole configuration is divided into 24 sections which are calculated in parallel. The flow conditions of the incoming airstream and the hydrogen jet are given in Table 1. In supersonic flow, the fluctuation in the downstream cannot propagate upstream. Therefore, the inflow boundary condition the boundary conditions for the hydrogen jet are fixed with the given parameters. The outflow is assumed to be supersonic such that its boundary condition is treated using extrapolation from the interior of the domain. According to the discussion of Génin et al. [14], the top and bottom walls of the chamber are assumed to be adiabatic and slip wall so as to avoid the interaction between the shock wave and the boundary layer. The surfaces of the strut are treated as no-slip and adiabatic walls. The fifth-order WENO (weighed essentially non-oscillation) scheme is used for the inviscid fluxes in order to improve the resolution of the flow field. The time integration is performed by means of a second-order accurate total-variation-diminishing (TVD) Runge-Kutta method. The Courant-Friedrichs-Lwey (CFL) number is fixed to be 0.5.
3. Results and discussions In our group, we have obtained many results about this combustor without and with combustion, and compared these results with experiment to test the validity and efficiency of our programming code for supersonic combustion [12]. In this essay, we
378
J. Wu et al. / Acta Astronautica 127 (2016) 375–383
Fig. 2. Perspective view in terms of contours of the axial momentum on the back plane ∼ ∼ ∼ and, an iso-surface of G = 0 colored by temperature and an iso-surface of Z = Zst .
only cared about the partially premixed combustion model we used and partially combustion phenomenon in this combustor, results without combustion will no longer be presented. Fig. 2 presents the perspective view of contours of the axial ∼ momentum, an iso-surface of G = 0 colored by temperature and an ∼ ∼ iso-surface of Z = Zst . It could be seen from Fig. 2 that the flow ∼ field has a great impact on the evolution of G equation, the interaction between jet and airstream in shear layer has wrinkled ∼ the iso-surface of G = 0 . The heat release in the recirculation zone behind the strut has increased the local pressure and expanded ∼ the recirculation zone, which also leads to the expansion of the G surface towards the airstream. Fig. 2 could also verify that the reinitialization and modified sub-cell-fix approach mentioned above ∼ for G equation are capable of capturing the surfaces of premixed front under complex computational model and complex flow condition. From the contour of the axial momentum on the back plane, the shock waves in the entire flow are well captured from ∼ ∼ calculation. It is shown in Fig. 3(a) that iso-surface of Z = Zst only covers a small region of the entire flow field, however, the iso∼ surface of G = 0 showed in Fig. 2 could extend to the entire region, most of the reactions occur on fuel lean condition. It could be seen ∼ in the Fig. 3(b) that the iso-surface of YOH = 0.1 indicates that there ∼ ∼ exists strong reaction zone in the vinicity of iso-surface of Z = Zst showed in Fig. 3(a). Besides, from the contours of the temperature on six cross planes showed in Fig. 3(c), we could see that the temperature in the second slice reaches maximum corresponding to the most intensive reaction zone compared with the other slices. The maximum temperature value is about 2500 K evaluated from Fig. 3(c). Fig. 4 further visualizes the instantaneous contour of tem∼ ∼ ∼ perature distribution together with Z = Zst line and G = 0 line at different spanwise locations. All high temperature zones at dif∼ ferent locations are surrounded G = 0 line. In the vicinity of the ∼ injector orifice, G = 0 lines drawn in the enlarged graph of Figs. 4 and 5 have engulfed the jet fluid and could extend to the thin mixing layer between the jet and the airstream. In the thin mixing layer after the strut, the mixture reaches the stoichiometric condition so that the turbulent velocity is large enough to overcome the flow speed in the shear layer and pushes the iso-surface of ∼ G = 0 toward to the strut edge. Due to the existence of fuel jet, iso∼ ∼ surface of G = 0 around the jet is broken into small G = 0 islands which will be transported to downstream as non-reacting circles as shown in both Figs. 4 and 5. As hydrogen jet develops downstream, mixing of the fuel and the air is promoted due to the largeeddy structures, mixture comes close to the stoichiometric state ∼ ∼ ( Z = Zst ) and mixing layers are broadened and the reaction zone moves from the upper and lower layers shown in Fig. 5 to the center of the flow field. It can be seen that the temperature in the center of the flow field achieves its peak value at the vicinity of
Fig. 3. Iso-surfaces and perspective view of computational results.(a) Iso-surface of ∼ ∼ ∼ Z = Zst (b) Iso-surface of YOH = 0.1 (c) Perspective view in terms of temperature on the six cross planes.
∼ ∼ line Z = Zst . In further downstream, combustion intensity decreases, reaction zone begins to shrink and breaks into small pieces and the temperature in the flow field decreases gradually due to continuous combustion consumption of fuel in the upstream. It also could be seen in the global graphs at different locations in Fig. 4 that the reaction wake flows are quite not the ∼ same due to the unsteadiness of the G = 0 surface, at Z/ZD ¼ 0.25, the reaction zone is close to the upper wall while at Z/ZD ¼ 0.75 the reaction zone is close to the lower wall. ∼ The turbulent burning velocity affected by Z is calculated using ∼ Eq. (2) and it combines with velocity to determine the evolution G and affects the combustion process indirectly. Fig. 6 gives the instantaneous contours of turbulent burning velocity and Mach number at different spanwise locations. It could be seen in Fig. 6 that the subsonic region downstream the strut mainly locates in the fuel rich region. Compared to white lines shown in Fig. 5, in the subsonic region, fuel is injected continuously to sustain the combustion and subsonic local speed condition is in favor of flame stabilization. It is illustrated in Fig. 6 that turbulent burning velocity differs from each other at different locations. Summarized from Fig. 6, the schematic of supersonic partially premixed combustion is also depicted in Fig. 7. We could see that turbulent burning velocity st1 (in shear layer near strut) are generally larger
J. Wu et al. / Acta Astronautica 127 (2016) 375–383
379
Fig. 4. Calculated front view of spanwise slice at Z/ZD ¼ 0.25, 0.5 and 0.75 (from up to down, ZD ¼22mm), instantaneous contour of temperature distribution together with ∼ ∼ ∼ Z = Zst line (white) and G = 0 line (black).
∼ Fig. 5. Snapshots of temperature contour evolvement together with G = 0 line (black) on central plane behind the strut at sixe different times ( Δt = 1.769 × 10−5s ).
than st2 (in subsonic zone) and st3 (close to the main flow). Under the interaction of turbulent burning velocity st1 and the recirculation zones generated near jet injection in the shear layer ∼ close to the strut, the iso-surface of G = 0 can propagate to the region very close to the strut edge. In subsonic zone, turbulent burning velocity st2 is small due to fuel rich condition, and the fuel ∼ injection velocity will push the iso-surface of G = 0 away from the strut. Together with the interaction of the recirculation zones ∼ generated by injection, the iso-surface of G = 0 is broken into very small iso-surfacs as shown in Figs. 4–6. In the region closes to the main flow, the turbulent burning velocity st3 is in competition with ∼ local velocity's projection on the iso-surface of G = 0 . Hot reaction zone almost exists in the region fuel/air mixture reaches
stoichiometry condition. Further downstream, the turbulent burning velocity decreases due to the consumption of the fuel and the large eddy structures generated have divided the iso-surface of ∼ G = 0 into small iso-surface islands. In these further regions, the flow field characteristics plays a more important role on the evo∼ lution of G equation. Fig. 8 presents time average temperature profiles using two different models compared with experiment at different locations. At x ¼78 mm, the peak value of temperature in the shear layer predicted by both models are higher than experimental data. This is also obtained by Berglund using two different flamelet models [13] and by Fureby et al. using different PaSR models [16]. There are two main factors that account for the over prediction of the
380
J. Wu et al. / Acta Astronautica 127 (2016) 375–383
∼ ∼ ∼ Fig. 6. Calculated instantaneous front view of spanwise slices at Z/ZD ¼ 0.25, 0.5 and 0.75 (from up to down, ZD ¼ 22 mm) together with Z = Zst line (black) and G = 0 line (purple).
temperature peak in the shear layer of the G/Z model. First, in the shear layer, the local equivalent ratio is close to stoichiometry equivalent ratio, so that the turbulent burning velocity is large ∼ enough to push the propagation of G = 0 surface forward to the edge of the strut as shown in Fig. 4. Then a more fierce combustion than experiment in the shear layer will be predicted. Second,
according to flamelet theory, if the local scalar dissipation rate is large enough, the flamelet solution will move close to the extinction. It is reasonable that the filtered scalar dissipation rate is under estimated so that the temperature peak obtained by flamelet model is larger than experiment. However, at this same location, the temperature profiles in the center zone after the strut predicted by the G/Z model are in better agreement with
J. Wu et al. / Acta Astronautica 127 (2016) 375–383
experiment than finite rate model. The finite rate model we used has over predicted the temperature which is similar with the results of Dharavath et al. [17]. In this zone, the combustion is partially premixed, the application of partially premixed model in this article could well reproduce this low temperature reaction zone by considering the interaction between G equation and Z equation. The simulation results obtained Berglund using two different
Fig. 7. Schematic of partially premixed supersonic combustion.
381
flamelet models have under predicted the temperature at this location. Further investigations by other authors show no better improvement [14–16], they have also under estimated the temperature in this partially premixed zone. At x¼ 125 mm, the width of the reaction zone predicted by G/Z model is better than finite rate model compared with experimental data. The results by finite rate model have broaden the reaction zone, this is also obtained by the work of Génin and Menon [14] and Dharavath et al. [17]. The peak value of the temperature predicted by both models shown in Fig. 8 are in good agreement with experimental data. At x¼ 233 mm, the width of the reaction zone and the peak value of the temperature predicted are smaller and lower than the experimental data compared with finite rate model. This will be explained in the following analysis. In order to analyze the influence of flow-field on the devel∼ opment of G and combustion, the average velocity and velocity fluctuation obtained by G/Z flamelet model and finite rate model are compared with experiment shown in Figs. 9–11. The average
Fig. 8. Average temperature compared with experiment using different combustion models at different axial locations.
Fig. 9. Average axial velocity compared with experiment using different combustion models at different axial locations.
Fig. 10. Average axial velocity fluctuation compared with experiment using different combustion models at different axial locations.
382
J. Wu et al. / Acta Astronautica 127 (2016) 375–383
Fig. 11. Average normal velocity fluctuation compared with experiment using different combustion models at different axial locations.
axial velocity predictions shown in Fig. 9 at x ¼ 78 mm in the shear layer are much smaller than the experiment by both models, even negative velocity values which could cause the generation of the recirculation are predicted. According to the ∼ ∼ properties of G equation, the surface of G = 0 could reach further upstream which causes the inaccurate prediction of reaction. This is consistent with the static temperature contours shown in Fig. 4. The average velocity fluctuation predictions in the shear layer are slightly over estimated by both models as shown in Figs. 10 and 11. The turbulent buring velocity are proportional to the velocity fluctuation. The over prediction of the velocity fluctuation could also contribute to the propagation of the surface of ∼ G = 0. At x ¼ 125 mm in Figs. 9–11, it is shown that deviations of average velocity and average velocity fluctuation compared with experiment by G/Z flamelet model are in more acceptable level than finite rate model so that G/Z flamelet model could well reproduce the overall combustion properties. Further downstream at x ¼233 mm, the average velocity fluctuation predicted by finite rate model is larger G/Z flamelet model, which may contribute to more intense reaction and higher temperature. The peak value of temperature and the width of the reaction zone predicted by G/Z flamelet model shown in Fig. 8 are smaller and lower than the experimental data. Besides the underestimation of velocity fluc∼ tuation, this may also be explained by the inaccuracy of G ∼ equation in the downstream for the surfaces of G = 0 are broken into pieces shown in Fig. 4 which could not sustain continous combustion.
4. Conclusions The partially premixed supersonic combustion in the DLR scramjet combustor is investigated numerically by hybrid RAN/LES using two different combustion models. Simulation results show that the combination of G equation method and flamelet libraries could capture the characteristics of partially premixed supersonic combustion. Further analyses of the comparsion between two models and experimental data show that the G/Z flamelet model used in this combustor could reproduce the low temperature reaction zone in the region after the strut. In this region, the predictions of temperature profiles by G/Z model are in better agreement with experiment than finite rate model and other simulation results, while the temperature in the shear layer are over predicted by the cause that filtered scalar dissipation rate might be under estimated. Although there exist some differences compared with experment in predicting the width of the reaction zone and temperature at further downstream, the deviations of average velocity and average velocity fluctuations compared with experiment are in acceptable level so that G/Z flamelet model could reproduce the overall combustion properties.
This paper validates the rationality and feasibility of the G/Z flamelet model in partially premixed supersonic combustion. On the one hand, it implies that the partially premixed combustion should be carefully dealt with. In order to reproduce the overall properties of the partially premixed combustion under supersonic condition, the modeling for numerical simulation is necessary. The research provides a guidance for the application of the G/Z flamelet model in partially premixed supersonic combustion. On the other hand, it points out there are some problems when applying the G/Z flamelet model in the shear layer near the strut and further downstream under supersonic condition. Some other modeling aspects should be taken into account for the quenching effect near the strut and elimating the strong dependence on G equation in further downstream.
Acknowledgments The research work was funded by the National Natural Science Foundation of China (Nos. 11472305, 11522222 and 91541101).
References [1] C.M. Muller, H. Breitbach, N. Peters, Partially premixed turbulent flame propagation in jet flames, Symposium (International) on Combustion, 25 (1994) 1099-1106. [2] P. Domingo, L. Vervisch, K. Bray, Partially premixed flamelets in LES of nonpremixed turbulent combustion, Combust. Theory Model. 6 (2002) 529–551. [3] C.D. Pierce, P. Moin, Progress-variable approach for large eddy simulation of non-premixed turbulent combustion, J. Fluid Mech. 504 (2004) 73–97. [4] H. Pitsch, M. Ihme, An unsteady/flamelet progress variable method for LES of nonpremixed turbulent combustion, in: Proceedings of the AIAA Paper 20052557, 2005. [5] H.B. Wang, Z.G. Wang, M.B. Sun, N. Qing, Combustion characteristics in a supersonic combustor with hydrogen injection upstream of cavity flameholder, Proc. Combust. Inst. (2013) 2073–2082. [6] M.B. Sun, H. Geng, J.H. Liang, Z.G. Wang, Mixing characteristics in a supersonic combustor with gaseous fuel injection upstream of a cavity flameholder, Flow., Turbul. Combust. 82 (2009) 271–286. [7] M.B. Sun, Z.G. Wang, J.H. Liang, H. Geng, Flame characteristics in a supersonic combustor with hydrogen injection upstream of a cavity flameholder, J. Propuls. Power 24 (2008) 688–696. [8] H.B. Wang, Z.G. Wang, M.B. Sun, N. Qin, Large eddy simulation of a hydrogenfueled scramjet combustor with dual cavity, Acta Astronaut. 108 (2012) 119–128. [9] H.B. Wang, Z.G. Wang, M.B. Sun, N. Qin, Large eddy simulation based studies of jet-cavity interactions in a supersonic flow, Acta Astronaut. 93 (2014) 182–192. [10] C.Y. Liu, Z.G. Wang, H.B. Wang, M.B. Sun, Numerical investigation on mixing and combustion of transverse hydrogen jet in a high-enthalpy supersonic crossflow, Acta Astronaut. 116 (2015) 93–105. [11] Y.X. Yang, Z.G. Wang, M.B. Sun, H.B. Wang, L. Li, Numerical and experimental study on flame structure characteristics in a supersonic combustor with dualcavity, Acta Astronaut. 117 (2015) 376–389. [12] Z.Q. Fan, W.D. Liu, M.B. Sun, Z.G. Wang, F.C. Zhuang, W.L. Luo, Theoretical analysis of flamelet model for supersonic turbulent combustion, Sci. China Tech. Sci. 55 (2012) 193–205. [13] M. Berglund, C. Fureby, LES of supersonic combustion in a scramjet engine
J. Wu et al. / Acta Astronautica 127 (2016) 375–383
model, Proc. Combust. Inst. 31 (2007) 2497–2504. [14] F. Génin, S. Menon, Simulation of turbulent mixing behind a strut injector in supersonic flow, AIAA J. 48 (2010) 526–539. [15] A.S. Potturi, J.R. Edwards, LES/RANS simulation of a supersonic combustion experiment, in: 50th AIAA Aerospace Sciences Meeting, AIAA Paper 2012-0611, 2012. [16] C. Fureby, E. Fedina, J. Tegner, A computational study of supersonic combustion behind a wedge-shaped flameholder, Shock. Waves 24 (2014) 41–50. [17] M. Dharavth, P. Manna, D. Chakraborty, Thermochemical exploration of hydrogen combustion in generic scramjet combustor, Aerosp. Sci. Technol. 24 (2013) 264–274. [18] N. Peters, Turbulent Combustion, Cambridge University Press, 2000. [19] M.B. Sun, Flow and flameholding mechanisms of cavity flameholders in supersonic flows (PhD thesis), National University of Defense Technology, 2008. [20] A. Pocheau, Scale invariance in turbulent front propagation, Phys. Rev. E 49 (1994) 1109–1122. [21] W.W. Kim, S. Menon, H. Mongia, Large eddy simulations of a gas turbine combustor flow, Combust. Sci. Technol. 143 (1999) 25–62. [22] A. Yoshizawa, K. Horiuti, A statistically-derived subgrid scale kinetic energy model for large-eddy simualtion of turbulent flows, J. Phys. Soc. Jpn. 54 (1985) 2834–2839. [23] F.R. Menter, Two equation eddy viscosity turbulence models for engineering applications, AIAA J. 32 (1994) 1589–1605. [24] R.A. Baurle, C.J. Tam, J.R. Edwards, H. Hassan, Hybrid simulation approach for cavity flows: blending, algorithm, and boundary treatment issues, AIAA J. 41 (2003) 1463–1484. [25] M. Sussman, P. Smereka, S. Osher, A level set approach for computing solutions to incompressible two-phase flow, J. Comp. Phys. 114 (1994) 146–159.
383
[26] M.B. Sun, Z.G. Wang, X.S. Bai, Assessment and modification of sub-cell-fix method for re-initialization of level-set distance function, Int. J. Numer. Methods Fluids 62 (2010) 211–236. [27] G. Russo, P. Smereka, A remark on computing distance functions, J. Comput. Phys. 163 (2000) 51–67. [28] D. Hartmann, M. Meinke, W. Schröder, Differential equation based constrained reinitialization for level set methods, J. Comput. Phys. 227 (2008) 6821–6846. [29] M. Ihme, Y.C. See, Large-eddy simulation of a turbulent lifted flame in a vitiated co-flow, in: 47th AIAA Aerospace Sciences Meeting 2009-2239, 2009. [30] F. di Mare, W.P. Jones, K.R. Menzies, Large eddy simulation of a model gas turbine combustor, Combust. Flame 137 (2004) 278–294. [31] H. Pitsch, 〈http://www.stanford.edu/group/pitsch/FlameMaster.htm〉, 2008. [32] U. Maas, J. Warnatz, Ignition processes in hydrogen-oxygen mixtures, Combust. Flame 74 (1988) 53–69. [33] V.E. Terrapon, F. Ham, R. Pecnik, H. Pitsch, A flamelet-based model for supersonic combustion, in, Center for Turbulence Research, Annu. Res. Briefs (2009). [34] K.J. Nogenmyr, P. Petersson, X.S. Bai, A. Nauert, J. Olofsson, C. Brackman, H. Seyfried, J. Zetterberg, Z.S. Li, M. Richter, A. Dreizler, M. Linne, M. Aldén, Large eddy simulation and experiments of stratified lean premixed methane/ air turbulent flames, Proc. Combust. Inst. 31 (2006) 1467–1475. [35] M. Oevermann, Numerical investigation of turbulent hydrogen combustion in a SCRAMJET using flamelet modeling, Aerosp. Sci. Technol. (2000) 463–480. [36] R. Guerra, W. Waidmann, C. Laible, An experimental investigation of the combustion of a hydrogen jet injected parallel in a supersonic air stream, in: 3rd International areospace planes conference, AIAA Paper 1991-5102, 1991.