The non-linear energy stability of Brinkman thermosolutal convection with reaction

We use the energy method to obtain the non-linear stability threshold for thermosolutal convection porous media of Brinkman type with reaction. The obtained non-linear boundaries for different values of the reaction terms are compared with the relevant linear instability boundaries obtained by Wang and Tan (Phys Lett A 373:776–780, 2009). Using the energy theory we obtain the non-linear stability threshold below which the solution is globally stable. The compound matrix numerical technique is implemented to solve the associated system of equations with the corresponding boundary conditions. Two systems are investigated, the heated below salted above case and the heated below salted below case. The effect of the reaction terms and Brinkman term on the Rayleigh number is discussed and presented graphically.


Introduction
Convection in porous media has attracted the attention of many researchers and has been an area of great interest in addition to its wide range of applications. Thermal convection in porous media and stability analysis returns back to Horton and Rogers [2], Lapwood [3] and Nield and Barletta [4]. The problem of double-diffusive convection in porous media is well investigated by Nield [5], Rudraiah et al. [6], Wollkind and Frisch [7,8], Nield and Bejan [9], Ingham and Pop [10,11], Vafai [12,13] and Vadasz [14]. Bdzil and Frisch [15] performed a linear stability analysis where the fluid catalysed at the lower boundary of the layer and they developed their work in Bdzil and Frisch [16] and a similar work carried by Gutkowicz-Krusin and Ross [17]. Many recent studies in double and multi-component convection are accomplished by Rionero [18][19][20][21]. The first study on the reactive convection in porous media was due to Steinberg and Brand [22,23]. More studies were carried out by Gatica et al. [24,25], Viljoen et al. [26] and Malashetty and Gaikwad [27]. Pritchard and Richardson [28] figured out a model similar to that of Steinberg and Brand [22,23]. They considered the Darcy model to study the onset of thermosolutal convection using linear instability technique. Wang and Tan [1] extended the previous work of Pritchard and Richardson [28] in which Wang and Tan [1] considered Darcy-Brinkman model and used normal mode analysis to carry out a linear instability analysis.
We are studying nonlinear stability using an energy stability method. This method is being used extensively by many leading mathematicians, see for example, Straughan [29,30], Rionero [31], Capone et al. [32], Straughan [33], Capone and De Luca [34], Rionero and Torcicollo [35], Capone et al. [36], Lombardo and Mulone [37], Rionero [38], De Luca [39] and De Luca and Rionero [40]. The work in this paper may be considered as an extension of Wang and Tan [1] and Pritchard and Richardson [28]. Al-Sulaimi [41] used the energy method to carry out a nonlinear stability analysis of Darcy thermosolutal convection with reaction. In this article, the energy stability of Brinkman thermosolutal convection with reaction is considered. The compound matrix numerical technique is used to solve the associated system of equations with the corresponding boundary conditions. Two systems are investigated separately, the heated below-salted above system and the heated below-salted below system. The energy stability boundaries obtained for different values of the reaction rates are compared with the relevant linear instability boundaries. Some linear instability boundaries are obtained by Wang and Tan [1], but they do not correspond directly to what we require and hence we recompute also the linear values using the D 2 Chebyshev tau method.
The aim of the study is to obtain the nonlinear stability boundaries below which the solution is globally stable by using the energy method and compare the nonlinear boundaries with the relevant linear instability boundaries obtained by Wang and Tan [1]. Considering a porous medium of Brinkman type occupying a bounded threedimensional domain, the variation of the onset of thermosolutal convection with the reaction rate and Brinkman coefficient is discussed.
Here v i , p, T, C are the velocity, pressure, temperature and salt concentration. K is the matrix permeability, μ is the fluid viscosity, ρ 0 is the fluid density. k C , k T are the molecular diffusivity of the solute through the fluid and the effective diffusivity of the heat through the saturated medium. M is the ratio of the heat capacity of the fluid to the heat capacity of the medium,φ is the matrix porosity,k is the reaction coefficient and f 0 + f 1 (T − T 0 ) = C eq (T ) in Pritchard and Richardson [28], where f 0 , f 1 and T 0 are constants. Moreover, g is the gravity, k = (0, 0, 1) and α T and α C are the thermal and solutal expansion coefficients respectively. The symbol is the Laplace operator. The Eq. (1) are taken in the domain where T L , T U , C L , C U all constants, with T L > T U since our systems are heated below. For the salted above porous medium C U > C L while for the salted below case C L > C U . In the steady state, we look for Assuming C eq (T (z)) =C(z) (see Pritchard and Richardson [28] and Al-Sulaimi [41]), we find the steady solution or the basic state to (1) which we are interested in studying its stability and which satisfies (2) as where β T = (T L − T U )/d and β C = (C L − C U )/d are the temperature and salt gradients respectively. To analyze the stability of the solutions (4) we define perturbations (u i , π, θ, φ) such that Using these perturbations in Eq. (1) we derive the equations governing (u i , π, θ, φ) as where w = u 3 . To non-dimensionalize the system (6), we define the length, time and velocity scales, L, τ and U , by L = d, τ = d/MU and U = k T /d. We introduce pressure, temperature and salt scales as where Le = k T /k C is the Lewis number. The temperature and salt Rayleigh numbers are defined as Then, the fully nonlinear, perturbed dimensionless form of (6) is where ε = M Le,γ = λK /μd 2 the Brinkman coefficient and h and η are the reaction terms Moreover, +R s is taken for the salted below system and −R s is taken for the salted above system. The corresponding boundary conditions are

Linear instability theory
To study the linear instability, we drop the nonlinear terms of (7) and take the double curl of equation (7) 1 and retaining only the third component of the resulting equation to reduce (7) to studying the system where * is the horizontal Laplacian. Assuming a normal mode representation for w, θ and φ of the form (see Straughan [29]) and a is a wave number. Using (10) and applying the normal mode representations to (9), we find where D = d/dz. This is an eigenvalue problem for σ to be solved subject to the boundary conditions System (11) with the corresponding boundary conditions (12) is solved using the D 2 Clebyshev tau method. Detailed numerical results for the heated below-salted above and heated below-salted below are reported separately in the subsections (6.1) and (6.2). We determine the critical Rayleigh number given by

Nonlinear energy stability theory
In order to study the nonlinear stability of the Brinkman model for the double diffusive convection, we consider the nonlinear system of equations in the dimensionless form (7) and the corresponding boundary conditions (8). Taking into consideration the periodicity of the system and the smoothness of the boundary to allow the application of the Divergence Theorem. Multiply Eq. (7) 1 by u i and integrate over V using integration by parts. Similarly, multiply Eq. (7) 3 by θ and Eq. (7) 4 by φ and integrate. The following system of energy equations is obtained Then we form the combination of the equations in system (13) as where λ a coupling parameter. This leads to the energy identity where Then is an energy inequality which follows from the energy identity, where H is the space of admissible solutions. Namely The nonlinear stability ensues when R E > 1 which implies that 1 − 1/R E > 0. By using the Poincaré inequality we can show where k = min 1 M Le , 1 . Then from (16) we may derive the inequality where the coefficient a 1 is defined by Upon integration we obtain Inequality (19) shows that under the condition R E > 1, E(t) → 0 as t → ∞. This result according to Eq. (15) 1 , proves that θ 2 → 0 and φ 2 → 0 as t → ∞.
To show the decay of u , we have to use the Poincaré inequality, the Arithmetic-Geometric Mean inequality and the fact w 2 ≤ u 2 in the energy equation (13) where α and β are constants to be chosen such that Rα + R s β = 1, which gives α = 1/2R and β = 1/2R s . This leads to Inequality (21) shows that R −1 E guarantees in addition to the decay of θ and φ , also decay of u .
Turning our attention to the maximization problem (17). We have to solve it by deriving the Euler-Lagrange equations. The maximum problem is Rescaling φ by puttingφ = √ λφ. Equation (22) will be where The Euler-Lagrange equations for this maximum are where P is a Lagrange multiplier. To remove the Lagrange multiplier, we take the double Curl of equation (24) 1 and retaining only the third component of the resulting equation to reduce (24) to where * = ∂ 2 /∂ x 2 + ∂ 2 /∂ y 2 is the horizontal Laplacian. Introducing the normal mode representation as presented in Sect. 3, system (25) becomes The Laplace operator is equivalent to = D 2 − a 2 , where D = ∂/∂z. The corresponding boundary conditions are We can determine the critical Rayleigh number given by Ra 2 E = max λ min a 2 R 2 (a 2 , λ), where for all R 2 < Ra 2 E the system is stable.

Numerical method
We have used the D 2 Chebyshev tau method (Dongarra et al. [42]) to find the bound for the linear instability theory, system (11) and the corresponding boundary conditions (12). For the energy theory we have used the compound matrix technique (Lindsay and Straughan [43]).

The compound matrix technique for the energy theory
To employ the compound matrix method (Lindsay and Straughan [43]), we have to write system (26) as The compound matrix for (33) works with the 4 × 4 minors of the 8 × 4 solution matrix formed from The solutions U i for i = 1, 2, 3, 4 are independent solutions to (33) implies that y 1 = W 1 W 2 W 3 W 4 + · · · , which gives 24 terms for y 1 . So, the idea is to derive y 2 , . . . , y 70 similarly and then obtain differential equations for the y i 's by differentiation. There is no need to write out the whole determinant each time. The first term, y 1 , suffices. By differentiating each y i and substituting from Eq. (33) we obtain differential equations for the y i , cf. Lindsay and Straughan [43] and chapter 19 of Straughan [29]. These equations are integrated numerically from 0 to 1. We keep the boundary conditions (27) at z = 0 and replace the ones at z = 1 by which using the y i 's yields the initial condition for the y i 's as Using y i 's, the final condition which satisfies (27) is seen to be The eigenvalue R is varied until (37) is satisfied to some pre-assigned tolerance.
6 Numerical results and conclusion

Heated below salted above system
The numerical integration is carried out for different values of the reaction rates, h and η and different values of the Brinkman coefficientγ . We found that when the layer is heated below and salted above in the case of no reaction i.e. h = η = 0 and when Brinkman coefficientγ = 1 that the numerical methods used give exactly the same values for Ra L and Ra E . The graphical representation of these values shows that the linear instability threshold coincide with the energy stability threshold as it is clear in Fig. 1a and that there is no region of subcritical instability. As we increase the values of the reaction rates h and η, the linear instability boundary starts to diverge from the energy stability boundary. Figure 1 shows the effect of increasing the values of the reaction rates, as we increase the values of h and η, the gap between the boundaries increases. Any point (Rs 2 , Ra 2 ) in the space above the linear instability boundary, the solid line Ra 2 L , represents a region where the system is unstable because the linear instability boundary guarantees instability. On the other hand, if (Rs 2 , Ra 2 ) lies below the energy stability boundary, the dashed line Ra 2 E , represents the space where the system is definitely stable. Note that as the reaction rates increase, the peak  of the linear instability curve moves to a higher position resulting in a wider region of possible subcritical instability between the energy stability threshold and the linear instability threshold. Moreover, there is a slight noticeable decrease in the energy stability threshold as the values of R s → +∞. Table 1 represents some numerical values obtained.
To study the effect of each one of h and η on the stability of the system, a bigger difference between their values is considered. It has been noticed that when h is bigger compared to η, the region of possible subcritical instability is wider and increasing the value of h implies more divergence of the linear instability boundary from the energy stability boundary and a movement of the peak value of the linear instability threshold to a higher position, as Fig. 2a shows. Compared to the case when η has a bigger value than h, the linear and energy boundaries coincide as shown in Fig. 2b and the linear boundary covers the content of stability. This is expected, as system (7) shows that hΘ is a destabilizing term while −ηΦ is a stabilizing term.
Examining the effect of different values of the Brinkman coefficient (effective viscosity term) on the stability boundaries, reveals that increasing the value ofγ results in a wider space of global stability below the energy stability threshold and a wider region of potential subcritical instability. The effect of different values of γ (= 0.5, 2) are presented graphically in Fig. 3.

Heated and salted Below system
It is instructive to write system (7) and the boundary conditions (8) for the salted below case as an abstract equation of form where u = (u 1 , u 2 , u 3 , θ, φ), N (u) represents the nonlinear terms in (7) so and L is the linear operator. In fact, the linear operator for (7) is We may split L into a symmetric plus skew-symmetric part as follows and For the salted above case, the previous subsection, L A would be zero and the analogous linear operator L would be symmetric. Even when h = 0 in the salted below case, we expect some problem with nonlinear energy stability theory since where (·, ·) is the inner product on (H 1 (V )) 5 with V being a period cell for the solution.
For the problem of this subsection, governed by Eqs. (7) and (8) for the salted below case, we have two sources of anti-symmetry, the R s term and the h term.
The numerical values are presented graphically for different values of the reaction rates h and η in Fig. 4. It has been noticed that as the reaction rate increases, the gap between the linear instability and energy stability boundaries increases due to the divergence of the linear threshold yielding a wider region of potential subcritical instability. Whereas, the energy stability threshold is approximately constant or more precisely it is decreasing unnoticeably as shown in Figs. 4 and 5. As expected from system (7) one sees that hΘ will destabilize the system while −ηΦ will stabilize the system which is clear and shown in Fig. 5 i.e, when the value of h is smaller compared to η the space of possible subcritical instability is less compared to the case when h is larger than η. The effect of changing the value ofγ can be noticed in Fig. 6 for γ = 0.5, 2. The gap between the boundaries increases and the space of global stability is wider asγ increases. The numerical values and their graphical representations show that the linear instability theory does not necessarily represent accurately the onset of convection and we may explain that this is due to the two sources of anti-symmetry the R s term and the h term. By this we mean that the linear instability boundary is definitely a threshold for instability, but in this case, it may be possible for instability to arise with a Rayleigh number below the linear instability boundary.