2 Problem formulation

Consider an infinite horizontal layer of porous material saturated with compressible PIP, heated from below and confined between two surfaces at $z=0$ and $z=d$. Let $\rho _n$ denote the density of the neutral components of PIP, $\rho _i$ the density of the ionized components of PIP and $K_1$ the permeability of the porous medium. The gravitational force, $\boldsymbol {g}=( 0,0,-g )$, acts downward. A geometrical representation of the problem is illustrated in figure 1.

Figure 1. Geometrical representation of the problem.

We assume that the solid phase of the porous medium and the ionized components of the PIP are not in LTE. The temperature of the solid phase, $T_s$, and the plasma phase, $T_i$, are maintained constant at the surfaces $z=0$ and $z=d$

(2.1a,b)\begin{equation} T_s=T_i=T_0, \quad \text{at} \ z=0 \quad \text{and} \quad T_s=T_i=T_d, \quad \text{at} \ z=d, \end{equation}

where $T_0>T_d$. A two-field model is employed to separately represent the temperature fields of the solid phase of the porous medium and the ionized components of the PIP in the energy equation.

In this analysis, the PIP is treated under the continuum hypothesis, behaving as a continuous fluid. Effects of LTNE, pressure, gravity and medium permeability on the neutral components of the PIP are considered negligibly small and are thus neglected.

Assuming adherence to the Boussinesq approximation (Spiegel & Veronis Reference Spiegel and Veronis1960), the governing equations are (cf. Sharma & Sharma Reference Sharma and Sharma1978; Kuznetsov Reference Kuznetsov1998; Straughan Reference Straughan2006; Nield & Bejan Reference Nield and Bejan2013)

(2.2)\begin{gather}\frac{1}{\epsilon} \frac{\partial \boldsymbol{q}}{\partial t} ={-}\frac{1}{{\rho_{m}}}\boldsymbol{\nabla} p + \boldsymbol{g} \left[ 1-\alpha_m (T-T_{m})+K_{m}\left(p-p_{m}\right)\right] \nonumber\\ -\frac{1}{\rho_{m}}\left[ \frac{\mu }{{{K}_{1}}}-\tilde{\mu }{{\nabla}^{2}} \right]\boldsymbol{q} + \frac{{{\rho }_{d}}{{\nu }_{c}}}{{{\rho }_{m}\epsilon}}\left( \boldsymbol{q}_{d}-\boldsymbol{q} \right), \end{gather}

(2.3)\begin{gather} \epsilon\frac{\partial\rho}{\partial t}+\left(\boldsymbol{q}\boldsymbol{\cdot}\boldsymbol{\nabla}\right)\rho+\rho\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{q}=0, \end{gather}

(2.4)\begin{gather} \epsilon \left(\rho c\right)_i \frac{\partial T_i}{\partial t} + \left(\rho c\right)_i \left( \boldsymbol{q}\boldsymbol{\cdot}\boldsymbol{\nabla}\right)T_i+ p \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{q}= \epsilon k_i \nabla^2 T_i +h\left(T_s-T_i\right), \end{gather}

(2.5)\begin{gather}\left(1-\epsilon\right) \left(\rho c\right)_s \frac{\partial T_s}{\partial t} = \left(1-\epsilon\right) k_s \nabla^2 T_s - h\left(T_s-T_i\right), \end{gather}

(2.6)\begin{gather}\frac{\partial \boldsymbol{q}_{d}}{\partial t}+\frac{1}{\epsilon}(\boldsymbol{q}_{d} \boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{q}_{d} ={-}\nu_{c}(\boldsymbol{q}_{d}-\boldsymbol{q}). \end{gather}

Here, $\epsilon$ is the porosity, $\boldsymbol {q}=(u, v, w)$ and $\boldsymbol {q}_{d}=(\ell, r, s)$ are the velocities of the ionized and neutral components of PIP, respectively, $t$ is time, $p$ is pressure, $\rho _m$, $T_m$ and $p_m$ are the constant space averages of density, temperature and pressure, respectively, $\tilde {\mu }$ is the effective viscosity of PIP, $\mu$ is the dynamic viscosity of PIP, $\nu _c$ is the collisional frequency between components of PIP, $k_i$ is the thermal conductivity of ionized components of PIP, $h$ is the inter-phase heat transfer coefficient, $c_i$ and $c_s$ are the specific heats of the ionized components of the PIP and the solid phase of the porous medium, respectively, and $k_s$ is the thermal conductivity of the solid phase of the porous medium. In (2.4) and (2.5), $(\rho c)_i = \rho _i c_i$ and $(\rho c)_s = \rho _s c_s$.

Additionally, $\alpha _m$ and $K_m$ are defined as

(2.7a,b)\begin{equation} \alpha_m={-}\left(\frac{1}{\rho}\frac{\partial\rho}{\partial t}\right)_m \left(=\alpha,say\right), \quad K_m=\left(\frac{1}{\rho}\frac{\partial\rho}{\partial p}\right)_m, \end{equation}

where $\alpha$ is the coefficient of thermal expansion.

3 Basic state and non-dimensionalized perturbed equations

The steady basic state (Sharma & Sunil Reference Sharma and Sunil1995) is defined as follows:

(3.1a-e)\begin{equation} \boldsymbol{q}={\boldsymbol{q}_{b}}=\textbf{0}, \quad p={{p}_{b}}\left( z \right), \quad T_i=T_{{i}_{b}}\left( z \right), \quad T_s=T_{{s}_{b}}\left( z \right) \quad \text{and} \quad {\boldsymbol{q}_{d}}={\boldsymbol{q}_{d}}_{_{b}}=\textbf{0}, \end{equation}

where

(3.2)\begin{equation} \left.\begin{gathered} \frac{1}{{\rho_{m}}}\boldsymbol{\nabla} {{p}_{b}}\left( z \right) = \boldsymbol{g} \left[ 1-\alpha_m (T-T_{m})+K_{m}\left(p-p_{m}\right)\right],\\ T_{{i}_{b}}\left( z \right)={-}\beta z+T_m=T_{{s}_{b}}\left( z \right). \end{gathered}\right\} \end{equation}

Here, the subscript ‘$b$’ denotes the basic state and $\beta (=(T_0-T_d)/d)$ is the temperature gradient.

Consider small disturbances in the form of

(3.3a-e)\begin{equation} \boldsymbol{{q}'}=\left(u',v',w'\right), \quad p', \quad {\boldsymbol{q}'_{d}}=\left(\ell',r',s'\right), \quad \theta, \quad \text{and} \quad \psi, \end{equation}

representing perturbations in $\boldsymbol {q}$, $p$, $\boldsymbol {q}^{}_{d}$, $T_i$ and $T_s$, respectively. The nonlinear perturbation equations are

(3.4)\begin{gather} \frac{1}{\epsilon} \frac{\partial \boldsymbol{{q}'}}{\partial t} ={-}\frac{1}{{\rho_{m}}}\boldsymbol{\nabla} p' - \boldsymbol{g} \alpha \theta -\frac{1}{\rho_{m}}\left[ \frac{\mu }{{{K}_{1}}}-\tilde{\mu }{{\nabla}^{2}} \right]\boldsymbol{q'} + \frac{{{\rho }_{d}}{{\nu }_{c}}}{{{\rho }_{m}\epsilon}}\left( {\boldsymbol{q}'_{d}}-\boldsymbol{{q}'} \right), \end{gather}

(3.5)\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{q'}=0, \end{gather}

(3.6)\begin{gather} \epsilon \left(\rho c\right)_i \frac{\partial \theta}{\partial t} + \left(\rho c\right)_i \left( \boldsymbol{q'}\boldsymbol{\cdot}\boldsymbol{\nabla}\right)\theta= \epsilon k_i \nabla^2 \theta +h\left(\psi-\theta\right)+\left(\rho c\right)_i \left(\beta -\frac{g}{c_i}\right) w', \end{gather}

(3.7)\begin{gather}\left(1-\epsilon\right) \left(\rho c\right)_s \frac{\partial \psi}{\partial t} = \left(1-\epsilon\right) k_s \nabla^2 \psi - h\left(\psi-\theta\right), \end{gather}

(3.8)\begin{gather} \frac{\partial {\boldsymbol{q}'_{d}}}{\partial t}+\frac{1}{\epsilon}\left({\boldsymbol{q}'_{d}} \boldsymbol{\cdot}\boldsymbol{\nabla}\right) {\boldsymbol{q}'_{d}} ={-}\nu_{c}\left({\boldsymbol{q}'_{d}}-\boldsymbol{{q}'}\right). \end{gather}

The perturbed equations (3.4)–(3.8) are non-dimensionalized using the following scales:

(3.9)\begin{equation} \left.\begin{gathered} t^*=\frac{k_i}{d^2 \left( \rho c \right)_i}t, \quad \boldsymbol{q^*}=\frac{\left(\rho c\right)_i d}{\epsilon k_i} \boldsymbol{{q'}}, \quad \boldsymbol{q}^{*}_{d}=\frac{\left(\rho c\right)_i d}{\epsilon k_i}{\boldsymbol{q}'_{d}}, \quad z^*=\frac{z}{d}, \\ p^*=\frac{K_1\left(\rho c\right)_i}{\mu \epsilon k_i}p^{*}, \quad \theta^*=\frac{\sqrt{Ra}}{\beta d}\theta, \quad \psi^*=\frac{\sqrt{Ra}}{\beta d}\psi. \end{gathered}\right\} \end{equation}

As a result, the non-dimensionalized perturbed equations (after removing the asterisk) are given as

(3.10)\begin{gather} \frac{1}{Va} \frac{\partial \boldsymbol{q}}{\partial t} ={-}\boldsymbol{\nabla} p + \sqrt{Ra}\theta \hat{\boldsymbol{k}} - \boldsymbol{q} + \widetilde{Da}\nabla^2 \boldsymbol{q} + \frac{FLPr}{Va}\left( \boldsymbol{q}_d - \boldsymbol{q}\right), \end{gather}

(3.11)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{q}=0, \end{gather}

(3.12)\begin{gather}\frac{\partial \theta}{\partial t}+ \left(\boldsymbol{q}\boldsymbol{\cdot}\boldsymbol{\nabla} \right)\theta= \nabla^2\theta + \mathcal{H}\left( \psi-\theta \right) + \sqrt{Ra}\left(1-\frac{1}{G}\right)w, \end{gather}

(3.13)\begin{gather} \mathcal{A} \frac{\partial \psi}{\partial t}= \nabla^2\psi-\mathcal{H}\gamma \left( \psi-\theta \right), \end{gather}

(3.14)\begin{gather}\frac{\partial \boldsymbol{q}_d}{\partial t}+ \left( \boldsymbol{q}_d\boldsymbol{\cdot}\boldsymbol{\nabla} \right)\boldsymbol{q}_d={-}L Pr \left( \boldsymbol{q}_d-\boldsymbol{q} \right). \end{gather}

Here, $Va=(\rho c)_i d^2 \mu \epsilon / k_i \rho _m K_1$ is the Vadasz number, $Ra=g\alpha \beta d^2K_1\rho _m(\rho c)_i/ \mu \epsilon k_i$ is the Rayleigh–Darcy number, $\vphantom{{1}^{2^{2^{2}}}}\widetilde {Da}=\tilde {\mu } K_1/ \mu d^2$ is the Darcy–Brinkman number, $L=\nu _c d^2 \rho _m/ \mu$ is the collisional frequency parameter of PIP, $Pr=\mu (\rho c)_i/ k_i \rho _m$ is the Prandtl number, $F=\rho _n/\rho _m$ is the ratio of densities of neutral to ionized components of PIP, $\mathcal {H}=hd^2/\epsilon k_i$ is the scaled inter-phase heat-transfer coefficient, $G=g\rho _m/ (\rho c)_i\beta$ is the compressibility parameter, $\mathcal {A}=(\rho c)_s k_i/ (\rho c)_i K_s$ is the diffusivity ratio and $\gamma =\epsilon k_i/ (1-\epsilon ) k_s$ is the porosity-modified conductivity ratio.

To investigate the impact of the collisional frequency of two components of PIP on the system's stability, we neglect the convective term in (3.14) as it is independent of the collisional frequency. Subsequently, using the normal mode analysis technique on (3.14), we obtain

(3.15)\begin{equation} {{\boldsymbol{q}}_{d}}=\left(\frac{Pr L}{n+L}\right)\boldsymbol{q}, \end{equation}

where $n (\equiv \partial /\partial t)$ is the frequency of harmonic disturbances (Sharma & Sharma Reference Sharma and Sharma1978). Substituting (3.15) into (3.10), we get

(3.16)\begin{equation} \frac{1}{Va} \left[1+\frac{FLPr}{n+LPr}\right] \frac{\partial \boldsymbol{q}}{\partial t} ={-}\boldsymbol{\nabla} p + \sqrt{Ra}\theta \hat{\boldsymbol{k}} - \boldsymbol{q} + \widetilde{Da}\nabla^2 \boldsymbol{q}. \end{equation}

The boundary conditions (BCs) associated with (3.11)–(3.13) and (3.16) are

(3.17)\begin{equation} \boldsymbol{q}=\boldsymbol{0}, \quad \theta =0, \quad\psi=0 \quad \text{on} \ z=0, 1, \end{equation}

and plane tiling periodicity in $x$ and $y$ is satisfied by $\boldsymbol {{q}}$, $\theta$ and $\psi$ (Straughan Reference Straughan2006).

4 Nonlinear analysis

To perform a nonlinear analysis, we employ the energy method. Let $\| \cdot \|$ denote the ${{L}^{2}}( V )$ norm, $\langle \cdot \rangle$ denote integration over $V$ and $V$ denote the three-dimensional periodicity cell. To begin, multiply (3.16) by $\boldsymbol {q}$, (3.12) by $\theta$ and (3.13) by $\psi /\gamma$. Integrating the resulting equations over $V$ and using (3.11) and the BCs, we obtain

(4.1)\begin{gather} \frac{1}{2Va} \left[1+\frac{FLPr}{n+LPr}\right] \frac{{\rm d}}{{\rm d} t}{{\left\| \boldsymbol{q} \right\|}^{2}}= \sqrt{Ra} \langle w\theta \rangle-\left\| \boldsymbol{q} \right\|^2 -\widetilde{Da} \left\| \boldsymbol{\nabla}\boldsymbol{q} \right\|^2, \end{gather}

(4.2)\begin{gather}\frac{1}{2} \frac{{\rm d}}{{\rm d}t} \left\| \theta \right\|^2 = \sqrt{Ra} \left( 1-\frac{1}{G} \right) \langle w\theta \rangle - \mathcal{H} \langle \theta \left( \theta -\psi \right) \rangle -\left\| \boldsymbol{\nabla} \theta \right\|^2 , \end{gather}

(4.3)\begin{gather} \frac{\mathcal{A}}{2\gamma} \frac{{\rm d}}{{\rm d}t} \left\| \psi \right\|^2 ={-} \mathcal{H} \langle \psi \left(\psi-\theta\right) \rangle -\frac{1}{\gamma} \left\| \boldsymbol{\nabla} \psi \right\|^2. \end{gather}

By adding (4.1)–(4.3), we get

(4.4)\begin{equation} \frac{{\rm d}E}{{\rm d}t}=I_0-D_0, \end{equation}

where $E(t)$ is energy, $I_0$ is the production term and $D_0$ is the dissipation term. Here,

(4.5)\begin{gather} E(t)=\frac{1}{2} \left\| \theta \right\|^2 + \frac{\lambda}{2Va} \left[1+\frac{FLPr}{n+LPr}\right] {{\left\| \boldsymbol{q} \right\|}^{2}} + \frac{\mathcal{A}}{2\gamma} \left\| \psi \right\|^2, \end{gather}

(4.6)\begin{gather} I_0=\lambda\sqrt{Ra} \langle w\theta \rangle + \sqrt{Ra} \left( 1-\frac{1}{G} \right) \langle w\theta \rangle, \end{gather}

(4.7)\begin{gather}D_0=\lambda \left( \widetilde{Da} \left\| \nabla\boldsymbol{q} \right\|^2 + \left\| \boldsymbol{q} \right\|^2 \right)+ \left\| \boldsymbol{\nabla} \theta \right\|^2 + \mathcal{H} \left\| \theta-\psi \right\|^2 + \frac{1}{\gamma} \left\| \boldsymbol{\nabla} \psi \right\|^2, \end{gather}

where $\lambda$ is the coupling parameter.

From (4.4), we obtain

(4.8)\begin{equation} \frac{{\rm d}E}{{\rm d}t} \leqslant \left(m-1 \right)D_0, \end{equation}

where

(4.9)\begin{equation} m=\underset{\mathcal{G}}{{\max}}\,\frac{I_0}{D_0}, \end{equation}

and $\mathcal {G}$ is the admissible solution space.

Using the Poincaré inequality, (4.7) becomes

(4.10)\begin{equation} D_0\geqslant K^{*} E(t), \end{equation}

where

(4.11)\begin{equation} K^{*}=2{\rm \pi}^2\min\left\{1,\frac{1}{\mathcal{A}},\frac{1}{Va} \left(\frac{1}{{\rm \pi}^2}+\widetilde{Da}\right)\left[1+\frac{FLPr}{n+LPr}\right]^{{-}1}\right\}. \end{equation}

Using inequality (4.10) in inequality (4.8), we get

(4.12)\begin{equation} \frac{{\rm d}E}{{\rm d}t} \leqslant -a_0K^{*} E(t), \end{equation}

where $a_0=1-m>0$.

By integrating inequality (4.12), we get

(4.13)\begin{equation} E(t) \leqslant E(0)\exp({-a_oK^*t}). \end{equation}

From inequality (4.13), it follows that $E( t )\to 0$ exponentially as $t\to \infty$. This ensures the stability for all values of $E( 0 )$. The collisional frequency plays a significant role in energy decay, particularly when

(4.14)\begin{equation} \left\{1, \frac{1}{\mathcal{A}}\right\} > \frac{1}{Va}\left(\frac{1}{{\rm \pi}^2}+\widetilde{Da}\right) \left[1+\frac{FLPr}{n+LPr}\right]^{{-}1}. \end{equation}

Similarly, the thermal diffusivity ratio is important if

(4.15)\begin{equation} \left\{1,\frac{1}{Va}\left(\frac{1}{{\rm \pi}^2}+\widetilde{Da}\right) \left[1+\frac{FLPr}{n+LPr}\right]^{{-}1} \right\}> \frac{1}{\mathcal{A}}. \end{equation}

These conditions define thresholds where the collisional frequency and thermal diffusivity ratio significantly influence the stability and energy dynamics of the system under consideration.

4.1 Eigenvalue problem for nonlinear analysis

The value of $Ra$ is calculated from the Euler–Lagrange equations derived from (4.9), namely

(4.16)\begin{gather} \frac{\sqrt{Ra}}{\sqrt{\lambda}} \left[ \lambda+\frac{G-1}{G} \right] \theta -2\boldsymbol{q} +2\widetilde{Da}\nabla^2\boldsymbol{q}=0, \end{gather}

(4.17)\begin{gather}\frac{\sqrt{Ra}}{\sqrt{\lambda}} \left[ \lambda+\frac{G-1}{G} \right] w + 2 \nabla^2\theta - 2\mathcal{H}\left( \theta -\psi \right)=0, \end{gather}

(4.18)\begin{gather}\frac{2}{\gamma} \nabla^2\psi + 2\mathcal{H} \left( \theta -\psi \right)=0. \end{gather}

Operating $\hat {\boldsymbol {k}}\boldsymbol {\cdot }\textrm {curlcurl}$ on (4.16), we get

(4.19)\begin{equation} \frac{\sqrt{Ra}}{\sqrt{\lambda}} \left[ \lambda+\frac{G-1}{G} \right] \nabla^{2}_{1}\theta - 2\nabla^2w +2\widetilde{Da}\nabla^4w=0, \end{equation}

where

(4.20)\begin{equation} \nabla^{2}_{1}\equiv\frac{{{\partial }^{2}}}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}}{\partial {{y}^{2}}}. \end{equation}

We assume a plane tiling of the form

(4.21)\begin{equation} \left\{ w, \theta, \psi \right\}= \left\{ W(z), \varTheta(z), \varPsi(z) \right\}f(x,y), \end{equation}

where $f$ is a planform satisfying $\nabla _{1}^{2} f +a^2 f =0$, $a$ being a wavenumber.

Using the plane form (4.21), (4.19), (4.17) and (4.18) become

(4.22)\begin{gather} \frac{\sqrt{Ra}}{\sqrt{\lambda}} \left[ \lambda+\frac{G-1}{G} \right] a^2\varTheta - 2\widetilde{Da}\left(D^2-a^2\right)^2W+ 2\left(D^2-a^2 \right)W =0, \end{gather}

(4.23)\begin{gather}\frac{\sqrt{Ra}}{\sqrt{\lambda}} \left[ \lambda+\frac{G-1}{G} \right] W + 2 \left( D^2-a^2 \right)\varTheta - 2\mathcal{H}\left( \varTheta -\varPsi \right)=0, \end{gather}

(4.24)\begin{gather}\frac{2}{\gamma} \left( D^2-a^2 \right)\varPsi + 2\mathcal{H} \left( \varTheta -\varPsi \right)=0, \end{gather}

where $D={\textrm {d}}/{\textrm {d}z}$. The BCs associated with (4.22)–(4.24) are

(4.25)\begin{equation} \left.\begin{gathered}W=\varPsi=\varTheta=0 \quad \text{for} \ z=0, 1,\\ \text{and} \quad DW=0 \quad \text{for rigid surface},\\ \text{and} \quad D^2W=0 \quad \text{for free surface}. \end{gathered}\right\} \end{equation}

The set of equations (4.22)–(4.24), together with the BCs (4.25), constitutes the eigenvalue problem for nonlinear analysis.

5 Linear analysis

To study the linear analysis, by considering the perturbation to be infinitesimally small, we ignore the nonlinear terms from the non-dimensional perturbed equations and obtain the linearized non-dimensional perturbed equations as

(5.1)\begin{gather} \frac{1}{Va} \left[1+\frac{FLPr}{n+LPr}\right] \frac{\partial \boldsymbol{q}}{\partial t} ={-}\boldsymbol{\nabla} p + \sqrt{Ra}\theta \hat{\boldsymbol{k}} - \boldsymbol{q} + \widetilde{Da}\nabla^2 \boldsymbol{q}, \end{gather}

(5.2)\begin{gather}\frac{\partial \theta}{\partial t}= \sqrt{Ra}\left(1-\frac{1}{G}\right)w + \mathcal{H}\left( \psi-\theta \right)+ \nabla^2\theta, \end{gather}

(5.3)\begin{gather}\mathcal{A} \frac{\partial \psi}{\partial t}={-}\mathcal{H}\gamma \left( \psi-\theta \right) + \nabla^2\psi. \end{gather}

Operating $\hat {\boldsymbol {k}}\boldsymbol {\cdot }\textrm {curlcurl}$ on (5.1), we get

(5.4)\begin{equation} \frac{1}{Va} \left[1+\frac{FLPr}{n+LPr}\right] \frac{\partial}{\partial t} \left( \nabla^2w \right) = \sqrt{Ra}\nabla^{2}_{1}\theta - \nabla^2w + \widetilde{Da}\nabla^4 w. \end{equation}

By applying the normal mode analysis method to equations (5.4), (5.2) and (5.3), where we assumed the perturbed quantities of the form

(5.5)\begin{equation} \left\{ w, \theta, \psi \right\}= \left\{ W(z), \varTheta(z), \varPsi(z) \right\} \exp\left\{ \iota \left( a_xx+a_yy \right) + nt \right\}. \end{equation}

Here, $\sqrt {a^2_x+a^2_y}=a$ is the wavenumber.

Using (5.5), equations (5.4), (5.2) and (5.3) become

(5.6)\begin{gather} \frac{1}{Va} \left[1+\frac{FLPr}{n+LPr}\right] n \left( D^2-a^2 \right)W ={-}\sqrt{Ra}a^2\varTheta - \left( D^2-a^2 \right)W + \widetilde{Da}\left( D^2-a^2 \right)^2W, \end{gather}

(5.7)\begin{gather}n\varTheta= \left( D^2-a^2 \right)\varTheta + \mathcal{H}\left( \varPsi-\varTheta \right) + \sqrt{Ra}\left(1-\frac{1}{G}\right)W, \end{gather}

(5.8)\begin{gather}\mathcal{A} n \varPsi= \left( D^2-a^2 \right)\varPsi-\mathcal{H}\gamma \left( \varPsi-\varTheta \right). \end{gather}

The BCs associated with (5.6)–(5.8) are same as (4.25).

5.1 Exchange of stabilities

Here, we prove that the principle of exchange of stabilities is valid, namely that marginally stable modes with $\sigma =0$ also have $\omega =0$, where $\sigma +\iota \omega =n$.

Multiply (5.6) by $W^*$ (complex conjugate of $W$), (5.7) by $\varTheta ^*$ (complex conjugate of $\varTheta$) and (5.8) by $\varPsi ^*$ (complex conjugate of $\varPsi$), integrating resulting equations over the range of $z$ and using BCs (4.25), we get

(5.9)\begin{gather} -\frac{n}{Va}\left[1+\frac{FLPr}{n+LPr}\right] \mathcal{I}_1 + \sqrt{Ra}a^2 \int_{0}^{1} W^* \varTheta \,{\rm d}z - \mathcal{I}_1-\widetilde{Da} \mathcal{I}_2=0, \end{gather}

(5.10)\begin{gather} n\mathcal{I}_3 + \mathcal{I}_4 - \mathcal{H} \int_{0}^{1} \varTheta^* \varPsi \,{\rm d}z + \mathcal{H} \mathcal{I}_3 - \sqrt{Ra}\left(1-\frac{1}{G}\right) \int_{0}^{1} W \varTheta^* \,{\rm d}z=0, \end{gather}

(5.11)\begin{gather}\mathcal{A}n\mathcal{I}_5 + \mathcal{I}_6 + \mathcal{H}\gamma\mathcal{I}_5-\mathcal{H}\gamma \int_{0}^{1} \varTheta \varPsi^* \,{\rm d}z, \end{gather}

where

(5.12)\begin{align} \left.\begin{gathered} \mathcal{I}_1=\int_{0}^{1} \left( \lvert DW \rvert^2+a^2\lvert W \rvert^2 \right) {\rm d}z, \quad \mathcal{I}_2=\int_{0}^{1} \left( \lvert D^2W \rvert^2+2a^2\lvert DW \rvert^2+a^4\lvert W \rvert^2 \right) {\rm d}z, \\ \mathcal{I}_3= \int_{0}^{1} \lvert \varTheta \rvert^2 \, {\rm d}z, \quad \mathcal{I}_4= \int_{0}^{1} \left( \lvert D\varTheta \rvert^2+a^2\lvert \varTheta \rvert^2 \right) {\rm d}z, \\ \mathcal{I}_5=\int_{0}^{1} \lvert \varPsi \rvert^2 \, {\rm d}z, \quad \mathcal{I}_6=\int_{0}^{1} \left( \lvert D\varPsi \rvert^2 +a^2 \lvert \varPsi \rvert^2 \right) {\rm d}z. \end{gathered}\right\} \end{align}

The integrals $\mathcal {I}_1$, $\mathcal {I}_2$, $\mathcal {I}_3$, $\mathcal {I}_4$, $\mathcal {I}_5$ and $\mathcal {I}_6$ are positive. So, subtracting the product of $a^2$ with (5.11) from the sum of the product of $(1-1/G)\gamma$ with (5.9) and the product of $a^2\gamma$ with complex conjugate of (5.10), we get

(5.13)\begin{align} & -\gamma \left(1-\frac{1}{G}\right) \left( \frac{n}{Va}\left[1+\frac{FLPr}{n+LPr}\right] +1 \right) \mathcal{I}_1 - \left(1-\frac{1}{G}\right)\gamma\widetilde{Da}\mathcal{I}_2 \nonumber\\ & \quad + a^2\gamma\left( n^*+\mathcal{H} \right)\mathcal{I}_3 + a^2\gamma\mathcal{I}_4 - a^2 \left( \mathcal{A}n+\gamma\mathcal{H} \right)\mathcal{I}_5 - a^2\mathcal{I}_6 = 0 , \end{align}

where $n^*$ is the complex conjugate of $n$.

The imaginary part of (5.13) gives

(5.14)\begin{equation} \omega \left\{ \left(1-\frac{1}{G}\right) \left( \frac{1}{Va} + \frac{FL^2(Pr)^2}{\left(\sigma+LPr\right)^2+\omega^2} \right)\gamma\mathcal{I}_1 + a^2\gamma\mathcal{I}_3 + a^2\mathcal{A}\mathcal{I}_5 \right\}=0. \end{equation}

The quantity inside the curly brackets is positive definite for $G$ is greater than $1$. Hence, (5.14) suggests that $\omega =0$. This establishes that the principle of exchange of stabilities is valid (Chandrasekhar Reference Chandrasekhar1981).

Now, keeping in mind the validation of the principle of exchange of stabilities, (5.6)–(5.8) can be rewritten as

(5.15)\begin{gather} \sqrt{Ra}a^2\varTheta + \left( D^2-a^2 \right)W - \widetilde{Da}\left( D^2-a^2 \right)^2W=0, \end{gather}

(5.16)\begin{gather}\left( D^2-a^2 \right)\varTheta + \mathcal{H}\left( \varPsi-\varTheta \right) + \sqrt{Ra}\left(1-\frac{1}{G}\right)W=0, \end{gather}

(5.17)\begin{gather}\left( D^2-a^2 \right)\varPsi-\mathcal{H}\gamma \left( \varPsi-\varTheta \right)=0. \end{gather}

The set of equations (5.15)–(5.17), together with the BCs (4.25), constitutes the eigenvalue problem for linear analysis.

6 Method of solution

The eigenvalue problems for both nonlinear as well as linear analyses have been solved using the single-term Galerkin method. In this method, the weighted functions are the same as the base (trial) functions (Yadav et al. Reference Yadav, Bhargava and Agarwal2013). Accordingly, we define $W$, $\varTheta$ and $\varPsi$ in the following form:

(6.1a-c)\begin{equation} W=\sum_{j=1}^{N}{{{A}_{j}}{{W}_{j}}}, \quad \varTheta=\sum_{j=1}^{N}{{{B}_{j}}{{\varTheta}_{j}}}, \quad \varPsi=\sum_{j=1}^{N}{{{C}_{j}}{{\varPsi}_{j}}}, \end{equation}

where ${{A}_{j}}$, ${{B}_{j}}$ and $C_j$ are unknown coefficients, $j= 1,2,3,\ldots, N$ and the base functions ${{W}_{j}}$, ${{\varTheta }_{j}}$ and $\varPsi _j$ are assumed to be in the following forms for free–free, rigid–free and rigid–rigid bounding surfaces (Chandrasekhar Reference Chandrasekhar1981; Yadav et al. Reference Yadav, Bhargava and Agarwal2013), respectively:

(6.2)\begin{gather} {{W}_{j}}=\sin (j{\rm \pi} z)={{\varTheta}_{j}}={{\varPsi }_{j}}, \end{gather}

(6.3a,b)\begin{gather}{{W}_{j}}=z^2\left(1-z\right)\left[ \left(j+2\right)-2z^j \right], \quad {{\varTheta}_{j}}=z^j-z^{j+1}={{\varPsi }_{j}}, \end{gather}

(6.4a,b)\begin{gather} {{W}_{j}}=z^{j+1}-2z^{j+2}+z^{j+3}, \quad {{\varTheta}_{j}}=z^j-z^{j+1}={{\varPsi }_{j}}, \end{gather}

such that $W_j$, $\varTheta _j$ and $\varPsi _j$ satisfy the corresponding BCs (4.25). By substituting the expressions for $W$, $\varTheta$ and $\varPsi$ into (4.22)–(4.24) and (5.15)–(5.17), and then multiplying the first equation by $W_j$, the second equation by $\varTheta _j$ and the third equation by $\varPsi _j$ and integrating the resulting equations over the interval from zero to unity, we obtain a set of linear homogeneous equations. This set of equations admits a non-trivial solution only if its determinant is equal to zero, which gives the characteristic equations of the eigenvalue problems in terms of the Rayleigh–Darcy number $Ra$.

6.1 Rayleigh–Darcy number for free–free bounding surfaces

For nonlinear analysis, $Ra$ has been found as a function of $a$, $\lambda$, $\mathcal {H}$, $\gamma$ and $G$. The expression for $Ra$ is given by

(6.5)\begin{equation} Ra=\frac{4G^2\left(a^2+{\rm \pi}^2\right)^2 \left(1+a^2\widetilde{Da}+{\rm \pi}^2\widetilde{Da}\right) \left( a^2+{\rm \pi}^2+\mathcal{H}+\gamma\mathcal{H} \right)\lambda}{a^2 \left( a^2+{\rm \pi}^2+\gamma\mathcal{H} \right)\left(G-1+G\lambda\right)^2}. \end{equation}

The optimal value of $\lambda$ has been found using the condition $\textrm {d}Ra/\textrm {d}\lambda =0$, which yields $\lambda =(G-1)/G$. Using this value, the expression for $Ra$ becomes

(6.6)\begin{equation} Ra=\left(\frac{G}{G-1}\right)\frac{\left(a^2+{\rm \pi}^2\right)^2 \left(1+a^2\widetilde{Da}+{\rm \pi}^2\widetilde{Da}\right) \left( a^2+{\rm \pi}^2+\mathcal{H}+\gamma\mathcal{H} \right)}{a^2 \left( a^2+{\rm \pi}^2+\gamma\mathcal{H} \right)}. \end{equation}

For linear analysis, the value of $Ra$ found by solving the eigenvalue problem (5.15)–(5.17) using the Galerkin method has been the same as (6.6).

6.2 Rayleigh–Darcy number for rigid–free bounding surfaces

For nonlinear analysis, $Ra$ has been found as a function of $a$, $\lambda$, $\mathcal {H}$, $\gamma$ and $G$. The expression for $Ra$ is given by

(6.7)\begin{align} Ra=\frac{112G^2\left(10+a^2\right) \left( 19a^4\widetilde{Da}+216\left(1+21\widetilde{Da}\right)+a^2\left(19+432\widetilde{Da}\right)\right) \left( 10+ a^2+\mathcal{H}+\gamma\mathcal{H} \right)\lambda}{507a^2 \left( 10+a^2+\gamma\mathcal{H} \right)\left(G-1+G\lambda\right)^2}. \end{align}

The optimal value of $\lambda$ has been found using the condition $dRa/d\lambda =0$, which yields $\lambda =(G-1)/G$. Using this value, the expression for $Ra$ becomes

(6.8)\begin{align} Ra=\frac{28G\left(10+a^2\right) \left( 19a^4\widetilde{Da}+216\left(1+21\widetilde{Da}\right)+a^2\left(19+432\widetilde{Da}\right)\right) \left( 10+ a^2+\mathcal{H}+\gamma\mathcal{H} \right)}{507a^2 \left(G-1\right)\left( 10+a^2+\gamma\mathcal{H} \right)}. \end{align}

For linear analysis, the value of $Ra$ found by solving the eigenvalue problem (5.15)–(5.17) using the Galerkin method has been the same as (6.8).

6.3 Rayleigh–Darcy number for rigid–rigid bounding surfaces

For nonlinear analysis, $Ra$ has been found as a function of $a$, $\lambda$, $\mathcal {H}$, $\gamma$ and $G$. The expression for $Ra$ is given by

(6.9)\begin{align} Ra=\frac{112G^2\left(10+a^2\right) \left( 12+504\widetilde{Da}+a^4\widetilde{Da}+a^2\left(1+24\widetilde{Da}\right)\right) \left( 10+ a^2+\mathcal{H}+\gamma\mathcal{H} \right)\lambda}{27a^2 \left( 10+a^2+\gamma\mathcal{H} \right)\left(G-1+G\lambda\right)^2}. \end{align}

The optimal value of $\lambda$ has been found using the condition $\textrm {d}Ra/\textrm {d}\lambda =0$, which yields $\lambda =(G-1)/G$. Using this value, the expression for $Ra$ becomes

(6.10)\begin{align} Ra=\frac{28G\left(10+a^2\right) \left( 12+504\widetilde{Da}+a^4\widetilde{Da}+a^2\left(1+24\widetilde{Da}\right)\right) \left( 10+ a^2+\mathcal{H}+\gamma\mathcal{H} \right)}{27a^2 \left(G-1\right)\left( 10+a^2+\gamma\mathcal{H} \right)}. \end{align}

For linear analysis, the value of $Ra$ found by solving the eigenvalue problem (5.15)–(5.17) using the Galerkin method has been the same as (6.10).

Here, $Ra$ has the same value for both nonlinear and linear analyses, indicating a lack of subcritical regions and demonstrating strong global stability. It has been observed that collisional effects contribute to energy decay, but do not affect the value of $Ra$ in either nonlinear or linear analyses.

8 Effect of magnetic field

To study the effect of the magnetic field on the problem of instability described in § 2, we consider a horizontal porous layer of PIP, subjected to a magnetic field $\boldsymbol {H}$. We have adopted the quasi-static magnetohydrodynamic approximation proposed by Galdi & Straughan (Reference Galdi and Straughan1985). The key modification is to include a term representing the Lorentz force ($\boldsymbol {\mathcal {L}}=\boldsymbol {j}\times \boldsymbol {B}_0$) in the equation of motion (2.3), where $\boldsymbol {j}$ is the current and $\boldsymbol {B}_0=(0,0, B_0)$ is the magnetic field with only vertical component (Chandrasekhar Reference Chandrasekhar1981). The non-dimensional perturbed equations (3.11)–(3.13) remains the same but (3.16) becomes

(8.1)\begin{gather} \frac{1}{Va} \left[1+\frac{FLPr}{n+LPr}\right] \frac{\partial \boldsymbol{q}}{\partial t} ={-}\nabla p + \sqrt{Ra}\theta \hat{\boldsymbol{k}} - \boldsymbol{q} + \widetilde{Da}\nabla^2 \boldsymbol{q}+ \frac{M^2\widetilde{Da}}{Vr} \left[\left(\boldsymbol{q}\times\hat{\boldsymbol{k}}\right)\times\hat{\boldsymbol{k}}\right], \end{gather}

where $M^2(=\sigma _1\boldsymbol {B}^{2}_{0}d^2/\epsilon \rho _m\nu )$ is the Hartmann number, which accounts for the effect of the magnetic field. Here, $\sigma _1$ is the electrical conductivity.

The eigenvalue problem for nonlinear analysis incorporating the effect of magnetic field is

(8.2)\begin{gather} \frac{\sqrt{Ra}}{\sqrt{\lambda}} \left[ \lambda+\frac{G-1}{G} \right] a^2\varTheta + 2\left(D^2-a^2 \right)W - 2\widetilde{Da}\left(D^2-a^2\right)^2W + \frac{2M^2\widetilde{Da}}{Vr}D^2W=0,\end{gather}

(8.3)\begin{gather}\frac{\sqrt{Ra}}{\sqrt{\lambda}} \left[ \lambda+\frac{G-1}{G} \right] W + 2 \left( D^2-a^2 \right)\varTheta - 2\mathcal{H}\left( \varTheta -\varPsi \right)=0, \end{gather}

(8.4)\begin{gather}\frac{2}{\gamma} \left( D^2-a^2 \right)\varPsi + 2\mathcal{H} \left( \varTheta -\varPsi \right)=0. \end{gather}

The eigenvalue problem for linear analysis incorporating the effect of magnetic field is

(8.5) \begin{gather} \sqrt{Ra}a^2\varTheta + \left( D^2-a^2 \right)W - \widetilde{Da}\left( D^2-a^2 \right)^2W+\frac{M^2\widetilde{Da}}{Vr}D^2W=0, \end{gather}

(8.6)\begin{gather}\left( D^2-a^2 \right)\varTheta + \mathcal{H}\left( \varPsi-\varTheta \right) + \sqrt{Ra}\left(1-\frac{1}{G}\right)W=0, \end{gather}

(8.7)\begin{gather}\left( D^2-a^2 \right)\varPsi-\mathcal{H}\gamma \left( \varPsi-\varTheta \right)=0. \end{gather}

The BCs associated with both eigenvalue problems remain the same as (4.25).

The Rayleigh–Darcy numbers for both nonlinear and linear analyses were determined by solving the eigenvalue problems using the Galerkin method as detailed in § 6. The value of $Ra$ for nonlinear analysis has been a function of $a$, $\lambda$, $\mathcal {H}$, $\gamma$, $M^2$ and $G$. The optimal value of $\lambda$ is $(G-1)/G$ for all combinations of bounding surfaces. Using this value of $\lambda$, we have found that $Ra$ for nonlinear analysis has been the same as $Ra$ for linear analysis. The Rayleigh–Darcy numbers ($Ra$) for distinct combinations of bounding surfaces are as follows.

For free–free bounding surfaces

(8.8)\begin{equation} Ra=\left(\frac{G}{G-1}\right)\frac{\left(a^2+{\rm \pi}^2\right) \left(a^2+{\rm \pi}^2+\left(a^2+{\rm \pi}^2\right)^2\widetilde{Da}+{\rm \pi}^2\widetilde{Da}M^2\right) \left( a^2+{\rm \pi}^2+\mathcal{H}+\gamma\mathcal{H} \right)}{a^2 \left( a^2+{\rm \pi}^2+\gamma\mathcal{H} \right)}. \end{equation}

For rigid–free bounding surfaces

(8.9)\begin{align} Ra& = \frac{28G\left(10+a^2\right) \left( 10+ a^2+\mathcal{H}+\gamma\mathcal{H} \right)}{507a^2 \left(G-1\right)\left( 10+a^2+\gamma\mathcal{H} \right)} \nonumber\\ & \quad \times\left[ 19a^4\widetilde{Da}+216\left(1+21\widetilde{Da}\right)+a^2\left(19+432\widetilde{Da}\right) +216\widetilde{Da}M^2\right]. \end{align}

For rigid–rigid bounding surfaces

(8.10)\begin{align} Ra& = \frac{28G\left(10+a^2\right) \left( 10+ a^2+\mathcal{H}+\gamma\mathcal{H} \right)}{27a^2 \left(G-1\right)\left( 10+a^2+\gamma\mathcal{H} \right)} \nonumber\\ & \quad \times \left[ 12+504\widetilde{Da}+a^4\widetilde{Da}+a^2\left(1+24\widetilde{Da}\right) +12\widetilde{Da}M^2\right]. \end{align}

For the limiting case of the magnetic field i.e. $M^2=0$, clearly the Rayleigh–Darcy number for the respective bounding surfaces has been the same as given in (6.6), (6.8) and (6.10).

Figure 6 shows the variation in $Ra_c$ with the square of the Hartmann number ($M^2$) for different boundary surface combinations. The figure shows that, with the increase in $M^2$, $Ra_c$ also increases, indicating that magnetic fields delay the onset of thermal convection, thus exerting a stabilizing influence. This occurs because the magnetic field impedes the motion within the PIP. When the PIP attempts convective motion, the Lorentz force generated by the magnetic field counteracts this movement, resisting the flow and suppressing convective instabilities, thereby delaying convection onset. This suppression of plasma motion is essential in postponing or preventing the onset of thermal convection. Additionally, the figure reveals that PIP confined within rigid–rigid bounding surfaces exhibits greater thermal stability compared with PIP confined within rigid–free or free–free bounding surfaces.

Figure 6. Variation of critical Rayleigh–Darcy number ($Ra_c$) with the square of Hartmann number ($M^2$) for the distinct combinations of bounding surfaces.

The stability analysis of PIP in a porous medium with LTNE effects has significant applications across various scientific and engineering domains. In astrophysics, this study enhances the understanding of stellar and planetary atmospheres, contributing to more accurate models of these complex systems. In the field of nuclear fusion, it offers insights into the stability of magnetic confinement systems, which is crucial for the development of efficient fusion reactors. The findings also aid in the design of thermal protection systems for spacecraft, ensuring their integrity during re-entry and other high-temperature conditions. Furthermore, the optimization of heat exchangers and combustion chambers in engineering applications benefits from the improved thermal stability insights provided by this study. Additionally, the study informs the development of plasma-based pollution control technologies and advances in plasma medicine, particularly in the design of stable plasma sources for medical treatments involving porous tissues or materials.