Heat Explosion In Porous Media Using Radial Basis Functions

The paper is devoted to the numerical investigation of the interaction between natural convection and heat explosion in a fluid-saturated porous media in a rectangular domain. The model consists of Darcy equations for an incompressible fluid in a porous medium coupled with the nonlinear heat equation. Numerical simulations are performed using the radial basis functions method (RBFs). We study the bifurcation of the periodic oscillation of the response born by Hopf bifurcation. First, a symmetry-breaking bifurcations observed; then is followed by successive period-doubling bifurcations leading to chaos.


Introduction
The study of the interaction between natural convection and heat explosion is important for the understanding and predicting thermal explosion.The theory of thermal explosion was started with the works by Semenov and Frank-Kamenetskii [1,2].In Semenov's theory, the temperature is considered to be uniform in the entire domain, while in Frank-Kamenetskii's model spatial temperature distribution is taken into account.The influence of natural convection on thermal explosion was first tackled in [3,4].It was shown that vigorous convection occurs when the value of the Frank-Kamenetskii parameter or of the Rayleigh number exceed some critical values.Dynamics of thermal explosion with natural convection described by the Navier-Stokes equations was studied in [5] for a square domain.The same issue was studied in [7] but for a rectangular domain.The authors of these two works showed how complex regimes appeared through successive bifurcations.In this paper, we continue the investigation of thermal explosion in the case of porous media.To this end, we consider nonlinear heat equation coupled with Darcy equations for fluid motion in porous media: where T denotes the temperature, v is the velocity, p is the pressure, κ the coefficient of thermal diffusivity, µ the kinematic viscosity, ρ the density, q the heat release, g is a e-mail: allali@hotmail.comb e-mail: joundy.youssef@yahoo.frc e-mail: taik ahmed@yahoo.frd e-mail: volpert@math.univ-lyon1.frthe acceleration due to gravity, γ a unit vector in the vertical direction, T 0 the characteristic value of the temperature and K p the permeability.Let us note that we consider unstationary form of the Darcy equation since conditions of heat explosion can be sensitive to the quasi-stationary approximation where the time derivative is cancelled.
In addition, the temperature dependence of the reaction rate K(T ) is given by the Arrhenius law: where E is the activation energy, R the universal gas constant and k 0 the pre-exponential factor.This system is considered in a 2D rectangular domain, 0 ≤ x ≤ n , 0 ≤ y ≤ .The boundary conditions will be specified below.
The paper is organized as follows.We first give the model in the stream function-vorticity formulation in the next section.The numerical method is described in Section 3. Section 4 deals with numerical simulation and Section 5 concludes the work.

The stream function-vorticity formulation
In order to write the dimensionless model, we introduce new spatial variables and keeping for convenience the same notation for the other variables, we obtain: We introduce the stream function ψ using incompressibility of the fluid, and the vorticity We apply the rotational operator to equations ( 6)-( 7) in order to eliminate the pressure.We get the equations with the following boundary condition: and some appropriate initial conditions on θ, ω , ψ. Here Eκµ is the Rayleigh number.Under the assumptions of large activation energy, RT 0 /E << 1, we can perform the Frank-Kamenetskii transform, so that the nonlinear reaction rate in the equation (11) is taken to be F K exp(θ) [2].Next, σ = 1 V a stands for the inverse of Vadasz number, Prandtl number and D a = K p 2 is the Darcy number.

Numerical method
Let us assume that (θ m , ω m , ψ m ) are known at time t = m∆t.They can be approximated as follows: ωm ψm We choose the Hardy's multiquadric or the MQ-RBF a given set of centres where X c i = (x c i , y c i ).For each collocation points X = (x, y) ∈ R 2 , the MQ-RBF is defined by ε is called a shape parameter.
Using the implicit scheme for the time derivative, the equations ( 11)-(13) become: where ⊗ is the element-by-element matrix product.

Periodic oscillations, symmetry-breaking and period-doubling
In our simulations, we fix the value of the Rayleigh number, R p = 10000, the geometry of the domain, n = 8, and σ = 0.01.We will vary the Frank-Kamenetskii parameter F K .Fig. 1 shows the convergence of the solution to a stationary state for the value of F K = 1.87.A small increase of this parameter leads to a Hopf bifurcation, after which the system oscillates.The corresponding limit cycle in the (ψ, θ) plane is an ellipse (Fig. 2(a)).When F K increased slightly, the system loses its symmetry (Fig. 2 Fig. 3(a) and 3(b) illustrate the first period-doubling where the system needs a 2T to browse the entire cycle.Further increase of this parameter to F K = 3.78 and to F K = 3.812 results in sequential bifurcations to the solutions with the periods 4T and 8T shown in Fig. 3(c) and 3(d), respectively.

Thermal explosion chaos
After successive period-doubling bifurcations, chaotic behavior of solutions is observed (Fig. 4(a)).Continuous frequency spectrum (Fig. 4(b)) is specific for chaotic oscillations.

Conclusion
Thermal explosion in porous media is modelled by a system coupling a non-linear heat equation with a system describing the fluid motion in rectangular domain with unsteady terms in the Darcy quation.The model is studied numerically using multi-quadratic radial basis functions.The Frank-Kamenetskii number F K was chosen as a bifurcation parameter.Increasing F K , we found transitions between periodic and chaotic oscillations through a sequence of period-doubling bifurcations.