© Hindawi Publishing Corp. STABILITY OF REACTION FRONTS IN THIN DOMAINS

The paper is devoted to the stability of reaction fronts in thin
domains. The influence of natural convection and of heat losses
through the walls of the reactor is studied numerically and
analytically. Critical conditions of stability of stationary
solutions are obtained.


Introduction.
The stability of reaction fronts is important for many technological devices and it is interesting from the point of view of nonlinear dynamics (see [6] and the references therein).It is known that thermal instability can be strongly influenced by heat loss [1] especially near the extinction limit.In this work, we study the influence of heat loss and of convection on thermal fronts in thin domains, where we can reduce the space dimension of the problem and in some cases simplify its analysis.
We consider the system of equations where the reaction-diffusion system is coupled with the Navier-Stokes equations under the Boussinesq approximation to describe natural convection which can occur because of the heat produced by the reaction: Here T is the temperature, α the depth of conversion (or the dimensionless concentration of the product of the reaction), v = (v x ,v y ,v z ) the velocity of the medium, p the pressure, κ the coefficient of thermal diffusivity, q the adiabatic heat release, ρ an average value of density, ν the coefficient of kinematic viscosity, g the gravity acceleration, β the coefficient of thermal expansion, γ the unit vector in the z-direction (upward), x, y, and z the spatial coordinates, (1.3) Here ∂T /∂n denotes the derivative in the direction of the outer normal vector.The boundary conditions at the lateral walls of the reactor correspond to heat loss through the boundaries.The free-surface boundary conditions for the velocity simplify the analysis of the problem.It corresponds to the case of the batch reactor.We will also consider the case where with a given velocity v 0 at the entrance and the exit of a continuous reactor.We introduce dimensionless spatial variables xc/κ, yc/κ, and zc/κ, time tc 2 /κ, velocity v/c, and pressure p/c 2 ρ.Here is the normal velocity of a condensed phase reaction front [4], where T b = T 0 + q is the adiabatic temperature.Denoting θ = (T − T b )/q and keeping for convenience the same notations for other variables, we rewrite system (1.1) in the form (1.6) Here P is the Prandtl number, The boundary conditions remain the same as above.
The contents of the paper are as follows.In the next section, we reduce the spatial dimension of the problem and introduce the heat loss in the heat equation instead of the boundary conditions.In Section 3, we formulate the interface problem using the infinitely narrow reaction zone method.Section 4 is devoted to the analysis of stationary thermal regimes and their stability.Finally in Section 5, we present numerical simulations of reaction fronts with and without convection.

Thin domains.
Without heat loss (σ = 0) and without convection (R = 0) problem (1.1), (1.3) has a one-dimensional (1D) stationary solution which depends on the zvariable only.The infinitely narrow reaction zone method (see Section 3) allows us to find it approximately analytically and then to use this explicit form of the solution to study its stability.If σ ≠ 0, the stationary solution is not one-dimensional any more and its analytical approximation cannot be found or becomes too complicated to be used.Therefore, the stability analysis cannot be carried out.To overcome this difficulty and to study reaction fronts with heat loss, we consider thin domains where the heat loss can be introduced in the equation instead of the boundary conditions.In this section, we suppose that the dimensionless domain is thin in the x-direction, that is, x ∈ ] − /2, /2[ with > 0 a small parameter.
Taking into account the dependence on , we rewrite system (1.6) in the form ) ) with ψ(θ) = Ze θ/(Z −1 +δθ) .We suppose that y ∈ ] − l y ,l y [, z ∈ ] − l z ,l z [.The boundary conditions are as follows (we suppose that the coefficient of heat loss at the boundary is of the order ): ) We make the usual change of variables x = x/ in (2.1), (2.2), (2.3), (2.4), (2.5), and (2.6) and we suppose the following formal expansions for the unknowns: Our goal is to prove, at least formally, that v 0 x = 0, θ 0 , α 0 , v 0 y , v 0 z , and p 0 are independent of x and that they satisfy the following system in (y, z) ∈ ]−l y ,l y [×]−l z ,l z [: with the boundary conditions (2.12) In this section we use the notation In order to obtain the limit problem, we substitute the above formal expansions in (2.1), (2.2), (2.3), (2.4), (2.5), and (2.6) and we equate the coefficients of k for any k ∈ Z.We first remark that the boundary conditions (2.12) are an immediate consequence of the corresponding boundary conditions from (2.5) and (2.6).
Equating the terms of the order O( −2 ) in (2.1) and O( −1 ) in the first equality of (2.5), we obtain which implies that θ 0 is independent of x.
Equating the terms of the order O( −1 ) in (2.1) and O( 0 ) in the first equality of (2.5), we obtain (2.14) Therefore, θ 1 is independent of x.
Equating the terms of the order O( −2 ) in the x-component of (2.3), we obtain Using also the condition v 0 x = 0 on x = ±1/2, we deduce Equating the terms of the order O( −2 ) in the y-component of (2.3), we obtain (2.17) Using also the condition ∂v 0 y /∂ x = 0 on x = ±1/2, we deduce that v 0 y is independent of x.
In the same manner, we obtain that v 0 z is independent of x.Equating the terms of the order O( 0 ) in (2.4), we obtain Equating the terms of the order O( 0 ) in (2.2), we obtain (2.9) with the help of (2.16) and (2.19).Taking into account the limit conditions, the function α 0 is the solution of a well-posed problem which does not depend on x, which implies that α 0 is independent of x.
Equating the terms of the order O( −1 ) in the x-component of (2.3), we deduce that p 0 is independent of x with the help of (2.16) and (2.19).
Equating the terms of the order O( 0 ) in (2.1), we obtain Equating the terms of the order O( ) in the first equality of (2.5), we obtain (2.21) Integrating (2.20) in x between −1/2 and 1/2 and using the above relations, we obtain (2.8).

Approximation of infinitely narrow reaction zone.
To study the problem analytically, we reduce it to a singular perturbation problem where the reaction zone is supposed to be infinitely narrow and the reaction term is neglected outside of it.It is a well-known approach for combustion problems [4,6,7].We fulfil a formal asymptotic analysis with = Z −1 = R 0 T 2 b /qE taken as a small parameter, and obtain a closed interface problem (i) when z ≠ ζ : The boundary conditions are the same as in Section 1.Here σ = 2σ * (see Section 2).The hat over σ is omitted below.Problem (3.1), (3.2) is coupled in the sense that it describes the thermal instability of the reaction front and the convective instability at the same time.There are different limiting cases here.For example, if the coefficient of thermal expansion β equals zero (i.e., R = 0), then this corresponds to a condensed medium since v ≡ 0. Another case is when we remove the thermal instability by decreasing the Zeldovich number Z. Linear stability analysis of problem (3.1), (3.2) without heat loss is carried out in [3].

Thermal regimes.
In this section, we consider a half-infinite continuous plug-flow reactor with a given velocity v of the medium along the axis of the reactor.We consider the 1D spatial case without hydrodynamics.It is a particular case of the complete problem (3.1), (3.2),where the width of the reactor is sufficiently small.The temperature distribution outside of the reaction zone satisfies the equation The jump conditions at the reaction zone z = ζ are The boundary condition is Moreover, we assume that the temperature is bounded.This problem is studied in [5] in the case σ = 0.

Stationary solution.
In this subsection, we find a stationary solution of problem (4.1), (4.2), and (4.3).We suppose that the reaction zone is located at z = l.Outside of it, the dimensionless temperature θ satisfies the following equation: We look for the solution of (4.4) in the form where From the boundary conditions and jump conditions at the interface and from the boundedness of the solution (c 3 = 0), we obtain where Solving this system, we find the temperature distribution .9) and the system of two equations with respect to the unknown temperature in the reaction zone θ f and the distance to the reaction zone l: In the particular case where σ = 0, µ 1 = v, and µ 2 = 0, If l = ∞, then θ f = 0 and we obtain the normal velocity of the front propagation v 2 n = I(0).If δ = 0, the integral can be found explicitly, v 2 n = 2, and we obtain the same expression for the distance to the reaction zone as in [5]:  Therefore, where (4.17) Linearizing the jump conditions at the reaction zone, taking into account the boundary condition at z = 0 and the fact that the perturbation is bounded at infinity, and introducing the notations we obtain the system of equations with respect to the unknown constants , c1 , c2 , and c4 .Here The condition of nontrivial solvability of this system gives the dispersion relation In the particular cases where there is no heat loss (σ = 0) and where l → ∞, this dispersion relation coincides with the relations obtained in [1,5].
At the oscillatory instability boundary, Ω = iφ, where φ is an unknown frequency.We can represent (4.21) as two real equations and find l from one of them and Z as a function of other parameters from the other one.The critical value of the Zeldovich number is shown in Figure 4.1.

Numerical simulations.
In this section, we use direct numerical simulations to study problem (1.3), (1.6) in two space dimensions.We consider the stream functionvorticity formulation of the Navier-Stokes equations and employ an alternative direction method to solve the finite-difference equations.
We begin with the problem without convection and analyze the convergence of solutions of the two-dimensional (2D) problem to the corresponding solution of the 1D problem as the width of the domain decreases.After that, we compute the complete problem with convection to study its influence on the thermal regimes.

Thermal regimes.
In this subsection, we consider the problem without convection (5.1) The corresponding 1D problem is where s = 2σ /L x .In Section 2, it is shown that the solution of problem (5.1) converges to the solution of problem (5.2) as σ → 0 and L x → 0 at any fixed time interval.The same method allows us to show the convergence of stationary solutions of these two problems.
We analyze the convergence of the stationary solutions numerically.We recall that there can exist high-temperature and low-temperature stationary regimes in plug-flow reactors [2].If the speed u of the medium is greater than the normal velocity of the front propagation, the high-temperature regime does not exist.
We vary two parameters σ and L x , keeping all other parameters constant.Let l 1 be such that α(l 1 ) = 0.5 for the stationary solution of problem (5.2).Denote by l(x) the function such that α(l(x), x) = 0.5 for the stationary solution of problem (5.1), (5.3) Figure 5.1 shows l 1 as a function of s and l 2 as a function of σ for L x fixed and as a function of L x for σ fixed.The values l 1 and l 2 are compared for the same values of s.We see that the results for two 2D simulations coincide.It is consistent with the fact that the 1D problem (5.2) depends on the ratio of these parameters only.So we can  expect a similar behavior for the 2D problem (5.1).Its solution converges to the solution of the 1D problem as s decreases.However, the difference between them increases as s approaches the extinction limit.We note that l d is very small compared with l max and l min (of the order 10 −3 ).So the 2D solution is practically independent of the x-variable.However, for s sufficiently large, it differs from the 1D solution.
According to the results of Section 2, l 2 should converge to l 1 as L x → 0 and σ → 0 for s fixed.Table 5.1 shows l 2 for s = 0.02 and for various L x and σ .We see that it has exactly the same values.This shows a very good convergence as L x and σ decrease, but to the value different from that of the 1D problem!The explanation of this paradox is connected with the numerical accuracy.The computations for L x = 0.5, σ = 0.005 and for L x = 0.1, σ = 0.001 are done with the same number m x of discretization points (Table 5.1, Section 1).The space step in the second case is 5 times less, the results coincide, and we conclude that the numerical accuracy is sufficient.However, this conclusion is wrong.Table 5.1 shows also the results for L x = 0.5, σ = 0.005, and for different numbers of discretization points in the xdirection (Sections 2 and 3).Increasing m x , we observe that l 2 approaches l 1 = 4.50.The conclusion of this analysis is rather unexpected: decreasing the width of the domain, we should increase the number of discretization points.If it is fixed, the numerical solution for the 2D problem does not converge to the numerical solution of the 1D problem.
If Z exceeds a critical value Z c , then the stationary solution loses its stability and periodic in time regimes appears as a result of a Hopf bifurcation.Figure 5.2 shows the amplitude of oscillations as a function of Z for the problem without heat loss (lower curve), for problem (5.1) with heat loss in the boundary conditions, and for problem (5.2) with heat loss in the equation (upper curve).We see that in agreement with the linear stability analysis, heat losses destabilize the front (Section 4.2) and that the 1D problem (5.2) provides a good approximation of the 2D problem (5.1).

Convection.
For the problem without convection (R = 0) and without heat loss (σ = 0), there can exist a stationary solution independent of x.This means that the reaction zone is horizontal and the temperature below it is greater than the temperature above it.Therefore, in the presence of gravity (R > 0), natural convection can appear.The critical Rayleigh number R c and the amplitude of convection depend on parameters.for a supercritical bifurcation.The maximum of the stream function as a function of L is well described by the square-root dependence x , (5.4) where the critical value L c x depends on the Zeldovich number and a = 1.7 appears to be the same for all Z. Increasing Z, we increase the maximal temperature gradient in the stationary temperature distribution, which is roughly proportional to qZ.This is why the convection becomes stronger.For a fixed width of the domain, the appearance of convection is determined by the value of the Rayleigh number.If R exceeds a critical value, convection appears and its amplitude grows together with R. As above, the critical value of R decreases with the increase of the Zeldovich number.If Z passes through the critical value where thermal oscillations appear, the stationary convective regime loses its stability and oscillating convective solutions are observed.

Call for Papers
This subject has been extensively studied in the past years for one-, two-, and three-dimensional space.Additionally, such dynamical systems can exhibit a very important and still unexplained phenomenon, called as the Fermi acceleration phenomenon.Basically, the phenomenon of Fermi acceleration (FA) is a process in which a classical particle can acquire unbounded energy from collisions with a heavy moving wall.This phenomenon was originally proposed by Enrico Fermi in 1949 as a possible explanation of the origin of the large energies of the cosmic particles.His original model was then modified and considered under different approaches and using many versions.Moreover, applications of FA have been of a large broad interest in many different fields of science including plasma physics, astrophysics, atomic physics, optics, and time-dependent billiard problems and they are useful for controlling chaos in Engineering and dynamical systems exhibiting chaos (both conservative and dissipative chaos).
We intend to publish in this special issue papers reporting research on time-dependent billiards.The topic includes both conservative and dissipative dynamics.Papers discussing dynamical properties, statistical and mathematical results, stability investigation of the phase space structure, the phenomenon of Fermi acceleration, conditions for having suppression of Fermi acceleration, and computational and numerical methods for exploring these structures and applications are welcome.
To be acceptable for publication in the special issue of Mathematical Problems in Engineering, papers must make significant, original, and correct contributions to one or more of the topics above mentioned.Mathematical papers regarding the topics above are also welcome.
Authors should follow the Mathematical Problems in Engineering manuscript format described at http://www .hindawi.com/journals/mpe/.Prospective authors should submit an electronic copy of their complete manuscript through the journal Manuscript Tracking System at http:// mts.hindawi.com/according to the following timetable:

5 Figure 4 . 1 .
Figure 4.1.Thermal instability boundary: critical value of Z as a function of the flow velocity.

Figure 5 . 1 .
Figure 5.1.Comparison of 1D and 2D simulations.The lower crosses represent 1D model; the upper crosses and the dashed line represent 2D model.

Figure 5 . 2 .
Figure 5.2.Amplitude of periodic oscillations as a function of Z. Lower crosses: without heat loss; upper crosses: with heat loss in the boundary condition and in the equation; and dashed curves: approximation by the squareroot formula.

Figure 5 .Figure 5 . 3 .
Figure 5.3.Maximum of the stream function as a function of the width of the domain for different Z. Dots: simulations, curves: approximation by the square-root formula.Lower curve: Z = 7.8, intermediate curve: Z = 7.9, and upper curve: Z = 8.0.

Table 5 .
1. Comparison of numerical results with different numbers of discretization points.