A DNS Study of Sensitivity of Scaling Exponents for Premixed Turbulent Consumption Velocity to Transient Effects

3D Direct Numerical Simulations of propagation of a single-reaction wave in forced, statistically stationary, homogeneous, isotropic, and constant-density turbulence, which is not affected by the wave, are performed in order to investigate the influence of the wave development on scaling (power) exponents for the turbulent consumption velocity UT as a function of the rms turbulent velocity u′, laminar wave speed SL, and a ratio L11/δF of the longitudinal turbulence length scale L11 to the laminar wave thickness δF. Fifteen cases characterized by u′/SL = 0.5,1.0,2.0,5.0, or 10.0 and L11/δF = 2.1, 3.7, or 6.7 are studied. Obtained results show that, while UT is well and unambiguously defined in the considered simplest case of a statistically 1D planar turbulent reaction wave, the wave development can significantly change the scaling exponents. Moreover, the scaling exponents depend on a method used to compare values of UT, i.e., the scaling exponents found by processing the DNS data obtained at the same normalized wave-development time may be substantially different from the scaling exponents found by processing the DNS data obtained at the same normalized wave size. These results imply that the scaling exponents obtained from premixed turbulent flames of different configurations may be different not only due to the well-known effects of the mean-flame-brush curvature and the mean flow non-uniformities, but also due to the flame development, even if the different flames are at the same stage of their development. The emphasized transient effects can, at least in part, explain significant scatter of the scaling exponents obtained by various research groups in different experiments, thus, implying that the scatter in itself is not sufficient to reject the notion of turbulent burning velocity.

D molecular diffusivity Da = τ T /τ F Damköhler number H height of a Bunsen flame Ka = τ F /τ η Karlovitz number L an integral length scale of turbulence L 11 longitudinal turbulence length scale Le Lewis number l 0 = /4 forcing scale N number of grid points in a transverse direction Pe= u L 11 / (S L δ F ) turbulent Péclet number Pr= ν/a Prandtl number p pressure q r scaling (power) exponent for a ratio of turbulence length scale to laminar flame thickness q s scaling (power) exponent for laminar flame speed q v scaling (power) exponent for rms turbulent velocity R radius of the nozzle of a burner R f mean radius of reaction-wave kernel Re T = u L 11 /ν u turbulent Reynolds number r f = R f /L 11 normalized mean radius of reaction-wave kernel S ij rate-of-strain tensor S L laminar flame speed S T turbulent flame (displacement) speed Sc= ν/D Schmidt number T Temperature t Time t * time instant for starting sampling statistics t d wave-development time U mean flow velocity at the nozzle of a burner U T turbulent burning (consumption) velocity u = {u 1 , u 2 , u 3 } velocity vector u rms turbulent velocity u t = U T /u normalized turrbulent consumption velocity W rate of product creation x = {x 1 , x 2 , x 3 } = {x, y, z} spatial coordinates x axis normal to mean flame brush Z distance from the nozzle of a burner Ze Zel'dovich number Greek symbols x numerical resolution δ F = D/S L laminar flame thickness δ T mean turbulent flame brush thickness ε rate of dissipation of turbulent kinetic energy θ = t d /τ T normalized wave-development time θ m maximal normalized wave-development time η Kolmogorov length scale width of computational domain ν kinematic viscosity ρ density σ = ρ u /ρ b density ratio τ F = δ F /S L laminar flame time scale τ T = L 11 /u turbulent time scale τ η = (ν u /ε) 1/2 Kolmogorov time scale

Subscripts and superscripts Q (x)
Reynolds-averaged value of a quantity Q, averaged over transverse plane and timē Q (x, t) value of a quantity Q, averaged over ensemble of transient fieldŝ Q (x, t) value of a quantity Q, averaged over transverse plane Q time and volume-averaged value of a quantity Q 0 turbulence characteristics at t =0 b burned

Introduction
When investigating premixed turbulent combustion, burning rate is commonly characterized using either turbulent flame (or displacement) speed S T , i.e., the speed of a mean flame surface with respect to the local mean flow, or burning (or consumption) velocity U T , i.e., the bulk mass rate of consumption of a reactant (or creation of a product), normalized using a mean flame surface area and the partial density of the reactant in the fresh mixture (or the partial density of the product in burned gas, respectively). Accordingly, as reviewed elsewhere [1][2][3], S T and U T were in the focus of experimental research into premixed combustion for many years. In particular, over the past decade, such measurements were conducted by various research groups, e.g., [4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21]. However, values of U T and S T obtained from different experiments under similar conditions, i.e., comparable rms turbulent velocities u , close unburned gas temperatures, the same pressure, the same fuel, and the same equivalence ratio, are strongly scattered, e.g., see Fig. 4.11 in Ref. [22]. Furthermore, even scaling (power) exponents q m in fits m to various experimental databases are scattered [3]. In these two expressions, M ≥ 3, Q 1 = u , Q 2 = S L , and other Q m are substituted with a ratio L/δ F of an integral length scale L of turbulence to the laminar flame thickness δ F = a/S L , turbulent Reynolds number Re t = u L/ν u , Damköhler number Da = τ T /τ F , or Karlovitz number Ka = τ F /τ η . Factors b 1 and b 2 are typically constant, but, in some fits, they depend on the density ratio σ = ρ u /ρ b , Prandtl number P r = ν/a, Schmidt number Sc = ν/D, or Lewis number Le = a/D. Here, τ T = L/u , τ η = (ν u /ε) 1/2 , and τ F = δ F /S L are eddyturn-over, Kolmogorov, and laminar-flame time scales, D is the molecular diffusivity of a deficient reactant in a mixture, ν and a are the kinematic viscosity and heat diffusivity of the mixture, subscripts u and b designate fresh reactants and equilibrium products, respectively, ε = 2νS ij S ij is the dissipation rate, S ij = 0.5 ∂u i /∂x j + ∂u j /∂x i is the rate-of-strain tensor, u i and x i are components of the velocity vector u and spatial coordinates x, respectively,Q is the Reynolds-averaged value of a quantity Q, and the summation convention applies for repeated indexes.
It is worth stressing that not only results of earlier measurements, processed in Ref. [3], but also recent experimental data on S T and U T show significant scatter of the aforementioned scaling exponents. For instance, the scaling exponent q v for S T or U T vs. the rms velocity u , e.g., S T ∝ u q v , was reported to be (i) less than unity [5, 7-14, 16, 19-21], e.g., q v = 0.49 [16], q v = 0.55 [21], q v = 0.63 [8], (ii) equal to unity [6], or (iii) even equal to two [17]. The fitted values of the scaling exponent q s for S T or U T vs. S L range from 0.37 [8] to 0.74 [10]. The scatter is much more pronounced for the scaling exponent q r for S T or U T vs. a ratio of L/δ L , with both negative, e.g., q r = −0.37 [8], and positive, e.g., q r = 0.25 [14], q r = 0.5 [13], or even q r = 1.35 [10], being reported.
There are two major, but fundamentally opposite standpoints regarding the significant scatter of experimental data on U T , S T , and the scaling exponents obtained from various flames using various methods. On the one hand, based on this scatter, the notions of U T and S T are sometimes put into question in the sense that they may not be extrapolated beyond a particular flame configuration used to evaluate them. Then, strictly speaking, research into turbulent burning velocity and flame speed is of minor fundamental value. On the other hand, while the scatter of measured U T , S T , and their scaling exponents seems to be a sufficient reason for expressing concern over the notions of these quantities, such a standpoint is not widely accepted and evaluation of U T and S T is still in the focus of experimental research into premixed combustion, e.g., see [4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21] and references quoted therein. Moreover, U T and S T are widely used in numerical simulations of turbulent flames, as reviewed elsewhere [3,[22][23][24][25]. To warrant such experimental and numerical investigations, it is necessary to assure that the aforementioned scatter does not prove that the notion of U T or S T should be limited to a particular flame configuration. Accordingly, it is necessary to reveal effects that (i) can cause the well-documented significant scatter of experimental data on U T , S T , and even the scaling exponents, but (ii) are consistent with a standpoint that the notion of U T or S T is of a fundamental value.
As far as the values of U T and S T are concerned, three such effects are well known. As shown [26][27][28][29][30][31] and reviewed [3,22,32] elsewhere, the quantitative scatter of available experimental data on U T and S T results, at least in part, from sensitivity of values of turbulent burning velocities and flame speeds to methods used to analyze raw experimental data and, in particular, to the choice of a mean flame surface within a thick mean flame brush. More specifically, first, U T is sensitive to such a choice, because (i) areas of different mean flame surfaces can significantly differ from each other if the mean flame brush is thick and curved [28,30,31] and (ii) U T is inversely proportional to the area used to normalize the bulk burning rate. Second, speeds of different mean flame surfaces can be substantially different due to non-uniformities of the mean flow of unburned gas within the thick flame brush, with such effects being well known in spatially diverging flows [26,27]. Third, S T characterizing different mean flame surfaces can be significantly different if the mean flame brush thickness δ T grows with time or distance from a flame-holder [28], with such a growth of δ T being well documented in a number of experiments with various premixed flames, as reviewed elsewhere [3,22,32], see also recent papers [10,16,33,34]. While the focus of studies [3,22,[26][27][28][29][30][31][32] was placed on the values of U T and S T , rather than the scaling exponents, a recent numerical work [35] has shown that the three aforementioned well-known effects can even cause substantial sensitivity of the scaling exponents q v , q s , and q r to methods used to evaluate U T or S T by processing raw experimental data.
These three effects are relevant, respectively, to (i) U T if the mean flame brush is curved, (ii) S T if the mean flow of fresh reactants is not spatially uniform, e.g., the mean flow diverges or converges, and (iii) S T if δ T grows. However, none of these three well-known effects is relevant to turbulent burning velocity obtained from a statistically 1D planar turbulent premixed flame. In this simplest case, U T is equal to properly normalized mean rate of a reactant consumption (or product creation), integrated along a straight line normal to a mean flame surface, with all mean flame surfaces being parallel to each other and having the same area. Therefore, U T is well and unambiguously defined. It is sensitive neither to the choice of a mean flame surface nor to the growth of the mean flame brush thickness, with the mean flow of unburned reactants being spatially uniform in the considered case. Nevertheless, even in this simplest case associated with negligible role played by the aforementioned three well-known effects, there is another effect, which can also contribute to significant scatter of the scaling exponents for U T . The major goal of present work is to emphasize substantial sensitivity of the scaling exponents to that effect.
More specifically, the paper aims at analyzing Direct Numerical Simulation (DNS) data obtained for various ratios of u /S L and L/δ F in the statistically 1D planar case in order to show that transient effects can yield a significant scatter of the scaling exponents q v , q s , and q r for the turbulent consumption velocity even in such a case. As noted above, consideration of the statistically 1D planar case offers an opportunity to place the focus of the study on the transient effects by eliminating the three other relevant well-known effects discussed earlier. Moreover, the paper aims at demonstrating that the transient effects can manifest themselves in apparent sensitivity of the scaling exponents to flame configuration and measurement method. Other related research issues such as comparison of computed scaling exponents with theoretical or experimental results are beyond the major scope of the paper.
In the next section, statement of the problem and DNS attributes are summarized. In the third section, results are discussed, followed by conclusions.

Statement of the problem
In order to place the focus of the study on the aforementioned transient effects and to isolate them from other phenomena as much as possible, the simplest relevant problem is addressed, i.e., we consider propagation of a statistically 1D planar (for the reasons discussed earlier) single-reaction wave in a constant-density turbulent flow, which is not affected by the wave. On the one hand, the invoked simplifications are very strong and make the studied problem substantially different from the problem of propagation of a premixed flame in a turbulent flow, because Lewis number and preferential diffusion, thermal expansion, complex chemistry and other effects can play an important role in the latter case, as reviewed, e.g., in Refs. [36][37][38][39] or shown in recent papers [40,41]. Accordingly, from purely numerical perspective, the present simulations are inferior to modern DNS studies that allow for both thermal expansion and complex chemistry in intense (u /S L 1) turbulence, e.g., [41][42][43][44][45][46][47]. However, the invoked assumptions are fully adequate to the qualitative goal of the present work. Indeed, if the transient effects are of significant importance under simple studied conditions, then, there are no solid reasons to disregard such effects under other, more realistic conditions.
On the other hand, the invoked simplifications allow us (i) to perform DNS by substantially varying both u /S L and L/δ F (in a range of L/δ F > 1), which is necessary for studying the scaling exponents, but is still not feasible in the case of a complex chemistry, and (ii) to obtain reliable statistics for both fully-developed and developing waves by significantly increasing the number of samples using a method valid solely in a constant-density case, as discussed in Section 2.2.
It is also worth noting that, first, turbulent flame speeds evaluated using single-step and multi-step chemistry were recently compared in two independent DNS studies [46,48]. Obtained results show that "mean turbulent flame properties such as burning velocity and fuel consumption can be predicted with the knowledge of only a few global laminar flame properties" [46, p.294] and "the global mechanism is adequate for predicting flame speed" [48, p.53]. Moreover, target-directed experiments [49] performed using the well-recognized Leeds fan-stirred bomb facility do not show a notable effect of combustion chemistry on turbulent flame speed either. Such effects are commonly expected to be of substantial importance when local combustion extinction occurs, but this phenomenon is beyond the scope of the present study.
Second, the assumption of a constant density, invoked to greatly improve statistical sampling, appears to be of minor importance for the goals of the present study for the following three reasons. (I) The vast majority of approximations of experimental data on U T or S T  do not invoke the density ratio σ , thus, implying a weak influence of σ on U T or S T . (II) Recent target-directed experiments [50], as well as earlier measurements [49], did not reveal a substantial influence of σ on U T or S T . (III) Recent DNS studies, e.g., [51, or [52, Fig. 2a], do not indicate such an influence either.
Thus, DNS of propagation of a statistically 1D, planar, single-reaction wave in forced, constant-density, homogeneous, and isotropic turbulence was performed by numerically solving the following 3D Navier-Stokes and reaction-diffusion equations, as well as ∇ · u = 0, in a fully periodic rectangular box of a size of x × × using a uniform cubic grid of N x ×N ×N cells (N x /N = x / = 4) and an in-house solver [53] developed for low-Mach-number reacting flows. Here, ∂ t designates partial derivative with respect to time, p is the pressure, a vector-function f is used to maintain turbulence intensity by applying energy forcing at low wavenumbers [54], c is the reaction progress variable (c = 0 and 1 in reactants and products, respectively), is the reaction rate, τ R is a constant reaction time scale, while parameters Ze = 6 and τ = 6 are counterparts of the Zel'dovich number Ze = a (T b − T u )/T 2 b and heat-release factor τ = σ −1, respectively. Substitution of c = (T − T u ) / (T b − T u ) into the exponent in Eq. 3 results in the classical Arrhenius law with the activation temperature a . Therefore, Eq. 3 offers an opportunity to mimic the behavior of the reaction rate in a flame by considering constant-density reacting flows.
As the state of the reacting mixture is characterized with a single scalar c, the simulated problem is associated with Le = 1 and Sc = P r.

DNS attributes
Because the DNS attributes are discussed in detail elsewhere [55][56][57][58][59], we will restrict ourselves to a very brief summary of the simulations.
The boundary conditions were periodic not only in transverse directions y and z, but also in direction x normal to the mean wave surface [56]. This was possible, because (i) the thickness of the entire wave brush was significantly smaller than the length x of the computational domain in each simulated case at each instant and (ii) the c-field did not affect the velocity and pressure fields in the studied constant-density cases. Accordingly, Eq. 2 was not solved in narrow layers upstream and downstream of the wave brush, but the reaction progress variable was equal to zero and unity, respectively, within these layers, e.g., Here, 1 and x le = x le (t) is the axial coordinate of the reactant boundary of the layer where (2) was solved. Due to the wave propagation, the time derivative dx le /dt was negative at all instances with the exception of instances when x le (t) = 0. At these instances, (i) the identical reaction wave entered the computational domain through its right boundary (x = x ) and (ii) c (x, t) dropped from unity to zero on the plane x = x (1 − ), but (2) was not solved in the vicinity of this plane. Such a method allowed us to strongly improve sampling statistics by simulating a number of cycles of wave propagation through the computational domain, but the method is justified only in the case of ρ =const and ν =const.
The initial turbulence field was generated by synthesizing prescribed Fourier waves [60] with an initial rms velocity u 0 and the forcing scale l 0 = /4. The initial turbulent Reynolds number Re 0 = u 0 l 0 /ν was changed by changing the domain width , with the numerical resolution x = y = z = x /N x = /N being the same in all cases. Subsequently, for each Re 0 , a forced incompressible turbulent field was simulated by integrating (1) with the same forcing vector-function f [61]. In all cases, the velocity, length, and time scales showed statistically stationary behavior at t > t * = 3.5τ 0 T = 3.5l 0 /u 0 [55,56], the turbulence achieved statistical homogeneity and isotropy over the entire domain [55,56], u ∼ = u 0 , and a ratio of L 11 / was about 0.12 [58]. Here, L 11 is the longitudinal length scale of the turbulence, evaluated at t > t * .
To obtain fully-developed statistics, i.e., mean characteristics of the statistically stationary stage of turbulent wave propagation, a planar wave c(x, 0) = c L (ξ ) was released at Here, ξ = x − x 0 and c L (ξ ) is a pre-computed laminar wave profile. Subsequently, evolution of a long-living field c(x, t) was simulated by solving (2). Sampling of the fully developed statistics was started at t = t * = 3.5τ 0 T and was performed over a time interval longer than 50τ 0 T . For that purpose, the time-dependent mean valueq(x, t) of a quantity q(x, t) was evaluated by averaging DNS data over transverse coordinates, followed by computing the fully-developed profile ofq(x) by averagingq(x, t) over time. Finally, x-coordinates were mapped tō c(x)-coordinates.
To study transient effects, the same pre- The transient simulations were run over 2τ 0 T before being reset. Subsequently, at t = t * + 2jτ 0 T , where j ≥ 1, the flow was again populated by M new profiles of c L (ξ ) and the transient simulations were repeated. Time-dependent mean quantities q t (x, t) were evaluated by averaging DNS data over transverse coordinates and over the entire ensemble [m = 1, . . . , M and various time intervals j (x, t). Then, x-coordinates were mapped to c t (xt)-coordinates, as discussed in detail elsewhere [56].
Such a method (i.e., simulations of M independent transient fields) significantly increased the sampling counts for calculating transient statistics. It is worth remembering, however, that such a method can only be used for simulating processes that do not affect the flow, e.g., a constant-density reaction wave addressed here.
Both fully-developed and transient bulk consumption velocities were calculated by integrating the mean reaction rate along the normal to the mean reaction wave brush, i.e., Various cases were set up by selecting a turbulent field and specifying the speed S L and thickness δ F = D/S L of the laminar reaction wave, with the required reaction time scale τ R in Eq. 3 being found in 1D pre-computations of the laminar wave. Because the reaction waves did not affect the flow, the choice of a turbulent field was independent of the choice of S L and δ F . Accordingly, three turbulence fields A (N x = 256, Re 0 = 50, η/ x = 0.68), B (N x = 512, Re 0 = 100, η/ x = 0.87), and C (N x = 1024, Re 0 = 200, η/ x = 1.07), characterized by different , l 0 , τ 0 T ,and L 11 , were generated and propagation of five different sets of reaction waves in each turbulence field was simulated. Within each set of waves, the speed S L was varied, but the thickness δ F retained the same value due to an appropriate adjustment of the Schmidt number Sc. Characteristics of all 15 cases are shown in Table 1, where η = ν 3 / ε 1/4 is the Kolmogorov length scale, P e = u L 11 / (S L δ F ) is turbulent Péclet number, ε is the dissipation rate averaged over the computational domain and time at t > t * . Moreover, extra cases were designed to show weak sensitivity of computed results to grid resolution, L 11 / , etc., but these results are discussed elsewhere [57,58].

Results and Discussion
To reveal significant influence of transient effects on the scaling exponents referred to, the DNS data on both fully-developedŪ T and developing U T (t) turbulent consumption velocities were processed using the same method. More specifically, data on Y = U T /u , U T /S L , (U T − S L ) /u , (U T − S L ) /S L vs. X = u /S L , Da, Ka, P e, u /S L (L 11 /δ F ) d were approximated with Y = aX b using least square fit. 1 In the case of X = u /S L (L 11 /δ F ) d , the power exponent d was varied from -4 to 4 with a step of 0.01 and the exponent that yielded the lowest value of = 1 − R 2 was selected. Here, is the coefficient of determination [62] and K = 15 is the number of studied cases. The aforementioned sets of expressions for X and Y were used to process the DNS results, because various fits to experimental data bases on U T (or S T ) can be found in the literature 32], e.g., (U T − S L ) /u = f (Da) [2], U T /u = f (Da) [3,14], U T /u = f (Ka) [3,4,9,11] The present DNS data on the fully-developed consumption velocities are best fitted, see triangles in Fig. 1, withŪ T /u ∝ u /S L −0.5 (L 11 /δ F ) 0.6 , i.e., the scatter = 1 − R 2 is the lowest in this case. In Ref. [58], a larger number of cases was analyzed andŪ T /S L ∝ P e 0.5 best fitted those data. For the present data, this fit is shown in circles in Fig. 1. Another fit to the present DNS data is plotted in squares in Fig. 1, because this fit is relevant to developing reaction waves, as will be discussed later. It is worth stressing that the DNS data on the fully-developed consumption velocity were already analyzed and compared with experimental data and theoretical results in Ref. [58]. Here, these DNS data are solely reported for comparison with the DNS data on the developing U T . Accordingly, the reader interested in a more detailed discussion of the behavior ofŪ T is referred to Ref. [58]. The DNS data on the developing consumption velocities were analyzed using three different methods. The simplest method consists in fitting the DNS data obtained at the same normalized wave-development time θ = t d /τ T . This method of the data processing mimics measurement of a turbulent consumption velocity, performed at certain distance z from flame-stabilization zone in a statistically stationary flame. Such turbulent consumption velocities were experimentally obtained, e.g., by Verbeek et al. [17] and by Sattler et al. [63] from V-shaped (rod-stabilized) flames. It is worth remembering that, as discussed in detail elsewhere [3,22], development of a statistically stationary premixed turbulent flame occurs during advection of the flame element by the mean flow, similarly to development (decay) of statistically stationary turbulence behind a grid or to development (growth) of a statistically stationary turbulent mixing layer. In such cases, the flame (turbulence, or mixing layer, respectively) development time can roughly be estimated as follows t d ≈ z/U , where U is the mean flow velocity in the z-direction. Accordingly, the normalized wave development time θ is equal to z/ (Uτ T ) in this case.
Certain representative results are shown in Fig. 2. At a low value of θ , see Fig. 2a, the data are very well fitted with U T /u ∝ u /S L −0.99 , see squares. Fitting with U T /u ∝ u /S L −0.99 (L 11 /δ F ) −0.1 does a little bit better job, see triangles, due to the use of the third fitting parameter. While circles in Fig. 2a might also look to be weakly scattered, this impression is wrong, as shown in insert. It indicates that a dependence of U T /S L on P e is weakly pronounced at θ = 0.2. Thus, at low θ , the consumption velocity is almost proportional to S L and almost independent of u . Such a result could be claimed to be expected, because U T (θ = 0) = S L , but, in any case, this result shows that transient effects can significantly change the scaling exponents. When the normalized time θ is increased, scatter of the DNS data around U T /u ∝ u /S L α and U T /S L ∝ P e β is increased and decreased, respectively, see squares and  Fig. 2b, c, and d, as well as circles and triangles in Fig. 3a. Due to the use of an extra parameter, U T /u ∝ u /S L γ (L 11 /δ F ) δ fits to the DNS data well at all θ , see squares in Fig. 3a, but the power exponents vary with increasing θ . In particular, triangles-down and squares in Fig. 3b show that the scaling exponents γ and δ, respectively, increase with increasing θ (because γ < 0, its magnitude decreases) and tend to values yielded by the best fit to the DNS data on the fully-developedŪ T , see dashed and double-dashed-dotted lines, respectively. Scaling exponents α in U T /u ∝ u /S L α and β in U T /S L ∝ P e β show the same trend, see circles and triangles-up, respectively, in Fig. 3b, but the scatter is increased and decreased, respectively, with increasing θ , see circles and triangles, respectively, in Fig. 3a.
Since the majority of measurements of turbulent burning velocity in V-shaped flames, e.g., [11,16,17], and, especially, in Bunsen (rim-stabilized) flames, e.g., [5,10,12,15,19,21,31], were performed by averaging experimental data over the entire flame volume or over a significant part of it, another method of processing the DNS data was also used in order to evaluate similarly averaged values of the developing U T . For this purpose, let us consider a simple mean-reaction-wave surface that has a conical shape, with the cone base radius and height being equal to R and H , respectively. Such a simple configuration is akin to the shape of the mean surface of a Bunsen flame stabilized on a nozzle rim of a radius R. Then, the volume-averaged consumption velocity U T is equal to where U is the mean flow velocity at the nozzle, i.e., at z = 0. Moreover, the bulk consumption velocity evaluated using the reactant flux through the nozzle, i.e., πR 2 U , should be equal to the bulk consumption velocity evaluated using the local value of U T (z) and integrated over distance z from the nozzle. Consequently, The volume-averaged U T (θ m ) evaluated by integrating U T (θ ) obtained in the DNS mimics turbulent burning velocity measured in experiments with statistically stationary Bunsen flames, e.g., [5,10,12,15,19,21,31]. If mixture composition, the rms turbulent velocity u , and the mean nozzle velocity U are varied in such experiments, but a ratio of u /U retains the same value, then, the maximal normalized flame-development time θ m = (H/U) / L 11 /u also retains the same value for the same H . Accordingly, evaluation of the scaling exponents for U T at various constant θ m is relevant to such experiments.
Results of processing the set of U T (θ m ), which were evaluated using Eq. 9 and U T (θ ) extracted from the DNS data in all 15 cases A1-A5, B1-B15, and C1-C5, are shown in Figs. 4 and 5. These results look similar to results obtained for the instantaneous developing U T (θ ), which are reported in Figs. 2 and 3. The major difference consists of U T (θ m ) < U T (θ ), which is obvious, because the function U T (θ ) integrated in Eq. 9 is an increasing function. Thus, while the two methods of processing the same DNS data on U T (θ ) clearly show substantial influence of the wave development on the considered scaling exponents, the simulated effect does not mimic an influence of measurement method on the scaling exponents. Nevertheless, application of one more method to processing the same DNS data revealed the latter influence also. The third method applied to process the same DNS data on U T (θ ) is associated with an experiment with an expanding statistically spherical premixed turbulent flame. In such experiments, e.g., [4,9,13,14,18,20], measurements are performed in a bounded range of the flame kernel radii R f , with this range being independent on mixture composition, turbulence characteristics, pressure, and temperature. To mimic such measurements, we numerically integrated the following equation Here, r f = R f /L 11 is the mean wave radius normalized using L 11 , and u t (θ ) = U T (θ )/u is taken from the DNS data. Subsequently, the DNS dependencies of U T (θ ) were transformed to U T r f (θ ) using results of integration of Eq. 10 and the fits discussed earlier were applied to sets of 15 values of U T r f (θ ) associated with the same wave "radius"r f in different cases A1-A5, B1-B5, and C1-C5.
Obtained results are reported in Figs. 6 and 7. Data on U T fitted (i) at the same normalized time, see Figs. 2, 3, 4 and 5, and (ii) at the same normalized wave radius r f , see Figs. 6 and 7, show similar qualitative trends, but, in the latter case, the difference between the developing U T and the fully-developedŪ T is less pronounced. Accordingly, the scaling exponents that yield the best fits to the DNS data depend not only on the wave-development time, but also on a method (whether U T is studied at the same θ or at the same r f ) used to compare data computed in different cases. These effects can cause significant variations of the scaling exponents q v , q s , and q r for turbulent consumption velocity as a function of u , S L , and L 11 /δ F , respectively, see Table 2 and note that values of 1 ≤ θ ≤ 2 or 1 ≤ r f ≤ 2 are typical for various experiments with statistically stationary or expanding premixed turbulent flames, respectively. In particular, Table 2 shows that the transient effects result in increasing q v and q r , but decreasing q s with wave-development time, see also the fully-developed scaling exponents reported in the right column titled with ∞.
More specifically, inspection of Table 2 shows that the scaling exponent q s is significantly larger for the volume-averaged consumption velocity U T (θ ), associated with Bunsen and V-shaped flames, when compared to the consumption velocity U T r f (θ ) associated with expanding statistically spherical flames. 2 The opposite trend is well pronounced for the scaling exponent q v . These results imply that development of premixed turbulent flames can affect the scaling exponents for burning velocity not only directly, but also indirectly. In other words, the scaling exponents obtained from premixed turbulent flames of different configurations may be different not only due to the well-known effects of the mean-flame-brush curvature and the mean flow non-uniformities, but also due to the flame development, even if the different flames are at the same stage of their development. 2 Direct comparison of the scaling exponents reported in Table 2 with fits to various measured data may be misleading not only due to (i) simplicity of the simulated problem, (ii) the significant scatter of the fits, and (iii) the lack of detailed information on flame-development time in many experimental papers, but also and mainly due to (iv) substantial sensitivity [35] of the scaling exponents to three other effects (the influence of the mean curvature of a flame brush on U T , the influence of the divergence of the mean flow of fresh reactants on S T , and the influence of the growth of mean flame brush thickness on S T ), which are not addressed in the present study, but play an important role in various laboratory premixed turbulent flames.   Such effects can, at least in part, explain significant scatter of the scaling exponents reported in different papers. The results presented above have the following implications. First, the transient effects should appropriately be taken into account when comparing experimental data reported by different research groups or when testing a model of premixed turbulent combustion. For instance, a widely accepted practice of assessment of such models consists in directly comparing the model expression for the fully-developed U T (or the scaling exponents in the expression) with a fit to experimental data on U T or S T . However, this method ignores the transient effects and, therefore, appears to be flawed. When using experimental data to test a numerical model, a method of evaluation of turbulent flame speed or burning velocity in simulations should be as close as possible to the method adapted to measure the data.
Second, while experimental data obtained from expanding premixed turbulent flames convey some information on the transient effects, e.g., see Refs. [3,9,13,14,18,20,49,50], development of U T or S T has yet been addressed in a few experimental investigations of statistically stationary premixed turbulent flames [17,63]. Definitely, the issue requires further target-directed research. In particular, dependencies of the local values of U T or/and S T on distance from flame-stabilization zone should be measured in statistically stationary premixed turbulent flames.
Third, a phenomenological method of processing experimental data on S T R f , obtained from expanding statistically spherical premixed turbulent flames, was proposed to be used in Ref. [64], where the method was also validated against seven experimental databases reported by six independent research groups. Subsequently, the method was theoretically substantiated [30] by assuming self-similarity of the mean structure of developing premixed turbulent flames [65,66], with that assumption being well supported by various experiments discussed in detail elsewhere [3,22,67,68]. On the contrary, a similar method has not yet been elaborated for statistically stationary premixed turbulent flames and this gap should be filled. To do so, a few model expressions for the transient U T (θ ), available in the literature, e.g., see Refs. [3,69], could be used and, eventually, further developed. However, as already noted in the previous paragraph, results of measurements of the local U T or S T at different distances from flame-holder are strongly required for this purpose.
Finally, it is worth noting that the aforementioned model expressions for the transient U T (θ ) also predict variations in the scaling exponents during flame development. For instance, the following equation [70] U T (θ ) U T (θ → ∞) = 1+ 1 θ e −θ −1 1/2 (11) which was developed within the framework of the Flame Speed Closure (FSC) model [22,25] and was validated in RANS simulations [24,25,35,71] of various experiments with statistically stationary premixed turbulent flames, predicts that q v (θ → ∞) = q v (θ 1) + 0.5 and q L (θ → ∞) = q L (θ 1) + 0.5 also. Here, q L is the scaling exponent with respect to an integral length scale L. In other words, both scaling exponents increase during flame development, with the total increase in q v or q L being equal to 0.5. The same increasing trend is shown for q v and q r in Table 2, with the magnitude of the total increase being sufficiently close to 0.5 at least for the latter scaling exponent (and for q v if U T (θ ) and U T (θ m ) are considered).

Conclusions
While turbulent consumption velocity U T is well and unambiguously defined in the simplest case of a statistically planar 1D turbulent reaction wave, the wave development can significantly affect scaling exponents for U T as a function of the rms turbulent velocity, laminar wave speed, and a ratio of the turbulence and laminar wave length scales. Moreover, the scaling exponents depend on a method used to compare values of U T , i.e., the scaling exponents found by processing data obtained at the same wave-development time may be substantially different from the scaling exponents found by processing data obtained at the same wave size.