Pressure e ﬀ ects on ﬂ ame structures and chemical pathways for lean premixed turbulent H 2 /air ﬂ ames: Three-dimensional DNS studies

This paper presents three-dimensional direct numerical simulations of lean premixed turbulent H 2 /air ﬂ ames over a range of pressures using a detailed chemical mechanism. E ﬀ ects of pressure on ﬂ ame front structures and heat release from pressure-dependent pathways are analysed. Under the same initial turbulence at di ﬀ erent pressures, the Kolmogorov length scale and local ﬂ ame thickness decrease with increasing pressure. Thinner and sharper structures are found on the ﬂ ame front at elevated pressures. As the pressure is increased from 1atm to 5atm, heat release is greatly enhanced at convex regions but weakened at concave regions of the ﬂ ame fronts, which indicates that the e ﬀ ect of Darrieus-Landau instability is becoming stronger. The correlation of heat release and fuel consumption is also strengthened as pressure is elevated. A main pressure-dependent heat release reaction, H+O 2 (+M)=HO 2 (+M), is found to contribute less to the total heat release with increasing pressure for turbulent ﬂ ames, which is contrary to the trend in laminar ﬂ ames. In the low temperature zones, this is due to the decreased H radical pool at elevated pressure. In the high temperature regions, the reaction is less competitive compared with H+OH+M=H 2 O+M, thereby reducing its contribution to the heat release.


Introduction
Flame/turbulence interactions are unsolved fundamental problems that have significant impact on the efficiency, emissions and overall performance of a wide range of combustion engines. The interactions are typically two-way: On the one hand, chemical reactions affect the flow properties and turbulence structures through gas expansion due to heat release [1,2]; on the other hand, turbulence affects chemical reactions both directly and indirectly [3]. The majority of studies on flame/turbulence interactions are conducted at atmospheric pressure. Nevertheless, the operation of combustion engines such as gas turbines and internal combustion engines is typically at elevated pressures. There has been a growing body of experimental work on combustion at elevated pressures [4][5][6][7][8][9][10][11][12]. However, it is still difficult and expensive to conduct experiments under high pressures.
With the availability of increasingly powerful supercomputers, direct numerical simulation (DNS) of turbulent reacting flow has become feasible and affordable. DNS with detailed chemical mechanisms provides a powerful tool for studying the fundamentals of turbulent combustion, and aids the development and validation of turbulent combustion models for engineering applications. Flame structures and propagation in premixed flames are widely investigated. Im and Chen [13] studied the effects of strain rate on local displacement speed of H 2 / air triple flames with detailed chemistries. Cecere et al. [14] simulated a turbulent slot flame and discussed the effects of H 2 enrichment on flame structures. The flames at different combustion regimes are also explored extensively. Aspden et al. [15] conducted a series of DNSs to study lean premixed H 2 /air combustion in a range of turbulence intensities and explored distributed burning. Carlsson et al. [16] investigated flame structures of lean premixed flames at high Karlovitz numbers and compared the influence of small-scale turbulence on H 2 / air and CH 4 /air flames. In addition, turbulence-chemistry interactions of lean premixed combustion have been addressed by some researchers. Aspden et al. [17] presented the decorrelation of fuel consumption and heat release at high turbulence intensity and Dasgupta et al. [18] discussed the modification of chemical pathways based on his results.
Most of the above work has been dedicated to flame structures and turbulence/chemistry interactions at the atmospheric pressure. An increasing effort of DNS studies of turbulent premixed flames under elevated pressures has been reported in the literature. Zhang et al. [19,20] conducted several two-dimensional (2D) DNSs to elucidate the complicated chemical and physical processes in a constant volume enclosure relevant to homogeneous charge compression ignition engines. Yu et al. [21] compared 2D and three-dimensional (3D) DNS results of the auto-ignition process of a lean H 2 /air mixture with temperature stratification in a constant volume enclosure. A delayed but more rapid ignition of the mixture has been found in the three-dimensional case than two-dimensional simulation. Most of them have employed periodic boundary conditions and are limited to 2D cases, due to the practical difficulties of performing DNS at elevated pressures and the excessive demand for computational resources as a result of very thin reaction layers that require extraordinarily fine spatial and temporal resolution. Dinesh et al. [22][23][24] extended 3D DNS to lean premixed expanding spherical flames at elevated pressures to address the effects of preferential diffusion, turbulence intensity and equivalence ratio on flame structure, propagation and instability. Although these studies provide much valuable information, DNS studies of transient flame characteristics with detailed chemical mechanisms at elevated pressures are far and few between. In particular, pressure-dependent reactions play an important role in heat release which is addressed for flames at atmospheric pressure [17,18,25], but the effects of pressure on theses reactions are not presented. It is thus critical to investigate the structures and chemical pathways of turbulent premixed flames at elevated pressures.
The aim of the present work is to enhance our fundamental understanding of the effects of pressure on turbulent flame structures and pressure-dependent reactions based on DNS of turbulent lean premixed H 2 /air flames. As following, the numerical method is described in section 2 and the computational cases at different pressures are given in section 3. DNS results are analysed in detail to elucidate the effects of pressure on turbulent flame structures and statistical properties in section 4, followed by conclusions in section 5.

Direct numerical simulation
The DNS is performed using the time-dependent compressible Navier-Stokes equations coupled with detailed transport properties and chemical kinetics [26]. The governing equations are: where ρ is the mixture density, u i velocity component in x i direction, p pressure and τ ij stress tensor. V k , Y k and ω̇are respectively diffusion velocity, mass fraction and reaction rate for specie k. E T is the total energy, q j jth-component of heat flux, R universal gas constant, T temperature and MW mixture molecular weight. Mixture averaged approximation is adopted for species diffusion, where Soret and buoyancy effects are neglected. Buoyancy has many effects on the propagation and instability of premixed and non-premixed flames [27,28]. However, thermo-diffusive and Darrieus-Landau instabilities play the most important roles in lean premixed combustion [23]. It has been demonstrated that buoyancy has a significant influence on laminar flames, but the effects on turbulent flames diminish with increased turbulence intensity [29] and the influence is limited for flames in small-scales [30]. Moreover, the characteristic time of flame propagation is much smaller than that of buoyant convection. As the considered cases are in small scales and under high turbulence intensities, the buoyancy effects should be small for the lean premixed H 2 /air flames in this study. Species diffusion velocity is expressed as [26]: where D k and X k are, respectively, the diffusion coefficient and mole fraction for species k.
The species reaction rate is given by [31]: where K fj and K rj are forward and reverse reaction rates for reaction j.
′ v kj and ″ v kj are forward and reverse stoichiometric coefficients for specie k in reaction j. The rate constant for pressure-dependent reactions (unimolecular/recombination fall-off reactions) is expressed as [32]: where K ∞ is the limit of high-pressure rate constant and P r is the pressure scale. F is determined with the TROE approach which is more sophisticated than the Lindemann form (where F is unity). The TROE approach is computationally more expensive but more realistic and compact enough in analysing large reaction systems [33,34].
In the present DNS of turbulent H 2 /air flames, three-dimensional, time-dependent, compressible Navier-Stokes equations and species transport equations are solved with a six-order explicit finite difference scheme. For time advancement, the low-storage third-order Runge-Kutta (RK3-2N) scheme is used for transport equations and third-order Livermore Solver for Ordinary Differential Equations (LSODE) is used for chemistry calculation. The time step in the simulations is controlled under 5 × 10 −9 s. This code was extensively used in astrophysics and magneto-hydrodynamics. Its application in premixed turbulent combustion could be found in [26,35].
The chemical mechanism for hydrogen oxidation adopted was developed by Li. et al. [36]. This mechanism contains 21 reversible reactions and 9 species (H 2 , O 2 , H 2 O, H, O, OH, HO 2 , H 2 O 2 , N 2 ). The chemistry implementation was validated quantitatively with Chemkin.

Simulation cases
Four 3D DNS cases were set up to simulate lean premixed H 2 /air flames under three different pressures in an inflow-outflow configuration with an aspect ratio 2:1:1. These cases have the same equivalence ratio = φ 0.6 and unburned gas temperature T u = 298 K. Table 1 lists the detailed parameters of the studied cases. Firstly, the domain was initialized by superimposing a one-dimensional laminar flame on a pregenerated turbulence. The flames were therefore developing in turbulence. The initial turbulence was isotropic and homogeneous which was generated with helical forcing functions [37]. Periodic boundary conditions were imposed in the span-wise and lateral directions. The Table 1 Parameters of simulation cases. For all the cases, the equivalence ratio is 0.6 and the temperature of the unburned gas is 298 K. .

X. Wang et al. Fuel 215 (2018) 320-329
Navier-Stokes Characteristics Boundary Conditions (NSCBC) were applied at the inlet/outlet to maintain a constant pressure. The domain length is 10 mm and width is 5 mm corresponding to 11 times the integral length scale l t in the crossflow directions. As shown in Fig. 1, four cases are within the thin reaction zone in the premixed combustion regime diagram. From Cases A to C, the pressure increases from 1 atm to 5 atm, while the root-mean-square turbulent fluctuation velocity ′ u is maintained the same. It is found from Table 1, when pressure is elevated, the laminar flame speed S L , thickness δ L and kinematic viscosity of mixture v decrease significantly. As a result, the initial Reynolds numbers increase from 153 to 763. However, the Karlovitz numbers are at the same level. Case D is set up with the same Reynolds number as Case A, while the corresponding Karlovitz number is much lower and the case is located close to the corrugated flamelets region.
The grids for each case are 1024 × 512 × 512 with uniform mesh resolution = x μ Δ 9. 77 m. Correspondingly, the number of grids across the flame thickness decreases from 38 to 8. In order to check that the cases are sufficiently resolved, a grid-independent-test case is performed. Case C with the thinnest flame is chosen as the baseline case. If the grid resolution is sufficient for Case C, then the resolution is good enough for the other three cases. All properties for the grid-independent-test case are the same as the baseline case except for resolution which decreases from 9.77 μm to 4.89 μm. Fig. 2 shows the comparison of the average heat release rate in the computational domain for the two cases. The two simulations agree with each other very well, which indicates that the mesh used here can produce grid-independent results.

Flame structures and local properties
Flame front is identified with the iso-surface which corresponds to the peak heat release rate in the turbulent flames. Fig. 3 shows instantaneous snapshots of flame fronts coloured by heat release rate (HRR) for the four cases. The convex structures (protruding to unburned mixture) are taken to be positive curvature and the concave structures (protruding to products) are taken to be negative one. When the flames are under the identical Reynolds number but different pressures (Case A and D), the flame structure shows significant differences. At p = 5 atm (Case D), the flame surface is slightly wrinkled, which corresponds well with Fig. 1 where Case D locates near the corrugated flamelets region. However, Case A shows a flame surface with considerable wrinkling. In this study, we will focus on the cases with the same initial velocity fluctuation but different pressures, which is representative of practical combustion [38,39]. Due to the hydrodynamic instability, the flames all display large scale wrinkling structures which are similar to those reported in [23]. It is noted that there are more small and sharp structures when pressure increases from 1 atm to 2 atm to 5 atm. This is due to the increased ′ u S / L ratio and decreased Kolmogorov length scale [39]. An interesting observation is that in high positive curvature regions the heat release rate is much lower than that in the surroundings at p = 1 atm, but the phenomenon is not obvious at p = 2 atm and 5 atm. In addition, high negative curvature areas show ultra-low heat release rate at elevated pressures.
Comparing concave and convex structures shown in Fig. 3, it is found that they have different characteristics. As the pressure increases, the concaves are becoming much thinner and sharper. To further investigate this phenomenon, probability density functions of mean flame surface curvature are plotted in Fig. 4. It is clear that there are more positive curvature areas and surface curvature is increasing when pressure is increased. For example, at p = 1 atm, the ratio of cumulative distribution functions (CDFs) for positive and negative curvature is 1.15:1, which is much smaller than 1.74:1 at p = 5 atm. This suggests more intermediates protrude into the unburned mixture at elevated pressures.
To further analyse the flame structures, two-dimensional slices of normalized heat release rate (NHRR, normalized by the peak heat release rate of corresponding 1D flame) taken from three-dimensional simulations for Case A, B and C are plotted in Fig. 5. Note only the regions that cover the flame fronts are plotted and local structures at p = 5 atm are also given. Fig. 5 also shows iso-lines of progress variable C T = 0.1-0.9 and NHRR = 35%NHRR peak which separately define the reaction zone and inner reaction zone [40]. The non-dimensional progress variable C T is given by: where T , T u and T p are the local, unburned and peak temperatures, respectively. It can be seen in Fig. 5 that the vorticity structures for the three cases are similar because these simulations are initialized with the same velocity field. However, the vorticity gradient is increasing with pressure due to the decreased kinematic viscosity, which tends to distort the flame front seriously. Even if the laminar flame thickness is decreased significantly with increased pressure, the flame front position variance is not obvious. It is also noted that, locally, the thinner reaction zone at elevated pressures makes it harder for small turbulence structures to penetrate into it. This observation confirms that the three cases are in the thin reaction zone as shown in Fig. 1.
A striking observation from Fig. 5 is that the flame front shows more finger structures with intense heat release rate when pressure is  elevated. In the low-pressure case, a cellular burning pattern is shown, with moderate heat release in regions of both positive and negative curvature. When pressure is increased to 5 atm, the cellular structures appear thinner and sharper; heat release in the convex regions is becoming quite different from that in the concave regions. It is clear from the zoomed-in views for p = 5 atm, heat release rate in regions of positive curvature is intense, but local extinction happens in regions of negative curvature. This is mainly due to the effect of Darrieus-Landau instability. In turbulent flames, any initial vortex serves as an initial disturbance for the onset of hydrodynamic instability, which leads to the wrinkling flame structure, together with turbulent eddies [41]. Thin flame reaction zone at elevated pressure will enhance the effect of Darrieus-Landau instability on local flame structure and increase topological complexity of the local flame. In convex regions, thermal Fig. 3. Instantaneous snapshots of flame front (T = 1000 K) coloured by heat release rate for Case A-D at three pressure values. The view for A1-D1 is protruding to unburned mixture (positive curvature regions) and the view for A2-D2 is protruding to products (negative curvature regions).
diffusion to the external streamlines is weaker than the chemical energy due to the deficient reactant transportation [42,43]. Local flame intensity will be enhanced. However, it is interesting to note that the peak heat release rate does not happen in regions of highest positive curvature, but in flank regions. This explains why the regions with highest positive curvature in Fig. 3 are related to lower heat release rate values. It can be seen from the iso-lines that when pressure is increased, the correlation between heat release rate and temperature is enhancing in regions with peak positive curvature; regions with peak temperature higher than burnt gases show in the centre of finger structures. These observations show the effects of pressure on flame structures and local heat release.

Variations of heat release and species distributions
It is known that flame wrinkling is affected by pressure and relationships among flame properties are changed by pressure. To quantitatively investigate the phenomenon, we look at the scatter plot of heat release rate and fuel consumption rate. The two quantities are normalized by the peak values of the corresponding 1D flames at different pressures. Fig. 6 shows the joint probability density function of temperature and normalized heat release rate at various pressures. The dashed lines provide references of 1D flames at the corresponding pressures. The distributions of 1D flames are close to normal distribution and the adiabatic temperatures are close to 1850 K, however the corresponding temperature of peak heat release region is increasing with pressure. For example, the peak heat release happens at a position near T = 1000 K at p = 1 atm. When pressure is elevated to 5 atm, the corresponding temperature region is at T = 1300 K. This observation is much more obvious for the 3D simulations in turbulent fields. In the two cases, the regions with T < 800 K are referred to as the low temperature zone of the reaction. There is noticeable heat release in the low temperature zone for the flame at p = 1 atm where the NHRR exceeds unity and even reaches a peak value as high as 2 in the domain. The intensified heat release in the low temperature zone is also observed in high Ka combustion [16] and this study suggests that the phenomenon is not unique for high Ka flames. It is interesting to note that when pressure is elevated, the heat release at low temperatures is decreased and the scatter distribution is enhanced towards the direction of deviating from 1D reference flames. It is reported in [25] that differential diffusion can enhance the heat release in the low temperature zone. This suggests that the influence of differential diffusion is obviously decreased for lean premixed flames when pressure is elevated. Fig. 7 shows the scatter plots and joint probability density function of temperature and local equivalence ratio for Cases A and C. The local equivalence ratio φ is given by:  where Y H and Y O are the molar mass of elements H and O, respectively. When the pressure is elevated to 5 atm, the distribution in the T − φ field becomes more scattered in the thin reaction zone due to the larger maximum curvature and stronger preferential diffusion for the tip shown in Fig. 4, but the data is more concentrated, towards the direction of yielding a unity effective Lewis number. However, findings in laminar flames show the opposite phenomenon: the deviation of local equivalence ratio from unburned mixture is stronger at p = 5 atm compared with that at p = 1 atm. It indicates that the rates of turbulent transport of energy and species are getting closer at the elevated pressure. Furthermore, at p = 5 atm, there is a high degree of scatter in the region T > 1850 K. This corresponds well with the high temperature regions in Fig. 5.
To examine the influence of pressure on fuel consumption and heat release, the scatter plots and joint probability density function of fuel consumption rate and heat release rate are shown in Fig. 8. The two variables are normalized by peak values of 1D flames, respectively. Owing to the decreased flame thickness, the distribution of scattered data changes from a triangle to an acute angle shape when pressure increases from 1 atm to 5 atm. Aspden et al. [17,44] observed a decorrelation of fuel consumption and heat release at high positive curvature regions and attributed it to the high mobility of atomic hydrogen at high turbulence intensity. In this study, when pressure is elevated, there are more regions of high curvature, but the decorrelation is clearly reduced. For Case A, more H radical diffuses from reaction zone to preheat zone. The consumption of H in the low temperature zone is associated with an exothermic reaction [25]: Owing to the above reaction, intensified heat release, which reaches peak values, is observed at low temperatures for Case A in Figs. 5 and 6. For Case C, the reduced decorrelation of heat release and fuel consumption is mainly due to the reduced mobility of atomic hydrogen, which is consistent with the analysis of Fig. 7. Comparing Case C and Case D, it is interesting to find that the correlation of heat release and fuel consumption shows a very similar trend. It implies the influence of turbulence on local flame chemistry is dramatically reduced when pressure is elevated. It becomes harder for turbulence to disrupt the flame inner structures.

The contribution of pressure-dependent reactions
In laminar flames, pressure-dependent reactions make a great contribution to total heat release. When the pressure is elevated, these reactions play a more important role. This section focuses on heat release from pressure-dependent reactions in turbulent flames. There are two pressure-dependent reactions in the Li et al. mechanism. Table 2 lists the main heat release reactions in which reaction 3 is a pressure-  dependent reaction. The heat release from another pressure-dependent reaction is very weak which is not discussed in this paper. The contribution of a given reaction to total heat release is computed as: where q i is heat release from reaction i and q tot is total heat release in the domain. The ratio of cumulative heat release from reaction i to total hear release is calculated as: where < < T T T u c is the temperature range in which heat release from reaction i is cumulated.
For a clear comparison, the heat release contribution of reaction 3 in laminar and turbulent flames is plotted in Fig. 9. A striking observation is that the contribution of reaction 3 shows different variation with pressure for laminar and turbulent flames. For the laminar flames, the contribution increases from 0.286 at 1 atm to 0.309 at 5 atm. A reverse trend is observed for the turbulent flames: the contribution drops from 0.304 to 0.282. This observation suggests a possible change in the local chemical pathways with pressure.
To investigate the reasons why the heat release contribution trend is reversed, we firstly plot the scatter distribution of temperature and heat   release rate of reaction 3 at 1 atm and 5 atm in Fig. 10. It is noted that heat release in the low temperature zone (T < 800 K) at p = 1 atm is very strong, and even reaches peak values. At P = 5 atm, the scatter distribution shifts towards the high temperature zone with the peak heat release value at around 1100 K. The reaction mainly occurs in high temperature regions and the distribution is similar to the total heat release shown in Fig. 6 because reaction 3 is the main exothermal reaction in the mechanism. Fig. 11 further shows the ratio of cumulative heat release to total heat release of two main H consumption reactions for both laminar and turbulent flames. For reference, Fig. 12 presents the scatter plot of H radical mass fraction at the two pressures for laminar and turbulent flames. For reaction 3 at p = 1 atm, the major increase of cumulative heat release ratio happens in the low temperature regions where the distribution of most H radical for the turbulent flame is above the level of the 1D laminar flame. The higher fraction of H can contribute to an intense heat release from reaction 3; this is the main factor responsible for the higher heat release contribution for the turbulent flame. Similarly, the heat release ratio is also increased for another H consumption reaction #2 in the low temperature zone. However, it is not enhanced to the same extent because there is excess O 2 in the low temperature regions and reaction 3 is more competitive. For reaction 3 at p = 5 atm, the contribution of heat release is apparently reduced in the low temperature zone. However, the ratio soars sharply when T > 1200 K for the laminar flame. In contrast, the ratio increases slowly for the turbulent flame, although a slight increase is observed in the low temperature regions. As can be seen in Fig. 12, the fraction of H radical stays at a very low level in low temperatures. It is enhanced for turbulent flames, but the scatter data density is much lower than that at p = 1 atm in the low temperature zone; most of the distribution is located in high temperature regions for the turbulent flame. There is also a high OH radical pool in the high temperature zone, which means reaction 2 is more competitive than reaction 3 at the high pressure. Correspondingly, the cumulative heat release ratio of the turbulent flame at 5 atm is enhanced greatly in high temperature regions. It is clear that the increased H distribution in the low temperature zone contributes to a higher heat release from reaction 3 at the low pressure and the enhanced H distribution at high temperatures results in more heat release from reaction 2 at the elevated pressure.

Conclusions
A series of three-dimensional direct numerical simulations have been performed for lean premixed turbulent H 2 /air flames in the thin reaction zone, with pressure ranging from 1 atm to 5 atm. Four flames are considered: Cases A-C are initialized with the same turbulence at three different pressures to study the influence of pressure on flame structures and chemical pathways; Case D is presented with the same Re number with Case A but a different pressure to isolate the effect of pressure on flames under the same turbulence intensity.
First, comparing Case D to Case A, an obvious decrease in flame surface areas and wrinkled structures is observed. Comparing Cases A-C where planar flames interact with the same turbulence, more sharp structures are observed at elevated pressures on the flame fronts. Heat  release is enhanced in concave regions, while reduced in convex regions, which can be attributed to the Darrieus-Landau instability. With increased pressure, the flame front area with positive curvature is increased. It is quite interesting to observe thin finger structures protruding into the unburned mixture under elevated pressures due to increased flame instability. It is also found that the peak heat release does not occur in fingertips but in flank regions. Second, the local heat release rate of turbulent flames is examined. It is observed that the peak value is over twice of that in the corresponding laminar flame. With increasing pressure, peaks are moving from low temperature regions to high temperature regions, which is more obvious for turbulent flames. Differential diffusion for turbulent flames is reduced with increased pressure, which results in an effective unity Lewis number. However, laminar flames experience an opposite trend. When pressure is elevated from 1 atm to 5 atm, the transport of H radical from high temperatures to low temperatures is reduced. The consequence is the enhanced correlation of heat release rate and fuel consumption rate.
Finally, contribution to heat release from the pressure-dependent reaction H + O 2 (+M) = HO 2 (+M) is analysed, considering it contributes around 30% of the total heat release in laminar flames at the atmosphere pressure. A striking observation is: the heat release contribution of this reaction increases with pressure for laminar flames, but decreases for turbulent flames. It is attributed to the decreased H radical pool in low temperature regions for turbulent flames at elevated pressures. In high temperature regions, the reaction is also less competitive compared with H + OH+M = H 2 O + M, thereby reducing its contribution to the total heat release. In summary, chemical pathways and heat release patterns are altered by pressure in turbulent flames compared with laminar flames.