- Email: [email protected]
Deterministic Multiple Mapping Conditioning (MMC) applied to a turbulent flame in Large Eddy Simulation (LES)
Available online at www.sciencedirect.com
Proceedings of the
Combustion Institute
Proceedings of the Combustion Institute 34 (2013) 1213–1221
www.elsevier.com/locate/proci
Deterministic Multiple Mapping Conditioning (MMC) applied to a turbulent flame in Large Eddy Simulation (LES) C.B. Devaud a,⇑, I. Stankovic´ b, B. Merci b a
Department of Mechanical and Mechatronics Engineering, University of Waterloo, 200 University Avenue West, Waterloo, ON, Canada N2L3G1 b Ghent University, Department of Flow, Heat and Combustion Mechanics, Belgium Available online 2 July 2012
Abstract The scalar dissipation rate is a key quantity in turbulent combustion modelling, in particular for Conditional Moment Closure (CMC). In Large Eddy Simulation (LES), a subgrid-scale (SGS) mixing model must be provided to determine the SGS mixture fraction variance needed for the filtered scalar dissipation rate and the filtered density function. Multiple Mapping Conditioning (MMC) is selected as a SGS mixing model and the deterministic form is implemented for the first time in LES. The objective of the present study is to assess the capabilities of MMC to model the turbulent mixing field and the conditionally filtered g needed in CMC. scalar dissipation rate, Njg The MMC transport equation is implemented in a LES code coupled with CMC to simulate a lifted jet flame in a vitiated coflow. The MMC formulation and computational details are presented. MMC predictions are compared with LES for the turbulent mixing field using instantaneous and time averaged data, and experimental values are included whenever available. Reasonable agreement between MMC, LES and experimental data is found. Best fit is obtained close to the fuel exit and some discrepancies can be seen further downstream. The SGS variance is underpredicted for many parts of the jet. A thorough analysis is performed to assess the modelling and numerical assumptions. The two proposed sources of discrepancy are the linear model and neglect of a transverse component g are significantly different from those of AMC and are for convection in MMC. The MMC profiles of Njg sensitive to the SGS variance. Future work is suggested. Ó 2012 The Combustion Institute. Published by Elsevier Inc. All rights reserved. Keywords: Multiple Mapping Conditioning; LES; Scalar dissipation rate; Conditional Moment Closure; Mixing
1. Introduction The accurate representation of the complex interactions between chemistry and turbulence is ⇑ Corresponding author. Fax: +1 519 885 5862.
E-mail Devaud).
address:
[email protected]
(C.B.
crucial for the development of industrial combustion devices. Turbulent mixing is of particular importance in non-premixed combustion where chemical reactions will only occur if fuel and oxidizer are mixed at the molecular level. The current computational tools for engineering applications rely on Reynolds Averaged Navier–Stokes (RANS) and Large Eddy Simulations (LES)
1540-7489/$ - see front matter Ó 2012 The Combustion Institute. Published by Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.proci.2012.06.076
1214
C.B. Devaud et al. / Proceedings of the Combustion Institute 34 (2013) 1213–1221
approaches where, in both cases, the turbulent small scales need to be modelled. The scalar dissipation rate, N, represents the rate of mixing at the molecular level and is proportional to the square gradient of the scalar, Z, such as N D$Z $Z, where D is the molecular diffusivity of Z. Thus, this is a key quantity in turbulent combustion modelling, in particular for flamelet, Conditional Moment Closure (CMC) and Probability Density Function (PDF) approaches [1]. Within the CMC and turbulent non-premixed combustion framework, the conditional average of N at a particular value, g, of mixture fraction, Z is of special interest with hNjgi = hD$Z $ZjZ = gi. The angular brackets denote a conditional average over an ensemble of realizations of the flow, subject to the condition to the right of the vertical bar. Evaluation of hNjgi is not straightforward and CMC further requires solution of the transport equations for the conditional averages to be consistent with that of the PDF transport equation [2]. Two of the most commonly used models in CMC are the presumed b-PDF model of Girimaji [3] and Amplitude Mapping Closure (AMC) [4]. Both models are derived assuming homogeneous turbulence and consequently, do not comply with the solution of the PDF transport equation in inhomogeneous turbulent flow. Other formulations for hNjgi have been proposed and can be divided into two groups. The first approach is based on integrating the PDF transport equation to enforce the required consistency in CMC [5–8]. This shows some success but suffers from potentially large numerical errors and inaccuracies due to the conditional velocity model [5]. The second strategy relies on Multiple Mapping Conditioning (MMC), initially proposed by Klimenko and Pope [9], in which hNjgi appears in closed form without homogeneous turbulence assumption or complex integration. Vogiatzaki et al. [10,11] implemented MMC with a single major scalar, mixture fraction, and calculated the unconditional and conditional scalar dissipation rate for simple jet flames. These results show good agreement with experimental data. All these past investigations were conducted in the context of RANS. In LES, the large scale turbulent motions are resolved but a subgrid-scale (SGS) mixing model must be provided to determine the SGS mixture fraction variance needed for the filtered scalar dissipation rate and the filtered density function (FDF) when a presumed form is considered for the FDF [12]. In LES–CMC, it is common practice to include the AMC model based on homogeneous turbulence to determine the conditionally filtered scalar g [13,14]. Alternative modelling dissipation Njg may use the conditional average of filtered values over one CMC cell [15]. Recent LES–CMC results for Sandia Flame F are found to be sensitive to g and a model constant in the modelling of Njg the SGS component of N was adjusted [16].
Encouraging results have been recently reported for Sandia Flame D using the stochastic formulation of MMC with a Lagrangian numerical scheme in LES [17]. Computational cost is minimized by using sparse Lagrangian technique. Good predictions for the filtered mixture fraction are obtained while the rms are overpredicted compared to the experimental data. Further, the results are shown to be sensitive to model constants in the particle mixing model and no density feedback is implemented. Cleary and Klimenko [18] further expanded the analysis for constant and variable density flows with density coupling. Comparisons were made between LES and MMC for mean mixture fraction and its rms. Similar to previous observations [17], good agreement is found for the mean mixture fraction but the rms is overpredicted. The reactive scalar predictions are shown to be sensitive to particle mixing length scales. The objective of the present investigation is to assess the capabilities of MMC to model the turg needed in CMC, in bulent mixing field and Njg LES. In contrast to previous LES–MMC studies, the deterministic form is implemented for the first time in LES and incorporates the approach initially developed in RANS. LES is used to determine the turbulent flow and mixing fields, CMC is included to simulate turbulent combustion and in parallel MMC is implemented to emulate the turbulent mixing field and determine the conditionally filtered scalar dissipation rate. For best g found in assessment of MMC performance, Njg MMC is not coupled back to CMC (CMC retains the AMC model). The current MMC investigation builds upon previous LES–CMC studies [14,19] applied to two experimental cases, hydrogen autoignition in a coflow of heated air [20] and the Cabra lifted hydrogen flame in a vitiated coflow [21]. The latter is kept for the current MMC analysis. Previous numerical work can also be found for the hydrogen Cabra flame using different combustion models in RANS [21–27] and LES [28–30]. Autoignition is shown to be the primary stabilization mechanism [19,22] and g has the potential of improved modelling of Njg strongly affecting the CMC predictions, in particular for the lift-off height. In the next sections, the MMC formulation and computational details are presented. The results focus on comparing MMC with LES for the turbulent mixing field using instantaneous and time averaged data, and experimental values are included whenever available. 2. Multiple Mapping Conditioning (MMC) MMC is derived from the idea that all turbulent fluctuations can be divided into major and minor groups and turbulent fluctuations of minor scalars are correlated with those of the major scalars
C.B. Devaud et al. / Proceedings of the Combustion Institute 34 (2013) 1213–1221
[9,31]. A reference variable is assigned to each major scalar. In the present work, one major scalar is selected, the conditioning variable used in CMC, namely mixture fraction. Thus, it is assumed that the fluctuations of all other quantities in MMC, such as temperature and species concentrations are correlated to those of the mixture fraction. This is the same idea behind the single conditioning CMC. The MMC equations solve for mapping functions which map between the reference space and the composition space following the same principle of mapping closure methods. The equations are derived for inhomogeneous turbulence and do not depend on the shape of the reference PDF. The reference variable is noted as n. Xz is the mapping function solved by the MMC equation and Xz maps between n and Z, mixture fraction. MMC offers additional benefits beyond the scope of the current investigation. It can be easily formulated for two or more conditioning variables [31], useful for the modelling of partially premixed turbulent combustion. The MMC equation for the spatial and temporal evolution of Xz is given by [9] @X z @X z @2X z þ UrX z þ A B ¼ 0; @t @n @n2
ð1Þ
where U, A and B correspond to velocity, drift and diffusion, respectively, and need to be closed. The present MMC formulation in LES is very similar to that in RANS, also consistent with similar form of the CMC equations in RANS and LES. Closure of the unclosed terms is obtained through consistency with the reference PDF transport equation. A spatially and temporally invariant Gaussian form with zero mean and unity variance is selected for the reference PDF, Pn, following common practice [9–11]. The coefficients can be determined in principle for any form of Pn. However, the choice of Pn has an impact on the resulting complexity for the expressions of U and A. All MMC coefficients must satisfy the transport equation for Pn. Different expressions may be obtained and depend on the assumptions made for U and B. In the present work, models suggested by Klimenko and Pope [9] and implemented in RANS by Vogiatzaki et al. [10,11] are used. A linear model is applied for U given by U = U(n, x, t) = U(0) + U(1)n, where U(0) corresponds to the filtered velocity vector solved in LES and U(1) is the velocity gradient. The expression for U(1) is determined such that the first moment of the filtered mixture fraction PDF is recovered following the procedure explained by Vogiatzaki et al. [11]: Uð1Þ hnX z i ¼ ðD þ Dt Þ$hX z i ;
ð2Þ
where D and Dt are the molecular and turbulent diffusivity, respectively. The terms with overtilde are filtered quantities and terms in angular brackets with a star are averages obtained by integrating
1215
the product of any function of n with Pn over n space. In Eq. (2), hXzi* is determined by R1 X ðnÞP z n dn. The linear model is commonly 1 used in RANS for CMC and MMC [5,10]. It is exact for jointly Gaussian turbulence and is considered a good approximation for flows within two to four standard deviations of the scalar under consideration [5]. However, its validity has not been demonstrated in LES [32]. The drift and diffusion coefficients are determined by 1 @X z @X z ð1Þ e: A ¼ Bn þ rqU ; B ¼N ð3Þ q @n @n The expression for B in Eq. (3) requires e obtained from LES and includes knowledge of N contribution of resolved and subgrid gradients: e e e ¼ m þ mt @ Z @ Z ; N ð4Þ Sc Sct @xk @xk where m is the molecular kinematic viscosity, mt the turbulent viscosity, Sc is the molecular Schmidt number and Sct the turbulent Schmidt number. g needed in CMC, does not appear explicitly Njg, in the MMC formulation. It can be retrieved through a transformation from the reference space, n, to the mixture fraction space, g, such that [10,11]: 2 g B @X z : Njg ð5Þ @n 1 z Pe ðgÞ may be obtained using Pe ðgÞ ¼ P n dX . dn 3. Computational details 3.1. Experimental case configuration The selected experimental case is a turbulent lifted jet flame of hydrogen diluted with nitrogen issuing into a wide co-flow of vitiated air [21]. The burner consists of a fuel jet nozzle (inner diameter, d = 4.57 mm) and a surrounding perforated plate disk (outer diameter 210 mm). The fuel stream consists of hydrogen diluted with nitrogen with exit velocity equal to 107 m/s at 305 K. The vitiated air consists of the products of a lean premixed hydrogen/air flame at 1045 K with a velocity of 3.5 m/s. The Reynolds number is equal to 23,600 for the fuel jet and 18,600 for the coflow. The stoichiometric mixture fraction is 0.474. 3.2. LES–CMC The Favre filtered equations for mass, momentum and mixture fraction, Z, are solved in three dimensions in Cartesian coordinates in the current LES code [14,19]. Convective transport is discretized with a second order central scheme for the momentum equations and a second order
1216
C.B. Devaud et al. / Proceedings of the Combustion Institute 34 (2013) 1213–1221
TVD (Total Variation Diminishing) [33] for the mixture fraction equation. Subgrid scale stresses are modelled using a standard Smagorinsky model with a model constant, C, equal to 0.1. Equal diffusivities, constant Schmidt number (Sc = Sct = 0.7) and unity Lewis number are applied. The SGS mixture fraction variance is obtained by a gradient type model, 2 @e g 002 Z @e Z Z ¼ CD @xk @xk , where D is the mesh dependent filter width. The expression for SGS variance assumes equilibrium between production and dissipation. Other approaches may be considered [32] but are beyond the scope of the present investigation. The computational domain extends 30d downstream from the jet inlet (approximately 137 mm) in the axial direction and 20d (91.4 mm) radially. The corresponding grid contains 192 48 48 cells. The grid is stretched smoothly towards the co-flow in the radial direction and is expanded smoothly in the axial direction. A digital filter is used to generate turbulence in the coflow with a turbulence lengthscale equal to 1.6 mm, the size of the holes in the outer disk. The turbulence intensity is set to 5% following experimental information [21] and previous simulations for the same flame [25]. At the inlet of the domain, Dirichlet boundary conditions are used for the velocity and mixture fraction. At the outlet, Neumann boundary conditions are applied for all quantities except for pressure which is imposed (atmospheric pressure). The implementation is in parallel with 4 blocks in the axial direction. Turbulent combustion is modelled using a first-order CMC approach with detailed mechanism including 9 species and 19 reversible reactions [34]. The equations for the conditionally filtered species mass fractions and temperature are solved. Integration with the density-weighted filtered b PDF, Pe ðgÞ over mixture fraction space yields the unconditional density, temperature and species mass fractions. The LES and CMC solvers are coupled and the information from CMC is transferred at every timestep of the simulation to the LES in order to update flow field. A coarser spatial mesh is used to solve the CMC equations due to much weaker spatial variations of conditional averages compared to those of Favre filtered scalars. The CMC equations are discretized in finite differences. The CMC grid consists of 80 5 5 cells and the mixture fraction space is discretized into 50 bins, clustered around the most reactive mixture fraction (equal to 0.053). Information is transferred from the finer LES spatial mesh to the coarser CMC grid using volume averaging [14,32]. The AMC model is applied directly to the CMC cells in order to g At the inlet of the CMC domain, inert obtain Njg. distributions for conditional averages are set. Zero gradient boundary condition is imposed at
the sides and the downstream edge of the domain in physical space. In g space, a Dirichlet boundary condition is used at g = 0 and g = 1. 3.3. MMC numerical implementation Due to increased memory allocation, capabilities of our current computational system and significant increase in computational time, the MMC equation, Eq. (1) is solved in two spatial dimensions: the axial (y axis) and transverse (z axis) directions are kept. Transport in the third (x, transverse) direction is ignored. The centre plane is selected from the LES 3D domain and the LES mesh resolution is maintained with the parallel implementation. Fifty three nodes are used to cover the reference space variable, n for the interval [4, 4], clustered around n = 0. A fractional step method is applied with a time step 100 times smaller than the LES time step. For the solver and further analysis purpose (see Section 4), Eq. (1) is expressed in the following form: @X z @X z @2X z ¼ UrX z A þB : |fflfflfflffl{zfflfflfflffl} @t @n @n2 Convection |fflfflfflffl{zfflfflfflffl} |fflfflfflfflffl{zfflfflfflfflffl} Drift
ð6Þ
Diffusion
The Xz profiles were shown to be insensitive to further reduction in the time step. Second order central differencing scheme is used for the discretization of diffusion terms, a second order TVD scheme [33] for the spatial convective term and hybrid scheme for convection in n space. At the inlet (y = 0), Dirichlet boundary conditions are used. At the outlet (y = 30d) and at both ends in the z direction, i.e. z = 10d and z = 10d, zero gradient is imposed. At n = 4 and n = 4 a zero gradient is imposed. Values of U, A and B are calculated explicitly based on the solution at the previous iteration. At every point in space, Xz are initialized in order to match the LES mixing e and its SGS variance. From the solufield for Z e and g tion of Eq. (1), Z Z 002 can be retrieved from Xz at any point in space and time using Z 1 e¼ Z X z P n dn; ð7Þ 1 Z 1 g ðX z hX z i Þ2 P n dn: ð8Þ Z 002 ¼ 1
LES and CMC are run for 80,000 time steps (total physical time of 20 ms) in order to have wellestablished burning conditions in the domain. Then, the mapping functions are initialized such that Eqs. (7) and (8) are satisfied. MMC is run as post-process onto the frozen turbulent flowfield, until steady-state is achieved in the mapping function profiles (see Sections 4.1 and 4.3). Starting from these converged MMC values, LES/CMC and MMC are run together for 10.25 ms to perform time averages (Section 4.2).
C.B. Devaud et al. / Proceedings of the Combustion Institute 34 (2013) 1213–1221
LES/CMC calculations are fully coupled, whereas MMC uses updated LES/CMC fields for inputs g to LES/CMC. only and does not return Njg Defining a flow reference time as 30d/ Ujet = 1.3 ms, 10.25 ms correspond to approximately 8 flow times. 4. Results The LES–CMC predicted mixture fraction, lift-off height, temperature and species concentrations are in good agreement with experimental data [19]. The present section focuses on comparing MMC statistics with LES results. Experimental data is also included for further comparison of g are examined and time averages. Results for Njg compared with those of AMC. 4.1. MMC results with LES–CMC frozen fields MMC is run alone starting from the initial conditions matching the turbulent mixing field predicted by LES. Thus, very small changes, if any, would be expected in Xz profiles and resulting e , and predictions of filtered mixture fraction, Z SGS variance, g Z 002 , if all terms in Eq. (1) are modelled and implemented correctly. Figure 1 shows e and SGS rms instantaneous radial profiles of Z obtained from MMC using Eqs. (7) and (8), respectively, compared with those from LES at three axial positions. For brevity, only three axial positions are selected for the paper, however, results at different locations lead to similar conclusions. Two sets of MMC predictions are included in Fig. 1: ‘regular MMC’ corresponds to MMC with the models presented in Section 2 and ‘scaled production MMC’ includes a scaling factor for production of SGS fluctuations only. Further discussion is to follow. As seen in Fig. 1, MMC e in excellent agreement with LES very predicts Z close to the inlet, at y/d = 1, also noticed for axial
e (left) and SGS Fig. 1. Instantaneous radial profiles of Z Z rms (right). Open symbols correspond to LES solutions, dashed line to regular MMC and solid line MMC with scaled production.
1217
positions y/d 6 5. Further downstream, at y/ d = 8, some discrepancies are seen between e values, in particular in the shear MMC and LES Z layer. MMC gives a centreline value approximately 10% larger compared to LES, but the radial distance at which half of the centreline is reached is underpredicted, approximately by 60%, yielding a much narrower profile than that of LES. At y/d = 26, much further downstream, close to the exit of the domain, MMC predictions e are in good agreement with LES with still a of Z sharper decrease in the shear layer compared to LES. The two sets of MMC values produce simie with a slight advantage to lar radial profiles for Z the ‘scaled production MMC’ for positions further downstream like at y/d = 26 when compared e is shown to LES. Overall, the centreline MMC Z to be in reasonable agreement with LES and MMC values decline faster in the radial direction compared to LES. At intermediate times before steady state, the SGS variance is shown to be much smaller than what is predicted by LES. The SGS fluctuations tend to zero in many locations when ‘regular MMC’ is used. This is not a realistic behaviour for the present LES, being far enough from Direct Numerical Simulation (DNS). When the SGS variance gets close to zero, Xz tends to become flat in e, n space at a constant value corresponding to Z consequently the filtered PDF becomes a delta g determined by Eq. (5) also tends function and Njg to zero. This is problematic for the calculation of g As steady-state is approached, the MMC preNjg. dicted SGS variance with ‘regular MMC’ have improved and the results are shown in Fig. 1 at the three selected axial positions. However, the SGS rms remains underpredicted for most parts of the jet flow using regular MMC, in particular at y/d = 26. The contributions of each term in Eq. (6) are investigated and it is found that the two major terms are convection and drift, with diffusion being much smaller. This is clearly demonstrated by the term budgets of Eq. (6) presented in Fig. 2. This is also qualitatively consistent with previous findings [11]. In deterministic MMC turbulent fluctuations are generated through convection, in particular U(1) determined by the linear model (Eq. (2)), while they are dissipated by the drift term. If initial high subgrid variances disappear too quickly, it means that either the turbulent subgrid fluctuations are not large enough or there is too much dissipation. In MMC, the drift is calculated in order to satisfy the transport equation of Pn, starting from models for the MMC velocity and diffusion coefficient, B. Thus, the source of the present discrepancies points towards the convection modelling and in particular, the linear model. Vaishnavi and Kronenburg [35] propose a different expression for U in order to correct previously observed discrepancies in the near field of
1218
C.B. Devaud et al. / Proceedings of the Combustion Institute 34 (2013) 1213–1221
Fig. 2. Term budget in Eq. (6) at y/d = 0.6 and z/ d = 0.6.
a jet flame. However, this is is still based on the linear model (in g space instead of n space) and was shown to bring improvement in the near-field region only. Several numerical tests have been performed (not shown here for brevity): the production was multiplied by a large constant value or the dissipation was divided by the same coefficient, each time leaving the other terms unchanged. The production (U(1)) is found to be the most influential term for the final MMC predictions. Thus, a multiplying coefficient equal to e g the ratio of LES values for ð Z Z 002 Þ to the MMC data is applied to the production term only to ensure that more plausible levels of SGS fluctuation are generated in Eq. (1), resulting in 002 e f Uð1Þ hnX z i ¼ ðD þ Dt Þ$hX z i ð Z LES Z LES Þ . All the f ðe Z MMC Z 002 MMC Þ other terms remain unchanged. The resulting MMC predictions are also shown in Fig. 1. The new scaling has a small impact on the values of e , except at y/d = 26. Some improvements can Z be seen for the SGS variances compared to regular MMC. This new multiplying coefficient brings the MMC predicted SGS rms closer to the LES results. A better multiplying coefficient may be found by trial and error, however, tuning constants in the linear model may not be the most suitable solution for future work. As shown in Section 2, the closure of the MMC coefficients must be consistent with the PDF equation. Thus, the present scaling does not satisfy this requirement anymore and this problem needs to addressed in future work. In Sections 4.2 and 4.3, the MMC results include this correction on the production term. Figure 2 presents the contribution of each term in Eq. (6) using the converged steady state solution at y/d = 0.6 and z/d = 0.6. When steady-state is reached, there is a balance mainly between drift and convection. Initially, convection is larger than drift and diffusion. Thus, convection is the term that creates the unbalance in the MMC transport
equation with (U(1)n) having a wide range of magnitudes depending on the location. At y/ d = 0.6 and z/d = 0.6, U(1) corresponds to 0.5% and 35% of the filtered velocity U(0) in the y and z direction, respectively. Thus, (U(1)n) is not negligible although its influence decreases as jnj tends to zero where the Pn peak is located. This may explain why the scaling of U(1) does not have a e , as seen in big impact on the predictions of Z Fig. 1. The linear model may still remain an issue for future work in terms of SGS variance predictions but may not be solely responsible for the dise . The crepancies in the radial profiles of Z remaining explanation may be the fact that MMC is run in 2D in space instead of 3D. Cartesian coordinates are used instead of axisymmetric coordinates, consequently spatial transport in the x-direction (one of the transverse directions) may not be small and should be included in the MMC calculations. Contributions of each convective component in x, y and z directions are determined using the LES velocity and mixing fields on the centre plane, i.e. same plane used for the present MMC calculations. Convection in the x direction cannot be considered negligible in many flow regions. In particular, convection in the y (axial) direction and that in the x direction are shown to have opposite signs at several points, which would reduce the net convection. Reduced magnitude of convection in MMC would help decrease the unbalance presently observed at early times and changes in Xz profiles from initial values. This would mean that neglecting convection in the x direction is not sufficiently accurate and causes discrepancies in the predictions of the mixture fraction field. There is no direct method to clearly quantify which source of modelling discrepancy between the LES and MMC results has the largest impact: the linear model or the neglect of spatial transport in the third direction in MMC. The numerical implementation of the transport equations in LES and MMC is similar: fractional step method, second-order scheme for spatial convection and diffusion and same spatial grid resolution in both approaches and MMC results were shown to be insensitive to further mesh refinement in space and n. Thus, no further significant numerical and modelling discrepancy is expected between LES and MMC. 4.2. Time averages Radial profiles of time averaged mixture fraction, Z, and corresponding resolved rms are shown in Fig. 3 at three axial positions. Each plot presents values obtained from MMC, LES (extracted from the centre plane on which Eq. (1) is solved) compared with the experimental data. Similar trends are seen between timeaveraged profiles of Z and instantaneous values shown in Fig. 1. Excellent agreement is obtained
C.B. Devaud et al. / Proceedings of the Combustion Institute 34 (2013) 1213–1221
between MMC, LES and the experiments close to the inlet (y/d = 1). At y/d = 8 and y/d = 26, MMC underpredicts the centreline values and radial spreading rates compared to LES and experimental values although the spreading rate is better captured in the time-averaged profiles than what is observed in Fig. 1. Again, the observed discrepancies may be explained by the fact that the current MMC formulation does not include spatial variations in the x-direction. MMC tends to overpredict the experimental rms in the jet core, in particular at y/d = 1 and y/ d = 8 but a reasonable match can be observed at y/d = 26. The current predictions are consistent with previous MMC studies, either in LES or RANS: in LES the mixture fraction rms are found to be overpredicted [18] and in RANS overpredicted in the near-field and underpredicted further downstream from the jet exit [11]. 4.3. Conditionally filtered scalar dissipation rate g obtained from Figure 4 presents profiles of Njg MMC (initial values and values after reaching steady-state) and AMC at two points, Point 1 and Point 2, along with the corresponding PDF from MMC and b PDF. Physical coordinates, vale and SGS variance for Point 1 and 2 are ues for Z given in Table 1. Point 1 is located close to the inlet and represents very rich mixture with nonnegligible SGS variance, whereas Point 2 is situated in a much leaner mixture. In Fig. 4a, the characteristic bell shaped profile of the AMC model can be seen, always centred at a mixture fraction of 0.5. The MMC profiles are significantly different from that of AMC. They are not e . The symmetrical and peak values occur near Z g MMC Njg values decrease to zero for g smaller than 0.4, consistent with what is seen in the PDF. MMC does not set a predefined shape for
e (left) and Z Fig. 3. Time-averaged radial profiles of Z rms (right). Open symbols correspond to experimental data, dashed line to LES and solid line MMC.
1219
g using AMC and MMC at Fig. 4. Left: profiles of Njg two locations. Right: corresponding profiles of PDF retrieved from MMC and calculated from beta function.
Table 1 Specifications of Point 1 and 2. y/d z/d e from LES Z e from MMC Z SGS rms from LES SGS rms from MMC
Point 1
Point 2
0.6 0.6 0.864 0.76 0.094 0.086
26 2.4 0.15 0.15 0.016 0.018
g and its shape closely follows the mixing state Njg described by the PDF. According to what is observed for the SGS rms in Fig. 1, the converged MMC values are often lower than the initial data e and g (initial MMC values for Z Z 002 are equal to LES data) due to lower SGS fluctuations. At Point 1 the converged MMC g Z 002 has decreased by 10%, resulting in a much larger change in the g approximately 50% decrease. magnitude of Njg, g from MMC is very sensitive This shows that Njg to the SGS variance. The peak of the corresponding MMC PDF has shifted to richer mixtures due e moving from 0.86 to to the MMC predicted Z 0.76. At Point 1, the MMC PDF is different from that of the b form. At these very rich mixtures far away from possible chemical reactions, the accuracy of the PDF is not as crucial as for mixtures closer to stoichiometry. Experimental or DNS data would be needed to determine what the correct profile is. At Point 2 (Fig. 4b), similar to what is found at Point 1, the MMC profiles (both coincide due to negligible change in the MMC values between inig are different in tial and converged data) for Njg g shape from AMC. The MMC Njg results go down
1220
C.B. Devaud et al. / Proceedings of the Combustion Institute 34 (2013) 1213–1221
to zero for g larger than 0.2 and smaller than 0.1, i.e. where the PDF also tends to zero. The b and MMC PDF are in close agreement, the b function only showing a slightly higher peak value, consistent e is with previous findings [11]. At both points, N R1 g recovered in MMC from 0 ð NjgÞP ðgÞdg, using g and P(g) (FDF) from MMC. No scaling is Njg e at Point 1 and Point 2, in conneeded to recover N trast to previous approaches based on PDF transport equation integration [5]. 5. Conclusion Deterministic MMC has been applied in LES for the first time to a lifted flame in a vitiated coflow. MMC has some interesting features like g no modelling assumption for the FDF and Njg, which makes it very attractive as a mixing model. Some modelling and numerical challenges, in particular in the context of LES, have been examined in the present paper. This first analysis includes MMC instantaneous results for mixture fraction, g and FDF starting from frozen SGS rms, Njg LES–CMC well-established conditions. Timeaveraged statistics are also collected for mixture fraction and resolved rms. Reasonable agreement e and resolved rms is found for radial profiles of Z between MMC, LES and experimental data. Best agreement is found close to the fuel exit and some discrepancies can be seen further downstream. It is believed that the current approximation of neglecting spatial transport in one of the transverse directions is mainly responsible for the observed discrepancies. Despite the reduced dimensionality implemented in the current MMC, it is very positive to see that deterministic MMC can reproduce the main turbulent mixing characteristics. Future work will focus on performing MMC in three dimensions. The resulting increased computational cost will be monitored and minimized. The second item for future investigation will be the linear model and alternatives to generate SGS fluctuations in MMC. For example, a gradient model may be implemented [31]. Accurately reproducing the level of SGS variance in MMC proves to be challenging but a key to accurate preg and the turbulent mixing field. Furdiction of Njg g are shown to be ther, the calculated values of Njg very sensitive to the level of SGS fluctuations. Thus, it is crucial to get good predictions for SGS variance. g are significantly difThe MMC profiles of Njg ferent from those of AMC, which certainly prompts further comparison with available experimental data and investigation of their impact on the CMC calculations for predicted lift-off height and species concentrations.
Acknowledgment The authors would like to thank Dr. K. Vogiatzaki for helpful discussions on MMC, as well as Dr. A. Klimenko. Funding from FWO project G.0079.07 is gratefully acknowledged. References [1] R. Fox, Computational Models for Turbulent Reacting Flows, Cambridge University Press, 2003. [2] A. Klimenko, R. Bilger, Prog. Energy Combust. Sci. 25 (6) (1999) 595–687. [3] S. Girimaji, Phys. Fluids A 4 (11) (1992) 2529–2537. [4] E. O’Brien, T. Jiang, Phys. Fluids A 3 (12) (1991) 3121–3123. [5] A. Milford, C. Devaud, Combust. Flame 157 (2010) 1467–1483. [6] C. Devaud, R. Bilger, T. Liu, Phys. Fluids 16 (6) (2004) 2004–2011. [7] M. Cleary, J. Kent, Combust. Flame 143 (2005) 357–368. [8] S. Sreedhara, Y. Lee, K.Y. Huh, D. Ahn, Combust. Flame 152 (2007) 282–286. [9] A.Y. Klimenko, S. Pope, Phys. Fluids (2003) 1907– 1925. [10] K. Vogiatzaki, A. Kronenburg, M.J. Cleary, J.H. Kent, Proc. Combust. Inst. 32 (2009) 1679–1685. [11] K. Vogiatzaki, M.J. Cleary, A. Kronenburg, J.H. Kent, Phys. Fluids 21 (2009) 025105. [12] A. Cook, J. Riley, Phys. Fluids 6 (1994) 2868–2870. [13] A. Triantafyllidis, E. Mastorakos, R. Eggels, Combust. Flame 156 (2009) 2328–2345. [14] I. Stankovic, A. Triantafyllidis, E. Mastorakos, C. Lacor, B. Merci, Flow Turbul. Combust. 86 (2011) 689–710. [15] S. Navarro-Martinez, S. Rigopoulos, Flow Turbul. Combust. 87 (2011) 407–423. [16] A. Garmory, E. Mastorakos, Proc. Combust. Inst. 33 (2011) 1673–1680. [17] M. Cleary, A. Klimenko, J. Janicka, M. Pfitzner, Proc. Combust. Inst. 32 (2009) 1499–1507. [18] M. Cleary, A. Klimenko, Phys. Fluids 23 (2011) 115102. [19] I. Stankovic, Numerical Simulations of Hydrogen Auto-ignition in Turbulent Flows, Ph.D. thesis, Ghent University, Belgium, 2011. [20] C. Markides, E. Mastorakos, Proc. Combust. Inst. 30 (2005) 883–891. [21] R. Cabra, T. Myhrvold, J. Chen, R. Dibble, A. Karpetis, R. Barlow, Proc. Combust. Inst. 29 (2002) 1881–1888. [22] R. Gordon, A. Masri, S. Pope, G. Goldin, Combust. Theory Model. 11 (2007) 351–376. [23] R. Cao, S. Pope, A. Masri, Combust. Flame 142 (2005) 438–454. [24] K. Gkagkas, R. Lindstedt, Combust. Theory Model. 13 (2009) 607–643. [25] A.R. Masri, R. Cao, S.B. Pope, G.M. Goldin, Combust. Theory Model. 8 (2004) 1–22. [26] T. Myhrvold, I. Ertesvag, I. Gran, R. Cabra, J. Chen, Combust. Sci. Tech. 178 (2009) 1001–1030. [27] S. Patwardhan, D. Santanu, K. Lakshmisha, Proc. Combust. Inst. 32 (2009) 1708–1712. [28] C. Duwig, L. Fuchs, Combust. Sci. Tech. 180 (2008) 453–480.
C.B. Devaud et al. / Proceedings of the Combustion Institute 34 (2013) 1213–1221 [29] W. Jones, S. Navarro-Martinez, Combust. Flame 150 (2007) 170–187. [30] S. Navarro-Martinez, A. Kronenburg, Flow Turbul. Combust. 87 (2011) 377–406. [31] M. Cleary, A. Klimenko, in: T. Echekki, E. Mastorakos (Eds.), Turbulent Combustion Modelling Advances, New Trends and Perspectives, Springer, 2011.
1221
[32] A. Triantafyllidis, E. Mastorakos, Flow Turbul. Combust. 84 (2009) 481–512. [33] B. Van Leer, J. Comput. Phys. 14 (1974) 361–370. [34] J. Li, Z. Zhao, A. Kazakov, F.L. Dryer, Int. J. Chem. Kinet. 36 (2004) 566–575. [35] P. Vaishnavi, A. Kronenburg, Combust. Flame 157 (2010) 1863–1865.
Proceedings of the
Combustion Institute
Proceedings of the Combustion Institute 34 (2013) 1213–1221
www.elsevier.com/locate/proci
Deterministic Multiple Mapping Conditioning (MMC) applied to a turbulent flame in Large Eddy Simulation (LES) C.B. Devaud a,⇑, I. Stankovic´ b, B. Merci b a
Department of Mechanical and Mechatronics Engineering, University of Waterloo, 200 University Avenue West, Waterloo, ON, Canada N2L3G1 b Ghent University, Department of Flow, Heat and Combustion Mechanics, Belgium Available online 2 July 2012
Abstract The scalar dissipation rate is a key quantity in turbulent combustion modelling, in particular for Conditional Moment Closure (CMC). In Large Eddy Simulation (LES), a subgrid-scale (SGS) mixing model must be provided to determine the SGS mixture fraction variance needed for the filtered scalar dissipation rate and the filtered density function. Multiple Mapping Conditioning (MMC) is selected as a SGS mixing model and the deterministic form is implemented for the first time in LES. The objective of the present study is to assess the capabilities of MMC to model the turbulent mixing field and the conditionally filtered g needed in CMC. scalar dissipation rate, Njg The MMC transport equation is implemented in a LES code coupled with CMC to simulate a lifted jet flame in a vitiated coflow. The MMC formulation and computational details are presented. MMC predictions are compared with LES for the turbulent mixing field using instantaneous and time averaged data, and experimental values are included whenever available. Reasonable agreement between MMC, LES and experimental data is found. Best fit is obtained close to the fuel exit and some discrepancies can be seen further downstream. The SGS variance is underpredicted for many parts of the jet. A thorough analysis is performed to assess the modelling and numerical assumptions. The two proposed sources of discrepancy are the linear model and neglect of a transverse component g are significantly different from those of AMC and are for convection in MMC. The MMC profiles of Njg sensitive to the SGS variance. Future work is suggested. Ó 2012 The Combustion Institute. Published by Elsevier Inc. All rights reserved. Keywords: Multiple Mapping Conditioning; LES; Scalar dissipation rate; Conditional Moment Closure; Mixing
1. Introduction The accurate representation of the complex interactions between chemistry and turbulence is ⇑ Corresponding author. Fax: +1 519 885 5862.
E-mail Devaud).
address:
[email protected]
(C.B.
crucial for the development of industrial combustion devices. Turbulent mixing is of particular importance in non-premixed combustion where chemical reactions will only occur if fuel and oxidizer are mixed at the molecular level. The current computational tools for engineering applications rely on Reynolds Averaged Navier–Stokes (RANS) and Large Eddy Simulations (LES)
1540-7489/$ - see front matter Ó 2012 The Combustion Institute. Published by Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.proci.2012.06.076
1214
C.B. Devaud et al. / Proceedings of the Combustion Institute 34 (2013) 1213–1221
approaches where, in both cases, the turbulent small scales need to be modelled. The scalar dissipation rate, N, represents the rate of mixing at the molecular level and is proportional to the square gradient of the scalar, Z, such as N D$Z $Z, where D is the molecular diffusivity of Z. Thus, this is a key quantity in turbulent combustion modelling, in particular for flamelet, Conditional Moment Closure (CMC) and Probability Density Function (PDF) approaches [1]. Within the CMC and turbulent non-premixed combustion framework, the conditional average of N at a particular value, g, of mixture fraction, Z is of special interest with hNjgi = hD$Z $ZjZ = gi. The angular brackets denote a conditional average over an ensemble of realizations of the flow, subject to the condition to the right of the vertical bar. Evaluation of hNjgi is not straightforward and CMC further requires solution of the transport equations for the conditional averages to be consistent with that of the PDF transport equation [2]. Two of the most commonly used models in CMC are the presumed b-PDF model of Girimaji [3] and Amplitude Mapping Closure (AMC) [4]. Both models are derived assuming homogeneous turbulence and consequently, do not comply with the solution of the PDF transport equation in inhomogeneous turbulent flow. Other formulations for hNjgi have been proposed and can be divided into two groups. The first approach is based on integrating the PDF transport equation to enforce the required consistency in CMC [5–8]. This shows some success but suffers from potentially large numerical errors and inaccuracies due to the conditional velocity model [5]. The second strategy relies on Multiple Mapping Conditioning (MMC), initially proposed by Klimenko and Pope [9], in which hNjgi appears in closed form without homogeneous turbulence assumption or complex integration. Vogiatzaki et al. [10,11] implemented MMC with a single major scalar, mixture fraction, and calculated the unconditional and conditional scalar dissipation rate for simple jet flames. These results show good agreement with experimental data. All these past investigations were conducted in the context of RANS. In LES, the large scale turbulent motions are resolved but a subgrid-scale (SGS) mixing model must be provided to determine the SGS mixture fraction variance needed for the filtered scalar dissipation rate and the filtered density function (FDF) when a presumed form is considered for the FDF [12]. In LES–CMC, it is common practice to include the AMC model based on homogeneous turbulence to determine the conditionally filtered scalar g [13,14]. Alternative modelling dissipation Njg may use the conditional average of filtered values over one CMC cell [15]. Recent LES–CMC results for Sandia Flame F are found to be sensitive to g and a model constant in the modelling of Njg the SGS component of N was adjusted [16].
Encouraging results have been recently reported for Sandia Flame D using the stochastic formulation of MMC with a Lagrangian numerical scheme in LES [17]. Computational cost is minimized by using sparse Lagrangian technique. Good predictions for the filtered mixture fraction are obtained while the rms are overpredicted compared to the experimental data. Further, the results are shown to be sensitive to model constants in the particle mixing model and no density feedback is implemented. Cleary and Klimenko [18] further expanded the analysis for constant and variable density flows with density coupling. Comparisons were made between LES and MMC for mean mixture fraction and its rms. Similar to previous observations [17], good agreement is found for the mean mixture fraction but the rms is overpredicted. The reactive scalar predictions are shown to be sensitive to particle mixing length scales. The objective of the present investigation is to assess the capabilities of MMC to model the turg needed in CMC, in bulent mixing field and Njg LES. In contrast to previous LES–MMC studies, the deterministic form is implemented for the first time in LES and incorporates the approach initially developed in RANS. LES is used to determine the turbulent flow and mixing fields, CMC is included to simulate turbulent combustion and in parallel MMC is implemented to emulate the turbulent mixing field and determine the conditionally filtered scalar dissipation rate. For best g found in assessment of MMC performance, Njg MMC is not coupled back to CMC (CMC retains the AMC model). The current MMC investigation builds upon previous LES–CMC studies [14,19] applied to two experimental cases, hydrogen autoignition in a coflow of heated air [20] and the Cabra lifted hydrogen flame in a vitiated coflow [21]. The latter is kept for the current MMC analysis. Previous numerical work can also be found for the hydrogen Cabra flame using different combustion models in RANS [21–27] and LES [28–30]. Autoignition is shown to be the primary stabilization mechanism [19,22] and g has the potential of improved modelling of Njg strongly affecting the CMC predictions, in particular for the lift-off height. In the next sections, the MMC formulation and computational details are presented. The results focus on comparing MMC with LES for the turbulent mixing field using instantaneous and time averaged data, and experimental values are included whenever available. 2. Multiple Mapping Conditioning (MMC) MMC is derived from the idea that all turbulent fluctuations can be divided into major and minor groups and turbulent fluctuations of minor scalars are correlated with those of the major scalars
C.B. Devaud et al. / Proceedings of the Combustion Institute 34 (2013) 1213–1221
[9,31]. A reference variable is assigned to each major scalar. In the present work, one major scalar is selected, the conditioning variable used in CMC, namely mixture fraction. Thus, it is assumed that the fluctuations of all other quantities in MMC, such as temperature and species concentrations are correlated to those of the mixture fraction. This is the same idea behind the single conditioning CMC. The MMC equations solve for mapping functions which map between the reference space and the composition space following the same principle of mapping closure methods. The equations are derived for inhomogeneous turbulence and do not depend on the shape of the reference PDF. The reference variable is noted as n. Xz is the mapping function solved by the MMC equation and Xz maps between n and Z, mixture fraction. MMC offers additional benefits beyond the scope of the current investigation. It can be easily formulated for two or more conditioning variables [31], useful for the modelling of partially premixed turbulent combustion. The MMC equation for the spatial and temporal evolution of Xz is given by [9] @X z @X z @2X z þ UrX z þ A B ¼ 0; @t @n @n2
ð1Þ
where U, A and B correspond to velocity, drift and diffusion, respectively, and need to be closed. The present MMC formulation in LES is very similar to that in RANS, also consistent with similar form of the CMC equations in RANS and LES. Closure of the unclosed terms is obtained through consistency with the reference PDF transport equation. A spatially and temporally invariant Gaussian form with zero mean and unity variance is selected for the reference PDF, Pn, following common practice [9–11]. The coefficients can be determined in principle for any form of Pn. However, the choice of Pn has an impact on the resulting complexity for the expressions of U and A. All MMC coefficients must satisfy the transport equation for Pn. Different expressions may be obtained and depend on the assumptions made for U and B. In the present work, models suggested by Klimenko and Pope [9] and implemented in RANS by Vogiatzaki et al. [10,11] are used. A linear model is applied for U given by U = U(n, x, t) = U(0) + U(1)n, where U(0) corresponds to the filtered velocity vector solved in LES and U(1) is the velocity gradient. The expression for U(1) is determined such that the first moment of the filtered mixture fraction PDF is recovered following the procedure explained by Vogiatzaki et al. [11]: Uð1Þ hnX z i ¼ ðD þ Dt Þ$hX z i ;
ð2Þ
where D and Dt are the molecular and turbulent diffusivity, respectively. The terms with overtilde are filtered quantities and terms in angular brackets with a star are averages obtained by integrating
1215
the product of any function of n with Pn over n space. In Eq. (2), hXzi* is determined by R1 X ðnÞP z n dn. The linear model is commonly 1 used in RANS for CMC and MMC [5,10]. It is exact for jointly Gaussian turbulence and is considered a good approximation for flows within two to four standard deviations of the scalar under consideration [5]. However, its validity has not been demonstrated in LES [32]. The drift and diffusion coefficients are determined by 1 @X z @X z ð1Þ e: A ¼ Bn þ rqU ; B ¼N ð3Þ q @n @n The expression for B in Eq. (3) requires e obtained from LES and includes knowledge of N contribution of resolved and subgrid gradients: e e e ¼ m þ mt @ Z @ Z ; N ð4Þ Sc Sct @xk @xk where m is the molecular kinematic viscosity, mt the turbulent viscosity, Sc is the molecular Schmidt number and Sct the turbulent Schmidt number. g needed in CMC, does not appear explicitly Njg, in the MMC formulation. It can be retrieved through a transformation from the reference space, n, to the mixture fraction space, g, such that [10,11]: 2 g B @X z : Njg ð5Þ @n 1 z Pe ðgÞ may be obtained using Pe ðgÞ ¼ P n dX . dn 3. Computational details 3.1. Experimental case configuration The selected experimental case is a turbulent lifted jet flame of hydrogen diluted with nitrogen issuing into a wide co-flow of vitiated air [21]. The burner consists of a fuel jet nozzle (inner diameter, d = 4.57 mm) and a surrounding perforated plate disk (outer diameter 210 mm). The fuel stream consists of hydrogen diluted with nitrogen with exit velocity equal to 107 m/s at 305 K. The vitiated air consists of the products of a lean premixed hydrogen/air flame at 1045 K with a velocity of 3.5 m/s. The Reynolds number is equal to 23,600 for the fuel jet and 18,600 for the coflow. The stoichiometric mixture fraction is 0.474. 3.2. LES–CMC The Favre filtered equations for mass, momentum and mixture fraction, Z, are solved in three dimensions in Cartesian coordinates in the current LES code [14,19]. Convective transport is discretized with a second order central scheme for the momentum equations and a second order
1216
C.B. Devaud et al. / Proceedings of the Combustion Institute 34 (2013) 1213–1221
TVD (Total Variation Diminishing) [33] for the mixture fraction equation. Subgrid scale stresses are modelled using a standard Smagorinsky model with a model constant, C, equal to 0.1. Equal diffusivities, constant Schmidt number (Sc = Sct = 0.7) and unity Lewis number are applied. The SGS mixture fraction variance is obtained by a gradient type model, 2 @e g 002 Z @e Z Z ¼ CD @xk @xk , where D is the mesh dependent filter width. The expression for SGS variance assumes equilibrium between production and dissipation. Other approaches may be considered [32] but are beyond the scope of the present investigation. The computational domain extends 30d downstream from the jet inlet (approximately 137 mm) in the axial direction and 20d (91.4 mm) radially. The corresponding grid contains 192 48 48 cells. The grid is stretched smoothly towards the co-flow in the radial direction and is expanded smoothly in the axial direction. A digital filter is used to generate turbulence in the coflow with a turbulence lengthscale equal to 1.6 mm, the size of the holes in the outer disk. The turbulence intensity is set to 5% following experimental information [21] and previous simulations for the same flame [25]. At the inlet of the domain, Dirichlet boundary conditions are used for the velocity and mixture fraction. At the outlet, Neumann boundary conditions are applied for all quantities except for pressure which is imposed (atmospheric pressure). The implementation is in parallel with 4 blocks in the axial direction. Turbulent combustion is modelled using a first-order CMC approach with detailed mechanism including 9 species and 19 reversible reactions [34]. The equations for the conditionally filtered species mass fractions and temperature are solved. Integration with the density-weighted filtered b PDF, Pe ðgÞ over mixture fraction space yields the unconditional density, temperature and species mass fractions. The LES and CMC solvers are coupled and the information from CMC is transferred at every timestep of the simulation to the LES in order to update flow field. A coarser spatial mesh is used to solve the CMC equations due to much weaker spatial variations of conditional averages compared to those of Favre filtered scalars. The CMC equations are discretized in finite differences. The CMC grid consists of 80 5 5 cells and the mixture fraction space is discretized into 50 bins, clustered around the most reactive mixture fraction (equal to 0.053). Information is transferred from the finer LES spatial mesh to the coarser CMC grid using volume averaging [14,32]. The AMC model is applied directly to the CMC cells in order to g At the inlet of the CMC domain, inert obtain Njg. distributions for conditional averages are set. Zero gradient boundary condition is imposed at
the sides and the downstream edge of the domain in physical space. In g space, a Dirichlet boundary condition is used at g = 0 and g = 1. 3.3. MMC numerical implementation Due to increased memory allocation, capabilities of our current computational system and significant increase in computational time, the MMC equation, Eq. (1) is solved in two spatial dimensions: the axial (y axis) and transverse (z axis) directions are kept. Transport in the third (x, transverse) direction is ignored. The centre plane is selected from the LES 3D domain and the LES mesh resolution is maintained with the parallel implementation. Fifty three nodes are used to cover the reference space variable, n for the interval [4, 4], clustered around n = 0. A fractional step method is applied with a time step 100 times smaller than the LES time step. For the solver and further analysis purpose (see Section 4), Eq. (1) is expressed in the following form: @X z @X z @2X z ¼ UrX z A þB : |fflfflfflffl{zfflfflfflffl} @t @n @n2 Convection |fflfflfflffl{zfflfflfflffl} |fflfflfflfflffl{zfflfflfflfflffl} Drift
ð6Þ
Diffusion
The Xz profiles were shown to be insensitive to further reduction in the time step. Second order central differencing scheme is used for the discretization of diffusion terms, a second order TVD scheme [33] for the spatial convective term and hybrid scheme for convection in n space. At the inlet (y = 0), Dirichlet boundary conditions are used. At the outlet (y = 30d) and at both ends in the z direction, i.e. z = 10d and z = 10d, zero gradient is imposed. At n = 4 and n = 4 a zero gradient is imposed. Values of U, A and B are calculated explicitly based on the solution at the previous iteration. At every point in space, Xz are initialized in order to match the LES mixing e and its SGS variance. From the solufield for Z e and g tion of Eq. (1), Z Z 002 can be retrieved from Xz at any point in space and time using Z 1 e¼ Z X z P n dn; ð7Þ 1 Z 1 g ðX z hX z i Þ2 P n dn: ð8Þ Z 002 ¼ 1
LES and CMC are run for 80,000 time steps (total physical time of 20 ms) in order to have wellestablished burning conditions in the domain. Then, the mapping functions are initialized such that Eqs. (7) and (8) are satisfied. MMC is run as post-process onto the frozen turbulent flowfield, until steady-state is achieved in the mapping function profiles (see Sections 4.1 and 4.3). Starting from these converged MMC values, LES/CMC and MMC are run together for 10.25 ms to perform time averages (Section 4.2).
C.B. Devaud et al. / Proceedings of the Combustion Institute 34 (2013) 1213–1221
LES/CMC calculations are fully coupled, whereas MMC uses updated LES/CMC fields for inputs g to LES/CMC. only and does not return Njg Defining a flow reference time as 30d/ Ujet = 1.3 ms, 10.25 ms correspond to approximately 8 flow times. 4. Results The LES–CMC predicted mixture fraction, lift-off height, temperature and species concentrations are in good agreement with experimental data [19]. The present section focuses on comparing MMC statistics with LES results. Experimental data is also included for further comparison of g are examined and time averages. Results for Njg compared with those of AMC. 4.1. MMC results with LES–CMC frozen fields MMC is run alone starting from the initial conditions matching the turbulent mixing field predicted by LES. Thus, very small changes, if any, would be expected in Xz profiles and resulting e , and predictions of filtered mixture fraction, Z SGS variance, g Z 002 , if all terms in Eq. (1) are modelled and implemented correctly. Figure 1 shows e and SGS rms instantaneous radial profiles of Z obtained from MMC using Eqs. (7) and (8), respectively, compared with those from LES at three axial positions. For brevity, only three axial positions are selected for the paper, however, results at different locations lead to similar conclusions. Two sets of MMC predictions are included in Fig. 1: ‘regular MMC’ corresponds to MMC with the models presented in Section 2 and ‘scaled production MMC’ includes a scaling factor for production of SGS fluctuations only. Further discussion is to follow. As seen in Fig. 1, MMC e in excellent agreement with LES very predicts Z close to the inlet, at y/d = 1, also noticed for axial
e (left) and SGS Fig. 1. Instantaneous radial profiles of Z Z rms (right). Open symbols correspond to LES solutions, dashed line to regular MMC and solid line MMC with scaled production.
1217
positions y/d 6 5. Further downstream, at y/ d = 8, some discrepancies are seen between e values, in particular in the shear MMC and LES Z layer. MMC gives a centreline value approximately 10% larger compared to LES, but the radial distance at which half of the centreline is reached is underpredicted, approximately by 60%, yielding a much narrower profile than that of LES. At y/d = 26, much further downstream, close to the exit of the domain, MMC predictions e are in good agreement with LES with still a of Z sharper decrease in the shear layer compared to LES. The two sets of MMC values produce simie with a slight advantage to lar radial profiles for Z the ‘scaled production MMC’ for positions further downstream like at y/d = 26 when compared e is shown to LES. Overall, the centreline MMC Z to be in reasonable agreement with LES and MMC values decline faster in the radial direction compared to LES. At intermediate times before steady state, the SGS variance is shown to be much smaller than what is predicted by LES. The SGS fluctuations tend to zero in many locations when ‘regular MMC’ is used. This is not a realistic behaviour for the present LES, being far enough from Direct Numerical Simulation (DNS). When the SGS variance gets close to zero, Xz tends to become flat in e, n space at a constant value corresponding to Z consequently the filtered PDF becomes a delta g determined by Eq. (5) also tends function and Njg to zero. This is problematic for the calculation of g As steady-state is approached, the MMC preNjg. dicted SGS variance with ‘regular MMC’ have improved and the results are shown in Fig. 1 at the three selected axial positions. However, the SGS rms remains underpredicted for most parts of the jet flow using regular MMC, in particular at y/d = 26. The contributions of each term in Eq. (6) are investigated and it is found that the two major terms are convection and drift, with diffusion being much smaller. This is clearly demonstrated by the term budgets of Eq. (6) presented in Fig. 2. This is also qualitatively consistent with previous findings [11]. In deterministic MMC turbulent fluctuations are generated through convection, in particular U(1) determined by the linear model (Eq. (2)), while they are dissipated by the drift term. If initial high subgrid variances disappear too quickly, it means that either the turbulent subgrid fluctuations are not large enough or there is too much dissipation. In MMC, the drift is calculated in order to satisfy the transport equation of Pn, starting from models for the MMC velocity and diffusion coefficient, B. Thus, the source of the present discrepancies points towards the convection modelling and in particular, the linear model. Vaishnavi and Kronenburg [35] propose a different expression for U in order to correct previously observed discrepancies in the near field of
1218
C.B. Devaud et al. / Proceedings of the Combustion Institute 34 (2013) 1213–1221
Fig. 2. Term budget in Eq. (6) at y/d = 0.6 and z/ d = 0.6.
a jet flame. However, this is is still based on the linear model (in g space instead of n space) and was shown to bring improvement in the near-field region only. Several numerical tests have been performed (not shown here for brevity): the production was multiplied by a large constant value or the dissipation was divided by the same coefficient, each time leaving the other terms unchanged. The production (U(1)) is found to be the most influential term for the final MMC predictions. Thus, a multiplying coefficient equal to e g the ratio of LES values for ð Z Z 002 Þ to the MMC data is applied to the production term only to ensure that more plausible levels of SGS fluctuation are generated in Eq. (1), resulting in 002 e f Uð1Þ hnX z i ¼ ðD þ Dt Þ$hX z i ð Z LES Z LES Þ . All the f ðe Z MMC Z 002 MMC Þ other terms remain unchanged. The resulting MMC predictions are also shown in Fig. 1. The new scaling has a small impact on the values of e , except at y/d = 26. Some improvements can Z be seen for the SGS variances compared to regular MMC. This new multiplying coefficient brings the MMC predicted SGS rms closer to the LES results. A better multiplying coefficient may be found by trial and error, however, tuning constants in the linear model may not be the most suitable solution for future work. As shown in Section 2, the closure of the MMC coefficients must be consistent with the PDF equation. Thus, the present scaling does not satisfy this requirement anymore and this problem needs to addressed in future work. In Sections 4.2 and 4.3, the MMC results include this correction on the production term. Figure 2 presents the contribution of each term in Eq. (6) using the converged steady state solution at y/d = 0.6 and z/d = 0.6. When steady-state is reached, there is a balance mainly between drift and convection. Initially, convection is larger than drift and diffusion. Thus, convection is the term that creates the unbalance in the MMC transport
equation with (U(1)n) having a wide range of magnitudes depending on the location. At y/ d = 0.6 and z/d = 0.6, U(1) corresponds to 0.5% and 35% of the filtered velocity U(0) in the y and z direction, respectively. Thus, (U(1)n) is not negligible although its influence decreases as jnj tends to zero where the Pn peak is located. This may explain why the scaling of U(1) does not have a e , as seen in big impact on the predictions of Z Fig. 1. The linear model may still remain an issue for future work in terms of SGS variance predictions but may not be solely responsible for the dise . The crepancies in the radial profiles of Z remaining explanation may be the fact that MMC is run in 2D in space instead of 3D. Cartesian coordinates are used instead of axisymmetric coordinates, consequently spatial transport in the x-direction (one of the transverse directions) may not be small and should be included in the MMC calculations. Contributions of each convective component in x, y and z directions are determined using the LES velocity and mixing fields on the centre plane, i.e. same plane used for the present MMC calculations. Convection in the x direction cannot be considered negligible in many flow regions. In particular, convection in the y (axial) direction and that in the x direction are shown to have opposite signs at several points, which would reduce the net convection. Reduced magnitude of convection in MMC would help decrease the unbalance presently observed at early times and changes in Xz profiles from initial values. This would mean that neglecting convection in the x direction is not sufficiently accurate and causes discrepancies in the predictions of the mixture fraction field. There is no direct method to clearly quantify which source of modelling discrepancy between the LES and MMC results has the largest impact: the linear model or the neglect of spatial transport in the third direction in MMC. The numerical implementation of the transport equations in LES and MMC is similar: fractional step method, second-order scheme for spatial convection and diffusion and same spatial grid resolution in both approaches and MMC results were shown to be insensitive to further mesh refinement in space and n. Thus, no further significant numerical and modelling discrepancy is expected between LES and MMC. 4.2. Time averages Radial profiles of time averaged mixture fraction, Z, and corresponding resolved rms are shown in Fig. 3 at three axial positions. Each plot presents values obtained from MMC, LES (extracted from the centre plane on which Eq. (1) is solved) compared with the experimental data. Similar trends are seen between timeaveraged profiles of Z and instantaneous values shown in Fig. 1. Excellent agreement is obtained
C.B. Devaud et al. / Proceedings of the Combustion Institute 34 (2013) 1213–1221
between MMC, LES and the experiments close to the inlet (y/d = 1). At y/d = 8 and y/d = 26, MMC underpredicts the centreline values and radial spreading rates compared to LES and experimental values although the spreading rate is better captured in the time-averaged profiles than what is observed in Fig. 1. Again, the observed discrepancies may be explained by the fact that the current MMC formulation does not include spatial variations in the x-direction. MMC tends to overpredict the experimental rms in the jet core, in particular at y/d = 1 and y/ d = 8 but a reasonable match can be observed at y/d = 26. The current predictions are consistent with previous MMC studies, either in LES or RANS: in LES the mixture fraction rms are found to be overpredicted [18] and in RANS overpredicted in the near-field and underpredicted further downstream from the jet exit [11]. 4.3. Conditionally filtered scalar dissipation rate g obtained from Figure 4 presents profiles of Njg MMC (initial values and values after reaching steady-state) and AMC at two points, Point 1 and Point 2, along with the corresponding PDF from MMC and b PDF. Physical coordinates, vale and SGS variance for Point 1 and 2 are ues for Z given in Table 1. Point 1 is located close to the inlet and represents very rich mixture with nonnegligible SGS variance, whereas Point 2 is situated in a much leaner mixture. In Fig. 4a, the characteristic bell shaped profile of the AMC model can be seen, always centred at a mixture fraction of 0.5. The MMC profiles are significantly different from that of AMC. They are not e . The symmetrical and peak values occur near Z g MMC Njg values decrease to zero for g smaller than 0.4, consistent with what is seen in the PDF. MMC does not set a predefined shape for
e (left) and Z Fig. 3. Time-averaged radial profiles of Z rms (right). Open symbols correspond to experimental data, dashed line to LES and solid line MMC.
1219
g using AMC and MMC at Fig. 4. Left: profiles of Njg two locations. Right: corresponding profiles of PDF retrieved from MMC and calculated from beta function.
Table 1 Specifications of Point 1 and 2. y/d z/d e from LES Z e from MMC Z SGS rms from LES SGS rms from MMC
Point 1
Point 2
0.6 0.6 0.864 0.76 0.094 0.086
26 2.4 0.15 0.15 0.016 0.018
g and its shape closely follows the mixing state Njg described by the PDF. According to what is observed for the SGS rms in Fig. 1, the converged MMC values are often lower than the initial data e and g (initial MMC values for Z Z 002 are equal to LES data) due to lower SGS fluctuations. At Point 1 the converged MMC g Z 002 has decreased by 10%, resulting in a much larger change in the g approximately 50% decrease. magnitude of Njg, g from MMC is very sensitive This shows that Njg to the SGS variance. The peak of the corresponding MMC PDF has shifted to richer mixtures due e moving from 0.86 to to the MMC predicted Z 0.76. At Point 1, the MMC PDF is different from that of the b form. At these very rich mixtures far away from possible chemical reactions, the accuracy of the PDF is not as crucial as for mixtures closer to stoichiometry. Experimental or DNS data would be needed to determine what the correct profile is. At Point 2 (Fig. 4b), similar to what is found at Point 1, the MMC profiles (both coincide due to negligible change in the MMC values between inig are different in tial and converged data) for Njg g shape from AMC. The MMC Njg results go down
1220
C.B. Devaud et al. / Proceedings of the Combustion Institute 34 (2013) 1213–1221
to zero for g larger than 0.2 and smaller than 0.1, i.e. where the PDF also tends to zero. The b and MMC PDF are in close agreement, the b function only showing a slightly higher peak value, consistent e is with previous findings [11]. At both points, N R1 g recovered in MMC from 0 ð NjgÞP ðgÞdg, using g and P(g) (FDF) from MMC. No scaling is Njg e at Point 1 and Point 2, in conneeded to recover N trast to previous approaches based on PDF transport equation integration [5]. 5. Conclusion Deterministic MMC has been applied in LES for the first time to a lifted flame in a vitiated coflow. MMC has some interesting features like g no modelling assumption for the FDF and Njg, which makes it very attractive as a mixing model. Some modelling and numerical challenges, in particular in the context of LES, have been examined in the present paper. This first analysis includes MMC instantaneous results for mixture fraction, g and FDF starting from frozen SGS rms, Njg LES–CMC well-established conditions. Timeaveraged statistics are also collected for mixture fraction and resolved rms. Reasonable agreement e and resolved rms is found for radial profiles of Z between MMC, LES and experimental data. Best agreement is found close to the fuel exit and some discrepancies can be seen further downstream. It is believed that the current approximation of neglecting spatial transport in one of the transverse directions is mainly responsible for the observed discrepancies. Despite the reduced dimensionality implemented in the current MMC, it is very positive to see that deterministic MMC can reproduce the main turbulent mixing characteristics. Future work will focus on performing MMC in three dimensions. The resulting increased computational cost will be monitored and minimized. The second item for future investigation will be the linear model and alternatives to generate SGS fluctuations in MMC. For example, a gradient model may be implemented [31]. Accurately reproducing the level of SGS variance in MMC proves to be challenging but a key to accurate preg and the turbulent mixing field. Furdiction of Njg g are shown to be ther, the calculated values of Njg very sensitive to the level of SGS fluctuations. Thus, it is crucial to get good predictions for SGS variance. g are significantly difThe MMC profiles of Njg ferent from those of AMC, which certainly prompts further comparison with available experimental data and investigation of their impact on the CMC calculations for predicted lift-off height and species concentrations.
Acknowledgment The authors would like to thank Dr. K. Vogiatzaki for helpful discussions on MMC, as well as Dr. A. Klimenko. Funding from FWO project G.0079.07 is gratefully acknowledged. References [1] R. Fox, Computational Models for Turbulent Reacting Flows, Cambridge University Press, 2003. [2] A. Klimenko, R. Bilger, Prog. Energy Combust. Sci. 25 (6) (1999) 595–687. [3] S. Girimaji, Phys. Fluids A 4 (11) (1992) 2529–2537. [4] E. O’Brien, T. Jiang, Phys. Fluids A 3 (12) (1991) 3121–3123. [5] A. Milford, C. Devaud, Combust. Flame 157 (2010) 1467–1483. [6] C. Devaud, R. Bilger, T. Liu, Phys. Fluids 16 (6) (2004) 2004–2011. [7] M. Cleary, J. Kent, Combust. Flame 143 (2005) 357–368. [8] S. Sreedhara, Y. Lee, K.Y. Huh, D. Ahn, Combust. Flame 152 (2007) 282–286. [9] A.Y. Klimenko, S. Pope, Phys. Fluids (2003) 1907– 1925. [10] K. Vogiatzaki, A. Kronenburg, M.J. Cleary, J.H. Kent, Proc. Combust. Inst. 32 (2009) 1679–1685. [11] K. Vogiatzaki, M.J. Cleary, A. Kronenburg, J.H. Kent, Phys. Fluids 21 (2009) 025105. [12] A. Cook, J. Riley, Phys. Fluids 6 (1994) 2868–2870. [13] A. Triantafyllidis, E. Mastorakos, R. Eggels, Combust. Flame 156 (2009) 2328–2345. [14] I. Stankovic, A. Triantafyllidis, E. Mastorakos, C. Lacor, B. Merci, Flow Turbul. Combust. 86 (2011) 689–710. [15] S. Navarro-Martinez, S. Rigopoulos, Flow Turbul. Combust. 87 (2011) 407–423. [16] A. Garmory, E. Mastorakos, Proc. Combust. Inst. 33 (2011) 1673–1680. [17] M. Cleary, A. Klimenko, J. Janicka, M. Pfitzner, Proc. Combust. Inst. 32 (2009) 1499–1507. [18] M. Cleary, A. Klimenko, Phys. Fluids 23 (2011) 115102. [19] I. Stankovic, Numerical Simulations of Hydrogen Auto-ignition in Turbulent Flows, Ph.D. thesis, Ghent University, Belgium, 2011. [20] C. Markides, E. Mastorakos, Proc. Combust. Inst. 30 (2005) 883–891. [21] R. Cabra, T. Myhrvold, J. Chen, R. Dibble, A. Karpetis, R. Barlow, Proc. Combust. Inst. 29 (2002) 1881–1888. [22] R. Gordon, A. Masri, S. Pope, G. Goldin, Combust. Theory Model. 11 (2007) 351–376. [23] R. Cao, S. Pope, A. Masri, Combust. Flame 142 (2005) 438–454. [24] K. Gkagkas, R. Lindstedt, Combust. Theory Model. 13 (2009) 607–643. [25] A.R. Masri, R. Cao, S.B. Pope, G.M. Goldin, Combust. Theory Model. 8 (2004) 1–22. [26] T. Myhrvold, I. Ertesvag, I. Gran, R. Cabra, J. Chen, Combust. Sci. Tech. 178 (2009) 1001–1030. [27] S. Patwardhan, D. Santanu, K. Lakshmisha, Proc. Combust. Inst. 32 (2009) 1708–1712. [28] C. Duwig, L. Fuchs, Combust. Sci. Tech. 180 (2008) 453–480.
C.B. Devaud et al. / Proceedings of the Combustion Institute 34 (2013) 1213–1221 [29] W. Jones, S. Navarro-Martinez, Combust. Flame 150 (2007) 170–187. [30] S. Navarro-Martinez, A. Kronenburg, Flow Turbul. Combust. 87 (2011) 377–406. [31] M. Cleary, A. Klimenko, in: T. Echekki, E. Mastorakos (Eds.), Turbulent Combustion Modelling Advances, New Trends and Perspectives, Springer, 2011.
1221
[32] A. Triantafyllidis, E. Mastorakos, Flow Turbul. Combust. 84 (2009) 481–512. [33] B. Van Leer, J. Comput. Phys. 14 (1974) 361–370. [34] J. Li, Z. Zhao, A. Kazakov, F.L. Dryer, Int. J. Chem. Kinet. 36 (2004) 566–575. [35] P. Vaishnavi, A. Kronenburg, Combust. Flame 157 (2010) 1863–1865.