A MESHFREE METHOD FOR HEAT EXPLOSION PROBLEMS WITH NATURAL CONVECTION IN INCLINED POROUS MEDIA

This paper investigates the interaction between natural convection and heat explosion in porous media. A meshless collocation method based on multiquadric radial basis functions has been applied to study the problem in an inclined two-dimensional porous media. The governing equations consist of coupling the Darcy equations in the Boussinesq approximation of low density variations to the heat equation with a nonlinear chemical source term. The numerical results obtained are in good agreement with some previous studies that consider the vertical direction. A complex behaviour of solutions is observed, including periodic and aperiodic oscillations. We have shown that a small inclination of the container stabilizes the reactive fluid and can prevent thermal explosion.


Introduction
Consider an enclosure filled with fluid subjected to a temperature difference on two opposite walls, while the rest of the walls are adiabatic. The fluid beside the hot wall will be heated and due to its decrease of density buoyancy will carry it upwards. An inverse phenomenon will occur along the cold wall causing the fluid to be colder and denser and thus it will travel downwards. This natural convection phenomenon appears frequently in nature and is present in many industrial applications. One key issue concerning the heat transfer process in fluids, is how efficient it is, if it varies depending on some critical parameters. In heat explosion domain, the development of natural convection can significantly alter the heat transfer from the reacting fluid to the environment and hence alter the balance between heat production and heat loss, which determines whether or not the system will explode. This work may be considered as a continuation of a series of investigations concerning the interaction between heat explosion and natural convection [1,6,7]. Most researches in this area, like [2,4], are concerned with the case when the fluid enclosure is oriented vertically or horizontally. In the present paper, we consider the case when the heat transfer process is produced in an inclined rectangular enclosure that is completely filled with a fluid-saturated porous medium. This is a complex process which involves various critical parameters like the Rayleigh number, the enclosure inclination angle and the Frank-Kamenetskii parameter. A range of numerical simulations is conducted by using the multiquadric radial basis functions (MQ-RBF) collocation method. a e-mail: Loubna.SALHI@um6p.ma The paper is organized as follows: The governing equations are presented in section 2. The stream functionvorticity formulation is described in section 3. A brief description of the numerical method used is given in section 4, and numerical results are presented in section 5.

Governing equations
The geometry considered is an inclined rectangular enclosure with an inclination angle α as shown in Figure  1. We consider the nonlinear heat equation coupled with Darcy equations for fluid motion in porous media. The fluid is Newtonian and all the thermophysical properties are supposed to be constant, except for the density in the buoyancy term that can be adequately modelled by the Boussinesq approximation, and that compression effects and viscous dissipation are neglected.
Under these assumptions, the problem is described in a two-dimensional space, 0 ≤ x ≤ nl, 0 ≤ y ≤ l, by the following equations: The boundary conditions will be specified below.
Here, T is the temperature, v = (u, v) the velocity, p the pressure, κ the coefficient of thermal diffusivity, q the adiabatic heat release, ρ the density, µ the coefficient of kinematic viscosity, g the gravity acceleration, β the coefficient of thermal expansion, T 0 is the characteristic value of the temperature and K p is the permeability. The function K(T ) describes the reaction rate where the temperature dependence is given by the Arrhenius law where E is the activation energy supposed to be large in this problem, R 0 is the universal gas constant and k 0 is the pre-exponential factor.

Stream function vorticity formulation
In order to obtain the dimensionless form, we introduce the following spatial variables: For convenience, we keep the same notation for all variables except for the temperature, so we may re-write system (1)-(3) as follows: With the boundary conditions: For incompressible flows with constant fluid properties, Eqs. (4)-(7) can be simplified by introducing the stream function Ψ and vorticity ω as dependent variables. We We apply the rotational operator to Eqs. (5)-(6) in order to eliminate the pressure. We get With the boundary conditions : and some appropriate initial conditions on θ, ω and Ψ . Here 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 Eq. (8) is taken to be F k exp(θ) as shown in [3]. Next, σ = 1 V a stands for the inverse of Vadasz number, V a = P r D a , P r = µ κ is the Prandtl number and D a = K p l 2 is the Darcy number.

MQ-RBF method for numerical interpolation
To understand the effect of natural convection phenomena on heat transfer in the case of quasi-steady flow, we use the multiquadric radial basis function (MQ-RBF) method because of its popularity in many applications and its good approximation properties [5,8]. The multiquadric is representive of the class of RBFs that are global, infinitely diffetentiable, and contains a shape parameter ε 0 which plays an important role for the accuracy of the method. The MQ-RBF is radially symmetric about its center, x ∈ R 2 , and its argument r = x 2 is dependent on the node location. It is defined by Now, let us assume that (θ n , ω n , Ψ n ) are known at time t = n∆t. They can be approximated as follows: where x c 1 ,x c 2 ,...,x c N ∈ Ω ⊂ R 2 is a given set of centres with x c i = (x c i , y c i ) and ϕ the MQ-RBF defined in (11).
By using the implicit scheme for the time derivative, Eqs. (8)-(10) become: where ⊗ is the element-by-element matrix product.

Numerical simulations
In all our simulations, we fix the value of the Rayleigh number, R p = 10000, the geometry of the domain, n = 8, σ = 0.01, F K = 3.8, and we vary the angle of inclination α.
We start with the case where the reservoir is placed on an horizontal foundation, α = 0 (no-inclination). Figure  2(a) shows the aperiodic behavior of the solution. Indeed, the system is aperiodic represented by a thick trajectory in the phase space of the mean of the temperature as a function of the mean of the stream function (see Figure  2(b)). A slight increase in the angle of inclination leads to a periodic regime. In Figure 3(a), we observe a periodic solution corresponding to the value of the inclination angle α = 0.01π as well as its phase space in Figure 3(b). If we continue to increase the angle α, the oscillations disappear and the solution converges to a steady state. Figure 4 shows the stationary solution for α = 0.011π. In order to highlight the way in which the aperiodic solution appears or disappears by varying the angle of inclination of the container we have drawn the bifurcation diagram (see Figure 5). Two types of behavior can be identifed.
In the first one, the attractor is represented by a set of points distributed over the vertical line (α ∈ [0; 0.0021π]), in the second one, the attractor is represented by a finite number of points attained successively, which corresponds to a periodic or steady behavior. We remark that aperiod motion disappears when the angle α exceeds the value of 0.0021π. We observe a reverse period doubling cascade near α = 0.007π. Figures 6(a) and 6(b) illustrate the solution obtained for the period 4 taking into consideration the inclination angle α = 0.00205π.

Conclusion
In this work, the problem of heat explosion with natural convection is studied in an inclined porous media. The inclination angle is considered to be the control parameter. The MQ-RBF collocation method has been used for numerical simulations. In the case of an inclined medium, the instability decreases considerably by slightly increasing the angle of inclination. We found transitions between aperiodic and periodic oscillations through a sequence of period-doubling bifurcations.