2 Governing equations in the isotropic case

A general theory for bidisperse porous media is presented in Franchi, Nibbi & Straughan (Reference Franchi, Nibbi and Straughan2017). These writers combine the model of Nield & Kuznetsov (Reference Nield and Kuznetsov2006a ) for a bidisperse porous medium together with the local thermal non-equilibrium equations for a monodisperse porous material of Nield (Reference Nield1998) and Banu & Rees (Reference Banu and Rees2002), to derive a theory which involves macro and micropore averaged velocities, an individual solid skeleton, macro and micro temperature fields. This may be seen as follows.

For an isotropic bidisperse porous medium the momentum and continuity equations of Nield & Kuznetsov (Reference Nield and Kuznetsov2006a ) may be written

(2.1) $$\begin{eqnarray}\left.\begin{array}{@{}ll@{}} & \displaystyle -{\displaystyle \frac{\unicode[STIX]{x1D707}}{K_{f}}}\,U_{i}^{f}+\tilde{\unicode[STIX]{x1D707}}\unicode[STIX]{x0394}U_{i}^{f}-\unicode[STIX]{x1D701}(U_{i}^{f}-U_{i}^{p})-p_{,i}^{f}\\[5.0pt] & \displaystyle \quad +\unicode[STIX]{x1D70C}_{F}g\unicode[STIX]{x1D6FC}k_{i}\left[{\displaystyle \frac{\unicode[STIX]{x1D719}}{\unicode[STIX]{x1D719}+\unicode[STIX]{x1D700}(1-\unicode[STIX]{x1D719})}}\,T^{f}+{\displaystyle \frac{\unicode[STIX]{x1D700}(1-\unicode[STIX]{x1D719})}{\unicode[STIX]{x1D719}+\unicode[STIX]{x1D700}(1-\unicode[STIX]{x1D719})}}\,T^{p}\right]=0,\\[5.0pt] & \displaystyle U_{i,i}^{f}=0,\\[5.0pt] & \displaystyle -{\displaystyle \frac{\unicode[STIX]{x1D707}}{K_{p}}}\,U_{i}^{p}-\unicode[STIX]{x1D701}(U_{i}^{p}-U_{i}^{f})-p_{,i}^{p}\\[5.0pt] & \displaystyle \quad +\unicode[STIX]{x1D70C}_{F}g\unicode[STIX]{x1D6FC}k_{i}\left[{\displaystyle \frac{\unicode[STIX]{x1D719}}{\unicode[STIX]{x1D719}+\unicode[STIX]{x1D700}(1-\unicode[STIX]{x1D719})}}\,T^{f}+{\displaystyle \frac{\unicode[STIX]{x1D700}(1-\unicode[STIX]{x1D719})}{\unicode[STIX]{x1D719}+\unicode[STIX]{x1D700}(1-\unicode[STIX]{x1D719})}}\,T^{p}\right]=0,\\[5.0pt] & U_{i,i}^{p}=0,\end{array}\right\}\end{eqnarray}$$

where $k_{i}=\unicode[STIX]{x1D6FF}_{i3}$ and where the divergence free conditions (continuity equations) express the fact that the fluid is incompressible. Nield & Kuznetsov (Reference Nield and Kuznetsov2006a ) additionally include a Laplacian (Brinkman) term in equation (2.1) $_{3}$ but we only use Darcy theory in the microphase. In these equations $U_{i}^{f},U_{i}^{p}$ , $p^{f},p^{p}$ , $T^{f}$ and $T^{p}$ are the pore averaged velocities, pressure and temperature fields in the macro (f) phase and micro (p) phase, respectively. The quantities $\unicode[STIX]{x1D707},\tilde{\unicode[STIX]{x1D707}},\unicode[STIX]{x1D70C},K_{f}$ and $K_{p}$ denote the dynamic viscosity of the fluid, the Brinkman viscosity, the fluid density, the macropermeability and the micropermeability. In addition, $g,\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D70C}_{F}$ denote gravity, thermal expansion coefficient of the fluid and a constant reference density, since in the buoyancy terms a Boussinesq approximation is employed to write the fluid density as a linear function of $T^{f}$ and $T^{p}$ . Finally, $\unicode[STIX]{x1D701}$ is a coefficient which represents the momentum transfer between the macro and microphases. All the variables and parameters defined above are listed in table 1.

In the solid skeleton, the macrophase and in the microphase the energy balances may be separately written to account for interactions between the temperatures of each phase. Taking into account the volume fraction occupied by each phase the energy balances are written as

(2.2) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D716}_{1}(\unicode[STIX]{x1D70C}c)_{s}T_{,t}^{s}=\unicode[STIX]{x1D716}_{1}\unicode[STIX]{x1D705}_{s}\unicode[STIX]{x0394}T^{s}+s_{1}(T^{f}-T^{s})+s_{2}(T^{p}-T^{s}),\\[5.0pt] \displaystyle \unicode[STIX]{x1D719}(\unicode[STIX]{x1D70C}c)_{f}T_{,t}^{f}+\unicode[STIX]{x1D719}(\unicode[STIX]{x1D70C}c)_{f}V_{i}^{f}T_{,i}^{f}=\unicode[STIX]{x1D719}\unicode[STIX]{x1D705}_{f}\unicode[STIX]{x0394}T^{f}+h(T^{p}-T^{f})+s_{1}(T^{s}-T^{f}),\\[5.0pt] \displaystyle \unicode[STIX]{x1D716}_{2}(\unicode[STIX]{x1D70C}c)_{p}T_{,t}^{p}+\unicode[STIX]{x1D716}_{2}(\unicode[STIX]{x1D70C}c)_{p}V_{i}^{p}T_{,i}^{p}=\unicode[STIX]{x1D716}_{2}\unicode[STIX]{x1D705}_{p}\unicode[STIX]{x0394}T^{p}+h(T^{f}-T^{p})+s_{2}(T^{s}-T^{p}),\end{array}\right\}\end{eqnarray}$$

where

$$\begin{eqnarray}\unicode[STIX]{x1D716}_{1}=(1-\unicode[STIX]{x1D700})(1-\unicode[STIX]{x1D719}),\quad \unicode[STIX]{x1D716}_{2}=(1-\unicode[STIX]{x1D719})\unicode[STIX]{x1D700},\end{eqnarray}$$

and where $s,f$ and $p$ denote solid skeleton, macrophase, microphase, $\unicode[STIX]{x1D70C}c$ denotes the product of the respective density and specific heat at constant pressure and $\unicode[STIX]{x1D705}$ denotes the thermal diffusivity. In (2.2) the velocities $V_{i}^{f}$ and $V_{i}^{p}$ are the actual velocities of the fluid which would be witnessed in an individual element in each phase. The coefficient $h$ is analogous to the $h$ of local non-thermal equilibrium theory as in Banu & Rees (Reference Banu and Rees2002), or in Rees (Reference Rees2009, Reference Rees2010). The coefficients $s_{1}$ and $s_{2}$ represent thermal interactions between the solid skeleton and the macro and microphases.

In the present work we deal with a local thermal equilibrium theory and take the temperatures $T^{s},T^{f}$ and $T^{p}$ to be the same, namely $T=T^{s}=T^{f}=T^{p}$ . In this way we may follow the procedure advocated by Joseph (Reference Joseph1976, p. 55) for a monodisperse porous material and we derive the governing system of equations from (2.1) and (2.2) as

(2.3) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \tilde{\unicode[STIX]{x1D707}}\unicode[STIX]{x0394}U_{i}^{f}-{\displaystyle \frac{\unicode[STIX]{x1D707}}{K_{f}}}U_{i}^{f}-\unicode[STIX]{x1D701}(U_{i}^{f}-U_{i}^{p})-p_{,i}^{f}+\unicode[STIX]{x1D70C}_{F}\unicode[STIX]{x1D6FC}gTk_{i}=0,\\[5.0pt] \displaystyle U_{i,i}^{f}=0,\\[5.0pt] \displaystyle -{\displaystyle \frac{\unicode[STIX]{x1D707}}{K_{p}}}U_{i}^{p}-\unicode[STIX]{x1D701}(U_{i}^{p}-U_{i}^{f})-p_{,i}^{p}+\unicode[STIX]{x1D70C}_{F}\unicode[STIX]{x1D6FC}gTk_{i}=0,\\[5.0pt] \displaystyle U_{i,i}^{p}=0,\\[5.0pt] \displaystyle (\unicode[STIX]{x1D70C}c)_{m}T_{,t}+(\unicode[STIX]{x1D70C}c)_{f}(U_{i}^{f}+U_{i}^{p})T_{,i}=\unicode[STIX]{x1D705}_{m}\unicode[STIX]{x0394}T.\end{array}\right\}\end{eqnarray}$$

Here, $U_{i}^{f}$ and $U_{i}^{p}$ are the pore averaged velocities, namely $U_{i}^{f}=\unicode[STIX]{x1D719}V_{i}^{f},$ $U_{i}^{p}=\unicode[STIX]{x1D700}(1-\unicode[STIX]{x1D719})V_{i}^{p}$ and $(\unicode[STIX]{x1D70C}c)_{m}$ and $\unicode[STIX]{x1D705}_{m}$ are given by

$$\begin{eqnarray}\displaystyle & \displaystyle (\unicode[STIX]{x1D70C}c)_{m}=(1-\unicode[STIX]{x1D719})(1-\unicode[STIX]{x1D716})(\unicode[STIX]{x1D70C}c)_{s}+\unicode[STIX]{x1D719}(\unicode[STIX]{x1D70C}c)_{f}+\unicode[STIX]{x1D716}(1-\unicode[STIX]{x1D719})(\unicode[STIX]{x1D70C}c)_{p}, & \displaystyle \nonumber\\ \displaystyle & \displaystyle \unicode[STIX]{x1D705}_{m}=(1-\unicode[STIX]{x1D719})(1-\unicode[STIX]{x1D716})\unicode[STIX]{x1D705}_{s}+\unicode[STIX]{x1D719}\unicode[STIX]{x1D705}_{f}+\unicode[STIX]{x1D716}(1-\unicode[STIX]{x1D719})\unicode[STIX]{x1D705}_{p}. & \displaystyle \nonumber\end{eqnarray}$$

As we remarked in the introduction, Imani & Hooman (Reference Imani and Hooman2017) have provided a justification for considering the temperatures to have the same value $T$ . The theory of Nield & Kuznetsov (Reference Nield and Kuznetsov2006a ) allows for different temperatures in the macro and micropores, but for problems involving the onset of thermal convection, since it is the same fluid in each type of pore, the thermal properties are the same and we believe this justifies the use of a single temperature. Note that we are assuming a relatively large porosity in the macrostructure and hence equation (2.3) $_{1}$ includes a Brinkman term, but we only require Darcy theory at the microlevel since the porosity there may be smaller. We believe that this is the first study of flow in bidisperse porous media with (2.3).

3 General theory for thermal convection

We suppose now the bidisperse porous medium is contained in the horizontal layer $0<z<d$ . In order to describe thermal convection we allow a generalization of (2.3) to the situation where the bidisperse porous medium may be anisotropic, which results in the permeabilities and mass transfer coefficients being second-order tensors. The problem of convection in an anisotropic bidisperse porous medium which employs Darcy theory in both the macro and microphases has led to some novel convection cell behaviour, see Straughan (Reference Straughan2018, Reference Straughan2019a ,Reference Straughan b ). Our first goal is to prove a general result showing coincidence of the linear instability and nonlinear stability boundaries and this we can do in a general anisotropic setting.

The relevant equations from (2.3) allowing for the solid skeleton to have an anisotropic structure have the form

(3.1) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \tilde{\unicode[STIX]{x1D707}}\unicode[STIX]{x0394}U_{i}^{f}-\unicode[STIX]{x1D614}_{ij}^{f}U_{j}^{f}-\unicode[STIX]{x1D701}_{ij}(U_{j}^{f}-U_{j}^{p})-p_{,i}^{f}+\unicode[STIX]{x1D70C}_{F}\unicode[STIX]{x1D6FC}gTk_{i}=0,\\[5.0pt] \displaystyle U_{i,i}^{f}=0,\\[5.0pt] \displaystyle -\unicode[STIX]{x1D614}_{ij}^{p}U_{j}^{p}-\unicode[STIX]{x1D701}_{ij}(U_{j}^{p}-U_{j}^{f})-p_{,i}^{p}+\unicode[STIX]{x1D70C}_{F}\unicode[STIX]{x1D6FC}gTk_{i}=0,\\[5.0pt] \displaystyle U_{i,i}^{p}=0,\\[5.0pt] \displaystyle (\unicode[STIX]{x1D70C}c)_{m}T_{,t}+(\unicode[STIX]{x1D70C}c)_{f}(U_{i}^{f}+U_{i}^{p})T_{,i}=\unicode[STIX]{x1D705}_{m}\unicode[STIX]{x0394}T.\end{array}\right\}\end{eqnarray}$$

The term $\tilde{\unicode[STIX]{x1D707}}$ is the Brinkman viscosity, $\unicode[STIX]{x1D6E5}$ is the Laplace operator and $\unicode[STIX]{x1D614}_{ij}^{f},\unicode[STIX]{x1D614}_{ij}^{p}$ are symmetric tensors given by

$$\begin{eqnarray}\unicode[STIX]{x1D614}_{ij}^{f}=\unicode[STIX]{x1D707}(\unicode[STIX]{x1D612}_{ij}^{f})^{-1},\quad \unicode[STIX]{x1D614}_{ij}^{p}=\unicode[STIX]{x1D707}(\unicode[STIX]{x1D612}_{ij}^{p})^{-1},\end{eqnarray}$$

$\unicode[STIX]{x1D707}$ being the dynamic viscosity of the saturating fluid, while $\unicode[STIX]{x1D612}_{ij}^{f}$ and $\unicode[STIX]{x1D612}_{ij}^{p}$ are the permeability tensors in the macro and microphases. The term $\unicode[STIX]{x1D701}_{ij}$ is a symmetric tensor representing the momentum transfer term but allowing for an anisotropic nature.

The boundary conditions imposed are that

(3.2a,b ) $$\begin{eqnarray}U_{i}^{f}=0,\quad U_{3}^{p}=0,\quad \text{on }z=0,d,\end{eqnarray}$$

and

(3.3) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D6FC}_{L}\left({\displaystyle \frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}z}}+\unicode[STIX]{x1D6FD}\right)d-(1-\unicode[STIX]{x1D6FC}_{L})T=-(1-\unicode[STIX]{x1D6FC}_{L})T_{L},\quad z=0,\\[10.0pt] \displaystyle \unicode[STIX]{x1D6FC}_{U}\left({\displaystyle \frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}z}}+\unicode[STIX]{x1D6FD}\right)d+(1-\unicode[STIX]{x1D6FC}_{U})T=(1-\unicode[STIX]{x1D6FC}_{U})T_{U},\quad z=d,\end{array}\right\}\end{eqnarray}$$

where $T_{L}>T_{U}>0$ are constants and $0\leqslant \unicode[STIX]{x1D6FC}_{L}<1$ , $0\leqslant \unicode[STIX]{x1D6FC}_{U}<1$ , are likewise constants. The parameter $\unicode[STIX]{x1D6FD}=(T_{L}-T_{U})/d$ . In § 6, where we present numerical results, we restrict our attention to $\unicode[STIX]{x1D6FC}_{L}=0=\unicode[STIX]{x1D6FC}_{U}$ , i.e. prescribed temperatures on $z=0$ and $d$ . However, for our general result on nonlinear stability it is expedient to allow the more general boundary conditions. The $\unicode[STIX]{x2202}T/\unicode[STIX]{x2202}z$ term is important when radiation boundary effects are acting, such as solar heating. Prescribed temperatures on the boundaries are likely under laboratory conditions or under cloudy skies, and (3.3) encompass the general situation since $\unicode[STIX]{x1D6FC}_{L}$ and $\unicode[STIX]{x1D6FC}_{U}$ may be taken in the range $[0,1)$ .

In order to investigate thermal convection we observe that (3.1), (3.2) and (3.3) admit the steady (conduction) solution

(3.4a-c ) $$\begin{eqnarray}\bar{U}_{i}^{f}\equiv 0,\quad \bar{U}_{i}^{p}\equiv 0,\quad \bar{T}=T_{L}-\unicode[STIX]{x1D6FD}z.\end{eqnarray}$$

To analyse the stability of solution (3.4) we let $u_{i}^{f},\,u_{i}^{p},\,\unicode[STIX]{x1D703},\,\unicode[STIX]{x03C0}^{f},\,\unicode[STIX]{x03C0}^{p}$ be perturbations to the steady solution. We then non-dimensionalize the resulting perturbation equations with the scalings

$$\begin{eqnarray}x_{i}=x_{i}^{\ast }d,\quad t=t^{\ast }{\mathcal{T}},\quad {\mathcal{T}}={\displaystyle \frac{(\unicode[STIX]{x1D70C}c)_{m}d^{2}}{\unicode[STIX]{x1D705}_{m}}},\quad U={\displaystyle \frac{\unicode[STIX]{x1D705}_{m}}{(\unicode[STIX]{x1D70C}c)_{f}d}},\quad P=dm_{11}U,\end{eqnarray}$$

where ${\mathcal{T}}$ , $U$ and $P$ are time, velocity and pressure scales. The number $m_{11}$ is given by $m_{11}=\unicode[STIX]{x1D614}_{11}^{f}$ and we take this out as a factor in $\unicode[STIX]{x1D614}_{ij}^{f}$ . We also write $\unicode[STIX]{x1D614}_{11}^{p}=\unicode[STIX]{x1D714}m_{11}$ and rescale $\unicode[STIX]{x1D614}_{ij}^{p}$ . Then put $\unicode[STIX]{x1D706}_{ij}=\unicode[STIX]{x1D701}_{ij}/m_{11}$ , $\unicode[STIX]{x1D706}=\tilde{\unicode[STIX]{x1D707}}/d^{2}m_{11}$ , define the temperature scale as

$$\begin{eqnarray}T^{\sharp }=\text{d}U\sqrt{{\displaystyle \frac{(\unicode[STIX]{x1D70C}c)_{f}\unicode[STIX]{x1D6FD}m_{11}}{\unicode[STIX]{x1D70C}_{F}\unicode[STIX]{x1D6FC}g\unicode[STIX]{x1D705}_{m}}}}\end{eqnarray}$$

and define the Rayleigh number $Ra=R^{2}$ by

(3.5) $$\begin{eqnarray}Ra={\displaystyle \frac{\unicode[STIX]{x1D70C}_{F}\unicode[STIX]{x1D6FC}g\unicode[STIX]{x1D6FD}d^{2}}{m_{11}[\unicode[STIX]{x1D705}_{m}/(\unicode[STIX]{x1D70C}c)_{f}]}}.\end{eqnarray}$$

The non-dimensionalization allows us to reduce the number of parameters in (3.1). For example, in the isotropic case there are in (2.3) the parameters $\unicode[STIX]{x1D707},\tilde{\unicode[STIX]{x1D707}},\mathit{K}_{f},\mathit{K}_{p},\unicode[STIX]{x1D701},\unicode[STIX]{x1D70C}_{F},g$ , $\unicode[STIX]{x1D6FC},(\unicode[STIX]{x1D70C}c)_{m},(\unicode[STIX]{x1D70C}c)_{f}$ and $\unicode[STIX]{x1D705}_{m}$ and the non-dimensionalization reduces this to four parameters, a Rayleigh number, and non-dimensional versions of $\unicode[STIX]{x1D706},\unicode[STIX]{x1D701}$ and $K_{r}=K_{f}/K_{p}$ the relative permeability, see § 6. We are able to reduce the behaviour of thermal convection in that case to a study of the dependence of the Rayleigh number, $Ra$ , and the wavenumber upon the three parameters $\unicode[STIX]{x1D706},\unicode[STIX]{x1D701}$ and $K_{r}$ .

The resulting perturbation equations arising from (3.1) may be written as

(3.6) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D706}\unicode[STIX]{x0394}u_{i}^{\,f}-\unicode[STIX]{x1D614}_{ij}^{f}u_{j}^{\,f}-\unicode[STIX]{x1D706}_{ij}(u_{i}^{\,f}-u_{i}^{p})-\unicode[STIX]{x03C0}_{,i}^{f}+Rk_{i}\unicode[STIX]{x1D703}=0,\\[5.0pt] u_{i,i}^{\,f}=0,\\[5.0pt] -\unicode[STIX]{x1D714}\unicode[STIX]{x1D614}_{ij}^{p}u_{j}^{p}+\unicode[STIX]{x1D706}_{ij}(u_{i}^{\,f}-u_{i}^{p})-\unicode[STIX]{x03C0}_{,i}^{p}+Rk_{i}\unicode[STIX]{x1D703}=0,\\[5.0pt] u_{i,i}^{p}=0,\\[5.0pt] \unicode[STIX]{x1D703}_{,t}+(u_{i}^{\,f}+u_{i}^{p})\unicode[STIX]{x1D703}_{,i}=R(w^{f}+w^{p})+\unicode[STIX]{x0394}\unicode[STIX]{x1D703},\end{array}\right\}\end{eqnarray}$$

where $\boldsymbol{u}^{f}=(u^{\,f},v^{f},w^{f})$ and $\boldsymbol{u}^{p}=(u^{p},v^{p},w^{p})$ . These equations hold in the domain $(x,y)\in \mathbb{R}^{2},$ $\{z\in (0,1)\},$ $t>0.$ The boundary conditions are

(3.7) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle u_{i}^{f}=0,\quad w^{p}=0,\quad z=0,1,\\[5.0pt] \displaystyle \unicode[STIX]{x1D703}_{z}-L_{L}\unicode[STIX]{x1D703}=0,\quad z=0,\\[5.0pt] \displaystyle \unicode[STIX]{x1D703}_{z}+L_{U}\unicode[STIX]{x1D703}=0,\quad z=1,\end{array}\right\}\end{eqnarray}$$

where $L_{L}=(1-\unicode[STIX]{x1D6FC}_{L})/\unicode[STIX]{x1D6FC}_{L}$ and $L_{U}=(1-\unicode[STIX]{x1D6FC}_{U})/\unicode[STIX]{x1D6FC}_{U}$ . In addition, we suppose $u_{i}^{f},\,u_{i}^{p},\,\unicode[STIX]{x1D703},\,\unicode[STIX]{x03C0}^{f}$ , $\unicode[STIX]{x03C0}^{p}$ satisfy a plane tiling shape in the $(x,y)$ plane. The periodic convection cell which arises will be denoted by $V.$

4 Instability and exchange of stabilities

To determine the linear instability boundary we discard the nonlinear terms in (3.6) and seek a solution in which $u_{i}^{f},$ $u_{i}^{p},$ $\unicode[STIX]{x1D703},\,\unicode[STIX]{x03C0}^{f},\,\unicode[STIX]{x03C0}^{f}$ have time dependence like $e^{\unicode[STIX]{x1D70E}t}.$ The system which remains is written as

(4.1) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D706}\unicode[STIX]{x0394}u_{i}^{\,f}-\unicode[STIX]{x1D614}_{ij}^{f}u_{j}^{\,f}-\unicode[STIX]{x1D706}_{ij}(u_{j}^{\,f}-u_{j}^{p})-\unicode[STIX]{x03C0}_{,i}^{f}+Rk_{i}\unicode[STIX]{x1D703}=0,\\[5.0pt] u_{i,i}^{\,f}=0,\\[5.0pt] -\unicode[STIX]{x1D714}\unicode[STIX]{x1D614}_{ij}^{p}u_{j}^{p}+\unicode[STIX]{x1D706}_{ij}(u_{j}^{\,f}-u_{j}^{p})-\unicode[STIX]{x03C0}_{,i}^{p}+Rk_{i}\unicode[STIX]{x1D703}=0,\\[5.0pt] u_{i,i}^{p}=0,\\[5.0pt] \unicode[STIX]{x1D70E}\unicode[STIX]{x1D703}=R(w^{f}+w^{p})+\unicode[STIX]{x0394}\unicode[STIX]{x1D703}.\end{array}\right\}\end{eqnarray}$$

The boundary conditions are (3.7), together with periodicity in $(x,y)$ . We now establish the strong form of the principle of exchange of stabilities. To achieve this, multiply equation (4.1) $_{1}$ by $u_{i}^{f\ast },$ the complex conjugate of $u_{i}^{f},$ and integrate over the period cell $V.$ Next, multiply (4.1) $_{2}$ by $u_{i}^{p\ast }$ and multiply (4.1) $_{3}$ by $\unicode[STIX]{x1D703}^{\ast }$ and likewise integrate each resulting equation over $V.$ Denote by $\langle \cdot \rangle$ integration over $V$ , and let $\unicode[STIX]{x1D6E4}_{1}$ be the boundary of $V$ which intersects the plane $z=0$ while $\unicode[STIX]{x1D6E4}_{2}$ is the boundary of $V$ which intersects the plane $z=1$ . Upon performing some integration by parts and using the boundary conditions one may derive the equations

(4.2) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}-\unicode[STIX]{x1D706}\langle u_{i,j}^{\,f}u_{i,j}^{\,f\ast }\rangle -\langle \unicode[STIX]{x1D614}_{ij}^{f}u_{j}^{\,f}u_{i}^{\,f\ast }\rangle \\[5.0pt] \quad \quad -\langle \unicode[STIX]{x1D706}_{ij}(u_{j}^{\,f}-u_{j}^{p})u_{i}^{\,f\ast }\rangle +R\langle \unicode[STIX]{x1D703}w^{\,f\ast }\rangle =0,\\[5.0pt] -\unicode[STIX]{x1D714}\langle \unicode[STIX]{x1D614}_{ij}^{p}u_{j}^{p}u_{i}^{p\ast }\rangle -\langle \unicode[STIX]{x1D706}_{ij}(u_{j}^{p}-u_{j}^{\,f})u_{i}^{p\ast }\rangle +R\langle \unicode[STIX]{x1D703}w^{p\ast }\rangle =0,\\[5.0pt] \unicode[STIX]{x1D70E}\langle \unicode[STIX]{x1D703}\unicode[STIX]{x1D703}^{\ast }\rangle =R[\langle w^{f}\unicode[STIX]{x1D703}^{\ast }\rangle +\langle w^{p}\unicode[STIX]{x1D703}^{\ast }\rangle ]-\langle \unicode[STIX]{x1D703}_{,i}\unicode[STIX]{x1D703}_{,i}^{\ast }\rangle \\[5.0pt] \phantom{\unicode[STIX]{x1D70E}\langle \unicode[STIX]{x1D703}\unicode[STIX]{x1D703}^{\ast }\rangle =}+\displaystyle \int _{\unicode[STIX]{x1D6E4}_{2}}\unicode[STIX]{x1D703}^{\ast }\,{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}z}}\,\,\text{d}A-\displaystyle \int _{\unicode[STIX]{x1D6E4}_{1}}\unicode[STIX]{x1D703}^{\ast }\,{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}z}}\,\text{d}A.\end{array}\right\}\end{eqnarray}$$

Next, add the above three equations and use the boundary conditions to obtain

(4.3) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D70E}\langle \unicode[STIX]{x1D703}\unicode[STIX]{x1D703}^{\ast }\rangle =-\unicode[STIX]{x1D706}\langle u_{i,j}^{\,f}u_{i,j}^{\,f\ast }\rangle -\langle \unicode[STIX]{x1D614}_{ij}^{f}u_{j}^{\,f}u_{i}^{\,f\ast }\rangle \\[5.0pt] \displaystyle -\unicode[STIX]{x1D714}\langle \unicode[STIX]{x1D614}_{ij}^{p}u_{j}^{p}u_{i}^{p\ast }\rangle -\langle \unicode[STIX]{x1D706}_{ij}(u_{j}^{\,f}-u_{j}^{p})(u_{i}^{\,f\ast }-u_{i}^{p\ast })\rangle \\[5.0pt] +R[\langle \unicode[STIX]{x1D703}w^{\,f\ast }\rangle +\langle w^{f}\unicode[STIX]{x1D703}^{\ast }\rangle +\langle \unicode[STIX]{x1D703}w^{p\ast }\rangle +\langle w^{p}\unicode[STIX]{x1D703}^{\ast }\rangle ]\\[5.0pt] \displaystyle -\langle \unicode[STIX]{x1D703}_{,i}\unicode[STIX]{x1D703}_{,i}^{\ast }\rangle -L_{L}\int _{\unicode[STIX]{x1D6E4}_{1}}\unicode[STIX]{x1D703}^{\ast }\unicode[STIX]{x1D703}\,\text{d}A-L_{U}\int _{\unicode[STIX]{x1D6E4}_{2}}\unicode[STIX]{x1D703}^{\ast }\unicode[STIX]{x1D703}\,\text{d}A.\end{array}\right\}\end{eqnarray}$$

One now writes each of $u_{i}^{\,f},u_{i}^{p},\unicode[STIX]{x1D703}$ as a sum of their real and imaginary parts and we put $\unicode[STIX]{x1D70E}=\unicode[STIX]{x1D70E}_{r}+i\unicode[STIX]{x1D70E}_{1}.$ The imaginary part of equation (4.3) yields the result

$$\begin{eqnarray}\unicode[STIX]{x1D70E}_{1}\langle \unicode[STIX]{x1D703}\unicode[STIX]{x1D703}^{\ast }\rangle =0.\end{eqnarray}$$

One requires $\langle \unicode[STIX]{x1D703}\unicode[STIX]{x1D703}^{\ast }\rangle \neq 0$ and so $\unicode[STIX]{x1D70E}_{1}=0.$ Thus, $\unicode[STIX]{x1D70E}\in \mathbb{R}$ and hence oscillatory convection does not occur. In conclusion, the strong form of the principle of exchange of stabilities is demonstrated. To find the linear instability boundary one may now investigate equations (4.1) with $\unicode[STIX]{x1D70E}=0$ .

We stress that in the two temperature case of Nield & Kuznetsov (Reference Nield and Kuznetsov2006a ) the exchange of stabilities has not been proved.

We return to calculate the instability threshold explicitly in § 6. There, we restrict attention to the case where $\unicode[STIX]{x1D6FC}_{U}=0,\unicode[STIX]{x1D6FC}_{L}=0$ , i.e. prescribed lower temperature $T_{L}$ , upper temperature $T_{U}$ . Detailed numerical results are presented.

5 Global nonlinear stability

Up to this point we have discussed only linear instability. We now establish an important result for global nonlinear stability (i.e. for all initial data). To achieve this result we employ the full nonlinear equations (3.6). We begin by multiplying (3.6) $_{1}$ by $u_{i}^{\,f}$ and integrating over  $V$ . Then, we multiply (3.6) $_{3}$ by $u_{i}^{p}$ and integrate over  $V$ , and finally we multiply (3.6) $_{5}$ by $\unicode[STIX]{x1D703}$ and integrate over $V$ . We add the results for the first two resulting equations and then after some integration by parts and use of the boundary conditions we obtain the following two results

(5.1) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}-\unicode[STIX]{x1D706}\langle u_{i,j}^{\,f}u_{i,j}^{\,f}\rangle -\langle \unicode[STIX]{x1D614}_{ij}^{f}u_{j}^{\,f}u_{i}^{\,f}\rangle -\unicode[STIX]{x1D714}\langle \unicode[STIX]{x1D614}_{ij}^{p}u_{j}^{p}u_{i}^{p}\rangle \\[5.0pt] \quad \quad -\langle \unicode[STIX]{x1D706}_{ij}(u_{j}^{\,f}-u_{j}^{p})(u_{i}^{\,f}-u_{i}^{p})\rangle +R\langle \unicode[STIX]{x1D703}(w^{f}+w^{p})\rangle =0,\\[5.0pt] {\displaystyle \frac{\text{d}}{\text{d}t}}{\displaystyle \frac{1}{2}}\langle \unicode[STIX]{x1D703}^{2}\rangle =R\langle \unicode[STIX]{x1D703}(w^{f}+w^{p})\rangle -\langle \unicode[STIX]{x1D703}_{,i}\unicode[STIX]{x1D703}_{,i}\rangle \\[5.0pt] \phantom{{\displaystyle \frac{\text{d}}{\text{d}t}}{\displaystyle \frac{1}{2}}\langle \unicode[STIX]{x1D703}^{2}\rangle =}-L_{L}\displaystyle \int _{\unicode[STIX]{x1D6E4}_{1}}\unicode[STIX]{x1D703}^{2}\,\text{d}A-L_{U}\displaystyle \int _{\unicode[STIX]{x1D6E4}_{2}}\unicode[STIX]{x1D703}^{2}\,\text{d}A.\end{array}\right\}\end{eqnarray}$$

Add the equations in (5.1) to obtain

(5.2) $$\begin{eqnarray}{\displaystyle \frac{\text{d}}{\text{d}t}}\,{\displaystyle \frac{1}{2}}\!\langle \unicode[STIX]{x1D703}^{2}\rangle =RI-D,\end{eqnarray}$$

where now the production term $I$ is

(5.3) $$\begin{eqnarray}I=2\langle (w^{f}+w^{p})\unicode[STIX]{x1D703}\rangle ,\end{eqnarray}$$

and the dissipation term $D$ is given by

(5.4) $$\begin{eqnarray}\left.\begin{array}{@{}rcl@{}}D\, & =\, & \unicode[STIX]{x1D706}\langle u_{i,j}^{\,f}u_{i,j}^{\,f}\rangle +\langle \unicode[STIX]{x1D614}_{ij}^{f}u_{i}^{\,f}u_{j}^{\,f}\rangle +\unicode[STIX]{x1D714}\langle \unicode[STIX]{x1D614}_{ij}^{p}u_{i}^{p}u_{j}^{p}\rangle \\[5.0pt] \, & \, & +\,\langle \unicode[STIX]{x1D706}_{ij}(u_{j}^{\,f}-u_{j}^{p})(u_{i}^{\,f}-u_{i}^{p})\rangle +\langle \unicode[STIX]{x1D703}_{,i}\unicode[STIX]{x1D703}_{,i}\rangle \\[5.0pt] \, & \, & +\,L_{U}\displaystyle \int _{\unicode[STIX]{x1D6E4}_{2}}\unicode[STIX]{x1D703}^{2}\,\text{d}A+L_{L}\displaystyle \int _{\unicode[STIX]{x1D6E4}_{1}}\unicode[STIX]{x1D703}^{2}\,\text{d}A.\end{array}\right\}\end{eqnarray}$$

Define the threshold $R_{E}$ by

(5.5) $$\begin{eqnarray}{\displaystyle \frac{1}{R_{E}}}=\max _{H}{\displaystyle \frac{I}{D}},\end{eqnarray}$$

where $H$ is the space of admissible functions. From the energy identity (5.2) we may see that

(5.6) $$\begin{eqnarray}{\displaystyle \frac{\text{d}}{\text{d}t}}\,{\displaystyle \frac{1}{2}}\,\langle \unicode[STIX]{x1D703}^{2}\rangle \,\leqslant \,-D\left(1-{\displaystyle \frac{R}{R_{E}}}\right).\end{eqnarray}$$

If now $R<R_{E},$ say $1-R/R_{E}=b>0,$ then by using Poincaré’s inequality in (5.6) one may obtain

(5.7) $$\begin{eqnarray}{\displaystyle \frac{\text{d}}{\text{d}t}}\,{\displaystyle \frac{1}{2}}\,\langle \unicode[STIX]{x1D703}^{2}\rangle \leqslant -b\unicode[STIX]{x03C0}^{2}\langle \unicode[STIX]{x1D703}^{2}\rangle .\end{eqnarray}$$

This may be integrated to see that

(5.8) $$\begin{eqnarray}\langle \unicode[STIX]{x1D703}(t)^{2}\rangle \leqslant \langle \unicode[STIX]{x1D703}(0)^{2}\rangle \,\exp (-2b\unicode[STIX]{x03C0}^{2}t).\end{eqnarray}$$

From inequality (5.8) one deduces that $\langle \unicode[STIX]{x1D703}(t)^{2}\rangle$ decays exponentially provided $R<R_{E}.$

To also establish decay of $u_{i}^{\,f}$ and $u_{i}^{p}$ we require $\unicode[STIX]{x1D614}_{ij}^{f}$ and $\unicode[STIX]{x1D614}_{ij}^{p}$ to be positive–definite and $\unicode[STIX]{x1D706}_{ij}$ to be non-negative. These are realistic physical assumptions. If

$$\begin{eqnarray}\unicode[STIX]{x1D614}_{ij}^{f}\unicode[STIX]{x1D709}_{i}\unicode[STIX]{x1D709}_{j}\geqslant k^{f}\unicode[STIX]{x1D709}_{i}\unicode[STIX]{x1D709}_{i},\quad \text{and}\quad \unicode[STIX]{x1D614}_{ij}^{p}\unicode[STIX]{x1D709}_{i}\unicode[STIX]{x1D709}_{j}\geqslant k^{p}\unicode[STIX]{x1D709}_{i}\unicode[STIX]{x1D709}_{i},\quad \text{for all }\unicode[STIX]{x1D709}_{i},\end{eqnarray}$$

$k^{f}>0,k^{p}>0$ , then from equation (5.1) $_{1}$ we may derive using the arithmetic–geometric mean inequality

(5.9) $$\begin{eqnarray}\displaystyle & & \displaystyle k^{f}\langle u_{i}^{\,f}u_{i}^{\,f}\rangle +\,k^{p}\unicode[STIX]{x1D714}\langle u_{i}^{p}u_{i}^{p}\rangle +\unicode[STIX]{x1D706}\langle u_{i,j}^{\,f}u_{i,j}^{\,f}\rangle \leqslant {\displaystyle \frac{R\unicode[STIX]{x1D6FE}_{1}}{2}}\langle w_{f}^{2}\rangle \nonumber\\ \displaystyle & & \displaystyle \quad +\,{\displaystyle \frac{R\unicode[STIX]{x1D6FE}_{2}}{2}}\langle w_{p}^{2}\rangle +{\displaystyle \frac{R}{2}}\left({\displaystyle \frac{1}{\unicode[STIX]{x1D6FE}_{1}}}+{\displaystyle \frac{1}{\unicode[STIX]{x1D6FE}_{2}}}\right)\langle \unicode[STIX]{x1D703}^{2}\rangle ,\end{eqnarray}$$

for constants $\unicode[STIX]{x1D6FE}_{1}>0,\unicode[STIX]{x1D6FE}_{2}>0$ . Choose $\unicode[STIX]{x1D6FE}_{1}=k^{f}/R$ and $\unicode[STIX]{x1D6FE}_{2}=k^{p}\unicode[STIX]{x1D714}/R$ . Then (5.9) allows us to derive

(5.10) $$\begin{eqnarray}k^{f}\langle u_{i}^{\,f}u_{i}^{\,f}\rangle +k^{p}\unicode[STIX]{x1D714}\langle u_{i}^{p}u_{i}^{p}\rangle +2\unicode[STIX]{x1D706}\langle u_{i,j}^{\,f}u_{i,j}^{\,f}\rangle \leqslant K\langle \unicode[STIX]{x1D703}^{2}\rangle ,\end{eqnarray}$$

where $K=R^{2}(1/k^{f}+1/k^{p}\unicode[STIX]{x1D714})$ . Since (5.8) shows $\langle \unicode[STIX]{x1D703}^{2}\rangle$ decays exponentially in time we may deduce $R<R_{E}$ is sufficient also for decay of $u_{i}^{\,f}$ and $u_{i}^{p}$ as shown in (5.10).

Hence, the number $R_{E}$ represents a global nonlinear stability threshold. We wish to compare this number to the analogous one of linear instability. To do this we derive the Euler–Lagrange equations for the maximum in (5.5). This follows a standard procedure whereby one replaces $u_{i}^{\,f}$ by $u_{i}^{\,f}+\unicode[STIX]{x1D716}\unicode[STIX]{x1D702}_{i}^{f}$ , $u_{i}^{p}$ by $u_{i}^{p}+\unicode[STIX]{x1D716}\unicode[STIX]{x1D702}_{i}^{p}$ , and $\unicode[STIX]{x1D703}$ by $\unicode[STIX]{x1D703}+\unicode[STIX]{x1D716}\unicode[STIX]{x1D702}$ and then one differentiates $I/D$ with respect to $\unicode[STIX]{x1D716}$ and takes the limit $\unicode[STIX]{x1D716}\rightarrow 0$ . Upon completing this procedure we find for Lagrange multipliers $\unicode[STIX]{x1D714}^{f}$ and $\unicode[STIX]{x1D714}^{p}$ the Euler–Lagrange equations are

(5.11) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D706}\unicode[STIX]{x0394}u_{i}^{\,f}-\unicode[STIX]{x1D614}_{ij}^{f}u_{j}^{\,f}-\unicode[STIX]{x1D706}_{ij}(u_{j}^{\,f}-u_{j}^{p})+R_{E}k_{i}\unicode[STIX]{x1D703}=\unicode[STIX]{x1D714}_{,i}^{f},\\[5.0pt] \displaystyle u_{i,i}^{\,f}=0,\\[5.0pt] \displaystyle -\unicode[STIX]{x1D714}\unicode[STIX]{x1D614}_{ij}^{p}u_{j}^{p}+\unicode[STIX]{x1D706}_{ij}(u_{j}^{\,f}-u_{j}^{p})+R_{E}k_{i}\unicode[STIX]{x1D703}=\unicode[STIX]{x1D714}_{,i}^{p},\\[5.0pt] \displaystyle u_{i,i}^{p}=0,\\[5.0pt] \displaystyle R_{E}(w^{f}+w^{p})+\unicode[STIX]{x0394}\unicode[STIX]{x1D703}=0,\end{array}\right\}\end{eqnarray}$$

together with the boundary conditions

(5.12a-c ) $$\begin{eqnarray}u_{i}^{\,f}=0,\quad w^{p}=0,\quad z=0,1,\end{eqnarray}$$

and

(5.13) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D703}_{z}+L_{U}\unicode[STIX]{x1D703}=0,\quad z=1,\\[5.0pt] \displaystyle \unicode[STIX]{x1D703}_{z}-L_{L}\unicode[STIX]{x1D703}=0,\quad z=0.\end{array}\right\}\end{eqnarray}$$

Since exchange of stabilities holds we may take $\unicode[STIX]{x1D70E}=0$ in (4.1) and we see that (5.11)–(5.13) are exactly the same as (3.7)–(4.1) for the linear instability threshold. Thus, we have shown that the linear instability problem yields the same Rayleigh number threshold as the fully nonlinear stability one. This is an optimal result which shows that the linear theory is capturing correctly the physics for the onset of convective motion. In § 8 we find the linear instability thresholds and we know these also represent the nonlinear stability ones.

We again stress that in the two temperature case exchange of stabilities has not been proved. Straughan (Reference Straughan2009, Reference Straughan2015) analyses linear instability and nonlinear energy stability results in detail for the two temperature situation where Brinkman theory is employed in both macro and micropores, and in the case where Darcy theory is employed in both. In the case of the Brinkman theory the energy stability limit is lower than the linear instability one. For the Darcy case the energy limit is much lower. This does suggest that there may well be oscillatory instability in the Darcy case. We should point out that the equivalence of linear instability and nonlinear stability proved in the present situation is a strong result which is not to be expected in general, see Xiong & Chen (Reference Xiong and Chen2019).

6 Heated below isotropic theory, two free surfaces

In this section we solve the linear instability problem for (4.1). We restrict attention to the case of isotropic permeability tensors and so $\unicode[STIX]{x1D614}_{ij}^{f}\equiv (\unicode[STIX]{x1D707}/\mathit{K}_{f})\unicode[STIX]{x1D6FF}_{ij}$ and $\unicode[STIX]{x1D614}_{ij}^{p}\equiv (\unicode[STIX]{x1D707}/\mathit{K}_{p})\unicode[STIX]{x1D6FF}_{ij}$ where $\mathit{K}_{f}$ and $\mathit{K}_{p}$ are the values of the permeability in the macro and microphases. We also suppose the interaction term is isotropic and so $\unicode[STIX]{x1D706}_{ij}\equiv \unicode[STIX]{x1D709}\unicode[STIX]{x1D6FF}_{ij}$ , where $\unicode[STIX]{x1D709}\geqslant 0$ , is the non-dimensional interaction coefficient. We also restrict our attention to the case where $\unicode[STIX]{x1D6FC}_{L}$ and $\unicode[STIX]{x1D6FC}_{U}$ are zero so that the boundary conditions on the temperature are that

$$\begin{eqnarray}T=T_{L},\quad z=0;\quad T=T_{U},\quad z=d,\end{eqnarray}$$

in the original problem.

Thus, the relevant system of equations becomes

(6.1) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D706}\unicode[STIX]{x0394}u_{i}^{\,f}-u_{i}^{\,f}-\unicode[STIX]{x1D709}(u_{i}^{\,f}-u_{i}^{p})-\unicode[STIX]{x03C0}_{,i}^{f}+Rk_{i}\unicode[STIX]{x1D703}=0,\\[5.0pt] \displaystyle u_{i,i}^{\,f}=0,\\[5.0pt] \displaystyle -\mathit{K}_{r}u_{i}^{p}-\unicode[STIX]{x1D709}(u_{i}^{p}-u_{i}^{\,f})-\unicode[STIX]{x03C0}_{,i}^{p}+Rk_{i}\unicode[STIX]{x1D703}=0,\\[5.0pt] \displaystyle u_{i,i}^{p}=0,\\[5.0pt] \displaystyle R(w^{f}+w^{p})+\unicode[STIX]{x0394}\unicode[STIX]{x1D703}=0,\end{array}\right\}\end{eqnarray}$$

where the Rayleigh number is now given by

(6.2) $$\begin{eqnarray}Ra=R^{2}={\displaystyle \frac{\unicode[STIX]{x1D6FD}d^{2}\unicode[STIX]{x1D70C}_{F}g\unicode[STIX]{x1D6FC}\mathit{K}_{f}}{\unicode[STIX]{x1D707}k}},\end{eqnarray}$$

where $k=\unicode[STIX]{x1D705}_{m}/(\unicode[STIX]{x1D70C}c)_{f}$ . The key parameters in (6.1) are $\unicode[STIX]{x1D706},\unicode[STIX]{x1D709}$ and $K_{r}$ and these are defined by

(6.3a-c ) $$\begin{eqnarray}\unicode[STIX]{x1D706}={\displaystyle \frac{\tilde{\unicode[STIX]{x1D707}}\mathit{K}_{f}}{d^{2}\unicode[STIX]{x1D707}}},\quad \unicode[STIX]{x1D709}={\displaystyle \frac{\unicode[STIX]{x1D701}\mathit{K}_{f}}{\unicode[STIX]{x1D707}}},\quad K_{r}={\displaystyle \frac{K_{f}}{K_{p}}}.\end{eqnarray}$$

The coefficient $\unicode[STIX]{x1D706}$ is essentially a non-dimensional measure of the Brinkman viscosity coefficient, or alternatively may be thought of as a non-dimensional version of the ratio of the Brinkman viscosity to the dynamic viscosity of the saturating fluid. The parameter $\unicode[STIX]{x1D709}$ is a non-dimensional version of the coefficient of momentum transfer between the macro and microphases. The coefficient $K_{r}$ is the ratio of the macropermeability to the micropermeability.

The boundary conditions on the perturbation variables are

(6.4a-c ) $$\begin{eqnarray}u_{i}^{\,f}=0,\quad w^{p}=0,\quad \unicode[STIX]{x1D703}=0,\quad z=0,1.\end{eqnarray}$$

Due to the presence of the Brinkman term we require a further boundary condition on $w^{f}$ . This depends on whether we treat stress free surfaces or fixed surfaces. The fixed surface problem, which involves more intricate numerical computation of eigenvalues of a singular problem, is addressed in § 7.

Next, take curlcurl of (6.1) $_{1,3}$ and retain the third components to derive the system of equations

(6.5) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}(1+\unicode[STIX]{x1D709})\unicode[STIX]{x0394}w^{f}-\unicode[STIX]{x1D709}\unicode[STIX]{x0394}w^{p}-R\unicode[STIX]{x1D6E5}^{\ast }\unicode[STIX]{x1D703}-\unicode[STIX]{x1D706}\unicode[STIX]{x1D6E5}^{2}w^{f}=0,\\[5.0pt] (K_{r}+\unicode[STIX]{x1D709})\unicode[STIX]{x0394}w^{p}-\unicode[STIX]{x1D709}\unicode[STIX]{x0394}w^{f}-R\unicode[STIX]{x1D6E5}^{\ast }\unicode[STIX]{x1D703}=0,\\[5.0pt] R(w^{f}+w^{p})+\unicode[STIX]{x0394}\unicode[STIX]{x1D703}=0,\end{array}\right\}\end{eqnarray}$$

where $\unicode[STIX]{x1D6E5}^{\ast }=\unicode[STIX]{x2202}^{2}/\unicode[STIX]{x2202}x^{2}+\unicode[STIX]{x2202}^{2}/\unicode[STIX]{x2202}y^{2}$ . We seek a solution to these equations of form $w^{f}=W^{f}(z)h(x,y)$ where $h$ is a planform which satisfies $\unicode[STIX]{x1D6E5}^{\ast }h=-a^{2}h$ , cf. Chandrasekhar (Reference Chandrasekhar1981, pp. 43–52), where $a$ is the wavenumber, with a similar representation for $w^{p}$ and $\unicode[STIX]{x1D703}$ . The functions $W^{f},W^{p}$ and $\unicode[STIX]{x1D6E9}$ are composed of terms like $\sin \,n\unicode[STIX]{x03C0}z$ , $n=1,2,\ldots$ This results in (6.5) leading to an expression for $R^{2}$ in terms of $\unicode[STIX]{x1D6EC}_{n}=n^{2}\unicode[STIX]{x03C0}^{2}+a^{2}$ . One may differentiate this equation in $n^{2}$ and conclude $n=1$ yields the minimum value and then the critical Rayleigh number is found by minimizing

(6.6) $$\begin{eqnarray}R^{2}={\displaystyle \frac{\unicode[STIX]{x1D706}(K_{r}+\unicode[STIX]{x1D709})\unicode[STIX]{x1D6EC}^{3}+(K_{r}+K_{r}\unicode[STIX]{x1D709}+\unicode[STIX]{x1D709})\unicode[STIX]{x1D6EC}^{2}}{a^{2}(\unicode[STIX]{x1D706}\unicode[STIX]{x1D6EC}+4\unicode[STIX]{x1D709}+K_{r}+1)}},\end{eqnarray}$$

in $a$ , where $\unicode[STIX]{x1D6EC}=\unicode[STIX]{x03C0}^{2}+a^{2}$ . Results are presented below for minimization of $R^{2}$ to yield the critical values of $Ra$ and $a$ , in terms of $\unicode[STIX]{x1D709},K_{r}$ and $\unicode[STIX]{x1D706}$ .

Two limit cases which arise from (6.6) are worth noting. We observe that as $\unicode[STIX]{x1D706}\rightarrow 0$ the critical value of $Ra$ is found as

(6.7) $$\begin{eqnarray}Ra=\left({\displaystyle \frac{K_{r}+K_{r}\unicode[STIX]{x1D709}+\unicode[STIX]{x1D709}}{1+4\unicode[STIX]{x1D709}+K_{r}}}\right)4\unicode[STIX]{x03C0}^{2},\end{eqnarray}$$

as obtained in Gentile & Straughan (Reference Gentile and Straughan2017a ). Furthermore, in the formal limit $K_{r}\rightarrow \infty$ , one finds

$$\begin{eqnarray}R^{2}\sim {\displaystyle \frac{\unicode[STIX]{x1D706}\unicode[STIX]{x1D6EC}^{3}+(1+\unicode[STIX]{x1D709})\unicode[STIX]{x1D6EC}^{2}}{a^{2}}},\end{eqnarray}$$

and then as $\unicode[STIX]{x1D709}\rightarrow 0$ we obtain the single-porosity case analysed by Rees (Reference Rees2002). For the Brinkman–Darcy theory studied here we find experimental values of $K_{f}$ and $K_{r}$ which suggest $K_{r}\gg 1$ .

To understand the behaviour of the critical Rayleigh number in both free and fixed surface cases it is helpful to recollect the energy equation (5.2) and the forms for $I$ and $D$ in (5.3) and (5.4). In the heated from below isotropic case under analysis in this section the function $I$ is as defined in (5.3). Note that none of the parameters $\unicode[STIX]{x1D706},K_{r},\unicode[STIX]{x1D709}$ appear in $I$ but they are all in $D$ , which here has the form

$$\begin{eqnarray}\displaystyle D & = & \displaystyle \unicode[STIX]{x1D706}\langle u_{i,j}^{\,f}u_{i,j}^{\,f}\rangle +\langle u_{i}^{\,f}u_{i}^{\,f}\rangle +K_{r}\langle u_{i}^{p}u_{i}^{p}\rangle \nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D709}\langle (u_{i}^{\,f}-u_{i}^{p})(u_{i}^{\,f}-u_{i}^{p})\rangle +\langle \unicode[STIX]{x1D703}_{,i}\unicode[STIX]{x1D703}_{,i}\rangle .\nonumber\end{eqnarray}$$

Thus, we expect $\unicode[STIX]{x1D706},K_{r}$ and $\unicode[STIX]{x1D709}$ to each exert a stabilizing influence (in that the critical Rayleigh number increases with each parameter increasing) and this is exactly what we witness numerically for either free or fixed boundary conditions.

7 Heated below isotropic theory, two fixed surfaces

In this section we solve the linear instability problem for (6.1), but we now analyse the fixed surface problem. Thus, we solve equations (6.5) with normal modes and employ the boundary conditions

(7.1a-e ) $$\begin{eqnarray}w\,^{f}=0,\quad w^{\,f\prime }=0,\quad w^{\,p}=0,\quad \unicode[STIX]{x1D703}=0,\quad z=0,1,\end{eqnarray}$$

where $w^{\,f\prime }=\unicode[STIX]{x2202}w\,^{f}/\unicode[STIX]{x2202}z$ . To solve this system numerically we use a modified $D^{2}$ Chebyshev tau method. The Chebyshev tau method is described in Orszag (Reference Orszag1971), while the $D^{2}$ version is described in Dongarra, Straughan & Walker (Reference Dongarra, Straughan and Walker1996). Thus, we write, as before, $w\,^{f}=W^{f}(z)h(x,y)$ , with a similar representation for $w^{p}$ and $\unicode[STIX]{x1D703}$ . Let $D=\text{d}/\text{d}z$ and let $a$ be the wavenumber. Then (6.5) reduce to

(7.2) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D706}(D^{2}-a^{2})^{2}W^{f}-(1+\unicode[STIX]{x1D709})(D^{2}-a^{2})W^{p}-R\,a^{2}\unicode[STIX]{x1D6E9}=0,\\[5.0pt] \displaystyle (K_{r}+\unicode[STIX]{x1D709})(D^{2}-a^{2})W^{p}-\unicode[STIX]{x1D709}(D^{2}-a^{2})W^{f}+R\,a^{2}\unicode[STIX]{x1D6E9}=0,\\[5.0pt] \displaystyle (D^{2}-a^{2})\unicode[STIX]{x1D6E9}+R(W^{f}+W^{p})=0,\end{array}\right\}\end{eqnarray}$$

for $0<z<1$ , which are to be solved together with the boundary conditions

(7.3a-e ) $$\begin{eqnarray}W^{f}=0,\quad DW^{f}=0,\quad W^{p}=0,\quad \unicode[STIX]{x1D6E9}=0,\quad z=0,1.\end{eqnarray}$$

The modified method we employ introduces variables $\unicode[STIX]{x1D712}$ and $\unicode[STIX]{x1D713}$ by

$$\begin{eqnarray}\unicode[STIX]{x1D712}=(D^{2}-a^{2})W^{f}\quad \text{and}\quad \unicode[STIX]{x1D713}=(D^{2}-a^{2})W^{p}.\end{eqnarray}$$

Equations (7.2) are then rewritten as the system

(7.4) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle (D^{2}-a^{2})W^{f}-\unicode[STIX]{x1D712}=0,\\[5.0pt] \displaystyle \left(D^{2}-a^{2}-{\displaystyle \frac{1+\unicode[STIX]{x1D709}}{\unicode[STIX]{x1D706}}}\right)\unicode[STIX]{x1D712}+{\displaystyle \frac{\unicode[STIX]{x1D709}}{\unicode[STIX]{x1D706}}}\unicode[STIX]{x1D713}={\displaystyle \frac{R}{\unicode[STIX]{x1D706}}}\,a^{2}\unicode[STIX]{x1D6E9},\\[10.0pt] \displaystyle (D^{2}-a^{2})W^{p}-\unicode[STIX]{x1D713}=0,\\[5.0pt] \displaystyle \unicode[STIX]{x1D713}-\left({\displaystyle \frac{\unicode[STIX]{x1D709}}{K_{r}+\unicode[STIX]{x1D709}}}\right)\unicode[STIX]{x1D712}=-{\displaystyle \frac{R}{(K_{r}+\unicode[STIX]{x1D709})}}\,a^{2}\unicode[STIX]{x1D6E9},\\[10.0pt] \displaystyle (D^{2}-a^{2})\unicode[STIX]{x1D6E9}=-R(W^{f}+W^{p}).\end{array}\right\}\end{eqnarray}$$

In terms of the Chebyshev polynomials $T_{n}(z)$ we now write

$$\begin{eqnarray}\displaystyle & \displaystyle W^{f}=\mathop{\sum }_{n=0}^{N}W_{n}^{f}T_{n}(z),\quad \unicode[STIX]{x1D712}=\mathop{\sum }_{n=0}^{N}\unicode[STIX]{x1D712}_{n}T_{n}(z), & \displaystyle \nonumber\\ \displaystyle & \displaystyle W^{p}=\mathop{\sum }_{n=0}^{N}W_{n}^{p}T_{n}(z),\quad \unicode[STIX]{x1D713}=\mathop{\sum }_{n=0}^{N}\unicode[STIX]{x1D713}_{n}T_{n}(z), & \displaystyle \nonumber\\ \displaystyle & \displaystyle \unicode[STIX]{x1D6E9}=\mathop{\sum }_{n=0}^{N}\unicode[STIX]{x1D6E9}_{n}T_{n}(z). & \displaystyle \nonumber\end{eqnarray}$$

In this way equations (7.4) reduce to solving the generalized eigenvalue problem

(7.5) $$\begin{eqnarray}\unicode[STIX]{x1D63C}\unicode[STIX]{x1D66D}=R\,\unicode[STIX]{x1D63D}\unicode[STIX]{x1D66D},\end{eqnarray}$$

where the matrices $\unicode[STIX]{x1D63C}$ and $\unicode[STIX]{x1D63D}$ are given by

$$\begin{eqnarray}\unicode[STIX]{x1D63C}=\left(\begin{array}{@{}ccccc@{}}D^{2}-a^{2}\unicode[STIX]{x1D644} & -I & 0 & 0 & 0\\ 0 & D^{2}-(a^{2}+(1+\unicode[STIX]{x1D709})\unicode[STIX]{x1D706}^{-1})I & 0 & (\unicode[STIX]{x1D709}/\unicode[STIX]{x1D706})I & 0\\ 0 & 0 & D^{2}-a^{2}I & -I & 0\\ 0 & -(\unicode[STIX]{x1D709}/(K_{r}+\unicode[STIX]{x1D709}))I & 0 & I & 0\\ 0 & 0 & 0 & 0 & D^{2}-a^{2}I\end{array}\right)\end{eqnarray}$$

and

$$\begin{eqnarray}\unicode[STIX]{x1D63D}=\left(\begin{array}{@{}ccccc@{}}0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & (a^{2}/\unicode[STIX]{x1D706})I\\ 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & -(a^{2}/(K_{r}+\unicode[STIX]{x1D709}))I\\ -I & 0 & -I & 0 & 0\end{array}\right),\end{eqnarray}$$

where $\unicode[STIX]{x1D644}$ is the $n\times n$ identity matrix and $0$ is the $n\times n$ zero matrix. The last two rows of each of the block matrices in rows 1, 2, 3 and 5 contain the boundary conditions. Details of how to do this are given in Dongarra et al. (Reference Dongarra, Straughan and Walker1996, pp. 406, 407). We make use of the fact that $T_{n}(\pm 1)=(\pm 1)^{n}$ and $T_{n}^{\prime }=(\pm 1)^{n+1}n^{2}$ . The last two rows of blocks 1, 1 and 2, 1 of $A$ contain the boundary conditions for $W_{f}$ while the last two rows of block 3, 3 contains the boundary conditions for $W_{p}$ . The last two rows of block 5, 5 contains the boundary conditions for $\unicode[STIX]{x1D6E9}$ . The matrices in block rows 4 contain the complete $n\times n$ identity matrix for both $A$ and $B$ .

The resulting finite dimensional generalized eigenvalue problem (7.5) is solved using the QZ algorithm of Moler & Stewart (Reference Moler and Stewart1971). We find this modified $D^{2}$ Chebyshev tau method works well and yields accurate results, although care has to be taken to avoid spurious eigenvalues.

8 Numerical results

We present numerical results for computations involving minimizing $Ra$ in (6.6) in $a$ for two free surfaces and likewise minimizing in $a$ for the solution from (7.4) or (7.5) when the surfaces are fixed. In both cases the Rayleigh number is a function of the wavenumber, $a$ , but also of the parameters $\unicode[STIX]{x1D706},\unicode[STIX]{x1D709}$ and $K_{r}$ defined in (6.3). We firstly assess what range of values these parameters may take in real bidisperse porous materials. In order to do this we make use of the articles of Givler & Altobelli (Reference Givler and Altobelli1994), and of Hooman et al. (Reference Hooman, Sauret and Dahari2015), and we employ the Carman–Kozeny relation for the permeability as given by Chen (Reference Chen1990) or Nield (Reference Nield2000).

In their analysis of an open cell rigid foam Givler & Altobelli (Reference Givler and Altobelli1994) deduce experimentally that

(8.1) $$\begin{eqnarray}5.1\leqslant {\displaystyle \frac{\tilde{\unicode[STIX]{x1D707}}}{\unicode[STIX]{x1D707}}}\leqslant 10.9\end{eqnarray}$$

and they also calculate $K_{f}=3.3\times 10^{-7}~\text{m}^{2}$ with an error

(8.2) $$\begin{eqnarray}0.0027.1\leqslant K_{f}\leqslant 0.0042~\text{cm}^{2}.\end{eqnarray}$$

For the momentum transfer coefficient, $\unicode[STIX]{x1D701}$ , we have found only one source of information and this is Hooman et al. (Reference Hooman, Sauret and Dahari2015) who report for a plate–fin heat exchanger

(8.3) $$\begin{eqnarray}\unicode[STIX]{x1D701}=63.3~\text{Pa}~\text{s}~\text{m}^{-2}.\end{eqnarray}$$

Hooman et al. (Reference Hooman, Sauret and Dahari2015) also supply values for the macro and micropermeabilities depending on the spacing between the plates and they report

(8.4) $$\begin{eqnarray}K_{f}=2.13\times 10^{-7}~\text{m}^{2},\end{eqnarray}$$

with (Hooman et al. (Reference Hooman, Sauret and Dahari2015) table 2, different values corresponding to different spacing)

(8.5) $$\begin{eqnarray}K_{r}={\displaystyle \frac{K_{f}}{K_{p}}}=701.75,456.14,263.16.\end{eqnarray}$$

Observe that the value of $K_{f}$ in (8.4) is in keeping with those of Givler & Altobelli (Reference Givler and Altobelli1994) in (8.2).

If we use the Carman–Kozeny relation of Chen (Reference Chen1990), then

(8.6) $$\begin{eqnarray}K_{f}={\displaystyle \frac{d_{f}^{2}}{172.8}}{\displaystyle \frac{\unicode[STIX]{x1D719}^{3}}{(1-\unicode[STIX]{x1D719})^{2}}},\end{eqnarray}$$

where $d_{f}$ is the diameter of the spheres comprising the porous medium and $\unicode[STIX]{x1D719}$ is the porosity. For a porosity of $\unicode[STIX]{x1D719}=0.8$ and 5 mm glass beads (8.6) yields a value for $K_{f}$ of

(8.7) $$\begin{eqnarray}K_{f}=1.85\times 10^{-6}~\text{m}^{2},\end{eqnarray}$$

which is consistent with (8.2) or (8.4).

If we use relation (8.6) to calculate $K_{r}$ then

(8.8) $$\begin{eqnarray}K_{r}={\displaystyle \frac{K_{f}}{K_{p}}}={\displaystyle \frac{d_{f}^{2}}{d_{p}^{2}}}{\displaystyle \frac{\unicode[STIX]{x1D719}^{3}}{(1-\unicode[STIX]{x1D719})^{2}}}{\displaystyle \frac{(1-\unicode[STIX]{x1D716})^{2}}{\unicode[STIX]{x1D716}^{3}}},\end{eqnarray}$$

where $d_{f}$ and $d_{p}$ denote the diameter of the spheres in the macro and microphases. With $\unicode[STIX]{x1D719}$ and $\unicode[STIX]{x1D716}$ the permeabilities, then for $d_{f}/d_{p}=10$ , $\unicode[STIX]{x1D719}=0.8$ and $\unicode[STIX]{x1D716}=0.3$ we find $K_{r}=2322$ .

We may also calculate the Brinkman parameter $\unicode[STIX]{x1D706}$ using equation (8.6) then

$$\begin{eqnarray}\unicode[STIX]{x1D706}={\displaystyle \frac{\tilde{\unicode[STIX]{x1D707}}}{\unicode[STIX]{x1D707}}}\,{\displaystyle \frac{1}{172.8}}\,{\displaystyle \frac{d_{f}^{2}}{d^{2}}}\,{\displaystyle \frac{\unicode[STIX]{x1D719}^{3}}{(1-\unicode[STIX]{x1D719})^{2}}},\end{eqnarray}$$

where $d$ is the depth of the convection layer. For a layer depth of 3 cm, with 3 mm glass beads, and assuming $\tilde{\unicode[STIX]{x1D707}}/\unicode[STIX]{x1D707}=7.5$ as in Givler & Altobelli (Reference Givler and Altobelli1994), we find

$$\begin{eqnarray}\unicode[STIX]{x1D706}\approx 5.56\times 10^{-3}.\end{eqnarray}$$

If we employ the Givler & Altobelli (Reference Givler and Altobelli1994) of $K_{f}=3.3\times 10^{-7}~\text{m}^{2}$ then for a 3 cm layer we obtain

$$\begin{eqnarray}\unicode[STIX]{x1D706}\approx 2.75\times 10^{-3}.\end{eqnarray}$$

If we take $d$ to be 1 cm, take $K_{f}=4.2\times 10^{-7}$ and $\tilde{\unicode[STIX]{x1D707}}/\unicode[STIX]{x1D707}=10.9$ , which are in the range of the Givler & Altobelli (Reference Givler and Altobelli1994) values, then we find

$$\begin{eqnarray}\unicode[STIX]{x1D706}={\displaystyle \frac{\tilde{\unicode[STIX]{x1D707}}}{\unicode[STIX]{x1D707}}}\,{\displaystyle \frac{K_{f}}{d^{2}}}\approx 4.578\times 10^{-2}.\end{eqnarray}$$

One may further employ equation (8.8) to see that if $d_{f}/d_{p}=5$ and we take $\unicode[STIX]{x1D719}=0.6$ , $\unicode[STIX]{x1D716}=0.6$ then $K_{r}=25$ whereas for $d_{f}/d_{p}=4$ , $\unicode[STIX]{x1D719}=0.8$ , $\unicode[STIX]{x1D716}=0.6$ then $K_{r}=151.7$ .

To calculate $\unicode[STIX]{x1D709}$ we use the relation

(8.9) $$\begin{eqnarray}\unicode[STIX]{x1D709}={\displaystyle \frac{\unicode[STIX]{x1D701}K_{f}}{\unicode[STIX]{x1D707}}},\end{eqnarray}$$

where $\unicode[STIX]{x1D709}$ is the momentum transfer coefficient which from Hooman et al. (Reference Hooman, Sauret and Dahari2015) is given by (8.3). We take the fluid to be water at $25\,^{\circ }\text{C}$ and then from the website engineersedge.com $\unicode[STIX]{x1D707}=8.9\times 10^{-4}~\text{Pa}~\text{s}$ . We have derived above four typical values for the macropermeability, $K_{f}$ , and these are

$$\begin{eqnarray}K_{f}=1.85\times 10^{-6}~\text{m}^{2},\;2.13\times 10^{-7}~\text{m}^{2},\;3.3\times 10^{-7}~\text{m}^{2},\;4.2\times 10^{-7}~\text{m}^{2},\end{eqnarray}$$

which, using (8.9), lead to values of $\unicode[STIX]{x1D709}$ as

(8.10) $$\begin{eqnarray}\unicode[STIX]{x1D709}=1.515\times 10^{-2},2.347\times 10^{-2},2.987\times 10^{-2},1.316\times 10^{-1}.\end{eqnarray}$$

The range of $\unicode[STIX]{x1D706}$ and $K_{r}$ values we have found are

(8.11) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D706}=2.75\times 10^{-3},5.56\times 10^{-3},4.578\times 10^{-2},\\[5.0pt] \displaystyle K_{r}=25,151.7,263.16,456.14,701.75,2322.\end{array}\right\}\end{eqnarray}$$

In our computations we mostly employ these values to determine the behaviour of the critical Rayleigh number, $Ra$ , and the critical wavenumber, $a$ , as functions of the non-dimensional parameters $\unicode[STIX]{x1D706},\unicode[STIX]{x1D709}$ and $K_{r}$ .

In both the fixed and free surface boundary condition cases the Rayleigh number is found to increase with variation in $K_{r},\unicode[STIX]{x1D706}$ or $\unicode[STIX]{x1D709}$ . The variation in $Ra$ with $K_{r}$ is shown in table 5 in the fixed case when $\unicode[STIX]{x1D706}$ and $\unicode[STIX]{x1D709}$ take possible realistic values. To interpret these we recall the definition of $Ra$ in (6.2). Since $K_{r}=K_{f}/K_{p}$ we may think of $K_{f}$ as fixed and then as $K_{r}$ increases $K_{p}$ decreases. This means the micropermeability decreases and so the fluid moves less easily in the micropores. In this case we expect convective motion to be less easy and so the system is more stable as is witnessed by increasing $Ra$ . The analogous variation in the critical value of the wavenumber, $a$ depends strongly on the boundary conditions. Since the wavenumber is inversely proportional to the cell width (aspect ratio) this means increasing $a$ leads to wider convection cells, whereas $a$ decreasing means narrower cells. In table 5 we see that, for fixed boundary conditions, $a$ increases as $K_{r}$ increases, although the increase in $a$ is not large. However, the free boundary case is very different, as seen in table 6. When $\unicode[STIX]{x1D709}=0.1316$ and $\unicode[STIX]{x1D706}=0.04578$ we find $a$ increases from the value 2.99266 when $K_{r}=0$ to a maximum of 2.99466 when $K_{r}=0.037$ and thereafter decreases to the value $a=2.64953$ when $K_{r}=10^{4}$ . For the case of $\unicode[STIX]{x1D709}=0.01515,\unicode[STIX]{x1D706}=2.75\times 10^{-3}$ , $a$ decreases from the value 3.13720 at $K_{r}=0$ to $a=3.06643$ when $K_{r}=10^{4}$ . The variation of $a^{2}$ against $K_{r}$ in the free case with $\unicode[STIX]{x1D706}=1,\unicode[STIX]{x1D709}=1$ is shown in figure 3.

Figure 1. Graph of $a^{2}$ versus $Kr$ with $Kr$ varying in the range 0 to 10, two free surfaces. Here, $\unicode[STIX]{x1D706}=1$ and $\unicode[STIX]{x1D709}=1$ . The values at the end points as shown are $(a^{2},Kr)=(8.491,0)$ and $(a^{2},Kr)=(7.220,10)$ . The maximum value is approximately at $(a^{2},Kr)=(8.502,0.21)$ .

Figure 2. Graph of $Ra$ versus $\unicode[STIX]{x1D706}$ with $\unicode[STIX]{x1D706}$ varying in the range 0 to 10, two free surfaces. Here $Kr=1$ and $\unicode[STIX]{x1D709}=1$ . The values at the endpoints as shown are $(Ra,\unicode[STIX]{x1D706})=(19.73921,0)$ and $(Ra2,\unicode[STIX]{x1D706})=(77.20063,10)$ .

Figure 3. Graph of $a^{2}$ versus $\unicode[STIX]{x1D709}$ with $\unicode[STIX]{x1D709}$ varying in the range 0 to 10, two free surfaces. Here $Kr=1$ and $\unicode[STIX]{x1D706}=1$ . The values at the endpoints as shown are $(a^{2},\unicode[STIX]{x1D709})=(9.448,0)$ and $(a^{2},\unicode[STIX]{x1D709})=(6.101,10)$ .

Figure 2 shows that $Ra$ increases rapidly as the Brinkman coefficient $\unicode[STIX]{x1D706}$ increases for $\unicode[STIX]{x1D706}$ small, and then the growth is less rapid for larger $\unicode[STIX]{x1D706}$ . This is in agreement with the case of a single-porosity Brinkman material. The variation in the critical wavenumber $a$ in figures 4–6 is interesting. In our notation, in the single-porosity case replacing $w\,^{f}$ and $w^{p}$ by $w$ , instead of (6.5) one has, cf. Rees (Reference Rees2002),

$$\begin{eqnarray}\displaystyle & \displaystyle -\unicode[STIX]{x1D706}\unicode[STIX]{x1D6E5}^{2}w+\unicode[STIX]{x0394}w-R\unicode[STIX]{x1D6E5}^{\ast }\unicode[STIX]{x1D703}=0, & \displaystyle \nonumber\\ \displaystyle & \displaystyle \unicode[STIX]{x0394}\unicode[STIX]{x1D703}+Rw=0, & \displaystyle \nonumber\end{eqnarray}$$

and instead of (6.6) one finds

(8.12) $$\begin{eqnarray}R^{2}={\displaystyle \frac{\unicode[STIX]{x1D706}\unicode[STIX]{x1D6EC}^{3}+\unicode[STIX]{x1D6EC}^{2}}{a^{2}}}.\end{eqnarray}$$

(We stress that Rees (Reference Rees2002) analyses the fixed surface problem.) The presence of the $\unicode[STIX]{x1D706}\unicode[STIX]{x1D6EC}$ term in the denominator of (6.6) changes the $a^{2}$ versus $\unicode[STIX]{x1D706}$ behaviour significantly between the one porosity and the bidispersive situation. In the one-porosity case we find for two stress free surfaces that $a^{2}$ decreases from $\unicode[STIX]{x03C0}^{2}$ to $\unicode[STIX]{x03C0}^{2}/2$ as $\unicode[STIX]{x1D706}$ increases from 0. In the one-porosity case the critical wavenumber may be obtained analytically and we find

$$\begin{eqnarray}a_{cr}^{2}=\{-(\unicode[STIX]{x1D706}\unicode[STIX]{x03C0}^{2}+1)+[(\unicode[STIX]{x1D706}\unicode[STIX]{x03C0}^{2}+1)^{2}+8\unicode[STIX]{x1D706}\unicode[STIX]{x03C0}^{2}(1+\unicode[STIX]{x1D706}\unicode[STIX]{x03C0}^{2})]^{1/2}\}/4\unicode[STIX]{x1D706}.\end{eqnarray}$$

Using Newton’s binomial expansion we find for $\unicode[STIX]{x1D706}$ small,

(8.13) $$\begin{eqnarray}a_{cr}^{2}=\unicode[STIX]{x03C0}^{2}-2\unicode[STIX]{x03C0}^{4}\unicode[STIX]{x1D706}+10\unicode[STIX]{x03C0}^{6}\unicode[STIX]{x1D706}^{2}+O(\unicode[STIX]{x1D706}^{3}),\end{eqnarray}$$

which confirms the decrease for $\unicode[STIX]{x1D706}\ll 1$ .

Figure 4. Graph of $a$ versus $\log _{10}\unicode[STIX]{x1D706}$ , for $\unicode[STIX]{x1D709}=1,K_{r}=1$ . The solid curve indicates two free surfaces whereas the solid dots are for two fixed surfaces. The minimum on the solid curve is where $a=2.7207,\unicode[STIX]{x1D706}=0.17$ , and the maximum value displayed for fixed surfaces is $a=3.2088,\unicode[STIX]{x1D706}=0.01$ .

Figure 5. Graph of $a$ versus $\log _{10}\unicode[STIX]{x1D706}$ , for $\unicode[STIX]{x1D709}=0.1316,K_{r}=2322$ . The solid curve indicates two free surfaces whereas the solid dots are for two fixed surfaces. The minimum on the solid curve is where $a=2.24586,\unicode[STIX]{x1D706}=3.46$ , and the maximum value displayed for fixed surfaces is $a=3.239,$ $\unicode[STIX]{x1D706}\in [0.005,0.006]$ .

Figure 6. Graph of $a$ versus $\log _{10}\unicode[STIX]{x1D706}$ , for $\unicode[STIX]{x1D709}=0.01515,K_{r}=25$ . The solid curve indicates two free surfaces whereas the solid dots are for two fixed surfaces. The minimum on the solid curve is where $a=2.43084,\unicode[STIX]{x1D706}=0.33$ , and the maximum value displayed for fixed surfaces is $a=3.234,$ $\unicode[STIX]{x1D706}\in [0.004,0.006]$ .

In the bidisperse case figures 4–6 show for two free surfaces a clear decrease then increase as $\unicode[STIX]{x1D706}$ increases for small $\unicode[STIX]{x1D706}$ . The initial decrease may be verified by an asymptotic analysis from (6.6). If we form $\unicode[STIX]{x2202}R^{2}/\unicode[STIX]{x2202}a^{2}$ from (6.6) and equate to zero and then solve the resulting equation by expanding $a^{2}$ as a power series in $\unicode[STIX]{x1D706}$ then for $0<\unicode[STIX]{x1D706}\ll 1$ we find

$$\begin{eqnarray}a^{2}=\unicode[STIX]{x03C0}^{2}-X_{1}\unicode[STIX]{x1D706}+O(\unicode[STIX]{x1D706}^{2}),\quad \text{where }X_{1}=2\unicode[STIX]{x03C0}^{6}\left[{\displaystyle \frac{K_{r}^{2}+4K_{r}\unicode[STIX]{x1D709}+4\unicode[STIX]{x1D709}^{2}}{(K_{r}+K_{r}\unicode[STIX]{x1D709}+\unicode[STIX]{x1D709})(1+K_{r}+4\unicode[STIX]{x1D709})}}\right].\end{eqnarray}$$

In fact, when $K_{r}=1,\unicode[STIX]{x1D709}=1$ the minimum value is $(a^{2},\unicode[STIX]{x1D706})=(7.402,0.17)$ whereas for $K_{r}=5,\unicode[STIX]{x1D709}=0.01$ the minimum is at $(a^{2},\unicode[STIX]{x1D706})=(6.953,0.15)$ . Thus the presence of micropores is changing the behaviour of the cell shape. When micropores are present and $\unicode[STIX]{x1D706}$ is small the cell shape increases in width (aspect ratio) as $\unicode[STIX]{x1D706}$ increases, reaching a maximum and then decreasing again. The wavenumber variation for $\unicode[STIX]{x1D706}$ small is, however, very different from the free–free case. We see that in tables 3 and 4 when the surfaces are fixed then the wavenumber has a maximum in $\unicode[STIX]{x1D706}$ . This is exactly the opposite behaviour to the free–free situation. This behaviour is, however, in agreement with the findings of Rees (Reference Rees2002) in the single-porosity case. He also predicts a rise then decrease in $a$ as $\unicode[STIX]{x1D706}$ decreases. We have checked the behaviour of the critical value of $a$ in the free–free single-porosity case and it is as in figures 4–6. We have computed the critical Rayleigh number and wavenumber for various other combinations of possible realistic values of parameter values and we always see the same type of behaviour as $\unicode[STIX]{x1D706}$ varies.

Table 2. Critical values of $Ra$ and $a$ for quoted values of $\unicode[STIX]{x1D706}$ , two fixed surfaces. Columns 2 and 3 are for $K_{r}=1$ , $\unicode[STIX]{x1D709}=1$ , while columns 4 and 5 are for $K_{r}=5$ , $\unicode[STIX]{x1D709}=0.01$ .

Table 3. Critical values of $Ra$ and $a$ for quoted values of $\unicode[STIX]{x1D706}$ , two fixed surfaces. Columns 2 and 3 are for $K_{r}=1$ , $\unicode[STIX]{x1D709}=1$ , while columns 4 and 5 are for $K_{r}=5$ , $\unicode[STIX]{x1D709}=0.01$ .

Table 4. Critical values of $Ra$ and $a$ for quoted values of $\unicode[STIX]{x1D706}$ , two fixed surfaces. Columns 2 and 3 are for $K_{r}=2322$ , $\unicode[STIX]{x1D709}=0.1316$ , while columns 4 and 5 are for $K_{r}=25$ , $\unicode[STIX]{x1D709}=0.01515$ .

Table 5. Critical values of $Ra$ and $a$ for quoted values of $Kr$ , two fixed surfaces. Columns 2 and 3 are for $\unicode[STIX]{x1D706}=0.04578$ , $\unicode[STIX]{x1D709}=0.1316$ , while columns 4 and 5 are for $\unicode[STIX]{x1D706}=0.00275$ , $\unicode[STIX]{x1D709}=0.01515$ .

Table 6. Critical values of the wavenumber $a$ for quoted values of $Kr$ , two free surfaces. Columns 2 and 3 are for $\unicode[STIX]{x1D706}=0.04578$ , $\unicode[STIX]{x1D709}=0.1316$ , while columns 4 and 5 are for $\unicode[STIX]{x1D706}=0.00275$ , $\unicode[STIX]{x1D709}=0.01515$ .

In figure 4 the maximum and minimum values for the fixed and free cases have reasonably close values of $\unicode[STIX]{x1D706}$ , but these are when $K_{r}=1,\unicode[STIX]{x1D709}=1$ . For $K_{r}$ and $\unicode[STIX]{x1D709}$ values which we have calculated for possible real bidisperse materials figures 5 and 6 show that the maximum in $\unicode[STIX]{x1D706}$ for the fixed surface problem is in the range of realistic values. The corresponding minima in $\unicode[STIX]{x1D706}$ for the free surface situation are for larger values of  $\unicode[STIX]{x1D706}$ .

The variations of $Ra$ and $a$ in $\unicode[STIX]{x1D709}$ are given in tables 7 and 8 with potentially realistic values of parameters. For fixed surfaces the situation of the behaviour of $a$ is not clear. In table 7, for $\unicode[STIX]{x1D706}=4.578\times 10^{-2}$ and $K_{r}=25$ or 2322 the critical wavenumber $a$ always increases as $\unicode[STIX]{x1D709}$ increases. However, from table 7, for $\unicode[STIX]{x1D706}=2.75\times 10^{-3}$ and $K_{r}=25$ or 2322, $a$ always decreases. For the free surface case $a$ is found to increase in all four cases as seen in table 9.

Table 7. Critical values of $Ra$ and $a$ for quoted values of $\unicode[STIX]{x1D709}$ , two fixed surfaces. Columns 2 and 3 are for $\unicode[STIX]{x1D706}=0.04578$ , $K_{r}=25$ , while columns 4 and 5 are for $\unicode[STIX]{x1D706}=0.04578$ , $Kr=2322$ .

Table 8. Critical values of $Ra$ and $a$ for quoted values of $\unicode[STIX]{x1D709}$ , two fixed surfaces. Columns 2 and 3 are for $\unicode[STIX]{x1D706}=0.00275$ , $K_{r}=25$ , while columns 4 and 5 are for $\unicode[STIX]{x1D706}=0.00275$ , $Kr=2322$ .

Table 9. Critical values of the wavenumber $a$ for quoted values of $\unicode[STIX]{x1D709}$ , two free surfaces. Column 2 is for $\unicode[STIX]{x1D706}=0.04578$ , $K_{r}=25$ , column 3 is for $\unicode[STIX]{x1D706}=0.04578$ , $Kr=2322$ , column 4 is for $\unicode[STIX]{x1D706}=0.00275$ , $K_{r}=25$ and column 5 is for $\unicode[STIX]{x1D706}=0.00275$ , $Kr=2322$ .