Non-linear stability bounds for a horizontal layer of a porous medium with an exothermic reaction on the lower boundary

We use the energy method to find regions of stability for a horizontal layer of a Darcy porous medium with an exothermic reaction on the lower layer. The results are compared to the linear instability results for this model found by Scott and Straughan [16]. It is shown that there is a region in which sub-critical instabilities may occur, but for small Lewis numbers, 0oLeo10, the non-linear stability boundary is reasonably close to the linear instability boundary. The effect of varying the parameters of the reaction on the stability curve is discussed. & 2013 The Author. Published by Elsevier Ltd. All rights reserved.


Introduction
For many convection problems linear studies have been used to demonstrate instability for Rayleigh numbers greater than a critical value. These studies include a number of works that consider systems containing a reaction such as Eltayeb et al. [5], Malashetty and Biradar [11] and McKay [12]. Although useful, this technique provides limited information about the system as it does not show that it is stable below this critical Rayleigh number. In order to prove stability for regions below the instability curve non-linear methods must be used; for example the variational methods used by Mulone and Rionero [14] and Hill et al. [6] and for non-constant boundary conditions by Capone and Rionero [2]. Particularly relevant is the energy method used by McTaggart and Straughan [13] to develop stability thresholds on a fluid with a reaction on the lower boundary. In the standard Bénard problem for a fluid it was shown by Joseph [7,8] that the linear instability and non-linear stability boundaries coincide. In this case the energy method is therefore of significant use as when combined with the linear method the stability of the system is fully described. It is shown in Chapter 4 of Straughan [17] that this is also the case for a Darcy porous medium and by Ahmadi [1] for a micropolar fluid layer heated from below. However, for many problems the stability boundary obtained by the energy method is below the linear instability boundary, in some instances by a large degree. Current methods are unable to show stability properties between these two boundaries and instead the equations must be solved numerically at each point using a three-dimensional computation. It is often found that sub-critical instabilities occur, for example Veronis [19] finds this is the case for a rotating fluid and Joseph and Shir [10] and Joseph and Carmi [9] demonstrate this for problems involving internal heating.
Postelnicu [15] and Scott and Straughan [16] linearly investigated a horizontal layer of a saturated porous medium with an exothermic reaction on the lower layer and discussed how the boundary reaction terms affect the instability boundary. The aim of the current work is to develop optimised stability boundaries for this problem using a fully non-linear energy method.

Non-linear perturbation equations
We begin by presenting the equations for our model as Contents lists available at ScienceDirect cf. Straughan [18], Scott and Straughan [16], where v, p, T and C are the velocity, pressure, temperature and reactant concentration, ρ 0 is the density at reference temperature T 0 , g is gravity, k ¼(0,0,1), μ is the dynamic viscosity of the fluid, α is the coefficient of thermal expansion, K is the permeability of the porous medium,φ is the porosity of the medium, M ¼ ðρ 0 c p Þ f =ðρ 0 cÞ m with ðρ 0 cÞ m ¼φ ðρ 0 c p Þ f þ ð1ÀφÞðρcÞ s and κ ¼ k m =ðρ 0 c p Þ f is the thermal diffusivity of the porous medium, k m being given by k m ¼ k s ð1ÀφÞ þ k fφ . Here the subscripts s and f represent the solid and saturating fluid, respectively.
On the upper boundary wall (z ¼h) the temperature and reactant concentration are at constant values, on the lower boundary wall (z¼ 0) an exothermic reaction occurs in which the reactant is converted to an inert product and there is no mass flux across either wall. The boundary conditions are then given by where n i is the unit normal to the boundary, h is the depth of the layer, k 0 is the rate constant, R n is the universal gas constant, E is the activation energy of the reaction, Q is the heat of the reaction, k T is the rate at which heat is conducted from the surface and k c is the reactant diffusivity.
Employing the non-dimensionalisations where Le ¼ κ=ðφk c Þ is the Lewis number. Application of the nondimensionalisations to the boundary conditions (2), (3) yields By now assuming that the concentration and temperature are dependent on only the vertical position, i.e. that C ¼ C(z) and T¼ T (z), we find that at the steady state, ðv; T ; C ; pÞ, the temperature, reactant concentration and pressure are Evaluating this steady state on the upper and lower boundaries gives the relations Given values of A, B and ξ, the relations (5) may be solved to find corresponding values for β 1 and β 3 . We introduce small perturbations u i , θ, ϕ, π from the steady state to Eqs. (4) such that and then subtract the steady state solution to obtain the nonlinear perturbation where w ¼ u 3 and the Rayleigh number is defined by Following this same process, the non-linear perturbation boundary conditions are Scott and Straughan [16] discuss how the values of A, B, ξ and Mφ influence the instability curve, however in the current work we are unable to consider the effect of ξ. It is not possible to carry out the following non-linear analysis with the conditions on the lower boundary in their current form. Instead we here consider, as a first approximation, the case E very small and consequently, we let ξ-0. In this limit the non-dimensional perturbation boundary conditions become

The energy method
It is now beneficial to rescale θ and ϕ byθ ¼ ffiffiffi where we have dropped the^symbols, R a ¼ ffiffiffi R p , and the boundary conditions (7) remain unchanged.
Next we define a general convection cell, V, in which 0 o z o 1 and the solution is periodic in x and y. We then multiply (8) 1 by u i , (8) 3 by λθ, (8) 4 byμϕ=Mφ, where λ;μ 4 0 are constants, and integrate over V, using the boundary conditions. We then add the results to find We define an energy E as in order to write (9) as where H is the set of admissible functions over which I and D are defined and we aim to find a maximum. By use of Poincaré's inequality we see that there exists a constant ξ such that D Z ξ 2 E and therefore if 1Àmax H I =D Z 0 we find We define R E by and see from (10) that if 1=R E r 1 the energy E decays exponentially and the system is stable, thus R E ¼1 on the stability boundary. Before continuing we now define the parameter μ to be μ ¼μ=Mφ 4 0: ð11Þ

Euler-Lagrange equations
In order to maximise I =D we must find the Euler-Lagrange equations using the calculus of variations, Section IV of Courant and Hilbert [3], Straughan [17], where the variables u; θ; ϕ are varied by arbitrary functions h; η 1 ; η 2 , respectively. We find that max H I =D occurs when cf. Chapter 2 of Straughan [17], where for the current problem δI and δD are defined as After differentiating and evaluating at ϵ ¼ 0 these become where the first three terms of δD have been integrated once using the boundary conditions and (since H is restricted to divergence free functions) the incompressibility condition u i;i ¼ 0 has been included in δI . We know that the stability boundary occurs when R E ¼1 and so from Eq. (12) we must solve δI ÀδD ¼ 0. As Eqs. (13) and (14) must hold for any arbitrary h; η 1 ; η 2 we find that the Euler-Lagrange equations are We also find that the natural boundary conditions that arise are Courant and Hilbert [3, pp. 208-211]. Next, we take curlcurl (15) 1 to eliminate the pressure term, then employ the Fourier transformation for a wave number k. Using similar expansions for θ and ϕ we obtain where D ¼d/dz, and the boundary conditions and It is now possible to solve Eqs. (16) with the boundary conditions (17) and (18) to find a value of R ¼ R 2 a for which the system is stable.

Numerical method and results
For each value of Le we initially fix λ and μ and use the Natural D Chebyshev-Tau method to solve for R a by treating it as an eigenvalue for the system, cf. Dongarra et al. [4]. To do this we use the boundary conditions (17) to define the natural variables χ 1 , χ 2 , χ 3 as Using these new variables, Eqs. (16) may be written as with the boundary conditions First, for each value of Le, with λ and μ remaining fixed, we minimise over the wave number k to find a value of R a for which we know that the system is stable. We then vary λ and μ to optimise the value of R a obtained. Defining λ s and μ s to be the values of λ and μ which maximise R a , k s to be the value of k for which this maximum occurs and R s to be the maximum Rayleigh number, R s ¼ R 2 a , we find the results given in Table 1 and Figs. 1-3. Here the linear instability results, R c , obtained by Scott and Straughan [16] are included for comparison. To the left of the peak in the R c curve stationary convection is dominant and to the right of this peak oscillatory convection dominates. The results for A¼0.5, B ¼0.5 are shown graphically in Fig. 1 and those for A¼0.5, B¼ 1 and A¼1, B¼ 1 are shown in Figs. 2 and 3, respectively.

Conclusions
The linear instability results obtained by Scott and Straughan [16] are useful as they provide a boundary above which we know that the system (6) with boundary conditions (7) is unstable, however they provide no information about stability or instability below this curve. The results in Table 1 and Figs. 1-3 provide stability curves, for values of R below these curves the system is stable. For the given choices of A and B we see that when the linear stationary convection occurs the two curves are close and as Table 1 For A¼ 0.5, B¼ 0.5, values of μ, λ, k for which the non-linear stability curve is optimised and the corresponding values of R s . The linear instability boundary R c is given for comparison.   oscillatory convection becomes dominant they differ more significantly. For the greater value of B the stability curve falls away from the instability curve more quickly. Scott and Straughan [16] found that when ξ ¼ 0 doubling A doubled R c , Figs. 2 and 3 show that this also doubles the non-linear stability boundary R s . For 0 o Le o 10 the stability curve remains relatively close to the instability curve, irrespective of which value of A or B we choose, and the non-linear energy analysis used here is therefore of use. From definition (11) we see that the value of Mφ does not influence the non-linear stability boundary, however it will affect the value of μ s for which this boundary occurs.
It is not possible to describe the stability properties of points that lie between the stability and instability curves the using current methods. To investigate the stability within this region the equations must be solved using a three-dimensional computation. As is usual with this type of problem we expect to find that subcritical instabilities occur.