Linear disturbance growth induced by viscous dissipation in Darcy–Bénard convection with throughflow

Abstract

Modal and non-modal linear stability analyses are employed to investigate the effect of internal and external heating on disturbance temporal growth for the Darcy–Bénard convection with throughflow. A matrix-forming approach is employed for both purposes, where the generalised eigenvalue problem is built using the generalised integral transform technique. Although the disturbance equations are not self-adjoint, the non-modal analysis indicates that there is no transient growth. Hence, any disturbance growth in time must be induced by modal mechanisms. An absolute instability analysis reveals that viscous dissipation has a destabilising effect and introduces new modes that are eventually destabilised by increasing the Péclet number. Beyond critical values of the Péclet number, where codimension-two absolutely unstable points exist, these modes become more unstable than the classical mode found in the absence of viscous dissipation, which is stabilised by an increasing Péclet number. This internal heating mechanism generated by viscous dissipation is so strong at high enough Péclet numbers that instability becomes possible through heating from above.

1. Introduction

Transition from stability to instability in fluid flows is a subject widely explored in the literature. Bénard (1901) was among the first to study it experimentally when he observed the appearance of convection cells in a thin layer of fluid after heating it from below beyond a certain critical temperature difference. The first explanation for this phenomenon was proposed by Rayleigh (1916), who used a linear stability analysis to introduce buoyancy as the driving mechanism, which is the reason why this became known as the Rayleigh–Bénard problem. Many decades later, however, Pearson (1958) also used linear stability analysis to propose a second explanation. He pointed out that the small film thickness employed by Lord Rayleigh essentially renders buoyancy effects negligible and, in turn, promotes surface tension gradients as the driving mechanism. The natural convection described by the Rayleigh–Bénard problem also occurs in a fluid saturated porous medium, in which case it is known as the Darcy–Bénard problem. Horton & Rogers (1945) as well as Lapwood (1948) were the first to investigate its linear instability. Prats (1966) was the first to extend it to include a horizontal throughflow and, hence, to consider mixed convection.

The above-mentioned studies focus on the asymptotic behaviour of individual normal modes in either time or space. This eigenvalue, or modal, based description is the traditional way linear stability analyses have been performed (Chandrasekhar 1961; Drazin & Reid 1981). Some key developments, however, have appeared since then. When the disturbance is local, i.e. it does not vary in the direction of the longitudinal flow, at least in an approximate sense, convective and absolute instability concepts took over the previous temporal/spatial understanding (Huerre & Monkewitz 1990). When the disturbance is no longer local, parabolised stability (Herbert 1997) and global stability (Theofilis 2011) concepts can be employed. Adjoint equations (Luchini & Bottaro 2014) also deserve special mention, since they have been connected to absolute instability in both discrete (Lesshafft & Marquet 2010) and continuous (Alves et al. 2019) senses. Many of these techniques have been applied to both natural and mixed convection in porous media. Most cases focus on convective instabilities (Nield & Bejan 2006), but there have been some recent attempts to investigate absolute instabilities as well (Barletta 2019; Barletta, Celli & Alves 2020; Schuabb, Alves & Hirata 2020).

On the other hand, the short time/space behaviour of superposed normal modes can also be relevant. This non-modal description is key in explaining transition in modally stable flows under subcritical conditions (Schmid 2007). In general, Rayleigh–Bénard-type problems are self-adjoint and non-modal growth is not possible. This can change, however, in the presence of throughflow. Biau & Bottaro (2004) studied the effect of stable thermal stratification in shear flows whereas Sameen & Govindarajan (2007) studied the effect of wall heating in channel flows from the perspective of both modal and non-modal linear growth. Finally, Jerome, Chomaz & Huerre (2012) studied non-modal growth in both Rayleigh–Bénard–Poiseuille and Rayleigh–Bénard–Couette problems. In the context of porous media flows, only a few studies have dealt with the non-modal linear growth of disturbances, but they focused on density-driven instability (Rapaka et al. 2008, 2009).

One of the main goals of the present paper is to fill this gap, investigating modal as well as non-modal mechanisms for linear disturbance temporal growth that might exist for flows in porous media. Another goal is to do so while also considering the influence of viscous dissipation effects, often neglected due to their small magnitude. Gebhart (1962) has identified, however, the parametric conditions under which viscous dissipation can be relevant for natural convection in pure fluids. Furthermore, mixed convection renders viscous dissipation effects even more important due to the added forced component. In the context of porous media, the former was shown to be true by Nakayama & Pop (1989), where Murthy (1998) extended this study to show that this is also true for the latter. According to Gebhart (1962), viscous dissipation effects can be dominant in many scenarios. He mentioned processes under a strong gravitational field (e.g. on larger planets), devices operating at high rotative speeds (e.g. internal cooling of turbine blades) as well as processes with large characteristic lengths (e.g. geophysical flows). Nield (2000) included particle bed nuclear reactors among possible interesting applications. Furthermore, Magyari, Rees & Keller (2005) pointed out that many natural convection processes could be qualitatively altered by viscous dissipation effects even when they appear negligible. The characteristic dimensionless parameter quantifying the strength of the viscous dissipation effect for buoyant flows is today known as Gebhart number, $Ge$, which can be interpreted as the ratio between the kinetic energy of the flow and the heat transferred to the fluid (Gebhart 1962). Gebhart (1962) and, after him, Turcotte et al. (1974) were the first to propose it, although they did so using a different name, i.e. the dissipation number. Nevertheless, they showed that such a parameter is usually small, but can achieve order unity under the aforementioned scenarios.

The effect of viscous dissipation on the onset of instability for several mixed convection problems in porous media was investigated in a series of studies. Barletta, Celli & Rees (2009a) and Storesletten & Barletta (2009) were the first to study the particular case where the internal heating generated by viscous dissipation is the sole cause of destabilisation. In their studies, no additional thermal forcing was imposed either internally or externally, namely through the walls or otherwise. The combined effect of internal heating through viscous dissipation and external heating through different thermal boundary conditions was studied by Barletta & Storesletten (2010) as well as Barletta, Celli & Nield (2010) and Nield, Barletta & Celli (2011). Nield & Barletta (2010) also explored two different models for the viscous dissipation effect. Viscous dissipation effects on non-Darcy models for flows in porous media were explored by Barletta, Celli & Rees (2009b) and Barletta, di Schio & Celli (2011b) whereas thermal non-equilibrium, heterogeneity and viscoelastic fluid models were explored by Barletta & Celli (2011), Barletta, Celli & Kuznetsov (2011a) and Alves et al. (2014), respectively. Roy & Murthy (2015) investigated the Soret effect on the double diffusive convection, where the convection is occurring just by means of viscous dissipation effects. Later on, Roy & Murthy (2017) studied the influence of viscous heating on the instability induced by an inclined temperature gradient. More recently, Barletta & Mulone (2021) showed that the classical problem studied by Horton & Rogers (1945) and Lapwood (1948) is conditionally stable from a nonlinear point of view in the presence of viscous dissipation.

These studies show that the presence of viscous dissipation has a significant impact on the onset of instability for flows in porous media when compared with the classical scenario where this effect is absent (Prats 1966). For instance, throughflow destabilises the onset of instability in the presence of viscous dissipation but it has no effect on this onset in the absence of viscous dissipation. This destabilising role of the viscous dissipation effects on the transition from stable conditions to convective instability also appears in the case of buoyant flows in clear fluids, as was demonstrated by Requilé, Hirata & Ouarzazi (2020). The transition from convective to absolute instability, on the other hand, is stabilised by throughflow in the absence of viscous dissipation (Hirata & Ouarzazi 2010). This is due to the fact that disturbances require more thermal energy from the base flow to be able to propagate upstream as the throughflow becomes stronger. The same is not true in the presence of viscous dissipation, although this is not yet entirely clear in this case (Brandão, Alves & Barletta 2014).

This literature review shows that linear disturbance temporal growth for convective porous media flows is not yet fully understood from an asymptotic (modal) perspective in the presence of viscous dissipation. Furthermore, this is not understood at all from a transient (non-modal) perspective, with or without viscous dissipation. Hence, the present paper explores both issues in detail in an attempt to fill these gaps. This is done here by investigating the possibility of transient disturbance growth, first by looking at eigenvector orthogonality, and then by performing a non-modal stability analysis based on optimal initial conditions. Then, an absolute instability analysis is performed in order to understand the time asymptotic disturbance behaviour. Section 2 shows the mathematical formulation of the physical problem, as well as the derivation of the linear disturbance equations. Section 3 shows the methodologies considered here to solve the eigenvalue problem, and also the particular issues of the non-modal and modal analyses. Section 4 discusses the results, while § 5 addresses the most relevant conclusions. Finally, the reader can find a convergence analysis in Appendix A and more details about the eigenvector matrix analysis in Appendix B.

2. Mathematical formulation

Fluid flow through a horizontal porous channel with a vertical temperature gradient induced by external heating from below and an internal heating induced by viscous dissipation is considered. The channel walls, located at $z=0,1$, are assumed impermeable with prescribed temperatures, where the lower boundary is hotter than the upper boundary. Momentum transfer is modelled by Darcy's law, where thermal equilibrium is assumed between solid and fluid phases for the local energy balance equation. Furthermore, viscous dissipation is taken into account and the Oberbeck–Boussinesq approximation is assumed valid. Therefore, the governing equations of the present problem can be written as

(2.1)\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0 , \end{gather}

(2.2)\begin{gather} \boldsymbol{u} = Ra T \hat{\boldsymbol{k}} - \boldsymbol{\nabla} P \end{gather}

and

(2.3)\begin{equation} \frac{\partial T}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} T = \nabla^2 T + \frac{Ge}{Ra} \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u} , \end{equation}

which is subject to the boundary conditions

(2.4a,b)\begin{equation} w = 0 \quad \text{and}\quad T = 1 \quad \text{at} \ z = 0 \end{equation}

and

(2.5a,b)\begin{equation} w = 0 \quad \text{and}\quad T = 0 \quad \text{at} \ z = 1 , \end{equation}

where the following dimensionless quantities are employed

(2.6a–c)\begin{gather} \boldsymbol{u} = \frac{\boldsymbol{u}^{*}}{\chi^{*}/h^{*}} ,\quad \boldsymbol{x}=\frac{\boldsymbol{x}}{h^{*}} ,\quad t=\frac{t^{*}}{\sigma {h^{*}}^2/\chi^{*}} , \end{gather}

(2.7a,b)\begin{gather} P=\frac{P^{*}}{\mu^{*}\chi^{*}/\kappa^{*}} \quad \mathrm{and} \quad T=\frac{T^{*}-T_0^{*}}{T_h^{*}-T_0^{*}} , \end{gather}

leading to the following definitions for the Rayleigh, Gebhart and Péclet numbers,

(2.8a–c)\begin{equation} Ra=\frac{\rho^{*} g^{*}\beta^{*}(T_0^{*}-T_h^{*})\kappa^{*}h^{*}}{\mu^{*}\chi^{*}} ,\quad Ge=\frac{g^{*}\beta^{*}h^{*}}{c^*} \quad \mathrm{and} \quad Pe=\frac{u_0^{*}h^{*}}{\chi^{*}} , \end{equation}

where the superscript asterisk denotes a dimensional quantity. In particular, $h^*$ is the distance between channel walls, $\boldsymbol {u}^{*}=\{u^{*},v^{*},w^{*}\}$ is the velocity vector, $\boldsymbol {x}^{*}=\{x^{*},y^{*},z^{*}\}$ is the coordinate vector, $\mu ^{*}$ is the dynamic viscosity, $\kappa ^{*}$ is the permeability, $t^{*}$ is the time coordinate, $T^{*}$ is the temperature, $\rho ^{*}$ is the fluid density at the reference temperature $T_0^{*}$, $c^{*}$ is the specific heat of the fluid, $\chi ^{*}$ is the effective thermal diffusivity, $P^{*}$ is the gauge pressure with respect to the hydrostatic pressure, $g^{*}$ is the gravity acceleration and $\beta ^{*}$ is the fluid thermal expansion coefficient. Furthermore, $\sigma$ is the dimensionless ratio between the volumetric heat capacity of the saturated porous medium and $\rho _h^{*}c^{*}$. Finally, the lower wall temperature is $T_h^{*}$, the upper wall temperature is $T_0^{*}$ and the uniform streamwise velocity is $u_0^{*}$. See Alves et al. (2014) for more details about the model.

2.1. Asymptotic expansion

The first step employed here to study the aforementioned problem is to assume an asymptotic expansion can be used to decompose its dependent variables as

(2.9)$$\begin{gather} \boldsymbol{u}(x,y,z,t) = \boldsymbol{u}_b(z) + \epsilon \boldsymbol{u}_d(x,y,z,t) + O(\epsilon^2) , \end{gather}$$

(2.10)$$\begin{gather}T(x,y,z,t) = T_b(z) + \epsilon T_d(x,y,z,t) + O(\epsilon^2) \quad \text{and} \end{gather}$$

(2.11)$$\begin{gather}P(x,y,z,t) = P_b(z) + \epsilon\,P_d(x,y,z,t) + O(\epsilon^2) , \end{gather}$$

where the subscripts $b$ and $d$ denote base flow and disturbances, respectively, while $\epsilon$ represents a dimensionless disturbance amplitude parameter. Two key assumptions are implicit to this expansion. First, $z$ is the only inhomogeneous coordinate. This implies that the base flow can be written as a steady state that depends on $z$ alone. Second, disturbance amplitudes are small, i.e. $\epsilon \ll 1$. This implies that all nonlinear terms, represented by the $O(\epsilon ^2)$ terms, are negligible.

2.2. Base flow

Equations (2.1)–(2.4a,b) have such a steady state and it is given by

(2.12)$$\begin{gather} \boldsymbol{u}_b(Z) = \{ Pe, 0, 0 \} , \end{gather}$$

(2.13)$$\begin{gather}T_b(z) = 1 - z + \frac{Ge Pe^2}{2 Ra}(1 - z)z \quad \text{and} \end{gather}$$

(2.14)$$\begin{gather}P_b(x,z) = P_0 - Pe\,x + \frac{Ra}{2}(2-z)z + \frac{Ge Pe^2}{12}(3-2z)z^2 , \end{gather}$$

where $P_0$ is a reference pressure. Two things are worth noting about the effect of viscous dissipation on this steady state. First, it can only act in the presence of throughflow ($Pe \ne 0$), because the steady state is a rest state otherwise. Second, it can create a stable temperature stratification near the hot wall, but only when the throughflow is strong enough ($Pe \gg 1$). Finally, the steady pressure field is non-local, i.e. it has a linear dependence on the $x$ coordinate. However, only the pressure gradient appears in the governing equations. Hence, this steady pressure gradient depends on $z$ alone, where its component in the $x$ direction is constant and responsible for the generation of a steady streamwise throughflow.

2.3. Linear disturbances

Substituting (2.9)–(2.11) into (2.1)–(2.4a,b), cancelling out the ${O}(\epsilon ^0)$ steady terms and neglecting the ${O}(\epsilon ^2)$ nonlinear terms, leaves the ${O}(\epsilon )$ linear terms that form the linear and homogeneous disturbance equations

(2.15)\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}_d = 0 , \end{gather}

(2.16)\begin{gather} \boldsymbol{u}_d = Ra T_d \hat{\boldsymbol{k}} - \boldsymbol{\nabla} P_d \end{gather}

and

(2.17)\begin{equation} \frac{\partial T_d}{\partial t} + \boldsymbol{u}_b \boldsymbol{\cdot} \boldsymbol{\nabla} T_d + \boldsymbol{u}_d \boldsymbol{\cdot} \boldsymbol{\nabla} T_b= \nabla^2 T_d + 2 \frac{Ge}{Ra} \boldsymbol{u}_b \boldsymbol{\cdot} \boldsymbol{u}_d , \end{equation}

which is subject to the linear and homogeneous boundary conditions

(2.18)\begin{equation} w_d = T_d = 0 \quad \text{at}\ z = 0 \end{equation}

and

(2.19)\begin{equation} w_d = T_d = 0 \quad \text{at}\ z = 1 , \end{equation}

where viscous dissipation has a direct effect on the linear disturbances through the last term in (2.17) when $Ge>0$, in addition to the indirect one through the base flow.

Disturbances are further assumed to behave as normal modes in the homogeneous directions, namely $x$, $y$ and with respect to time $t$. In other words,

(2.20)\begin{equation} \{\boldsymbol{u}_d,T_d,P_d\} = \{\hat{\boldsymbol{u}}(z),\hat{T}(z),\hat{P}(z)\}\exp\left[{\rm i}(\alpha x + \beta y - \omega t) \right]+ \mathrm{c.c.}, \end{equation}

where c.c. denotes complex conjugate. Furthermore, $\alpha$ and $\beta$ are the streamwise and spanwise complex wavenumbers, respectively. Their real parts can be used to calculate the real wavelengths whereas their imaginary parts are the spatial damping rates. Finally, $\omega$ is the complex angular frequency. Its real part is the real angular frequency whereas its imaginary part is the temporal growth rate.

Substituting (2.20) for the normal modes into (2.15)–(2.19) for the linear disturbances leads to the differential eigenvalue problem

(2.21)\begin{gather} {\rm i}\,\alpha\hat{u}(z) + {\rm i}\,\beta\hat{v}(z) + \hat{w}'(z) = 0 , \end{gather}

(2.22a,b)\begin{gather} {\rm i}\,\alpha\hat{P}(z) + \hat{u}(z) = 0 , \quad {\rm i}\,\beta\hat{P}(z) + \hat{v}(z) = 0 , \quad \hat{P}'(z) + \hat{w}(z) = Ra\,\hat{T}(z) , \end{gather}

(2.23)\begin{gather} {\rm i}(Pe\,\alpha-\omega)\hat{T}(z)+T_b'(z)\,\hat{w}(z)=\hat{T}''(z)-(\alpha^2 + \beta^2)\hat{T}(z)+2 \frac{Ge Pe}{Ra}\hat{u}(z) , \end{gather}

which is subject to the following normal mode boundary conditions:

(2.24)\begin{gather} \hat{w} = \hat{T} = 0 \quad \text{at}\ z = 0 , \end{gather}

(2.25)\begin{gather} \hat{w} = \hat{T} = 0 \quad \text{at}\ z = 1 . \end{gather}

It is convenient to rearrange this system into a single ordinary differential equation, which can be written in terms of the normal disturbance velocity as

(2.26)$$\begin{align}& \hat{w}''''(z) - {\rm i}(Pe\,\alpha-2\,{\rm i}(\alpha^2+\beta^2)-\omega)\hat{w}''(z)- 2\,{\rm i}\,\alpha\,Ge\,Pe\,\hat{w}'(z) \nonumber\\ &\quad + (\alpha^2+\beta^2)({\rm i}\,Pe\,\alpha+(\alpha^2+\beta^2) - {\rm i}\,\omega+Ra\,T_b'(z))\,\hat{w}(z)=0 , \end{align}$$

which is a fourth-order ordinary differential equation. Hence, it requires two additional boundary conditions. They can be obtained from the relation between temperature and normal velocity disturbances, derived from (2.21) and (2.22a,b) and given by

(2.27)\begin{equation} \hat{T}(z)=\frac{(\alpha^2+\beta^2)\hat{w}(z)-\hat{w}''(z)}{Ra(\alpha^2+\beta^2)} , \end{equation}

and, hence, (2.26) is subject to the following boundary conditions:

(2.28a,b)\begin{gather} \hat{w} = 0 \quad \text{and}\quad \hat{w}'' = 0 \quad \text{at}\ z = 0 , \end{gather}

(2.29a,b)\begin{gather} \hat{w} = 0 \quad \text{and}\quad \hat{w}'' = 0 \quad \text{at}\ z = 1 . \end{gather}

3. Analysis methodology

Equation (2.26) and its boundary conditions (2.28a,b) and (2.28a,b) are solved here using a matrix-forming approach. It transforms the ordinary differential equation into a system of algebraic equations that is recast as a generalised eigenvalue problem. The main difference between modal and non-modal analyses lies on which information they extract from this eigensystem to model linear disturbance growth. On one hand, a modal analysis evaluates the eigenvalues of the generalised eigenvalue problem. On the other hand, a non-modal analysis evaluates their respective eigenvectors. When doing so, the former models the behaviour of a single disturbance in the case of convective instability, and of a disturbance wavepacket in the case of absolute instability. On the same token, the latter models the behaviour of superposed disturbances.

One consequence of this process is that a modal analysis provides the time asymptotic behaviour of linear disturbances. In the case of a convectively unstable flow, an excitation source is capable of promoting the spatial growth of the targeted disturbance downstream of its location. In the case of an absolutely unstable flow, a disturbance present in the initial condition grows spatially both downstream and upstream of its original location. In other words, the former displays the extrinsic dynamics typical of a noise amplifier whereas the latter displays the intrinsic dynamics typical of an oscillator. Another consequence of this process is that a non-modal analysis provides initial transient behaviour induced by weighted disturbance superposition. Commonly known as transient growth, it can predict a temporary disturbance superposition energy growth even when each individual disturbance being superposed is time asymptotically stable. All the steps employed by each approach to the eigenvalue problem are described next.

3.1. Integral transform pair

A truncated series solution for the vertical velocity disturbance is first proposed in the form of the inverse function

(3.1)\begin{equation} \hat{w}(z)=\sum_{m=1}^{N_t}\tilde{w}_m \tilde{\psi}_m(z) , \end{equation}

where the number of terms $N_t$ in this summation series must be chosen high enough to guarantee a user-prescribed tolerance.

The orthogonal basis function is obtained from

(3.2)\begin{equation} \psi_m''''(z)=\lambda_m^4 \psi_m(z) , \end{equation}

which is a Sturm–Liouville-type problem subject to boundary conditions

(3.3)\begin{equation} \psi_m = \psi_m'' = 0 \quad \text{at} \ z = 0 \end{equation}

and

(3.4)\begin{equation} \psi_m = \psi_m'' = 0 \quad \text{at} \ z = 1 . \end{equation}

One can then define the orthonormal basis function

(3.5)\begin{equation} \tilde{\psi}_m(z)=\frac{\psi_m(z)}{\sqrt{N_m}}={-}\frac{\sinh(\lambda_m)\sin(\lambda_m z)}{\sqrt{N_m}} , \end{equation}

which satisfies the normalised versions of (3.2) and (3.3), as well as eigenvalues

(3.6)\begin{equation} \lambda_m=m{\rm \pi} , \end{equation}

which allows (3.5) to also satisfy the normalised version of (3.4), where the normalisation function is defined as

(3.7)\begin{equation} N_m=\frac{\cos(2\lambda_m)+\cosh(2\lambda_m)}{4} - \frac{1}{2} , \end{equation}

in order to guarantee (3.5) is orthonormal, i.e.

(3.8)\begin{equation} \int_0^1\tilde{\psi}_m(z)\tilde{\psi}_n(z)\,{\rm d}z=\delta_{m,n} , \end{equation}

where $\delta _{m,n}$ is the Kronecker delta.

Finally, multiplying (3.1) by the orthonormal eigenfunction in (3.5), integrating the result over the domain length and using (3.8) leads to

(3.9)\begin{equation} \tilde{w}_m=\int_0^1 \tilde{\psi}_m(z) \hat{w}(z) \,\mathrm{d}z , \end{equation}

which defines the integral transformed vertical velocity disturbance.

3.2. Matrix-forming approach

The procedure that transformed (3.1) into (3.9) can be applied to (2.26) as well. In other words, multiplying it by the orthonormal eigenfunction, integrating the result over the domain length and using the orthonormality condition leads to

(3.10) $$\begin{gather} \int_0^1 \tilde{\psi}_m(z)\hat{w}''''(z) \,{\rm d}z - {\rm i}\,\left(Pe \, \alpha - 2{\rm i}(\alpha^2+\beta^2)-\omega\right)\int_0^1 \tilde{\psi}_m(z)\hat{w}''(z) \,{\rm d}z \nonumber\\ -\,2{\rm i}\,\alpha \,Ge \,Pe \int_0^1 \tilde{\psi}_m(z)\hat{w}'(z)\,{\rm d}z - (\alpha^2+\beta^2)Ge\, Pe^2 \int_0^1 z\tilde{\psi}_m(z)\hat{w}(z)\,{\rm d}z \nonumber\\ +\,\tfrac{1}{2}(\alpha^2+\beta^2)\left(Ge \,Pe^2+2{\rm i}\left(Pe \, \alpha - \omega - {\rm i}(\alpha^2+\beta^2-Ra)\right)\right)\int_0^1 \tilde{\psi}_m(z)\hat{w}(z)\,{\rm d}z=0 , \end{gather}$$

after replacing the base flow terms, defined in (2.12)–(2.14), where the first term can be further simplified to yield

(3.11)\begin{equation} \int_0^1 \tilde{\psi}_m(z)\hat{w}''''(z)\,{\rm d}z=\int_0^1 \tilde{\psi}_m''''(z)\hat{w}(z)\,{\rm d}z=\lambda_m^4 \int_0^1 \tilde{\psi}_m(z)\hat{w}(z)\,{\rm d}z=\lambda_m^4 \tilde{w}_m , \end{equation}

after using the boundary conditions in (2.28a,b), (2.28a,b), (3.3) and (3.4) when integrating by parts, (3.2) to eliminate the eigenfunction derivative and (3.9). Equation (3.10) can now be simplified using (3.11), as well as the inverse/transform pair given by (3.1) and (3.9), respectively, to yield

(3.12)\begin{equation} \sum_{n=1}^{N_t} A_{m,n}\hat{w}_n = 0 \quad \mbox{or} \quad \boldsymbol{\mathsf{A}} \boldsymbol{\cdot} \hat{\boldsymbol{w}} = 0 , \end{equation}

where the coefficients $A_{m,n}$ of the matrix $\boldsymbol{\mathsf{A}}$ formed are defined as

(3.13) $$\begin{align} A_{m,n} &= \left(\lambda_m^4 + \tfrac{1}{2}(\alpha^2 + \beta^2)\left(Ge Pe^2 + 2{\rm i}\left(Pe \, \alpha - \omega - {\rm i}(\alpha^2 + \beta^2 - Ra)\right)\right) \right) \delta_{m,n} \nonumber\\ &\quad - (\alpha^2+\beta^2)\,Ge \,Pe^2 \, A_{m,n}^{(1)} - 2\,{\rm i}\,\alpha \,Ge \,Pe\, A_{m,n}^{(2)} - {\rm i}\,(Pe \, \alpha - 2\, {\rm i} \, (\alpha^2+\beta^2)-\omega)\,A_{m,n}^{(3)} , \end{align}$$

which depends on the integral transformed coefficient matrices

(3.14)$$\begin{gather} A_{m,n}^{(1)}=\int_0^1 z \hat{\psi}_m(z)\tilde{\psi}_n(z)\,\mathrm{d}z , \end{gather}$$

(3.15)$$\begin{gather}A_{m,n}^{(2)}=\int_0^1 \hat{\psi}_m(z)\tilde{\psi}_n'(z)\,\mathrm{d}z \quad \text{and} \end{gather}$$

(3.16)$$\begin{gather}A_{m,n}^{(3)}=\int_0^1 \hat{\psi}_m(z)\tilde{\psi}_n''(z)\,\mathrm{d}z , \end{gather}$$

whose integrals can be obtained analytically.

3.3. Non-modal analysis

Extending the inner product between two functions $u(z)$ and $v(z)$ used earlier to define the real orthogonal eigenfunctions in (3.8) towards complex functions, i.e.

(3.17)\begin{equation} \langle u,v \rangle = \int_0^1 u\, v^{{\ast}} \, \mathrm{d}z , \end{equation}

where $\ast$ denotes the complex conjugate, enables one to write

(3.18)\begin{equation} \langle L w,\xi\rangle = \langle w,\mathcal{L}\xi \rangle , \end{equation}

which defines $\xi$, the adjoint of $w$. In the above integral relation, $L$ is the linear operator associated with (2.26) and, hence, $\mathcal {L}$ is the respective adjoint linear operator. It is obtained using integration by parts to derive the right-hand side of (3.18) from its left-hand side while imposing appropriate boundary conditions to maintain homogeneity. Doing so yields

(3.19)$$\begin{align} &\hat{\xi}''''(z) - {\rm i}\{Pe\,\alpha-2\,{\rm i}(\alpha^2+\beta^2)-\omega\}\hat{\xi}''(z)+ 2{\rm i}\,\alpha\,Ge\,Pe\,\hat{\xi}'(z) \nonumber\\ &\quad + (\alpha^2+\beta^2)\{{\rm i}\,Pe\,\alpha+(\alpha^2+\beta^2) - {\rm i}\,\omega+Ra\,T_b'(z)\}\,\hat{\xi}(z)=0 , \end{align}$$

which means operators $L$ and $\mathcal {L}$ are identical, except for their third term having opposite signs. In other words, (2.26) is no longer self-adjoint when both viscous dissipation ($Ge>0$) and throughflow ($Pe>0$) are present. This implies that transient growth is indeed possible.

In order to quantify transient growth, an energy metric must be defined. The most common choice for incompressible isothermal flows is the kinetic energy. Otherwise, when temperature gradients become relevant, one can use instead

(3.20)\begin{equation} E(t) = \int_0^1 \left(\sigma(|\hat{u}|^2 + |\hat{v}|^2 + |\hat{w}|^2) + \gamma|\hat{T}|^2 \right)\, \mathrm{d}z , \end{equation}

where $\sigma$ and $\gamma$ are arbitrary positive scalars. Nevertheless, it is always possible to prescribe one of the constants, e.g. $\sigma =1$, since only relative growth measures are important. Furthermore, even though $\gamma$ has a quantitative impact on the energy metric, such an impact is usually not significant (Hanifi, Schmid & Henningson 1996; Biau & Bottaro 2004; Sameen & Govindarajan 2007). The same lack of sensitivity with respect to $\gamma$ was noted in the present problem, so the results shown here use $\gamma =1$ as well. Finally, the inner product defined in (3.17) can be associated with an energy norm as follows:

(3.21)\begin{equation} E(t)=\langle \boldsymbol{q}, \boldsymbol{q} \rangle=\|\boldsymbol{q}(z,t)\|_E^2 , \end{equation}

where $\boldsymbol {q}(z,t)=\{\hat {u},\hat {v},\hat {w},\hat {T}\}^{tr}$ and superscript $tr$ represents the transpose.

It is now possible to define the gain as

(3.22)\begin{equation} G(t) = \max_{\boldsymbol{q}_0 \neq 0} \left(\frac{E(t)}{E(0)}\right) =\max_{\boldsymbol{q}_0 \neq 0} \frac{\|\boldsymbol{q}(z,t)\|_E^2}{\|\boldsymbol{q}(z,0)\|_E^2} , \end{equation}

which represents the maximum possible growth at a given time $t$ over all possible initial conditions $\boldsymbol {q}_0 = \boldsymbol {q}(0)$. The state vector $\boldsymbol {q}$ is defined as a linear combination of infinite eigenvectors. However, this infinite summation has to be truncated for numerical reasons. Hence, the state vector $\boldsymbol {q}$ must be approximated by

(3.23)\begin{equation} \boldsymbol{q}(z,t) \simeq \boldsymbol{\tilde{q}}_i(z) \, k_i(t) , \end{equation}

using Einstein summation notation for a repeated index, where $i=1, 2, 3, \ldots, N$, with $N$ characterising the truncation of the infinite summation. Furthermore, $k_i(t)=k_i(0)\mathrm {e}^{-{\rm i} \omega _i t}$ represents the weight of each eigenvector $\boldsymbol {\tilde {q}}_i$ on the final composition of the state vector $\boldsymbol {q}$, as well as the time evolution of each eigenvector $\boldsymbol {\tilde {q}}_i$ given by its associated eigenvalue $\omega _i$. Substituting (3.23) into (3.21) yields

(3.24)$$\begin{align} \|\boldsymbol{q}(z,t)\|_E^2 &= \int \boldsymbol{q}^{{\ast}}(z,t) \boldsymbol{q}(z,t) \, \mathrm{d}z = \int \left(\boldsymbol{\tilde{q}}_j^{{\ast}}(z) k_j^{{\ast}}(t) \right) \left( \boldsymbol{\tilde{q}}_k(z) \, k_k(t) \right) \, \mathrm{d}z \nonumber\\ &= k_j^{{\ast}}(t) M_{j,k} k_k(t) = \boldsymbol{k}^{{\ast}}(t) \boldsymbol{M} \boldsymbol{k}(t) , \end{align}$$

where the matrix $\boldsymbol {M}$ coefficients are given by

(3.25)\begin{equation} M_{j,k} = \int \boldsymbol{\tilde{q}}_j^{{\ast}}(z) \boldsymbol{\tilde{q}}_k(z) \, \mathrm{d}z , \end{equation}

and $\boldsymbol {k} = \{ k_1, k_2,\ldots,k_N \}$. Since $\boldsymbol {M}$ is a positive-definite Hermitian matrix, it can be decomposed as $\boldsymbol {M}=\boldsymbol {F}^{{\dagger} } \boldsymbol {F}$, where $\boldsymbol {F}^{{\dagger} }$ is the adjoint (conjugate transpose) of $\boldsymbol {F}$. It is possible to transform the energy norm in (3.21) into

(3.26)\begin{equation} \|\boldsymbol{q}(z,t)\|_E^2 = \boldsymbol{k}^{{\ast}}(t) \boldsymbol{M} \boldsymbol{k}(t) = \boldsymbol{k}^{{\ast}}(t) \boldsymbol{F}^{{{\dagger}}} \boldsymbol{F} \boldsymbol{k}(t) = \|\boldsymbol{F} \boldsymbol{k}(t)\|^2_2 , \end{equation}

namely an Euclidean norm, leads to the new expression for the gain

(3.27)\begin{equation} G(t) = \max_{\boldsymbol{k}_0\neq0} \frac{\|\boldsymbol{F} \boldsymbol{k}(t)\|^2_2}{\|\boldsymbol{F} \boldsymbol{k}(0)\|^2_2} = \max_{\boldsymbol{k}_0\neq0} \frac{\|\boldsymbol{F} \varLambda \boldsymbol{k}(0)\|^2_2}{\|\boldsymbol{F} \boldsymbol{k}(0)\|^2_2} = \max_{\boldsymbol{k}_0\neq0} \frac{\|\boldsymbol{F} \varLambda \boldsymbol{F}^{{-}1}\boldsymbol{F} \boldsymbol{k}(0)\|^2_2}{\|\boldsymbol{F} \boldsymbol{k}(0)\|^2_2}, \end{equation}

where $\boldsymbol {k}_0 = \boldsymbol {k}(0)$ and $\varLambda = \mathrm {diag}({\rm e}^{-{\rm i} \omega _1 t},{\rm e}^{-{\rm i} \omega _2 t},\ldots,{\rm e}^{-{\rm i} \omega _N t})$. Hence, the gain can be optimised over all initial conditions at each time $t$ by solving the matrix norm

(3.28)\begin{equation} G(t) = \|\boldsymbol{F} \varLambda \boldsymbol{F}^{{-}1}\|^2_2 , \end{equation}

where the superscript $-1$ means inverse. In other words, (3.28) provides the maximum energy growth at a given time $t$ for any given pair $\alpha$ and $\beta$. An important feature of the formula given by (3.28) is that it can be easily determined by means of a singular value decomposition (SVD), as it is always true for the Euclidean norm of a matrix. If this gain is large enough for a given initial disturbance amplitude, it will likely trigger a subcritical transition towards a more complex flow pattern.

Since the eigenvectors become orthogonal in the absence of viscous dissipation and throughflow, it is interesting to note that (3.28) reduces to

(3.29)\begin{equation} G(t) = \| \varLambda \|^2_2 = {\rm e}^{2 \mathrm{Im}\left[\omega\right]_{max} t} , \end{equation}

because $\boldsymbol {M}$ and $\boldsymbol {F}$ become diagonal matrices, where $\mathrm {Im}[\omega ]_{max}$ is the imaginary part of the least stable (or most unstable) eigenvalue.

3.4. Absolute instability analysis

In order to identify the transition to absolute instability, one must investigate the behaviour of a disturbance wavepacket in the limit of very large times ($t \to \infty$). If the analysis is restricted to a two-dimensional wavepacket, one must evaluate an integral on $\alpha$ over a path $\gamma$, which coincides with the infinite real domain $\alpha \in (-\infty,\infty )$. Its temporal behaviour for $t \to \infty$ can be given by the largest growth rate on the saddle point of $\omega$ on $\alpha _0$ ($\partial \omega / \partial \alpha =0$ at $\alpha = \alpha _0$). This conclusion cames from the steepest-descent approximation, which requires the wavepacket to be holomorphic and the paths $\gamma$, coincident to the real axis of $\alpha$, and $\gamma ^*$, crossing the saddle point $\alpha _0$, to be homotopic. In other words, it must be possible to continuously deform $\gamma$ into $\gamma ^*$ in order to apply this approximation. As already pointed out by Barletta (2019), if there are multiple saddle points $\alpha _0$, the steepest-descent approximation just keeps those with largest Im$[\omega ]$. In the case such saddle points share the same value of Im$[\omega ]$, their contributions have to be summed up in order to apply the steepest-descent approximation. The aforementioned considerations are based on a two-dimensional wavepacket. However, as pointed out by Brevdo (1991), the same should be true for three-dimensional wavepackets, namely the asymptotic behaviour of the wavepacket can be given by looking at Im$[\omega ]$ on the saddle points of $\omega$ on $\alpha _0$ ($\partial \omega / \partial \alpha =0$ at $\alpha = \alpha _0$) and $\beta _0$ ($\partial \omega / \partial \beta = 0$ at $\beta = \beta _0$).

Identifying the transition to absolute instability can be computationally quite intensive when employing classical techniques, e.g. finding the steepest descent curve or verifying the collision criterion, unless saddle points can be cheaply calculated a priori (Alves et al. 2019). In the present case, this can be done by applying the zero group velocity conditions to the dispersion relation, coupling it with auxiliary dispersion relations that can identify saddle points. They are, however, a necessary but not sufficient condition for absolute instability. Once these points have been found, one must either obtain a steepest descent curve or verify the collision criterion in order to make sure they are, in fact, pinching points (Barletta 2019).

In order to find the aforementioned saddle points, one must first note that (3.12) only has non-trivial solutions when

(3.30)\begin{equation} \det \left(\boldsymbol{\mathsf{A}} \right)=0 , \end{equation}

for a fixed value of $N_t$. This equation is, in fact, the dispersion relation for this problem. It must then be coupled with the additional equations

(3.31a,b)\begin{equation} \frac{\partial \det \left( \boldsymbol{\mathsf{A}} \right)}{\partial \alpha}=0 \quad\mathrm{and}\quad \frac{\partial \det \left( \boldsymbol{\mathsf{A}} \right)}{\partial \beta}=0 , \end{equation}

respectively, restricted by the zero group velocity conditions

(3.32a,b)\begin{equation} \frac{\partial \omega}{\partial \alpha}=0\quad\mathrm{and}\quad \frac{\partial \omega}{\partial \beta}=0 , \end{equation}

to provide the auxiliary dispersion relations for this problem. Together, they form a set of three complex equations that yields the saddle points in the complex $\alpha$ and $\beta$ planes and their complex frequency $\omega$ for a set of prescribed control parameter values. The reader is referred to recent in depth reviews for more information about absolute instability calculations (Alves et al. 2019; Barletta 2019).

4. Results and discussion

4.1. Code verification

Before modal and non-modal linear stability analyses are employed to investigate the effects of viscous dissipation on the asymptotic and transient disturbance growth in time, respectively, the codes developed for these analyses are verified under four different scenarios. The first two involve the modal onset of (convective) instability. In the absence of viscous dissipation, which is the first scenario, the present model reduces to the classical problem originally solved by Prats (1966). Furthermore, the classical critical Rayleigh number for the onset of instability when $Pe=0$, i.e. $Ra=4 {\rm \pi}^2$, is recovered even when truncating (3.12) with $N_t=1$. This simplified dispersion relation is given by $4 \beta (2 {\rm i} {\rm \pi}^2 - {\rm i} Ra - \alpha Pe + 2 {\rm i} (\alpha ^2 + \beta ^2) + \omega ) = 0$ when $Pe\neq 0$, which suggests that (convective) instability first appears with $\beta =0$ unless $\omega =\alpha \,Pe$, since $\alpha _i=\beta _i=\omega _i=0$. Numerical evaluations of the converged ($N_t\gg 1$) dispersion relation indicate that this is indeed the case. It turns out this is also true in the presence of viscous dissipation, which is the second scenario. Figure 1(a) shows the critical Rayleigh number for the onset of (convective) instability as a function of the Péclet number obtained from the present code (line) when $Ge=0.1$ and from Nield et al. (2011) (points). There is a very good agreement between both sets of results. The third scenario considered for verification purposes is the modal onset of absolute instability in the absence of viscous dissipation. Figure 1(b) shows the critical Rayleigh number for the onset of absolute instability as a function of the Péclet number obtained from the present code (line) when $Ge=0.0$ and from Hirata & Ouarzazi (2010) (points). Once again, a very good agreement is observed between both sets of results. Furthermore, this onset also occurs for $\beta =0$. The fourth and final verification scenario considers non-modal growth. Since direct and adjoint linear disturbance governing (2.26) and (3.19) are identical in the absence of viscous dissipation, i.e. when $Ge=0$, transient growth cannot occur. Figure 2 compares the gain temporal behaviour calculated by two metrics for $Pe=100$, $Ra=39$, $\alpha ={\rm \pi}$ and $\beta =0$ when $Ge=0$. One metric is given by (3.29), which calculates the linear gain based on the dominant eigenvalue (red dashed lines). The other metric is given by (3.28), which calculates the gain associated with the entire linear dynamics, i.e. based on both eigenvalues and eigenvectors (black dotted lines). The modal behaviour is manifested in the exponential temporal evolution of the gain, which is recovered as a straight line in semi-logarithmic plots. Since this figure shows that both curves are essentially identical straight lines, it indicates that transient growth indeed does not occur. Finally, summation series convergence studies are presented in Appendix A.

Figure 1. Critical Rayleigh number for the onset of convective (a) and absolute (b) instability as a function of the Péclet number when $Ge=0.1$ and $0.0$, respectively. Lines represent current results whereas points represent results from the literature for the former (Nield et al. 2011) and latter (Hirata & Ouarzazi 2010) cases, respectively.

Figure 2. Gain as a function of time for $Pe=100$, $Ra=39$, $\alpha ={\rm \pi}$, $\beta =0$ and $Ge=0$, obtained from (red dashed) (3.29) and (black dotted) (3.28).

4.2. Transient disturbance growth

Non-modal growth in the presence of viscous dissipation is considered next. Since direct and adjoint linear disturbance governing (2.26) and (3.19) are not identical when $Ge>0$ and $Pe>0$, transient growth is possible. Gain calculations under several different parametric conditions, however, do not reveal any significant transient growth. Figure 3 provides evidence in favour of this claim by showing the same results from figure 2 but for the particular cases where $Pe=200$, $Ra=-800$ and $Ge=0.1$ with $\alpha =1$, $2$, $3$ and $4$. Differences between results from (3.28) and (3.29) are barely noticeable. Hence, the behaviour of all four cases is essentially modal, with the former three ($\alpha =1$, $2$ and $3$) being stable and the latter ($\alpha =4$) being unstable. The condition number of the eigenvector matrix can quantify this dominance of the modal behaviour, or weakness of the non-modal one, since this number is equal to one (infinity) when the eigenvectors are orthogonal (parallel). Present calculations show that its value is indeed 1 when $Ge=0$ and remains at ${O}(1)$ when $Ge>0$. Table 1 provides similar evidence through the condition number of the eigenvector matrix for different Rayleigh, Péclet and streamwise wavenumbers when $Ge=0.1$. Even though the aforementioned results were obtained for $\beta =0$, the same trends were observed for non-zero $\beta$ as well. Hence, one can conclude that transient growth is not relevant in this problem.

Figure 3. Gain vs time for $Pe=200$, $Ra=-800$ and $Ge=0.1$ for different $\alpha$ with $\beta =0$, obtained from (red dashed) (3.29) and (black dotted) (3.28).

Table 1. Condition number of the eigenvector matrix for $\beta =0$ and $Ge=0.1$.

4.3. Asymptotic disturbance growth

4.3.1. Without internal heating

Since transient growth cannot induce a disturbance amplitude increase in time that is strong enough to cause a transition from stability to instability, the only other known linear mechanism that can do so is asymptotic growth. Hence, this section presents results regarding the influence of viscous dissipation on the transition from convective to absolute instability. The onset of the former has been investigated for mixed convection within a porous medium in the presence of viscous dissipation (Nield et al. 2011), but the onset of the latter has only been investigated for this problem in the absence of viscous dissipation (Hirata & Ouarzazi 2010). Nevertheless, before doing so, it is important to review some key aspects of this asymptotic growth when $Ge=0$. Figure 4(a) shows the critical Rayleigh numbers for the onset of convective (dashed line) and absolute (solid line) instability as a function of the Péclet number obtained using the present code. The former is Péclet independent, as expected. The latter, on the other hand, shows that an increasing Péclet number has a stabilising effect. This is a typical result (Barletta 2019). More energy must be provided to the wavepacket for it to overcome a convectively stronger base flow and propagate upstream, which is obtained from the external heat source by increasing the Rayleigh number. Figure 4(b) shows the collision criterion for the particular case with $Pe=50$, showing that the saddle point found using (3.30)–(3.31a,b) is indeed a pinching point. Although not shown here, the same trends were observed at several other Péclet numbers. Hence, there is enough evidence supporting the claim that all saddle points shown in figure 4(a) are indeed pinching points. Two final remarks deserve special mention for their relevance to subsequent results. First, extensive numerical simulations found no additional saddle points competing for the role of a pinching point. Second, all pinching points found have $\beta =0$.

Figure 4. (a) Critical Rayleigh numbers for the onset of convective (dashed line) and absolute (solid line) instability as a function of the Péclet number for $\beta =0$ and $Ge=0$. (b) Collision criterion (coloured lines) just before (blue), at (orange) and just after (green) the pinching point (large dot) for $\beta =0$, $Ge=0$, $Pe=50$ and $\omega _i=0$, with $(\alpha _r,\alpha _i) = (2.88674,-3.12176)$, where the dashed line represents the steepest ascent path followed by the pinching point since the onset of convective instability.

4.3.2. With internal heating: small Péclet numbers

The influence of viscous dissipation on the asymptotic disturbance behaviour is now discussed, focusing first on small Péclet numbers, namely $Pe\leq 50$. Table 2 presents critical Rayleigh number, real frequency and complex wavenumber for a few given Péclet and Gebhart numbers at the onset of absolute instability. On one hand, throughflow still has a stabilising effect on the transition to absolute instability in the presence of viscous dissipation. On the other hand, viscous dissipation has a destabilising effect on this transition in the presence of throughflow. This is due to the fact that it acts as an internal heating mechanism and, hence, less external heating is required. Further insight can be gained using (3.13). According to this equation, throughflow effects are approximately $O(Pe)$ whereas viscous dissipation effects are approximately $O(Ge\,Pe^2)$, which explains why viscous dissipation effects only become relevant at large $Ge$ for these moderate $Pe$. Such a dimensional analysis also indicates that throughflow will eventually have a destabilising effect when it is large enough.

Table 2. Throughflow ($Pe$) and viscous dissipation ($Ge$) influence on the transition to absolute instability ($\omega _i=0$).

Figure 5. Onsets of convective (small dashed line) and absolute instabilities for $Ge=0.1$ and $Pe\leq 450$, showing saddle points (long dashed lines) competing for the role of pinching point (solid line). Red dots labelled (a–f) are evaluated at the critical points $(Pe,Ra)=(10,57.7822)$, $(100,413.001)$, $(206.2,857.197)$, $(300,728.532)$, $(373.2,402.244)$ and $(420,-58.0117)$, respectively.

Before proceeding any further, it is important to remind the reader about the Gebhart number magnitude. In most engineering applications, $Ge\ll 1$. However, it is possible to reach $Ge \sim O(1)$ in some geophysical flows. This is the reason why such a range was used in table 2. Nevertheless, for the purposes of the present study, $Ge=0.1$ is assumed to be a reasonable upper bound in the following discussion.

4.3.3. With internal heating: large Péclet numbers

Focus is now switched to a larger Péclet number range, i.e. $50 \le Pe \leq 450$. It is important to remind the reader that the Péclet number is the product of the Prandtl and Reynolds numbers. Since the former can have quite high values even for moderate values of the latter, such a large Péclet number range is indeed not uncommon. Furthermore, the influence of viscous dissipation on the time asymptotic disturbance behaviour is discussed in detail from this point for $Ge=0.1$, which has already been established as a reasonable upper bound for this parameter. Nevertheless, the disturbance behaviour at smaller $Ge$ values is summarised at the end of this subsection.

Figure 5 shows the onsets of convective (small dashed line) and absolute (solid line) instabilities over the aforementioned $(Ge,Pe)$ parametric region, where three different saddle points (long dashed lines) compete for the role of pinching point. Their respective collision criterions (solid lines) are provided in figure 6 for the six $(Pe,Ra)$ critical points (red dots) shown in figure 5. The paths taken by all three saddle points in the complex streamwise wavenumber plane as the Péclet number varies between $0\leq Pe \leq 450$ are also shown (non-solid lines) in figure 6. As already illustrated in figure 1, it is important to note that the onset of convective instability is always destabilised by throughflow, so much so that the flow becomes convectively unstable even when heated from above for $Pe\geq 68.6575$ (Nield et al. 2011). The effect of throughflow on the onset of absolute instability, on the other hand, is significantly less straightforward. A single downstream propagating branch $\alpha ^+$ is involved in the collision that forms the pinching point for all Péclet numbers within $0\leq Pe \leq 450$. This collision, however, occurs against one of three different upstream propagating branches, namely $\alpha _1^{-}$, $\alpha _2^{-}$ and $\alpha _3^{-}$, within the same Péclet number range. Consider first scenario (a) from figures 5 and 6. This particular pinching point is formed by a collision between branches $\alpha ^+$ and $\alpha _1^{-}$ and called here mode 1. The other two saddle points are formed by collisions between upstream propagating branches and, hence, do not have any physical meaning. Mode 1 is the same transition mechanism observed in the absence of viscous dissipation, as shown in figure 4, where the latter two saddle points do not exist. A similar transition mechanism is observed in scenario (b) from figures 5 and 6. The only qualitative difference observed at this higher Péclet number is related to the wavenumber of the saddle points induced by viscous dissipation, whose real part is no longer zero. A further increase in the throughflow magnitude to $Pe=206.2$, however, leads to the emergence of a new transition mechanism. Scenario (c) from figures 5 and 6 show a codimension-two point where the traditional (convection dominated) and novel (viscous dissipation dominated) saddle points are equally unstable. In other words, branch $\alpha ^+$ collides simultaneously with branches $\alpha _1^{-}$ and $\alpha _2^{-}$. The latter collision, however, is called mode 2 here and is the one that leads to the pinching point at higher throughflow, as shown in scenario (d) from figures 5 and 6. Increasing the throughflow even further to $Pe\simeq 373.2$, on the other hand, leads to a new codimension-two point, but now between both novel viscous dissipation dominated saddle points, as shown in scenario (e) from figures 5 and 6. In other words, branch $\alpha ^+$ collides simultaneously with branches $\alpha _2^{-}$ and $\alpha _3^{-}$. This second viscous dissipation related collision is called here mode 3 and leads to the pinching point at higher throughflow, as shown in scenario ( f) from figures 5 and 6. There are more viscous dissipation related saddle points, but they only become relevant at higher Gebhart numbers, which are non-physical. Hence, they are not discussed. Finally, it is important to note that a linear stability analysis cannot provide additional information about either codimension-two point. Under either parametric condition, a nonlinear analysis is required to clarify the dynamical system behaviour.

Figure 6. Collision checks for the parametric conditions shown in figure 5 at (a) $Pe=10$, where $(\alpha _r,\alpha _i) = (3.39332,-1.89382)$, (b) $Pe=100$, between branches $\alpha ^+$ and $\alpha _1^{-}$, where $(\alpha _r,\alpha _i) = (2.54770,-3.41712)$, (c) $Pe\simeq 206.2$ between branches $\alpha ^+$ and $\alpha _1^{-}$, where $(\alpha _r,\alpha _i) = (1.92771,-3.81869)$ and branches $\alpha ^+$ and $\alpha _2^{-}$, where $(\alpha _r,\alpha _i) = (3.14560,-6.65050)$, (d) $Pe=300$, between branches $\alpha ^+$ and $\alpha _2^{-}$, where $(\alpha _r,\alpha _i) = (2.68895,-6.76400)$, (e) $Pe\simeq 373.2$, between branches $\alpha ^+$ and $\alpha _2^{-}$, where $(\alpha _r,\alpha _i) = (2.15204,-7.06621)$, and branches $\alpha ^+$ and $\alpha _3^{-}$, where $(\alpha _r,\alpha _i) = (2.70854,-8.96509)$, and ( f) $Pe=420$, between branches $\alpha ^+$ and $\alpha _3^{-}$, where $(\alpha _r,\alpha _i) = (2.76478,-8.93637)$. The non-solid curves represent the paths taken by each saddle point (i.e. modes 1, 2 and 3) as the Péclet number varies.

Figures 7 and 8 show the disturbance streamlines and isotherms at the codimension-two points $Pe \simeq 206.2$ and $\simeq 373.2$, respectively. Dashed (solid) lines stand for negative (positive) values in the streamline plots. Modes 1, 2 and 3 are clearly distinct from one another. Similarly to the typical cell pattern behaviour observed in convectively unstable conditions, higher order modes have a larger number of cells. They are concentrated in the downstream end of both figures because all three modes grow in space at the onset of absolute instability. Finally, it is also interesting to note that the cell pattern is closer to horizontal for mode 1 but closer to vertical for mode 3. This is related to the fact that the spatial growth rates are different for each mode, as shown in figure 6. In figure 7, $\alpha _i = -3.81869$ for mode 1 but $\alpha _i = -6.65050$ for mode 2. In figure 8, $\alpha _i = -7.06621$ for mode 2 and $\alpha _i = -8.93637$ for mode 3.

Figure 7. Stream function (a,c) and temperature (b,d) isolines for mode 1 (a,b) and mode 2 (c,d) at $Pe\simeq 206.2$, i.e. at the first codimension-two point.

Figure 8. Stream function (a,c) and temperature (b,d) isolines for mode 2 (a,b) and mode 3 (c,d) at $Pe\simeq 373.2$, i.e. at the second codimension-two point.

When $Ge=0.1$, as shown in figure 5, the transition from convective to absolute instability for $Pe_C \simeq 414.651$ occurs when $Ra=0$. This means that transition is induced by internal heating alone, i.e. without external heating. In other words, viscous dissipation alone, without the influence of buoyancy effects, is responsible for the onset of absolute instability. Furthermore, when $Pe \ge Pe_C$, absolute instability occurs even in the presence of negative Rayleigh numbers. Under these conditions, a stable temperature stratification would be induced with $Ge=0$. Hence, the internal heating mechanism created by viscous dissipation is capable of inducing transition even in the presence of the stabilising external heating mechanism created by buoyancy. It turns out this critical Péclet number is induced by mode 3 and depends on the Gebhart number according to $Pe_C \simeq 41.9479/Ge$, since the mode 2 dependence is given by $Pe_C \simeq 43.3805/Ge$. Both correlations were obtained using nonlinear regression based on the data provided in tables 3 and 4, which also provides eigenvalues. The standard errors of both parameters are $2.32605\times 10^{-3}$ and $1.35113\times 10^{-3}$, respectively. These correlations imply that $Pe_C \rightarrow \infty \, (0)$ when $Ge \rightarrow 0 \, (\infty )$, which is expected since more (less) throughflow is required to make viscous dissipation important when the Gebhart number decreases (increases). Finally, this correlation also implies that the qualitative trends shown in figures 5–8 would be the same for different Gebhart numbers, although occurring at different Rayleigh and Péclet numbers, which has been confirmed but is not shown.

Table 3. Saddle point data at the onset of absolute instability when $Ra=0$ as a function of the Gebhart number for mode 2.

Table 4. Same as table 3 but for mode 3.

5. Conclusions

The present paper has investigated the appearance of natural convection, as induced by temporal disturbance growth, in the otherwise forced convection in porous media. At any given throughflow (Péclet number), instability could be induced by internal (Gebhart number) and/or external (Rayleigh number) heating. Both modal (asymptotic) and non-modal (transient) linear mechanisms have been considered. A matrix-forming approach based on a continuous spectral method has been employed to solve the differential eigenvalue problem. On one hand, the modal analysis has been performed using the dispersion relation obtained from the determinant of the resulting matrix. Auxiliary dispersion relations were then obtained by applying the zero group velocity conditions to verify the possibility of modal growth in time through absolute instability. On the other hand, the non-modal analysis has been performed using both the SVD of this matrix as well as the condition number of the respective eigenvector matrix. They allowed us to verify the possibility of transient growth. The major findings of our study are summarised as follows:

1. (i) In the absence of viscous dissipation, i.e. when the Gebhart number is zero, the differential eigenvalue problem is self-adjoint. However, this is no longer true for positive Gebhart numbers. Hence, non-modal growth is possible.

2. (ii) Non-modal growth was found negligible for a wide range of Péclet, Gebhart and Rayleigh numbers. Hence, it is possible to infer that there is no transient growth. This is true for two- and three-dimensional disturbances.

3. (iii) In the absence of internal heating, modal growth does occur for strong enough external heating, although throughflow has a stabilising effect. In the presence of internal heating, however, this is only true for a weak enough throughflow.

4. (iv) For a moderate throughflow, internal heating drives modal growth, where throughflow has a destabilising effect. Absolute instability is possible even in the absence of external heating, i.e. at zero Rayleigh number.

5. (v) For strong throughflows, internal heating is capable of inducing modal growth even in the presence of external cooling from below. In other words, absolute instability occurs for negative Rayleigh numbers.

6. (vi) In the three cases mentioned above, namely weak, moderate and strong throughflows, modal growth is controlled by mode 1, 2 and 3, respectively.

7. (vii) Two codimension-two points separate these three cases when $Ge=0.1$. The first one, separating weak and moderate throughflows, is located at $Pe \simeq 206.2$ whereas the second one, separating moderate and strong throughflows, is located at $Pe \simeq 373.2$.

8. (viii) Internal heating due to viscous dissipation effects is responsible for the appearance of modes 2 and 3, becoming the dominant mechanism over external heating responsible for the onset of absolute instability at moderate and strong throughflows.

9. (ix) Transverse modes are responsible for all modal growths discussed here. This is true when either internal or external heating acts as the dominant mechanism.

Current research is simulating the fully nonlinear system of governing equations to investigate both codimension-two points, which mark the switch between external and internal heating dominated absolute instabilities, as well as the switch between both internal heating dominated absolute instabilities, as throughflow increases. Furthermore, the possibility of a spatial non-modal growth is also being investigated.