A priori investigation of subgrid correlation of mixture fraction and progress variable in partially premixed flames

Subgrid correlation of mixture fraction, Z, and progress variable, c, is investigated using direct numerical dimulation (DNS) data of a hydrogen lifted jet flame. Joint subgrid behaviour of these two scalars are obtained using a Gaussian-type filter for a broad range of filter sizes. A joint probability density function (JPDF) constructed using single-snapshot DNS data is compared qualitatively with that computed using two independent β-PDFs and a copula method. Strong negative correlation observed at different streamwise locations in the flame is captured well by the copula method. The subgrid contribution to the Z–c correlation becomes important if the filter is of the size of the laminar flame thickness or larger. A priori assessment for the filtered reaction rate using the flamelet approach with independent β-PDFs and correlated JPDF is then performed. Comparison with the DNS data shows that both models provide reasonably good results for a range of filter sizes. However, the reaction rate computed using copula JPDF is found to have a better agreement with the DNS data for large filter sizes because the subgrid Z–c correlation effect is included.


Introduction
Partially premixed combustion is ubiquitous in many combustion applications where perfect premixing of fuel and oxidiser is generally difficult to achieve [1][2][3]. In typical experimental and numerical combustion studies, it is common to use a mixture fraction, Z, to describe scalar mixing and a progress variable, c, to denote the progress of chemical reactions. The partially premixed combustion mode involves flame propagation in mixtures with varying equivalence ratios, and turbulent fluctuations of the mixture fraction and the progress variable can mutually influence one another resulting in a significant cross correlation [4][5][6].
In the Reynolds-averaged Navier-Stokes (RANS) context, the effect of this correlation has been assessed previously by Ruan et al. [7] and Chen et al. [8] using a presumed joint probability density function (JPDF) approach for turbulent lifted flames. The statistical correlation was included in the JPDF using a copula method [7,9]. It was found that the effect of Z-c correlation is quite significant and it must be taken into consideration in order to predict the correct flame lift-off height. The importance of Z-c correlation and its significance for subgrid reaction rate closure in large eddy simulation (LES) are unclear and have not been investigated.
As LES directly solves large-scale unsteady fluid motions, the Favre cross correlation (or covariance) can be decomposed into its resolved and subgrid parts written as where { · } denotes a density-weighted (or Favre) time-averaging operation and [ Z c ] sgs = ( Zc − Z c ) is the unresolved subgrid covariance requiring a closure model. The tilde here represents Favre filtering. Note that the two additional residual terms involving cross correlation between the temporal and spatial averaging procedures [10] are considered to be negligible for this study as in [11]. Previous RANS studies [7,8,12,13] showed that the total time-averaged correlation, {Z c }, is quite influential in partially premixed lifted flames. This was further confirmed in a recent experimental study by Barlow et al. [5] showing a strong Z-c correlation in a partially premixed piloted jet flame. However, in LES this correlation is partly resolved by the numerical grid and it is of interest to understand the effect of the remaining unresolved part, [ Z c ] sgs , in Equation (1) and its behaviour with the filter size. In many existing LES models for partially premixed combustion such as the flamelet-based approaches [14][15][16][17][18][19][20][21] and conditional moment closure (CMC) with doubleconditioning [22,23], a presumed JPDF is used to account for the subgrid fluctuations of mixture fraction and progress variable. This joint PDF is computed locally using the transported scalars, usually the Favre-filtered values and subgrid variances of Z and c. The filtered thermo-chemical quantities can then be obtained from the flamelet manifolds using where ξ and ζ are the sample space variables for Z and c, respectively. The flamelet value for Φ is F Φ and the joint PDF is generally taken to be the product of two marginal PDFs [14,21,24,25] as P (ξ, ζ ) ≈ P (ξ ; Z, Z 2 sgs ) P (ζ ; c, c 2 sgs ) after dropping t and x for conciseness. This assumes statistical independence between the subgrid fluctuations of Z and c, i.e. [ Z c ] sgs = 0, and its validity is an open question. However, the unresolved scales of Z and c interact at the subgrid level and this interaction, depending on the filter size, can have mutual influence leading to their correlation which can affect the filtered reaction rate. To account for this Z-c correlation, the subgrid covariance can be considered while modelling the JPDF as P (ξ, ζ ) = P (ξ, ζ ; Z, Z 2 sgs , c, c 2 sgs , [ Z c ] sgs ) using the copula method [7][8][9]. However, the superiority of this approach over the independent JPDF is yet to be assessed for LES since modelling [ Z c ] sgs would involve further complexity and additional computational cost.
It is worth noting that the function P (ξ, ζ ; x, t) in Equation (2), referred here as subgrid JPDF, takes into account the subgrid-scale fluctuations of Z and c, and it can be obtained by collecting the samples in a given subgrid space (with spatial filtering) at a given time over many ensembles having the same resolved fields, as first performed experimentally by Tong [26]. This differs from the filtered density function (FDF), which could be obtained formally by filtering fine-grain function for a single ensemble (or realisation) as shown in past studies [27][28][29][30][31][32][33]. This concept was introduced by Pope [27] and the FDF was shown to have all the properties of PDF when the filtering kernel is positive by Gao and O'Brien [28]. However, the FDF has some randomness associated with it and one needs to be cautious while evaluating the subgrid statistics using FDF obtained from single-realisation experimental or numerical data [33,34]. Moreover, the filter width, , also influences the subgrid fluctuations and thus the PDF depends on further to the five quantities noted above. This subgrid PDF becomes the one-point one-time PDF in the classical sense when → 0 [27].
Recently, direct numerical simulation (DNS) has become a useful tool for model assessment through a priori analysis. In the context of partially premixed flames, it has been utilised to study the effect of mixture fraction gradients on flame propagation speed in turbulent mixing layers [35], droplet combustion [36], and slot-jet flames with varying equivalence ratio inlets [37]. Indeed, these studies have provided many useful insights for partially premixed flames but the joint behaviour of Z and c were not examined from a modelling perspective. In relevance to LES, the DNS data can be filtered in physical space to obtain the subgrid information at each time instant. It helps to understand the joint behaviour of scalars at the subgrid level and thus the validity of combustion models for LES can be assessed. This technique has been used to test various filtered reaction rate closure models for purely premixed or stratified flames [38][39][40] but it is scarce for partially premixed combustion.
The objective of this paper is to investigate the relative significance of the subgrid Z-c correlation, [ Z c ] sgs , at different filter sizes using the DNS data of a partially premixed lifted hydrogen jet flame [41,42]. It is of interest to assess the validity of the commonly used statistical independence assumption for the subgrid JPDF. The correlated JPDF computed using the copula approach is examined for this. The filtered reaction rate of progress variable is computed using a partially premixed flamelet model [24] with the JPDF obtained using the independent and correlated presumed-shapes. These modelled reaction rates are then compared with the filtered reaction rate from the DNS to elucidate the importance of [ Z c ] sgs . This paper is organised as follows. Firstly, a brief description of the DNS dataset is given along with the post-processing methodology used here. The joint statistics of Z and c are then discussed for a particular time instant along with time-averaged results. Subsequently, the filtered reaction rate closures are assessed through a priori analysis. Both conditional averages based on single snapshot data and time averages of reaction rate are investigated. Finally, conclusions are summarised.

DNS data and post-processing methodology
The present DNS data have been discussed extensively in [4,[41][42][43] and used in a previous RANS modelling study [7]. Thus, only a brief description is presented here. The configuration comprises a pure hydrogen jet issuing from a round nozzle into quiescent air at 280 K and atmospheric pressure, establishing a turbulent lifted flame in the downstream. The nozzle diameter, D, is 2 mm and the bulk mean jet exit velocity, U j , is 680 m/s, corresponding to a jet Reynolds number of about 13,600 and a Mach number of 0.52. The computational domain starts from 2D below the jet exit plane to 20D in the streamwise direction, and extends to ±12.5D in the transverse directions. A uniform grid spacing of 50 μm is used in the entire flame region, which is about 1/10 of the laminar premixed flame thermal thickness for a hydrogen-air mixture, [δ 0 L ] st = 0.44 mm. This grid resolution resolves the flame and scalar statistics quite well, as shown in previous studies [4,41], and thus the DNS data used for the purpose of the present study is considered to have been validated. Instantaneous DNS solutions are saved for 1000 time snapshots spanning a physical time of 1.2 ms. This covers about six characteristic laminar flame timescales, τ L = δ 0 L /S 0 L , based on a premixed stoichiometric H 2 -air mixture. Detailed transport properties and chemical kinetics are considered in the DNS using a mechanism involving 9 species and 17 reactions [44]. Figure 1 shows the instantaneous and filtered progress variable fields along with the stoichiometric mixture fraction contour. The middle-plane snapshot showing the flame evolution in the streamwise direction is presented in Figure 1(a) (the right half is the filtered version of the left half) and cross-plane contours are shown in Figures 1(b) and 1(c) for unfiltered and filtered c fields, respectively. The progress variable is defined using the H 2 O mass fraction normalised by its equilibrium value for the local mixture fraction as in [4,7,43]: The mixture fraction is calculated using Bilger's [45] formulation written as where the subscripts 1 and 2 denote fuel and oxidiser streams, respectively, and W i is the molar mass for element i. The stoichiometric value is Z st = 0.03. A typical laminar flame thickness mentioned earlier, [δ 0 L ] st , is taken as a reference length hereafter and the normalised filter size is calculated as where F is the filter width. A filter size of + = 1 is used in Figure 1 to obtain the filtered field of the progress variable, c, and the stoichiometric mixture fraction, Z st . The filter size is marked in Figure 1(b) for illustration. Three streamwise locations, x/D = 6.5, 10, and 15, marked as A, B, and C, respectively, are chosen to cover the streamwise variation of the reaction zone for later investigation. It is shown in Figure 1 that the filtering operation smooths out the sharp scalar gradients and fine turbulent structures in the flame as one would expect, and the methodology used for this operation is described next. The various quantities of interest from the DNS data are filtered with density weighting and the filtering is performed within a 2 F × 2 F × 2 F filter box, centred at a filtering point, x. For example, the Favre-filtered mixture fraction is computed as where ρ and ρ are the unfiltered and filtered mixture densities, respectively. The coordinate vector, x , corresponds to the sample point in the filter sub-space. The Gaussian filter kernel, G , is given as [46] The Favre subgrid variance of Z can then be obtained using The values of c and c 2 sgs are computed in a similar manner. The Favre subgrid covariance, which signifies the interaction of the mixture fraction and progress variable fluctuations, is Note that the scalar variance and covariance are often defined using different expressions while deriving their transport equations required for LES to avoid the additional correlation terms, for instance, [ Z c ] sgs = Zc − Z c. This expression is found to provide very similar results as Equation (9) and the difference observed is less than 1% in this a priori study. These spatial filtering techniques are applied for every two instantaneous snapshots of the DNS data and, in total, 500 instants are processed to obtain time-averaged statistics. This sampling period, τ = 1.2 ms, is about 20 flow-through times, which is computed as L/U j , where L is the computational domain length. A simple time averaging procedure is applied over the entire sampling period, τ , for the filtered quantities. For example, the mean mixture fraction is calculated as The corresponding total variance is obtained using with < Z 2 res > being the resolved variance. The total covariance is computed in a similar manner For the subgrid JPDF analysis, spatial samples are collected inside the filter at a given location from a single time instant of the DNS data. Different filter sizes are used to investigate the sensitivity of the subgrid scalar statistics and these details are summarised in Table 1. Note that, due to the high computational cost, the time-averaging procedure is performed only for a number of representative filter sizes. For small filter sizes, the number of sample points is insufficient to show small-scale fluctuations in space. This could be overcome by combining the samples from multiple snapshots within a time window smaller than the typical LES time-step size, which is typically two or three orders of magnitude larger than the time step used for the DNS and a typical LES time step is of the order of 10 −6 s. Unfortunately, the DNS solution was saved every 1.2 × 10 −6 s and thus the saved data do not contain small-scale fluctuations in time information. This common practice for DNS is due to data storage limitation, but it would be useful to have time-resolved data (sets of a few time steps with each set uncorrelated in time and space) archived in future DNS studies.
Using the first and second moments of a given random variable, a common approach to model its PDF is to assume a β-distribution shape [47]. For instance, the Favre subgrid PDF of the mixture fraction is calculated as where is the gamma-function. A similar procedure is used to obtain P β (ζ ; c, c 2 sgs ) for the progress variable. Assuming statistical independence between Z and c, their Favre subgrid joint PDF is simply the product of these two marginal PDFs: To construct the correlated joint PDF, the subgrid covariance, [ Z c ] sgs , is used in the copula method to couple the univariate marginal distributions, P β (ξ ) and P β (ζ ). In a brief description, the correlated JPDF for non-zero values of [ Z c ] sgs is calculated as [7,9] with where C β is the cumulative distribution function (CDF) for the β-distribution, and θ is the odds ratio calculated using [ Z c ] sgs and two sets of independent random variables. Elaborate details of this method are described in [7] and [9].
To visualise the effect of Z-c correlation on the filtered reaction rate in the look-up table approach, Figure 2 shows typical variation of the filtered reaction rate as a function of Z and c. Normalised variances, g Z = Z 2 sgs / ( Z(1 − Z )) and g c = c 2 sgs / ( c (1 − c )), are used as control parameters, in addition to Z, c, and g Zc , for the flamelet library [8] and typical values of g Z = 0.06 and g c = 0.45 are used here for illustration purpose. The correlation coefficient, g Zc , is defined as By definition, the value of g Zc varies from −1 to 1 and g Zc = 0 indicates no correlation. It is seen that both positive and negative correlations result in a higher peak reaction rate compared to the uncorrelated value. Another correlation effect is that the reaction rate contour becomes more asymmetric having a negative slope when the correlation is negative as shown in Figure 2(a) and vice versa in Figure 2(c). These effects may be reflected in the LES through the source terms of appropriate transport equations leading to different flame behaviour when the correlation is included in the modelling.

Joint subgrid Z -c statistics 3.1. Analysis of single snapshot data
The three streamwise locations marked in Figure 1 are chosen to study the joint subgrid statistics. The results in this subsection are shown for a representative time instant. It is well known that strong partial premixing exists in the lifted flame base region forming a tripleflame configuration [48,49]. This was also discussed in the original DNS investigation [41] of the present flame. Thus, for this study, subgrid statistics are collected for three different radial positions at the streamwise location A, which is close to the flame base. These three positions are designated as A-l, A-st, and A-r, corresponding to lean, stoichiometric and rich mixture fractions of about Z = 0.015, 0.03 and 0.07, respectively. For the downstream locations B and C, as the flame is burning mainly in rich mixtures [41], only filtering points with about Z = 0.1 are considered. All these locations are selected with c ≈ 0.6 for a direct comparison. Figure 3 shows the subgrid Favre JPDFs of Z and c obtained directly from the DNS data and those computed using the two independent β-PDFs and the copula method. The results are shown for five filtering points, A-l, A-st, A-r, B and C marked in Figure 1. The filter size used is + = 1. Note that density-weighted JPDF from the DNS is used here because typically the Favre-filtered quantities such as Z and c are transported in LES of reacting flows. This JPDF is obtained using where the unweighted JPDF, f Z,c (x , t), is given by (19) for x ∈ [x − F , x + F ] and P is the probability function. Thus, the quantity f Z,c (x, t) is the fraction of the fluid mass at t and around x weighted by G having Z and c in the infinitesimal ranges noted in Equation (19) given above [27]. The first and second moment values obtained from the DNS data are listed in the first frame of each row in Figure 3 and the notation for subgrid quantities, "sgs", is omitted for clarity. The correlation coefficient is also given in the copula JPDF frames along with the streamwise and radial locations of the filtering point. The number of sampling bins used is 200 in both Z and c spaces and, because single snapshot is used, it is not possible to construct the averaged subgrid JPDF. However, this subgrid quantity constructed from the DNS data is adequate to show the subgrid correlation between Z and c existing for all five probed locations. The correlation is positive for the location A-l with g Zc = 0.84. For the stoichiometric (A-st) and rich (A-r) cases, the correlation becomes negative but with a similar magnitude of g Zc being 0.89 and 0.87 respectively. This is consistent with the previous finding based on RANS methodologies [4,7,8], where the change of sign was discussed in detail on a physical basis and thus is not repeated here. This correlated behaviour and its variation with mixture fraction at the flame base is well captured by the copula approach, whereas the independent β-PDFs fail to predict both the shape and the peak. As one moves downstream, stronger negative correlation is observed with the g Zc value close to −1. As a result, the copula JPDF shapes agree with the DNS data. However, the peak values are overpredicted and concentrated in narrower regions compared to the DNS results. These results suggest that the commonly assumed statistical independence of subgrid Z and c fluctuations is questionable for the chosen filter size of + = 1, which is close to the laminar flame thickness.
Before investigating the influence of + on f Z,c (x, t), it is worth making the following remarks. The analysis discussed above was also performed using DNS snapshots separated by an interval larger than the typical turbulence integral and flame timescales to ensure that these snapshots are uncorrelated. Results similar to those in Figure 3 were observed for the same conditions (axial position, Z and c values) and thus they are not shown here. This suggests that the results obtained using a single snapshot data is sufficient to show typical Z-c correlation behaviour of this subgrid filtered PDF. Strictly, one must consider a large number of DNS runs to construct f Z,c (x, t) which can then can be averaged to get the subgrid PDF. However, this would be a very costly exercise. Alternatively, one could consider sets of snapshots from a very long DNS run with the sets separated by the largest characteristic timescale (the large-scale convection time, as suggested by an anonymous reviewer) of the flow. This is not possible with the DNS data used for this study which has a limited total duration. However, the correlation between Z and c is driven by combustion although turbulence produces Z and c . Combustion timescales are typically shorter compared to those for the large flow scales and thus some insights into the subgrid Z-c correlation can be gained by analysing DNS snapshots separated by a few flame timescales, which has been done as noted above. Thus, single snapshot data is used here to show the correlation and to reflect on the flamelet modelling discussed in Section 4.
To investigate the dependence of Z-c correlation on the filter size, Figure 4 shows the JPDFs for + = 0.6, 1.5, and 3 at the filtering point A-st. The results for + = 1 are shown earlier in the second row of Figure 3. As one would expect, the spread of the JPDF in both Z and c spaces increases as the filter size increases. This is because when the filter is large compared to the flame thickness the scalar fluctuations become less resolved by the numerical grid. A significant increase of about four times is observed for Z 2 sgs , c 2 sgs and the magnitude of [ Z c ] sgs moving from + = 0.6 to 3. The magnitude of g Zc , however, drops from 0.9 to 0.68 and it is due to the faster increase of the denominator than the numerator in Equation (17). Although it seems that the value of [ Z c ] sgs compared to Z 2 sgs and c 2 sgs is smaller for larger + , the Z-c correlation effect on combustion is more important as the filter size increases because the subgrid scalar fluctuations and their interaction become more influential and they require closure models. It can be seen in Figure 4 that the JPDF computed using the copula method captures the Z-c correlation quite well, giving a closer agreement with the DNS result compared to that computed using two independent β-PDFs for a range of filter sizes. To visualise the spatial distribution of the Z-c correlation, Figure 5 shows the instantaneous g Zc field along with Z contours for three different filter sizes in the mid-plane. It is more evident in this figure that positive Z-c correlation mainly appears in very lean mixtures as observed earlier in Figure 3. However, the negative correlation is quite strong along the Z st iso-line at the flame base and then extends to rich mixtures in the downstream. At the flame base, the magnitude of peak positive g Zc increases with the filter size while it slightly decreases for the negative one around stoichiometry due to the increase of Z 2 sgs and c 2 sgs . This is in line with the JPDFs shown in Figure 4, where the magnitude of g Zc is found to be smaller for larger filter sizes.

Significance of subgrid covariance
Since large-scale scalar fluctuations are partly resolved in LES, only the subgrid scale Z-c correlation is of interest for modelling. Figure 6 presents radial profiles of the resolved covariance, [ Z c ] res , and time-averaged subgrid covariance, [ Z c ] sgs for filter sizes of + = 0.5, 1, and 1.5 at three streamwise locations x/D = 6, 10, and 15. These secondorder statistics are obtained in a similar manner as in Equation (11). If we consider the subgrid-to-total ratio for the covariance, R = [ Z c ] sgs / Z c , it can be seen that the subgrid contribution of covariance is negligible with R < 0.01 when + = 0.5 is applied for all streamwise locations. For a larger filter size of + = 1, the value of R at the peaking radial location becomes about 0.25, 0.17, 0.09 for x/D = 6, 10, and 15, respectively. These numbers further increase to 0.5, 0.29, and 0.13 for the + = 2 case. The relative contribution of the subgrid Z-c correlation increases with the filter size as one would expect and is more significant in the upstream regions possibly due to higher level of turbulence. Thus, for an LES grid with a filter size similar to or larger than the typical laminar flame thickness, which is common in LES practice, the subgrid Z-c correlation can have a substantial effect on time-averaged statistics. This is consistent with the importance of Z-c correlation observed in previous RANS studies [7][8][9]12,13]. This is specifically so in flame stabilisation regions as can be seen in Figure 6 for x/D = 6.

Effect on filtered reaction rate modelling
Reaction rate closure for the progress variable transport equation is a central aspect of turbulent combustion modelling. For fully premixed flames, it is straightforward to derive this equation from the species conservation equation based on the definition of c given in Equation (3) because the equilibrium mass fraction is a constant. However, for partially premixed flames, this equilibrium value varies with mixture fraction and thus one needs to be cautious while deriving the convection-diffusion-reaction form as has been elaborated by Bray et al. [50]. The resulting reaction rate term after filtering is written as [7,8,19] where the asterisk * appearing inω * c denotes the partially premixed reaction rate, andω c signifies the contribution of premixed mode combustion with mixture fraction stratifications whereasω np denotes the contributions from non-premixed mode. Note that here the cross dissipation contribution is considered to be small and thus neglected following previous studies [7,24]. By considering these multi-mode combustion effects, this modelling approach has been shown to be adequate to capture the triple flame structures which are present in the flame base region of lifted jet flames using both RANS [7,8,13] and LES [24].
As in common LES flamelet approaches, the premixed term is modelled using the form of Equation (2), written aṡ where the flamelet reaction rate,ω c , is obtained from unstrained laminar planar premixed H 2 -air flame calculations with the number of equivalence ratios varying from the lean to rich flammability limit. The second term in Equation (20),ω np , denoting contributions of non-premixed mode combustion is modelled as [7,24] ω np ρ c χ Z where χ Z is the filtered scalar dissipation rate including both resolved and subgrid parts. In this section, the performance of this model forω * c with the JPDF obtained from modelling in Equations (14) and (15) is investigated by comparing with the filtered reaction rate computed from the DNS data. Figure 7 compares the instantaneous spatial variations of filtered reaction rate contours computed using the DNS data and the flamelet model with the two presumed JPDFs. Four rows correspond to the filter sizes of + = 0.6, 1, 2, and 3, respectively. It can be seen that the overall reaction zone shape is captured well by the independent β-PDFs and copula JPDF models. It is also observed that, as the filter size increases, the flame becomes less wrinkled and the maximum reaction rate decreases at the flame base. A considerable overestimation is observed in the downstream, x/D > 7, for the independent-β and copula models. This overestimation was also found in [7] and may be due to the overprediction of laminar flame reaction rate in the premixed flamelets of a rich mixture. It seems in Figure 7 that no apparent difference is observed between the reaction rate contours computed using the independent β-PDFs and copula JPDF, even though the latter showed a clear improvement when compared to the DNS JPDF earlier in Figures 3 and 4. To investigate this further, more quantitative comparison of the filtered reaction rate is required as will be seen later in this section. Figure 8 presents the premixed and non-premixed mode contributions toω * c computed using Equations (21) and (22) with the copula method. The contour ofω c is very similar to that ofω * c for all filter sizes used suggesting that premixed mode combustion is dominant. The contribution from non-premixed mode is negative as has been observed in previous studies [4,7,8] and this effect is more evident in the downstream of the flame base which is around x/D = 6.5. The value ofω np is an order of magnitude smaller than that ofω c in the flame base region, whereas they become comparable in the downstream. This is consistent with previous findings [1,48] arguing that partial premixing plays an important role in the lifted jet flame stabilisation mechanism at the flame base and that the downstream combustion is diffusion controlled.
Scalar dissipation rate (SDR) of mixture fraction is a key parameter in diffusion flames and its effect is included in Equation (22) through the filtered SDR, χ Z = ( D Z |∇ Z| 2 + χ Z,sgs ), where D Z is the filtered diffusivity of the mixture fraction. Figure 9 compares the resolved and subgrid parts of χ Z at streamwise locations of x/D = 6, 10, and 15 for filter sizes ranging from + = 0.6 to 3. Here the subgrid SDR is obtained by subtracting the resolved part from the filtered total SDR computed from the DNS data. It is seen that for x/D = 6 the SDR is largely under-resolved due to the high turbulence in the jet near-field, and the values of SDR decrease as one moves downstream. For the downstream locations, the resolved part is substantially larger than the subgrid one for filter sizes + ≤ 1 and vice versa. This suggests that both the resolved and subgrid parts need to be included in the modelling ofω np in Equation (22).
In order to isolate the mixture stratification effect on the reaction rate, Figure 10 compares the mixture fraction conditioned average of the filtered reaction rate between the DNS data and modelled values. The samples used are obtained from cross-sectional slices of single-time snapshot data at different streamwise locations. Reasonably good agreement is obtained for both models using the independent β-PDFs and copula JPDF for the four filter sizes and three streamwise locations considered. At the flame base location, x/D = 6, the two presumed-shape models show quite good agreement with the DNS data for the conditional average value ofω * c . Moving to the downstream locations, at x/D = 10 and 15 the previously observed overprediction of filtered reaction rate in rich mixtures ( Z > 0.03) is more evident here for both models and this overprediction increases with the filter size. A quite substantial improvement is seen for the copula JPDF compared to the independent one and it increases from about 10% for + = 1 to 30% for + = 3. This difference can be quite influential during transient processes such as the spark ignition of jet flames -in a previous study [24] using the independent β-PDFs model, the flame leading edge propagation speed was overestimated. Figure 11 compares the time-averaged reaction rate, ω * c , computed using the independent β-PDFs and copula approaches, with the filtered DNS data for + = 1, 2, and 3. The results for + = 0.6 are not shown because they are very similar to those for + = 1 as in Figure 10. The modelled results presented here were obtained by applying the reaction rate models for 500 DNS snapshots individually and the time-averaged values were then averaged along the azimuthal direction. It is shown that the reaction rate is predicted very well by both models at the near-field location, x/D = 6, for all three filter sizes applied. This is consistent with the observations in Figure 10. Again, the reaction rate is overpredicted in the downstream locations for large filter sizes and the copula approach slightly improves the agreement with the DNS data especially for large filter sizes. This is because, as a flame becomes largely under-resolved, i.e. + > 1, the subgrid interaction between mixture fraction and progress variable fluctuations increases resulting in a stronger statistical correlation between these two variables. It was also shown in a previous LES study [24] that the flame leading edge propagation speed in the spark ignition process of a methane-air jet flame was overpredicted by the independent β-PDF approach. Therefore, it may be worthwhile to include this correlation in the modelling approach. However, this would require modelling of the subgrid covariance, which is not straightforward.
In previous RANS studies [7,8,13], a transport equation for the covariance was used with a simple linear relaxation model for the cross dissipation rate closure which may not be appropriate for LES. This is because the cross dissipation rate is partly resolved by the numerical grid and the subgrid part is more likely to be dominated by combustion, not turbulence effects. Therefore, further investigation is needed for the modelling of covariance and cross dissipation rate. One way to do this is to analyse the balance of the various terms in the covariance transport equation using DNS data, which will be explored in a future work.

Conclusions
In the present work, subgrid correlation of mixture fraction and progress variable is investigated using DNS data of a lifted turbulent hydrogen jet flame. The DNS data is processed using a low-pass Gaussian filter with a number of filter sizes covering the commonly used range in LES. Both single and multiple snapshots of data are processed to study the subgrid correlation. The probability distribution constructed using a single snapshot is very similar to that obtained using multiple snapshots separated by a few (stoichiometric) flame timescales. This is because the subgrid Z-c correlation is driven predominantly by combustion. This correlation is observed to be pronounced for various streamwise positions using filter sizes ranging from 0.6 to 3 times a typical laminar flame thickness. Although one strictly requires a large ensemble (many DNS runs or a single run for a very long duration with data saved for many large-scale turbulence and convection timescales) to draw statistical inferences required for subgrid PDF modelling, some useful insights are obtained by using limited available DNS data. The joint subgrid PDF computed using a copula approach agrees quite well with the DNS result compared to that for the commonly used two β-PDFs method, which assumes a statistical independence between Z and c. The results of time-averaged covariance show that the Z-c correlation is predominantly negative and its subgrid component becomes influential when the filter size is similar to or larger than the flame thickness. It would be more useful to construct the covariance conditioned on the filtered mixture fraction and progress variable for a statistical inference from this correlation, and this conditional averaging can only be meaningful using a large ensemble, which can be addressed in the future.
A priori analysis of the filtered reaction rate model is conducted for the independent β-PDFs and copula approaches, and it is found that both can provide reasonably good results compared to the DNS data. The effect of subgrid Z-c correlation is more evident in the instantaneous aspect and the copula method shows a substantial improvement up to 30% for large filter sizes. This is because the joint behaviour of Z and c at the subgrid level are captured well. The improvement is found to be less obvious for the time-averaged values of reaction rate. However, it should be noted that these time-averaged results are based on DNS snapshots that are post-processed independently, and hence the transient effects of the subgrid correlation cannot be seen in such an a priori study. These transient effects can play quite important roles in a posteriori LES, which will be addressed in future studies.