Effect of small heat release and viscosity on thermal-diffusive instability

The linear stability of a thermal reaction front has been investigated based on the thermal-diffusive model proposed by Zel’dovich and Frank-Kamenetskii, which is called ZFK model. In the framework of ZFK model, heat-conduction and mass-diffusion equations are treated without the effect of hydrodynamic flow. Then, two types of instability appear: cellular and oscillatory instabilities. The cellular instability has only positive real part of growth rate, while the oscillatory instability is accompanied with non-zero imaginary part. In this study, the effect of heat release and viscosity on both instabilities is investigated asymptotically and numerically. This is achieved by coupling mass-conservation and Navier–Stokes equations with the ZFK model for small heat release. Then, the stable range of Lewis number, where the real part of growth rate is negative, is widened by non-zero values of heat release for any wavenumber. The increase of Prandtl number also brings the stabilization effect on the oscillatory instability. However, as for the cellular instability, the viscosity leads to the destabilization effect for small wavenumbers, opposed to its stabilization effect for moderate values of wavenumber. Under the limit of small wavenumber, the property of viscosity becomes clear in view of cut-off wavenumber, which makes the real part of growth rate zero.

www.nature.com/scientificreports/ flame front is treated as a zero-thickness interface and the basic flow, which is assumed to be steady one-dimensional flow, is described by the constant density and velocity, which are discontinuous across the front. Such a discontinuity is the intrinsic property of hydrodynamic instability, or the Darrieus-Landau instability (DLI): the DLI claims that the growth rate of perturbations is an linear function of wavenumber, whose coefficient is positive due to the heat release, and a planar flame front is unstable for any wavenumber. As a consequence of the jump of normal velocity across a flame front, the tangential components of velocity field induced by small perturbations are also discontinuous. This fact is essential for the existence of DLI. However, our study is performed in the scale of thermal thickness of a flame front, or in the framework of thermal-diffusive model. Then, the boundary conditions across the reaction front, which have been derived based on the large activation energy asymptotics, are described by the continuity of all quantities except for the pressure, whose discontinuity is supported by the viscosity. Due to the continuity of normal velocity, the tangential components of velocity becomes continuous for small perturbations. For a basic flow, the reaction front is considered as a density continuous interface, whereas the derivative of density is discontinuous. As a result, the DLI does not appear in our treatment.
In this study, we consider the case of small heat release. Then, we need to deal with not only the thermal-diffusive model but also hydrodynamic one. So far, the coupling of these two models has been achieved by the multi-scale analysis based on the thermal thickness of a flame front, or by introducing the preheat zone inside the front [23][24][25][26] . This approach successfully combines the DLI [27][28][29] with the Markstein effect which reflects the thermal and mass diffusivities [30][31][32] . The numerical study to understand the dynamics of a flame front has been performed for cellular and non-cellular fronts 33,34 and for cellular front with heat loss [35][36][37] . Despite of many works on the instability of a flame front, the influence of heat release and viscosity on the stability boundaries of cellular and oscillatory instabilities has not ever been declared clearly. Besides, the effect of higher order terms of wavenumber, such as third and forth order terms, on the growth rate has not been discussed in detail. The existence of third order term of wavenumber has the possibility of leading to two cut-off wavenumbers, opposed to the case of zero heat release, where the cut-off wavenumber is determined uniquely: the cut-off wavenumber means the value of wavenumber at which the growth rate of perturbations changes its sign. In order to reveal these influence, we directly examine the linear stability of reaction front for infinitesimal perturbations with the heat release being a small parameter. Besides, the calculation of numerical solutions is performed based on the TVD-MacCormack scheme [38][39][40][41][42][43][44] .

Methods
We employ the following non-dimensional governing equations 13,25 where the non-dimensional variables, ρ , V = (U, V , W) , T, P M and Y are the density of mixture, velocity field, temperature, pressure and mass fraction of a deficient species which is consumed through the exothermic chemical reaction. The parameters γ , Ma, Pr, q and Le are the specific heat ratio, Mach number, Prandtl number, heat release and Lewis number, respectively. The dependent variables are made dimensionless by use of those of fresh mixture at a position far from a flame front, where the flow field is assumed to be steady planar, denoted by ρ −∞ , S L , T −∞ , P −∞ and Ỹ −∞ . Especially, S L is called the laminar flame speed, or the burning velocity, which is the propagation velocity of a planar flame front into a quiescent fuel mixture 13 . Hereafter, the variables with tilde symbol stand for dimensional ones.
The parameters l d , D th , D , μ , c p and ˜ are the thermal thickness of a flame front, thermal diffusivity, diffusion coefficient of a deficient species, coefficient of viscosity, specific heat at constant pressure and thermal conductivity. We assume that D , μ , c p and ˜ are all constant. The adiabatic sound speed is defined by c s = (γP −∞ /ρ −∞ ) 1/2 . We normalize the coordinates x and the time t , by use of l d and S L , as Boundary conditions. In the far field from a flame front on the unburned side, where the flow is assumed to be steady planar, all quantities become unity due to the definition (6). www.nature.com/scientificreports/ The deficient component of combustible mixture is assumed to be depleted during the exothermic chemical reaction inside the reaction front. Then, the mass fraction Y is equivalently zero on the burned side.
In this study, we employ the zero-Mach-number approximation, Ma ≪ 1 . Under this limit, as readily found from (2), the pressure is assumed to have the following asymptotic form 22 .
In the subsequent discussion, we consider the variation of pressure at O(Ma 2 ). In order to know the profile of physical quantities outside the reaction front, we need the boundary conditions on it. In the limit of large activation energy, ε(= q/β) ≪ 1 , the following jump conditions across the reaction front have been derived 15,16 . where the brackets for arbitrary function �(x, y, ξ , t) means the difference of its values on the burned and unburned sides, In the subsequent calculations, the boundary conditions (20), (21), (23)- (27) are imposed to determine the integral constants of unperturbed and perturbed solutions of (1)-(5) for a weakly corrugated front, whose amplitude is infinitesimally small.
Basic flow. As a basic flow without perturbations, we consider a steady planar flow. The notation of over-bar is used to denote this flow field. For example, the velocity of basic flow is denoted by W =W(ξ ) . In this case, the jump conditions (24)-(26) are given by We note that (23) and (27) are automatically satisfied by the solutions obtained by use of (28). For a steady planar flow, the mass conservation Eq. (11) becomes From (29), by use of (20) and (28), we readily find that the steady planar mass flux is unity on both sides of reaction front.

Formulation of perturbations.
We consider the situation where the reaction front is slightly deviated from its unperturbed state. We assume that the heat release is sufficiently small.
Based on (36), the amplitude of weakly corrugated front is given in the form of harmonic perturbation by where ω is the growth rate and the wavenumber k is defined by k = (k 2 x + k 2 y ) 1/2 . The non-dimensionalization is made as ω =ωl d /S L and k =kl d with ω and k dimensional quantities.
We impose small perturbations on the basic flow as follows.
where the order of each perturbation is assumed to be Upon (37)-(40), the perturbed governing equations have the following form with O(q 4 ) terms omitted due to (36).
As is readily confirmed from (35), the derivatives of ρ , P , W and T are O(q) and that of Ȳ is O (1). Unlike the ordinary thermal-diffusive model, the influence of flow field appears in (45) due to the existence of small heat release. In this case, the coupling of hydrodynamic equations and thermal-diffusive ones is possible because the heat equation remains unchanged from the ordinary one and solvable. The detail of calculation is given in the succeeding subsection.

Calculation of perturbations.
We calculate perturbations (39) and (40) by solving perturbed governing Eqs. (41)- (46). From (46), we find that once the perturbation of temperature is calculated, that of density is obtained.
Assuming Re[ω] > 0 , the requirement that perturbations are bounded in space leads to the solution of (48) as where We note that the signs are in the same order. The integral constants θ m and θ p are determined in the succeeding subsection.
Next, in order to calculate P p , we transform the perturbed mass conservation Eq. (41) into where (35) and (46) are used and higher order terms with respect to q are omitted due to (36). Then, taking the transverse derivative ∇ t · of (42) and differentiating (43) with respect to ξ , and combining them with the help of (51), we obtain By use of (35), (37) and (39), the combined equation (52) becomes the following differential equation for P p .
where and u m , u p , v m , v p , w m and w p are integral constants.
The perturbation of mass fraction is calculated by integrating (45) by use of (35), (37), (40), (47), (49) and (57). www.nature.com/scientificreports/ where and c m is an integral constant. Unlike the classical thermal-diffusive model, the solution of perturbed mass fraction (59) includes the terms reflecting the effect of hydrodynamic flow due to non-zero values of heat release q. Thanks to these terms, it becomes possible to incorporate the influence of heat release and viscosity into the thermal-diffusive instability. At last, by substituting solutions (49), (55)-(57) into the mass conservation equation (51), we obtain the following relations.
Dispersion relation. In this subsection, we calculate the dispersion relation of thermal-diffusive instability, which represents the relation between the growth rate and wavenumber for perturbations, by substituting the solutions obtained above into the boundary conditions across the reaction front. For the perturbed quantities, the boundary conditions (23)- (27) are written down as We note that the continuity condition of mass flux in (24) leads to the continuity of density perturbation due to (64), but this gives no information. In fact, this is the same condition as the first one in (66) because of perturbed equation of state (46). From (61)-(67), we determine integral constants of perturbations. After that, (68) is used to calculate the dispersion relation of thermal-diffusive instability including the effect of small heat release accompanied with the viscosity.
At first, the integral constants for thermal-diffusive quantities are easily determined. Indeed, substituting (49) into (66), we obtain Then, integral constants θ m and θ p of temperature perturbation are determined by solving (69) as The expression of c m is obtained from (67), by use of (59), as If q = 0 , the above solutions (70) and (71) coincides with those in the classical thermal-diffusive model 22,32 . Next, the integral constants for hydrodynamic flow field are determined. The substitution of (55) and (56) into (63) gives  (61) and (72), u m , u p , v m and v p are eliminated to give the following relation.
The continuity of normal velocity (64), with (57) substituted, leads to Solving (73) and (74), we find By use of solutions (54) and (57) The numerical and analytical studies of (81) are given below.

Asymptotic results
The case of q = 0. In the pure thermal-diffusive model, or the ZFK model 1,2 , the effect of hydrodynamic flow was effectively decoupled from the heat and mass-diffusion equations. This situation is easily recovered by taking q = 0 in (81). It may be helpful to recall the original dispersion relation for thermal-diffusive instability 3,32 , which is given by www.nature.com/scientificreports/ In addition, resorting to the assumption that the Lewis number is nearly equal to unity, which is acceptable for real gases 6,13 , we introduce O(1) constant le as under the limit of large Zel' dovich number, β ≫ 1 . By use of (83), the dispersion relation (82) is transformed into the following form.
The numerical calculation of (84) has been shown and summarized by many works 3,22,32 . The cellular instability is observed for some range of Lewis number less than unity and the oscillatory one appears for large values of Lewis number. The cellular instability is characterized by the positive values of real part of growth rate without imaginary part and the oscillatory one is by those with imaginary part. This tendency remains unchanged for the case of non-zero values of heat release as discussed below.
Computation of dispersion relation. We consider the case of q = 0 . In order to compare our result with the previous research, we apply the assumption (83) to the dispersion relation (81) under β ≫ 1 . Then, we obtain the reduced dispersion relation, with qβ = O(1) and O(q) terms omitted, as follows.  . The stable range of Lewis number, where the real part of growth rate takes negative values, is greatly altered from that in the previous case of q = 0 due to the heat release accompanied with the viscous effect. At first, from Figs. 1 and 3, the inclusion of heat release is found to extend the stable domain for both cellular and oscillatory instabilities. In addition, Fig. 4 shows that the non-zero values of Prandtl number, or the viscous effect, also broadens the stable range of Lewis number for the oscillatory instability. Contrary to such stabilization effect, it should be emphasized that the viscosity has the destabilizing effect on the cellular instability if the wavenumber is small, as shown in Fig. 2. For some moderate values of wavenumber, the cellular instability is stabilized by the viscosity.
Cellular instability for small wavenumber and growth rate. As discussed above, we find that the viscosity brings the destabilization effect on the cellular instability for small values of wavenumber, opposed to the stabilizing one for finite values. Therefore, it may be meaningful to investigate the dispersion relation (81) in detail for small wavenumbers to reveal the destabilization effect of viscosity. For this purpose, we consider two cases: ω ∼ k 2 ≪ 1 and ω ∼ k 4 ≪ 1 . Especially, in the latter case, our result contains the new effect of third order term with respect to wavenumber in the reduced dispersion relation.
The case of ω ∼ k 2 ≪ 1. At first, the effect of wavenumber is examined up to second order term, which corresponds to the Markstein effect reflecting the curvature effect 30,31,45 . By assuming that ω ∼ k 2 ≪ 1 , the dispersion relation (81) is reduced to For q = 0 , the above Eq. (86) is consistent with that obtained in the previous work 45 . In addition, assuming (83) for β ≫ 1 with qβ = O(1) , the reduced dispersion relation (86) gives the following growth rate.  www.nature.com/scientificreports/ The growth rate (87) clearly shows how the variation of le in the range of small wavenumber influences the cellular instability with small heat release. The stability condition of (87), or ω < 0 , is given by le > −2 − qβ(5 + 3Pr)/(2 + 2Pr) . The stability boundary le = −2 , which is obtained in the previous research 3,32,45 , is extended by the inclusion of heat release q accompanied with the viscous effect represented by Prandtl number Pr. The increase of q leads to the relaxation of lower bound of le. On the other hand, the rise of Pr narrows the stable range of le. This result explains the destabilizing effect of viscosity for small wavenumber as implied by the numerical computation in Fig. 2.
The case of ω ∼ k 4 ≪ 1. Next, we proceed to higher order terms of wavenumber, O(k 3 ) and O(k 4 ) terms. We note that O(k 3 ) term is intrinsic to the effect of heat release and not included in the original dispersion relation of thermal-diffusive model 3,32 . In the present case, the cut-off wavenumber is no longer determined uniquely and there is the possibility of two cut-off wavenumbers depending on values of Lewis number.
Based on the assumption of ω ∼ k 4 ≪ 1 , the dispersion relation (81) is transformed into Solving the reduced dispersion relation (88), by use of the assumption (83) for β ≫ 1 with qβ = O(1) , we find with We note that the Prandtl number Pr typically takes positive and small values: for example, Pr ≈ 3/4 in air 13 . For q = 0 , the O(k 4 ) term in (89) is consistent with previous work 45 , which treated the case of le ≈ −2 , that is ω 4 k 4 /2 ≈ −4k 4 . For non-zero values of heat release, as confirmed from ω 2 and ω 3 , the rise of Prandtl number increases the growth rate at second order of wavenumber and decreases it at third order. Because the wavenumber is assumed to be small, the viscous effect is not so important at the forth order of wavenumber. By setting ω = 0 in (89), the cut-off wavenumber is obtained as In the case of q = 0 , the cut-off wavenumber is uniquely determined as k c = √ (le + 2)/(3le − 2) . The set up of cellular instability is given by le = −2 or 2/3. The positive value of growth rate for small wavenumbers, or long wavelengths, is observed for le < −2 and that for large wavenumber, or short wavelength, is observed for le > 2/3 . For −2 < le < 2/3 , the growth rate is negative for arbitrary wavenumber.
Although (90) implies the existence of instability at some finite value of wavenumber for le > 2/3 , it should be omitted because the wavenumber is now assumed to be small. In fact, even if the value of le is sufficiently large, k c is not so small: k c → 1/ √ 3 + 0 as le → +∞. The new tendency of cellular instability due to the existence of heat release is characterized by ω 3 in (89). This term promotes the cellular instability accompanied with two cut-off wavenumbers for some range of le. The stability boundary of le for the cellular instability is obtained from (90) as follows. , (91) le 2 < le < le 3 (no cut-off wavenumber), (92) le < le 1 , le > le 4 (one cut-off wavenumber), (93) le 1 < le < le 2 , le 3 < le < le 4 (two cut-off wavenumbers), www.nature.com/scientificreports/ where le 2 and le 3 are the solutions of ω 2 3 − 4ω 2 ω 4 = 0 and determine the cut-off wavenumber which induces the cellular instability. We note that le 1 = le 2 = −2 and le 3 = le 4 = 2/3 for q = 0.
(95)   www.nature.com/scientificreports/ In the stable range of le given by (91), the dispersion relation (89) shows negative values of ω for arbitrary wavenumber and the perturbations disappear in time. On the other hand, in the case of (92), the positive value of growth rate for small wavenumbers is observed for le < le 1 and that for large wavenumbers is observed for le > le 4 . Besides, the new feature of instability due to the existence of heat release is observed in the range of (93). The increase of Pr narrows the range of (93), where two cut-off wavenumbers exist, and decreases the values of k c . The cases of le > le 4 and le 3 < le < le 4 are omitted because of the assumption k 4 ≪ 1 . In these cases, the cut-off wavenumbers does not take small values.
The growth rate in the case of (93) is shown in Figs. 5 and 6. The existence of viscosity raises the value of le 2 , which implies the Lewis number approaches unity. Therefore, the cellular instability tends to be promoted by the viscosity for small wavenumbers. This fact is consistent with the result obtained numerically in Fig. 2. In the case of (92), it is clear that the viscous effect has two aspects on the instability as plotted in Fig. 7. Although the viscosity enhances the growth rate for small wavenumbers, it shrinks the range of wavenumber where the instability exists. For k ≈ 0 , by setting ω = 0 in (89), the stability boundary with respect to the values of Lewis number and Prandtl number is given by ω 2 = 0 and shown in Fig. 8. The destabilization due to the viscosity for small wavenumbers is easily comprehended from Fig. 8. Finally, the dependence of cut-off wavenumber on the Prandtl number is shown in Fig. 9. It is apparent that, among two cut-off wavenumbers, the smaller one is reduced due to the viscosity. In contrast, as for the larger cut-off wavenumber, the rise of Pr brings either increase or decrease of it subject to the values of le.

Numerical study
Discretization. In this section, the numerical solutions of governing Eqs. (1)-(5) are computed by use of the TVD-MacCormack scheme [38][39][40][41][42][43][44] . For simplicity, the two-dimensional flow field in x-z plane is considered, where the reaction front propagates in the negative z-direction. In order to employ the MacCormack scheme 39,40 , (1)-(5) are rewritten into the conservation form: where the flux vectors are  www.nature.com/scientificreports/ Recalling (8) and (9), the coefficient of reaction term is expressed by A = q −2 β 2 � exp (1 + q −1 )β . In addition, utilizing the result from burning-rate eigenvalue problem under large-activation-energy asymptotics 3,12,13 , we may estimate � ≈ q 2 T b /(2Le) . Therefore, the coefficient A is approximated by In this study, the time-splitting method 38 is used to calculate numerical solutions of (98). Then, (98) is divided into two one-dimensional problems corresponding to x-and z-directions: The one-dimensional explicit MacCormack scheme is applied to solve each equation in (102), respectively. The MacCormack scheme has second-order accuracy in both time and space. The procedure of discretization is expressed by finite-difference operators L x (�t) and L z (�t) [46][47][48] . Each operator consists of two steps: the predictor and corrector steps. For each equation in (102), the predictor step is and the corrector step is where E p and F p are calculated by use of U p . The superscript n implies the time step, and the subscripts i and j imply the spatial steps in x-and z-directions, respectively.
Each component of flux vectors is discretized in space as follows.
(103) www.nature.com/scientificreports/ For a flux vector E , the backward and central differences are applied to the partial derivatives with respect to x and z, respectively. On the other hand, for F , the central and backward differences are applied to the partial derivatives with respect to x and z, respectively.
To avoid the oscillation of numerical solutions, we introduce a total variation diminishing (TVD) numerical scheme in the corrector step [41][42][43][44] . We employ the Causon's model 44  where (6) and (7) are used. Accounting for the above formulae, the finite-difference operators of TVD-MacCormack scheme are given by (106) L x (�t)U n : Settings. The initial conditions are given as follows. The development of perturbed reaction front is examined by initially imposing small perturbations on a planar front. We use the sinusoidal displacement periodic in x-direction, which is given by where a is the initial amplitude of displacement and k = 2π/ is a wavenumber with the wavelength of perturbations. The initial solutions are given by those of a steady planar flow in (35). Therefore, the initial values of a flux vector U consists of the following components.
The perturbations are imposed on the planar reaction front initially located at z = 9.
As for the boundary conditions, the spatially periodic conditions are employed in the x-direction. In the z-direction, the far-field condition (20) is used at the upstream edge and zero-gradient conditions are made at both upstream and downstream edges.
We consider the case of k = 0.5 , correspondingly = 4π . The parameters are set as γ = 1.4 , Ma = 0.01 , β = 10 and a = 0.1 . The Lewis number Le is calculated by use of (83). The computational domain is arranged for one wavelength of a perturbed front, , in x-direction. In the propagation direction, we consider 0 ≤ z ≤ 10.
We use the Grid Convergence Index (GCI) method for estimation of discretization error [51][52][53] . A representative grid size h is given by h = (�x�z) 1/2 because we consider the case of uniform grid. We select three different sets of grid so as to h 1 < h 2 < h 3 . The grid refinement factor is defined by the ratio of grid size, which should be greater than 1.3. Let r 21 = h 2 /h 1 and r 32 = h 3 /h 2 . Then, the apparent order p of the numerical method is calculated by where (115) L z (�t)U n : (117) d = a sin(kx), (120) p = 1 ln(r 21 ) ln(ε 32 /ε 21 ) + q(p) ,  www.nature.com/scientificreports/ In this study, the key variable φ of simulation is set as ρY , which is useful to identify the location of reaction front because Y = 0 on the burned side. The percentage of oscillatory convergence is identified by the number of s < 0 among entire lattice points. The approximate relative error e 21 a and the grid-fine convergence index GCI 21 fine are given by The numerical uncertainty at each lattice point is estimated by use of GCI 21 fine (p ave ) and indicated by error bars in plotting fine-grid solution φ 1 , where p ave is an average value of p and measures the global order of accuracy. Table 1. For 65 × 513 and 65 × 1025 grids, the finegrid solution φ 1 is plotted in Figs. 10 and 11 with error bars calculated by GCI 21 fine (p ave ) . In view of GCI, a 65 × 1025 grid looks like better than others, though the oscillatory convergence occurs at some points and the average apparent order is not near 2, which is a theoretically desirable order of accuracy in the case of MacCormack scheme. By adopting a 65 × 1025 grid as the spatial resolution, the increments in x-and z-directions are �x = 4π/64 ≈ 0.2 and �z = 10/1024 ≈ 0.01 . The time increment is set as t = 3 × 10 −5 . The temperature distribution is shown in Figs. 12 and 13 for two values of Lewis number, which are chosen so as to show both unstable and stable cases as guessed from Fig. 1. It is possible to estimate the numerical growth rate by considering the front displacement A(t), which is calculated from the difference of maximum and minimum positions of the front. The front is identified by an isoline ρY = 0.1 . Figures 14 and 15 show the calculated front displacement in time and the estimated exponential curves, −0.005 + 0.195 exp(−0.035t) and −0.01 + 0.195 exp(0.055t) , which mean the numerical growth rates are −0.035 for le = −4 and 0.055 for le = −6 . There is a concern on the decrease of A(t) at the early stage for le = −6 . This might indicate that other grids, which would be finer in z-direction, are demanded. Besides, the case of small wavenumber has not been examined yet, which should be compared with the analytical result.

Discussion
We remark that our results do not contain the Darrieus-Landau instability (DLI), or the hydrodynamic instability. This is because we have employed the boundary conditions across the reaction front, which have been derived systematically under the large activation energy asymptotics, and focused on the pure thermal-diffusive instability. In our case, the density, or the normal velocity, is continuous across the reaction front even for a basic flow, whereas their gradients are discontinuous. Then, the x-y components of velocity field are also continuous across the front, opposed to their discontinuity in the hydrodynamic model. The discontinuity of density is essential for the existence of DLI. The boundary conditions across a flame front, which reflect the density jump, have been studied in detail based on the multi-scale analysis by introducing the preheat zone inside a flame front [23][24][25][26] .
We roughly examine why the DLI drops out in the present study by removing the continuity condition of density, or that of normal velocity, which is expected to lead to the appearance of O(k) term in the dispersion relation corresponding to the DLI. In order to confirm this, we consider the following solutions of a basic flow, instead of (35), in accordance with those of hydrodynamic model 27,28 .
In this case, boundary conditions (62) and (63) are changed into Figure 11. The same as Fig. 10  www.nature.com/scientificreports/ We employ solutions (47), (49), (54)-(57) and (59), which remain unchanged in order to couple the thermaldiffusive instability with the hydrodynamic one. Substituting these solutions into boundary conditions (64)-(68), (124) and (125), the dispersion relation is calculated, under the assumption of ω ∼ k 2 ≪ 1 and with O(q 2 ) terms omitted, as www.nature.com/scientificreports/ The dispersion relation (126) gives the following growth rate for small wavenumbers, which possesses O(k) term, which reflects the influence of hydrodynamic instability. Although this result has been obtained in very rude manner, the intrinsic importance of density discontinuity for the existence of DLI may become apparent. Finally, we note that the present study is restricted to the case of continuous density across the reaction front. Therefore, the DLI does not appear.

Conclusions
We have investigated how the thermal-diffusive instability is affected by the heat release and viscosity. Instead of introducing the effect of thermal thickness of a flame front, or the preheat zone, we directly dealt with the governing equations in the framework of thermal-diffusive model, or ZFK model. This approach has a merit of (127) ω = k 2 − 1 β(1 + q) + 1 2Le + LePr 2 1 + LePr k 2 , Figure 13. The same as Fig. 12 but for le = −6. www.nature.com/scientificreports/ understanding the viscous effect on the stability boundary of cellular and oscillatory instabilities in a simple manner: the stability boundary is described in the plane with respect to Lewis number and wavenumber. Our main assumption is the small values of heat release with the continuity condition of density across the reaction front, which is valid in the framework of large activation energy asymptotics. Thanks to this assumption, the coupling of hydrodynamic and thermal-diffusive models becomes possible. As a result, the effect of viscosity is incorporated into the dispersion relation for the pure thermal-diffusive instability, though the DLI does not appear. Due to the existence of heat release and viscosity, the stability boundary of thermal-diffusive instability is largely changed from that obtained in the previous case of no heat release. The heat release has a stabilizing effect on both of cellular and oscillatory instabilities. In other words, for arbitrary wavenumber, the stable range of Lewis number, where the real part of growth rate of perturbations is negative, is widened as the heat release increases. Such a stabilizing effect is also brought by the rise of Prandtl number, which originates from the Navier-Stokes equation, on the oscillatory instability. However, the effect of viscosity has two opposite aspects on the cellular instability. For some moderate values of wavenumber, the stable range of Lewis number for the cellular instability is extended. On the other hand, for small wavenumbers, the stable range is narrowed by the increase of Prandtl number.
The growth rate of cellular instability has been studied in detail by restricting ourselves to small values of growth rate and wavenumber. In this case, we find that the third order term with respect to wavenumber appears in the reduced dispersion relation due to the existence of heat release. This term brings the new band of Lewis   Table 1. The resolution test is examined for several grids N x × N z , where N x and N z are the number of lattice points in x-and z-directions, respectively. All estimation is made at y = 2π and t = 0.075 for every lattice points in z-direction. www.nature.com/scientificreports/ number, where two cut-off wavenumbers exist: in other words, there are two wavenumbers at which the growth rate of perturbations changes its sign. The relation between the Lewis number and the cut-off wavenumber is explicitly derived with the viscous effect included. While the value of smaller cut-off wavenumber is reduced as the Prandtl number increases, the behavior of larger one shows both of the increase and decrease depending on values of Lewis number. The numerical calculation is also made by use of TVD-MacCormack scheme. The numerical growth rate is estimated for different values of Lewis number with small heat release. For moderate value of wavenumber, the numerical result seems to be consistent with the stability boundary obtained in the asymptotic analysis. However, the case of small wavenumbers has not been checked yet. We need to compare this case with the analytical result.

Data availibility
No datasets were generated or analysed during the current study.