The onset of thermal convection in anisotropic and rotating bidisperse porous media

Abstract

The onset of thermal convection in anisotropic rotating bidisperse porous media is investigated. The optimal result concerning the coincidence between linear instability and nonlinear stability thresholds with respect the $L^2$-norm is obtained.

1 Introduction

The onset of thermal convection for both clear fluids (see [7, 20]) and fluids saturating porous materials (see [19, 25]) has been widely studied by many scientists—in the past as nowadays—due to its relevant applications in different fields, such as geophysics (geothermal reservoirs, geological storage of carbon dioxide), astrophysics (pore water convection within carbonaceous chondrite parent bodies), engineering and industrial process (water treatment process, nuclear waste disposal, chemical reactor engineering, and the storage of heat-generating materials such as grain and coal) (see [14] and references therein). In [1, 19] numerical studies of fluid flows through porous media have been performed; then in [5, 6, 12] the onset of convection in anisotropic porous media under various assumptions has been analyzed. Moreover, [13, 31] investigated clear fluid flows with pressure dependent viscosity, while [21,22,23] deal with the same problem but through porous media. Finally, [11, 20, 24, 30] analyze the onset of convection of a binary mixture in a porous medium, of a non-Newtonian fluid between two vertical plate, of a fluid saturating a porous medium under LTNE assumption and of fluids saturating rotating porous media, respectively.

However, in recent times, double-porosity materials have attracted the attention of a large number of researchers. A double-porosity material is also referred to as bidisperse porous medium (BDPM): it has the normal pore structure, but the solid skeleton has cracks or fissures in it. In particular, a BDPM is composed of clusters of large particles that are agglomerations of small particles, there are macropores between the clusters and micropores within them, and in particular the macropores are referred to as f-phase, while the remainder of the structure is referred to as p-phase. The reason of this nomenclature is that one can think of the f-phase as being a fracture phase, the p-phase as a porous phase [17].

Let $\varphi$ be the porosity of the macropores and $\epsilon$ be the porosity of the micropores, thus $(1- \varphi ) \epsilon$ is the fraction of volume occupied by the micropores, $\varphi + (1- \varphi ) \epsilon$ is the fraction of volume occupied by the fluid, $(1- \epsilon )(1-\varphi )$ is the fraction of volume occupied by the solid skeleton.

The fundamental theory for thermal convection in bidisperse porous media can be found in [15,16,17].

The study of bidisperse convection has a large number of practical applications, such as industrial ones, for example in order to design heat pipes (as reported in [14], since the bidisperse wick structure significantly increases the area available for liquid film evaporation, it has been proposed for use in the evaporator of heat pipes), or medical ones, in fact brain and human bones may be modeled as bidisperse porous media (in particular, the analysis of the onset convection in a BDPM allows to understand and model interstitial fluid flow in bones or to develop tissue engineering strategy for bone defects [25, 26]). Since in bidisperse porous media the convection occurs at higher Rayleigh numbers, the heat transfer in BDPM due to the convective movement is delayed, and for this reason BDPM are significantly better used in thermal insulation problems and thermal management problem (such as cooling of data centers) than simple porous materials [8, 9, 17].

As reported in [27], bidisperse porous media are increasingly important in the chemical engineering field and regarding anisotropic materials, while anisotropic single porosity media have been widely studied by several authors (see for example [5, 12, 18]), anisotropic bidisperse porous materials may have much more potentials, since they offer many possibilities to design manmade materials for heat transfer or insulation problems, for oil recovery from underground reservoir, for nuclear waste recovery and so on (see [3, 10, 28, 29] and references therein). Therefore, in the present, we allow fully anisotropic permeabilities in both f-phase and p-phase.

The key role of rotation on the onset of thermal convection in porous materials has been pointed out by many authors in view of its practical applications in geophysics and in engineering (food process industry, chemical process industry, centrifugal filtration processes, rotating machinery) [30], hence the study of thermal convection in rotating BDPM may be necessary and useful as well, see for instance [2, 4], which deal with the onset of convection in a horizontal layer of rotating isotropic BDPM, according to Darcy’s law and taking into account the Vadasz term, respectively, or [3], which analyzes the onset of convection in a horizontally isotropic rotating BDPM.

Envisaging a rotating machinery constituted by an engineered fully anisotropic bidisperse porous material, the aim of this paper is to analyze the onset of thermal convection in an anisotropic bidisperse porous medium uniformly rotating about a vertical axis, through linear and nonlinear stability theory. The plan of the paper is as follows. The mathematical model and the associated perturbation equations are introduced in Sect. 2. In Sect. 3 the strong version of the principle of exchange of stabilities is proved and the linear instability analysis of the thermal conduction solution is performed, to determine the steady instability threshold. In Sect. 4, the nonlinear stability analysis of the thermal conduction solution is performed, proving the coincidence between the linear and the nonlinear stability thresholds with respect to the $L^2$-norm. In Sect. 5, in order to analyze the influence of rotation and of anisotropic permeability on the onset of convection, numerical simulations are presented. The paper ends with a conclusions section, in which all the results are collected.

2 The mathematical model

Let Oxyz be a reference frame with fundamental unit vectors $\mathbf{i},\mathbf{j},\mathbf{k}$ ($\mathbf{k}$ pointing vertically upward) and let L be a layer of an anisotropic bidisperse porous medium, saturated by an homogeneous incompressible fluid heated from below. Let us assume that the layer L—of thickness d—rotates about the vertical axis z, with constant angular velocity ${\varvec{\Omega }}=\Omega \mathbf{k}$ and that the temperature in the macropores $(T_f)$ and the temperature in the micropores $(T_p)$ are the same, i.e., $T^f=T^p=T$.

Let us assume that the axes (x, y, z) are the principal axes of the permeability, so the permeability tensors in the macropores and in the micropores may be written as

$$\begin{aligned}&\begin{aligned} \mathbf{K}^f&= \text {diag}(K^f_x,K^f_y,K^f_z)=K^f_z \ \mathbf{K}^{f*}, \\ \mathbf{K}^p&= \text {diag}(K^p_x,K^p_y,K^p_z)=K^p_z \ \mathbf{K}^{p*}, \end{aligned} \\&\mathbf{K}^{f*}=\text {diag}(k_1,k_2,1), \quad \mathbf{K}^{p*}=\text {diag}(h_1,h_2,1), \end{aligned}$$

where

$$\begin{aligned} \begin{array}{l} k_1=\dfrac{K^f_x}{K^f_z}, \quad k_2=\dfrac{K^f_y}{K^f_z}, \\ \\ h_1=\dfrac{K^p_x}{K^p_z}, \quad h_2=\dfrac{K^p_y}{K^p_z}. \end{array} \end{aligned}$$

In the Oberbeque–Boussinesq approximation and extending the Darcy’s Law in order to include the Coriolis term in the momentum equations for the micropores and the macropores, the governing equations for thermal convection are [2, 27]:

$$\begin{aligned} {\left\{ \begin{array}{ll} \mathbf{v}^f = \dfrac{1}{\mu } \mathbf{K}^f \cdot \Bigl [ - \zeta (\mathbf{v}^f - \mathbf{v}^p) - \nabla p^f + \varrho _F \alpha g T \mathbf{k} - \dfrac{2 \varrho _F \Omega }{\varphi } \mathbf{k} \times \mathbf{v}^f \Bigr ], \\ \mathbf{v}^p = \dfrac{1}{\mu } \mathbf{K}^p \cdot \Bigl [ - \zeta (\mathbf{v}^p - \mathbf{v}^f) - \nabla p^p + \varrho _F \alpha g T \mathbf{k} - \dfrac{2 \varrho _F \Omega }{\epsilon } \mathbf{k} \times \mathbf{v}^p \Bigr ], \\ \nabla \cdot \mathbf{v}^f = 0, \\ \nabla \cdot \mathbf{v}^p = 0, \\ (\varrho c)_m \dfrac{\partial T}{\partial t} + (\varrho c)_f (\mathbf{v}^f + \mathbf{v}^p) \cdot \nabla T = k_m \Delta T. \end{array}\right. } \end{aligned}$$

(1)

where

$$\begin{aligned} p^s=P^s-\frac{\varrho _F}{2} \vert {\varvec{\Omega }} \times \mathbf{x} \vert ^2, \quad s=\{f,p\} \end{aligned}$$

are the reduced pressures, $\mathbf{x}=(x,y,z)$, $\mathbf{v}^s$ = seepage velocity for $s=\{f,p\}$, $\zeta$ = interaction coefficient between the f-phase and the p-phase, $\mathbf{g}=-g \mathbf{k}$ = gravity, $\mu$ = fluid viscosity, $\varrho _F$ = reference constant density, $\alpha$ = thermal expansion coefficient, c = specific heat, $c_p$ = specific heat at a constant pressure, $(\varrho c)_m=(1-\varphi )(1-\epsilon )(\varrho c)_{sol}+\varphi (\varrho c)_f+\epsilon (1-\varphi )(\varrho c)_p$, $k_m=(1-\varphi )(1-\epsilon )k_{sol}+\varphi k_f+\epsilon (1-\varphi )k_p$ = thermal conductivity (the subscript sol is referred to the solid skeleton).

To (1), the following boundary conditions are appended

$$\begin{aligned} \begin{aligned} \mathbf{v}^s \cdot \mathbf{n} = 0&\, \text {on} \ z=0,d, \ s=\{f,p\} \\ T=T_L&\, \text {on} \ z=0,\\ T=T_U&\, \text {on} \ z=d, \end{aligned} \end{aligned}$$

(2)

where $\mathbf{n}$ is the unit outward normal to the impermeable horizontal planes delimiting the layer and $T_L>T_U$.

The problem (1)-(2) admits the steady state (conduction solution):

$$\begin{aligned} \mathbf{{\overline{v}}}^f=\mathbf{0}, \ \mathbf{{\overline{v}}}^p=\mathbf{0}, \ {\overline{T}}=- \beta z + T_L, \ \beta =\dfrac{T_L-T_U}{d}. \end{aligned}$$

Defining $\{ \mathbf{u}^f, \mathbf{u}^p, \theta , \pi ^f, \pi ^p \}$ a perturbation to the steady solution, the evolution equations for the perturbation fields are

$$\begin{aligned} {\left\{ \begin{array}{ll} \mathbf {u}^f = \dfrac{1}{\mu } \mathbf {K}^f \cdot \Bigl [ - \zeta (\mathbf {u}^f - \mathbf {u}^p) - \nabla \pi _f + \varrho _F \alpha g \theta \mathbf {k} - \dfrac{2 \varrho _F \Omega }{\varphi } \mathbf {k} \times \mathbf {u}^f \Bigr ], \\ \mathbf {u}^p = \dfrac{1}{\mu } \mathbf {K}^p \cdot \Bigl [ - \zeta (\mathbf {u}^p - \mathbf {u}^f) - \nabla \pi _p + \varrho _F \alpha g \theta \mathbf {k} - \dfrac{2 \varrho _F \Omega }{\epsilon } \mathbf {k} \times \mathbf {u}^p \Bigr ], \\ \nabla \cdot \mathbf {u}^f = 0, \\ \nabla \cdot \mathbf {u}^p = 0, \\ (\varrho c)_m \dfrac{\partial \theta }{\partial t} + (\varrho c)_f (\mathbf {u}^f + \mathbf {u}^p) \cdot \nabla \theta = (\varrho c)_f \beta (w^f+w^p) + k_m \Delta \theta . \end{array}\right. } \end{aligned}$$

(3)

where $\mathbf{u}^f=(u^f,v^f,w^f), \ \mathbf{u}^p=(u^p,v^p,w^p)$. Using the following non-dimensional parameters

$$\begin{aligned} \begin{aligned}&\mathbf{x}^{*} = \frac{\mathbf{x}}{d}, \ t^{*}=\frac{t}{{\tilde{t}}}, \ \theta ^{*} = \frac{\theta }{{\tilde{T}}}, \\&\mathbf{u}^{s*} = \frac{\mathbf{u}^s}{{\tilde{u}}}, \ \pi ^{s*} = \frac{\pi ^s}{{\tilde{P}}}, \ \text {for} \ s=\{f,p\}, \\&\eta =\frac{\varphi }{\epsilon }, \ \gamma _1=\frac{\mu }{K^f_z \zeta }, \ \gamma _2=\frac{\mu }{K^p_z \zeta }, \end{aligned} \end{aligned}$$

where the scales are given by

$$\begin{aligned} \small { {\tilde{u}}=\frac{k_m}{(\varrho c)_f d}, \ {\tilde{t}} = \frac{d^2 (\varrho c)_m}{k_m}, \ {\tilde{P}} = \frac{\zeta k_m}{(\varrho c)_f}, \ {\tilde{T}} = \sqrt{\frac{\beta k_m \zeta }{(\varrho c)_f \varrho _F \alpha g}} } \end{aligned}$$

and introducing the Taylor number ${\mathcal {T}}$ and the thermal Rayleigh number R, respectively, given by

$$\begin{aligned} {\mathcal {T}}= \frac{2 \varrho _F \Omega K^f_z}{\varphi \mu }, \qquad R = \sqrt{\frac{\beta d^2 (\varrho c)_f \varrho _F \alpha g}{k_m \zeta }}, \end{aligned}$$

the resulting non-dimensional perturbation equations, omitting all the asterisks, are

$$\begin{aligned} {\left\{ \begin{array}{ll} \gamma _1 (\mathbf{K}^f)^{-1} \mathbf{u}^f + (\mathbf{u}^f - \mathbf{u}^p) = - \nabla \pi ^f + R \theta \mathbf{k} - \gamma _1 {\mathcal {T}} \mathbf{k} \times \mathbf{u}^f, \\ \gamma _2 (\mathbf{K}^p)^{-1} \mathbf{u}^p - (\mathbf{u}^f - \mathbf{u}^p) = - \nabla \pi ^p + R \theta \mathbf{k} - \eta \gamma _1 {\mathcal {T}} \mathbf{k} \times \mathbf{u}^p, \\ \nabla \cdot \mathbf{u}^f = 0, \\ \nabla \cdot \mathbf{u}^p = 0, \\ \dfrac{\partial \theta }{\partial t} + (\mathbf{u}^f + \mathbf{u}^p) \cdot \nabla \theta = R (w^f+w^p) +\Delta \theta , \end{array}\right. } \end{aligned}$$

(4)

under the initial conditions

$$\begin{aligned} \mathbf{u}^s(\mathbf{x},0)=\mathbf{u}^s_0(\mathbf{x}), \quad \pi ^s(\mathbf{x},0)=\pi ^s_0(\mathbf{x}), \quad \theta (\mathbf{x},0)=\theta _0(\mathbf{x}), \end{aligned}$$

with $\nabla \cdot \mathbf{u}^s_0=0$, for $s=\{f,p\}$, and the boundary conditions

$$\begin{aligned} w^f=w^p= \theta =0 \quad \text {on} \ z=0,1. \end{aligned}$$

(5)

The above equations are defined in $\{ (x,y,z,t) \in {\mathbb {R}}^4 | z \in (0,1), t>0 \}$.

In the sequel, we will suppose that the perturbation fields are periodic in the x and y directions of period $\displaystyle {\frac{2 \pi }{l}}$ and $\displaystyle {\frac{2 \pi }{m}}$, respectively, and we will denote by

$$\begin{aligned} V=\Big [ 0,\frac{2 \pi }{l} \Big ] \times \Big [ 0,\frac{2 \pi }{m} \Big ] \times [0,1] \end{aligned}$$

the periodicity cell.

3 Instability analysis

In this section we will perform linear instability analysis of the basic solution, to this aim let us first linearise the perturbation Eq. (4), i.e.,

$$\begin{aligned} {\left\{ \begin{array}{ll} \gamma _1 (\mathbf{K}^f)^{-1} \mathbf{u}^f + (\mathbf{u}^f - \mathbf{u}^p) = - \nabla \pi ^f + R \theta \mathbf{k} - \gamma _1 {\mathcal {T}} \mathbf{k} \times \mathbf{u}^f, \\ \gamma _2 (\mathbf{K}^p)^{-1} \mathbf{u}^p - (\mathbf{u}^f - \mathbf{u}^p) = - \nabla \pi ^p + R \theta \mathbf{k} - \eta \gamma _1 {\mathcal {T}} \mathbf{k} \times \mathbf{u}^p, \\ \nabla \cdot \mathbf{u}^f = 0, \\ \nabla \cdot \mathbf{u}^p = 0, \\ \dfrac{\partial \theta }{\partial t} = R (w^f+w^p) +\Delta \theta , \end{array}\right. } \end{aligned}$$

(6)

Since the system (6) is autonomous, we seek solutions which have time-dependence like $e^{\sigma t}$, i.e., solutions of form

$$\begin{aligned} \begin{aligned} \mathbf{u}^s(t,\mathbf{x})&= e^{\sigma t} \mathbf{u}^s(\mathbf{x}), \\ \theta (t,\mathbf{x})&= e^{\sigma t} \theta (\mathbf{x}), \\ \pi ^s(t,\mathbf{x})&= e^{\sigma t} \pi ^s(\mathbf{x}), \end{aligned} \end{aligned}$$

(7)

with $\sigma \in {\mathbb {C}}$ and $s=\{f,p\}$. By virtue of (7), (6) becomes

$$\begin{aligned} {\left\{ \begin{array}{ll} \gamma _1 (\mathbf{K}^f)^{-1} \mathbf{u}^f + (\mathbf{u}^f - \mathbf{u}^p) = - \nabla \pi ^f + R \theta \mathbf{k} - \gamma _1 {\mathcal {T}} \mathbf{k} \times \mathbf{u}^f, \\ \gamma _2 (\mathbf{K}^p)^{-1} \mathbf{u}^p - (\mathbf{u}^f - \mathbf{u}^p) = - \nabla \pi ^p + R \theta \mathbf{k} - \eta \gamma _1 {\mathcal {T}} \mathbf{k} \times \mathbf{u}^p, \\ \sigma \theta = R (w^f+w^p) +\Delta \theta . \end{array}\right. } \end{aligned}$$

(8)

Let $^*$ be the complex conjugate of a field. Multiplying $(8)_1$ by $\mathbf{u}^{f*}$, $(8)_2$ by $\mathbf{u}^{p*}$ and $(8)_3$ by $\theta ^{*}$, integrating each equation over the periodic cell V and adding the resulting equations, one obtains

$$\begin{aligned} \begin{aligned} \sigma \Vert \theta \Vert ^2&= - \gamma _1 (\mathbf{M}^f \mathbf{u}^f,\mathbf{u}^{f*}) - \gamma _2 (\mathbf{M}^p \mathbf{u}^p,\mathbf{u}^{p*}) \\&\quad - \Vert \mathbf{u}^f - \mathbf{u}^p \Vert ^2 - \Vert \nabla \theta \Vert ^2 \\&\quad + R[(\theta ,w^{f*}+w^{p*}) + (w^f+w^p,\theta ^{*})] \\&\quad - \gamma _1 {\mathcal {T}}(\mathbf{k} \times \mathbf{u}^f, \mathbf{u}^{f*}) - \eta \gamma _1 {\mathcal {T}}(\mathbf{k} \times \mathbf{u}^p,\mathbf{u}^{p*}). \end{aligned} \end{aligned}$$

(9)

where $\mathbf{M}^f=(\mathbf{K}^f)^{-1}$ and $\mathbf{M}^p=(\mathbf{K}^p)^{-1}$, while $(\cdot , \cdot )$ and $\Vert \cdot \Vert$ are inner product and norm on the Hilbert space $L^2(V)$, respectively. Setting $\sigma =\sigma _r+i \sigma _i$, the imaginary part of equation (9) is

$$\begin{aligned} \sigma _i \Vert \theta \Vert ^2 = - \gamma _1 {\mathcal {T}}(\mathbf{k} \times \mathbf{u}^f, \mathbf{u}^{f*}) - \eta \gamma _1 {\mathcal {T}}(\mathbf{k} \times \mathbf{u}^p, \mathbf{u}^{p*}). \end{aligned}$$

(10)

Applying the same procedure to the complex conjugate of (8), multiplying by $\mathbf{u}^f, \mathbf{u}^p, \theta$ one gets

$$\begin{aligned} \sigma _i \Vert \theta \Vert ^2 = - \gamma _1 {\mathcal {T}}(\mathbf{k} \times \mathbf{u}^{f*}, \mathbf{u}^f) - \eta \gamma _1 {\mathcal {T}}(\mathbf{k} \times \mathbf{u}^{p*}, \mathbf{u}^p). \end{aligned}$$

(11)

Adding (10) and (11), one obtains

$$\begin{aligned} 2 \sigma _i \Vert \theta \Vert ^2 = 0 \ \Rightarrow \ \sigma _i=0 \ \Rightarrow \ \sigma \in {\mathbb {R}}\end{aligned}$$

and hence the strong version of the principle of exchange of stabilities holds: if the convection sets in, it arises necessary via a stationary motion (steady convection).

To determine the instability threshold for the onset of stationary convection, let us consider system (8) with $\sigma =0$, i.e.,

$$\begin{aligned} \left\{ \begin{aligned}&\frac{\gamma _1}{k_1} u^f+u^f-u^p = - \pi ^f_x + \gamma _1 {\mathcal {T}}v^f, \\&\frac{\gamma _1}{k_2} v^f+v^f-v^p = - \pi ^f_y - \gamma _1 {\mathcal {T}}u^f, \\&\gamma _1 w^f+w^f-w^p = - \pi ^f_z +R \theta , \\&\frac{\gamma _2}{h_1} u^p+u^p-u^f = - \pi ^p_x + \eta \gamma _1 {\mathcal {T}}v^p, \\&\frac{\gamma _2}{h_2} v^p+v^p-v^f = - \pi ^p_y - \eta \gamma _1 {\mathcal {T}}u^p, \\&\gamma _2 w^p+w^p-w^f = - \pi ^p_z +R \theta , \\&u^f_x+v^f_y+w^f_z = 0, \\&u^p_x+v^p_y+w^p_z = 0, \\&Rw^f+Rw^p+\Delta \theta =0. \end{aligned} \right. \end{aligned}$$

(12)

From system (1), one obtains

$$\begin{aligned} \left\{ \begin{aligned} u^f =&\frac{1}{D} \Bigl \{ H \pi ^p_x + \gamma _1 {\mathcal {T}}(L+\eta ) \pi ^f_y + [L(b+1)-(c+1)] \pi ^f_x + \gamma _1 {\mathcal {T}}N \pi ^p_y \Bigr \}, \\ u^p =&\frac{1}{D} \Bigl \{ H \pi ^f_x + \gamma _1 {\mathcal {T}}(\eta K+1) \pi ^p_y + [K(d+1)-(a+1)] \pi ^p_x + \gamma _1 {\mathcal {T}}M \pi ^f_y \Bigr \}, \\ v^f =&\frac{1}{D} \Bigl \{ G \pi ^p_y - \gamma _1 {\mathcal {T}}(L+\eta ) \pi ^f_x + [L(a+1)-(d+1)] \pi ^f_y -\gamma _1 {\mathcal {T}}M \pi ^p_x \Bigr \}, \\ v^p =&\frac{1}{D} \Bigl \{ G \pi ^f_y - \gamma _1 {\mathcal {T}}(\eta K +1) \pi ^p_x + [K(c+1)-(b+1)] \pi ^p_y - \gamma _1 {\mathcal {T}}N \pi ^f_x \Bigr \}, \end{aligned} \right. \end{aligned}$$

(13)

where A is the coefficients matrix of system $(12)_{1,2}$–$(12)_{4,5}$ and

$$\begin{aligned} \begin{aligned}a&=\frac{\gamma _1}{k_1}, \ b=\frac{\gamma _1}{k_2}, \ c=\frac{\gamma _2}{h_1}, \ d=\frac{\gamma _2}{h_2},\\ H&= (b+1)(d+1)-1-\eta (\gamma _1 {\mathcal {T}})^2, \\ G&= (a+1)(c+1)-1-\eta (\gamma _1 {\mathcal {T}})^2, \\ K&= (a+1)(b+1)+(\gamma _1 {\mathcal {T}})^2, \\ L&= (d+1)(c+1)+\eta ^2 (\gamma _1 {\mathcal {T}})^2, \\ M&= (d+1) + \eta (a+1), \\ N&= (c+1) + \eta (b+1), \\ D&= \text {det}(A)\,. \end{aligned} \end{aligned}$$

Hence, differentiating equations (13) with respect to z and equations $(12)_3$–$(12)_6$ with respect to x and y, by virtue of the incompressible conditions $(12)_{7}$–$(12)_{8}$, differentiating with respect to z, i.e.,

$$\begin{aligned} \begin{aligned} u^f_{xz}+v^f_{yz}&= -w^f_{zz}, \\ u^p_{xz}+v^p_{yz}&= -w^p_{zz}, \end{aligned} \end{aligned}$$

(14)

one obtains the following system

$$\begin{aligned} {\left\{ \begin{array}{ll} a_1 w^f_{xx} + a_2 w^f_{yy} + D w^f_{zz} + a_3 w^p_{xx} + a_4 w^p_{yy} + c_1 w^f_{xy} \\ \quad - c_1 \hat{\gamma _2} w^p_{xy} + R a_5 \theta _{xx} + R a_6 \theta _{yy} + R c_1 \theta _{xy} = 0, \\ \\ b_1 w^f_{xx} + b_2 w^f_{yy} + b_3 w^p_{xx} + b_4 w^p_{yy} + D w^p_{zz} + c_1 \hat{\gamma _1} w^f_{xy} \\ \quad - c_1 w^p_{xy} + R b_5 \theta _{xx} + R b_6 \theta _{yy} - R c_1 \theta _{xy} = 0, \\ \\ R w^f + R w^p + \Delta \theta =0. \end{array}\right. } \end{aligned}$$

(15)

where $\hat{\gamma _s}=\gamma _s+1$, for $s=1,2$ and

$$\begin{aligned} \begin{aligned}&a_1=H-\hat{\gamma _1}[L(b+1)-(c+1)], \\&a_2=G-\hat{\gamma _1}[L(a+1)-(d+1)], \\&a_3=-\hat{\gamma _2}H+L(b+1)-(c+1), \\&a_4=-\hat{\gamma _2}G+L(a+1)-(d+1), \\&a_5=H+L(b+1)-(c+1), \\&a_6=G+L(a+1)-(d+1), \\&b_1=- \hat{\gamma _1}H+K(d+1)-(a+1), \\&b_2=-\hat{\gamma _1}G+K(c+1)-(b+1), \\&b_3=H-\hat{\gamma _2}[K(d+1)-(a+1)], \\&b_4=G-\hat{\gamma _2}[K(c+1)-(b+1)], \\&b_5=H+K(d+1)-(a+1), \\&b_6=G+K(c+1)-(b+1), \\&c_1=\gamma _1 {\mathcal {T}}(N-M). \end{aligned} \end{aligned}$$

By virtue of periodicity of perturbation fields in the horizontal directions x and y, taking into account the boundary conditions (5), since $\{ \sin (n \pi z) \}_{n \in {\mathcal {N}}}$ is a complete orthogonal system for $L^2([0,1])$, employing normal modes solutions [7]:

$$\begin{aligned} \begin{aligned} w^f&= W^f_0 \sin (n \pi z) e^{i(lx+my)}, \\ w^p&= W^p_0 \sin (n \pi z) e^{i(lx+my)}, \\ \theta&= \Theta _0 \sin (n \pi z) e^{i(lx+my)}, \end{aligned} \end{aligned}$$

(16)

from (15) one obtains

$$\begin{aligned} {\left\{ \begin{array}{ll} h_{11} W^f_0 +h_{12} W^p_0 + R h_{13} \Theta _0 = 0, \\ h_{21} W^f_0 +h_{22} W^p_0 + R h_{23} \Theta _0 = 0, \\ R W^f_0 + R W^p_0 - \Lambda _n \Theta _0 = 0. \end{array}\right. } \end{aligned}$$

(17)

where $\Lambda _n=n^2 \pi ^2 + l^2 + m^2$ and

$$\begin{aligned} \begin{aligned} h_{11}&= a_1 l^2 + a_2 m^2 +D n^2 \pi ^2 + c_1 lm, \\ h_{12}&= a_3 l^2 +a_4 m^2 -c_1 \hat{\gamma _2} lm, \\ h_{13}&= a_5 l^2 + a_6 m^2 +c_1 lm, \\ h_{21}&= b_1 l^2 +b_2 m^2 +c_1 \hat{\gamma _1} lm, \\ h_{22}&= b_3 l^2 + b_4 m^2 +D n^2 \pi ^2 -c_1 lm, \\ h_{23}&= b_5 l^2 +b_6 m^2 -c_1 lm. \end{aligned} \end{aligned}$$

Finally, requiring a zero determinant for (17), the linear instability threshold for the onset of stationary convection is:

$$\begin{aligned} R^2_L=\min _{n,l,m} \dfrac{\Lambda _n (h_{11} h_{22} - h_{12} h_{21})}{h_{12} h_{23} - h_{13} h_{22} - h_{11} h_{23} + h_{21} h_{13}} \end{aligned}$$

(18)

where both numerator and denominator of the right hand side of (18) are positive. The minimization (18) is numerically analyzed in Sect. 5.

Let us point out that

• (i) assuming $h_1=h_2=h$, $k_1=k_2=k$ (horizontally isotropic case) and for ${\mathcal {T}}\rightarrow 0$, i.e., in the absence of rotation, we get

  $$\begin{aligned} R_L^2=\min _{n,a^2} \dfrac{{\hat{\Gamma }} \Lambda _n}{a^2} \dfrac{a^4 \Gamma {\hat{\Gamma }}^{-1} + n^4 \pi ^4 + a^2 n^2 \pi ^2 {\hat{K}} {\hat{\Gamma }}^{-1} }{a^2 (4+\gamma _1+\gamma _2) + n^2 \pi ^2 (\frac{\gamma _1}{k} + \frac{\gamma _2}{h} +4)} \end{aligned}$$

  where $a^2=l^2+m^2$ and

  $$\begin{aligned} \begin{aligned} {\hat{K}}&= \Bigl [ \gamma _1+\gamma _2+\frac{\gamma _1}{k}+\frac{\gamma _2}{h} + \gamma _1 \gamma _2 \Big ( \frac{1}{k} + \frac{1}{h} \Big ) \Bigr ], \\ \Gamma&= \gamma _1 \gamma _2 + \gamma _1 + \gamma _2, \\ {\hat{\Gamma }}&= \frac{\gamma _1}{k} + \frac{\gamma _2}{h} + \frac{\gamma _1}{k} \frac{\gamma _2}{h}, \end{aligned} \end{aligned}$$

  that coincides with the instability threshold found in [29];

• (ii) the case of a non-rotating layer of isotropic bidisperse porous medium (assuming $h_s=k_s=1$ for $s=1,2$ and as ${\mathcal {T}}\rightarrow 0$) leads to

  $$\begin{aligned} R^2_L= \min _{n,a^2} \dfrac{\Lambda ^2}{a^2} \dfrac{\gamma _1 \gamma _2 + \gamma _1 + \gamma _2}{\gamma _1 + \gamma _2 +4} \end{aligned}$$

  that is the same threshold found in [8].

4 Nonlinear stability

In order to study the influence of rotation on the nonlinear stability of the conduction solution, since the Coriolis terms in momentum equations are antisymmetric, instead of applying the standard energy method, let us apply the differential constraint approach (see [5, 11, 24]).

To this end, let us set

$$\begin{aligned} \begin{aligned} E(t)&= \frac{1}{2} \Vert \theta \Vert ^2, \\ I(t)&= (w^f+w^p,\theta ), \\ D(t)&= \Vert \nabla \theta \Vert ^2, \end{aligned} \end{aligned}$$

(19)

and by virtue of $(4)_5$, one obtains

$$\begin{aligned} \dfrac{dE}{dt} = \Bigl ( R\frac{I}{D} -1 \Bigr )D\,. \end{aligned}$$

(20)

Setting

$$\begin{aligned} \frac{1}{R_E} = \max _{{\mathcal {H}}^{*}} \frac{I}{D} \end{aligned}$$

(21)

with

$$\begin{aligned} \begin{aligned} {\mathcal {H}}^{*} = \{&(w^f,w^p,\theta ) \in (H^1)^3 | w^f=w^p=\theta =0 \ \text {on} \ z=0,1; \\&\text {periodic in} \ x,y \ \text {with periods} \ 2 \pi /l, 2 \pi /m;D<\infty ; \\&\text {verifying} \ (15)_{1,2} \} \end{aligned} \end{aligned}$$

the space of kinematically admissible solutions. The variational problem (21) is equivalent to the following variational problem:

$$\begin{aligned} \frac{1}{R_E}= \max _{{\mathcal {H}}} \frac{I + \int \limits _{V} \lambda g_1 \ dV + \int \limits _{V} \psi g_2 \ dV }{D}, \end{aligned}$$

(22)

where $\lambda (\mathbf{x})$ and $\psi (\mathbf{x})$ are Lagrange multipliers and

$$\begin{aligned} \begin{aligned} g_1 \equiv&R^{-1}(a_1 w^f_{xx} +a_2 w^f_{yy} + D w^f_{zz} + a_3 w^p_{xx}+ a_4 w^p_{yy} + c_1 w^f_{xy} - c_1 \hat{\gamma _2} w^p_{xy})+ a_5 \theta _{xx} + a_6 \theta _{yy} + c_1 \theta _{xy}, \\ g_2 \equiv&R^{-1}(b_1 w^f_{xx} + b_2 w^f_{yy} + b_3 w^p_{xx} + b_4 w^p_{yy}+ D w^p_{zz} + c_1 \hat{\gamma _1} w^f_{xy} - c_1 w^p_{xy}) + b_5 \theta _{xx}+ b_6 \theta _{yy} - c_1 \theta _{xy}, \end{aligned} \end{aligned}$$

(23)

$$\begin{aligned} \begin{aligned} {\mathcal {H}} = \{&(w^f,w^p,\theta ) \in (H^1)^3 | w^f=w^p=\theta =0 \ \text {on} \ z=0,1; \\&\text {periodic in} \ x,y \ \text {with periods} \ 2 \pi /l, 2 \pi /m;D< \infty \}. \end{aligned} \end{aligned}$$

By virtue of Poincaré inequality, since

$$\begin{aligned} D \ge \pi ^2 \Vert \theta \Vert ^2, \end{aligned}$$

from (20) one obtains that condition $R<R_E$ guarantees the global nonlinear stability of the conduction solution with respect to the $L^2$-norm, according to the following inequality

$$\begin{aligned} E(t) \le E(0) \exp \Bigl [ \frac{R-R_E}{R_E} \ t \Bigr ]. \end{aligned}$$

Remark 1

Multiplying $(8)_1$ by $\mathbf{u}^f$, $(8)_2$ by $\mathbf{u}^p$, integrating over the period cell V and adding the resulting equations, one finds

$$\begin{aligned} \begin{aligned}&\gamma _1 \int \limits _{V} \Bigl [ \frac{1}{k_1} (u^{f})^2 + \frac{1}{k_2} (v^{f})^2 + (w^{f})^2 \Bigr ] \ dV + \gamma _2 \int \limits _{V} \Bigl [ \frac{1}{h_1} (u^{p})^2 + \frac{1}{h_2} (v^{p})^2 + (w^{p})^2 \Bigr ] \ dV \\&\qquad + \Vert \mathbf{u}^f - \mathbf{u}^p \Vert ^2 = R (\theta , w^f+w^p). \end{aligned} \end{aligned}$$

(24)

Setting ${\hat{k}}=\max (k_1,k_2,1)$ and ${\hat{h}}=\max (h_1,h_2,1)$ and using the generalized Cauchy inequality on the right-hand side of (24), one obtains

$$\begin{aligned} \frac{\gamma _1}{{\hat{k}}} \Vert \mathbf{u}^f \Vert ^2 + \frac{\gamma _2}{{\hat{h}}} \Vert \mathbf{u}^p \Vert ^2 \le R^2 \Big ( \frac{{\hat{k}}}{\gamma _1} + \frac{{\hat{h}}}{\gamma _2} \Big ) \Vert \theta \Vert ^2 \end{aligned}$$

(25)

and hence condition $R<R_E$ guarantees that $\Vert \mathbf{u}^f \Vert ^2 \rightarrow 0$ and $\Vert \mathbf{u}^p \Vert ^2 \rightarrow 0$ as $t \rightarrow \infty$, too.

In order to solve the variational problem (22), let us consider the associated Euler–Lagrange equations:

$$\begin{aligned} {\left\{ \begin{array}{ll} R_E (w^f+w^p) + R_E (a_5 \lambda _{xx} + a_6 \lambda _{yy} + c_1 \lambda _{xy} + b_5 \psi _{xx} + b_6 \psi _{yy} - c_1 \psi _{xy} ) +2 \Delta \theta =0, \\ R_E \theta + a_1 \lambda _{xx} + a_2 \lambda _{yy} + D \lambda _{zz} + c_1 \lambda _{xy} + b_1 \psi _{xx} + b_2 \psi _{yy} + c_1 \hat{\gamma _1} \psi _{xy} =0, \\ R_E \theta + a_3 \lambda _{xx} + a_4 \lambda _{yy} - c_1 \hat{\gamma _2} \lambda _{xy} + b_3 \psi _{xx} +b_4 \psi _{yy} + D \psi _{zz} - c_1 \psi _{xy} =0, \\ a_1 w^f_{xx} + a_2 w^f_{yy} + D w^f_{zz} + a_3 w^p_{xx} + a_4 w^p_{yy} + c_1 w^f_{xy} - c_1 \hat{\gamma _2} w^p_{xy} + R_E( a_5 \theta _{xx} + a_6 \theta _{yy} + c_1 \theta _{xy}) =0, \\ b_1 w^f_{xx} + b_2 w^f_{yy} + b_3 w^p_{xx} + b_4 w^p_{yy} + D w^p_{zz} + c_1 \hat{\gamma _1} w^f_{xy} - c_1 w^p_{xy} + R_E (b_5 \theta _{xx} + b_6 \theta _{yy} - c_1 \theta _{xy}) =0. \end{array}\right. } \end{aligned}$$

(26)

Defining the operators

$$\begin{aligned} \begin{aligned} \Delta _f&\equiv a_1 \partial _{xx}+a_2 \partial _{yy}+D \partial _{zz}, \\ \Delta _p^*&\equiv a_3 \partial _{xx}+a_4 \partial _{yy}, \\ \Delta _f^*&\equiv b_1 \partial _{xx}+b_2 \partial _{yy}, \\ \Delta _p&\equiv b_3 \partial _{xx}+b_4 \partial _{yy}+D \partial _{zz}, \\ {\mathcal {L}}_1&\equiv a_5 \partial _{xx}+a_6 \partial _{yy}, \\ {\mathcal {L}}_2&\equiv b_5 \partial _{xx}+b_6 \partial _{yy}. \end{aligned} \end{aligned}$$

and taking $2 \Delta$ of $(26)_{2,3,4,5}$, the Euler–Lagrange equations become

$$\begin{aligned} {\left\{ \begin{array}{ll} R_E(w^f+w^p)+R_E({\mathcal {L}}_1 \lambda + c_1 \lambda _{xy} + {\mathcal {L}}_2 \psi - c_1 \psi _{xy}) = -2 \Delta \theta , \\ 2 \Delta \Delta _f \lambda + 2 c_1 \Delta \lambda _{xy} + 2 \Delta \Delta _f^* \psi + 2 c_1 \hat{\gamma _1} \Delta \psi _{xy} = -2 R_E \Delta \theta , \\ 2 \Delta \Delta _p^* \lambda - 2 c_1 \hat{\gamma _2} \Delta \lambda _{xy} + 2 \Delta \Delta _p \psi - 2 c_1 \Delta \psi _{xy} = -2 R_E \Delta \theta , \\ 2 \Delta \Delta _f w^f + 2 c_1 \Delta w^f_{xy} + 2 \Delta \Delta _p^* w^p - 2 c_1 \hat{\gamma _2} \Delta w^p_{xy} \\ \quad + 2 R_E {\mathcal {L}}_1 \Delta \theta + 2 R_E c_1 \Delta \theta _{xy}=0, \\ 2 \Delta \Delta _f^* w^f + 2 c_1 \hat{\gamma _1} \Delta w^f_{xy} + 2 \Delta \Delta _p w^p - 2 c_1 \Delta w^p_{xy} \\ \quad + 2 R_E {\mathcal {L}}_2 \Delta \theta - 2 R_E c_1 \Delta \theta _{xy}=0 \end{array}\right. } \end{aligned}$$

(27)

Eliminating variable $\theta$ and setting

$$\begin{aligned} \begin{aligned} {\mathcal {M}}_1&\equiv 2 \Delta \Delta _f + 2 c_1 \Delta \partial _{xy} -R_E^2 {\mathcal {L}}_1 - R_E^2 c_1 \partial _{xy}, \\ {\mathcal {M}}_2&\equiv 2 \Delta \Delta _f^* + 2 c_1 \hat{\gamma _1} \Delta \partial _{xy} -R_E^2 {\mathcal {L}}_2 + R_E^2 c_1 \partial _{xy}, \\ {\mathcal {M}}_3&\equiv 2 \Delta \Delta _p^* - 2 c_1 \hat{\gamma _2} \Delta \partial _{xy} -R_E^2 {\mathcal {L}}_1 - R_E^2 c_1 \partial _{xy}, \\ {\mathcal {M}}_4&\equiv 2 \Delta \Delta _p - 2 c_1 \Delta \partial _{xy} -R_E^2 {\mathcal {L}}_2 + R_E^2 c_1 \partial _{xy}, \\ {\mathcal {N}}_1&\equiv -R_E( {\mathcal {L}}_1+c_1 \partial _{xy})^2, \\ {\mathcal {N}}_2&\equiv -R_E^2( {\mathcal {L}}_1+c_1 \partial _{xy})({\mathcal {L}}_2-c_1 \partial _{xy}), \\ {\mathcal {N}}_3&\equiv -R_E^2({\mathcal {L}}_2-c_1 \partial _{xy})^2, \end{aligned} \end{aligned}$$

one obtains

$$\begin{aligned} {\left\{ \begin{array}{ll} - R_E^2 w^f - R_E^2 w^p +{\mathcal {M}}_1 \lambda +{\mathcal {M}}_2 \psi =0, \\ - R_E^2 w^f - R_E^2 w^p +{\mathcal {M}}_3 \lambda +{\mathcal {M}}_4 \psi =0, \\ {\mathcal {M}}_1 w^f + {\mathcal {M}}_3 w^p + {\mathcal {N}}_1 \lambda + {\mathcal {N}}_2 \psi =0, \\ {\mathcal {M}}_2 w^f + {\mathcal {M}}_4 w^p + {\mathcal {N}}_2 \lambda + {\mathcal {N}}_3 \psi =0. \end{array}\right. } \end{aligned}$$

(28)

By employing normal modes

$$\begin{aligned} \begin{aligned} w^f&= W^f_0 \sin (n \pi z) e^{i(lx+my)}, \\ w^p&= W^p_0 \sin (n \pi z) e^{i(lx+my)}, \end{aligned} \end{aligned}$$

(29)

and choosing [6, 11]

$$\begin{aligned} \begin{aligned} \lambda&= \lambda _0 \sin (n \pi z) e^{i(lx+my)}, \\ \psi&= \psi _0 \sin (n \pi z) e^{i(lx+my)}, \end{aligned} \end{aligned}$$

(30)

from (28) one obtains

$$\begin{aligned} {\left\{ \begin{array}{ll} -R_E^2 W^f_0 - R_E^2 W^p_0 +(2 \Lambda _n h_{11} + R_E^2 h_{13}) \lambda _0 +(2 \Lambda _n h_{21}+R_E^2 h_{23}) \psi _0=0, \\ -R_E^2 W^f_0 - R_E^2 W^p_0 +(2 \Lambda _n h_{12} + R_E^2 h_{13}) \lambda _0 +(2 \Lambda _n h_{22}+R_E^2 h_{23}) \psi _0=0, \\ (2 \Lambda _n h_{11}+R_E^2 h_{13}) W^f_0 + (2 \Lambda _n h_{12} +R_E^2 h_{13}) W^p_0 -R_E^2 h_{13}^2 \lambda _0 -R_E^2 h_{13} h_{23} \psi _0=0, \\ (2 \Lambda _n h_{21}+R_E^2 h_{23}) W^f_0 + (2 \Lambda _n h_{22} +R_E^2 h_{23}) W^p_0 -R_E^2 h_{13} h_{23} \lambda _0 -R_E^2 h_{23}^2 \psi _0=0. \end{array}\right. } \end{aligned}$$

(31)

Requiring a zero determinant for (31) we find

$$\begin{aligned} R_E^2=R_L^2, \end{aligned}$$

and hence the global nonlinear stability threshold and the linear instability threshold coincide and subcritical instabilities do not exist.

5 Numerical simulations

We now present numerical results to solve (18), in order to analyze the asymptotic behavior of $R_L^2$ with respect to ${\mathcal {T}}$, $h_i$, $k_i$, for $i=1,2$, i.e., to study the influence of rotation and anisotropic permeability on the onset of convection. As regards the physical parameters, in all numerical simulations, we have chosen a set of values analogous to those ones fixed in [27], in order to compare our results with those ones obtained in [27], to stress the influence of rotation and anisotropy on the onset of convection.

In all the computations, we have performed, the minimum of $R_L^2$ with respect to n is attained at $n=1$. Each of the following tables and figures show the stabilizing effect of rotation on the onset of convection.

Table 1 Critical $R^2_L,l,m$ for increasing ${\mathcal {T}}$ and $h_1=10,h_2=1,k_1=0.1,k_2=1,\eta =0.2,\gamma _1=10,\gamma _2=50$

Table 2 Critical $R^2_L,l,m$ for increasing ${\mathcal {T}}$ and $h_1=0.1,h_2=1,k_1=10,k_2=1,\eta =0.2,\gamma _1=10,\gamma _2=50$

Fig. 1

a Critical Rayleigh number $R^2_L$ as function of the Taylor number ${\mathcal {T}}$ for $h_1=10,h_2=1,k_1=0.1,k_2=1,\eta =0.2,\gamma _1=10,\gamma _2=50$. b Critical Rayleigh number $R^2_L$ as function of the Taylor number ${\mathcal {T}}$ for $h_1=0.1,h_2=1,k_1=10,k_2=1,\eta =0.2,\gamma _1=10,\gamma _2=50$

Table 1 shows that for large values of the Taylor number ${\mathcal {T}}$ and when $h_1>>k_1$, m becomes zero, this means that, as rotation increases, the convection cells become rolls with the axis in the y-direction. Table 2 shows a transition from convection patterns as rolls along y-axis ($m=0$ for very small ${\mathcal {T}}$) to convection patterns as rolls along x-axis ($l=0$), as the rotation increases and for $h_1<<k_1$. For these physical values, the asymptotic behavior of $R_L^2$ with respect to ${\mathcal {T}}$ is shown in Fig. 1. We can also observe that, as ${\mathcal {T}}$ increases, $R_L^2$ increases more slowly when $h_1<<k_1$ then $h_1>>k_1$.

Let us point out that bi-dimensional convection cells (rolls along x-axis for $l=0$ and rolls along y-axis for $m=0$) were already found in [27] as an effect of anisotropic macropermeability and micropermeability in absence of rotation.

Table 3 (a): critical $R^2_L,l,m$ for increasing ${\mathcal {T}}$. (b): critical Rayleigh number $R^2_L$ as function of the Taylor number ${\mathcal {T}}$. For $h_1=1,h_2=0.1,k_1=1,k_2=10,\eta =0.2,\gamma _1=2,\gamma _2=0.2$

From Table 3, we numerically find out that for parameters $\{ h_1=1,h_2=0.1,k_1=1,k_2=10,\eta =0.2,\gamma _1=2,\gamma _2=0.2 \}$ the critical value of m is mainly zero, except for very small values of the Taylor number ${\mathcal {T}}\in [0,3)$, for which l and m are both nonzero, i.e., for very little rotation of the layer, three-dimensional convection cells are expected.

As a matter of fact, the wavelengths in the x and y directions are ${\hat{x}}=\frac{2 \pi }{l}$ and ${\hat{y}}=\frac{2 \pi }{m}$. The condition ${\hat{y}} / {\hat{x}}=0$ implies $l=0$, this means that the convective fluid motion occurs in the y and z directions (the solution is a function of y and z), i.e., the convection cells are rolls in the x-direction. Instead, the condition ${\hat{y}} / {\hat{x}} \rightarrow \infty$ implies $m=0$ and the convective fluid motion occurs in the x and z directions, so the cells are rolls in the y-direction [27].

Table 4 Critical $R^2_L,l,m$ for quoted values of $h_1$ (a) and for quoted values of $h_2$ (b). Table a: $h_2=0.9,k_1=0.2,k_2=1.1,\eta =0.2,\gamma _1=0.9,\gamma _2=1.8,{\mathcal {T}}=10$. Table b: $h_1=3.3,k_1=0.2,k_2=1.1,\eta =0.2,\gamma _1=0.9,\gamma _2=1.8,{\mathcal {T}}=10$

Table 5 Critical $R^2_L,l,m$ for quoted values of $k_1$ (a) and for quoted values of $k_2$ (b). Table a: $h_1=3.3,h_2=0.9,k_2=1.1,\eta =0.2,\gamma _1=0.9,\gamma _2=1.8,{\mathcal {T}}=10$. Table b: $h_1=3.3,h_2=0.9,k_1=0.2,\eta =0.2,\gamma _1=0.9,\gamma _2=1.8,{\mathcal {T}}=10$

Tables 4 and 5 exhibit the influence of anisotropy parameters for both macropores and micropores on the onset of convection, and the values of $h_1,h_2,k_1,k_2$ are fixed such that the permeability ratios in the macropores and micropores are different, in particular we set $\{ h_1=3.3,h_2=0.9,k_1=0.2,k_2=1.1 \}$ (see [27]) and we vary $h_s,k_s$ for $s=1,2$ in turn to see how each parameter affects the Rayleigh number. As in [27], we numerically find out a very complex relationship between the macro and micro permeability parameters and the critical Rayleigh and wave numbers. For increasing $h_1,h_2,k_1,k_2$, we can see a similar trend, i.e., $R_L^2$ increases up to a maximum before decreasing. From 4(a) and from 5(a), we can see a first transition from rolls along x-axis to three-dimensional cells and then another transition to rolls along y-axis, while 4(b) and 5(b) displays a mirror behavior with respect to 4(a) and 5(a) , respectively.

In Fig. 2 the critical Rayleigh number $R_L^2$ is represented as function of the Taylor number ${\mathcal {T}}$ for $h_1=0.1,1,5,10$ and the others parameters are fixed as $h_2=0.9,k_1=0.2,k_2=1.1,\eta =0.2,\gamma _1=0.9,\gamma _2=1.8$, with the aim to graphically analyze the values shown in Table 4a.

Fig. 2

Critical Rayleigh number $R_L^2$ as function of the Taylor number ${\mathcal {T}}$ for $h_1=0.1,1,5,10$ and $h_2=0.9,k_1=0.2,k_2=1.1,\eta =0.2,\gamma _1=0.9,\gamma _2=1.8$

6 Conclusions

The onset of thermal convection in a horizontal layer of anisotropic BDPM, uniformly rotating about a vertical axis and uniformly heated from below, has been analyzed, according to Darcy’s law in both micropores and macropores. In particular, it has been proved that:

• the strong version of the principle of exchange of stabilities holds, and hence, when the convection arises, it sets in through a stationary motion;

• the linear instability threshold and the global nonlinear stability threshold in the $L^2-$norm coincide: this is an optimal result since the stability threshold furnishes a necessary and sufficient condition to guarantee the global (i.e., for all initial data) nonlinear stability.

Moreover, it has been numerically analyzed the influence of the rotation and the influence of the anisotropy on the onset of convection.