Instability of Vertical Throughflows in Bidisperse Porous Media

Abstract

In this paper, the instability of a vertical fluid motion, or throughflow, is investigated in a horizontal bidisperse porous layer that is uniformly heated from below. By means of the order-1 Galerkin approximation method, the critical Darcy–Rayleigh number for the onset of steady instability is determined in closed form. The coincidence between the linear instability threshold and the global nonlinear stability threshold, in the energy norm, is shown.

1. Introduction

In recent years, thermal convection in bidisperse porous media has attracted the interest of many researchers due to the numerous applications in which they are used. Bidisperse porous media (BDPM) were defined in [1] as media composed of clusters of large particles that are agglomerations of small particles. In this way, BDPM can be regarded as regular porous media containing fissures or cracks. There are macropores between the clusters and micropores within them. The macropores are referred to as f-phase, while the remainder of the structure is referred to as p-phase. Artificial BDPM can be constructed in order to inhibit or promote the onset of convection.

The onset of thermal convection in BDPM is widely studied. In [2,3,4,5], a theoretical fundation for fluid motion in BDPM can be found. A mathematical model incorporating two velocities, two pressures, and two temperatures for macro and micro phases was introduced.

In [6], the local thermal equilibrium between the macro and micropores was assumed, and hence a mathematical model with independent velocity and pressure has been introduced. This model has been widely used to investigate the onset of thermal convection in BDPM incorporating various effects. In particular, the effect of the uniform rotation about a vertical axis has been investigated in [7,8,9,10,11] for isotropic and anisotropic BDPM with or without the inertia term; the double-diffusive thermal convection has been analyzed in [12,13,14,15,16].

When the rest state loses its stability and steady convection occurs, a secondary steady motion, or throughflow, is observable. The stability of fluid motion in porous media finds relevant applications in industrial processes, geophysics, and astronomy and has been analyzed in many papers (see, for example [17,18,19,20,21,22,23,24,25,26] and the references therein).

In [17,18,19,20,21,22,23,24,25,26], the stability of vertical constant throughflows has been performed incorporating various effects, such as viscous dissipation [17]; chemicals dissolved in the fluid [19,22,25,26]; and an external magnetic field acting on an electrically conducting fluid [23]. The stability of non-constant throughflows was performed in [27].

In this paper, stability analysis of a vertical constant throughflow, saturating a bidisperse porous medium heated from below, is performed. The paper is organized as follows. Section 2 is devoted to the mathematical model and to the determination of a constant steady state solution that is different from the rest state (i.e., different from the conduction solution). Section 3 deals with the linear stability analysis and, by using an order-1 Galerkin approximation method, the critical Darcy–Rayleigh number for the onset of steady instability is determined in closed form. Nonlinear stability analysis in the energy norm is investigated in Section 4 showing the absence of subcritical instabilities. The paper ends with a final Section summarizing the results obtained.

2. Mathematical Model and Preliminaries

Let $\mathtt{L}$ be a horizontal layer of depth d filled by a bidisperse porous medium. Assume that the layer is uniformly heated from below and is filled by a Newtonian, homogeneous, incompressible fluid moving vertically. Introducing a reference frame $Oxyz$ with fundamental unit vectors $\mathbf{i},\mathbf{j},\mathbf{k}$ ($\mathbf{k}$ pointing vertically upward) so that $\mathtt{L}={\mathbb{R}}^2\times \left[0,d\right]$, the equations governing the fluid motion in $\mathtt{L}$, on assuming the Oberbeck–Boussinesq approximation, are [6]:

$$\left\{\begin{array}{c}{\displaystyle -\frac{\mu }{k_f}{\mathbf{v}}^f-\zeta ({\mathbf{v}}^f-{\mathbf{v}}^p)-\nabla p^f+{\varrho }_F\alpha gT\mathbf{k}=\mathbf{0},}\hfill \\ {\displaystyle -\frac{\mu }{k_p}{\mathbf{v}}^p-\zeta ({\mathbf{v}}^p-{\mathbf{v}}^f)-\nabla p^p+{\varrho }_F\alpha gT\mathbf{k}=\mathbf{0},}\hfill \\ \nabla ·{\mathbf{v}}^f=0,\hfill \\ \nabla ·{\mathbf{v}}^p=0,\hfill \\ {\left(\varrho c\right)}_m{T}_{,t}+{\left(\varrho c\right)}_f({\mathbf{v}}^f+{\mathbf{v}}^p)·\nabla T=k_m\mathsf{\Delta }T,\hfill \end{array}\right.$$

(1)

where $\mathbf{x}=(x,y,z)$, subscript ”${}_{,}t$” denotes the time derivative, ${\mathbf{v}}^s$ is the seepage velocity, $p^s$ is the pressure, T is the temperature, $\varrho$ is the density, $\zeta$ is the interaction coefficient between the f-phase and the p-phase, $\mathbf{g}=-g\mathbf{k}$ is the gravity, $\mu$ is the fluid viscosity, ${\varrho }_F$ is the reference constant density, $\alpha$ is the thermal expansion coefficient, c is the specific heat, $c_p$ is the specific heat at a constant pressure, ${\left(\varrho c\right)}_m=(1-\phi )(1-ϵ){\left(\varrho c\right)}_{\mathrm{sol}}+\phi {\left(\varrho c\right)}_f+ϵ(1-\phi ){\left(\varrho c\right)}_p$, and $k_m=(1-\phi )(1-ϵ)k_{\mathrm{sol}}+\phi k_f+ϵ(1-\phi )k_p$ is the thermal conductivity (the subscript ”$\mathrm{sol}$” refers to the solid skeleton, and $s=\{f,p\}$ for f-phase and p-phase, respectively).

The followng boundary conditions are applied to the system (1):

$$\begin{array}{c}T\equiv T_L,\mathrm{at}z=0;T\equiv T_U\mathrm{at}z=d,\hfill \end{array}$$

(2)

where $T_L>T_U$.

The problem (1)–(2) admits the stationary solution (vertical constant throughflow) $m_t=\{{\overline{\mathbf{v}}}^f,{\overline{\mathbf{v}}}^p,{\overline{p}}^f,{\overline{p}}^p,\overline{T}\}$:

$$\begin{array}{c}{\displaystyle {\overline{\mathbf{v}}}^f=Q^f\mathbf{k},{\overline{\mathbf{v}}}^p=Q^p\mathbf{k},\overline{T}=\frac{T_U-T_L{e}^{\overline{Q}d/k}+(T_L-T_U)e^{\overline{Q}z/k}}{1-e^{\overline{Q}d/k}},\overline{Q}=Q^f+Q^p,}\hfill \\ {\overline{p}}^f={\overline{p}}_0^f+{\varrho }_F\alpha g{\displaystyle \frac{T_L-T_U}{1-e^{\overline{Q}d/k}}}{\displaystyle \frac{k}{\overline{Q}}}(e^{\overline{Q}z/k}-1)+\left({\varrho }_F\alpha g{\displaystyle \frac{T_U-T_L}{1-e^{\overline{Q}d/k}}}{e}^{\overline{Q}d/k}-{\displaystyle \frac{\mu }{k_f}}{Q}^f-\zeta (Q^f-Q^p)\right)z,\hfill \\ {\overline{p}}^p={\overline{p}}_0^p+{\varrho }_F\alpha g{\displaystyle \frac{T_L-T_U}{1-e^{\overline{Q}d/k}}}{\displaystyle \frac{k}{\overline{Q}}}(e^{\overline{Q}z/k}-1)+\left({\varrho }_F\alpha g{\displaystyle \frac{T_U-T_L}{1-e^{\overline{Q}d/k}}}{e}^{\overline{Q}d/k}-{\displaystyle \frac{\mu }{k_f}}{Q}^p-\zeta (Q^p-Q^f)\right)z,\hfill \end{array}$$

with $\overline{Q}=Q^f+Q^p,k={\displaystyle \frac{k_m}{{\left(\varrho c\right)}_f}}$(thermal diffusivity). Setting

$${\mathbf{u}}^f={\mathbf{v}}^f-{\overline{\mathbf{v}}}^f,{\mathbf{u}}^p={\mathbf{v}}^p-{\overline{\mathbf{v}}}^p,{\mathsf{\Pi }}^f=p^f-{\overline{p}}^f,{\mathsf{\Pi }}^p=p^p-{\overline{p}}^p,\theta =T-\overline{T},$$

(3)

the system governing the evolution of the perturbation fields is:

$$\left\{\begin{array}{c}-{\displaystyle \frac{\mu }{k_f}}{\mathbf{u}}^f-\zeta ({\mathbf{u}}^f-{\mathbf{u}}^p)-\nabla {\mathsf{\Pi }}^f+{\varrho }_F\alpha g\theta \mathbf{k}=\mathbf{0},\hfill \\ -{\displaystyle \frac{\mu }{k_p}}{\mathbf{u}}^p-\zeta ({\mathbf{u}}^p-{\mathbf{u}}^f)-\nabla {\mathsf{\Pi }}^p+{\varrho }_F\alpha g\theta \mathbf{k}=\mathbf{0},\hfill \\ \nabla ·{\mathbf{u}}^f=0,\hfill \\ \nabla ·{\mathbf{u}}^p=0,\hfill \\ {\displaystyle \frac{{\left(\varrho c\right)}_m}{{\left(\varrho c\right)}_f}}{\theta }_{,t}+({\mathbf{u}}^f+{\mathbf{u}}^p)·\nabla \theta =-(w^f+w^p){\overline{T}}_{,z}-\overline{Q}{\theta }_{,z}+k\mathsf{\Delta }\theta ,\hfill \end{array}\right.$$

(4)

where ${\mathbf{u}}^f=(u^f,v^f,w^f),{\mathbf{u}}^p=(u^p,v^p,w^p)$, under the boundary conditions,

$${\mathbf{u}}^f·\mathbf{n}={\mathbf{u}}^p·\mathbf{n}=\theta =0\mathrm{at}z=0,d$$

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

Introducing the non-dimensional parameters,

$${\displaystyle {\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}},{\mathsf{\Pi }}^{s*}=\frac{{\mathsf{\Pi }}^s}{\tilde{P}},\mathrm{for}s=\{f,p\},}$$

$${\displaystyle {\gamma }_1=\frac{\mu }{k_f\zeta },{\gamma }_2=\frac{\mu }{k_p\zeta },}$$

where the scales are given by

$$\tilde{u}=\frac{k_m}{{\left(\varrho c\right)}_{f}d},\tilde{t}=\frac{d^2{\left(\varrho c\right)}_m}{k_m},\tilde{P}=\frac{\zeta k_m}{{\left(\varrho c\right)}_f},\tilde{T}=\sqrt{\frac{\beta k_m\zeta }{{\left(\varrho c\right)}_f{\varrho }_F\alpha g}},$$

the system (4) becomes (dropping all the asterisks):

$$\left\{\begin{array}{c}-{\gamma }_1{\mathbf{u}}^f-({\mathbf{u}}^f-{\mathbf{u}}^p)-\nabla {\mathsf{\Pi }}^f+\mathrm{Ra}\theta \mathbf{k}=\mathbf{0},\hfill \\ -{\gamma }_2{\mathbf{u}}^p-({\mathbf{u}}^p-{\mathbf{u}}^f)-\nabla {\mathsf{\Pi }}^p+\mathrm{Ra}\theta \mathbf{k}=\mathbf{0},\hfill \\ \nabla ·{\mathbf{u}}^f=0,\hfill \\ \nabla ·{\mathbf{u}}^p=0,\hfill \\ {\theta }_{,t}+({\mathbf{u}}^f+{\mathbf{u}}^p)·\nabla \theta =-\mathrm{Ra}\tilde{f}\left(z\right)(w^f+w^p)-\mathrm{Pe}{\theta }_{,z}+\mathsf{\Delta }\theta ,\hfill \end{array}\right.$$

(5)

where

$$\mathrm{Ra}=\sqrt{{\displaystyle \frac{\beta d^2{\left(\varrho c\right)}_f{\varrho }_F\alpha g}{k_m\zeta }}}$$

is the Darcy–Rayleigh thermal number and

$$\mathrm{Pe}={\displaystyle \frac{\overline{Q}d}{k}}(\mathrm{Péclet} \mathrm{number}),\tilde{f}\left(z\right)={\displaystyle \frac{\mathrm{Pe}{e}^{\mathrm{Pe}z}}{1-e^{\mathrm{Pe}}}}<0,\forall z\in \left[0,1\right],\forall \mathrm{Pe}.$$

The initial boundary conditions,

$$\begin{array}{c}{\mathbf{u}}^s(\mathbf{x},0)\equiv {\mathbf{u}}_0^s\left(\mathbf{x}\right),{\mathsf{\Pi }}^s(\mathbf{x},0)\equiv {\mathsf{\Pi }}_0\left(\mathbf{x}\right),\theta (\mathbf{x},0)\equiv {\theta }_0\left(\mathbf{x}\right),\hfill \\ w^f=w^p=\theta =0\mathrm{at}z=0,1,\hfill \end{array}$$

(6)

are appended to the system (5), with $\nabla ·{\mathbf{u}}_0^s=0$.

In the sequel, it is assumed the perturbation fields are periodic in the horizontal directions x and y of periods ${\displaystyle \frac{2\pi }{a_x}}$ and ${\displaystyle \frac{2\pi }{a_y}}$, respectively, and the periodicity cell is denoted by

$$V=\left[0,\frac{2\pi }{a_x}\right]\times \left[0,\frac{2\pi }{a_y}\right]\times [0,1].$$

3. Linear Instability

Let $\{{\widehat{\mathbf{u}}}^f,{\widehat{\mathbf{u}}}^p,{\widehat{\mathsf{\Pi }}}^f,{\widehat{\mathsf{\Pi }}}^p,\widehat{\theta }\}$ be the solution of the linearized version of the system (5), i.e.,

$$\left\{\begin{array}{c}-{\gamma }_1{\widehat{\mathbf{u}}}^f-({\widehat{\mathbf{u}}}^f-{\widehat{\mathbf{u}}}^p)-\nabla {\widehat{\mathsf{\Pi }}}^f+\mathrm{Ra}\widehat{\theta }\mathbf{k}=\mathbf{0},\hfill \\ -{\gamma }_2{\widehat{\mathbf{u}}}^p-({\widehat{\mathbf{u}}}^p-{\widehat{\mathbf{u}}}^f)-\nabla {\widehat{\mathsf{\Pi }}}^p+\mathrm{Ra}\widehat{\theta }\mathbf{k}=\mathbf{0},\hfill \\ \nabla ·{\widehat{\mathbf{u}}}^f=0,\hfill \\ \nabla ·{\widehat{\mathbf{u}}}^p=0,\hfill \\ {\widehat{\theta }}_{,t}=-\mathrm{Ra}\tilde{f}\left(z\right)({\widehat{w}}^f+{\widehat{w}}^p)-\mathrm{Pe}{\widehat{\theta }}_{,z}+\mathsf{\Delta }\widehat{\theta },\hfill \end{array}\right.$$

(7)

under the initial-boundary conditions

$$\begin{array}{c}{\widehat{\mathbf{u}}}^s(\mathbf{x},0)={\widehat{\mathbf{u}}}_0^s\left(\mathbf{x}\right),{\widehat{\mathsf{\Pi }}}^s(\mathbf{x},0)={\widehat{\mathsf{\Pi }}}_0\left(\mathbf{x}\right),\widehat{\theta }(\mathbf{x},0)={\widehat{\theta }}_0\left(\mathbf{x}\right),\hfill \\ {\widehat{w}}^f={\widehat{w}}^p=\widehat{\theta }=0\mathrm{at}z=0,1.\hfill \end{array}$$

(8)

The third components of the double curl of the first two equations of the system (7) along with the last equation of the system constitute a linear system governing the evolution of the three independent fields ${\widehat{w}}^f,{\widehat{w}}^p,\widehat{\theta }$:

$$\left\{\begin{array}{c}(1+{\gamma }_1)\mathsf{\Delta }{\widehat{w}}^f-\mathsf{\Delta }{\widehat{w}}^p-\mathrm{Ra}{\mathsf{\Delta }}_1\widehat{\theta }=0,\hfill \\ -\mathsf{\Delta }{\widehat{w}}^f+(1+{\gamma }_2)\mathsf{\Delta }{\widehat{w}}^p-\mathrm{Ra}{\mathsf{\Delta }}_1\widehat{\theta }=0,\hfill \\ {\widehat{\theta }}_{,t}=-\mathrm{Ra}\tilde{f}\left(z\right)({\widehat{w}}^f+{\widehat{w}}^p)-\mathrm{Pe}{\widehat{\theta }}_{,z}+\mathsf{\Delta }\widehat{\theta }.\hfill \end{array}\right.$$

(9)

Let us look for solutions of normal modes type

$$\left(\begin{array}{c}{\widehat{w}}^f(x,y,z,t)\\ {\widehat{w}}^p(x,y,z,t)\\ \widehat{\theta }(x,y,z,t)\end{array}\right)=\left(\begin{array}{c}{\overline{w}}^f\left(z\right)\\ {\overline{w}}^p\left(z\right)\\ \overline{\theta }\left(z\right)\end{array}\right)exp[-\sigma t+i(a_{x}x+a_{y}y)],$$

(10)

where $\sigma \in \mathbb{C}$. Setting

$$a^2=a_x^2+a_y^2,D\equiv {\displaystyle \frac{d}{dz}},$$

$\forall \phi \in \{{\widehat{w}}^f,{\widehat{w}}^p,\widehat{\theta }\}$, one has that

$${\mathsf{\Delta }}_1\phi =-a^2\phi ,\mathsf{\Delta }\phi =(D^2-a^2)\phi .$$

(11)

Then, the system (9) reads:

$$\left\{\begin{array}{c}(1+{\gamma }_1)(D^2-a^2){\overline{w}}^f-(D^2-a^2){\overline{w}}^p+\mathrm{Ra}{a}^2\overline{\theta }=0,\hfill \\ -(D^2-a^2){\overline{w}}^f+(1+{\gamma }_2)(D^2-a^2){\overline{w}}^p+\mathrm{Ra}{a}^2\overline{\theta }=0,\hfill \\ -\sigma \overline{\theta }=-\mathrm{Ra}\tilde{f}\left(z\right)({\overline{w}}^f+{\overline{w}}^p)-\mathrm{Pe}D\overline{\theta }+(D^2-a^2)\overline{\theta },\hfill \end{array}\right.$$

(12)

under the boundary conditions,

$${\overline{w}}^f={\overline{w}}^p=\overline{\theta }=0\mathrm{at}z=0,1.$$

(13)

To determine an approximation of the critical Darcy–Rayleigh number for the onset of steady instability, let us employ the order-1 Galerkin weighted residual method [24]. To this end, let us choose, as trial functions satisfying the boundary conditions (13),

$${\overline{w}}^f=C_{1}sin\pi z,{\overline{w}}^p=C_{2}sin\pi z,\overline{\theta }=C_{3}sin\pi z,$$

(14)

where $C_i$ are constants $i\in \{1,2,3\}$. Substituting the trial functions (14) in the system (12), one obtains three residuals. Making these residuals orthogonal to the trial functions over the range $0\le z\le 1$, one obtains the following system of three linear algebraic equations in the three unknowns, $C_1,C_2,C_3$:

$$\left\{\begin{array}{c}-(1+{\gamma }_1)(a^2+{\pi }^2)C_1+(a^2+{\pi }^2)C_2+\mathrm{Ra}{a}^2{C}_3=0,\hfill \\ (a^2+{\pi }^2)C_1-(1+{\gamma }_2)(a^2+{\pi }^2)C_2+\mathrm{Ra}{a}^2{C}_3=0,\hfill \\ {\displaystyle \frac{4{\pi }^2\mathrm{Ra}}{4{\pi }^2+{\mathrm{Pe}}^2}}{C}_1+{\displaystyle \frac{4{\pi }^2\mathrm{Ra}}{4{\pi }^2+{\mathrm{Pe}}^2}}{C}_2-\left[(a^2+{\pi }^2)-\sigma \right]C_3=0.\hfill \end{array}\right.$$

(15)

Requiring the vanishing of the determinant of the system (15), one has that

$${\mathrm{Ra}}^2={\displaystyle \frac{({\gamma }_1+{\gamma }_2+{\gamma }_1{\gamma }_2)(4{\pi }^2+{\mathrm{Pe}}^2)(a^2+{\pi }^2)(a^2+{\pi }^2-\sigma )}{4{\pi }^2{a}^2(4+{\gamma }_1+{\gamma }_2)}}.$$

(16)

From Equation (16), ${\mathrm{Ra}}^2$ is a real number if and only if $\sigma \in \mathbb{R}$, i.e., the principle of exchange of stability holds and instability can arise only via a steady motion. Hence, setting $\sigma =0$ in Equation (16), one obtains that the critical Darcy–Rayleigh number for the onset of steady instability is:

$${\displaystyle {\mathrm{Ra}}^2={\mathrm{Ra}}^{\left(s\right)}=\underset{a^2\in {\mathbb{R}}^{+}}{min}{\displaystyle \frac{({\gamma }_1+{\gamma }_2+{\gamma }_1{\gamma }_2)(4{\pi }^2+{\mathrm{Pe}}^2){(a^2+{\pi }^2)}^2}{4{\pi }^2{a}^2(4+{\gamma }_1+{\gamma }_2)}}.}$$

(17)

Simple calculations show that the minimum is reached for $a^2=a_s^2={\pi }^2$ and is given by

$${\mathrm{Ra}}^{\left(s\right)}={\displaystyle \frac{({\gamma }_1+{\gamma }_2+{\gamma }_1{\gamma }_2)(4{\pi }^2+{\mathrm{Pe}}^2)}{4+{\gamma }_1+{\gamma }_2}}.$$

(18)

Let us remark that ${\mathrm{Ra}}^{\left(s\right)}$ increases with $\mathrm{Pe}$, and

$${\displaystyle \underset{\mathrm{Pe}\to 0}{lim}{\mathrm{Ra}}^{\left(s\right)}={\displaystyle \frac{4{\pi }^2({\gamma }_1+{\gamma }_2+{\gamma }_1{\gamma }_2)}{4+{\gamma }_1+{\gamma }_2}},}$$

(19)

which is the critical Darcy–Rayleigh thermal number for the onset of steady convection found in [28].

4. Nonlinear Stability

Let us consider the nonlinear system,

$$\left\{\begin{array}{c}(1+{\gamma }_1)\mathsf{\Delta }{w}^f-\mathsf{\Delta }{w}^p-\mathrm{Ra}{\mathsf{\Delta }}_1\theta =0,\hfill \\ -\mathsf{\Delta }{w}^f+(1+{\gamma }_2)\mathsf{\Delta }{w}^p-\mathrm{Ra}{\mathsf{\Delta }}_1\theta =0,\hfill \\ {\theta }_{,t}+({\mathbf{u}}^f+{\mathbf{u}}^p)·\nabla \theta =-\mathrm{Ra}\tilde{f}\left(z\right)(w^f+w^p)-\mathrm{Pe}{\theta }_{,z}+\mathsf{\Delta }\theta ,\hfill \\ \nabla ·{\mathbf{u}}^f=0,\nabla ·{\mathbf{u}}^p=0,\hfill \end{array}\right.$$

(20)

under the boundary conditions

$$w^f=w^p=\theta =0\mathrm{at}z=0,1.$$

(21)

Let us denote by $\parallel ·\parallel$ and $\langle ·,·,\rangle$ the $L^2\left(V\right)$-norm and scalar product, respectively. Multiplying the first equation of the sytem (20) by $w^f$, the second equation by $w^p$, and the third equatuon by $\theta$; adding the resulting equations and integrating over V, results in

$${\displaystyle \frac{1}{2}}{\displaystyle \frac{d}{dt}}{\parallel \theta \parallel }^2=\mathrm{Ra}I-D,$$

(22)

with

$$\begin{array}{c}I=-\langle \tilde{f}\left(z\right)(w^f+w^p),\theta \rangle +\lambda \langle {\nabla }_1\theta ,{\nabla }_1(w^f+w^p)\rangle ,\hfill \\ D={\parallel \nabla \theta \parallel }^2+\lambda \left\{{\gamma }_1\parallel \nabla w^f{\parallel }^2+{\gamma }_2\parallel \nabla w^p{\parallel }^2+{\parallel \nabla w^f-\nabla w^p\parallel }^2\right\}\hfill \end{array}$$

(23)

with ${\nabla }_1$ as the horizontal gradient and $\lambda \in \mathbb{R}$ a coupling parameter to be suitably chosen later. Defining

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

(24)

with $\mathcal{H}$ being the class of the kinematically admissible perturbations, i.e.,

$$\begin{array}{cc}\hfill \mathcal{H}=\{& (w^f,w^p,\theta )\in {\left(H^1\right)}^3|w^f=w^p=\theta =0\mathrm{on}z=0,1;\mathrm{periodic}\mathrm{in}x,y\hfill \\ & \mathrm{with}\mathrm{periods}2\pi /a_x,2\pi /a_y;D<\infty ;\mathrm{verifying}\nabla ·{\mathbf{u}}^s=0,s=\{f,p\}\}.\hfill \end{array}$$

$\mathrm{Ra}<R_E$ implies the nonlinear, global, asymptotic, exponential stability in the energy norm. The Euler–Lagrange equations are:

$$\left\{\begin{array}{c}{R}_E[\tilde{f}\left(z\right)\theta +\lambda {\mathsf{\Delta }}_1\theta ]+\lambda [-2(1+{\gamma }_1)\mathsf{\Delta }{w}^f+2\mathsf{\Delta }{w}^p]=0,\hfill \\ R_E[\tilde{f}\left(z\right)\theta +\lambda {\mathsf{\Delta }}_1\theta ]+\lambda [2\mathsf{\Delta }{w}^f-2(1+{\gamma }_2)\mathsf{\Delta }{w}^p]=0,\hfill \\ R_E[\tilde{f}\left(z\right)(w^f+w^p)+\lambda {\mathsf{\Delta }}_1(w^f+w^p)]-2\mathsf{\Delta }\theta =0.\hfill \end{array}\right.$$

(25)

By using the order-1 Galerkin residual method, one obtains an approximation of the solution of the system (25) given by

$$R_E^2={\displaystyle \frac{4\lambda ({\gamma }_1+{\gamma }_2+{\gamma }_1{\gamma }_2){(a^2+{\pi }^2)}^2{(4{\pi }^2+{\mathrm{Pe}}^2)}^2}{(4+{\gamma }_1+{\gamma }_2){\left[4{\pi }^2(1+\lambda a^2)+a^2\lambda {\mathrm{Pe}}^2\right]}^2}}.$$

(26)

Choosing

$$\lambda ={\displaystyle \frac{4{\pi }^2}{a^2(4{\pi }^2+{\mathrm{Pe}}^2)}},$$

(27)

to maximize $R_E^2$, one obtains:

$$R_E^2={\displaystyle \frac{({\gamma }_1+{\gamma }_2+{\gamma }_1{\gamma }_2){(a^2+{\pi }^2)}^2(4{\pi }^2+{\mathrm{Pe}}^2)}{4{\pi }^2{a}^2(4+{\gamma }_1+{\gamma }_2)}}.$$

(28)

The minimum—with respect to $a^2\in {\mathbb{R}}^{+}$—is given by ${\mathrm{Ra}}^{\left(s\right)}$, i.e., there is coincidence between the linear instability threshold and the global nonlinear stability threshold in the $L^2\left(V\right)-$norm.

5. Results and Conclusions

In this paper, the instability of vertical constant throughflows saturating a horizontal BDPM that is uniformly heated from below is analyzed. An approximation of the critical Darcy–Rayleigh thermal number for the onset of steady instability was determined in algebraic closed form by using the order-1 Galerkin approximation method. The coincidence between the linear instability threshold and the global nonlinear stability threshold in the energy norm was proven. It is found that instability set in when the Darcy–Rayleigh number reached the threshold in Equation (18). Then, from Equation (18), it follows that:

(i)

In the absence of throughflow (i.e., when the horizontal bidisperse porous medium is filled by a fluid at the rest state), ${\mathrm{Ra}}^{\left(s\right)}$ reverts to the critical Darcy–Rayleigh thermal number $R_S$ at which steady convection sets in.

(ii)

Since ${\mathrm{Ra}}^{\left(s\right)}>R_S,\forall Pe$, the throughflow has a stabilizing effect in the sense that it loses its stability for a higher Darcy–Rayleigh number, compared to that related to thermal conduction solution.

(iii)

In order to compare the result obtained here with the case of a monodispersive porous layer, let us define the classical Darcy–Rayleigh number ${\mathrm{Ra}}_{\mathrm{cl}}$ by

$${\mathrm{Ra}}_{\mathrm{cl}}^2={\displaystyle \frac{{\varrho }_{F}g\alpha \beta d^2{k}_f{\left(\varrho c\right)}_f}{k_m\mu }}.$$

(29)

Hence,

$${\mathrm{Ra}}^2={\gamma }_1{\mathrm{Ra}}_{\mathrm{cl}}^2.$$

(30)

Then, the critical Darcy–Rayleigh number for the onset of steady instability for a vertical throughflow saturating a monodispersive layer is obtained by substituting Equation (30) into Equation (18) and letting $\zeta \to 0$, i.e., it is given by

$${\mathrm{Ra}}_{\mathrm{cl},s}^2=4{\pi }^2+{\mathrm{Pe}}^2,$$

(31)

with $\mathrm{Pe}=Qd/k$. The threshold (31) coincides with the one found in [26] in the absence of chemicals dissolved in the fluid.

(iv)

Comparing Equation (18) with Equation (31), it turns out that (see Figure 1)

−

if ${\gamma }_1{\gamma }_2<4$, then ${\mathrm{Ra}}^{\left(s\right)}<{\mathrm{Ra}}_{\mathrm{cl},s}^2$, i.e., the double porosity has a destabilizing effect; and

−

if ${\gamma }_1{\gamma }_2>4$, then ${\mathrm{Ra}}^{\left(s\right)}>{\mathrm{Ra}}_{\mathrm{cl},s}^2$, i.e., the double porosity has a stabilizing effect.