1. Introduction
Convective instability induced by thermal gradients is a subject widely explored in the literature. The typical set-up giving rise to unstable behaviour is one where the vertical temperature gradient is directed downward. Such a configuration may pertain to a motionless basic state, as in the classical Rayleigh–Bénard problem and its many variants (Drazin & Reid Reference Drazin and Reid2004), as well as to basic forced or mixed convection flow states.
Much research has been done over the past sixty years to investigate thermal instability in fluid-saturated porous media. Surveys of current knowledge in this field are available in Nield & Bejan (Reference Nield and Bejan2013), as well as in Rees (Reference Rees, Vafai and Hadim2000), Tyvand (Reference Tyvand, Ingham and Pop2002) and Barletta (Reference Barletta, Öchsner and Murch2011). The thermoconvective instability of a basic motionless state, much like a porous medium version of the Rayleigh–Bénard problem, was first studied by Horton & Rogers (Reference Horton and Rogers1945) and Lapwood (Reference Lapwood1948). These studies were related to a horizontal layer with impermeable and isothermal walls kept at different temperatures, and they defined what is now well known as the Darcy–Bénard problem.
A direct extension of the Darcy–Bénard problem arises when the plane porous layer is inclined to the horizontal. The pioneering studies on this subject are the papers by Bories & Combarnous (Reference Bories and Combarnous1973), Weber (Reference Weber1975) and Caltagirone & Bories (Reference Caltagirone and Bories1985). The basic set-up considered by these authors is a porous medium modelled by Darcy’s law, and confined between a pair of impermeable walls kept at different uniform temperatures. The main consequence of the layer inclination is that the basic state cannot be motionless, but it is given by a stationary and parallel buoyant flow with a zero mass flow rate. Further developments have been achieved more recently by several authors (Karimi-Fard, Charrier-Mojtabi & Mojtabi Reference Karimi-Fard, Charrier-Mojtabi and Mojtabi1999; Storesletten & Tveitereid Reference Storesletten and Tveitereid1999; Rees & Bassom Reference Rees and Bassom2000; Rees & Postelnicu Reference Rees and Postelnicu2001; Rees, Storesletten & Postelnicu Reference Rees, Storesletten and Postelnicu2006; Barletta & Storesletten Reference Barletta and Storesletten2011; Nield, Barletta & Celli Reference Nield, Barletta and Celli2011; Rees & Barletta Reference Rees and Barletta2011; Barletta & Rees Reference Barletta and Rees2012). Many different physical effects have been studied, including the anisotropy of the porous medium, the viscous dissipation, and the modelling of the layer walls as isoflux, instead of isothermal, boundaries. In all these studies, with the single exception of Barletta & Rees (Reference Barletta and Rees2012), the thermal boundary conditions were such that the temperature in the basic state is streamwise-invariant, with no net heating or cooling effects.
Linear stability analysis has been employed as a useful methodology in determining the onset of convective stability for many of the above-mentioned studies. This methodology generally leads to a complex differential eigenvalue problem that needs to be solved for the Rayleigh number and wavenumber values, in order to determine the condition when a given flow becomes unstable. While a variety of techniques are available for this purpose, many of which are of a numerical nature, a hybrid analytical–numerical methodology known as the generalised integral transform technique (GITT) (Cotta Reference Cotta1990, Reference Cotta1993, Reference Cotta1994) can be used as a powerful tool for solving differential eigenvalue problems that arise in linear stability analysis. The solution of differential eigenproblems via the GITT, as described in Mikhailov & Cotta (Reference Mikhailov and Cotta1994), has been employed for different convection and diffusion problems, such as multidimensional heat diffusion in irregular geometries (Sphaier & Cotta Reference Sphaier and Cotta2000, Reference Sphaier and Cotta2002), diffusion in heterogeneous media (Naveira-Cotta et al. Reference Naveira-Cotta, Cotta, Orlande and Fudym2009), conjugate convection–conduction in a single channel (Knupp, Naveira-Cotta & Cotta Reference Knupp, Naveira-Cotta and Cotta2012), and conjugate convection–conduction in multiple micro-channels (Knupp, Naveira-Cotta & Cotta Reference Knupp, Naveira-Cotta and Cotta2014).
Figure 1. Sketch of the system.
The aim of this paper is to go beyond the analysis of the Darcy–Bénard problem in an inclined porous channel by devising a set-up where both walls are impermeable and symmetrically heated or cooled. The existence of a stationary solution is ensured, in this case, only if a net mass flow rate is allowed, so that the excess heat is removed or supplied by convection along the streamwise direction. Hence the basic state is one of mixed convection. The analysis to be carried out is an extension of the work by Barletta (Reference Barletta2012, Reference Barletta2013) and Sphaier & Barletta (Reference Sphaier and Barletta2014), with reference to the special cases of a horizontal and a vertical layer. The results obtained herein reveal an important feature peculiar to the inclined-channel situation: two different regimes occur, and one of them is unstable for any finite inclination angle, no matter how small. This is an important finding for experimental research involving mixed convection in horizontal channels, as it shows that the control of the channel inclination must be done in a very careful fashion. In fact, the usual scheme supporting a linear stability analysis is based, according to Lyapunov’s theory, on a test of the effects of small deviations from the initial conditions. However, vanishingly small deviations from the prescribed boundary conditions and from the geometrical layout can also be crucial when assessing the critical conditions for the onset of instability. This novel finding may cast new light on future investigations concerning the stability of convective flow systems.
2. Mathematical model
The problem considered in this study is that of a plane porous layer bounded by inclined impermeable walls separated by a distance
$H$
(see figure 1). A uniform wall heat flux,
$q_{0}$
, is prescribed both on the lower boundary plane
$z=0$
and on the upper boundary plane
$z=H$
. The buoyancy-driven flow in the porous layer is investigated under the following assumptions:
-
(i) the Oberbeck–Boussinesq approximation holds;
-
(ii) the saturated porous medium is considered as homogeneous, isotropic and with constant properties, except for the density in the gravitational body force;
-
(iii) the momentum transfer is described by Darcy’s law;
-
(iv) local thermal equilibrium between the solid phase and the fluid phase holds;
-
(v) viscous dissipation is negligible and no internal heat generation effects are present.
The dimensionless quantities are defined by rescaling the dimensional quantities as follows:
$$\begin{eqnarray}\displaystyle & \displaystyle (x,y,z)\,\frac{1}{H}\rightarrow (x,y,z),\quad t\,\frac{{\it\alpha}}{{\it\sigma}\,H^{2}}\rightarrow t, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle (u,v,w)\,\frac{H}{{\it\alpha}}\rightarrow (u,v,w),\quad (T-T_{0})\,\frac{k}{q_{0}\,H}\rightarrow T, & \displaystyle\end{eqnarray}$$
$(x,y,z)$
are the Cartesian coordinates,
$t$
is the time,
$\boldsymbol{u}=(u,v,w)$
is the velocity,
$T$
is the temperature and
$T_{0}$
is the constant reference temperature,
${\it\alpha}$
is the average thermal diffusivity of the saturated porous medium,
${\it\sigma}$
is the ratio between the volumetric heat capacity of the saturated medium and that of the fluid,
$g$
is the modulus of the gravitational acceleration
$\boldsymbol{g}$
,
${\it\beta}$
is the thermal expansion coefficient of the fluid,
$K$
is the permeability and
$k$
is the average thermal conductivity of the saturated porous medium.
With the above assumptions, the governing local balance equations can be written in dimensionless form as
$$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}\,\boldsymbol{\cdot }\,\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}\,\times \,\boldsymbol{u}=\mathit{Ra}\,\boldsymbol{{\rm\nabla}}\,\times \,[T\,(\sin {\it\phi}\,\hat{\boldsymbol{e}}_{x}+\cos {\it\phi}\,\hat{\boldsymbol{e}}_{z})], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial T}{\partial t}+\boldsymbol{u}\,\boldsymbol{\cdot }\,\boldsymbol{{\rm\nabla}}T={\rm\nabla}^{2}T, & \displaystyle\end{eqnarray}$$
where
$\hat{\boldsymbol{e}}_{x}$
,
$\hat{\boldsymbol{e}}_{z}$
are unit vectors along the
$x$
and
$z$
directions. The following boundary conditions are prescribed:
$$\begin{eqnarray}\displaystyle & \displaystyle z=0:\quad w=0,\quad -\frac{\partial T}{\partial z}=1, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle z=1:\quad w=0,\quad \frac{\partial T}{\partial z}=1. & \displaystyle\end{eqnarray}$$
As is shown in figure 1,
${\it\phi}$
is the channel inclination angle to the horizontal. The Darcy–Rayleigh number
$\mathit{Ra}$
is defined in terms of the heat flux:
$$\begin{eqnarray}\mathit{Ra}=\frac{g{\it\beta}q_{0}KH^{2}}{{\it\nu}{\it\alpha}k},\end{eqnarray}$$
where
${\it\nu}$
is the kinematic viscosity of the fluid. Either fluid cooling,
$\mathit{Ra}<0$
, or fluid heating,
$\mathit{Ra}>0$
, are allowed.
A time-independent solution of (2.2) is found by assuming a net horizontal flow rate along the
$x$
direction,
$$\begin{eqnarray}\displaystyle & \boldsymbol{u}_{b}=F(z)\,\hat{\boldsymbol{e}}_{x}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\rm\nabla}}T_{b}=\frac{2}{\mathit{Pe}}~\hat{\boldsymbol{e}}_{x}+\frac{1}{\mathit{Pe}}\,G(z)\,\hat{\boldsymbol{e}}_{z}, & \displaystyle\end{eqnarray}$$
where the functions
$F(z)$
and
$G(z)$
are given by
$$\begin{eqnarray}\displaystyle F(z) & = & \displaystyle \frac{{\it\Gamma}}{2}\left[\mathit{Pe}\,\coth \left(\frac{{\it\Gamma}}{2}\right)+2\tanh \left(\frac{{\it\Gamma}}{2}\right)\cot {\it\phi}\right]\cosh ({\it\Gamma}\,z)\nonumber\\ \displaystyle & & \displaystyle -\,\frac{{\it\Gamma}}{2}\,\left(\mathit{Pe}+2\,\cot {\it\phi}\right)\,\sinh ({\it\Gamma}\,z),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle G(z) & = & \displaystyle 2\cot {\it\phi}+\left[\mathit{Pe}\,\coth \left(\frac{{\it\Gamma}}{2}\right)+2\tanh \left(\frac{{\it\Gamma}}{2}\right)\cot {\it\phi}\right]\sinh ({\it\Gamma}\,z)\nonumber\\ \displaystyle & & \displaystyle -\,(\mathit{Pe}+2\cot {\it\phi})\,\cosh ({\it\Gamma}\,z),\end{eqnarray}$$
where the parameter
${\it\Gamma}$
is defined as
$$\begin{eqnarray}{\it\Gamma}=\sqrt{\frac{2\,\mathit{Ra }\,\sin {\it\phi}}{\mathit{Pe}}},\end{eqnarray}$$
and the Péclet number is the average dimensionless velocity in the porous channel:
$$\begin{eqnarray}\mathit{Pe}=\int _{0}^{1}F(z)\hspace{2.22198pt}\text{d}z.\end{eqnarray}$$
Two flow regimes exist: buoyancy-assisted flow
$(\mathit{Ra}/\mathit{Pe}>0)$
, and buoyancy-opposed flow
$(\mathit{Ra}/\mathit{Pe}<0)$
.
As noted from (2.4), positive values of
$\mathit{Pe}$
lead to an increasing temperature in the positive
$x$
direction, while a negative value naturally yields the opposite behaviour. Moreover, it is evident from these equations that the basic solution is left invariant by the transformation
$$\begin{eqnarray}\mathit{Ra}\rightarrow -\mathit{Ra},\quad \mathit{Pe}\rightarrow -\mathit{Pe},\quad {\it\phi}\rightarrow {\rm\pi}-{\it\phi},\quad x\rightarrow -x.\end{eqnarray}$$
Hence, it is not restrictive to focus the forthcoming stability analysis to cases where
${\it\phi}\in [0^{\circ },90^{\circ }]$
,
$\mathit{Ra}$
is positive, while
$\mathit{Pe}$
is either positive or negative. Further cases, not explicitly accounted for, can be easily recovered by applying the transformation (2.5).
3. Formulation of the stability problem
In order to study the linear stability of the basic flow, the basic solution is perturbed with very small disturbances:
$$\begin{eqnarray}\boldsymbol{u}=\boldsymbol{u}_{b}+{\it\epsilon}\,\boldsymbol{U},\quad T=T_{b}+{\it\epsilon}{\it\Theta},\end{eqnarray}$$
where
${\it\epsilon}$
is a real positive parameter, such that
${\it\epsilon}\ll 1$
.
Substituting (3.1) into (2.2), taking into account (2.4), and finally dropping terms
$O({\it\epsilon}^{2})$
yields
$$\begin{eqnarray}\displaystyle & \displaystyle {\rm\nabla}^{2}{\it\Psi}=\mathit{Ra}\left(\sin {\it\phi}\,\frac{\partial {\it\Theta}}{\partial x}+\cos {\it\phi}\,\frac{\partial {\it\Theta}}{\partial z}\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\Theta}}{\partial t}+F(z)\,\frac{\partial {\it\Theta}}{\partial x}+\frac{2}{\mathit{Pe}}\left(\!\mathit{Ra}\,{\it\Theta}\,\sin {\it\phi}-\frac{\partial {\it\Psi}}{\partial x}\!\right)+\frac{1}{\mathit{Pe}}~G(z)\left(\!\mathit{Ra}\,{\it\Theta}\,\cos {\it\phi}-\frac{\partial {\it\Psi}}{\partial z}\!\right)={\rm\nabla}^{2}{\it\Theta}, & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle z=0,1:\quad \frac{\partial {\it\Psi}}{\partial z}=\mathit{Ra}\,{\it\Theta}\,\cos {\it\phi},\quad \frac{\partial {\it\Theta}}{\partial z}=0, & \displaystyle\end{eqnarray}$$
${\it\Psi}$
was introduced so that the perturbation
$\boldsymbol{U}$
can be written as 
$$\begin{eqnarray}\boldsymbol{U}=\mathit{Ra}\left[{\it\Theta}\left(\sin {\it\phi}\,\hat{\boldsymbol{e}}_{x}+\cos {\it\phi}\,\hat{\boldsymbol{e}}_{z}\right)\right]-\boldsymbol{{\rm\nabla}}{\it\Psi}.\end{eqnarray}$$
3.1. Normal mode analysis
The auxiliary function
${\it\Psi}$
and the dimensionless temperature
${\it\Theta}$
are represented as plane waves propagating along an arbitrary direction in the
$(x,y)$
plane:
$$\begin{eqnarray}{\it\Psi}={\it\psi}(z)\,\text{e}^{\text{i}\,(a_{x}\,x+a_{y}\,y-{\it\omega}\,t)},\quad {\it\Theta}={\it\theta}(z)\,\text{e}^{\text{i}\,(a_{x}\,x+a_{y}\,y-{\it\omega}\,t)},\end{eqnarray}$$
where
$a_{x}$
and
$a_{y}$
are the real-valued components of the wavevector,
$a=(a_{x}^{2}+a_{y}^{2})^{1/2}$
is the wavenumber, and
${\it\omega}$
is the complex-valued frequency. Moreover, functions
${\it\psi}(z)$
and
${\it\theta}(z)$
are the amplitudes of the disturbances. The real part of
${\it\omega}$
defines the angular frequency, whereas the imaginary part determines the unstable (if positive) or stable (if negative) nature of the flow. Since this study focuses on the neutral stability analysis, one considers
$\text{Im}({\it\omega})=0$
, where
$\text{Im}$
denotes the imaginary part of a complex number.
Substitution of (3.4) into system (3.2) leads to the following eigenproblem:
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\psi}^{\prime \prime }(z)-a^{2}\,{\it\psi}(z)-\mathit{Ra}[\text{i}\,a_{x}\,{\it\theta}(z)\,\sin {\it\phi}+{\it\theta}^{\prime }(z)\,\cos {\it\phi}]=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\theta}^{\prime \prime }(z)-\left[a^{2}-\text{i}\,{\it\omega}+\text{i}\,a_{x}\,F(z)+\frac{2\,\mathit{Ra}}{\mathit{Pe}}\,\sin {\it\phi}+\frac{\mathit{Ra}}{\mathit{Pe}}\,G(z)\,\cos {\it\phi}\right]\,{\it\theta}(z) & \displaystyle \nonumber\\ \displaystyle & \displaystyle \hspace{-20.00003pt}+\,\text{i}\,a_{x}\,\frac{2}{\mathit{Pe}}\,{\it\psi}(z)+\frac{1}{\mathit{Pe}}\,G(z)\,{\it\psi}^{\prime }(z)=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\psi}^{\prime }(0)-\mathit{Ra}\,{\it\theta}(0)\,\cos {\it\phi}={\it\theta}^{\prime }(0)=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\psi}^{\prime }(1)-\mathit{Ra}\,{\it\theta}(1)\,\cos {\it\phi}={\it\theta}^{\prime }(1)=0, & \displaystyle\end{eqnarray}$$
$z$
. In order to have an eigenvalue problem with simpler boundary conditions, two new functions are introduced, 
$$\begin{eqnarray}f(z)={\it\psi}^{\prime }(z)-\mathit{Ra}\,{\it\theta}(z)\,\cos {\it\phi},\quad h(z)=a\,{\it\theta}(z),\end{eqnarray}$$
together with a new parameter
${\it\gamma}\in [0,1]$
, defined as
$$\begin{eqnarray}{\it\gamma}=\frac{a_{x}}{a}.\end{eqnarray}$$
The special cases
${\it\gamma}=0$
and
${\it\gamma}=1$
lead to longitudinal and transverse rolls, respectively, while
$0<{\it\gamma}<1$
defines the oblique rolls.
Using the new definitions (3.6), the eigenvalue problem can be rewritten as
$$\begin{eqnarray}\displaystyle & \displaystyle f^{\prime \prime }(z)-a^{2}\,f(z)-\mathit{Ra}\,[a\,h(z)\,\cos {\it\phi}+\text{i}\,{\it\gamma}\,h^{\prime }(z)\,\sin {\it\phi}]=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle h^{\prime \prime }(z)-\left[a^{2}-\text{i}\,{\it\omega}+\text{i}\,{\it\gamma}\,a\,F(z)+\frac{2\,\mathit{Ra}}{\mathit{Pe}}(1-{\it\gamma}^{2})\,\sin {\it\phi}\right]h(z) & \displaystyle \nonumber\\ \displaystyle & \displaystyle \hspace{30.00005pt}+\,\frac{2\,\text{i}\,{\it\gamma}}{\mathit{Pe}}\,f^{\prime }(z)+\frac{a}{\mathit{Pe}}\,G(z)\,f(z)=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle f(0)=h^{\prime }(0)=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle f(1)=h^{\prime }(1)=0. & \displaystyle\end{eqnarray}$$
$f(z)=h(z)=0$
, due to its homogeneous nature. However, the existing non-trivial solutions are the ones that are sought. Such solutions can be found only for specific values of
$\mathit{Ra}$
and
${\it\omega}$
, once (
$a,{\it\gamma},{\it\phi},\mathit{Pe}$
) are prescribed.
3.2. Splitting up into two coupled real eigenproblems
In order to analyse oblique and transverse rolls (i.e.
${\it\gamma}\neq 0$
), the eigenfunctions
$f$
and
$h$
are rewritten in terms of real and imaginary components:
$$\begin{eqnarray}f(z)=f_{\text{Re}}(z)+\text{i}\,f_{\text{Im}}(z),\quad h(z)=h_{\text{Re}}(z)+\text{i}\,h_{\text{Im}}(z).\end{eqnarray}$$
Thus, the original system (3.8) is rewritten as
$$\begin{eqnarray}\displaystyle & \displaystyle f_{\text{Re}}^{\prime \prime }(z)-a^{2}\,f_{\text{Re}}(z)-\mathit{Ra}\,[a\,h_{\text{Re}}(z)\,\cos {\it\phi}-{\it\gamma}\,h_{\text{Im}}^{\prime }(z)\,\sin {\it\phi}]=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle f_{\text{Im}}^{\prime \prime }(z)-a^{2}\,f_{\text{Im}}(z)-\mathit{Ra}\,[a\,h_{\text{Im}}(z)\,\cos {\it\phi}+{\it\gamma}\,h_{\text{Re}}^{\prime }(z)\,\sin {\it\phi}]=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle h_{\text{Re}}^{\prime \prime }(z)-\left[a^{2}+\frac{2\,\mathit{Ra}}{\mathit{Pe}}(1-{\it\gamma}^{2})\,\sin {\it\phi}\right]h_{\text{Re}}(z)-[{\it\omega}-{\it\gamma}\,a\,F(z)]h_{\text{Im}}(z) & \displaystyle \nonumber\\ \displaystyle \hspace{20.00003pt} & \displaystyle -\,\frac{2\,{\it\gamma}}{\mathit{Pe}}\,f_{\text{Im}}^{\prime }(z)+\frac{a}{\mathit{Pe}}\,G(z)\,f_{\text{Re}}(z)=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle h_{\text{Im}}^{\prime \prime }(z)-\left[a^{2}+\frac{2\,\mathit{Ra}}{\mathit{Pe}}(1-{\it\gamma}^{2})\,\sin {\it\phi}\right]h_{\text{Im}}(z)+[{\it\omega}-{\it\gamma}\,a\,F(z)]h_{\text{Re}}(z) & \displaystyle \nonumber\\ \displaystyle \hspace{20.00003pt} & \displaystyle +\,\frac{2\,{\it\gamma}}{\mathit{Pe}}\,f_{\text{Re}}^{\prime }(z)+\frac{a}{\mathit{Pe}}\,G(z)\,f_{\text{Im}}(z)=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle f_{\text{Re}}(0)=f_{\text{Im}}(0)=h_{\text{Re}}^{\prime }(0)=h_{\text{Im}}^{\prime }(0)=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle f_{\text{Re}}(1)=f_{\text{Im}}(1)=h_{\text{Re}}^{\prime }(1)=h_{\text{Im}}^{\prime }(1)=0. & \displaystyle\end{eqnarray}$$
$f$
and
$h$
, namely
$f_{\text{Re}}$
,
$f_{\text{Im}}$
,
$h_{\text{Re}}$
and
$h_{\text{Im}}$
, as well as for the eigenvalues
$\mathit{Ra}$
and
${\it\omega}$
.
4. Integral transform solution of eigenvalue problem
4.1. Integral transform pairs and eigenfunction bases
A solution to the eigenproblem (3.10a – f ) is now developed using the GITT (Cotta Reference Cotta1990, Reference Cotta1993, Reference Cotta1994). The solution process starts with the definition of the following transformation pairs:
$$\begin{eqnarray}\displaystyle & \displaystyle f~\text{inverse}~\Rightarrow ~f_{\text{Re}}(z)+\text{i}\,f_{\text{Im}}(z)=\mathop{\sum }_{n=1}^{\infty }\frac{\bar{f}_{n}^{\text{Re}}\,{\it\Omega}_{n}(z)}{N_{n}}+\text{i}\,\mathop{\sum }_{n=1}^{\infty }\frac{\bar{f}_{n}^{\text{Im}}\,{\it\Omega}_{n}(z)}{N_{n}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle f~\text{transform}~\Rightarrow ~\bar{f}_{n}^{\text{Re}}+\text{i}\,\bar{f}_{n}^{\text{Im}}=\int _{0}^{1}f_{\text{Re}}(z)\,{\it\Omega}_{n}(z)\hspace{2.22198pt}\text{d}z+\text{i}\,\int _{0}^{1}f_{\text{Im}}(z)\,{\it\Omega}_{n}(z)\hspace{2.22198pt}\text{d}z, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle h~\text{inverse}~\Rightarrow ~h_{\text{Re}}(z)+\text{i}\,h_{\text{Im}}(z)=\mathop{\sum }_{n=0}^{\infty }\frac{\bar{h}_{n}^{\text{Re}}\,{\it\Lambda}_{n}(z)}{\hat{N}_{n}}+\text{i}\,\mathop{\sum }_{n=0}^{\infty }\frac{\bar{h}_{n}^{\text{Im}}\,{\it\Lambda}_{n}(z)}{\hat{N}_{n}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle h~\text{transform}~\Rightarrow ~\bar{h}_{n}^{\text{Re}}+\text{i}\,\bar{h}_{n}^{\text{Im}}=\int _{0}^{1}h_{\text{Re}}(z)\,{\it\Lambda}_{n}(z)\hspace{2.22198pt}\text{d}z+\text{i}\,\int _{0}^{1}h_{\text{Im}}(z)\,{\it\Lambda}_{n}(z)\hspace{2.22198pt}\text{d}z, & \displaystyle\end{eqnarray}$$
${\it\Omega}_{n}$
and
${\it\Lambda}_{n}$
are infinite sets of orthogonal solutions of Sturm–Liouville-type eigenvalue problems, while
$N_{n}$
and
$\hat{N}_{n}$
are the corresponding
$\mathscr{L}^{2}$
-norms. For this particular application, Helmholtz problems that satisfy the same boundary conditions as the original eigensystem are selected: 
$$\begin{eqnarray}\displaystyle & {\it\Omega}_{n}^{\prime \prime }+{\it\mu}_{n}^{2}\,{\it\Omega}_{n}=0,\quad {\it\Omega}_{n}(0)={\it\Omega}_{n}(1)=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & {\it\Lambda}_{n}^{\prime \prime }+{\it\lambda}_{n}^{2}\,{\it\Lambda}_{n}=0,\quad {\it\Lambda}_{n}^{\prime }(0)={\it\Lambda}_{n}^{\prime }(1)=0. & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & {\it\Omega}_{n}=\sin ({\it\mu}_{n}\,z),\quad {\it\mu}_{n}=n\,{\rm\pi},~\text{for}~n=1,2,\dots , & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & {\it\Lambda}_{n}=\cos ({\it\lambda}_{n}\,z),\quad {\it\lambda}_{n}=n\,{\rm\pi},~\text{for}~n=0,1,2,\dots , & \displaystyle\end{eqnarray}$$
where the associated norms are given by
$$\begin{eqnarray}\displaystyle & \displaystyle N_{n}=\int _{0}^{1}{\it\Omega}_{n}^{2}\hspace{2.22198pt}\text{d}z=\frac{1}{2}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \hat{N}_{n}=\int _{0}^{1}{\it\Lambda}_{n}^{2}\hspace{2.22198pt}\text{d}z=\left\{\begin{array}{@{}ll@{}}1\quad & \text{for}~n=0,\\ {\textstyle \frac{1}{2}}\quad & \text{for}~n>0.\end{array}\right. & \displaystyle\end{eqnarray}$$
4.2. Integral transformation
The transformation of the problem is carried out by multiplying (3.10a
) and (3.10b
) by
${\it\Omega}_{n}$
and (3.10c
) and (3.10d
) by
${\it\Lambda}_{n}$
and integrating the resulting equations within the problem domain. By applying the boundary conditions, substituting the inversion formulas into the non-transformable terms, and simplifying, one obtains
$$\begin{eqnarray}\displaystyle & & \displaystyle -{\it\mu}_{n}^{2}\,\bar{f}_{n}^{\text{Re}}-a^{2}\,\bar{f}_{n}^{\text{Re}}-a\,\mathit{Ra}\,\cos {\it\phi}\mathop{\sum }_{m=0}^{\infty }A_{n,m}\bar{h}_{m}^{\text{Re}}+\mathit{Ra}\,{\it\gamma}\,\sin {\it\phi}\mathop{\sum }_{m=0}^{\infty }R_{n,m}\bar{h}_{m}^{\text{Im}}=0,\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle -{\it\mu}_{n}^{2}\,\bar{f}_{n}^{\text{Im}}-a^{2}\,\bar{f}_{n}^{\text{Im}}-a\,\mathit{Ra}\,\cos {\it\phi}\mathop{\sum }_{m=0}^{\infty }A_{n,m}\bar{h}_{m}^{\text{Im}}-\mathit{Ra}\,{\it\gamma}\,\sin {\it\phi}\mathop{\sum }_{m=0}^{\infty }R_{n,m}\bar{h}_{m}^{\text{Re}}=0,\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle -{\it\lambda}_{n}^{2}\,\bar{h}_{n}^{\text{Re}}\,-\left[a^{2}\,+\frac{2\,\mathit{Ra}}{\mathit{Pe}}(1-{\it\gamma}^{2})\,\sin {\it\phi}\right]\bar{h}_{n}^{\text{Re}}-{\it\omega}\,\bar{h}_{n}^{\text{Im}}+{\it\gamma}\,a\,\mathop{\sum }_{m=0}^{\infty }B_{n,m}\,\bar{h}_{m}^{\text{Im}}\nonumber\\ \displaystyle & & \displaystyle \hspace{60.00009pt}-\,2\,\frac{{\it\gamma}}{\mathit{Pe}}\,\mathop{\sum }_{m=1}^{\infty }S_{n,m}\,\bar{f}_{m}^{\text{Im}}+\frac{a}{\mathit{Pe}}\mathop{\sum }_{m=1}^{\infty }C_{n,m}\,\bar{f}_{m}^{\text{Re}}=0,\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle -{\it\lambda}_{n}^{2}\,\bar{h}_{n}^{\text{Im}}\,-\left[a^{2}\,+\frac{2\,\mathit{Ra}}{\mathit{Pe}}(1-{\it\gamma}^{2})\,\sin {\it\phi}\right]\bar{h}_{n}^{\text{Im}}+{\it\omega}\,\bar{h}_{n}^{\text{Re}}-{\it\gamma}\,a\,\mathop{\sum }_{m=0}^{\infty }B_{n,m}\,\bar{h}_{m}^{\text{Re}}\nonumber\\ \displaystyle & & \displaystyle \hspace{60.00009pt}+\,2\,\frac{{\it\gamma}}{\mathit{Pe}}\,\mathop{\sum }_{m=1}^{\infty }S_{n,m}\,\bar{f}_{m}^{\text{Re}}+\frac{a}{\mathit{Pe}}\mathop{\sum }_{m=1}^{\infty }C_{n,m}\,\bar{f}_{m}^{\text{Im}}=0,\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle A_{n,m}=\frac{1}{\hat{N}_{m}}\int _{0}^{1}{\it\Lambda}_{m}\,{\it\Omega}_{n}\hspace{2.22198pt}\text{d}z, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle B_{n,m}=\frac{1}{\hat{N}_{m}}\int _{0}^{1}F(z)\,{\it\Lambda}_{m}\,{\it\Lambda}_{n}\hspace{2.22198pt}\text{d}z, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle C_{n,m}=\frac{1}{N_{m}}\int _{0}^{1}G(z)\,{\it\Omega}_{m}\,{\it\Lambda}_{n}\hspace{2.22198pt}\text{d}z, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle S_{n,m}=\frac{1}{N_{m}}\int _{0}^{1}{\it\Omega}_{m}^{\prime }\,{\it\Lambda}_{n}\hspace{2.22198pt}\text{d}z, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle R_{n,m}=\frac{1}{\hat{N}_{m}}\int _{0}^{1}{\it\Lambda}_{m}^{\prime }\,{\it\Omega}_{n}\hspace{2.22198pt}\text{d}z. & \displaystyle\end{eqnarray}$$
$n=1,2,\dots ,\infty$
, while (4.4c
,
d
) are valid for
$n=0,1,\dots ,\infty$
.
System (4.4) is then truncated to a finite order
$n_{max}$
, and rearranged in a compact vector form, which is obtained by eliminating the vectors
$\bar{\boldsymbol{f}}_{\text{Re}}$
and
$\bar{\boldsymbol{f}}_{\text{Im}}$
from the truncated algebraic system:
$$\begin{eqnarray}\displaystyle & \displaystyle -\unicode[STIX]{x1D640}\,\bar{\boldsymbol{h}}_{\text{Re}}-\unicode[STIX]{x1D646}\,\bar{\boldsymbol{h}}_{\text{Im}}+\frac{\mathit{Ra}}{\mathit{Pe }}\left(\unicode[STIX]{x1D642}\,\bar{\boldsymbol{h}}_{\text{Re}}+\unicode[STIX]{x1D641}\,\bar{\boldsymbol{h}}_{\text{Im}}\right)=\mathbf{0}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle -\unicode[STIX]{x1D640}\,\bar{\boldsymbol{h}}_{\text{Im}}+\unicode[STIX]{x1D646}\,\bar{\boldsymbol{h}}_{\text{Re}}+\frac{\mathit{Ra}}{\mathit{Pe }}\left(\unicode[STIX]{x1D642}\,\bar{\boldsymbol{h}}_{\text{Im}}-\unicode[STIX]{x1D641}\,\bar{\boldsymbol{h}}_{\text{Re}}\right)=\mathbf{0}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D641}=a\,{\it\gamma}\,\left[2\,\cos {\it\phi}\,\unicode[STIX]{x1D64E}\,\unicode[STIX]{x1D63F}\,\unicode[STIX]{x1D63C}+\sin {\it\phi}\,\unicode[STIX]{x1D63E}\,\unicode[STIX]{x1D63F}^{-1}\unicode[STIX]{x1D64D}\right], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D642}=2\,\sin {\it\phi}\left[{\it\gamma}^{2}\,\unicode[STIX]{x1D64E}\,\unicode[STIX]{x1D63F}\,\unicode[STIX]{x1D64D}-\,(1-{\it\gamma}^{2})\,\,\unicode[STIX]{x1D644}\right]-a^{2}\,\cos {\it\phi}\,\unicode[STIX]{x1D63E}\,\unicode[STIX]{x1D63F}\,\unicode[STIX]{x1D63C}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D646}=({\it\omega}\,\unicode[STIX]{x1D644}-{\it\gamma}\,a\,\unicode[STIX]{x1D63D}) & \displaystyle\end{eqnarray}$$
and
$\unicode[STIX]{x1D63F}$
and
$\unicode[STIX]{x1D640}$
are diagonal matrices given by the following coefficients:
$$\begin{eqnarray}\unicode[STIX]{x1D60B}_{n,n}=({\it\mu}_{n}^{2}+a^{2})^{-1},\quad \unicode[STIX]{x1D60C}_{n,n}=({\it\lambda}_{n}^{2}+a^{2}).\end{eqnarray}$$
System (4.6) is a nonlinear eigenvalue problem if
$\mathit{Ra}$
is to be calculated in terms of the remaining parameters, while it is a linear eigenvalue problem if the angular frequency
${\it\omega}$
is to be calculated in terms of the remaining parameters. Since the computation of
${\it\omega}$
in terms of the other parameters requires considerably less CPU time, system (4.6) is rewritten in the form
$$\begin{eqnarray}\unicode[STIX]{x1D643}\,\boldsymbol{y}={\it\omega}\,\unicode[STIX]{x1D648}\,\boldsymbol{y},\end{eqnarray}$$
in which
$\unicode[STIX]{x1D643}$
and
$\unicode[STIX]{x1D648}$
are block matrices:
The new vector
$\boldsymbol{y}$
is defined as
$$\begin{eqnarray}\boldsymbol{y}=(\bar{h}_{0}^{\text{Re}},\bar{h}_{1}^{\text{Re}},\dots ,\,\bar{h}_{n_{max}}^{\text{Re}},\bar{h}_{0}^{\text{Im}},\bar{h}_{1}^{\text{Im}},\dots ,\,\bar{h}_{n_{max}}^{\text{Im}}).\end{eqnarray}$$
System (4.8) constitutes a generalised linear algebraic eigenvalue problem, which is finally solved numerically using a single run of Mathematica (© Wolfram Research) function Eigenvalues.
4.3. Simplification for large Péclet numbers
The limiting behaviour of
$F(z)$
and
$G(z)$
for
$\mathit{Pe}\gg 1$
is such that
$$\begin{eqnarray}\lim _{\mathit{Pe}\rightarrow \infty }{\displaystyle \frac{F(z)}{\mathit{Pe}}}=1,\quad \lim _{\mathit{Pe}\rightarrow \infty }{\displaystyle \frac{G(z)}{\mathit{Pe}}}=2\,z-1.\end{eqnarray}$$
Thus, (3.8b ) is simplified to yield
$$\begin{eqnarray}h^{\prime \prime }(z)-(a^{2}-\text{i}\,{\it\chi})h(z)+a\,(2\,z-1)\,f(z)=0,\end{eqnarray}$$
where
${\it\chi}={\it\omega}-{\it\gamma}\,a\,\mathit{Pe}$
is assumed to be
$O(1)$
when
$Pe\rightarrow \infty$
.
With this simplification, the truncated transformed system (4.6) is reduced to
$$\begin{eqnarray}\displaystyle & -\unicode[STIX]{x1D640}\,\bar{\boldsymbol{h}}_{\text{Re}}-{\it\chi}\,\bar{\boldsymbol{h}}_{\text{Im}}+\mathit{Ra}\left(\unicode[STIX]{x1D642}^{+}\,\bar{\boldsymbol{h}}_{\text{Re}}+\unicode[STIX]{x1D641}^{+}\,\bar{\boldsymbol{h}}_{\text{Im}}\right)=\mathbf{0}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & -\unicode[STIX]{x1D640}\,\bar{\boldsymbol{h}}_{\text{Im}}+{\it\chi}\,\bar{\boldsymbol{h}}_{\text{Re}}+\mathit{Ra}\left(\unicode[STIX]{x1D642}^{+}\,\bar{\boldsymbol{h}}_{\text{Im}}-\unicode[STIX]{x1D641}^{+}\,\bar{\boldsymbol{h}}_{\text{Re}}\right)=\mathbf{0}, & \displaystyle\end{eqnarray}$$
where
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D641}^{+}=a\,{\it\gamma}\,\sin {\it\phi}\,\unicode[STIX]{x1D63E}^{+}\,\unicode[STIX]{x1D63F}\,\unicode[STIX]{x1D64D}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D642}^{+}=-a^{2}\,\cos {\it\phi}\,\unicode[STIX]{x1D63E}^{+}\,\unicode[STIX]{x1D63F}\,\unicode[STIX]{x1D63C}, & \displaystyle\end{eqnarray}$$
and
$\unicode[STIX]{x1D63E}^{+}$
is a matrix generated by the coefficients
$$\begin{eqnarray}\unicode[STIX]{x1D60A}_{n,m}^{+}=\frac{1}{N_{m}}\int _{0}^{1}(2\,z-1)\,{\it\Omega}_{m}\,{\it\Lambda}_{n}\hspace{2.22198pt}\text{d}z.\end{eqnarray}$$
Introducing the
$\boldsymbol{y}$
vector, as defined by (4.9), equations (4.12) can be written in the form
$$\begin{eqnarray}\unicode[STIX]{x1D643}^{+}\,\boldsymbol{y}={\it\chi}\,\unicode[STIX]{x1D648}\,\boldsymbol{y},\end{eqnarray}$$
in which
$\unicode[STIX]{x1D643}^{+}$
is the following block matrix:
Finally, the solution of system (4.14) is directly obtained by using Mathematica function Eigenvalues, as similarly employed for solving system (4.8).
4.4. Calculation of
$\mathit{Ra}$
and
$a$
in terms of
${\it\omega}$
or
${\it\chi}$
The previous sections described how to obtain the angular frequency
${\it\omega}$
, or its modified counterpart
${\it\chi}$
, using a single numerical stage which involves the calculation of matrix eigenvalues from given values of
$a$
and
$\mathit{Ra}$
, as well as
$\mathit{Pe}$
,
${\it\phi}$
and
${\it\gamma}$
. Nevertheless, in order to draw neutral stability curves and to determine critical Rayleigh number values, an additional step must be incorporated. This involves calculating
$a$
and
$\mathit{Ra}$
(and
${\it\omega}$
or
${\it\chi}$
as well) from a given configuration, that is, from given
$\mathit{Pe}$
,
${\it\phi}$
and
${\it\gamma}$
values. This step is carried out by enforcing a condition that the angular frequency must be real, thus satisfying the neutral stability requirement. This is done by calculating values of
$\mathit{Ra}$
and
$a$
that yield
$\text{Im}({\it\omega})=0$
. In the Mathematica implementation employed, this is accomplished by using the FindRoot routine (Newton and secant-type methods) to find the values of
$\mathit{Ra}$
and
$a$
that satisfy
$\text{Im}({\it\omega})=0$
.
5. Asymptotic solution for
$a\rightarrow 0$
Let us expand eigenfunctions
$f$
and
$h$
, as well as eigenvalues
$\mathit{Ra}$
and
${\it\omega}$
, in powers of
$a$
:
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}f(z)=f_{0}(z)+f_{1}(z)\,a+f_{2}(z)\,a^{2}+O(a^{3}),\\ h(z)=h_{0}(z)+h_{1}(z)\,a+h_{2}(z)\,a^{2}+O(a^{3}),\\ \mathit{Ra}=\mathit{Ra}_{0}+\mathit{Ra}_{1}\,a+\mathit{Ra}_{2}\,a^{2}+O(a^{3}),\\ {\it\omega}={\it\omega}_{0}+{\it\omega}_{1}\,a+{\it\omega}_{2}\,a^{2}+O(a^{3}).\end{array}\right\}\end{eqnarray}$$
By keeping only terms of order 0 in
$a$
, the eigenproblem (3.8) for arbitrary oblique rolls yields
$$\begin{eqnarray}\displaystyle & f_{0}^{\prime \prime }(z)-\text{i}\,\mathit{Ra}_{0}\,{\it\gamma}\,h_{0}^{\prime }(z)\,\sin {\it\phi}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle h_{0}^{\prime \prime }(z)+\left[\text{i}\,{\it\omega}_{0}-\frac{2\,\mathit{Ra}_{0}}{\mathit{Pe}}(1-{\it\gamma}^{2})\sin {\it\phi}\right]h_{0}(z)+\frac{2\,\text{i}\,{\it\gamma}}{\mathit{Pe}}\,f_{0}^{\prime }(z)=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & f_{0}(0)=h_{0}^{\prime }(0)=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & f_{0}(1)=h_{0}^{\prime }(1)=0. & \displaystyle\end{eqnarray}$$
$z$
, so that we get 
$$\begin{eqnarray}\displaystyle & \displaystyle h_{0}^{\prime \prime \prime }+\left(\text{i}\,{\it\omega}_{0}-\frac{2\,\mathit{Ra}_{0}}{\mathit{Pe}}\,\sin {\it\phi}\right)h_{0}^{\prime }=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & h_{0}^{\prime }(0)=h_{0}^{\prime }(1)=0. & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}h_{0}=\cos (n\,{\rm\pi}\,z)+C,\end{eqnarray}$$
where
$n=0,1,2,\dots$
and
$C$
is a real constant, provided that the dispersion relation,
$$\begin{eqnarray}n^{2}{\rm\pi}^{2}=\text{i}\,{\it\omega}_{0}-\frac{2\,\mathit{Ra}_{0}}{\mathit{Pe}}\,\sin {\it\phi},\end{eqnarray}$$
holds. Equation (5.5) allows one to conclude that
${\it\omega}_{0}=0$
and that
$$\begin{eqnarray}\mathit{Ra}_{0}\cos {\it\phi}=-\frac{\mathit{Pe}}{2}\,n^{2}\,{\rm\pi}^{2}\cot {\it\phi}.\end{eqnarray}$$
On account of (5.6), for any buoyancy-opposed basic flow (
$\mathit{Ra}/\mathit{Pe}<0$
) and for any non-vanishing inclination angle
${\it\phi}$
, one may conclude that the neutral stability curves in the
$(a,\mathit{Ra}\,\cos {\it\phi})$
-plane display an infinite sequence of points at
$a=0$
. The lowest is
$\mathit{Ra}\cos {\it\phi}=0$
, and this proves that even a vanishingly small inclination to the horizontal, producing a buoyancy-opposed flow, is capable of yielding instability. Higher values are
$$\begin{eqnarray}\mathit{Ra}\,\cos {\it\phi}=-\frac{\mathit{Pe}\,{\rm\pi}^{2}}{2}\,\cot {\it\phi},-2\,\mathit{Pe}\,{\rm\pi}^{2}\cot {\it\phi},\dots .\end{eqnarray}$$
This information is used in the numerical calculation of the higher neutral stability curves that occur for buoyancy-opposed flow.
6. Results and discussion
In this section, numerical results are presented and discussed. The first data are displayed for the validation of the hybrid methodology employed to obtain the solution of the differential eigenvalue problem (3.8). The results obtained with this methodology are compared with the data computed using a numerical shooting scheme based on a fourth-order Runge–Kutta (RK4) solver. More details on this numerical technique can be found in Barletta (Reference Barletta2012). Tables 1 and 2, relating to
$\mathit{Pe}\rightarrow \infty$
and
$\mathit{Pe}=50$
respectively, display the critical values of the rescaled Rayleigh number
$\mathit{Ra}\,\cos {\it\phi}$
, wavenumber
$a$
and modified angular frequency
${\it\chi}$
calculated with various truncation orders. The results are presented for longitudinal rolls (
${\it\gamma}=0$
) and transverse rolls (
${\it\gamma}=1$
), with reference to inclination angles
${\it\phi}=1^{\circ },10^{\circ }$
.
Table 1. Critical Rayleigh number, wavenumber and modified angular frequency for different truncation orders, for
$\mathit{Pe}\rightarrow \infty$
, and different values of
${\it\phi}$
and
${\it\gamma}$
.
Table 2. Critical Rayleigh number, wavenumber and modified angular frequency for different truncation orders, for
$\mathit{Pe}=50$
, and different values of
${\it\phi}$
and
${\it\gamma}$
.
As can be seen from the results, there is perfect agreement between the integral transform implementation and the solution obtained using the fully numerical RK4-shooting scheme. All values agree within six significant digits for a truncation order of 80. However, if four-digit precision is desired, working with
$n_{max}=20$
is sufficient, while 10 terms give an error of less than 1 %.
After examining the convergence behaviour of the solution, the next results are focused on the stability analysis itself. As mentioned earlier, the effect of the channel inclination to the horizontal is different depending on the sign of the Péclet number. On assuming positive values of
$\mathit{Ra}$
, a positive
$\mathit{Pe}$
yields a buoyancy-assisted flow, whereas a negative
$\mathit{Pe}$
yields a buoyancy-opposed flow. Figure 2 shows neutral stability curves for longitudinal rolls in the plane
$(a,\mathit{Ra}\cos {\it\phi})$
with different inclination angles and different Péclet numbers. This figure displays neutral stability curves for the buoyancy-assisted regime (
$\mathit{Pe}>0$
) including the limiting case with infinite
$\mathit{Pe}$
. As one can observe from these curves, when the Péclet number is very large, the neutral stability curves reach an asymptotic regime where the onset of convection is practically independent of the inclination angle. We note that the case with
$\mathit{Pe}=500$
almost matches that with
$\mathit{Pe}\rightarrow \infty$
. Conversely, for smaller values of
$\mathit{Pe}$
, the neutral stability curves move upward as
${\it\phi}$
increases, thus describing more stable conditions. When
$\mathit{Pe}$
becomes as small as 20, the neutral stability curves assume the shape of closed loops that shrink in size and, as
$\mathit{Pe}$
further decreases, they collapse and eventually disappear. This phenomenon is first displayed with small inclinations and then it also appears for larger values of
${\it\phi}$
. As already observed, the horizontal case, i.e.
${\it\phi}=0^{\circ }$
, is always the most unstable case. This implies, for the buoyancy-assisted regime, that horizontal channels are more prone to the development of instability than their inclined counterparts.
Figure 2. Buoyancy-assisted flow: neutral stability curves in the plane
$(a,\mathit{Ra}\cos {\it\phi})$
with different Péclet numbers and inclination angles for longitudinal rolls (
${\it\gamma}=0$
): (a)
$Pe=\infty$
, (b) 500, (c) 100, (d) 50, (e) 25, (f) 20.
Figure 3 displays neutral stability curves for longitudinal rolls in the plane
$(a,\mathit{Ra}\cos {\it\phi})$
with different inclination angles and different Péclet numbers, for the buoyancy-opposed flow regime. With the exception of the horizontal channel, whose basic state is independent of the flow direction, very different behaviour is seen in comparison with the buoyancy-assisted flow. The main feature of figure 3 is that the neutral stability curves display their minimum at
$\mathit{Ra}=0$
and
$a=0$
, as predicted by the asymptotic analysis for small wavenumbers (see § 5). Moreover, as predicted by the same analysis, the curves intersect the vertical axis at multiple positions. Higher branches of neutral stability are also displayed in figure 3.
Figure 3. Buoyancy-opposed flow: neutral stability curves in the plane
$(a,\mathit{Ra}\cos {\it\phi})$
with different Péclet numbers and inclination angles for longitudinal rolls (
${\it\gamma}=0$
): (a)
$Pe=-500$
, (b)
$-100$
, (c)
$-50$
, (d)
$-25$
, (e)
$-10$
, (f)
$-1$
.
Figure 4 refers to oblique rolls and displays neutral stability curves in the plane
$(a,\mathit{Ra}\cos {\it\phi})$
for the case of buoyancy-assisted flow with
${\it\phi}=5^{\circ }$
and different values of
${\it\gamma}$
and
$\mathit{Pe}$
. This figure shows that the longitudinal rolls are the most unstable, as the lowest curves are always those with
${\it\gamma}=0$
. It is also interesting to note that, for small inclination angles, the dependence on
${\it\gamma}$
is minor, which is further reduced for higher values of
$\mathit{Pe}$
. In fact, in the limiting situation of an infinite Péclet number, there is neither dependence on
${\it\gamma}$
nor on the inclination angle. Results for a null inclination angle are not presented, since the dependence on
${\it\gamma}$
is even smaller than for
${\it\phi}=5^{\circ }$
, being noticeable only for smaller values of
$\mathit{Pe}$
, as previously documented in Barletta (Reference Barletta2012). Following the previous results, higher inclinations are now analysed. Figures 5 and 6 display the neutral stability curves relative to oblique rolls for
${\it\phi}=15^{\circ }$
and
${\it\phi}=30^{\circ }$
, respectively. As one can infer from these figures, the longitudinal rolls are still the most unstable. In contrast to what was discussed for small inclination angles, the neutral stability curves are much more dependent on the oblique roll parameter
${\it\gamma}$
when larger inclination angles are considered. If the inclination is sufficiently high (e.g. 30°), the curves for larger values of
${\it\gamma}$
move upwards to such an extent that they are no longer visible within the vertical range of the plots. As well as this phenomenon, for sufficiently small Péclet numbers, some neutral stability curves become a closed loop, which shrinks and collapses to a point as the value of
$\mathit{Pe}$
is further reduced, similarly to what happens for longitudinal rolls when
$\mathit{Pe}$
approaches its lowest bound for instability. This means that unstable oblique rolls with larger values of
${\it\gamma}$
disappear for a sufficiently high inclination. This occurrence is further intensified as
$\mathit{Pe}$
decreases so that, for sufficiently small Péclet numbers and a sufficiently high channel inclination, there are no unstable transverse rolls. In order to further clarify what happens with the neutral stability curves for
$\mathit{Pe}=500$
and
${\it\phi}=30^{\circ }$
, as
${\it\gamma}$
is increased past
$5/8$
, additional curves for this case are displayed in figure 7. As can be seen from this figure, the point on the curve for
${\it\gamma}=5/8$
near
$a=6$
and
$\mathit{Ra}\,\cos ({\it\phi})=600$
rapidly moves upwards, splitting the curves in two parts, as
${\it\gamma}$
is increased from
$5/8$
to
$3/4$
. The split occurs just after
${\it\gamma}=21/32$
; for
${\it\gamma}=169/256$
the curve is already separated into two branches. As
${\it\gamma}$
is further increased the left branch of neutral stability shrinks and moves upwards at a higher rate than the branch on the right, until it ceases to exist, as was observed in figure 6 for
${\it\gamma}=7/8$
.
Figure 4. Buoyancy-assisted flow: neutral stability curves in the plane
$(a,\mathit{Ra}\cos {\it\phi})$
with different Péclet numbers and values of
${\it\gamma}$
for
${\it\phi}=5^{\circ }$
: (a)
$Pe=500$
, (b) 100, (c) 50, (d) 25.
Figure 5. Buoyancy-assisted flow: neutral stability curves in the plane
$(a,\mathit{Ra}\cos {\it\phi})$
with different Péclet numbers and values of
${\it\gamma}$
for
${\it\phi}=15^{\circ }$
: (a)
$Pe=500$
, (b) 100, (c) 50, (d) 25.
Figure 6. Buoyancy-assisted flow: neutral stability curves in the plane
$(a,\mathit{Ra}\cos {\it\phi})$
with different Péclet numbers and values of
${\it\gamma}$
for
${\it\phi}=30^{\circ }$
: (a)
$Pe=500$
, (b) 100, (c) 50, (d) 25.
Figure 7. Buoyancy-opposed flow: neutral stability curves in the plane
$(a,\mathit{Ra}\cos {\it\phi})$
with
$\mathit{Pe}=500$
and
${\it\phi}=30^{\circ }$
and different values of
${\it\gamma}$
.
Figure 8 displays neutral stability curves for oblique rolls, ranging from longitudinal (
${\it\gamma}=0$
) to transverse (
${\it\gamma}=1$
), under different channel inclination angles for cases of buoyancy-opposed flow. This figure shows data for a channel inclination angle
${\it\phi}=5^{\circ }$
. As one can infer from the plots, for most of the wavenumber range presented there is a small variation in the curves with
${\it\gamma}$
, especially for larger values of the Péclet number, which was similarly observed for buoyancy-assisted flows. However, as smaller values of
$a$
are involved (i.e. near zero), a notable difference between the curves is seen, where curves with higher
${\it\gamma}$
values move above the local maximum that occurs in the range
$0<a<1$
. This behaviour is observed for the cases with
$\mathit{Pe}=-500$
,
$-100$
and
$-50$
. In the case
$\mathit{Pe}=-25$
, there is no local maximum as the neutral stability curves for all values of
${\it\gamma}$
are separated into two disconnected branches. As already pointed out in relation to the previous figures, figure 8 confirms that longitudinal rolls are more unstable than oblique and transverse rolls.
Figure 8. Buoyancy-opposed flow: neutral stability curves in the plane
$(a,\mathit{Ra}\cos {\it\phi})$
with different Péclet numbers and values of
${\it\gamma}$
for
${\it\phi}=5^{\circ }$
: (a)
$Pe=-500$
, (b)
$-100$
, (c)
$-50$
, (d)
$-25$
.
Figures 9 and 10 show similar results for 15° and 30°, respectively. In these higher channel inclination angles, there is a larger difference between curves with different values of
${\it\gamma}$
throughout the wavenumber spectrum, again with the minimum for a given value of
$a$
occurring for
${\it\gamma}=0$
. For cases where
${\it\gamma}$
is sufficiently large, the curves again separate into two regions, where the one to the right assumes the shape of an inverted C or becomes a closed loop. Although this can easily be seen for
$\mathit{Pe}=-25$
and
${\it\phi}=15^{\circ }$
, it occurs also in other instances. Most of them are not displayed in these figures as it would require much larger vertical and horizontal scales.
Figure 9. Buoyancy-opposed flow: neutral stability curves in the plane
$(a,\mathit{Ra}\cos {\it\phi})$
with different Péclet numbers and values of
${\it\gamma}$
for
${\it\phi}=15^{\circ }$
: (a)
$Pe=-500$
, (b)
$-100$
, (c)
$-50$
, (d)
$-25$
.
Figure 10. Buoyancy-assisted flow: neutral stability curves in the plane
$(a,\mathit{Ra}\cos {\it\phi})$
with different Péclet numbers and values of
${\it\gamma}$
for
${\it\phi}=30^{\circ }$
: (a)
$Pe=-500$
, (b)
$-100$
, (c)
$-50$
, (d)
$-25$
.
By performing an overall view of figures 8–10, a few additional points can be highlighted. The portion of the neutral stability curves that connects the two intercepts with the vertical axis is a perfect straight line for transverse rolls, which can be seen in all of the previous figures. Regardless of this, transverse rolls can also lead to a closed region on the right. With respect to the closed detached region on the right, one can also see that curves for a given value of
${\it\gamma}$
are always enclosed by curves with smaller
${\it\gamma}$
. For some situations, as the Péclet number decreases or the inclination increases, the closed portion of the curve on the right collapses and eventually disappears.
In figure 11, the behaviour of the critical Rayleigh number versus the Péclet number is displayed. Obviously, only the results for buoyancy-assisted flow are displayed, as the buoyancy-opposed flow regime leads to
$\mathit{Ra}_{c}=0$
for all cases with a finite Péclet number. The results confirm the previous observations about the existence of a minimum Péclet number value that may lead to instability. This effect was demonstrated in Barletta (Reference Barletta2012), for a horizontal channel, and it is corroborated here for
${\it\phi}=0^{\circ }$
. In fact, the case with
${\it\phi}=0^{\circ }$
yields the lowest
$\mathit{Pe}$
compatible with the onset of convective instability. As the channel inclination increases, the lowest unstable
$\mathit{Pe}$
becomes higher, reflecting the fact that increasing the inclination of the channel has a stabilising effect on the flow.
Figure 11. Critical Rayleigh number,
$\mathit{Ra}_{c}$
, versus
$\mathit{Pe}$
with different inclination angles for: (a) longitudinal (
${\it\gamma}=0$
), (b) oblique (
${\it\gamma}=1/3$
), (c) oblique (
$2/3$
) and (d) transverse (
${\it\gamma}=1$
) rolls in buoyancy-assisted flow.
7. Summary and conclusions
A linear stability analysis of the mixed convection in a porous channel inclined to the horizontal was carried out. Both channel walls were assumed to be impermeable and uniformly heated or cooled by a uniform flux. The stability of the basic stationary and parallel flow was tested in a traditional way, by perturbing the basic solution with small-amplitude plane waves propagating along an arbitrary direction parallel to the bounding walls. This formulation led to a complex differential eigenvalue problem where the perturbation amplitudes were the desired eigenfunctions, and where the Rayleigh number,
$\mathit{Ra}$
, the wavenumber,
$a$
, and the angular frequency,
${\it\omega}$
, were the eigenvalues. The computation was performed for different Péclet numbers,
$\mathit{Pe}$
, channel inclinations to the horizontal,
${\it\phi}$
, and disturbance wavevector orientations, defined by the parameter
${\it\gamma}$
. A symmetry of the basic solution was determined, allowing the analysis to be focused on cases where
$\mathit{Ra}>0$
and
${\it\phi}\in [0^{\circ },90^{\circ }]$
, with a negative or positive
$\mathit{Pe}$
.
The stability eigenvalue problem was solved by means of the GITT framework, namely a hybrid analytical–numerical method. By employing the GITT, the differential eigenproblem was transformed into an algebraic eigenvalue problem, which was solved numerically using a matrix eigenvalue calculation routine. Since this problem is nonlinear in
$\mathit{Ra}$
and
$a$
, but linear in the angular frequency
${\it\omega}$
, this computational algorithm was applied such that
${\it\omega}$
could be calculated from any combination of input data
$(a,\mathit{Ra},\mathit{Pe},{\it\gamma},{\it\phi})$
. Then, by constraining the imaginary part of
${\it\omega}$
to be zero, the neutral stability curves in the plane
$(a,\mathit{Ra}\cos {\it\phi})$
were determined. Compared to traditional numerical methods (e.g. Runge–Kutta) for the calculation of the eigenvalues involved in the stability analysis, the methodology presented herein has the advantage of computing multiple eigenvalues at the same time. Furthermore, the semi-analytical nature of the integral transform methodology allows very efficient control of the global solution error. In addition, it is important to emphasise that the solution algorithm can be quickly modified to solve other eigenvalue problems that arise in stability analysis. The main changes would be different expressions for the integral coefficients and a different form for the block matrices that appear in the generalised algebraic eigenvalue problem (equations (4.14)).
The methodology was verified by performing a convergence analysis of the solution, and a comparison with a fully numerical solution based on the Runge–Kutta method and the shooting method. We have shown that a 1 % error can be ensured with 10 terms in the truncated series and four-digit precision can be obtained if the truncation is raised to 20 terms. All figures were plotted using 20 terms in the series to ensure graphical convergence for the presented results. As well as this verification, an asymptotic analysis for small wavenumbers showed that there can be multiple intercepts at
$a=0$
for the buoyancy-opposed regime, where
$\mathit{Pe}<0$
.
Neutral stability curves were analysed for different channel inclination angles and different Péclet numbers, as well as for longitudinal, oblique and transverse rolls. The results have been presented for two distinct types of flow: the buoyancy-assisted regime (positive
$\mathit{Pe}$
), in which the buoyancy force has a component in the direction of the through-flow, and a buoyancy-opposed regime (negative
$\mathit{Pe}$
), in which the same force has a component opposing the through-flow. These two regimes naturally have no meaning when the channel is perfectly horizontal, but also when
$\mathit{Pe}$
approaches infinity, as the magnitude of the buoyant component in the flow direction has no effect on the forced convection flow. For finite Péclet numbers, the buoyancy-assisted regime was shown to have a stabilising effect on the flow, which is greater as the channel inclination increases. In this regime, higher critical Rayleigh numbers are required for the onset of convective instability when larger inclination angles are present. There exists a minimum Péclet number required to achieve an instability of the basic flow. This feature was previously documented for horizontal channels. The current work showed that this minimum value increases with the channel inclination angle,
${\it\phi}$
.
While the buoyancy-assisted basic flow becomes more and more stable as
${\it\phi}$
increases, the buoyancy-opposed flow is always unstable, as the neutral stability curves have their minimum at
$\mathit{Ra}=0$
with
$a\rightarrow 0$
. This feature has an important implication for the design of experiments dealing with horizontal channels. In fact, even a minimal systematic error in setting the horizontal may cause an undesired buoyancy-opposed configuration with a completely altered transition to instability, even with extremely small Rayleigh numbers.
The analysis of stability carried out under the assumption of a perfectly horizontal channel leads to the conclusion that instability takes place if
$Ra$
exceeds its critical value, which depends on
$Pe$
and which cannot be smaller than 119.064 (Barletta Reference Barletta2012). Relaxing the assumption of a perfectly horizontal channel implies that even a vanishingly small inclination could be sufficient for the onset of instability for every value of
$Ra$
, however small. This conclusion casts new light on the linear stability analysis of a base flow. In fact, the linear stability analysis is usually based on the idea of perturbing the initial conditions while the boundary conditions and the geometrical layout of the system are kept constrained. However, in an actual experimental set-up, nothing is perfectly constrained. Hence, prediction of the effects of small, possibly random, perturbations of the initial conditions, but also of the boundary conditions and of the geometrical layout, is definitely important. This kind of analysis represents an essential extension of the usual scope of a linear stability analysis. The result described above is just an example of a fresh view which can be extended to a plethora of different cases.
