Abstract
The statistics of mean fluid velocity components conditional in unburned reactants and fully burned products in the context of Reynolds Averaged Navier Stokes (RANS) simulations have been studied using a Direct Numerical Simulation database of statistically planar turbulent premixed flame representing the corrugated flamelets regime combustion. Expressions for conditional mean velocity and conditional velocity correlations which are derived based on a presumed bimodal probability density function of reaction progress variable for unity Lewis number flames are assessed in this study with respect to the corresponding quantities extracted from DNS data. In particular, conditional surface averaged velocities and the velocity correlations in the unburned reactants are demonstrated to be effectively modelled by the unconditional velocities and velocity correlations , respectively, for the major part of turbulent flame brush with the exception of the leading edge. By contrast, conditional surface averaged velocities and the velocity correlations in fully burned products are shown to be markedly different from the unconditional velocities and velocity correlations , respectively.
1. Introduction
Mean fluid velocities conditional in reactants, products, and instantaneous flame surface often play a key role in the theoretical analysis and modelling of turbulent premixed combustion [1–7]. Theoretical analysis of Bray et al. [1] demonstrated the relation between the turbulent scalar fluxes of reaction progress variable and the difference in the mean conditional velocities in products and reactants , where the overbar signifies a Reynolds averaging operation, , , and are the Reynolds fluctuation, Favre fluctuation, and Favre mean of a general quantity , is the ith component of velocity, is the gas density, and the subscripts and are used to refer to conditional mean quantities in reactants and products, respectively. Domingo and Bray [4] showed the relation between the slip velocity and the pressure transport terms in the Reynolds stress transport equation. In addition to the conditional mean velocities (i.e. and ), the surface averaged quantities play a pivotal role in the closure of turbulent premixed flames using the Flame Surface Density (FSD) [3, 5, 8, 9] or conditional mean equations [6, 7, 10–13].
In the case of an asymptotically high Damköhler number , which characterises a ratio of the large-scale turbulent time scale to the chemical time scale, the conditional mean velocities may be extracted from the Favre mean velocities, which are readily available in conventional Reynolds Averaged Navier Stokes (RANS) simulations. Under , the probability density function (pdf) of reaction progress variable could be approximated by a bimodal distribution and analytical expressions for conditional-velocity-related quantities could easily be derived based on a presumed bimodal pdf approach [1, 2]. In many practical flames, the underlying combustion process can be characterised by a high but finite within the corrugated flamelets regime on the regime diagram [14, 15], where the flame thickness remains smaller than the Kolmogorov length scale and the pdf of reaction progress variable is expected to be bimodal.
However, in spite of the key importance of conditional velocities in turbulent premixed combustion modelling, there is yet to be a detailed Direct Numerical Simulation (DNS) based assessment of the validity of the expressions of the mean conditional velocities, which can be derived based on a presumed bi-modal pdf of , in the corrugated flamelets regime combustion.
In this respect the main objectives of this study are as follows.(1)To analyse the statistical behaviours of the conditional mean velocity and their relation to the Favre mean velocity in the context of RANS for the corrugated flamelets regime combustion.(2)To analyse the statistical behaviours of the conditional surface averaged velocities and their relation to the conditional mean velocities for the corrugated flamelets regime combustion.
The rest of the paper will be organised as follows. The necessary mathematical background will be provided in the next section of the paper. This will be followed by a brief discussion on numerical implementation. Following this, results will be presented and subsequently discussed. Finally the main findings will be summarised and conclusions will be drawn.
2. Mathematical Background
DNS of turbulent combustion should account for both three-dimensionality of turbulence and detailed chemical mechanism. However, the limitation of computer storage capacity until recently restricted DNS either to two dimensions with detailed chemistry or to three dimensions with simplified chemistry. As turbulent velocity field is inherently three-dimensional in nature and vortex stretching mechanism is absent in two dimensions, the second approach takes precedence in the present study, which is based on a single-step irreversible Arrhenius-type chemistry.
In premixed combustion the species field is often represented in terms of a reaction progress variable which rises monotonically from zero in pure reactants to unity in fully burned products. The reaction progress variable can be defined in terms of the mass fraction of a suitable reactant in the following manner: where subscripts 0 and indicate values in pure reactants and in fully burned products, respectively. According to BML analysis [1, 2], the pdf of reaction progress variable and the joint pdf of velocity vector and are given by where , and are the coefficients for the progress variable pdf, and refer to pdfs of in unburned gases and fully burned products, respectively, and the functions and originate due to burning gases at the interior of the flame. According to BML analysis [1, 2], it is assumed that , , and are normalised in such a manner that they individually integrate to unity. The last term on the right-hand side of (2a) and (2b) refers to the contribution of burning fluid. Based on this pdf, one obtains where is the mean burned gas density, is the unburned gas density, and is the mean value of density conditional on . For low Mach number adiabatic unity Lewis number flames and are given by where is the heat release parameter with and being the adiabatic flame temperature and unburned gas temperature respectively. Under high Damköhler number (i.e., ) the contribution of becomes negligible in (3), which yields the following expressions of and : Based on (2b) and (5) one obtains the following relation for the Favre mean velocity components [1, 2]: The conditional mean values in unburned reactants and fully burned products (i.e., and ) for a general quantity dependent on are evaluated aswhere is a small number (i.e., ). It is possible to obtain the following relation for turbulent scalar flux using (2b), (5), and (6a) [1, 2]: Using (6a)-(6b) and (7) one obtains the following relations for and : Using (2b) it is possible to express and in the following manner: The expressions for Reynolds stresses conditional in reactants and products (i.e., and ) can be obtained using (7) and (9) in the following manner:The surface averaged value of fluid velocity component based on the generalised Flame Surface Density (i.e. ) [8, 9] can be expressed in the following manner: As proposed in [10] based on the DNS data by Im et al. [6], the simplest possible model for the surface averaged velocity components conditional on unburned gas side of flamelets (i.e. ) can be given as follows: In the reference frame attached to the flame, the conditional mean velocity approaches the turbulent flame speed (i.e., at the leading edge for a hypothetical fully developed, constant-density, statistically planar, one-dimensional flame that propagates from right to left. Under this condition, the surface averaged velocity component approaches (i.e., ) at the leading edge because the instantaneous propagation velocity of a flamelet should be equal to zero in the selected reference frame under steady state (otherwise the flamelet would cross the leading edge, and that is impossible by definition). Therefore, (11b) yields wrong result at in this simple case. For these reasons, the following linear interpolation has also been proposed in [10]: Here, and in which is either equal to or equal to . This model effectively ensures that when either or . Similarly when either or according to (11c). In deriving (11c), it is assumed [10] that, at the leading and trailing edges, (i) the flamelet is parallel to the flame brush (i.e., ), (ii) the velocity component tangential to the flamelet is assumed to be unaffected (i.e. ), and (iii) the velocity component normal to the flamelet is taken to be equal to laminar burning velocity (i.e., and ).
The simplest possible model for the conditional surface-weighted mean of the correlation between velocity components (i.e., ) can be proposed as [10] whereas the following linear interpolation [10] ensures the correct behaviour of at for a constant density “flame”:
3. Numerical Implementation
For the present study, a decaying turbulence DNS database of statistically planar freely propagating turbulent premixed flame representing the corrugated flamelets regime combustion has been considered. The initial values of normalised root-mean-square value of turbulent velocity fluctuation , integral length scale normalised by the thermal flame thickness , heat release parameter , global Lewis number , Damköhler number , Karlovitz number , and are given in Table 1, where is the unstrained planar laminar burning velocity, and is the thermal laminar flame thickness, and are the unburned gas density and viscosity respectively. The subscript is used to denote the quantities in the unstrained planar steady laminar flame. It is evident from Table 1 that the Karlovitz number remains smaller than unity in the case considered here which suggests that combustion in this case belongs to the corrugated flamelets regime [14, 15]. Standard values are chosen for Zel’dovich number , Prandtl number , and the ratio of specific heats (i.e.,) for the case considered in the present study, where is the activation temperature.
The DNS database considered here is taken from [16] where the simulation is carried out using low Mach number assumption. A domain of size is discretised by a Cartesian grid of with uniform grid spacing in each direction where is the thermal diffusivity in the unburned gas. In this case, inlet and outlet boundaries are specified in the mean direction of flame propagation. Realistic turbulent velocity fluctuations are specified at the inlet boundary using Taylor’s hypothesis by scanning a plane through a precomputed box of frozen turbulence. The inlet mean flow velocity is taken to be equal to . The outlet boundary is taken to be partially nonreflecting in nature and is specified according to Navier Stokes Characteristic Boundary Condition (NSCBC) formulation [17]. The transverse directions are assumed to be periodic.
A sixth-order finite difference scheme is used for evaluating spatial derivatives in the direction of flame propagation. The spatial derivatives in the transverse direction are evaluated using a pseudospectral method. The density change is accounted for by the relation between density and reaction progress variable according to the BML formulation for the unity Lewis number flames [1, 2]. In this case, the initial velocity field is specified by an initially homogeneous isotropic turbulence using a pseudo-spectral method following the Batchelor-Townsend spectrum [18]. The time advancement for viscous and diffusive terms in this case is carried out using an implicit solver based on Crank-Nicholson scheme whereas the convective terms are time advanced with the help of a third-order low storage Runge-Kutta method [19]. The flame is initialised using a steady unstrained planar laminar flame solution. The grid resolution is determined by the resolution of the flame structure. About 10 grid points are placed within the thermal flame thickness .
In this case, flame-turbulence interaction takes place under decaying turbulence. The simulation in this case was run for about 4 initial eddy turn over times which is about 27-fold greater than the chemical time scale . It is recognised that the simulation time remains small but previous studies with similar or smaller simulation time provided useful information in the past [8, 9, 20–26]. The values of in the unburned reactants ahead of the flame at the time when the statistics were extracted decreased by about 78%. The value of has increased from its initial value by a factor of about 1.62, but there are still enough turbulent eddies on each side of the computational domain. For the final values of and the underlying combustion situation belongs to the corrugated flamelets regime according to the combustion regime diagram [14, 15].
For the purpose of postprocessing DNS data, the mean quantities are assumed to be a function of the coordinate in the direction of mean flame propagation ( direction in the present case). The Reynolds/Favre averaged quantities are evaluated based on averaging the relevant quantities in transverse directions (). The statistical convergence of the Reynolds/Favre averaged quantities is checked by comparing the corresponding values obtained using half of the sample size in the transverse directions with those obtained based on full available sample size. Both the qualitative and quantitative agreements between the values obtained based on full and half sample sizes are found to be satisfactory. In the next section the results obtained based on full sample size available in transverse directions will be presented for the sake of brevity.
In the present study the conditional mean quantities are evaluated using (6b), where is taken to be 0.1 following the previous analysis by Domingo and Bray [4]. This essentially suggests that the mean quantities conditional in reactants are evaluated using the data corresponding to . Similarly the mean quantities conditional in products are evaluated using the data corresponding to . A smaller value of yields qualitatively similar results to those obtained using but the sample size for evaluating the conditional mean values decreases with decreasing . By contrast, increasing value (e.g., ) gives rise to similar qualitative trends to obtained using but increasing the value of by a large margin is likely to lead to and values which are not representative of conditional means in reactants and products, respectively. The quantitative agreement between the conditional velocity statistics obtained for is indeed found to be excellent.
4. Results and Discussion
The contours of the reaction progress variable for the flame considered here are plotted in Figure 1(a) which shows that contours remain mostly parallel to each other because the flame thickness remains smaller than the Kolmogorov length scale in this case and thus turbulent eddies cannot penetrate into the flame structure while background turbulence only wrinkles the flame. This is consistent with the behaviour expected in the corrugated flamelets regime [15]. However, a careful examination of Figure 1(a) reveals that the contours of progress variables for small values of (i.e., ) are not perfectly parallel to the contours representing the reaction zone (i.e. ) for the present thermochemistry. This suggests that the preheat zone is relatively more affected by turbulence than the reaction zone, indicating that some eddies of the order of the Kolmogorov length scale may perturb the preheat zone which is likely to give rise to minor departure from pure bimodal distribution of and .
(a)
(b)
(c)
The nature of pdf of can be characterised in terms of the variance of reaction progress variable . For a bimodal distribution of , the variance of reaction progress variable is given by The variations of with are shown in Figure 1(b) which indicates that remains close to for this case. This implies that the assumption of bimodal distribution of , as done in BML analysis [1, 2], holds reasonably well for the case considered here. The pdf of at the location corresponding to is shown in Figure 1(c), which clearly demonstrates that the bimodal distribution of satisfactorily approximates the pdf of in the flame considered here. However, Figure 1(b) clearly indicates that remains smaller than , which demonstrates that contribution is small but not negligible in this case. This behaviour is expected as the reaction progress variable contours representing the leading edge of the preheat zone are affected by the turbulent eddies and introduce nonnegligible contributions.
The variation of with is shown in Figure 2(a), which demonstrates that increases from unburned gas side to the burned gas side of the flame brush. The increase in due to thermal expansion and the effects of flame propagation are responsible for the increase in with increasing . However, the comparison between and in Figure 2(a) reveals that there is a noticeable difference between and . The variations of and are not shown because these components are identically equal to zero for statistically planar flames. It is evident from Figure 2(a) that the relation given by (8a) satisfactorily captures the variation of with for the flame considered here. As (8a) is derived based on the assumption of the bimodal pdf of , this expression satisfactorily captures for this flame where is roughly bimodal in nature.
(a)
(b)
(c)
(d)
The variation of with is shown in Figure 2(b). In this case also increases from unburned gas side to the burned gas side of the flame brush and there is a noticeable difference between and in this flame. The effects of increase in due to thermal expansion and flame propagation are also responsible for the increase of from unburned gas side to the burned gas side of the flame brush. The variations of and are not shown because these components are identically equal to zero for statistically planar flames. In this case the pdf of remains roughly bimodal and thus (8b) satisfactorily predicts obtained in the DNS (see Figure 2(b)).
Comparing Figures 2(a) and 2(b) it can be seen that remains greater than for this case. This behaviour can be observed clearly from the variation of with across the flame brush as shown in Figure 2(c). According to (7) a positive (negative) value of leads to a countergradient (gradient) transport of turbulent scalar flux , which can be substantiated from the variations of shown in Figure 2(d), where a positive (negative) value implies countergradient (gradient) type transport. In the present case, the sign of the normalised slip velocity indicates the direction of turbulent scalar flux. A similar observation was made earlier by Veynante et al. [24] and Chakraborty and Cant [27]. The modelling of has been discussed in [24, 27–29] and thus will not be further addressed in this paper for the sake of conciseness. It can be seen from Figure 2(c) that the value of remains close to, but smaller than, for a major portion of the flame brush (i.e., ). This behaviour is found to be consistent with the experimental and DNS based findings reviewed in [30, Section 3.3]. Veynante et al. [24] proposed a simple estimate for the normalised slip velocity as , where is an appropriate efficiency function which assumes a magnitude of the order of unity. The second term in the square bracket on the right-hand side is associated with the velocity jump induced by the pressure gradient along the normal to the laminar flame. For the present database the velocity jump remains much greater than (i.e., ), which ultimately gives rise to for the major part of the flame brush (see Figure 2(c)). Because , the fact that the ratio of is smaller than unity in the present case is associated with due to random orientations of the instantaneous flame normal vectors with respect to the -axis. As remains positive for the major part of the flame brush (see Figure 2(c)), assumes smaller (greater) values than in the current case (see Figures 2(a) and 2(b)).
The variations of with are shown in Figure 3(a). It is evident from Figure 3(a) that decays from unburned gas side to burned gas side in the present case. The variations of and the prediction of (10a) are also shown in Figure 3(a). The difference between and is significant for this case, which suggests that the contribution of cannot be ignored while evaluating . The variations of and with are shown in Figure 3(b). As and are statistically similar to and , respectively, the variations of and are not shown here for the sake of conciseness.
(a)
(b)
The evolution of and within a statistically planar one-dimensional turbulent flame brush (that propagates from right to left) is affected by several physical mechanisms [30–32], which can be evident from the following transport equations: The transport equation for is not shown here because it takes the same form as that of (14b) because and directions are statistically similar. The statistical behaviours of the various terms of (14a) and (14b) for this database have been discussed elsewhere [31, 32] in detail and will be not be repeated here for the sake of brevity. The main behavioural trends will be summarised here. The mean dilatation term is a sink in (14a) and acts to reduce but this term is identically equal to zero for (14b) in the context of statistically planar flames. The mean pressure gradient acts as a source term in (14a) for the present database due to countergradient transport (), which serves to increase . However, the term vanishes for (14b) in the case of statistically planar flames. The pressure dilatation term acts as a source term in both (14a) and (14b). The term acts as a sink due to viscous dissipation which serves to decrease and . The terms and denote the dispersion due to pressure fluctuation and velocity fluctuation, respectively, and these terms remain small in comparison to the magnitude of the contributions of , and ( and ) in the context of () transport. The generation mechanisms associated with and in the context of transport dominate over the dissipative actions of and in the middle of the flame brush, whereas the contributions of and dominate at the leading and trailing edges. As a result of this, decreases from the leading edge of the flame brush before increasing significantly at the middle of the flame brush, and it eventually decreases as the burned gas side is approached. By contrast, and roughly balance each other for the major part of the flame brush in the context of transport, but the dissipative action of dominates over the contribution of towards the leading and trailing edges of the flame brush. This gives rise to the decay of from the leading edge of the flame brush but it does not change significantly for the major part of the flame brush before decaying again towards the burned gas side. The smaller values of and in comparison to and , respectively, are associated with the influence of the second and third terms on the right-hand side of (10a). For instance, the third term on (10a) (i.e., is negative and its magnitude increases with when (10a) is written for . However, the contribution of the third term on the right-hand side (i.e. vanishes when (10a) is written for and the difference between and in Figure 3(b) is significantly smaller than the difference between and in Figure 3(a).
It is surprising that poor agreement (at ) between the DNS data and (10a) is observed in this case that is, under the approximate conditions for that the equation was derived. It is evident from Figure 3(a) that (10a) in this case significantly overpredicts the magnitude of towards the burned gas side of the flame brush. Although the pdf of remains roughly bimodal (see Figure 1(c)), the departure from bimodal distribution can be substantiated from Figure 1(b) which shows that remains slightly smaller than . The departure from bimodal distribution seems to have major influences on the predictive capability of (10a) because the reactive contribution to is not accurately represented by when departs from the bimodal distribution.
The variations of and with are shown in Figure 4(a) along with the prediction of (10b). It can be seen from Figure 4(a) that (10b) overpredicts the magnitude of for the major portion of the flame brush. Figure 4(a) also shows that the variation of remains markedly different from for the case considered here due to the contribution of . The variations of and with are shown in Figure 4(b). The variations of and with are both qualitatively and quantitatively similar to and , respectively, and thus are not shown here for the sake of brevity.
(a)
(b)
Figures 3(b) and 4(b) indicate that the predictions of both (10a) and (10b) agree reasonably well with the corresponding quantities obtained from DNS data for the major part of the flame brush. Although the differences and are of the same order, the relative difference is significantly greater than the relative difference . In an actual RANS simulation the quantities and are modelled and the accuracy of their modelling is also likely to affect the performance of the models given by (10a) and (10b). Modelling of and has been discussed elsewhere [1, 24, 27, 30, 33, 34] and thus will not be repeated here for the sake of brevity.
The variations of and with along with the predictions of (11c) according to and are shown in Figure 5(a). It is worth stressing that these velocities were evaluated in the coordinate framework selected so that the mean velocity of unburned gas upstream the flame brush is equal to . It is evident from Figure 5(a) that the simplest model given by (11b) [10] captures the behaviour of for the major portion of the flame brush, which is in agreement with earlier DNS results [6, 7]. Figure 5(a) show that the performance of (11c) is comparable for and and the model given by (11c) underpredicts for the case considered here. It is worth reminding that (11c) is a linear interpolation between two limit cases, that is, at (or ) to at (or ). Figure 5(a) indicates that the linear interpolation does not capture the behaviour of in the major part of the flame brush, but, at the leading edge of the flame brush (), (11c) agrees with the DNS data in a better manner than the simple model given by (11c).
(a)
(b)
(c)
Lee and Huh [7] proposed a model for in the following manner: where is a tuning parameter that is proportional to but significantly larger than the kinematic eddy viscosity conditional in unburned gas [7] where is a model constant. Lee and Huh [7] found that the ratio was equal to 6.7 and 20 in the two flames studied by them. The prediction of (17) for is also plotted in Figure 5(a), which shows that remains much smaller than and and the prediction of (15) remains comparable to (11c) for the major part of turbulent flame brush with the exception of the leading edge. The variations of and with are plotted in Figure 5(b) which indicates that satisfactorily captures the qualitative behaviour of . The quantitative agreement between and remain satisfactory towards the unburned gas side of the flame brush but the quantitative agreement is not satisfactory towards the burned gas side of the flame brush.
The variations of and with across the flame brush are shown in Figure 5(c) which demonstrates that remains close to but smaller than , which is also consistent with DNS results of Im et al. [6]. A physical mechanism that makes greater than is associated with the acceleration of burned gas under the action of the mean pressure gradient from the flame surface to the burned gas side of the turbulent flame brush. This physical mechanism has been discussed in detail elsewhere [12] and will therefore not be elaborated here for the sake of brevity.
The variations of and with are shown in Figure 6(a) along with the predictions of (12a) and (12b) (for and ). The corresponding variations for and with are shown in Figure 6(b). It is evident from both Figures 6(a) and 6(b) that and are appropriately captured by and , respectively. A comparison between Figure 2(a) and Figure 3(a) shows that the magnitude of is much larger than the magnitude of . Therefore, the behaviour of appears to be well captured by .
(a)
(b)
It is evident from Figure 6(a) that (12b) for both and overpredicts the magnitude of for the major part of the flame brush. According to (12b) the quantity is bound between at (or ) to at (or ), and, as is greater than in this case because of greater values of and than and , respectively (see Figure 2 and compare Figures 3 and 4), the model given by (12b) overpredicts the magnitude of for the major portion of the flame brush. Moreover, in this case the value of remains greater than , which leads to an overprediction for the model given by (12b) whereas (12a) satisfactorily predicts throughout the flame brush.
The variations of and with are shown in Figures 7(a) and 7(b), respectively, along with the variations of and . In the present case overpredicts whereas underpredicts . The former trend is consistent with the trend demonstrated in Figure 5(c), that is, , and is explained when discussing Figure 5(c).
(a)
(b)
It is worth noting that the above analysis is carried out using a simplified chemistry-based DNS database representing the corrugated flamelets regime combustion for a moderate value of turbulent Reynolds number following several previous studies [16, 24, 25]. More analysis is needed based on detailed chemistry-based DNS and experimental data for higher values of turbulent Reynolds number to develop more efficient models for , , and .
5. Conclusions
The statistics of mean fluid velocities conditional in reactants and products in the context of Reynolds Averaged Navier Stokes simulations are analysed in detail using a DNS database of statistically planar turbulent premixed flame representing the corrugated flamelets regime of premixed turbulent combustion. It has been found that the conditional mean velocities and and the velocity fluctuation correlations and differ substantially from the corresponding Favre averaged quantities.
While the conditional velocities and evaluated from DNS data agree well with the predictions obtained from the BML model, the conditional Reynolds stresses and are not appropriately predicted by the expressions derived based on a presumed bi-modal distribution although the segregation factor remains close to unity. It is also demonstrated that the conditional surface averaged velocities and the velocity correlations can be effectively modelled by and , respectively, for the major part of the flame brush with the exception of the leading edge. However, and cannot be modelled by and , respectively.
As the present analysis has been carried out for moderate value of turbulent Reynolds number and simplified chemistry, detailed chemistry-based DNS and experimental data for higher values of will be necessary for more comprehensive physical understanding and the development of more efficient models. Moreover, it remains to be seen how these models perform for nonunity Lewis number and low Damköhler number conditions where the pdf of is unlikely to be bimodal in nature. Some of these issues will form the basis of future investigations.
Acknowledgments
N. Chakraborty gratefully acknowledges the financial support provided by EPSRC, UK. A. N. Lipatnikov gratefully acknowledges the financial support provided by Swedish Energy Agency and by the Chalmers Combustion Engine Research Centre (CERC).