Transient Response of a Laminar Premixed Flame to a Radially Diverging/Converging Flow

Laminar flamelets are often used to model premixed turbulent combustion. The libraries of rates of conversion
from chemical to thermal enthalpies used for flamelets are typically based on counter-flow, stained laminar planar
flames under steady conditions. The current research seeks further understanding of the effect of stretch on
premixed flames by considering laminar flame dynamics in a cylindrically-symmetric outward radial flow geometry (i.e., inwardly propagating flame). This numerical model was designed to study the flame response when the flow and scalar fields align (i.e., no tangential strain on the flame) while the flame either expands (positive stretch) or contracts (negative stretch, which is a case that has been seldom explored) radially. The transient response of a laminar premixed flame has been investigated by applying a sinusoidal variation of mass flow rate at the inlet boundary with different frequencies to compare key characteristics of a steady unstretched flame to the dynamics of an unsteady stretched flame. An energy index (EI), which is the integration of the source term in the energy equation over all control volumes in the computational domain, was selected for the comparison. The transient response of laminar premixed flames, when subjected to positive and negative stretch, results in amplitude decrease and phase shift increase with increasing frequency. Other characteristics, such as the deviation of the EI at the mean mass flow rate between when the flame is expanding and contracting, are nonmonotonic with frequency. Also, the response of fuel lean flames is more sensitive to the frequency of the periodic stretching compared to a stoichiometric flame. An analysis to seek universality of transient flame responses across lean methane-air flames of different equivalence ratios (i.e., 1.0 to 0.7) using Damköhler Numbers (i.e., the ratio of a flow to chemical time scales) had limited success.


Introduction
Numerical simulation of turbulent combustion has significance in practical devices, such as industrial furnaces, stationary gas turbines, aero-combustors or internal combustion engines. In particular, premixed turbulent combustion is a canonical case for the reactant state of some practical combustion systems, such as HCCI (Homogeneous Charge Compression Ignition) engines. The modelling of such systems is a challenging task due to the unsteady, multi-component and multidimensional nature, and the large range of length and time scales of these flows. These complexities have been discussed in the literature as part of the accurate determination of the rate of conversion of reactants to products, and the sensitivity of the turbulent flame speed to the geometry of the flame [1].
To decouple the small-scale chemistry from large-scale flow features, premixed turbulent flames have often been viewed as ensembles of premixed laminar stretched flames that are wrinkled and potentially torn by the turbulent flow field, in the so-called 'flamelet' approach [2]. The local burning characteristics have been calculated from strained laminar flames as a function of equivalence ratio, pressure, temperature, and strain rate [3], but without consideration of potential impacts of spatial gradients of these quantities over the flame surface or any temporal gradient. The effect of strain rate on the laminar premixed flame has been examined by many studies (e.g., [4]) and it has been shown how different phenomena can affect local burning rates when small sections of the turbulent flame experience a range of strain rates and curvatures locally [5].
Stretch rate (κ) as a mechanism to affect the rate of combustion was first introduced by Karlovitz et al. [6] as the Lagrangian time (t) derivative of an element of the flame surface area (A) as in Eq. (1).
In that work, a planar combustion wave, exposed to a velocity gradient, was analyzed, while the impact of the motion of curved flame fronts normal to itself was first discussed by Markstein [7]. The combination of these two effects; namely, the underlying hydrodynamic strain and the flame surface curvature effects on an interface have been expressed as stretch rate [8,9]. Matalon [8] introduced an expression (Eq. (2)) for calculating stretch rate of a flame in an arbitrary shape on the flame surface indicating by the function ( , ) = 0.
where is the fluid velocity, is the local normal to the flame front and parallel to the gradient in the scalar field (local normal is considered positive when pointing from reactants to products), ! is the speed of an identifiable flame surface feature and thus can be the local flame speed in laboratory coordinates. The first term on the right hand side of Eq. (2) represents the flame curvature effects and results from the divergent/convergent flow field and motion of the flame in that flow field. The second term illustrates the tangential strain rate and results from a non-uniform flow field across a scalar field of the flame. The summation of these two terms creates an equivalence between the strain and curvature effects, but that is strictly true only for an interface, while flames of a finite thickness and multiple scalar species add uncertainty to this equivalence and introduce some unique challenges regarding which interface to consider [10,11].
In counter-flow configurations, due to the planar flame shape, the divergence of the scalar field (∇. = ) becomes zero and subsequently the first term in Eq. (2) is eliminated; however, misalignment of flow and scalar fields ( × ≠ 0) results in a non-zero value for the second term. Therefore, the laminar premixed flame in this geometry is strained. A time-dependent investigation in this geometry was performed by Saitoh and Otsuka [12] for both premixed and diffusion flames numerically and experimentally. They varied the velocity normal to the stagnation plane sinusoidally around its mean.
They concluded that the positional amplitude of temperature and concentration fluctuations decreased with increasing oscillation frequency. In a similar work, Stahl and Warnatz [3] carried out a numerical investigation on the transient response of strained flamelets. They studied the influence of timedependent sinusoidal change of strain rate on the flame front behaviour. The dependence of flame oscillation amplitude and phase shift with frequency of the strain rate was analyzed and it was concluded that the amplitude of flame position oscillation decreases and the phase shift increases with an increase in frequency.
These studies showed that in a counter-flow configuration, the response of a strained flame to a periodic change in flow velocity and strain rate is not instantaneous, implying limitations in validity of the quasisteady laminar flamelet assumption in premixed turbulent modeling. As part of this endeavour, Petrov and Ghoniem [13] revisited the validity of this assumption by analyzing the transient response of The above investigations were concerned with transient response of strained flames, while fewer studies have been concerned with the effects of curvature. Giannakopoulos et al. [14] numerically studied the effect of flame curvature by using spherically outwardly propagating flames. This geometry has also been used extensively to study the influence of flow rate on stabilized [15] and flame-stretch interactions in premixed flames (e.g. [16]). In this geometry, due to the alignment of flow and scalar fields, the tangential strain rate (the second term in Eq. (2)) is null, but the flame is stretched because of the motion of the curved flame. In this condition, flame stretch varies with time as the flame propagates outward and has a positive value at all times.
In premixed turbulent combustion, the local flame experiences both positive and negative stretch.
Kostiuk and Bray [5] studied the distribution of stretch rates on the flamelet surface, it has been shown that in an outwardly propagating spherical flame, 30% to 50% of the flame is under compression (i.e., negative stretch rates) at any time. Therefore, modeling the region of negative stretch rates in terms of consumption velocity and analyzing the mean effects of stretch rate on conversion from reactants to products is of importance. However, as mentioned previously, the geometry considered in the literature to analyze the flamelets, namely, counter-flow flame configurations, only involves positive stretch.
In this paper, through the consideration of an inwardly-propagating premixed cylindrical flame that is forced to radially expand and contract, a numerical investigation has been made to revisit the quasisteady assumptions for laminar flamelet models when the flame is subjected to both positive and negative stretch. Furthermore, since the flame stretch in this work is due to the motion of a curved flame, the second term in Eq. (2) is zero and the effects of tangential strain are eliminated. Finally, in order to study the effect of equivalence ratio on transient response of a laminar flame to a periodic flow, two fuel lean flames are compared with the stoichiometric condition.
In premixed turbulent combustion, the time-varying vortices stretch the flamelets and cause small curvature on flame fronts; therefore, generating libraries to describe the flamelets based on steady laminar stagnation point flames (planar strained flames) has been challenged in some ways [17]. One of these challenges is related to the instantaneous structure of a turbulent flame where the flame front includes positive stretch and negative stretch at certain locations. The planar strained flames, which are only positively stretched, are not capable of describing the concave curved flame fronts towards the fresh gases. As it has been mentioned earlier, most premixed turbulent flames involve negative stretch and highly curved flame fronts, therefore, to generate a more rigorous flamelets library to model a premixed turbulent flame, in addition to a counter-flow configuration, one needs to include the timedependent motion of a curved flame. The model that has been proposed and developed in this paper enables us to study this aspect of flamelets and include the negative stretch in the computations of a laminar premixed flame.

Model Description
It was desired to create one-dimensional radially (cylindrical) outward flow of a pre-mixture of CH 4 and air (O 2 and N 2 ) of specifiable composition and temperature, as depicted in Fig. 1. The mass flow rate of reactants ( ) was imagined to occur along the axial coordinate ( Fig. 1-a) and then be distributed radially outward through a porous cylinder with characteristics to create a uniform diverging flow. It is at the exit plane of the porous cylinder that the computational domain begins. The computational domain ends at the inflow surface of a larger concentric porous cylinder where the identical mass flow is extracted.
The system is assumed to be of sufficient axial length and neglecting of gravity that flow that is onedimension in the radial direction is obtained, and independent of the axial or azimuthal coordinates.
cylindrical flame (red dotted lines in Fig. 1-b) forms between the inlet and outlet boundaries. Although this flame has not been studied experimentally due to difficulties in buoyant stability, it represents unique features for theoretical and numerical studies of laminar premixed flames for two main reasons.
First, similar to a stationary spherical flame [18], the flow field in the domain is radially diverging. In such a condition, the strain of the flow is balanced by the curvature of the flame, and thus the total stretch rate equals zero, which is essential to analyze the separate effects of curvature and strain rate.
Second, oscillating flow at the inlet boundary enables the study of flame creation (expanding) and destruction (contracting).

Governing Equations
The governing equations comprise a set of coupled, nonlinear partial differential equations and algebraic constraints to impose conservation of mass, momentum, energy and the evolution of individual species mass. For the application presented here, the low Mach number approximation is applied. This problem was formulated in terms of primitive variables -pressure, velocity, temperature, and species mass fractions. Since their derivation is presented in several references (e.g., [19]), they are just shown here in vector form: Conservation of total mass: where is the fluid velocity vector; , density; and , time.
Balance of momentum: In Eq. (4), is the static pressure; , the dynamic viscosity of the mixture, and is the unit tensor. As is common in the study of low speed flows of dilute gases, the bulk viscosity is assumed to be negligible [20].

Conservation of energy: [21]
In this equation, is the species index which ranges over the !" species; , temperature; ! , mass fraction of species ; ! and !,! are the specific heats at constant pressure of the mixture and species , respectively; the mixture coefficient of thermal conductivity; ℎ ! the specific enthalpy of species ; ! the molecular weight of species k; ! the molar production rate of species ; and ! diffusion velocity.
Viscous effects and the material derivative of pressure were neglected [22].
Balance of individual species mass: where ! is the volumetric mass production rate of species due to chemical reaction.
The ideal gas equation of state completes the system: where is the mean molecular weight of the mixture (defined as , and is the universal gas constant. The first order nature of the continuity equation makes it stable only neutrally. Numerical stability can be maintained by adding a first order artificial damping term ( ∆ ! ! ) to the continuity equation.
The value of should be selected sufficiently small (10 !! ) to make sure that it does not affect the results. [23] In this work, GRI-Mech 3.0 1 [24] is used within CHEMKIN 2 [25] subroutines to model combustion without nitrogen chemistry. This mechanism includes 36 species and 219 reactions. Mixture averaged diffusion is used to calculate diffusion velocities: where !" is the mixture-averaged diffusion coefficient of the species .

Formulation and Numerical Method
A standard solution approach for such flow problems is to use staggered grids in the context of finitevolume methods [26]. Therefore, a finite volume formulation was derived for the one-dimensional total mass, momentum, energy, and species equations. The diffusive terms were discretized using a secondorder central difference scheme while the convective terms were discretized using the exact exponential solution to approximate the values between the centers of two control volumes [27]. Once the coupled, nonlinear system of equations had been discretized, the system was solved iteratively using Newton's method at each time step and continues until the maximum norm of the solution vector was reduced to within a specified tolerance. Since the convergence of Newton's method requires a sufficiently good initial guess and many convergence problems with this method can be eliminated by applying time-stepping, a time-dependent approach is used to solve this one-dimensional flame problem instead of solving a steady-state boundary value problem directly. This method has been chosen mainly because it is robust and it converges for sufficiently small time steps. Newton's method has been used to solve the system of nonlinear governing equations for each time step. Typically, the relative (scaled) tolerance should be in the range of 10 -3 to 10 -4 [28]. The linear system formed in Newton iteration was solved using a Bi-CGSTAB linear solver [29].
Since the problem was solved with dimensional variables, as opposed to their non-dimensional counterparts, it is necessary to discuss the actual physical size of the computational domain. The naturally occurring length scale of this problem is the thickness of the flame, which is expected to be of the order of 3 mm. Consequently, the calculated results could be expected to depend on the radius of curvature of the flame if it were of the same order of magnitude as the flame thickness. To avoid these size dependencies, the inflow boundary was set to be 60 mm such that the radius of curvature of the flame was at least 20 times that of its thickness.

Steady State Condition
For ignition, the reactant mixture at 300 K was introduced into pure nitrogen at 1800 K, which initially filled the computational domain. The temperature and species fields evolve until they reach steady state.
It is worth pointing out that this stationary flame was solved by time stepping of the transient model with steady boundary conditions. Figure 2 shows the species mass fractions and temperature profile of the stationary flame at the stoichiometric condition. To ensure the final product condition was met, the solution was tested for domain length independence by doubling its size and verifying solution consistency. In order to achieve a mesh independent solution and also ensure that the resolution is sufficient to capture the thin reaction layers of minor species, a non-uniform grid was used with spacings varying from 0.02 cm to 0.005 cm. This stationary flame was used as the initial condition of the timedependent boundary condition flame and also to compare some of its characteristics with expanding and contracting flames at various frequencies.
Prior to presenting the periodic results, steady computational results were verified against appropriate conditions in the literature. Currently, there are no experimental results for this exact geometry, so the results were compared to those of others in terms of the flame velocity, temperature, and species mass fractions. The first test was a comparison with the equilibrium state temperature and composition. The Chemical Equilibrium Calculation Spreadsheet [30] which invokes JANAF thermodynamical tables [31] was used to estimate the equilibrium state, and compared in Table 1 to the burned gas composition calculated herein at the end of the computational domain where all gradients were zero, and showed very good agreement. As discussed previously, even though this stationary flame has a finite curvature, it is not stretched. To Also shown in Fig. 3 as horizontal dashed lines, are the estimated planar unstretched laminar burning velocities of a stoichiometric and a fuel lean methane-air flame at atmospheric pressure reported in [32] which has been selected here for comparison as it is among the more classical studies in flame velocity measurements. There is excellent agreement between the experiments and the current simulation, demonstrating that the steady state solution computed here captures the basic physical processes of a methane-air flame.

4.2.Transient Periodic Flow
For transient flow cases, the inlet boundary condition of mass flow per unit axial length changes sinusoidally, given by: In this equation, B was set equal to 10% of ′ ! ( ′ ! = 0.2 / . ) and represents the amplitude of the oscillation and is the frequency, which is varied between 5 Hz and 2000 Hz. This lower limit was needed because the effects of unsteadiness for low equivalence ratio mixtures were observed even at these conditions. The reason to increase the frequency up to 2000 Hz is for the relevance to turbulent flows. For example, in the middle of a premixed turbulent duct flame stabilized by a backward facing step, mean frequency of the flame motion across a fixed point was observed to be in this range [33].
Experiments show that wrinkling process of flamelets depends on the geometry of turbulent flame [34].
Therefore, it has been suggested to study premixed turbulent flames in different geometries (i.e., where ′ is referred to here as the energy index (per unit axial length), and ′ represents the volume of cell (per unit axial length). Also, by comparing this quantity at the same point in consecutive cycles, it was used to determine if the solution has reached the fully-periodic state from its initial conditions.
In order to have a better comparison between different equivalence ratios, the flame speeds ( ! ) and the steady values of ′ at the mean mass flow rate (steady state case which has been used as the reference conditions) are summarized in Table 2: where is the frequency. As illustrated in Fig. 4, each data set decays to negligible variation after 9-12 cycles. Based upon these results, the comparisons have been made using the computational data in the tenth cycle.
In Table 3, ′ has been calculated and compared for different conditions; namely, steady state, expanding, and contracting flames, when the value of instantaneous mass flow rate was equal to that of the steady state ( ′ has the units J s per unit axial length, but has been normalized here by the steady state value). Since ′ is an integral value of heat release due to chemical reaction over the whole domain, it is related to the amount of mass entering the flame area. It has been discussed previously that due to the difficulty and ambiguity of defining a leading flame surface to relate the ′ to flame radius, the frequency response of the flame has been plotted versus the mass flow rate, and comparisons have been made using these integral values at the same mass flow rate, rather than at the same flame radius. At each frequency, the contracting flame has a value higher than unity and the expanding flame has a value lower than unity. This result is expected since the expanding flame has a positive stretch and the contracting flame has a negative stretch. The essential point being that at these instants for the two separate equivalence ratios their imposed inlet upstream hydrodynamic state is identical, but extra information is required regarding whether this instant is part of an expanding, contracting, or stationary flame in order to accurately estimate the instantaneous rate of energy conversion. While the differences between stationary, expanding, and contracting at this instant appear to be modest (+/-1 to 3% for these frequencies, which is linked to the relatively small fluctuations in mass flow rate), it is interesting to note that this characteristic is not monotonic with frequency for either equivalence ratio. Therefore, it is worth following the flame response throughout a complete cycle. The small variations observed in Table   3 are mainly because of two reasons. First, the methane-air flame with Le~1 is generally stretch insensitive suggesting that the transient flame responses should not deviate substantially from the instantaneous response. It should be noted that this mixture has been selected to keep consistency with similar works in other configurations. Second, in order to stay relatively close to the linear response limit, the amplitude change of mass flow rate is kept small at +/-10% throughout this study. With the quasi-steady assumption (solid line), changing the mass flow rate sinusoidally with time at the inlet affects the value of EI in the flame. Therefore, the flame responds to the change in mass flow rate instantaneously with no phase lag or amplitude change. In this regard, each point on the solid line corresponds to its own steady radius. As it can be seen in Fig. 5-a, the amplitude of the oscillations of ′ diminish and its phase shift increase with an increase in frequency to as much as ~π/4 at 2000 Hz.
Although the results of premixed methane-air flames can be compared qualitatively to those of the counter-flow configuration, the significant distinguishing feature of this study compared to previous works is the unstrained nature of the flame in the new geometry. Therefore, the present work helps to fill In this equation, is the frequency of the mass flow rate at the inlet boundary, is the thermal diffusivity and ! ! is the unstretched laminar flame velocity. This definition has the advantage of being evaluated on quantities that are specified (imposed) on the flow, or known a priori, such as the material properties. The Damköhler number could also be defined by considering the flame's response to the changing hydrodynamics. In the previous section, with increasing frequency it was observed that ′ had a diminishing amplitude and increasing phase lag, which aligns with observations of diminishing flame movement and increasing phase lag [3,12,13], meaning that flame response may be important to consider when comparing flames of different equivalence ratios. Essentially, variations in the flow were imposed on the flame, but the stretch rate experienced by the flames in terms of their changing positions and thermo-chemical structures are responding variables. Therefore, a second Damköhler number is also considered for which the flow time scale is based on the variation in stretch rate that includes the frequency and flame response, though not explicitly altering the chemical time scale due to varying chemical structures in the flame.
To calculate a characteristic stretch rate (stretch rate amplitude defined as the difference between maximum and minimum stretch rate in a cycle) in each case, equation (2) was used. As mentioned before, due to the particular geometry of the flame, the stretch rate equation is simplified to = ! ∇. !!! . Considering flame stretch rate, the flame displacement speed is used to analyze the dynamics of premixed flames. This quantity is generally ambiguous except for with a steadily propagating planar flame. The main reason is that the mass flow rate through the combustion region varies with distance through the flame. This variation makes it difficult to select an iso-surface to represent the flame surface. In the present work, as ! represents local flame speed, selecting an isosurface is challenging. Therefore, this quantity is sensitive to iso-surface selection. In order to apply this equation to a flame, an arbitrary isotherm needs to be selected where the velocity and divergent to the normal of the flame front is calculated. As suggested in [18], an iso-surface close to the hot side of the flame is more insensitive to the magnitude of the calculated stretch rate. Therefore, the selected isotherm is where temperature rise reaches 90% of the total temperature increase. Fig. 8 illustrates how the radius of this isotherm and its corresponding stretch rate change with time in an oscillating flame with equivalence ratio of 0.8 and frequency of 20 Hz. Comparing the two plots, there are two locations that flame stretch rate becomes zero (where the flame stops moving and changes direction); maximum and minimum stretch rates occur when the flame is expanding and contracting while passing through the mid-point (mean mass flow rate), respectively. Fig. 9 shows the stretch rate amplitudes that characterize the flame responses for the different equivalence ratios over the range of frequencies tested, and shows the connection between the two time scales. As expected for a constant equivalence ratio, as the frequency rises, the stretch rate amplitude increases due to the higher time rate of change of the iso-surface radius, though not linearly because of the lagging response of the flames. Also, at the same frequency, lower equivalence ratio flames have smaller stretch rate amplitude. The main reason is that the oscillating flame amplitude is smaller in lower equivalence ratios due to the greater phase lag discussed earlier. Therefore, the amplitude of the iso-surface becomes smaller with a constant rate of change of radius. the Damköhler Number based on the stretch rate is smaller than that based on the frequency, while the opposite is true for the higher frequencies when the flame cannot respond to such rapid changes. The characteristic stretch rate considered in ! calculations is defined as the difference between maximum and minimum stretch rates (stretch rate amplitude). As illustrated in Fig. 9, the stretch rate amplitude increases with frequency at each equivalence ratio in the same order of magnitude. This increase suggests that at the limit of very high frequencies, the stretch rate amplitude is also increased to a range that makes ! very small (close to zero). Due to the similarity of the two plots, the characteristic Damköhler Number based on imposed frequency was used in the remainder of this study.  Figure 12 shows the normalized EI amplitude at the mean mass flow rate (i.e., minor axis of ellipses in Fig. 7), plotted with respect to Damköhler number (Eq. 12). As can be observed, there is a maximum point in all three plots. As the equivalence ratio decreases, this maxima move to higher Damköhler Numbers (i.e., lower frequencies).
In this figure, ′ amplitude at the mean mass flow rate is the difference between maximum and minimum values of ′ at each frequency, which is then normalized by the ′ amplitude in the quasisteady condition. The responding behaviour of the different equivalence ratio flames is not well captured by the Damköhler Number, with the data for the equivalence ratio of 0.7 being dramatically different from the other two cases.

Conclusions
A model was developed in order to revisit the quasi-steady assumptions that are used in current approaches of premixed laminar flamelet models. This model is based on a cylindrically-symmetric radial outward flow geometry with an applied sinusoidal variation of mass flow rate that was only 10% It has been concluded that similar to planar flames exposed to oscillating strain rates where the flame stretch always has a positive value, the transient response of laminar premixed flames results in decreasing amplitudes of motion and burning rates and increasing phase lag with increasing frequency.
This changing ′amplitude and phase lag resulted in a elliptical shape when it was plotted with respect to the instantaneous mass flow rate, which showed the different burning rates between positively and negatively stretched flames with the same imposed hydrodynamic state, as well as how that differed from the quasi-steady state.
A comparison was also made between stoichiometric and two fuel lean flames to study the effect of equivalence ratio on the transient response of laminar premixed flames to various frequencies. It was concluded that as the equivalence ratio was made leaner, so that the chemical time scales became longer,      Normalized EI amplitude at mean mass flow rate Da phi=0.8 Fig. 11. Normalized EI amplitude with respect to Damköhler number ( ! ) for three equivalence ratios