Transition from pulled to pushed fronts in premixed turbulent combustion: Theoretical and numerical study

Combustion and Flame 162 (2015) 2893–2903

Contents lists available at ScienceDirect

Combustion and Flame j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m b u s t fl a m e

Transition from pulled to pushed fronts in premixed turbulent combustion: Theoretical and numerical study Vladimir A. Sabelnikov a,⇑, Andrei N. Lipatnikov b a b

ONERA – The French Aerospace Lab., F-91761 Palaiseau, France Department of Applied Mechanics, Chalmers University of Technology, Gothenburg 412-96, Sweden

a r t i c l e

i n f o

Article history: Received 8 January 2015 Received in revised form 16 March 2015 Accepted 26 March 2015 Available online 29 April 2015 Keywords: Premixed turbulent combustion Turbulent flame speed Countergradient transport Theory Pulled and pushed traveling waves

a b s t r a c t This paper extends a previous theoretical study (Sabelnikov and Lipatnikov, 2013) of the influence of countergradient transport (CGT) on the speed of a statistically stationary, planar, 1D premixed flame that passively propagates in homogenous turbulence in the form of a traveling wave, i.e. retains its mean thickness and structure. While two particular models of the mean rate of product creation were addressed in the previous article, with the shape of the rate as a function of the Favre-averaged combustion progress variable being concave in both cases, the present paper deals with a more general model that subsumes both concave functions and functions with an inflection point, i.e. a point where the function changes from being concave to convex or vice versa. In this more general case, transition from pulled (flame speed is controlled by processes localized to the flame leading edge) to pushed (flame speed is controlled by processes within the entire flame brush) flames can occur both due to interplay of the nonlinear reaction term and a nonlinear convection term associated with CGT and due to the change of the shape of the reaction term in the absence of CGT. Explicit pushed traveling wave solutions to the studied problem are theoretically derived and conditions under that developing flames approach either pushed or pulled traveling wave solution are obtained by analyzing the governing equations at the flame leading edge and invoking the steepness selection criterion which highlights traveling wave with the steepest profile at the leading edge. Other analytical results include conditions for transition from pulled to pushed premixed turbulent flames, dependence of flame speed on the magnitude of the CGT term and the shape of the mean reaction rate, analytical expressions for the mean thickness of the pushed flames and turbulent scalar flux within the pushed flames. All these theoretical findings are validated by results of unsteady numerical simulations of the initial boundary value problem with steep initial wave profiles. Ó 2015 The Combustion Institute. Published by Elsevier Inc. All rights reserved.

1. Introduction Due to extreme complexity of various highly nonlinear and multiscale phenomena associated with premixed turbulent combustion, its theoretical analysis is commonly simplified to a study of propagation of a statistically stationary, planar, 1D premixed flame in homogeneous turbulence which is not affected by the flame. Moreover, in such theoretical studies, the state of the mixture is often characterized with a single scalar variable c called combustion progress variable (c = 0 and c = 1 in unburned mixture and combustion products, respectively) and the flame propagation is modeled with a single balance equation [1]

q

@ ~c @ ~c @ qu00 c00 ~ þ u þq ¼ Wð~cÞ; @t @x @x

which involves two unclosed terms; turbulent scalar flux F t ¼ qu00 c00 and the mean mass rate W of product creation. Here, t is the time, x and u are a spatial coordinate and the flow velocity, respectively, q ~ ¼ qq=q ~ designate the Favre-av and q00 ¼ q  q is the gas density, q eraged and fluctuating quantities, respectively, with the Reynolds . averages being denoted with over-bars, e.g. q If the turbulent scalar flux in Eq. (1) is closed by invoking the gradient transport (GT) approximation

 Dt F t ¼ qu00 c00 ¼ q ⇑ Corresponding author. E-mail addresses: [email protected] (V.A. Sabelnikov), lipatn@ chalmers.se (A.N. Lipatnikov).

ð1Þ

@ ~c ; @x

ð2Þ

where Dt > 0 is a turbulent diffusion coefficient, then, Eq. (1) belongs to a wide class of partial differential equations, which are called convection–diffusion–reaction (CDR) equations. Such

http://dx.doi.org/10.1016/j.combustflame.2015.03.016 0010-2180/Ó 2015 The Combustion Institute. Published by Elsevier Inc. All rights reserved.

2894

V.A. Sabelnikov, A.N. Lipatnikov / Combustion and Flame 162 (2015) 2893–2903

equations are widely studied in various fields of science starting from the pioneering work by Fisher [2] and Kolmogorov et al. [3] who developed the well-known KPP theory of such equations. More precisely, Kolmogorov et al. [3] considered a diffusion–reaction (DR) equation, which results from Eqs. (1) and (2) in the case ~ ¼ 0, a constant Dt, and a constant density. They also assumed a of u concave, i.e. W½ð~c1 þ ~c2 Þ=2 > ½Wð~c1 Þ þ Wð~c2 Þ=2 for any 0 6 ~c1 < ~c2 6 1, ‘‘hump’’-like shape of the source term in the DR equation, i.e. Wð~cÞ in Eq. (1), with the highest slope dW=d~c being reached at ~c ¼ 0. Moreover, Wð0Þ ¼ Wð1Þ ¼ 0 and Wð~cÞ > 0 if 0 < ~c < 1. Logistic expression ~cð1  ~cÞ and term ~c  ~cq with q > 1 are wellknown simple examples of such a concave source term multiplied with a proper time scale. Kolmogorov et al. [3] sought for permanent monotonous traveling wave (TW) solutions ~cðx; tÞ ¼ ~cðXÞ, where X = x + u0t, to the studied DR equation with the boundary conditions of ~cð1Þ ¼ 0 and ~cðþ1Þ ¼ 1. Such a TW solution satisfies the following ordinary differential equation (ODE)

u0 = u0,min is basically different, i.e. it has the steepest profile of ~cðXÞ at ~c ! 0, with the speed u0,min being controlled by processes in the entire wave. Such a TW solution is called pushed wave [4,5]. Examples of pushed waves that are modeled by CDR equa~ ¼ 0, and source terms that are convex sometions with q = const, u where within the interval of 0 < ~c < 1 can be found e.g. in Refs. [16–18]. Aronson and Weinberger [15] have also solved the speed selection problem and have shown that a solution to the IBVP associated with Eqs. (1) and (2) approaches the TW solution characterized by the slowest propagation speed u0,min provided that the initial wave profile is steep enough, e.g. the step function. If applied to premixed turbulent combustion, the general mathematical result by Aronson and Weinberger [15] indicates that the KPP approach cannot be used to determine Ut if an invoked closure relation for Wð~cÞ has an inflection point in the interval of 0 < ~c < 1. A laminar premixed flame modeled with the following CDR equation

2 d~c d ~c u0 ¼ Dt 2 þ Wð~cÞ; dX dX

q

ð3Þ

retains its shape ~cðXÞ, and moves at a constant speed u0. Kolmogorov et al. [3] have analytically solved this nonlinear boundary value problem (BVP) and have proved that it does not have a unique solution, i.e. there is a family of TW solutions and a continuous spectrum of eigenvalues, i.e. the TW speeds u0, qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi bounded with a minimum (slowest) speed u0;KPP ¼ 2Dt W 0 ð0Þ where W 0 ð~cÞ ¼ dW=d~c. The constraint of u0 P u0;KPP has been obtained by linearizing the ODE at the leading edge ~c ! 0 (an unstable state) of a TW. Moreover, Kolmogorov et al. [3] analytically studied the initial boundary value problem (IBVP) with localized initial conditions, i.e. ~cðx; t ¼ 0Þ ¼ ~c0 ðxÞ, with ~c0 ðxÞ being a monotonously increasing function in an interval of x1 < x < x2, but ~c0 ðxÞ ¼ 0 if x < x1 and ~c0 ðxÞ ¼ 1 if x > x2 > x1. A step function H(x  x1) is a particular case of such an initial profile ~c0 ðxÞ with x1 = x2. Finally, Kolmogorov et al. [3] have solved the speed selection problem and have proved that solutions to the IBVP approach the TW solution characterized by the slowest propagation speed provided that the initial wave profiles are steep enough, e.g. the step function. TW solutions found by Kolmogorov et al. [3] are called ‘‘pulled’’ waves [4,5] in order to stress that the speed of such a solution is determined from the linear analysis of the problem at the unstable state ~c  1, i.e. the TW is pulled by its leading edge. In the combustion literature, the KPP approach was widely applied to the CDR equation that results from Eqs. (1) and (2) in order to determine fully developed turbulent burning velocity Ut [6–14], which was associated with the slowest TW speed in the cited papers. Various expressions for Ut were obtained in Refs. [6–14] depending on invoked closure relation for the mean rate Wð~cÞ, but Ut was controlled by Dt and the slope of Wð~cÞ at ~c ! 0 in all these studies. It is worth noting, however, that the source term in the DR equation is not necessary to be a concave function in a general case. The maximum slope of Wð~cÞ can be reached inside the interval of 0 < ~c < 1 and the source term can have an inflection point, where the dependence of W on ~c changes from being concave to for any convex, i.e. W½ð~c1 þ ~c2 Þ=2 < ½Wð~c1 Þ þ Wð~c2 Þ=2 0 6 ~c1 < ~c2 6 1, or vice versa. Such a case was studied by Aronson and Weinberger [15] who have proved that the BVP for qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Eq. (3) has a solution for any u0 P u0;min P u0;KPP ¼ 2Dt W 0 ð0Þ, but, contrary to the KPP case, the slowest speed u0,min cannot be found by linearizing Eq. (3) at ~c  1. They have proved also that TW solutions characterized by a continuous spectrum of u0 > u0,min are pulled, whereas, if u0,min > u0,KPP, then, the TW solution with

  @c @c @ @c qUðcÞ þ qu ¼ qD þ @t @x @x @x s0   Zeð1  cÞ ;  exp  1  ð1  r1 Þð1  cÞ

ð4Þ

where U(c) = 0 at 0 6 c 6 c1 < 1 and Ze  1, is another well-known example of a pushed TW, as shown by Zeldovich and FrankKamenetsky [19–21]. Here, c = (T  Tu)/(Tb  Tu) is the normalized temperature, Ze ¼ ðT b  T u ÞE=ðR0 T 2b Þ is the Zeldovich number, Tu and Tb designate temperatures of unburned and burned gases, respectively, E is the activation temperature of a single global reaction that combustion chemistry is reduced to, R0 is the universal gas constant, D is molecular diffusivity, s0 is a reaction time scale, U(c) is a weakly non-linear (when compared to the Arrhenius exponential term) function such that U(c1 < c < 1) > 0 and U(1) = 0, and r = Tb/Tu. The KPP theory is not applicable to this problem, because dU/dc = 0 at the leading edge, and there is a single TW solution [19– 21]. In the case of Uð0 6 c < 1Þ > 0, Eq. (4) can have an intermediate asymptotic solution [21–23], which retains its shape and speed over a long time interval and is associated with a pushed or pulled wave if r  1 is much larger than or on the order of R0Tb/E, respectively, with the speed of the pulled flame being consistent with the KPP theory. As far as premixed turbulent combustion is concerned, transition from pulled to pushed TWs is straightforwardly relevant to the following long-standing basic issue. On the one hand, production of flamelet surface area within a premixed turbulent flame brush is often considered to control its speed including the speed of its leading edge [24,25]. On the other hand, within the framework of the leading point concept pioneered by Zeldovich [26], see also pp. 450–452 in Ref. [21], the mean turbulent flame speed is hypothesized to be controlled by processes localized to the leading edge of the mean flame brush, whereas the production of flamelet surface area within the flame brush adjusts itself to the speed of the leading edge. The readers interested in further discussion of the leading point concept and facts that indirectly support it are referred to books [27,28], a review paper [29], and recent contributions by Lieuwen et al. [30–34]. To the best of the present authors knowledge, TW solutions to Eqs. (1) and (2) have not yet been investigated analytically in the case of Wð~cÞ with an inflection point at 0 < ~c < 1 in the turbulent combustion literature. The major goal of the present work is to fill this gap. Moreover, the following issue specific to turbulent flames will also be addressed. As predicted by Clavin and Williams [35] and Bray and Libby [36], experimentally discovered by Moss [37] and by Yanagi and Mimura [38], and documented in a number of subsequent experimental and direct numerical simulation (DNS) studies reviewed

V.A. Sabelnikov, A.N. Lipatnikov / Combustion and Flame 162 (2015) 2893–2903

elsewhere [39], Eq. (2) yields wrong direction of turbulent scalar flux in many weakly turbulent flames, where the gradient of ~c and the turbulent flux of ~c have the same sign. This phenomenon known as countergradient diffusion or countergradient transport (CGT) impedes straightforwardly applying the KPP approach to weakly turbulent premixed combustion. In a recent paper [40], the present authors invoked a simple model to close the countergradient contribution to the flux Ft in Eq. (1) and investigated that CDR equation following earlier mathematical studies by Kolmogorov et al. [3], Murray [41,42], and von Saarloos et al. [43–45]. Speeds of TW solutions to the considered CDR equation were analytically found in two cases associated with two widely used models of the mean rate of product creation, with the shape of Wð~cÞ being concave in both cases. In the present work, the previous analysis is extended to a more general case of the source term Wð~cÞ with an inflection point in order to investigate transition from pulled to pushed TW solutions not only due to interplay of the nonlinear source term and nonlinear convection term associated with CGT, but also due to the change of the shape of the nonlinear source term. Thus, the major difference between the present study and our previous investigation [40] consists of invoking a more general expression for Wð~cÞ. The paper is organized as follows. Section 2 describes a mathematical model addressed by us. In Section 3, the boundary value problem (BVP) for finding TW solutions is stated and the KPP linear analysis is applied to find a continuous spectrum of the TW speeds. Explicit pushed TW solutions to the problem are analytically derived in Section 4. In Section 5, the speed selection problem is addressed and physically realizable TW solutions to the considered IBVP are selected. In Section 6, the obtained analytical results and, in particular, a method invoked to select the physically realizable TW solutions are validated in numerical simulations. Conclusions are drawn in Section 7. 2. Mathematical formulation and governing equations

~  @q u @q þ ¼ 0: @x @t

ð5Þ

If the combustion progress variable c is associated with the normalized temperature (T  Tu)/(Tb  Tu) and qT = quTu, then, q = qu/ (1 + sc) and the mean density is determined by the following equations

q ð~cÞ ¼ qu Rð~cÞ; R ¼

r¼

1 < 1; 1 þ s~c

s¼

qu  1 ¼ r  1; qb

T b qu ¼ ; T u qb

ð6Þ

which can be obtained by applying the definition of the Favre aver~ ¼ qq=q  to q = 1/q. Moreover, in order to link the Reynolds c age q and Favre ~c averages of the combustion progress variable, we will use the well-known Bray-Moss-Libby (BML) model [1,8], which is based on a hypothesis that probability of finding intermediate (between unburned and burned) states of a reacting mixture, i.e. 0 < c < 1, is much less than unity everywhere within a premixed turbulent brush. Under this assumption,

q ~c ¼ qb c; c ¼ rRð~cÞ~c;

ð7Þ

and Eqs. (6) and (7) yield

RðcÞ ¼ 1  ac 6 1;

0 6 a ¼ 1  r1 < 1:

Based on earlier studies [46–51], the turbulent scalar flux is closed as follows

 Dt F t  qu00 c00 ¼ q

ð8Þ

@ ~c  U B ~cð1  ~cÞ: þq @x

ð9Þ

This closure relation is based on (i) separately modeling gradient and countergradient contribution to the flux [46,48], i.e.

F t  qu00 c00 ¼ F GT þ F CGT ;

ð10Þ

(ii) gradient diffusion closure of FGT given by Eq. (2), and (iii) the following simple model of the countergradient contribution

 ~cð1  ~cÞ F CGT ¼ sSL q

ð11Þ

introduced by Veynante et al. [46]. Here, SL is the laminar flame speed. For the sake of generality, the velocity U B P 0 in Eq. (9) is not associated with sSL in the present work, but is considered to be an input parameter of the problem, required to characterize the magnitude of the CGT flux. In particular, such an approach allows us to subsume not only the model by Veynante et al. [46], who, following Bray [52], highlighted the difference between velocities conditioned on unburned and burned gas due to pressure drop across inherently laminar flamelets, but also a model by Zimont and Biagioli [53,54], who highlighted the difference between the conditioned velocities due to different accelerations of the unburned and burned gas by the combustion-induced mean pressure gradient within the mean flame brush. Recent models [49–51,55,56] address both mechanisms. For instance, Eq. (9) is supported by Eqs (29) and (31) from Ref. [50]. It is worth also noting that the following analysis holds even if Dt in Eq. (9) is associated with the sum of turbulent and molecular diffusivities, but the turbulent contribution dominates under typical conditions, i.e. it is larger by a factor of the order of turbulent Reynolds number. The following algebraic expression

W¼

Let us consider a statistically planar 1D premixed flame that propagates from right to left in statistically stationary homogeneous turbulence. The flame is modeled by Eq. (1) supplemented with the Favre-averaged continuity equation

2895

xð~c; rÞ qu ~cð1  ~cÞ  qu ~cð1  ~cÞ ~ ¼ f ðc Þ ¼ f 1 ðc; rÞ; sf sf sf

ð12Þ

is invoked to close the mean rate of product creation in the present study. Various expressions for a time scale sf as a function of the rms turbulent velocity u0 , integral length scale L, laminar flame speed SL and thickness dL, etc. can be found in the literature, e.g. see Table 1 in Ref. [57]. For the goals of the present work, a particular expression for sf is not required, provided that this time scale is independent of t, x, or ~c (this assumption is required in order to derive exact analytical solutions in Section 4, whereas conclusion regarding transition from pulled to pushed TW solutions at a sufficiently large U B P 0, which will be drawn later, holds even if the time scale depends on ~c, as discussed in the end of Section 3). Without loss of generality, we will assume that the time scale sf is set so that, if the derivative dW=d~c is finite at ~c ! 0, then, dW=d~c ! s1 f . The function f ðcÞ is given by the following relation

f ðcÞ ¼ 1 þ 2e1 c þ e2 c2 :

ð13Þ

It involves two input parameters e1 and e2, which should be set so that constraint of f ðcÞ > 0 is satisfied in an interval of 0 < c < 1. If e1 = e2 = 0, then, f ðcÞ ¼ 1 and the source term given by Eq. (12) has a concave shape, with the maximum slope being reached at ~c ¼ 0, in line with the KPP constraints. The function f ðcÞ is introduced in order to study not only concave functions Wð~cÞ, but also functions with an inflection point. Dependencies of the normalized mean rate Wð~cÞ= maxfWð~cÞg on Favre ~c and Reynolds c averages of the combustion progress variable, calculated using Eqs. (12) and (13), various e1, and e2 = e1a (in this particular case, the problem admits an analytical

2896

V.A. Sabelnikov, A.N. Lipatnikov / Combustion and Flame 162 (2015) 2893–2903

2 -0.875 -0.5 -0.25 0 0.25 0.5 1

0.8 1

first derivative

normalized source term

1

0.6 -0.875 -0.75 -0.5 -0.25 0 1 4

0.4

0.2

0

0

0.2

0

-1

0.4

0.6

0.8

1

-2

Reynolds-averaged c

0

0.2

(a)

normalized source term

0.6

0.8

1

Fig. 2. Dependencies of the derivative dðW sf Þ=d~c yielded by Eqs. (12) and (13) on the Favre-averaged combustion progress variable, computed for various e1 specified in legends and e2 = e1a. Filled circles show inflection points. Density ratio is equal to 7.5.

1

0.8

~c ¼ 0. However, this KPP constraint is not satisfied in three other cases, i.e. e1 = 0.25, 0.5, and 1.0 and the KPP results could not be valid even if UB = 0. To simplify the subsequent analysis, let us rewrite Eqs. (1), (9), (12) and (13) in terms of the Reynolds average c

0.6 -0.875 -0.75 -0.5 -0.25 0 1 4

0.4

0.2

0

0.4

Favre-averaged c

0

0.2

0.4

0.6

0.8

1

Reynolds-averaged c

(b) Fig. 1. Dependencies of the normalized reaction rate Wð~cÞ= maxfWð~cÞg given by Eqs. (12) and (13) on the (a) Favre- and (b) Reynolds-averaged combustion progress variable, computed for various e1 specified in legends and e2 = e1a. Density ratio is equal to 7.5.

solution, see Eq. (47) in Section 4) are shown in Fig. 1. When e1 is increased, the values of ~cm and cm associated with peaks of Wð~cÞ and WðcÞ, respectively, are also increased. Accordingly, the use of Eq. (13) allows us to study functions Wð~cÞ and WðcÞ characterized by various ~cm and cm , respectively. It is worth noting that recent DNS results indicate that ~cm (or cm ) associated with the maximum flame surface density is increased by u0 /SL, see Fig. 3 in Ref. [58], Fig. 2 in Ref. [59], Fig. 7 in Ref. [60], or Fig. 7 in Ref. [61], thus, implying an increase in e1 by u0 /SL. As far as eventual dependence of e1 on the density ratio is concerned, DNS data by Nishiki et al. [62,63] indicate a weak increase in ~cm when r is decreased, see Fig. 3 in Ref. [64], thus, implying reduction in e1 by r, but the issue definitely needs a target-directed study. It is worth remembering that the two closure relations invoked in our previous paper [40], i.e. W sf ¼ ~cð1  ~cÞ and W sf ¼ R2 ~cð1  ~cÞ ¼ cð1  cÞ=r, are characterized by constant ~cm ¼ 0:5 and cm ¼ 0:5, respectively. Note that these two simple closure relations, which are widely used in the literature, e.g. see Ref. [9], results from Eqs. (12) and (13) if e1 = e2 = 0 and 2 2 e1 ¼ q1 u ðqu  qb Þ ¼ a, e2 ¼ e1 ¼ a , respectively. Dependencies of dðW sf Þ=d~c on the Favre-averaged combustion progress variable, calculated, see curves, for various e1 specified in legends and e2 = e1a, are plotted in Fig. 2, with filled circles showing inflection points for e1 = r/(1 + r), 0.25, 0.5, and 1.0. In the case of e1 = r/(1 + r), i.e. the lowest value of e1 consistent with f ð0 < ~c < 1Þ > 0, there is an inflection point, but the KPP analysis still holds if UB = 0, because the derivative dðW sf Þ=d~c peaks at

q q @ c q q @c q u 2 b þ q u~ u 2 b q @t q @x   @ q q q q @ c u b þ q U B  2 cð1  cÞ  q Dt u 2 b @x q q @x qu qb qu cð1  cÞ  f ðcÞ; ¼ 2 q sf

ð14Þ

where we used the following well-known relations

qq qu ð1  cÞ ¼ q ð1  ~cÞ; qb c ¼ q ~c; ~cð1  ~cÞ ¼ u 2 b cð1  cÞ; q q ¼ qu þ ðqb  qu Þc @ ~c qu qb @ c ¼ 2 ; @s q @s

s ¼ t;

s ¼ x:

ð15Þ

Eq. (14) reads

   @ c @ c @ 1 @ c q þ q u~ þ q 2 U B cð1  cÞ  Dt  @x q @t @x @x   q cð1  cÞ  f ðcÞ: ¼ u

sf

ð16Þ

Normalization of Eqs. (5) and (16) results in

@R @Rv þ ¼0 @h @n

ð17Þ

and

R

   @ c @ c @ 1 @ c ¼ cð1  cÞf ðcÞ; NB cð1  cÞ  þ Rv þ R2 @n R @h @n @n

ð18Þ

respectively, where

h¼

t

sf

;

x n ¼ pffiffiffiffiffiffiffiffiffiffi ; D t sf

~ u D t = sf

v ¼ pffiffiffiffiffiffiffiffiffiffiffiffi ;

UB NB ¼ pffiffiffiffiffiffiffiffiffiffiffiffi ; Dt =sf

ð19Þ

are the dimensionless time, distance, the Favre-averaged axial velocity, and a number basically similar to the Bray number [46,52], respectively. In the present work, NB is considered to be an input parameter. The normalized turbulent scalar flux pffiffiffiffiffiffiffiffiffiffiffiffi f t  qu00 c00 =ðqu Dt =sf Þ is equal to

2897

V.A. Sabelnikov, A.N. Lipatnikov / Combustion and Flame 162 (2015) 2893–2903

  @ ~c 1 @ c f t ¼ N B R~cð1  ~cÞ  R ¼ NB cð1  cÞ  @n rR @n

ð20Þ

by virtue of Eqs. (9) and (15). The use of constant (independent of time and distance) time sf, pffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffi length Dt sf , and velocity Dt =sf scales in the transformation of Eqs. (5), (9) and (16) to Eqs. (17), (18), and (20), respectively, implies that the turbulence is stationary and homogeneous, i.e. combustion does not affect it. Such a simplification is typical for many theoretical studies of premixed turbulent flames [6– 14,40,65] and appears to be necessary for obtaining analytical results. Nevertheless it is worth stressing that even if sf, Dt, and UB depend on c either straightforwardly or due to dependence of local turbulence characteristics on c, then, conclusion regarding transition from pulled to pushed TW solutions at a sufficiently large U B P 0, which is drawn in the end of Section 3, see discussion of Eq. (39), still hold provided that the governing equations are normalized using sf ðc ! 0Þ and Dt ðc ! 0Þ. Boundary conditions are as follows

cð1; hÞ ¼ 0;

cð1; hÞ ¼ 1:

ð21Þ

The initial boundary value problem for the nonlinear CDR Eq. (16) with boundary conditions given by Eq. (21) and initial condition of cðn; h ¼ 0Þ ¼ c0 ðnÞ describes a combustion front that propagates into an unstable state, where c ¼ 0 and dW=dc > 0, from a stable state, where c ¼ 1 and dW=dc < 0. It is a common practice to consider the IBVP with localized initial conditions, i.e. c0 ðnÞ ¼ 0 if n < n1 and c0 ðnÞ ¼ 1 if n > n2 > n1, with c0 ðnÞ being a monotonously increasing function in an interval of n1 < n < n2. A step function H(n  n1) is a particular case of such an initial profile c0 ðnÞ with n1 = n2. Typically, a solution to the IBVP tends to a TW solution cðfÞ with time and we will study first the latter solution. Here, f = n + Kh is the wave variable. In the coordinate framework attached to the unburned mixture, K denotes the normalized TW speed. Eq. (18) reduces to the KPP DR equation with a concave nonlinear source term provided that v = 0, R = 1, NB = 0, e1 = e2 = 0. In a general case of R ¼ RðcÞ 6 1 and NB – 0, Eq. (18) involves two more nonlinear terms in addition to the nonlinear source term. First, the mean density depends on c, thus, making the diffusion term in Eq. (18) nonlinear. Second, a nonlinear convection term induced by the CGT appears in Eq. (18), see the first term in square brackets. In the subsequent sections, we will show that the two extra nonlinearities and the convex shape of the nonlinear source term restrict the domain of applicability of the KPP method (linear analysis at the leading edge) for determining the TW speed.

In the coordinate framework attached to the unburned gas,

v(1) = 0 and the continuity equation (22) reads RðK þ v Þ ¼ K:

ð25Þ

Substitution of this result into Eq. (23) yields

K

   dc d 1 dc þ R2 NB cð1  cÞ  df R df df ¼ cð1  cÞð1 þ 2e1 c þ e2 c2 Þ:

It is worth stressing that the value of K should be positive. Indeed, let us divide Eq. (26) by R2

K

   1 dc d 1 dc   þ N c ð1  c Þ  B df R2 df df R 1 ¼ 2 cð1  cÞð1 þ 2e1 c þ e2 c2 Þ: R

Substitution of Rðn; sÞ ¼ R½cðfÞ, cðn; sÞ ¼ cðfÞ, and v(n, s) = v(f), where f = n + Kh, into Eqs. (17) and (18) results in the following ordinary differential equations (ODE)

ðK þ v Þ

dR dv d þR ¼ ½RðK þ v Þ ¼ 0; df df df

RðK þ v Þ

ð22Þ

   dc d 1 dc NB cð1  cÞ  þ R2 df R df df

¼ cð1  cÞð1 þ 2e1 c þ e2 c2 Þ;

ð23Þ

respectively. The boundary conditions are given by

cð1Þ ¼ 0;

cð1Þ ¼ 1:

ð24Þ

ð27Þ

This equation can be transformed to

K

   1 1 dR d 1 dc NB cð1  cÞ  þ df a df df R 1 ¼ 2 cð1  cÞð1 þ 2e1 c þ e2 c2 Þ R

ð28Þ

using the chain rule 1 1 dR dR dc ¼ ; df dc df

1

dR a a ¼ ¼ dc ð1  acÞ2 R2

ð29Þ

and Eq. (8). Integration of Eq. (28) from f = 1 to f = 1 yields

K

 Z 1 1 1 cð1  cÞð1 þ 2e1 c þ e2 c2 Þdf: 1 ¼ 2 a 1a 1 R 1



ð30Þ

Thus,

K ¼ ð1  aÞ

Z

1

1

1

R2

cð1  cÞf ðcÞdf > 0;

ð31Þ

because the integral in Eq. (31) is positive for f ðcÞ > 0 and 0 6 a < 1. Eq. (26) reads

R

2 d c 2

 ½K þ aNB cð1  cÞ þ RNB ð1  2cÞ

df þ cð1  cÞð1 þ 2e1 c þ e2 c2 Þ ¼ 0:

dc dc dc þa df df df ð32Þ

The second order OD Eq. (32) can be reduced to the following nonlinear first order ODE

Rpp0 ¼ ap2 þ ½K þ aNB cð1  cÞ þ RNB ð1  2cÞp  cð1  cÞð1 þ 2e1 c þ e2 c2 Þ

3. Traveling wave equation. Boundary value problem and linear analysis

ð26Þ

ð33Þ

for the derivative p  dc=df considered to be a function of c, i.e. 2

p0  dp=dc, as d c=df2 ¼ dp=df ¼ ðdp=dcÞðdc=dfÞ ¼ p0 p in this case. Eq. (33) describes curves in the phase plane ðc; pÞ. Boundary conditions for Eq. (33) read

pðc ¼ 0Þ ¼ pðc ¼ 1Þ ¼ 0:

ð34Þ

Either the BVP given by Eqs. (23) and (24) or the BVP given by Eqs. (33) and (34) is an eigenvalue problem. We have to find positive eigenvalues K such that the obtained trajectory pðcÞ > 0, 0 < c < 1 connects two singular points (1,0) and (0,0) in the phase space. The complete solution to the BVP requires considering the global behavior of trajectories pðcÞ in the phase space, e.g. see Ref. [3]. However, the necessary conditions for the existence of such trajectories can be found by linearizing the problem with respect to c  1 at the leading edge of a TW, i.e. in the vicinity of the unstable equilibrium point (0, 0), similarly to the KPP analysis [3].

2898

V.A. Sabelnikov, A.N. Lipatnikov / Combustion and Flame 162 (2015) 2893–2903

Substitution of

cðf ! 1Þ ! 0;

4. Pushed TW. Exact solution

R ! 1;

p ! jc;

j>0

ð35Þ

into Eq. (33) yields the following dispersion relation

j2  ðK þ NB Þj þ 1 ¼ 0;

ð36Þ

which links the TW speed K and the decay rate (wave steepness)

j > 0 of the profile of cðfÞ near its leading edge. Note that this dispersion relation holds not only for FCGT and f ð~cÞ given by Eqs. (11) and (13), respectively, but also for any flux FCGT and any source term xð~c; rÞ, see Eq. (12), that scale as ~c when ~c ! 0. The quadratic dispersion relation given by Eq. (36) has two branches

K þ NB jþ ¼ þ 2

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðK þ N B Þ2  1; 4

K þ NB j ¼  2

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðK þ N B Þ2  1; 4

ð37Þ

ð38Þ

with j+j = 1. Because cðfÞ is a monotonous positive function, both j+ and j should be real. Therefore TW solutions can exist if and only if the spectrum of the TW speeds is bounded from below by a minimum speed Kmin, i.e. a necessary condition for the existence of a TW solution is as follows

K P Kmin P KKPP ¼ 2  N B :

p0 ¼ kð1  2cÞ;

p ¼ kcð1  cÞ;

c ¼ 1=½1 þ expðkfÞ:

ð40Þ

if K = KKPP. In the KPP case (v = 0, R = 1, NB = 0, e1 = e2 = 0), Kmin = KKPP = 2 and the decay rates of TWs with K > KKPP belong to the j-branch given by Eq. (38), as proved by Kolmogorov et al. [3]. If NB > 2, then, the value of KKPP given by Eq. (39) is negative, whereas eigenvalues of the BVP should be positive by virtue of Eq. (31). Therefore, the linear analysis does not allow us to find a TW speed if N B P N cr B and pulled TWs do not exist. The value of N cr B < 2 cannot be determined from the linear analysis and will be found in Section 5, see Eq. (57). To solve the eigenvalue BVP in this case, the behavior of trajectories pðcÞ in the entire phase space should be analyzed. Accordingly, if N B P N cr B , then, the smallest speed of the wave Kmin > KKPP is controlled by its interior profile, i.e. the wave is pushed. Therefore, transition from pulled to pushed TW occurs at N B ¼ N cr B . It is worth noting that, the transition cannot be revealed using the linear analysis at the leading edge in the simplified case of NB = 0, v = 0, R = 1, e1 > 0, and e2 = 0. Nevertheless, in this particular case, Hadeler and Rothe [16] have proved that pulled TWs exist only if 0 6 e1 < 1, while TWs are pushed TW if 1 6 e1 . All conclusions drawn in the previous paragraph hold even in a more general case, i.e. for any flux FCGT such that F CGT ð~c ! 0Þ ! U B ~c and for any source term xð~c; rÞ such that xð~c  1; rÞ 6 ~c, because linearization of this more general BVP at ~c  1 results in Eqs. (36)– (38), with unity being substituted with @ x=@ ~cð~c ! 0Þ 6 1 in each of this equation. Furthermore, even if sf, Dt, and UB vary within the flame brush, then, we can arrive at Eqs. (36)–(39) and again draw the same conclusions by normalizing the governing equations using sf ðc ! 0Þ and Dt ðc ! 0Þ, with NB designating pffiffiffiffiffiffiffiffiffiffiffiffi U B = Dt =sf evaluated at c ! 0 in this case. However, in order to obtain exact analytical solutions discussed in the next section, sf, Dt, and UB should be independent of ~c.

ð41Þ

where k > 0 in order to satisfy the boundary conditions. Substitution of Eq. (41) into Eq. (33) yields, after some algebra, the following equality 1

1

ðK þ N B  k  k Þ  2ðN B  k þ k

e1 Þc þ ½ðNB  kÞa  k1 e2 c2 ¼ 0: ð42Þ

ð39Þ

The last equality, i.e. KKPP = 2  NB, results from the constraint of vanishing discriminant D = (K + NB)2  4 = 0 of the quadratic Eq. (36). Note that

jþ ¼ j ¼ jKPP ¼ 1

The mathematical theory and various methods of finding pushed TW solutions to DR equations are discussed in detail elsewhere [4,5,15]. In particular, the wave steepness (decay rate) of the pushed TW with slowest speed Kmin has been shown to belong to the j+-branch given by Eq. (37). It should be stressed that this result is valid only for the TW with Kmin, i.e. the spectrum K(j) for pushed TWs consists only of an isolated discrete point K = Kmin(j+). The decay rates of TWs with K > Kmin > KKPP are associated with the j-branch given by Eq. (38), similarly to the pulled TW in the KPP case. In the present section, we will use a single method of finding pushed TW solutions, which consists of searching for explicit analytical solutions to the BVP given by Eq. (32) or (33). To the best of the present authors’ knowledge, the method was pioneered by Murray [41,42] and was subsequently applied e.g. by von Saarloos [43] or Benguria and Depassier [66]. Let us seek a solution pðcÞ in a form of a quadratic function

In order for this equality to hold for all c from an interval of 0 6 c 6 1, the following three constraints 1

K þ NB  k  k ðNB  kÞ þ k

1

¼ 0;

ð43Þ

e1 ¼ 0;

ð44Þ

1

ðNB  kÞa  k

e2 ¼ 0

ð45Þ

should be satisfied. Note that Eq. (43) is identical to the dispersion relation given by Eq. (36). This is not surprising, because the dispersion Eq. (36) should be satisfied at cðfÞ  1. Eqs. (43) and (44) result in

K ¼ k1 ð1 þ e1 Þ;

ð46Þ

while Eqs. (43) and (45) yield the following relation

e2 ¼ e1 a

ð47Þ

that links two parameters e1 and e2. Therefore, the explicit analytical solution given by Eq. (41) can exist only for the following function

f ðcÞ ¼ ð1 þ 2e1 c  e1 ac2 Þ;

ð48Þ

where

e1 P emin ¼ 

1 1þs qu ¼ ¼ > 1; 2a 2þs qu þ qb

ð49Þ

in order for f ðcÞ to be positive in the interval of 0 < c < 1. In the rest of the present paper, we will assume that Eqs. (47) and (49) hold. Substitution of Eq. (46) into Eq. (43) yields 2

k  NB k  e1 ¼ 0:

ð50Þ

If e1 P 0, then, the positive solution to the quadratic Eq. (50) is as follows

NB k ¼ kþ ¼ þ 2

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi NB þ N2B þ 4e1 N2B þ e1 ¼ ; 4 2

ð51Þ

and exists if NB P 0 (NB > 0 provided that e1 = 0). For negative

e1 P r=ð1 þ rÞ, a solution to Eq. (50) exists only if

2899

V.A. Sabelnikov, A.N. Lipatnikov / Combustion and Flame 162 (2015) 2893–2903

pffiffiffiffiffiffiffiffiffi NB P 2 e1 :

ð52Þ

In this case, there are two real positive solutions, but k+ given by Eq. (51) is larger and k+k = e1 < 1. Therefore, k < 1. Bearing in mind the theoretical result cited in the beginning of this section (the pushed TW is characterized by a larger j+ > 1), we will restrict ourselves to the k+ -solution given by Eq. (51). 5. Selection problem As shown above, for given e1 and NB, the considered problem admits a continuous spectrum of pulled TW solutions and a discrete spectrum of pushed TW solutions, similar to the case studied by Aronson and Weinberger [15]. As discussed in Introduction, the selection problem consists of finding the sole TW solution that is approached by a solution to the IBVP, given by Eqs. (13), (17), (18) and (21), provided that the initial profile of cðn; h ¼ 0Þ is sufficiently steep. This TW solution is commonly considered to be physically realizable. To solve the selection problem, we will apply the steepness criterion [43–45], which states that a physically realizable TW solution is the TW solution with the maximum decay rate j. Accordingly, an inequality of k+ > jKPP = 1 should be satisfied in order for the explicit pushed TW solution found in the previous section to be physically realizable. If k+ < jKPP = 1, then, the pulled solution with the KPP decay rate jKPP = 1, which results from the linear analysis, is physically realizable. It is worth remembering, however, that the steepness criterion has not yet been proved rigorously in the mathematical literature [4,5] and the present work aimed, in particular, at assessing this criterion. For this purpose, a numerical study of the considered IBVP was performed, as discussed in the next section. Comparison of the decay rate k of the exact solution, given by Eq. (51) with jKPP = 1 yields the following results. First, if e1 < 1 (and e1 P r=ð1 þ rÞ in order for the source term to be positive at 0 < c < 1), then, the speeds of the physically realizable TW solutions and their decay rates are equal to

(

K¼

KKPP ¼ 2  NB 2ð1þe1 Þ pffiffiffiffiffiffiffiffiffiffiffi ffi 2 NB þ

N B þ4e1

if if

0 6 NB < 1  e1 NB > 1  e1

ð53Þ

and

(

j¼

jKPP ¼ 1 NB þ

pffiffiffiffiffiffiffiffiffiffiffi ffi 2 N B þ4e1 2

if > 1 if

0 6 NB < 1  e1 NB > 1  e1

ð54Þ

ð55Þ

and decay rate

jþ ¼

NB þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi N2B þ 4e1 2

NB ¼ Ncr B ¼ maxf0; 1  e1 g:

ð57Þ

The pushed TWs are physically realizable if N B > N cr B. In the particular cases of (i) e1 = 0, e2 = 0, f ðcÞ ¼ 1 and (ii) e1 = a, e2 = a2, f ðcÞ ¼ 1  2ac þ a2 c2 ¼ ð1  acÞ2 ¼ R2 ; Eqs. (53) and (54) reduce to results obtained in Ref. [40], see cases n = 0 and n = 2, respectively. In the particular case of e1 > 0, NB = 0, and q = const, Eqs. (53) and (54) reduce to results derived by Hadeler and Rothe [16]. Eqs. (53)–(56) show that both the flame speed and thickness, which is inversely proportional to the decay rate for the pushed TWs, see Eq. (41) which yields maxfpðcÞg ¼ kþ =4, are decreased when the magnitude NB of the countergradient contribution to the scalar flux is increased, in line with common expectations. Moreover, the speed of the physically realizable pushed TW is increased by e1, whereas its thickness is decreased when e1 is increased. Finally, the speed does not depend on the density ratio if both e1 and NB are kept constant, but the speed of the pushed TW that is physically realizable at NB = 1  emin depends on the density ratio due to the dependence of emin on s, given by Eq. (49). It is also worth remembering that NB depends on the density ratio and both the time scale sf in Eq. (12) and parameter e1 could also depend on it. Accordingly, the obtained results should not be considered to show that turbulent flame speed is independent of the density ratio. The speed can depend on r due to dependence of NB, sf, or e1 on it. When applied to a TW, Eq. (20), which models the entire scalar flux, reads

ft ¼

  1 dc : NB cð1  cÞ  df rR

ð58Þ

For pushed TWs, the flux can be determined analytically. Indeed, substitution of Eqs. (41) and (51) into Eq. (58) yields

respectively. Note that if e1 < 0 and NB > 1  e1, then, Eq. (52) is satpffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffi isfied, because 1  e1  2 e1 ¼ 1 þ je1 j  2 je1 j ¼ ð1  je1 j2 Þ > 0. If e1 < 0 and NB < 1  e1, then, the decay rate k+ of the pushed TW, given by Eq. (51) provided that Eq. (52) holds, is lower than jKPP = 1 and the pulled TW with the KPP decay rate (jKPP = 1) is selected. Accordingly, the constraint given by Eq. (52) is of minor importance if e1 < 0 and NB < 1  e1. It is also worth noting that, if NB = 0 and 1/2r < e1 < 1, then, the KPP results for the TW speed K and decay rate j hold in spite of existence of an inflection point and a peak of Wð~cÞ at ~c > 0, e.g. see curves e1 = 0.25, 0.5, and 1.0 in Fig. 2. Second, if e1 P 1, then, 1  e1 6 0 and, therefore, the pushed TWs with the following speed

K ¼ j1 þ ð1 þ e1 Þ

are physically realizable independently of N B P 0. In particular, if NB = 0, then, transition from pulled to pushed TW solutions is still possible and occurs at e1 = 1. Therefore, this transition can be controlled not only by the interplay of the nonlinear convection and the nonlinear reaction term, as found in Ref. [40] for two particular concave source terms Wð~cÞ, but also solely by the nonlinear reaction term provided that it has an inflection point. In summary, transition from pulled to pushed TW solution occurs at

ð56Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi NB  N2B þ 4e1 1 cð1  cÞ: ft ¼ ðN  k Þcð1  cÞ ¼ 2rR rR B þ

ð59Þ

Thus, the flux is countergradient, i.e. ft > 0, (gradient, ft < 0) for negative (positive) e1, with the direction of the flux being the same in the entire flame brush independently of c. Eq. (59) shows that a single Bray number is not sufficient to characterize the direction of turbulent scalar flux. In the studied case, the flux direction is solely controlled by the dependence of the mean rate of product creation on the mean combustion progress variable provided that the flame is pushed, i.e. NB > 1  e1. Certainly, an increase in NB is required in order for the pushed flame to be physically realizable and the flux to be countergradient in the case of e1 < 0. To the contrary, the flux is always gradient independently of N B P 0 if e1 > 1. For the pulled TWs, the flux cannot be found analytically due to lack of an explicit solution. Nevertheless, Eqs. (35), (40) and (58) allow us to determine the direction of turbulent scalar flux at the LE of the mean flame brush. Indeed, substitution of Eqs. (35) and (40) into Eq. (58) yields

ft !

1

r

ðNB c  jKPP cÞ ¼

1

r

ðNB  1Þc;

at c ! 0;

ð60Þ

V.A. Sabelnikov, A.N. Lipatnikov / Combustion and Flame 162 (2015) 2893–2903

with the pulled TWs being physically realizable provided that NB < 1  e1. The flux is countergradient if NB > 1, which is possible only for negative e1. Otherwise, the flux is gradient at the leading edge of the pulled TW, but can change direction within the flame brush. Inability of a single Bray number to characterize the direction of turbulent scalar flux was already noted in a couple of papers reviewed elsewhere [39]. For instance, Robin et al. [50] pointed out that the magnitude of countergradient turbulent scalar flux should depend on the wrinkling factor or available amount of the flame surface, thus, implying the straightforward dependence of the flux on Wð~cÞ.

2

normalized flame speed

2900

-0.5 -0.5 -0.2 -0.2 0 0 0.2 0.2

1.5

1

0.5

0

0

1

2

3

4

5

NB

6. Numerical simulations

(a) 2.5

normalized flame speed

Because the analytical results obtained in the previous section are strongly based on the steepness criterion [43–45], which is widely accepted, but has not yet been proved rigorously, the above analysis requires validation in numerical simulations. In particular, it should be shown that solutions to the initial boundary value problem given by Eqs. (13), (17), (18) and (21) tend to the TW solutions obtained in the previous sections provided that the initial profile of ~cðh ¼ 0; nÞ is sufficiently steep. For this purpose, the initial boundary value problem was numerically solved. In the simulations, the initial profile of ~cðh ¼ 0; nÞ varied linearly from zero at n = n1 to unity at n = n2 > n1, with independence of computed TW speed and thickness on Dn = n2  n1 < < 1 being checked. Typically, n1 was about 0.1 of the width of the computational ~ ðh ¼ 0; nÞ was domain and Dn = 0.4. The initial velocity v

0.5 0.5 1.0 1.0 2.0 2.0 4.0 4.0

2

1.5

1

0.5

0

0

1

2

3

4

5

NB

normalized flame speed

2

(b)

-0.5 -0.5 -0.2 -0.2 0 0 0.2 0.2

1.5

1

Fig. 4. Normalized speeds of fully-developed turbulent flames, computed numerically, see symbols, and calculated using Eqs. (53) and (55), see lines, for various e1 specified in legends. Constant density.

0.5

0

0

1

2

3

5

4

NB

(a)

normalized flame speed

2.5

0.5 0.5 1.0 1.0 2.0 2.0 4.0 4.0

2

1.5

1

0.5

0

0

1

2

3

4

5

NB

(b) Fig. 3. Normalized speeds of fully-developed turbulent flames, computed numerically, see symbols, and calculated using Eqs. (53) and (55), see lines, for various e1 specified in legends. Density ratio is equal to 7.5.

independent of the spatial coordinate n. In order to keep the flame in the middle of the computational domain, the inlet velocity v(h, n = 0) of the unburned mixture far ahead the flame was set equal to K determined by Eq. (53) or (55) for e1 < 1 or e1 > 1, respectively. Simulations were performed utilizing an in-house code exploited in previous studies [40,47]. For various NB and e1, Eqs. (17) and (18) were discretized by spatially integrating them using a uniform staggered grid, with the power law scheme being applied jointly to the convection and diffusion term [67]. The nonlinear convection term associated with CGT was discretized using the first-order upwind scheme. The CGT and reaction terms were linearized following general recommendations by Patankar [67]. All the aforementioned schemes were implicit. The unsteady term was discretized using the first-order scheme, with the time step being evaluated from a constraint on the maximum allowed change Dc 6 d  1 of c during a single time step, with weak sensitivity of computed results to variations in d being checked. The magnitude of d was varied from 5  105 to 5  106. Because the thickness of a numerical TW solution depended strongly on NB and e1, e.g. see Fig. 5a in Ref. [40], the width of the computational domain was varied from 25 to 500 depending on NB and e1, but regions of c  1 and 1  c  1 were well inside the computational domain in each case. Number of grid points was varied from 16,000 to 64,000, with weak sensitivity of computed results to halving this number being checked. Convergence of each unsteady numerical solution to a TW was checked by monitoring evolution of mean burning velocity, flame speed, and flame brush thickness, see Eqs. (61)–(63), respectively. Time required to reach the TW

2901

V.A. Sabelnikov, A.N. Lipatnikov / Combustion and Flame 162 (2015) 2893–2903

solution depended strongly on NB and e1. For instance, it was of unity order for NB = 5 and e1 = 4, but about 200 for NB = 0 and e1 = 0. In the case of NB = 0 and e1 = e2 = 0, numerical TW solutions obtained for various density ratios were validated by comparing computed TW speeds with the theoretical expressions resulting from the KPP analysis [3]. In the case of e1 = e2 = 0 and a constant density, numerical TW solutions obtained for various NB were validated by comparing computed TW speeds with the analytical results of Murray [41,42] and the numerical results of Méndez and Fort [68]. In the present paper, we restrict ourselves to reporting results computed for fully-developed flames associated with the TW solutions to Eqs. (17) and (18). The normalized turbulent burning velocity Ut, flame speed St, and thickness dt were determined as follows

Ut ¼

Z

normalized flame brush thickness

4

-0.5 -0.5 -0.2 -0.2 0 0 0.2 0.2

3

2

1 1

2

3

4

5

NB

(a)

nmax

~cð1  ~cÞf ðcÞdn;

ð61Þ

St ¼ v ðh; n ¼ 0Þ 

dnf ; dh

normalized flame brush thickness

0

ð62Þ

1 dt ¼ : maxf@ c=@ng

ð63Þ

normalized flame brush thickness

where the coordinate nf(t) was evaluated from the constraint of cðh; nf Þ ¼ 0:5. In all simulations, a constraint of Ut = St was satisfied. 4

-0.5 -0.5 -0.2 -0.2 0 0 0.2 0.2

3

3

4

normalized flame brush thickness

(a) 0.5 0.5 1.0 1.0 2.0 2.0 4.0 4.0

1 1

2

3

4

2

3

4

5

Figures 3 and 4 show physically realizable TW speeds and validate the theoretical Eqs. (53) and (55), including independence of the speed on the density ratio, cf. Figs. 3 and 4, provided than NB and e1 are kept constant. Accordingly, Figs. 3 and 4 further support the steepness criterion invoked by us to select physically realizable TW solutions. It is worth remembering that the steepness criterion has not yet been proved rigorously in the mathematical literature [4,5]. The exact solution given by Eqs. (41) and (51) offers an opportunity to evaluate not only the flame speed, but also its thickness

dt ¼

2

0

1

Fig. 6. Normalized mean thicknesses of fully-developed turbulent flames, computed numerically, see symbols, and calculated using Eq. (64), see lines, for various e1 specified in legends. Constant density.

5

3

1

(b)

NB

4

2

NB

1 2

3

0

2

1

0.5 0.5 1.0 1.0 2.0 2.0 4.0 4.0

4

5

NB

4 ¼ kþ

8 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi NB þ N2B þ 4e1

ð64Þ

provided that N B P Ncr B . If e1 < 1, then, the maximum normalized thickness yielded by Eq. (64) is equal to 4 and is reached at NB ¼ Ncr B ¼ 1  e1 . If e1 P 1, then, the maximum normalized thickpffiffiffiffiffi ness is reached at NB = 0 and is equal to 4= e1 . The thickness does not depend on the density ratio provided than NB and e1 are kept constant. All these theoretical results (curves) are validated by numerical data (symbols) shown in Figs. 5 and 6, including the independence of the thickness on the density ratio at constant NB and e1. 7. Conclusions

(b) Fig. 5. Normalized mean thicknesses of fully-developed turbulent flames, computed numerically, see symbols, and calculated using Eq. (64), see lines, for various e1 specified in legends. Density ratio is equal to 7.5.

In the case of a sufficiently general closure relation for the mean rate of product creation Wð~cÞ, which subsumes various dependencies of the rate on the mean combustion progress variable and

2902

V.A. Sabelnikov, A.N. Lipatnikov / Combustion and Flame 162 (2015) 2893–2903

Table 1 Summary of analytical results. Source term Bray number Flame speed

 1þrr 6 e1 < 0 1 < NB < 1 - e1

Flame type Flame thickness

0 < e1 < 1 0 6 N B < 1  e1

NB < 1 K = 2  NB

 1þrr 6 e1 < 0 1  e1 < NB

K¼

Pulled No analytical expression

0 < e1 < 1 2ð1þe1 Þ

pffiffiffiffiffiffiffiffiffiffiffiffi ffi 2

NB þ

CGT at LE

GT at LE

allows both concave functions and functions with an inflection point, the following results were theoretically obtained and numerically validated. An analytical expression, see Eq. (53), for the normalized speed of the pulled flame as a function of the Bray number NB, which characterizes the magnitude of the countergradient contribution to the turbulent scalar flux. The flame speed is decreased when NB is increased, but, if NB is kept constant, the flame speed  ! qu at the leadis independent of the density ratio, because q ing edge of the flame brush. Analytical expressions, see Eqs. (53) and (54), for normalized speed of the pushed flame as a function of the Bray number NB and a parameter e1, which characterizes the shape of Wð~cÞ, i.e. transition from the concave shape to the shape with an inflection point occurs when e1 is increased. The flame speed is decreased when NB is increased, but is increased by e1. Even in the case of the pushed flame, its speed does not depend on the density ratio provided that both e1 and NB are kept constant. It is worth noting, however, that the speed of the pushed TW that is physically realizable at NB = 1  emin depends on the density ratio due to the dependence of emin on s, given by Eq. (49). Moreover, the flame speed can depend on the density ratio e.g. due to dependence of NB on it. Analytical pushed TW solutions, see Eqs. (41), (46) and (51), and expression for the pushed flame thickness, see Eq. (64), and turbulent scalar flux, see Eq. (59), within the pushed flames. The thickness is decreased when either NB or e1 is increased. Depending on e1, the flux can be either gradient (e1 > 0) or countergradient (e1 < 0) in the entire pushed flame brush, including its leading edge. This analytical result proves that a single Bray number is not sufficient to characterize the direction of turbulent scalar flux, which is also controlled by the dependence of the mean rate of product creation on the mean combustion progress variable. Criterion for transition from pulled to pushed flames, see Eq. (57). The transition can occur not only due to interplay of the nonlinear reaction term and a nonlinear convection term associated with CGT (an increase in NB promotes the transition to the pushed flames), but also due to the change of the shape of the reaction term even in the absence of CGT, i.e. the pushed flames are physically realizable if Wð~cÞ has an inflection point and e1 > 1. The obtained analytical results are summarized in Table 1. Because (i) TW solutions that should be approached during flame development were selected by analyzing the governing equations at the flame leading edge and invoking the steepness criterion and (ii) the obtained theoretical results were numerically validated, the present work further supports the steepness selection criterion. Acknowledgments The first author (VS) was supported by ONERA and the second author (AL) was supported by the Swedish Energy Agency and by

N B þ4e1

Pushed p8ffiffiffiffiffiffiffiffiffiffiffiffi ffi dt ¼ 2 NB þ

Flux

1 < e1 0 6 NB

CGT

N B þ4e1

GT

the Chalmers Combustion Engine Research Center (CERC). Helpful discussion with A. Mura and his comments are gratefully acknowledged. The authors express their gratitude to the reviewers of this manuscript for their valuable comments and suggestions. References [1] K.N.C. Bray, J.B. Moss, Acta Astronaut. 4 (1977) 291–319. [2] R.A. Fisher, Ann. Eugenics 7 (1937) 355–369. [3] A.N. Kolmogorov, I. Petrovsky, N. Piskounov, Bjull. MGU A 1(6) (1937) 1–26; English translation in: P. Pelcé (Ed.), Dynamics of Curved Fronts, Academic Press, Boston, USA, 1988, pp. 105–130. [4] U. Ebert, W. van Saarloos, Physica D 146 (2000) 1–99. [5] W. van Saarloos, Phys. Rep. 386 (2003) 29–222. [6] K.N.C. Bray, P.A. Libby, Phys. Fluids 19 (1976) 1687–1701. [7] B. Hakberg, A.D. Gosman, Proc. Combust. Inst. 20 (1984) 225–232. [8] P.A. Libby, Prog. Energy Combust. Sci. 11 (1985) 83–96. [9] K.N.C. Bray, Proc. R. Soc. London A431 (1990) 315–335. [10] C.A. Catlin, R.P. Lindstedt, Combust. Flame 85 (1991) 427–439. [11] J. Duclos, D. Veynante, T. Poinsot, Combust. Flame 95 (1993) 101–118. [12] C.R. Choi, K.Y. Huh, Combust. Flame 114 (1998) 336–348. [13] H. Kolla, J.W. Rogerson, N. Chakraborty, N. Swaminathan, Combust. Sci. Technol. 181 (2009) 518–535. [14] H. Kolla, J.W. Rogerson, N. Swaminathan, Combust. Sci. Technol. 182 (2010) 284–308. [15] D.G. Aronson, H.F. Weinberger, in: J.A. Goldstein (Ed.), Partial Differential Equations and Related Topics, vol. 446, Springer, Berlin, 1975, pp. 5–49. [16] K.P. Hadeler, F. Rothe, J. Math. Biol. 2 (1975) 251–263. [17] A.P. Aldushin, S.I. Khudyaev, Ya.B. Zeldovich, Arch. Combust. 1 (1981) 9–21. [18] P. Clavin, A. Linãn, in: M.G. Velarde (Ed.), Nonequilibrium Cooperative Phenomena in Physics and Related Fields, NATO Science Series B: Advanced Science Institutes Series, vol. 116, Springer, Berlin, 1984, pp. 291–338. [19] Ya.B. Zeldovich, D.A. Frank-Kamenetsky, Doklady Acad. Sci. (USSR) 19 (1938) 693–695 (in Russian). [20] Ya.B. Zeldovich, D.A. Frank-Kamenetsky, Zh. Fiz. Khim. 12 (1938) 100–105 (in Russian). [21] Ya.B. Zeldovich, G.I. Barenblatt, V.B. Librovich, G.M. Makhviladze, The Mathematical Theory of Combustion and Explosions, Plenum Publ. Corp., New York, 1985. [22] Ya.B. Zeldovich, Combust. Flame 39 (1980) 211–214. [23] Ya.B. Zeldovich, Combust. Flame 39 (1980) 219–224. [24] T. Poinsot, D. Veynante, Theoretical and Numerical Combustion, second ed., Edwards, Philadelphia, PA, 2005. [25] D. Veynante, L. Vervisch, Prog. Energy Combust. Sci. 28 (2002) 193–266. [26] Ya.B. Zeldovich, D.A. Frank-Kamenetsky, Turbulent and Heterogeneous Combustion, MMI, Moscow, 1947 (in Russian). [27] V.R. Kuznetsov, V.A. Sabelnikov, Turbulence and Combustion, Hemisphere, New York, 1990. [28] A. Lipatnikov, Fundamentals of Premixed Turbulent Combustion, CRC Press, Boca Raton, FL, 2012. [29] A.N. Lipatnikov, J. Chomiak, Prog. Energy Combust. Sci. 31 (2005) 1–73. [30] P. Venkateswaran, A. Marshall, D.H. Shin, D. Noble, J. Seitzman, T. Lieuwen, Combust. Flame 158 (2011) 1602–1614. [31] P. Venkateswaran, A. Marshall, J. Seitzman, T. Lieuwen, Proc. Combust. Inst. 34 (2013) 1527–1535. [32] A. Amato, T. Lieuwen, Combust. Flame 161 (2014) 1337–1347. [33] A. Amato, M. Day, R.K. Cheng, J. Bell, T. Lieuwen, Proc. Combust. Inst. 35 (2015) 1313–1320. [34] A. Marshall, J. Lundrigan, P. Venkateswaran, J. Seitzman, T. Lieuwen, Proc. Combust. Inst. 35 (2015) 1417–1424. [35] P. Clavin, F.A. Williams, J. Fluid Mech. 90 (1979) 589–604. [36] P.A. Libby, K.N.C. Bray, AIAA J. 19 (1981) 205–213. [37] J.B. Moss, Combust. Sci. Technol. 22 (1980) 119–129. [38] T. Yanagi, Y. Mimura, Proc. Combust. Inst. 18 (1981) 1031–1039. [39] A.N. Lipatnikov, J. Chomiak, Prog. Energy Combust. Sci. 36 (2010) 1–102. [40] V.A. Sabelnikov, A.N. Lipatnikov, Combust. Theory Modell. 17 (2013) 1154– 1175. [41] J.D. Murray, Lectures on Nonlinear-Differential-Equation Models in Biology, Clarendon Press, Oxford, UK, 1977. [42] J.D. Murray, Mathematical Biology 1. An Introduction, third ed., Springer, Berlin, 2002.

V.A. Sabelnikov, A.N. Lipatnikov / Combustion and Flame 162 (2015) 2893–2903 [43] [44] [45] [46] [47] [48] [49] [50] [51] [52] [53] [54] [55] [56] [57]

W. van Saarloos, Phys. Rev. A 37 (1988) 211–229. G.T. Dee, W. van Saarloos, Phys. Rev. Lett. 60 (1988) 2641–2644. W. van Saarloos, Phys. Rev. A 39 (1989) 6367–6390. D. Veynante, A. Trouvè, K.N.C. Bray, T. Mantel, J. Fluid Mech. 332 (1997) 263– 293. A.N. Lipatnikov, J. Chomiak, Combust. Theory Modell. 8 (2004) 211–225. J. Chomiak, J.R. Nisbet, Combust. Flame 102 (1995) 371–386. V. Robin, A. Mura, M. Champion, T. Hasegawa, Combust. Sci. Technol. 182 (2010) 449–464. V. Robin, A. Mura, M. Champion, J. Fluid Mech. 689 (2011) 149–182. V. Robin, A. Mura, M. Champion, Combust. Sci. Technol. 184 (2012) 1718– 1742. K.N.C. Bray, Proc. R. Soc. London A 451 (1995) 231–256. V.L. Zimont, F. Biagioli, Combust. Theory Modell. 6 (2002) 79–101. F. Biagioli, V.L. Zimont, Proc. Combust. Inst. 29 (2012) 2087–2095. A.N. Lipatnikov, Flow Turbul. Combust. 86 (2011) 609–637. A.N. Lipatnikov, Proc. Combust. Inst. 33 (2011) 1497–1504. A.N. Lipatnikov, Int. J. Spray Combust. Dynam. 1 (2009) 39–66.

2903

[58] J. Hult, S. Gashi, N. Chakraborty, M. Klein, K.W. Jenkins, R.S. Cant, C.F. Kaminski, Proc. Combust. Inst. 31 (2007) 1319–1326. [59] N. Chakraborty, J.W. Rogerson, N. Swaminathan, Phys. Fluids 20 (2008) 045106. [60] I. Han, K.Y. Huh, Proc. Combust. Inst. 32 (2009) 1419–1425. [61] E. Lee, K.Y. Huh, Flow Turbul. Combust. 84 (2010) 339–356. [62] S. Nishiki, T. Hasegawa, R. Borghi, R. Himeno, Proc. Combust. Inst. 29 (2002) 2017–2022. [63] S. Nishiki, T. Hasegawa, R. Borghi, R. Himeno, Combust. Theory Modell. 10 (2006) 39–55. [64] A.N. Lipatnikov, V.A. Sabelnikov, S. Nishiki, T. Hasegawa, in: 6th European Combustion Meeting, Lund, Sweden, 25–28 June, 2013, CD, 6pp. [65] V.A. Sabelnikov, P. Bruel, J. Combust. 2011 (2011) 821358. [66] R.D. Benguria, M.C. Depassier, Phys. Rev. E 50 (1994) 3701–3704. [67] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere, New York, 1980. [68] V. Méndez, J. Fort, Phys. Rev. E 64 (2001) 011105.