Thermal Convection in a Rotating Anisotropic Fluid Saturated Darcy Porous Medium

Abstract

The stability of the thermal convection in a fluid-saturated rotating anisotropic porous material is investigated. We take into account the rotation of a layer of saturated porous medium about an axis orthogonal to the planes bounding the layer. The permeability is allowed to be an anisotropic tensor. In particular, we restrict our attention to the case where the permeability in the vertical direction is different to that in the horizontal plane. The linear instability and nonlinear stability analysis, in the case where the inertial term vanishes, are performed. It is shown, by using an energy method, that the nonlinear critical Rayleigh numbers coincide with those of the linear analysis. The results reveal that the system becomes more stable when the rotation is present.

1. Introduction

Thermal convection in a rotating porous medium is an active topic of research since it has many applications including geophysics, chemical engineering, food process industry, binary alloy solidification, cooling of electronics equipment, solidification and centrifugal casting of metals and rotating machinery (see e.g., Vadasz [1], Vadasz [2], Vadasz and Govender [3], Govender [4] Nield and Bejan [5], Ingham and Pop [6], and the references therein). Indeed, thermal convection involving the rotation of the layer of saturated porous medium is a subject receiving attention and is being studied extensively by many researchers, such as Vadasz [1,7,8,9]. In particular, he investigated the effect of the Coriolis force on thermal convection when the Darcy model is extended by including the time-derivative term in the momentum equation [9]. A comprehensive review of thermal convection in a rotating porous medium are given by Vadasz [10]. Palm and Tyvand [11] showed that the results of thermal instability in a rotating porous layer are equivalent to those of non-rotating anisotropic porous media. Vadasz and Govender [3] also considered the influence of gravity and centrifugal forces on the onset of convection in a rotating porous layer. Straughan [12] presented an analysis of the nonlinear stability problem for convection in a rotating isotropic porous medium. He showed that the global nonlinear stability boundary is exactly the same as the linear instability.

It is important to note that the above-mentioned studies considered assumed saturated that the porous medium is isotropic. However, the effect of anisotropy combined with the rotation effect on thermal instability has been the contribution of Alex and Patil [18], who investigated thermal instability subject to the centrifugal acceleration and the anisotropy effect as in the case of both the Darcy and Brinkman models. Govender [4] considered the Vadasz paper 1994 [7] but included the anisotropy effects for both permeability and thermal diffusivity. Later, Malashetty and Swamy [19] also performed linear instability and weakly nonlinear theory to investigate the anisotropy effects on the onset of convection in a rotating porous medium. They found that increasing an anisotropy parameter for both permeability and thermal diffusivity leads to advancing oscillatory convection. The same authors in [20] employed linear instability theory to investigate the effect of both thermal modulation and rotation on the onset of the stationary convection. Govender and Vadasz [21] also deal with the effect that thermal diffusivity and permeability anisotropy have on the thermal convection in a rotating porous medium with a thermal non-equilibrium model.

Recently, Vanishree and Siddheshwar [22] performed linear instability for an anisotropic porous medium with a temperature dependent viscosity. They also found that the onset of convection in a rotating porous medium is qualitatively similar to that in a non-rotating one. Additionally, the linear instability and nonlinear stability in an anisotropic porous medium were adopted by Kumar and Bhadauria [23] who considered viscoelastic fluid in a rotating anisotropic porous medium. Saravanan and Brindha [24] deal with the onset of centrifugal convection in the Brinkman model, and Gaikwad and Begum [25] considered the onset of double-diffusive reaction convection in an anisotropic porous medium.

In this article, we consider that the system of equations is essentially the same as that given in Vadasz [9], but we allow for the symmetric permeability tensor to be anisotropic. In particular, we consider the case where the permeability in the vertical direction is different to that in the horizontal plane. In fact, we consider the case of the inverse of the permeability tensor $M=diag\left\{1/k_x,1/k_x,1/k_z\right\}$. The goal of this article is to investigate the effect of anisotropy with rotation on the stability thresholds using linear instability and nonlinear stability methods. Here, we will ignore the inertia term in the momentum equation. More precisely, we consider cases where the values of Vadasz number tend to be large [9]. We show that the critical Rayleigh number of the linear theory is the same as the critical Rayleigh number of the nonlinear theory. We observe that energy methods are very much in vogue in the current hydrodynamical stability literature cf. Rionero [13], Capone and Rionero [14], Hill and Carr [15,16], and Hill and Malashetty [17].

2. Governing Equations

Consider a layer of porous medium heated from below and bounded by two horizontal planes $z=0$ and $z=d$, with gravity acting in the vertical direction of the z-axis. We assume that an incompressible Newtonian fluid saturates the porous layer and occupies the spatial domain $\left\{(x,y)\in {\mathbb{R}}^2\right\}\times \left\{z\in (0,d)\right\}$. Furthermore, we suppose that the layer rotates about the z-axis. The Boussineq approximation is assumed to be valid.

The governing equations incorporating fluid inertia for thermal convection in an anisotropic rotating porous media of Darcy type may be written as, cf. Malashetty, and Swamy [19],

$$\begin{array}{c}{a}_0{v}_{j,t}=-p_{,j}-\mu M_{ij}{v}_i,+g{\rho }_0\alpha T{k}_j-\frac{2}{\phi }{(\mathsf{\Omega }\times v)}_j\end{array}$$

(1)

$$\begin{array}{c}{v}_{i,i}=0,\end{array}$$

(2)

$$\begin{array}{c}{T}_{,t}+v_i{T}_{,i}=\kappa \Delta T.\end{array}$$

(3)

Here, $\mathbf{v}$, t, p, T are the velocity field, time, pressure, and temperature, respectively, and $\mu$, $\kappa$, g, $\alpha$, ${\rho }_0$, $\phi$ are dynamic viscosity, thermal diffusivity, gravity, thermal expansion coefficient of the fluid, constant density coefficient, and porosity, respectively, and $\mathsf{\Omega }$ is the angular velocity vector with $k=(0,0,1)$, $a_0=\widehat{a}/\phi$ is an inertia coefficient, $\widehat{a}$ is constant and $\phi$ is porosity. Standard indicial notation is employed throughout.

The inverse of the permeability tensor is assumed to be of the form

$$M_{ij}=diag\left\{1/k_x,1/k_x,1/k_z\right\},$$

where $k_x$, $k_z$ are constants. The boundary conditions for the problem are

$$\begin{array}{c}\mathbf{n}·\mathbf{v}=0,\mathrm{on}z=0,d,\hfill \\ T=T_L,z=0,T=T_U,z=d,\hfill \end{array}$$

(4)

where $T_L$, $T_U$ are constants with $T_L>T_U$, and $\mathbf{n}$ is the unit outward normal to the boundary, so $\mathbf{n}=(0,0,1)$ on $z=d$ and $\mathbf{n}=(0,0,-1)$ on $z=0$.

When no motion occurs and the temperature gradient is constant throughout the layer, the basic steady state solution $(\overline{v},\overline{p},\overline{T})$ whose stability is under investigation is

$$\begin{array}{c}{\overline{v}}_i\equiv 0,\overline{T}=-\beta z+T_L,\hfill \\ \overline{p}=p_0-g{\rho }_{0}z-\frac{1}{2}\alpha \beta g{\rho }_0{z}^2,\hfill \end{array}$$

(5)

with $p_0$ is the pressure at the surface $z=0$, $\beta =(T_L-T_U)/d.$

Letting $v_i={\overline{v}}_i+u_i$, $T=\overline{T}+\vartheta$, $\overline{p}=\overline{p}+\pi$, the nonlinear perturbation equations arising from Equations (1)–(3), are

$$\begin{array}{c}{a}_0{u}_{j,t}=-{\pi }_{,j}-\mu M_{ij}{u}_i+k_{j}g{\rho }_0\alpha \vartheta -\frac{2}{\phi }{(\mathsf{\Omega }\times u)}_j,\hfill \\ u_{i,i}=0,\hfill \\ {\vartheta }_{,t}+u_i{\vartheta }_{,i}=\beta w+\kappa \Delta \vartheta ,\hfill \end{array}$$

(6)

where $w=u_3$.

The perturbation equations are non-dimensionalised with the following scalings as follows:

$$x_i=d{x}_i^{*},u_i=U{u}_i^{*},t=\mathcal{T}{t}^{*},\pi =P{\pi }^{*},\vartheta =T^{\#}{\vartheta }^{*},$$

$$U=\frac{\kappa }{d},\mathcal{T}=\frac{d^2}{\kappa },P=\frac{d\mu U}{k_x},T^{\#}=U\sqrt{\frac{d^2\beta \mu }{\kappa {\rho }_{0}g\alpha k_x}},$$

$${\widehat{V}}_a=\frac{\phi Pr}{\widehat{a}Da},R=\sqrt{\frac{d^2{\rho }_{0}g\alpha \beta k_x}{\mu \kappa }},\tilde{T}=\frac{2\mathsf{\Omega }{k}_x}{\mu \phi },$$

where $Ra=R^2$ is the Rayleigh number, $Ta={\tilde{T}}^2$ is the Taylor number, ${\widehat{V}}_a$ is the Vadasz number, and $Pr$ is the Prandtl number, with $Da=k_x/d^2$ being the Darcy number.

Omitting all stars, the nonlinear non-dimensional perturbation equations are

$$\begin{array}{c}{\displaystyle \frac{1}{\widehat{V}a}{u}_{i,t}=-{\pi }_{,i}+R{k}_i\vartheta -\tilde{T}{(\mathbf{k}\times \mathbf{u})}_i-m_{ij}{u}_j,}\hfill \\ u_{i,i}=0,\hfill \\ {\vartheta }_{,t}+u_i{\vartheta }_{,i}=Rw+\Delta \vartheta .\hfill \end{array}$$

(7)

Here, $m_{ij}=diag\left\{1,1,\xi \right\}$, where $\xi =k_x/k_z$ is the anisotropy parameter.

The corresponding boundary conditions are

$$n_i{u}_i=\vartheta =0,z=0,1,$$

(8)

with $\left\{u_i,\vartheta ,\pi \right\}$ satisfying a plane tiling periodicity in $\left(x,y\right)$.

3. The Principle of Exchange of Stabilities Ignoring Inertia Term

As stated in the Vadasz paper [9], the values of Vadasz for many porous media applications in a real life are large. To this end, we let ${\widehat{V}}_a\to \infty$ in the Equation (7) be

$$\begin{array}{c}-{\pi }_{,i}+R{k}_i\vartheta -\tilde{T}{(\mathbf{k}\times \mathbf{u})}_i-m_{ij}{u}_j=0,\hfill \end{array}$$

(9a)

$$\begin{array}{c}{u}_{i,i}=0,\hfill \end{array}$$

(9b)

$$\begin{array}{c}{\vartheta }_{,t}+u_i{\vartheta }_{,i}=Rw+\Delta \vartheta .\hfill \end{array}$$

(9c)

We now take curl of Equation (9a) and curlcurl of the same equation to find

$$\begin{array}{c}R({\vartheta }_{,y}{\delta }_{i1}-{\vartheta }_{,x}{\delta }_{i2})+\tilde{T}\frac{\partial u_i}{\partial z}-{\epsilon }_{ijk}{m}_{kq}{u}_{q,j}=0,\hfill \end{array}$$

(10)

and

$$\begin{array}{c}{m}_{ir}\Delta u_r-m_{jr}{u}_{r,ji}+\tilde{T}\frac{\partial {\omega }_i}{\partial z}=R(k_i{\Delta }^{*}\vartheta -{\vartheta }_{,xz}{\delta }_{i1}-{\vartheta }_{,yz}{\delta }_{i2}),\hfill \end{array}$$

(11)

where ${\Delta }^{*}={\partial }^2/\partial x^2+{\partial }^2/\partial y^2$ is the horizontal Laplacian operator, and ${\omega }_i$ is the vorticity.

Upon taking the third component of the foregoing equations, we obtain

$$\begin{array}{c}{m}_{3r}\Delta u_r-m_{jr}{u}_{r,j3}+\tilde{T}{\omega }_{3,z}=R{\Delta }^{*}\vartheta ,\hfill \end{array}$$

(12a)

$$\begin{array}{c}\tilde{T}{w}_{,z}-{\epsilon }_{3jk}{m}_{kq}{u}_{q,j}=0,\hfill \end{array}$$

(12b)

$$\begin{array}{c}{\vartheta }_{,t}+u_i{\vartheta }_{,i}=Rw+\Delta \vartheta .\hfill \end{array}$$

(12c)

We now consider the linearised Equation of (12) by removing the nonlinear term of Equation (12c), and therefore we seek for solutions of the form

$$\mathbf{u}\left(\mathbf{x},t\right)=\mathbf{u}\left(\mathbf{x}\right)e^{\sigma t},\vartheta \left(\mathbf{x},t\right)=\vartheta \left(\mathbf{x}\right)e^{\sigma t}.$$

By substituting into Equation (12) and removal of exponential parts, we have to solve the system

$$\begin{array}{c}{m}_{3r}\Delta u_r-m_{jr}{u}_{r,j3}+\tilde{T}{\omega }_{3,z}=R{\Delta }^{*}\vartheta ,\hfill \end{array}$$

(13a)

$$\begin{array}{c}\tilde{T}{w}_{,z}-{\epsilon }_{3jk}{m}_{kq}{u}_{q,j}=0,\hfill \end{array}$$

(13b)

$$\begin{array}{c}\sigma \vartheta =Rw+\Delta \vartheta .\hfill \end{array}$$

(13c)

The corresponding boundary conditions are

$$\begin{array}{c}w=\vartheta =0,z=0,1.\hfill \end{array}$$

(14)

In order to show that $\sigma \in \mathbb{R}$, and that the principle of exchange of stabilities holds, we consider a three-dimensional periodic cell V for solution to Equation (13) and assume momentarily that $\sigma$, $u_i$, and $\vartheta$ are complex. Then, we multiply Equation (13a) by $w^{*}$ (the complex conjugate of w)and integrate over V to obtain

$$\begin{array}{c}{\int }_V\left(m_{3r}\Delta u_r-m_{jr}{u}_{r,j3}\right)w^{*}dV+{\int }_V\tilde{T}{\omega }_{3,z}{w}^{*}dV={\int }_{V}R{\Delta }^{*}\vartheta w^{*}dV,\hfill \end{array}$$

(15)

since $m_{ij}=diag\left\{1,1,\xi \right\}$, so one may rewrite the first term in Equation (15) as shown below

$$\begin{array}{ccc}\hfill m_{3r}\Delta u_r-m_{jr}{u}_{r,j3}& =& m_{33}\Delta u_3-m_{11}{u}_{1,13}-m_{22}{u}_{2,23}-m_{33}{u}_{3,33}\hfill \\ & =& \xi (u_{3,11}+u_{3,22}+u_{3,33})-u_{1,13}-u_{2,23}-\xi u_{3,33}\hfill \\ & =& \xi (u_{3,11}+u_{3,22})-{\left(u_{1,1}+u_{2,2}\right)}_{,3}.\hfill \end{array}$$

Recalling $u_{1,1}+u_{2,2}=-u_{3,3}$, we have

$$\begin{array}{c}{m}_{3r}\Delta u_r-m_{jr}{u}_{r,j3}=\xi {\Delta }^{*}w+w_{,zz}.\hfill \end{array}$$

(16)

Making use of Equation (13b)

$$\begin{array}{ccc}\hfill \tilde{T}{w}_{,z}& =& {\epsilon }_{3jk}{m}_{kq}{u}_{q,j}\hfill \\ & =& {\epsilon }_{321}{m}_{11}{u}_{1,2}+{\epsilon }_{312}{m}_{22}{u}_{2,1}\hfill \\ & =& -u_{,y}+v_{,x}.\hfill \end{array}$$

(17)

Furthermore, we now make use of vorticity equation

$$\begin{array}{c}\hfill {\omega }_i=\nabla \times u_i={\epsilon }_{ijk}{u}_{k,j}\equiv \left(w_{,y}-v_{,z},u_{,z}-w_{,x},v_{,x}-u_{,y}\right),\end{array}$$

(18)

and so

$$\begin{array}{c}{\omega }_3=v_{,x}-u_{,y}.\hfill \end{array}$$

(19)

Then, we form the combination of Equations (17) and (19) to find

$$\begin{array}{c}{\omega }_3=\tilde{T}{w}_{,z}.\hfill \end{array}$$

(20)

After differentiating Equation (20) with respect to z and expressing $w_{,zz}=\Delta w-{\Delta }^{*}w$, we employ the results and Equation (16) into Equation (15) to obtain

$$\begin{array}{c}{\int }_V\left(1+{\tilde{T}}^2\right)\Delta w{w}^{*}dV-{\int }_V\left(1-\xi +{\tilde{T}}^2\right){\Delta }^{*}w{w}^{*}dV={\int }_{V}R{\Delta }^{*}\vartheta w^{*}dV,\hfill \end{array}$$

and hence we arrive at

$$\begin{array}{c}-\left(1+{\tilde{T}}^2\right){∥\nabla w∥}^2+\left(1-\xi +{\tilde{T}}^2\right){∥{\nabla }^{*}w∥}^2=-R\left({\nabla }^{*}w,{\nabla }^{*}\vartheta \right),\hfill \end{array}$$

(21)

where ${\nabla }^{*}\equiv \left(\partial /\partial x,\partial /\partial y,0\right)$, $(., .)$ and $∥·∥$ denote the inner product and norm on the complex Hilbert space $L^2(V)$.

By applying the horizontal Laplacian operator ${\Delta }^{*}$ to Equation (13c), multiplying by ${\vartheta }^{*}$ (the complex conjugate of $\vartheta$) and again integrating, we find

$$\begin{array}{c}\sigma {∥{\nabla }^{*}\vartheta ∥}^2=R\left({\nabla }^{*}w,{\nabla }^{*}\vartheta \right)+{∥{\nabla }^{*}\nabla \vartheta ∥}^2.\hfill \end{array}$$

(22)

Next, the addition of Equations (21) and (22) yields

$$\begin{array}{c}\sigma {∥{\nabla }^{*}\vartheta ∥}^2=\left(1+{\tilde{T}}^2\right){∥\nabla w∥}^2-\left(1-\xi +{\tilde{T}}^2\right){∥{\nabla }^{*}w∥}^2+{∥{\nabla }^{*}\nabla \vartheta ∥}^2.\hfill \end{array}$$

(23)

Since $\sigma ={\sigma }_r+i{\sigma }_i$, the equating the imaginary parts of Equation (23) yields

$${\sigma }_i{∥{\nabla }^{*}\vartheta ∥}^2=0.$$

Thus, ${\sigma }_i=0$ and so $\sigma \in \mathbb{R}$, which implies that the linearized Equation (13) satisfy the strong principle of exchange of stabilities. As such, the instability set in as stationary convection.

4. Linear Instability Analysis

In this section, we seek to find the critical Rayleigh number of linear theory and we follow the work of Chandrasekhar [26]. To this end, we set $\sigma =0$ into Equation (13). We further employ Equation (16), and the governing system can be reduced to

$$\begin{array}{c}\xi {\Delta }^{*}w+w_{,zz}+\tilde{T}{\omega }_{3,z}=R{\Delta }^{*}\vartheta ,\hfill \end{array}$$

(24a)

$$\begin{array}{c}{\omega }_{3,z}-\tilde{T}{w}_{,zz}=0,\hfill \end{array}$$

(24b)

$$\begin{array}{c}Rw+\Delta \vartheta =0,\hfill \end{array}$$

(24c)

where Equation (20) has been differentiated with respect to z.

We now eliminate ${\omega }_{3,z}$ from Equation (24a,b), and therefore system (24) can be written as follows:

$$\begin{array}{c}\xi {\Delta }^{*}w+\left(1+{\tilde{T}}^2\right)w_{,zz}=R{\Delta }^{*}\vartheta ,\hfill \\ Rw+\Delta \vartheta =0.\hfill \end{array}$$

(25)

To proceed, we assume a normal mode representation for w, and $\vartheta$ of the form

$$\vartheta =\mathsf{\Theta }(z)f(x,y),w=W(z)f(x,y),$$

where $f(x,y)$ is the horizontal planform that satisfies ${\Delta }^{*}f=-a^{2}f$, a being a wave number. With $D=d/dz$, we arrive at the following system

$$\begin{array}{c}\left[(1+{\tilde{T}}^2)D^2-\xi a^2\right]W=-a^{2}R\mathsf{\Theta },\hfill \end{array}$$

(26a)

$$\begin{array}{c}\left(D^2-a^2\right)\mathsf{\Theta }=-RW.\hfill \end{array}$$

(26b)

The corresponding boundary conditions are

$$\begin{array}{c}W=\mathsf{\Theta }=0,z=0,1.\hfill \end{array}$$

(27)

The variable $\mathsf{\Theta }$ is eliminated from Equation (26) to yield the fourth order differential equation

$$\begin{array}{c}\left[(1+{\tilde{T}}^2)(D^2-a^2)D^2-\xi a^2(D^2-a^2)\right]W=a^2{R}^{2}W.\hfill \end{array}$$

(28)

In view of the boundary conditions (27) and from Equation (26a), we obtain

$$D^{2}W=0,z=0,1.$$

Applying these boundary conditions to Equation (28), it turns out that

$$\begin{array}{c}{D}^{4}W=0,z=0,1.\hfill \end{array}$$

(29)

Further differentiation of Equation (28) yields

$$D^{(2n)}W=0,\mathrm{on}z=0,1,\mathrm{for}n=0,1,2,\dots$$

Thus, we may select $W= \sin n\pi z,$ for $n\in \mathbb{N}$. Upon substituting in Equation (28), we have

$$\left[(1+{\tilde{T}}^2)(n^2{\pi }^2+a^2)n^2{\pi }^2+\xi a^2(n^2{\pi }^2+a^2)\right]=a^2{R}^2,$$

which leads to

$$\begin{array}{c}{R}_L^2=\frac{\left(1+{\tilde{T}}^2\right){\pi }^2{n}^2{\mathsf{\Lambda }}_n}{a^2}+\xi {\mathsf{\Lambda }}_n,\hfill \end{array}$$

(30)

where ${\mathsf{\Lambda }}_n=n^2{\pi }^2+a^2$. Minimizing over n yields $n=1$. Then, differentiating $R^2$ with respect to $a^2$ yields the stationary convection boundary

$$\begin{array}{c}\hfill R_{L(sc)}^2={\pi }^2{\left(\sqrt{\xi }+\sqrt{1+{\tilde{T}}^2}\right)}^2,\end{array}$$

(31)

and the corresponding critical wave number $a_{L(c)}$ is given by

$$\begin{array}{c}\hfill a_{L(c)}^2={\pi }^2\sqrt{\frac{1+{\tilde{T}}^2}{\xi }}.\end{array}$$

(32)

It is worth observing that as ${\tilde{T}}^2=0$, and $\xi =1$, we recover the result for the isotropic problem [27]

$$a_{L(c)}^2={\pi }^2,R_{L(sc)}^2=4{\pi }^2.$$

5. Nonlinear Stability Analysis

In this section, we commence with the derivation of further boundary conditions that will be used to continue with the nonlinear stability analysis. To obtain these, we observe from Equations (10) and (18),

$$\begin{array}{c}{\omega }_1=\left(1-\xi \right)w_{,y}+\tilde{T}{u}_{,z}+R{\vartheta }_{,y},{\omega }_2=\left(\xi -1\right)w_{,x}+\tilde{T}{v}_{,z}-R{\vartheta }_{,x}.\hfill \end{array}$$

(33)

One may then deduce from the boundary conditions (14),

$$\begin{array}{c}{\omega }_1=\tilde{T}{u}_{,z},{\omega }_2=\tilde{T}{v}_{,z},z=0,1.\hfill \end{array}$$

(34)

In addition, from Equation (18), we also find on the boundaries

$$\begin{array}{c}{\omega }_1=-v_{,z},{\omega }_2=u_{,z},z=0,1.\hfill \end{array}$$

(35)

One may then deduce from Equations (34) and (35),

$$\begin{array}{c}{u}_{,z}=v_{,z}=0,z=0,1,\hfill \end{array}$$

(36)

and hence

$$\begin{array}{c}{\omega }_1={\omega }_2=0,z=0,1.\hfill \end{array}$$

(37)

Furthermore, we find from Equations (19) and (20) that

$$\begin{array}{c}\tilde{T}{w}_{,zz}=v_{,xz}-u_{,yz}.\hfill \end{array}$$

(38)

It follows from Equation (36) that

$$\begin{array}{c}{w}_{,zz}=0,z=0,1.\hfill \end{array}$$

(39)

Since $w\equiv \vartheta \equiv 0$ on $z=0,1$, we obtain from Equation (9c)

$$\begin{array}{c}{\vartheta }_{,zz}=0,z=0,1.\hfill \end{array}$$

(40)

Then differentiating Equation (9c) $2n$ times with respect to z, we find

$${\vartheta }_{,t}^{(2n)}+\sum _{s=0}^{2n}\left(\genfrac{}{}{0pt}{}{2n}{s}\right)u_i^{(s)}{\vartheta }_{,i}^{(2n-s)}=R{w}^{(2n)}+\Delta {\vartheta }^{(2n)},$$

where we have used the General Leibniz Rule.

Furthermore, we may rewrite the foregoing equation as shown below:

$$\begin{array}{c}\hfill {\vartheta }_{,t}^{(2n)}+\sum _{s=0}^{2n}\left(\genfrac{}{}{0pt}{}{2n}{s}\right)\left[u^{(s)}{\vartheta }_{,x}^{(2n-s)}+v^{(s)}{\vartheta }_{,y}^{(2n-s)}+w^{(s)}{\vartheta }^{(2n-s+1)}\right]\\ \hfill =R{w}^{(2n)}+{\Delta }^{*}{\vartheta }^{(2n)}+{\vartheta }^{(2n+2)}.\end{array}$$

Now, upon setting $n=1$, we have

$$\begin{array}{c}\hfill {\vartheta }_{,t}^{(2)}+\sum _{s=0}^2\left(\genfrac{}{}{0pt}{}{2}{s}\right)\left[u^{(s)}{\vartheta }_{,x}^{(2-s)}+v^{(s)}{\vartheta }_{,y}^{(2-s)}+w^{(s)}{\vartheta }^{(3-s)}\right]\\ \hfill =R{w}^{(2)}+{\Delta }^{*}{\vartheta }^{(2)}+{\vartheta }^{(4)}.\end{array}$$

(41)

Thus, employing Equations (14), (36), (39) and (40) yields

$$\begin{array}{c}{\vartheta }^{(4)}=0,z=0,1.\hfill \end{array}$$

(42)

We next differentiate Equation (33), an even number of times with respect to z, to find

$$\begin{array}{c}{\omega }_{1,zz}=\tilde{T}{u}_{,zzz},{\omega }_{2,zz}=\tilde{T}{v}_{,zzz},z=0,1.\hfill \end{array}$$

(43)

In addition, we also differentiate Equation (35) an even number of times with respect to z, we have

$$\begin{array}{c}{\omega }_{1,zz}=-v_{,zzz},{\omega }_{2,zz}=u_{,zzz},z=0,1.\hfill \end{array}$$

(44)

Therefore, from Equations (43) and (44), we obtain

$$\begin{array}{c}{u}_{,zzz}=0,v_{,zzz}=0,z=0,1.\hfill \end{array}$$

(45)

By further differentiation of Equation (38) an even number of times with respect to z, we find

$$\begin{array}{c}{w}^{(4)}=0,z=0,1.\hfill \end{array}$$

(46)

The above process may be repeated to derive the general boundary conditions

$$\begin{array}{c}{w}^{(2n)}=0,{\vartheta }^{(2n)}=0,z=0,1,\mathrm{for}n\in \mathbb{N},\hfill \end{array}$$

(47)

which hold for the solution of the nonlinear problem.

We aim now to study nonlinear energy stability and find a stability threshold. Again, we let V be a periodic cell for a disturbance to Equation (9), and let $∥.∥$ and $(., .)$ be the norm and inner product on $L^2(V)$. The energy identities are derived by multiplying the vertical component of Equation (11) by w, upon use of Equations (16) and (20) with $i=3$, and also use some integrations by parts, with the aid of boundary conditions, one may show that

$$\begin{array}{c}\xi {∥{\nabla }^{*}w∥}^2+\left(1+{\tilde{T}}^2\right){∥w_{,z}∥}^2=R\left({\nabla }^{*}\vartheta ,{\nabla }^{*}w\right).\hfill \end{array}$$

(48)

Next, multiply Equation (9c) by $\vartheta$ and integrate over V to find

$$\begin{array}{c}\frac{1}{2}\frac{d}{dt}{∥\vartheta ∥}^2=R\left(w,\vartheta \right)-{∥\nabla \vartheta ∥}^2.\hfill \end{array}$$

(49)

By adding $\lambda$ (48) to (49), for $\lambda >0$ a parameter to be chosen, we may derive an energy identity of form

$$\begin{array}{c}\frac{dE}{dt}=RI-D,\hfill \end{array}$$

(50)

where

$$\begin{array}{c}{\displaystyle E(t)=\frac{1}{2}{∥\vartheta ∥}^2,}\end{array}$$

(51)

$$\begin{array}{c}I=\left(w,\vartheta \right)+\lambda \left({\nabla }^{*}\vartheta ,{\nabla }^{*}w\right),\end{array}$$

(52)

$$\begin{array}{c}D={∥\nabla \vartheta ∥}^2+\lambda \left(\xi {∥{\nabla }^{*}w∥}^2+\left(1+{\tilde{T}}^2\right){∥w_{,z}∥}^2\right).\end{array}$$

(53)

Define $R_E$ by

$$\frac{1}{R_E}=\underset{\mathcal{H}}{max}\frac{I}{D},$$

(54)

where $\mathcal{H}$ is the space of admissible functions given by

$\mathcal{H}=\{u_i,\vartheta |u_i\in L^2(V),\vartheta \in H^1(V),u_{i,i}=0,u_i,\vartheta$ are periodic in $x,y\}$.

Therefore, from Equation (50), we deduce

$$\frac{dE}{dt}\le -D\left(\frac{R_E-R}{R_E}\right).$$

(55)

Then, from the Poincaré’s inequality on D, we have

$$D\ge {\pi }^2{∥\vartheta ∥}^2.$$

Provided $R<R_E$, put $c=1-R/R_E>0$ and then, from Equation (55), we have

$$\frac{dE}{dt}\le -2{\pi }^{2}cE(t).$$

This yields

$$E(t)\le E(0)e^{-2{\pi }^{2}ct}.$$

Thus, $E(t)$ tends to 0 as $t\to \infty$ at least exponentially. Therefore, $∥\vartheta (t)∥\to 0$ at least exponentially.

To obtain the decay of $\mathbf{u}$, we multiply Equation (9) by $u_i$ and integrate over V to obtain

$$\begin{array}{c}\left(m_{ij}{u}_j,u_i\right)=R\left(\vartheta ,w\right).\hfill \end{array}$$

(56)

We may observe that

$$\left(m_{ij}{u}_j,u_i\right)\ge \widehat{\mu }{∥\mathbf{u}∥}^2,$$

where

$$\widehat{\mu }=min\left\{1,\xi \right\}.$$

From Equation (56), it now follows that

$$\widehat{\mu }{∥\mathbf{u}∥}^2\le R\left(\vartheta ,w\right),$$

and then with use of the arithmetic geometric mean inequality, one shows

$$\begin{array}{ccc}\hfill \widehat{\mu }{∥\mathbf{u}∥}^2& \le & \frac{R}{2\widehat{\alpha }}{∥\vartheta ∥}^2+\frac{R\widehat{\alpha }}{2}{∥w∥}^2\hfill \\ & \le & \frac{R}{2\widehat{\alpha }}{∥\vartheta ∥}^2+\frac{R\widehat{\alpha }}{2}{∥\mathbf{u}∥}^2,\hfill \end{array}$$

for $\widehat{\alpha }>0$ to be chosen.

If we now pick $\widehat{\alpha }=\widehat{\mu }/R$, then we show

$$0<{∥\mathbf{u}∥}^2\le \frac{R^2}{{\widehat{\mu }}^2}{∥\vartheta ∥}^2,$$

which implies ${∥\mathbf{u}∥}^2$ must also decay at least exponentially. Hence, the global nonlinear stability criterion is determined by Equation (54).

In order to determine $R_E$, we have to derive the Euler–Lagrange equations and maximise in the coupling parameter $\lambda$. To do this, we must find the stationary point of $I/D$, by using the calculus of the variations technique, the Euler–Lagrange equations arising from Equation (54) are determined from

$$\begin{array}{c}{R}_E\delta I-\delta D=0,\hfill \end{array}$$

(57)

for all $h_i\in \mathcal{H}$, and $\eta \in \mathcal{H}$. We have that

$$\begin{array}{c}\delta D=\frac{d}{d\epsilon }{\left.\underset{V}{\int }\left[{\left(\nabla \left(\vartheta +\eta \epsilon \right)\right)}^2+\lambda \xi {\left({\nabla }^{*}\left(w+h_3\epsilon \right)\right)}^2+\lambda \left(1+{\tilde{T}}^2\right){\left(w_{,z}+h_{3,z}\epsilon \right)}^2\right]dV\right|}_{\epsilon =0},\hfill \end{array}$$

$$\begin{array}{c}={\left.\underset{V}{\int }\left[2\nabla \left(\vartheta +\eta \epsilon \right)\nabla \eta +2\lambda \xi {\nabla }^{*}\left(w+h_3\epsilon \right){\nabla }^{*}{h}_3+2\lambda \left(1+{\tilde{T}}^2\right)\left(w_{,z}+h_{3,z}\epsilon \right)h_{3,z}\right]dV\right|}_{\epsilon =0},\hfill \end{array}$$

and

$$\begin{array}{c}\delta I=\frac{d}{d\epsilon }{\left.\underset{V}{\int }\left[\left(w+\epsilon h_3\right)\left(\vartheta +\epsilon \eta \right)+\lambda {{\nabla }^{*}\left(\vartheta +\eta \epsilon \right)\nabla }^{*}\left(w+h_3\epsilon \right)-\left(u_{i,i}+\epsilon h_{i,i}\right)\pi (x)\right]dV\right|}_{\epsilon =0},\hfill \end{array}$$

$$\begin{array}{c}={\left.\underset{V}{\int }\left[\left(w+\epsilon h_3\right)\eta +h_3\left(\vartheta +\epsilon \eta \right)+\lambda {{\nabla }^{*}\left(\vartheta +\eta \epsilon \right)\nabla }^{*}{h}_3+{\nabla }^{*}\left(w+h_3\epsilon \right){\nabla }^{*}\eta -h_{i,i}\pi (x)\right]dV\right|}_{\epsilon =0},\hfill \end{array}$$

where we have included the constraint $u_{i,i}=0$ by way of a Lagrange multiplier $2\pi (x)$, and $\epsilon$ is a positive constant.

Furthermore, after some integrations by parts and using the boundary conditions, we find that

$$\begin{array}{ccc}\hfill \delta D& =\hfill & \hfill \underset{V}{\int }\left[-2\eta \Delta \vartheta +2\lambda h_3\left(-\xi {\Delta }^{*}w-\left(1+{\tilde{T}}^2\right)w_{,zz}\right)\right]dV,\end{array}$$

$$\begin{array}{ccc}\hfill \delta I& =\hfill & \hfill \underset{V}{\int }\left[\eta \left(w-\lambda {\Delta }^{*}w\right)+h_i\left({\delta }_{i3}\left(\vartheta -\lambda {\Delta }^{*}\vartheta \right)-{\pi }_{,i}\right)\right]dV.\end{array}$$

Since $h_i$ and $\eta$ were chosen arbitrary functions, from Equation (57), we obtain the Euler–Lagrange equations

$$\begin{array}{c}{R}_E\left(\vartheta -\lambda {\Delta }^{*}\vartheta \right)+2\lambda \left(\xi {\Delta }^{*}w+\left(1+{\tilde{T}}^2\right)w_{,zz}\right)={\pi }_{,i},\end{array}$$

(58)

$$\begin{array}{c}{R}_E\left(w-\lambda {\Delta }^{*}w\right)+2\Delta \vartheta =0,\end{array}$$

(59)

where $\pi (\mathbf{x})$ is now a Lagrange multiplier. Applying the horizontal Laplacian operator to Equation (58), we obtain

$$\begin{array}{c}{R}_E\left(\lambda {\Delta }^{*}-1\right){\Delta }^{*}\vartheta -2\lambda {\Delta }^{*}\left(\xi {\Delta }^{*}w+\left(1+{\tilde{T}}^2\right)w_{,zz}\right)=0,\end{array}$$

(60)

$$\begin{array}{c}{R}_E\left(w-\lambda {\Delta }^{*}w\right)+2\Delta \vartheta =0.\end{array}$$

(61)

We again use a normal mode representation, as for the linear stability analysis, $\vartheta =\mathsf{\Theta }(z)f(x,y)$, $w=W(z)f(x,y)$. This leaves us to solve the eigenvalue problem

$$\begin{array}{c}{R}_E\left(1+\lambda a^2\right)\mathsf{\Theta }+2\lambda \left[\left(1+{\tilde{T}}^2\right)D^2-\xi a^2\right]W=0,\hfill \\ R_E\left(1+\lambda a^2\right)W+2\left(D^2-a^2\right)\mathsf{\Theta }=0.\hfill \end{array}$$

(62)

This system would have to be solved for $R_E$ subject to the boundary conditions Equation (27). Furthermore, we observe that W and $\mathsf{\Theta }$ satisfy the boundary conditions

$$\begin{array}{c}{W}^{(2n)}=0,{\mathsf{\Theta }}^{(2n)}=0,z=0,1,forn\in \mathbb{N}.\hfill \end{array}$$

(63)

By eliminating $\mathsf{\Theta }$, we obtain a fourth order equation in W,

$$\begin{array}{c}4\lambda \left(1+{\tilde{T}}^2\right)\left(D^2-a^2\right)D^{2}W-4\lambda \xi a^2\left(D^2-a^2\right)W=R_E^2{\left(1+\lambda a^2\right)}^{2}W.\hfill \end{array}$$

(64)

Hence, $W(z)$ may be written in the form

$$W= \sin n\pi z,n=1,2,\cdots .$$

After some calculations, following the method in Section 4, one may find

$$\begin{array}{c}{R}_E^2=4\lambda \frac{{\pi }^2{n}^2(1+{\tilde{T}}^2)({\pi }^2{n}^2+a^2)+\xi a^2({\pi }^2{n}^2+a^2)}{{(1+\lambda a^2)}^2}.\hfill \end{array}$$

(65)

For any fixed wave number $a^2$, the minimum with respect to $n^2$ of $R_E^2(a^2,n^2)$ is obtained for $n=1$. Then,

$$\begin{array}{c}{R}_E^2=4\lambda \frac{{\pi }^2(1+{\tilde{T}}^2)({\pi }^2+a^2)+\xi a^2({\pi }^2+a^2)}{{(1+\lambda a^2)}^2}.\hfill \end{array}$$

(66)

Let us now select $\lambda =1/a^2$, and then

$$\begin{array}{c}{R}_E^2=\frac{{\pi }^2({\pi }^2+a^2)(1+{\tilde{T}}^2)}{a^2}+\xi ({\pi }^2+a^2).\hfill \end{array}$$

(67)

This is exactly the same Equation (30) with $n=1$ for linear instability problem. This is, in a sense, the best possible threshold for the onset of linear unconditional stability. Thus, the minimum of $R_E^2$ with respect to $a^2$ is identical to the minimum of $R_L^2$ with respect to $a^2$, and hence no subcritical instabilities can arise. This result is undoubtedly due to the fact that the operator attached to the linear theory is symmetric in this case (see Straughan [27] and Falsaperla et al. [28]).

6. Numerical Results

The aim of this paper was to investigate how the inclusion of the Taylor number ${\tilde{T}}^2$ affects the thermal instability threshold in an anisotropic porous medium. The results of different values of the anisotropy parameter $\xi$ and the Taylor number ${\tilde{T}}^2$ are presented in

Table 1. Critical values of Rayleigh number $R_c$, vs. $\xi$, for ${\tilde{T}}^2=5,10,15,20,25$.

	$R_c$
$\mathit{\xi }$	${\tilde{\mathit{T}}}^{\mathbf{2}}=\mathbf{0}$	${\tilde{\mathit{T}}}^{\mathbf{2}}=\mathbf{5}$	${\tilde{\mathit{T}}}^{\mathbf{2}}=\mathbf{10}$	${\tilde{\mathit{T}}}^{\mathbf{2}}=\mathbf{15}$	${\tilde{\mathit{T}}}^{\mathbf{2}}=\mathbf{20}$	${\tilde{\mathit{T}}}^{\mathbf{2}}=\mathbf{25}$	${\tilde{\mathit{T}}}^{\mathbf{2}}=\mathbf{100}$
1	39.478	117.438	183.903	246.740	307.588	367.130	1205.076
2	57.524	147.335	220.890	289.315	354.926	418.690	1297.116
3	73.668	172.573	251.568	324.280	393.546	460.550	1370.037
4	88.826	195.398	278.979	355.306	427.653	497.389	1433.062
5	103.356	216.682	304.304	383.815	458.876	531.019	1489.762
6	117.438	236.871	328.145	410.535	488.051	562.370	1541.969
7	131.182	256.230	350.864	435.901	515.674	591.993	1590.772
8	144.657	274.932	372.693	460.194	542.068	620.249	1636.881
9	157.914	293.097	393.795	483.611	567.457	647.388	1680.786
10	170.987	310.813	414.288	506.293	592.006	673.591	1722.848

and

Table 2. Critical values of wave number $a_c$, vs. $\xi$, for ${\tilde{T}}^2=0,5,10,15,20,25,100$.

	$a_c$
$\mathit{\xi }$	${\tilde{\mathit{T}}}^{\mathbf{2}}=\mathbf{0}$	${\tilde{\mathit{T}}}^{\mathbf{2}}=\mathbf{5}$	${\tilde{\mathit{T}}}^{\mathbf{2}}=\mathbf{10}$	${\tilde{\mathit{T}}}^{\mathbf{2}}=\mathbf{15}$	${\tilde{\mathit{T}}}^{\mathbf{2}}=\mathbf{20}$	${\tilde{\mathit{T}}}^{\mathbf{2}}=\mathbf{25}$	${\tilde{\mathit{T}}}^{\mathbf{2}}=\mathbf{100}$
1	3.142	4.917	5.721	6.283	6.725	7.094	9.959
2	2.642	4.135	4.811	5.284	5.655	5.965	8.375
3	2.387	3.736	4.347	4.774	5.110	5.390	7.567
4	2.221	3.477	4.046	4.443	4.755	5.016	7.042
5	2.101	3.288	3.826	4.202	4.497	4.744	6.660
6	2.007	3.142	3.656	4.015	4.297	4.533	6.363
7	1.931	3.023	3.517	3.863	4.135	4.361	6.123
8	1.868	2.924	3.402	3.736	3.999	4.218	5.922
9	1.814	2.839	3.303	3.628	3.883	4.096	5.750
10	1.767	2.765	3.217	3.533	3.782	3.989	5.601

, and are presented graphically in Figure 1, Figure 2 and Figure 3.

Table 1 and Figure 1 present the values of $R_{L(sc)}^2=R_c$, the critical Rayleigh number for both the onset of linear instability and for the nonlinear stability. This shows that the effect of increasing the Taylor number ${\tilde{T}}^2$ always results in an increase in the critical Rayleigh number $R_c$, so that rotation stabilizes the system. Furthermore, the effect of increasing the anisotropy parameter $\xi$ is seen also to increase the critical Rayleigh number $R_c$. This means that, when the rate of rotation and $\xi$ increase, the stability becomes more pronounced, i.e., $R_c$ increases. For example, for $\xi =3$ and ${\tilde{T}}^2=5$, we see from Table 1 that the critical Rayleigh number is $R_c=172.573$, whereas, when $\xi =10$ and ${\tilde{T}}^2=25,$ the critical Rayleigh number is $R_c=673.591$.

From Figure 3, it is evident that, when there is no rotation, ${\tilde{T}}^2=0$, the instability curve starts at $R_c=4{\pi }^2$ when $\xi =1$ and increases when rotation is included. Note that the critical Rayleigh number $R_c$ became significantly higher at ${\tilde{T}}^2=100$, which leads to stabilize the system. One can see that, when $\xi$ increases to $\xi =5$ in case ${\tilde{T}}^2=0$, the instability curve starts at $R_c=103.356$. This means that the effect of increasing the anisotropy parameter is to delay the onset of convection in a fluid layer. Again, we observe that increasing the Taylor number ${\tilde{T}}^2$ leads to an increase in the critical Rayleigh number $R_c$. Thus, an increase in the anisotropy parameter $\xi$ in the vertical direction with an increase in the Taylor number ${\tilde{T}}^2$ has the effect of stabilizing the system. We can, therefore, conclude that the effect of rotation is to enhance the stability of the system. In addition, these results are reinforcing the fact that the linear instability analysis is accurately capturing the physics of the onset of convection.

Table 2 and Figure 2 present the values of $a_{L(c)}=a_c$, the critical wave number for both the onset of linear instability and for the nonlinear stability. It can be observed that, for a fixed value of the anisotropy parameter $\xi$, the effect of increasing the Taylor number ${\tilde{T}}^2$ is to increase the wave number. For example, for $\xi =3$ and ${\tilde{T}}^2=5$, we see from Table 2 that the critical wave number is $a_c=3.736$, whereas, when ${\tilde{T}}^2=25$ for the same anisotropy parameter $\xi =3$, the critical wave number is $a_c=5.390$. It is also observed that increasing the anisotropy parameter $\xi$ had the effect of decreasing the value of the wave number. However, as soon as the value of the Taylor number ${\tilde{T}}^2$ increases, one can observe the critical wave number also increases, which corresponds to the narrower convection cells. These results indicate the effect of incorporating rotation in an anisotropic porous medium.

7. Conclusions

In this article, we investigated the combined effects of the Taylor number ${\tilde{T}}^2$ and the anisotropy parameter $\xi$ on the stability threshold for the thermal convection problem. We have studied a model of the thermal convection in a fluid saturated rotating anisotropic Darcy medium allowing the Vadasz number ${\widehat{V}}_a$ to be infinite. The validity of the linear instability is tested and the nonlinear analysis has performed to confirm the validity of linear instability. Our analysis emphasized that the subcritical instabilities are not possible when the inertia term is neglected and hence the linear instability analysis is accurately capturing the physics of the onset of convection. These results showed that the effect of rotation is to enhance the stability of the system.