Evaluation of scalar dissipation rate sub-models for modeling unsteady reacting jets in engines

Chemical Engineering Science 127 (2015) 334–343

Contents lists available at ScienceDirect

Chemical Engineering Science journal homepage: www.elsevier.com/locate/ces

Evaluation of scalar dissipation rate sub-models for modeling unsteady reacting jets in engines M.M. Ameen n, J. Abraham School of Mechanical Engineering, Purdue University, West Lafayette, IN 47906-2088, United States

H I G H L I G H T S







Evaluation of filtered scalar dissipation rate models for reacting LES. Strain rate tensor model is the most suitable model. Exponential and log‐normal PDFs for the scalar dissipation rate are evaluated. A model equation is formulated for the scalar dissipation rate variance.

art ic l e i nf o

a b s t r a c t

Article history: Received 15 December 2014 Received in revised form 13 January 2015 Accepted 24 January 2015 Available online 4 February 2015

The filtered scalar dissipation rate is an important variable required in many turbulent combustion models used in large eddy simulations (LES). In the present study, direct numerical simulations (DNS) of reacting and non-reacting turbulent mixing layers are used to evaluate different models for the filtered scalar dissipation rate. This study is conducted at elevated temperature and pressure conditions relevant to engine applications. The models evaluated are the turbulent diffusivity model, k–ε model, the strain rate tensor (SRT) model and the subfilter kinetic energy (SKE) model. It is found that the SRT model is the best choice when considering the performance and ease of implementation in LES codes. DNS results are also used to evaluate the marginal PDF for scalar dissipation rate. It is found that an exponential PDF works well for low values of scalar dissipation rate and smaller filter sizes whereas a lognormal PDF works well for larger values of scalar dissipation rate and larger filter sizes. A model is also formulated for the variance of the scalar dissipation rate to be used when employing the lognormal PDF by relating it to the mean and variance of the mixture fraction. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Modeling reacting jets Engine combustion modeling Scalar dissipation rate Flamelet model

1. Introduction With recent advances in computational capabilities, large eddy simulation of turbulent reacting flows has started to become feasible for conditions where the pressure and temperature are elevated such as in engine applications. The sub-grid scale modeling of turbulence–chemistry interactions is, however, a continuing subject of inquiry. For turbulent non-premixed combustion, flamelet-based models are attractive when the combustion occurs in reaction zones that are smaller than the Kolmogorov scale. This class of models includes the unsteady flamelet model (Pitsch et al., 1998), steady flamelet progress variable model (Pierce and Moin, 2004) and unsteady flamelet progress variable model (Ihme et al.,

n Correspondence to: School of Mechanical Engineering, Purdue University, 585 Purdue Mall, West Lafayette, IN 47906-2088, United States. Tel.: þ 1 765 494 5660. E-mail address: [email protected] (M.M. Ameen).

http://dx.doi.org/10.1016/j.ces.2015.01.055 0009-2509/& 2015 Elsevier Ltd. All rights reserved.

2005; Ihme and Pitsch, 2008; Ihme and See, 2010). When the reactions occur in thin zones, the effect of the turbulence is to locally stretch as well as extinguish these reaction zones. Under these conditions, the flame can be assumed to be locally onedimensional in the mixture fraction (Z) space. Within these reaction zones, which are called as flamelets, the evolution of the species mass fractions are governed by the unsteady flamelet equations given by ∂ϕ χ ∂2 ϕ _ ϕ; ¼ þω ∂t 2 ∂Z 2

ð1Þ

where ϕ is a vector representing the set of all reactive scalars, which includes the temperature and mass fractions of all the _ ϕ is the corresponding source term due to chemical species, and ω reactions. The symbol χ is the instantaneous scalar dissipation rate defined as

χ ¼ 2Dj∇Z j2 ;

ð2Þ

M.M. Ameen, J. Abraham / Chemical Engineering Science 127 (2015) 334–343

where D is the molecular diffusivity. Notice that this fundamental definition cannot be employed in LES (or RANS) because only filtered (or Reynolds-averaged) values of Z are available. In a mixing layer, the functional form of the dependence of χ on Z is typically assumed to follow an error function profile (Peters, 2000; Mukhopadhyay and Abraham, 2012). n
2 o exp 2 erf c  1 ð2Z Þ ð3Þ χ ¼ χ st n
 2 o: exp 2 erf c  1 ð2Z st Þ By using this assumption, the value of the scalar dissipation rate at any Z can be related to its value at the stoichiometric mixture fraction, χst. The unsteady flamelet equations (Eq. (1)) can be solved for different values of χst and the solution is usually tabulated as a function of 3 independent variables Z, χst and Λ, where Λ is called the progress variable and it is an indicator of how much the reactions have progressed in the flamelet (Pierce and Moin, 2004; Ihme et al., 2005; Ihme and Pitsch, 2008; Ihme and See, 2010). These models are only valid if a single flamelet is present in an LES computational cell. Typically, every computational cell contains multiple flamelets, and the filtered reaction rate in a computational cell is given by averaging the reaction rate over all these flamelets. Z Z Z Z Z Z       g _ϕ ¼ ω ω_ ϕ P Z; χ st ; Λ dZdχ st dΛ ¼ ω_ ϕ P ðZ ÞP χ st P Λ ð4Þ   In the above equation, P Z; χ st ; Λ is the joint-PDF of Z, χst and Λ, and P(Z), P(χst) and P(Λ) are the marginal PDFs of Z, χst and Λ respectively. The implicit assumption is made that these 3 variables are statistically independent so that the joint PDF can be written as the product of the marginal PDFs. Mukhopadhyay and Abraham (2012) have verified this assumption for reacting mixing layers. In this study, models for the filtered scalar dissipation rate, χf st , and its marginal PDF are studied using DNS. The usual reference in the literature to DNS is to computations that resolve the flow scales, but not necessarily the flame scales. In the present study, both flow and flame scales are resolved. This makes the computations with complex kinetics computationally intensive. The complex kinetics is required because the interest is in autoigniting mixing layers of relevance to fuel-injected compression-ignition engines. Furthermore, the computations are carried out under high-pressure and high-temperature conditions relevant to engines. For this reason, the DNS simulations performed in this work are 2D. In prior 3D DNS studies, single-step and two-step kinetics, or other simplifications using tabulated chemistry have been employed (van Oijen et al., 2007). 3D DNS of reacting flows in which the flame scales are also resolved and multi-step kinetics is employed are still rare (Yoo et al., 2011, 2013). They are generally restricted to atmospheric conditions. There have been a few other studies reported in the literature, which also evaluate models for the filtered scalar dissipation rate. Girimaji and Zhou (1996) performed DNS of isotropic non-reacting turbulence and used the results to analyze and model the effect of subgrid scales on the filtered scalar dissipation rate by accounting for two physical phenomena: (i) the interaction among the subgrid scales and (ii) the interaction between the resolved and unresolved scales. Using spectral arguments, they derived a model, which, to the leading order, involves enhancing the molecular diffusivity with a turbulent diffusivity. They also included a model for “backscatter” of energy from the small scales to the large scales. Although this model performed well when compared to the DNS results, its performance for reacting simulations have not been evaluated. The implementation of this model for LES is challenging, as there is a requirement to solve additional transport equations. Jimenez et al. (2001) performed DNS of inert scalar mixing in homogenous turbulence to model the filtered

335

scalar dissipation rate using a characteristic time-scale model. Model parameters were obtained by filtering the DNS data and the model was shown to predict the decay of the scalar variance accurately. Balarac et al. (2008) performed DNS of forced isotropic turbulence and used the principle of optimal estimators to evaluate different models for filtered scalar dissipation rate. They found that a model based on the ratio of turbulent-to-scalar time scale, where the turbulent time scale is evaluated from the subfilter kinetic energy gives the lowest irreducible error. These studies have all been performed using nonreacting DNS and it is not known whether their conclusions can be extended to reacting flows. A recent study by Mukhopadhyay and Abraham (2012) has shown that combustion and the associated heat release modifies the scalar dissipation rate significantly, especially near the flame front. So, the conclusions made regarding the models using non-reacting DNS may not hold when these models are applied for reacting flows. To the authors' best knowledge, the only DNS study for modeling the filtered scalar dissipation rate in reacting flows is the one performed by Knudsen et al. (2012). They performed DNS of autoigniting jets and used the results to evaluate an algebraic model and a transport equation model for the filtered scalar dissipation rate and showed that the transport equation model performed superior to the algebraic model. They only evaluated one algebraic model – the commonly used turbulent diffusivity model. Although the transport equation model performed well, its computational overhead is greater as there is a need to solve additional transport equations. It is the purpose of the present study to evaluate other algebraic closure models for the filtered scalar dissipation rate and test their performance for a variety of conditions. In this study, we assess the performance of the scalar dissipation rate models by using data from DNS of both nonreacting and reacting flows. We also analyze how the model performance changes with changes in filter sizes and turbulence intensities. Note that the reacting flows of specific interest in this work are those generated through autoignition of fuel/air mixing layers at elevated temperature and pressure conditions with application to compression-ignited fuel-injected engines. The outline for the rest of the paper is as follows. In the next section, the details about the computational setup and model are described. The different models for the filtered scalar dissipation rate evaluated in this study are discussed in Section 3. The models are assessed using the DNS results for a variety of conditions. In Section 4, two models for the PDF of the scalar dissipation rate are described and the conditions under which each PDF works well are examined. The variance of the scalar dissipation rate is an important parameter that is required to determine the PDF. In Section 5, a simple model for the scalar dissipation rate variance is derived and it is shown to be accurate for a range of conditions. The paper closes with summary and conclusions.

2. Computational setup The numerical code employed in this work solves the compressible form of the Navier–Stokes equations for multicomponent gaseous mixtures with chemical reactions. The sixth-order compact finite-difference scheme of Lele (1992) is employed to spatially discretize the governing equations. The resulting discretized equations are solved using the tridiagonal matrix algorithm. A fourth order Runge–Kutta scheme is employed to perform the timeintegration. Boundary conditions are implemented using the Navier–Stokes characteristic boundary conditions method of Poinsot and Lele (1992), which is extended to account for multicomponent transport (Anders et al., 2007). Chemical kinetic source terms are computed through an interface with CHEMKIN-like subroutines. The effective binary diffusion coefficient model for

336

M.M. Ameen, J. Abraham / Chemical Engineering Science 127 (2015) 334–343

computing species diffusion, using the method of Bird et al. (1960), is implemented. The code is written in FORTRAN 90 and parallelized using the message passing interface library. The accuracy of the code has been assessed in several prior studies (Anders et al., 2007; Venugopal and Abraham, 2008; Mukhopadhyay and Abraham, 2011). The method that is employed to generate turbulence is the one developed by Fathali et al. (2007). This method generates a homogenous decaying turbulent flow field with a specified set of Reynolds stresses and length scales. In this method, a complex flow field is considered as the linear superposition of multiple simpler white noise fields. When combining the simpler flow fields, the weight assigned to each simple flow fields is determined by constraints imposed by the prescribed set of Reynolds stresses and length scales. The generated turbulent flow field is employed as an initial condition for DNS simulations. In this work, we select the Reynolds stresses to be equal in the two directions. The numerical grid is smaller than the Kolmogorov scale. In prior work, it is reported that the turbulent integral length scale, l0, should be less than about 0.3 of the length of the domain to avoid the influence of boundary conditions. In our simulations, the largest turbulent length scale l0 is selected to be 500 μm, which is 0.1 times the length of the domain in either x or y direction. Initial turbulence intensities u0 range from 0.5 to 2.5 m/s. The turbulent Reynolds number Ret defined as Ret ¼ u0 l0/ν, ranges from 50 to 250 in our simulations. The Kolmogorov scale lk which is estimated as lk ¼ l0 Ret  3=4 varies from 8 to 27 μm. Therefore, numerical grid resolution is selected such that the smallest length scale is resolved in our turbulent simulations. Initialization of the turbulent velocity is done as described above. We employ periodic boundary condition and adiabatic slip boundary condition in x and y directions, respectively. N-heptane is used as the fuel in our simulations. Fig. 1 shows the initial mixture fraction field in the computational domain. Table 1 shows the list of simulation parameters, which are used to generate the DNS database. Parametric variations of the mixing layer thickness and the turbulence intensity are also performed. The case with u0 ¼1.0 m/s and δ ¼120 μm is considered as the baseline case in the discussion. Both non-reacting and reacting simulations are performed for each case. The 37-species mechanism developed by Peters et al. (2002) is used as the chemical kinetic mechanism. This mechanism has been employed in several

Table 1 List of simulation parameters employed in this study. u0 (m/s)

l0 (μm)

δ (μm)

0.5 1.0 2.5

500 500 500

90, 120, 240, 480 90, 120, 240, 480 90, 120, 240

prior studies to study autoignition and flame development in engines (Bajaj et al., 2013). 3. Modeling the filtered scalar dissipation rate Fig. 2 shows the instantaneous profile of the scalar dissipation rate at a time of 0.7 ms for the (a) non-reacting and (b) reacting mixing layers. The regions of high χ are localized to small regions in the computational domain. It is seen that the peak values of scalar dissipation rate are approximately two times higher for the non-reacting mixing layer. This is expected as chemical reactions cause an increase in temperature leading to local expansion and thus reduction in the gradients of the mixture fraction. This directly corresponds to a reduction in the scalar dissipation rate. This is clearer in Fig. 3 which compares the conditionally-averaged scalar dissipation rates, 〈χ j Z〉, for the non-reacting and reacting mixing layers at the same time. It can be seen that 〈χ j Z〉 is higher for the non-reacting mixing layer for all value of Z. For the reacting mixing layer, it is also important to look at the evolution of χ during the ignition process. Fig. 4 shows contour plots of χ at 0.20 and 0.30 ms after start of computation. The heat release due to chemical reactions cause the scalar dissipation rate profiles to change from a bimodal distribution at 0.2 ms to an approximately unimodal distribution at 0.3 ms. This is made clearer in Fig. 5 which shows the conditionally-averaged scalar dissipation rate profiles as a function of mixture fraction at 0.2 and 0.3 ms. The conditionally-averaged temperature profiles are also shown. Consistent with Fig. 4(a), the scalar dissipation rate shows two peaks at 0.2 ms. One peak at Z ¼0.5 corresponds to the original peak in the χ profile. The other peak near Z¼0.06 corresponds to the location of peak temperature. Recall from Eq. (2) that χ depends on molecular diffusivity and mixture fraction gradient |∇Z|. When the temperature rises, D rises along with it. The heat release also causes the reacting mixture to expand which leads to a reduction in |∇Z|. This, in turn, leads to a reduction in χ. This expansion due to heat release is, however, not instantaneous and thus the reduction in χ due to heat release is only felt after a finite time. For this reason, although the temperature does not change noticeably between 0.2 and 0.3 ms, the local peak in the 〈χ|Z〉 profile at ZE 0.06 decreases from 60 s  1 to 10 s  1. As expected, the χ distribution will also depend on the value of the turbulence intensity u0 . Fig. 6 shows the distribution of χ at a time of 0.5 ms for three values of u0 , 0.5, 1.0 and 2.5 m/s, when the mixture is reacting. Not surprisingly, as u0 increases, the mixing layer becomes increasingly stretched, leading to larger values of χ. The quantitative increase in χ with increase in turbulence intensity is made clearer in Fig. 7, which shows the conditionally-averaged scalar dissipation rate profiles as a function of Z for these three cases. At this point, the DNS database will be employed to assess the accuracy of scalar dissipation rate sub-models. Starting with Eq. (2), the filtered scalar dissipation rate is defined as e Ze j2 þ χ χe ¼ 2Dg j∇Z j2 ¼ 2Dj∇ model ;

Fig. 1. The initial mixture fraction field in the computational domain, when δ ¼120 μm.

ð5Þ

where χ model is the term to be modeled, and the tilde denotes filtered variables. In this study, several models proposed in the literature for

M.M. Ameen, J. Abraham / Chemical Engineering Science 127 (2015) 334–343

337

Fig. 2. Instantaneous scalar dissipation rate contours at 0.7 ms for (a) non-reacting mixing layer and (b) reacting mixing layer.

simulations when the k–ε model is employed for turbulence, a model can be formulated as proposed by Jimenez et al. (2001) where

εe χ model ¼ χ kε ¼ C kε Z v ;

ð8Þ

e k

 e 1=2ug ee where i ui   k ¼g  ui ui is the sub-grid scale kinetic energy, εe ¼ ν ∂ui =∂xj U ∂ui =∂xj is the filtered kinetic-energy dissipation, ν is the viscosity, Zv is the variance of mixture fraction and Ckε is a model parameter to be determined. ej , The turbulent time-scale can also be defined as τ ¼ 1=j S where j e S j is the magnitude of the large-scale strain rate tensor,  1=2 e f f j S j ¼ 2S with ij Sij  e e ∂ u 1 ∂ u j i f þ : ð9Þ S ij ¼ 2 ∂xj ∂xi

The strain-rate tensor (SRT) model (Balarac et al., 2008) then is defined as Fig. 3. Variation of the conditionally-averaged scalar dissipation rate with Z for the non-reacting and reacting baseline cases at 0.7 ms.

the filtered scalar dissipation rate are evaluated by comparing their performance against DNS results. The models are described next. The first model that is evaluated is the turbulent diffusivity (TD) model (Girimaji and Zhou, 1996; Pierce and Moin, 2004), given by

χ model ¼ χ TD ¼ 2DT j∇Ze j2 ;

ð6Þ

where DT is the sub-grid scale turbulent diffusivity in LES. In RANS models, the scalar dissipation rate is usually linked algebraically to the variance of the mixture fraction, Zv (Sanders and Gokalp, 1998). A similar approach can be employed for LES as well. Here, the sub-filter scalar mixing time, Zv/χ model , is assumed to be proportional to the sub-filter turbulent timescale, τ, i.e. Z

χ model ¼ C v : τ

ð7Þ

Based on the choice for the definition of the time-scale τ, different models for the scalar dissipation rate can be derived. Following a method that is commonly employed in RANS

χ model ¼ χ SRT ¼ C SRT Z v j Se j ;

ð10Þ

where C SRT is a model parameter to be determined. 1=2 The turbulent time scale can also be defined as τ ¼ k=Δ , where Δ is the filter size (which is assumed to be equal to the LES grid size). This leads to the subfilter kinetic energy (SKE) model (Schmidt and Schumann, 1989; Balarac et al., 2008) given by

χ model ¼ χ SKE ¼ C SKE

Zvk

1=2

Δ

;

ð11Þ

where again CSKE is a model parameter to be determined. Table 2 lists the different models for the filtered scalar dissipation rate that are evaluated in this study. The DNS database is used to obtain the filtered scalar dissipation rate by explicitly filtering the DNS results using a box filter with different filter sizes ranging from 50 to 500 μm. The TD model, k–ε model, SRT model and SKE model are then assessed in detail by comparing their accuracy with that of the DNS values for nonreacting and reacting mixing layers. Among these models, the TD and SRT models are the easiest to implement in LES computations, as there is no need to solve additional transport equations for k and ε. Note that with the exception of the TD model, an equation for the mixture fraction variance has to be solved. To compare the performance of the models, a normalized error E is defined similar

338

M.M. Ameen, J. Abraham / Chemical Engineering Science 127 (2015) 334–343

Fig. 4. Instantaneous scalar dissipation rate contours for the reacting baseline case during the ignition process at (a) t¼ 0.20 ms and (b) t¼ 0.30 ms.

Fig. 5. Conditionally-averaged temperature and scalar dissipation rate profiles as a function of Z at (a) t¼ 0.20 ms and (b) t ¼0.30 ms.

to the one used by Balarac et al. (2008) and given by

  2  e j ∇Ze j 2 χ model  χe  2D

 ; E¼ 2  e j ∇Ze j 2 χe 2D

ð12Þ

where the angular brackets denote an ensemble averaging over all the e j ∇Z e j 2 Þ2 is the actual quantity e  2D filtered cells. Note that the term ðχ

estimated from the DNS database. This error norm is a measure of the lack of correlation between the model and the DNS results. A value of E¼0 implies perfect correlation and as the correlation reduces, E increases. The value of the error norm, of course, depends on the model constants. The model constants are selected by minimizing the

error norm for a baseline case (u0 ¼1.0 m/s and δ ¼120 μm) over a range of filter sizes and then averaging it across the filter sizes. For the non-reacting case, the model parameters were found to be Ckε ¼0.7, CSRT ¼0.1 and CSKE ¼0.57. Fig. 8 compares the value of E for the baseline non-reacting case for filter sizes (Δ) varying from 50 to 500 μm. Fig. 8 shows that the error generally increases with increasing Δ for all the

models. The error is negligible when Δ ¼ 50 μm, but increases with increasing Δ. The SKE model performs the best among the 4 models across the range of filter sizes. The performance of the TD model is relatively poor for all the filter sizes. The normalized error for the SKE model is about 0.15 for Δ ¼500 μm but about 0.6 for the TD model. In fact, the error norm for the TD model has values that are not very different from using no model. Examining Eqs. (5) and (6), this suggests that the contribution of the sub-grid scale term (Eq. (6)) is relatively small in the TD model. This also suggests that the turbulent diffusivity DT, which is obtained from the Smagorinsky model in this study, is underpredicted. Notice that the difference between the two curves (no model vs TD) increases with increasing filter size suggesting that the sub-grid scale contribution does increase as expected. This conclusion is similar to the one made by Balarac et al. (2008). It is also seen that the SRT model performs relatively well under all conditions. Fig. 9 makes the same comparison for the reacting baseline case. It is seen that the TD model again performs poorly, showing negligible improvement over using no model. The TD model is used in many LES computations (Ihme et al., 2005; Ihme and Pitsch, 2008; Ihme and See, 2010). When flamelet models are employed (see Eq. (4)) together with

M.M. Ameen, J. Abraham / Chemical Engineering Science 127 (2015) 334–343

339

Fig. 6. Effect of turbulence intensity on the instantaneous scalar dissipation rate contours for the turbulent reacting mixing layers at t ¼ 0.5 ms for (a) u0 ¼0.5 m/s, (b) u0 ¼1.0 m/s and (c) u0 ¼2.5 m/s.

the TD model, these results suggest that there can be large errors. The SRT and SKE models have the minimum errors across the range of filter widths. Recall that the SKE model had the minimum error for the nonreacting case. For the reacting case, the model parameters were found to be Ckε ¼0.11, CSRT ¼0.11 and CSKE ¼0.45. While the model parameters for the SRT and SKE model for the reacting cases are within 10% and 20%, respectively, of the values obtained for the non-reacting simulations, the optimum model parameter for the k–ε model is a factor of about seven lower. The reason for this is that in the reacting mixing layer, the increase in temperature leads to higher values of kinematic viscosity ν and thus to higher values for ε. Since the filtered scalar dissipation rate predicted by the k–ε model is directly proportional to ε (see Eq. (8)), Ckε has to be lowered to adjust for the rise in temperature. Note that a temperature-dependent property of the fluid, e.g. ν, appears only in the k–ε model. From this point on, only the SRT and SKE model will be considered for further assessment. The constants derived for the reacting cases are now used to perform further parametric studies by varying mixing layer thickness and turbulence intensity. Figs. 10 and 11 show the effect of the mixing layer thickness on the performance of the SRT and SKE models for the non-reacting

and reacting simulations, respectively, for a filter size of 200 μm using the constants determined from the baseline case. Also shown for comparison are the error values obtained without using any model. The general trend from Fig. 10 is that as the mixing layer thickness increases, the error E reduces. This is expected, as the increase in mixing layer thickness leads to reduction in the gradients in Z and thus lowers the values of the scalar dissipation rate. A similar trend is observed for Fig. 11. It can be concluded that both the SRT and SKE models perform equally well for a range of mixing layer thicknesses. Figs. 12 and 13 show the effect of the turbulence intensity on the performance of the models for non-reacting and reacting mixing layers, respectively. As the turbulence intensity u0 increases, the mixing layer becomes increasingly stretched, leading to larger values of χ. The error increases as the turbulence intensity increases with and without models; but, the SRT and SKE models reduce the error to about 50% of that when using no model. It can be concluded that the SRT and SKE models are able to model the filtered scalar dissipation rate relatively well for the wide range of conditions selected in this study. Use of the SKE model, however, requires that a transport equation for the subgrid turbulent kinetic

340

M.M. Ameen, J. Abraham / Chemical Engineering Science 127 (2015) 334–343

Fig. 9. Comparison of the different models for filtered scalar dissipation rate for the baseline reacting case.

Fig. 7. Conditionally-averaged scalar dissipation rate profiles as a function of Z at 0.5 ms for different turbulence intensities.

Table 2 List of models for filtered scalar dissipation rate. Model

χ model

1

TD model

2 χ model ¼ χ TD ¼ 2DT ∇Ze

2

k–ε model

3

SRT model

χ model ¼ χ SRT ¼ C SRT Z v j Se j

4

SKE model

χ model ¼ χ SKE ¼ C SKE Z v kΔ

χ model ¼ χ kε ¼ C kεeε Z v e k

1=2

Fig. 10. Effect of mixing layer thickness (δ) on the performance of the SRT and SKE models for non-reacting mixing layers.

Fig. 8. Comparison of the different models for filtered scalar dissipation rate for the baseline non-reacting case (u0 ¼ 1.0 m/s and δ ¼ 120 μm).

energy has to be solved. Since the differences in performances between these two models are not very significant, the SRT model is the recommended model to be used for LES computations.

4. Modeling the PDF of filtered scalar dissipation rate As shown in Eq. (4), in addition to the accurate modeling of the e , the marginal PDF of χ is also filtered scalar dissipation rate χ required for use in flamelet models. Figs. 14 and 15 show the PDFs obtained from the DNS simulations for the baseline reacting case.

Fig. 11. Effect of mixing layer thickness (δ) on the performance of the SRT and SKE models for reacting mixing layers.

When the filter size Δ is small, the entire filtered cell can be expected to have an almost uniform value of χ. In other words, the variance of χ can be expected to be very small. Under these conditions, the PDF of χ is expected to behave like a δ-function e . Fig. 14 shows that for Δ ¼100 μm, the shape of centered on χ ¼ χ the PDF of χ is similar to a δ-function. It can be expected that when Δ is large, the regions with high χ are localized to a small fraction of the filter size, and the remaining fraction has lower χ. Hence,

M.M. Ameen, J. Abraham / Chemical Engineering Science 127 (2015) 334–343

341

Fig. 12. Effect of turbulence intensity on the performance of SRT and SKE model for non-reacting mixing layers. Filter size ¼200 μm.

Fig. 14. Comparison of marginal PDFs for the filtered scalar dissipation rate with the actual PDF for a filter size of 100 μm: (a) χe ¼ 1.1 s  1, and (b) χe ¼ 76.3 s  1.

Fig. 13. Effect of turbulence intensity on the performance of SRT and SKE model for reacting mixing layers. Filter size ¼ 200 μm.

e and a larger the PDF of χ is expected to show a small peak near χ peak at a value of χ ¼0. Fig. 15 shows that this is indeed the case for Δ ¼500 μm. Based on the general shapes of these PDFs, two functional forms are used to model the PDF – the exponential PDF and the lognormal PDF. The exponential PDF is given by    1 χ PDFexp χ ¼ exp  : ð13Þ χe χe The lognormal PDF is given by   PDFLN χ ¼

1

χχ 1=2 var ð2π Þ

exp  1=2



e ln χ  χ 2χ var

1=2 !

;

ð14Þ

where χ var is the variance of the scalar dissipation rate. Fig. 14 compares the actual PDF with the exponential and lognormal PDF for a filter size of 100 μm. It is seen that the exponential PDF e whereas the lognormal PDF performs well for low values of χ e . The reason for the poor performs well for larger values of χ e is that the performance of the lognormal PDF at low values of χ lognormal PDF always has a value of 0 at χ ¼ 0, whereas the discussion above showed that the actual PDF tends to have a none . Fig. 15 makes the same zero value at χ ¼0 for low values of χ comparison for a filter size of 500 μm and the results show that the exponential PDF does a better job at this filter size because the

e . In other larger value of the filter size results in a lower value of χ words, the localized region where χ is high is compensated by a larger region with lower χ. An error norm EPDF can be defined to quantify the performance of the model PDFs as EPDF ¼

Z

ðP ðZ Þ  P model ðZ ÞÞ2 dZ;

ð15Þ

where P(Z) is the actual PDF from the DNS database and Pmodel(Z) is the model PDF. The error norm was calculated for the e in this set of cases shown in Table 1. The range of values of χ simulations is from 5 to 50 s  1. A normalized scalar dissipation rate is defined to be

χe χ norm ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; u02 =Δδ

ð16Þ

where u0 is the turbulence intensity, Δ is the filter size, and δ is the mixing layer thickness. Because the performance of the model e , insight can be gained by plotting PDF depends on the value of χ the error norm as a function of χnorm. Fig. 16 compares the error norms for the exponential and lognormal PDFs. It is seen that the exponential PDF works better when χnorm is lower than 0.6, and the lognormal PDF works better for larger values. It can be e , and large values of concluded that for low values of Δ, δ and χ u0 , the exponential PDF is the better choice.

342

M.M. Ameen, J. Abraham / Chemical Engineering Science 127 (2015) 334–343

Fig. 16. Comparison of errors for exponential and lognormal PDFs. The vertical line shows the value of χnorm above which the lognormal PDF is more accurate.

where g e 2 ð∇Ze U ∇Z″Þ2  16D e 2 j∇Ze j ∇Z sd 2 ; χ var;1 ¼ 16D 2

g

4

χ var;2 ¼ 4D″2 j∇Ze j ;

ð19Þ ð20Þ

and g

2

e χ var;3 ¼ 16DD''ð∇ Ze :∇Z″Þ j∇Ze j :

ð21Þ

In Eq. (19), Zsd is the standard deviation of the mixture fraction, which is the square root of the mixture fraction variance. In Eq. (19), the approximation is made that the gradient of the mean mixture fraction and the gradient of the standard deviation of the mixture fraction are parallel vectors and hence their dot product can the product of their magnitudes. By examining Eqs. (19)–(21), it can be seen that χvar,1 is more dominant as compared to χvar,2 and χvar,3. So, as a first approximation, Eq. (18) can be rewritten as

χ var;model ¼ K χ var;1 ; Fig. 15. Comparison of marginal PDFs for the filtered scalar dissipation rate with the actual PDF for a filter size of 500 μm: (a) χe ¼2.5 s  1, and (b) χe ¼ 20.1 s  1.

5. Modeling scalar dissipation rate variance When using the SRT or SKE model, a transport equation for the variance of mixture fraction, i.e. Z v , has to be solved. Knowing Z v , the variance of the scalar dissipation rate can be modeled as shown below. The variance of the scalar dissipation rate is needed to employ the lognormal PDF. In this section, a simple model is derived for the variance, which is found to be applicable for all the cases considered here. Using Eq. (2), the variance of the scalar dissipation rate can be expressed as g

e Þ2 χ var ¼ χg ''2 ¼ ðχ  χ

ð17Þ

Eq. (17) can be expanded to express χvar as a sum of three terms as follows:

χ var ¼ χ var;1 þ χ var;2 þ χ var;3 ;

where K is a model parameter. In the above expression, the variance of the scalar dissipation rate can be obtained solely based on the values of the filtered mixture fraction, the variance of the mixture fraction and the filtered diffusivity. Figs. 17 and 18 confirm the validity of Eq. (22) for filter sizes of 100 μm and 200 μm respectively. The model parameter K can be obtained as the slope of the best-fit line. K is found to be in the range of 2–2.3 for all the cases considered in this study. It is worth mentioning that although the slope of K obtained from the linear fit is along expected lines, there are significant departures from the linear correlation especially for large filter sizes as Fig. 18 shows. These departures are attributed to the approximations made in Eqs. (19) and (22).

6. Summary and conclusions

g e þ D″Þ ∇ðZ þ Z''Þ 2  2ðD e þ D″Þg ∇ðZ þ Z″Þ 2 Þ ¼ ð2ðD g e Ze d∇Z″Þ þ2D''j∇Ze j2 Þ2 : ¼ ð4Dð∇

ð22Þ

ð18Þ

In this study, direct numerical simulations of non-reacting and reacting mixing layers have been carried out to generate databases, which are then employed to assess the accuracy of four models for the filtered scalar dissipation rate in LES. The pressure and temperature conditions selected are relevant to compression-ignition combustion engines. N-heptane, often used as a surrogate for diesel fuel, is used as the fuel. The four models assessed are the turbulent diffusivity model, the k–ε model, the strain rate tensor (SRT) model and the subfilter

M.M. Ameen, J. Abraham / Chemical Engineering Science 127 (2015) 334–343

343

variance and PDF of the scalar dissipation rate are expected to improve predictions when employed as subgrid-scale models in LES. Future studies will evaluate these models through LES of reacting jets.

Acknowledgments The authors would like to thank the US National Institute of Computational Sciences (NICS), eResearch SA (eRSA) and the Australian National Computational Infrastructure (NCI) for providing the computing resources for this work, and Caterpillar Inc. (grant number 204421) and Purdue Research Foundation (grant number 204533) for financial support of this work. The contributions of Professor Vinicio Magi in the development of the numerical codes are gratefully acknowledged. References

Fig. 17. Determining the validity of χvar,model (refer Eq. (22)) for filter size of 100 μm. The slope of the best-fit line gives the value of K to be 2.17.

Fig. 18. Determining the validity of χvar,model (refer Eq. (22)) for filter size of 200 μm. The slope of the best-fit line gives the value of K to be 2.28.

kinetic energy (SKE) model. An error norm, E, is employed to quantify the differences between the model results and the DNS results. The assessment is carried out for a range of initial mixing layer thicknesses (90–480 μm), turbulence intensities (0.5–2.5 m/s), and filter widths (50–500 μm). Based on the values of the error norm, it is found that the SRT and SKE models perform the best among the models considered for the range of conditions considered. The SRT model is recommended because, unlike the SKE model, it can be used without solving additional transport equations (for the subfilter turbulent kinetic energy). To employ the model, the probability density function (PDF) of the scalar dissipation rate is required. It is shown that the choice of PDF depended on the value of the filtered scalar dissipation rate – for lower values, an exponential PDF performs well, whereas a lognormal PDF performs better for larger values. A model was also introduced that relates the variance of the scalar dissipation rate to the mean and variance of the mixture fraction. This model is shown to give satisfactory agreement with the DNS results. These improved models for the mean,

Anders, J.W., Magi, V., Abraham, J., 2007. Large-eddy simulation in the near-field of a transient multi-component gas jet with density gradients. Comput. Fluids 36, 1609–1620. Bajaj, C., Ameen, M., Abraham, J., 2013. Evaluation of an unsteady flamelet progress variable model for autoignition and flame lift-off in diesel jets. Combust. Sci. Technol. 185, 454–472. Balarac, G., Raman, V., Pitsch, H., 2008. Modeling of the subfilter scalar dissipation rate using the concept of optimal estimators. Phys. Fluids 20, 091701. Bird, R.B., Stewart, W.E., Lightfoot, E.N., 1960. Transport Phenomena. Wiley, New York. Fathali, M., Klein, M., Broeckhoven, T., Lacor, C., Baelmans, M., 2007. Generation of turbulent inflow and initial conditions based on multi-correlated random fields. Int. J. Numer. Methods Fluids 57, 93–117. Girimaji, S.S., Zhou, Y., 1996. Analysis and modeling of subgrid scalar mixing using numerical data. Phys. Fluids 8, 1224. Ihme, M., Cha, C.M., Pitsch, H., 2005. Prediction of local extinction and re-ignition effects in non-premixed turbulent combustion using a flamelet/progress variable approach. Proc. Combust. Inst. 30, 793–800. Ihme, M., Pitsch, H., 2008. Prediction of extinction and re-ignition in nonpremixed turbulent flames using a flamelet/progress variable model 1. A priori study and presumed PDF closure. Combust. Flame 155, 70–89. Ihme, M., See, Y.C., 2010. Prediction of autoignition in a lifted methane/air flame using an unsteady flamelet/progress variable model. Combust. Flame 157, 1850–1862. Jimenez, C., Ducros, F., Cuenot, B., Bedat, B., 2001. Subgrid scale variance and dissipation of a scalar field in large eddy simulations. Phys. Fluids 13, 1748. Knudsen, E., Richardson, E.S., Doran, E.M., Pitsch, H., Chen, J.H., 2012. Modeling scalar dissipation and scalar variance in large eddy simulation: algebraic and transport equation closures. Phys. Fluids 24, 055103. Lele, S.K., 1992. Compact finite difference schemes with spectral-like resolution. J. Comput. Phys. 103, 16–42. Mukhopadhyay, S., Abraham, J., 2011. Influence of compositional and thermal stratification on autoignition in n-heptane/air mixtures. Combust. Flame 158, 1064–1075. Mukhopadhyay, S., Abraham, J., 2012. Evaluation of an unsteady flamelet progress variable model for autoignition and flame development in compositionally stratified mixtures. Phys. Fluids 24, 075115. Peters, N., 2000. Turbulent Combustion. Cambridge University Press; Cambridge. Peters, N., Paczko, G., Seiser, R., Seshadri, K., 2002. Temperature cross-over and nonthermal runaway at two-stage ignition of n-heptane. Combust. Flame 128, 38–59. Pierce, C.D., Moin, P., 2004. Progress-variable approach for large-eddy simulation of non-premixed turbulent combustion. J. Fluid Mech. 504, 73–97. Pitsch, H., Chen, M., Peters, N., 1998. Unsteady flamelet modeling of turbulent hydrogen–air diffusion flames. Symp. (Int.) Combust. 27, 1057. Poinsot, T.J., Lele, S.K., 1992. Boundary conditions for direct simulations of compressible viscous flows. J. Comput. Phys. 103, 104–129. Sanders, J.P.H., Gokalp, I., 1998. Scalar dissipation rate modelling in variable density turbulent axisymmetric jets and diffusion flames. Phys. Fluids 10, 938. Schmidt, H., Schumann, U., 1989. Coherent structure of the convective boundary layer derived from large-eddy simulations. J. Fluid Mech. 200, 511–562. van Oijen, J.A., Bastiaans, R.J.M., de Goey, L.P.H., 2007. Low-dimensional manifolds in direct numerical simulations of premixed turbulent flames. Proc. Combust. Inst. 31, 1377–1384. Venugopal, R., Abraham, J., 2008. A 2-D DNS investigation of extinction and reignition dynamics in nonpremixed flame vortex interactions. Combust. Flame 153, 442–464. Yoo, C.S., Richardson, E.S., Sankaran, R., Chen, J.H., 2011. A DNS study on the stabilization mechanism of a turbulent lifted ethylene jet flame in highlyheated coflow. Proc. Combust. Inst. 33, 1619–1627. Yoo, C.S., Luo, Z., Lu, T., Kim, H., Chen, J.H., 2013. A DNS study of ignition characteristics of a lean iso-octane/air mixture under HCCI and SACI conditions. Proc. Combust. Inst. 34, 2985–2993.