Numerical Methods for Solving Some Hydrodynamic Stability Problems

Abstract

In this paper we solve the stability problem of the thermal convection in a thin fluid layer with free-free, slip-slip, and fixed-slip boundary conditions. We use different numerical methods to check their flexibility and accuracy where we use the following numerical methods: finite difference, High order finite difference, $p$ order finite element, Chebyshev collocation-1, Chebyshev collocation-2 and Chebyshev tau. Moreover, we solve the equations of Poiseuille flow in a Brinkman porous media with slip-slip boundary conditions using $p$ order finite element, Chebyshev collocation-1 and Chebyshev collocation-2 methods. The advantage and disadvantage of each method have been discussed. According to our result, we believe that the finite difference and finite element methods are very flexible methods and we can apply them to solve any problem easily. However, the accuracy of the the finite difference and finite element methods is less than the accuracy of the Chebyshev methods.

Introduction

Conventionally, stability calculations involve determining eigenvalues and eigenfunctions, with few of the associated eigenvalue problems solvable analytically. Two powerful existing techniques for finding eigenvalues and eigenfunctions numerically are the compound matrix (cf. [3, 30–32]) and the Chebyshev tau method (cf. [5, 35]). These numerical methods are used in the linear and nonlinear analyses to yield generalised eigenvalue problems of the form,

$$\begin{aligned} A x=\sigma B x \end{aligned}$$

where $A$ and $B$ are matrices and $x$ is some vector, all of which depend upon the system under consideration. The compound matrix method, which belongs to the family of shooting techniques, performs competently for stiff differential equations, with the specific purpose of reducing rounding error, as explored in [11, 41], see also the references therein. The Chebyshev tau technique is a spectral method. This method calculates as many eigenvalues as required as opposed to just one at a time as is done in the compound matrix method.

Orszag [35] introduced the earlier numerical study to the Orr–Sommerfeld problem. He solved this problem by using a direct (tau) spectral method. The other important papers was introduced by [4, 41]. In [4], they analyse the Orr–Sommerfeld equation supplied with homogeneous boundary conditions which contain derivative up to the first order. Straughan and Walker [41] applied these two techniques to linear and nonlinear stability problems for convection in porous media. Their paper provides an excellent summary of the two aforementioned methods. They compare the techniques and highlight the advantages and disadvantages of both when investigating stability problems. Hill and Straughan [22] applied high accuracy Legendre spectral element method for solving second order as well as higher order (especially fourth order) differential eigenvalue problems. Gheorghiu and Pop [10] built up a Chebyshev tau method in order to solve a hydrodynamic stability problem connected with the Marangoni–Plateau–Gibbs convection. The eigenvalue problem was a non-standard one, i.e., two out of the four boundary conditions attached to the Orr–Sommerfeld equation contained derivatives of order two and three. Gheorghiu and Dragomirescu [9] used weighted residuals Galerkin method, weighted residuals Petrov–Galerkin method and the Chebyshev collocation method to solve the linear hydrodynamic stability problem of a convective flow in varying gravity field. However, the main difficulties with spectral methods are the imposing of boundary conditions which involve derivatives of order higher than one and sometimes significant rounding errors appear in the computational results. One of the technique which was used to avoid these difficulties is to choose the spaces of test and trial functions such that these spaces satisfy as many as possible boundary conditions [9, 10].

Recently, Straughan and Harfash introduce important numerical techniques to solve some hydrodynamic stability problems. In [12], we studied the problem of double-diffusive convection in a reacting fluid with a concentration and magnetic field effect based internal heat source. We performed a linear instability analysis and nonlinear stability analysis and then we used the FE method of $p$ order to compute the numerical results. Harfash and Straughan [21] studied the problem of convective movement of a reacting solute in a viscous incompressible fluid occupying a plane layer and subjected to a vertical magnetic field. They solved the linear instability and global nonlinear energy stability systems using the finite differences method. Straughan and Harfash [40] studied a model for Poiseuille flow instability in a porous medium of Brinkman type. In this study, two Chebyshev collocation methods and finite element method have been applied to incorporate the boundary conditions and to provide an independent check. Recently, we used the above numerical schemes to solve the eigenvalue systems of linear instability and nonlinear stability theories for varied convection problems (cf. [13–16]).

In this paper, we focus in two kind of problems. Firstly, we apply the finite difference, High order finite difference, $p$ order finite element, Chebyshev collocation-1, Chebyshev collocation-2 and Chebyshev tau methods to solve the eigenvalue systems of standard thermal convection with free-free, slip-slip, and fixed-slip boundary conditions. [36] showed that, in the case of free-free boundary conditions, we may obtain the analytical result of the Rayleigh number $Ra_{crit} = 27 \pi ^{4}/4$ (see also [1, 3]), thus, we select free-free boundary conditions to check the accuracy of numerical methods. However, slip-slip, and fixed-slip boundary conditions have been selected to check the flexibility of the numerical methods in dealing with these boundary conditions.

Secondly, the focus of attention is on the problem of Poiseuille flow in a channel which is filled with a porous medium saturated with a linear viscous fluid. Due to the invention of man-made materials such as metallic foams which have porosity close to one, and are of much use in the heat transfer industry, the porous-Poiseuille problem is of practical interest.

The classical hydrodynamic problem of stability of Poiseuille flow in a channel is a major one in fluid dynamics, see e.g. [26], chapter 3, [38], chapter 8. There are two major problems correspond with obtaining a meaningful nonlinear energy stability theory for such flows since the nonlinear energy stability threshold is inevitably far away from the linear instability one. Moreover, the Numerical solution of the eigenvalue systems of the equations of the Poiseuille flow stability problem is very difficult, see e.g. [4, 44]. However, this area has attracted considerable interest during the last 40 years. Particularly, we mention to the nonlinear energy stability work of [25], the analyses of time-dependent and time-periodic Poiseuille flows by [6, 7], the work on the problem with slip boundary conditions by [43], the problem of Poiseuille flow of a fluid overlying a porous medium, see [2, 23], and the Poiseuille problem for flow in a channel with one fluid overlying another, see [44]. In this section, we study the linear instability of Poiseuille flow in a porous medium of Brinkman type. The equivalent of the Orr–Sommerfeld eigenvalue problem is solved numerically. Moreover, we discuss the difficulties which corresponding with obtaining the spectrum of the porous Orr–Sommerfeld equation.

Since solving the eigenvalue systems of the equations of the Poiseuille flow stability problem is a difficult numerical problem, we adopt the Chebyshev collocation method to solve the eigenvalue systems and applied two variations of the Chebyshev collocation method in order to incorporate the boundary conditions and also to provide an independent check. As an additional check, we have also employed a finite element method. We have found the finite element method to be more unstable numerically than either of the two collocation methods reported above. However, we include the numerical results of finite element method to compare the performance of the collocation and finite element methods.

The Effect of Boundary Conditions on Convective Instability

Penetrative convection refers to a convective motion, where part of the layer has a tendency to be unstable, which will then induces instability in the rest of the layer (cf. Nield and Bejan [34]). One of the most widely employed mechanisms to describe penetrative convection is internal heating.Aninternal heat source (or sink) can give rise to a situation where one part of a layer is naturally convecting while the other remains stable, allowing penetrative convection to occur. Recently, mathematical models to describe penetrative convection have been developed and analyzed intensely by Harfash [17–20].

Let $x = (x, y, z)$ denote Cartesian coordinates in $\mathbb {R}^{3}$. We consider a fluid contained in the region $\Omega \subset \mathbb {R}^{3}$, which is the infinite layer defined by $\Omega = (- \infty , \infty ) \times (- \infty , \infty ) \times [0, d]$. The behaviour of this fluid is described by the Boussinesq Eqs. (1)–(3), which comprise the Navier–Stokes equations and an energy balance equation:

$$\begin{aligned} \rho \bigg (\frac{\partial \mathbf {v}}{\partial t}+(\mathbf {v} \cdot \nabla ) \mathbf {v} \bigg )= & {} - \nabla p+\mu \, \nabla ^{2} \mathbf {v}-\rho \mathbf {k} \mathrm {g} [1-\alpha (T-T_{m})],\end{aligned}$$

(1)

$$\begin{aligned} \nabla \cdot \mathbf {v}= & {} 0, \end{aligned}$$

(2)

$$\begin{aligned} \frac{\partial T}{\partial t}+\mathbf {v}. \nabla T= & {} \kappa \nabla ^{2} T, \end{aligned}$$

(3)

where $\mathbf {v}=(u, v, w), p$ and $T$ are velocity field, pressure, and temperature, respectively. Additionally, $\alpha$ is the thermal expansion coefficient, $\mu$ is the dynamic viscosity, $\mathrm {g}$ is the acceleration due to gravity, $\rho$ is the density at the reference temperature $T_{m}, \kappa$ is the thermometric conductivity and $\mathbf {k}=(0, 0, 1)$.

$\Omega$ is bounded above by the plane $z = d$ and below by the plane $z = 0$. The temperature at the upper and lower surfaces is kept constant

$$\begin{aligned} T|_{z=0}=T_{L},\,\,\,\,\,\,\,\,\,\,T|_{z=d}=T_{U}, \end{aligned}$$

(4)

for constants $T_{L}>T_{U}$, and thus the layer is heated from below.

Navier [29] proposed a linear boundary condition relating $\mathbf {v}$ to the shear rate, which has become standard in the study of boundary slip problems. Letting the surface $\partial \Omega$ have unit normal $n(x)$ directed out of the fluid, and $\hat{t}(x)$ be either of the vectors tangent to $\partial \Omega$ at $x \in \partial \Omega$, this boundary condition can be expressed as

$$\begin{aligned} v_{i}n_{i}|_{\partial \Omega }= & {} V_{i}n_{i},|_{\partial \Omega } \end{aligned}$$

(5)

$$\begin{aligned} v_{i} \hat{t}_{i}|_{\partial \Omega }= & {} (V_{i}-\lambda \, \epsilon _{ij} \,n_{j})\, \hat{t}_{i},|_{\partial \Omega } \end{aligned}$$

(6)

where $\epsilon = \epsilon (\mathbf {v})$ is the shear strain tensor, and $V_{i} = V_{i}(\partial \Omega )$ is the $i$th component of the local surface velocity. The model is essentially to set the component of $\mathbf {v}$ normal to $\partial \Omega$ to be zero, thus imposing a condition of zero flux across the surface, while setting the two tangential components of $\mathbf {v}$ proportional to the corresponding components of shear stress. We denote the constant of proportionality $\lambda \ge 0$, which has the dimension of length, and it can be seen that $\lambda = 0$ in (5) and (6) recovers the no-slip boundary condition.

We now apply the boundary conditions (5)–(6) to the Boussinesq model. Since the fluid is confined to $\Omega$, from (5) we impose,

$$\begin{aligned} w=0,\quad \mathrm {at}\quad z=0, d, \end{aligned}$$

(7)

and we note that since there is no variation of $w$ in the surfaces $\partial \Omega _{L}$ and $\partial \Omega _{U}$, we must have

$$\begin{aligned} w_{,x}=w_{,y}=0,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\mathrm {at}\,\,\,\,\,\,\,\,\,\,\,\,z=0, d. \end{aligned}$$

(8)

Let $\lambda _{L}$ be the slip length associated with the fluid-solid interface at $\partial \Omega _{L}$, and define $\lambda _{U}$ similarly. Then, from (6) we have

$$\begin{aligned} u-\lambda _{L}\,u_{,z}=0,\quad v-\lambda _{L}\,v_{,z}=0,\quad \mathrm {at} \quad z=0, \end{aligned}$$

(9)

$$\begin{aligned} u-\lambda _{U}\,u_{,z}=0,\quad v-\lambda _{U}\,v_{,z}=0,\quad \mathrm {at} \quad z=d. \end{aligned}$$

(10)

We note that these boundary conditions allow the zero solution $\mathbf {v} = 0$, which represents a fluid at rest. Let us now consider the basic steady state solution $(\bar{\mathbf {v}}, \bar{p}, \bar{T})$ of the system, where, as there is no fluid flow, $\bar{\mathbf {v}} \equiv 0$. Utilizing the boundary conditions and assuming that the basic steady state solutions are functions of $z$ only

$$\begin{aligned} \bar{T}=-\beta z+T_{L}, \end{aligned}$$

(11)

where $\beta =(T_{L}-T_{U})/d$. The steady pressure $\bar{p}$ may then be found from (1) which reduces to

$$\begin{aligned} - \frac{1}{\rho } \frac{d \bar{p}}{d z }= \mathrm {g} [1-\alpha (T-T_{m})]. \end{aligned}$$

(12)

To study the stability of (1)–(3), we introduce a perturbation $(\mathbf {u}, \pi , \theta )$ to the steady state solution $(\bar{\mathbf {v}}, \bar{p}, \bar{T})$, where

$$\begin{aligned} v_{i}=\bar{v}_{i}+u_{i},\,\,\,\,\,\,\,\,\,\,\,\, p=\bar{p}+ \pi ,\,\,\,\,\,\,\,\,\,\,\,\, T=\bar{T}+\theta . \end{aligned}$$

Using (11), (12) the perturbed system is

$$\begin{aligned} \frac{\partial \mathbf {u}}{\partial t}+(\mathbf {u} \cdot \nabla ) \mathbf {u}= & {} - \frac{1}{\rho } \nabla \pi +\nu \nabla ^{2} \mathbf {u}+\mathrm {g} \alpha \mathbf {k} \theta , \end{aligned}$$

(13)

$$\begin{aligned} \frac{\partial \theta }{\partial t}+\mathbf {u} \cdot \nabla \theta= & {} \beta w+\kappa \nabla ^{2} \theta , \end{aligned}$$

(14)

where $u_{i}$ is solenoidal, i.e. $\nabla \cdot \mathbf {u}=0$ and $\nu =\mu /\rho$.

We now introduce non-dimensionalised variable with scaling of

$$\begin{aligned} \mathbf x= & {} \mathbf x ^{*} d,\,\,\,\,\ t=t^{*}\frac{d^{2}}{\nu },\,\,\,\,\ \mathbf u =U \mathbf u ^{*},\,\,\,\,\ \theta =T^{\sharp } \theta ^{*},\,\,\,\,\ \pi =P \pi ^{*},\,\,\,\,\ U=\frac{\nu }{d},\,\,\,\ P=\frac{\rho \, \nu ^{2}}{d^{2}},\\ T^{\sharp }= & {} U \sqrt{\frac{\nu \, \beta }{\kappa \alpha \mathrm {g}}},\,\,\,\,\ R=\sqrt{\frac{\alpha \mathrm {g}\, d^{4} \, \beta }{\kappa \, \nu }},\,\,\,\,\ P_{r}=\frac{\nu }{\kappa }. \end{aligned}$$

Here $Pr$ is the Prandtl number and $Ra=R^{2}$ is the Rayleigh number. With this scaling the non-dimensional form of (13)–(14) becomes (omitting the stars for case of notation)

$$\begin{aligned} \frac{\partial \mathbf {u}}{\partial t}+(\mathbf {u} \cdot \nabla ) \mathbf {u}= & {} - \nabla \pi + \nabla ^{2} \mathbf {u}+ \mathbf {k} R \theta , \end{aligned}$$

(15)

$$\begin{aligned} \nabla \cdot \mathbf {u}= & {} 0, \end{aligned}$$

(16)

$$\begin{aligned} P_{r}\left( \frac{\partial \theta }{\partial t}+\mathbf {u} \cdot \nabla \theta \right)= & {} R w+ \nabla ^{2} \theta . \end{aligned}$$

(17)

The spatial domain is now $\{(x,y) \in \mathbb {R}^{2}\} \times \{z \in (0,1)\}$. The perturbed boundary conditions are given by

$$\begin{aligned}&u-\lambda _{L}\,u_{,z}=0,\qquad v-\lambda _{L}\,v_{,z}=0,\qquad \mathrm {at}\qquad z=0, \end{aligned}$$

(18)

$$\begin{aligned}&u-\lambda _{U}\,u_{,z}=0,\qquad v-\lambda _{U}\,v_{,z}=0,\qquad \mathrm {at}\qquad z=1, \end{aligned}$$

(19)

$$\begin{aligned}&\quad \theta =0 ,\,\,\,\,\,\,\, \mathrm {on} \,\,\,\,\,\,\, z=0,1. \end{aligned}$$

(20)

We have reduced the problem of finding conditions for the onset of convection in our fluid layer to that of investigating the stability of the basic steady state solutions with respect to perturbations $\mathbf {u}, \pi , \theta$ as defined above. In this way we aim to find, for fixed $\lambda _{L}$ and $\lambda _{U}$, the critical Rayleigh number $R^{2}_{crit}(\lambda _{L}, \lambda _{U})$ such that solutions to (15)–(20) decay over time for $R <R_{crit}$ and grow for $R > R_{crit}$, regardless of the initial data .

We do this by showing that there exists $R_{L}$ such that thermal instability will occur for $R > R_{L}$, and $R_{E}$ such that $R < R_{E}$ guarantees stability of the the basic steady state solutions. It has been shown that $R_{E} = R_{L} = R_{crit}$ for system (15)–(17) for no-slip boundary conditions (see [27, 28]) and for slip boundary conditions (see [42]), thus, it is enough to solve the linear system to find the linear and nonlinear threshold and hence we can solve the system with solutions of the form of single Fourier modes.

The linearised equations are obtained from (15)–(17) by omitting a nonlinear terms. The resulting linearized equations possess solutions of type

$$\begin{aligned} u_{i}(\mathbf x ,t)=u_{i}(\mathbf x ) e^{\sigma t},\quad \theta (\mathbf x ,t)=\theta (\mathbf x ) e^{\sigma t},\quad \pi (\mathbf x ,t)=\pi (\mathbf x ) e^{\sigma t}, \end{aligned}$$

where $\sigma$ is the growth rate and a complex constant. $u_{i}(\mathbf x ), \phi (\mathbf x ), \pi (\mathbf x )$ then satisfy

$$\begin{aligned} -\pi _{,\,i}+ \Delta u_{i}+k_{i} R\, \theta= & {} \sigma u_{i}, \end{aligned}$$

(21)

$$\begin{aligned} R w+\Delta \theta= & {} \sigma P_{r} \,\theta . \end{aligned}$$

(22)

Taking the double $curl$ of (21), using the third component, (and the fact that $u$ is solenoidal) we have

$$\begin{aligned} \Delta ^{2} w+R \, \Delta ^{*} \theta = \sigma \Delta w, \end{aligned}$$

(23)

where $\Delta ^{*}={\partial ^{2}}/{\partial x^{2}}+ {\partial ^{2}}/{\partial y^{2}}, D=d/dz$. We now introduce normal modes of the form $w=W(z) f(x,y)$, and $\theta =\Theta (z) f(x,y)$ where $f(x,y)$ is a plan-form which tiles the plane $(x,y)$ with

$$\begin{aligned} \Delta ^{*} f = -a^{2} f. \end{aligned}$$

(24)

The plan-forms represent the horizontal shape of the convection cells formed at the onset of instability. These cells from a regular horizontal pattern tiling the $(x,y)$ plane, where the wavenumber $a$ is a measure of the width of the convection cell. Using (24), and applying the normal mode representations to (22) and (23) we find

$$\begin{aligned} (D^{2}-a^{2})^{2}W-a^{2} R \, \Theta= & {} \sigma (D^{2}-a^{2})W, \end{aligned}$$

(25)

$$\begin{aligned} (D^{2}-a^{2})\Theta + R W= & {} \sigma P_{r} \Theta , \end{aligned}$$

(26)

where $D=d/dz,$ and $z\in (0, 1)$. It is easy to show that $\sigma \in \mathbb {R}$, and therefore the principle of “exchange of stabilities” applies to the linearized system, and thus it is enough to solve system (25)–(26) with $\sigma =0$ i.e we shall solve the following system

$$\begin{aligned} (D^{2}-a^{2})^{2}W= & {} a^{2} R \, \Theta , \end{aligned}$$

(27)

$$\begin{aligned} (D^{2}-a^{2})\Theta= & {} - R W. \end{aligned}$$

(28)

The boundary conditions which $w$ must satisfy are derived from the conditions (18) and (19) and the incompressibility condition $u_{,x}+v_{,y}+w_{,z}=0,$ are

$$\begin{aligned} w(0)=w(1)=0,\,\,\,\,\,\,\,\,\,\lambda _{L} w{^{\prime \prime }}(0)-w{^\prime }(0)=0,\,\,\,\,\,\,\,\,\,\lambda _{U} w{^{\prime \prime }}(1)+w{^\prime }(1)=0, \end{aligned}$$

(29)

and the boundary conditions for $\Theta$ are

$$\begin{aligned} \Theta (0)=\Theta (1)=0. \end{aligned}$$

(30)

Equations (27) and (28) are the classic stability equations for the Bénard problem. We note that in the limit $\lambda _{U} \rightarrow 0$ we obtain from (29) the no-slip boundary condition at $z = 1,$ while for $\lambda _{U} \rightarrow \infty$ we recover the free boundary condition (and similarly for $\lambda _{L}$ at z $=$ 0). In the next section, six numerical methods are used to solve this system (27) and (28).

Numerical Methods for Thermal Convection Eigenvalue System

Chebyshev Tau

Equation (27) has a fourth order derivative. Dongarra et al. [4] show that high order differentiation matrices, for instance in this case the $D^{4}$ matrix, can introduce significant round off errors. Therefore we use what is described in the literature as a $D^{2}$ method, and make the substitution. Letting $\lambda _{U}$ and $\lambda _{L} \rightarrow \infty$, then we obtain from (29) the free boundary condition at $z = 0$ and $z = 1$. We introduce new function $A=(D^{2}-a^{2})W$. Rewriting Eqs. (27) and (28) in terms of the new variable $A$, we require to solve

$$\begin{aligned} (D^{2}-a^{2})W= & {} A, \end{aligned}$$

(31)

$$\begin{aligned} (D^{2}-a^{2})A= & {} a^{2} R \, \Theta , \end{aligned}$$

(32)

$$\begin{aligned} (D^{2}-a^{2})\Theta= & {} -R W , \end{aligned}$$

(33)

with boundary conditions

$$\begin{aligned} \Theta = W=A= 0, \,\,\,\mathrm {at} \,\,z=0,1. \end{aligned}$$

(34)

To employ the Chebyshev tau technique, system (31)–(34) is converted to the Chebyshev domain $(-1,1)$, and then $W, A$ and $\Phi$ are written as a finite series of Chebyshev polynomials

$$\begin{aligned} W=\sum _{k=0}^{N+2}W_{k} T_{k}(z), \,\,\,\,\,\,\,\, A=\sum _{k=0}^{N+2}A_{k} T_{k}(z), \,\,\,\,\,\,\,\, \Theta =\sum _{k=0}^{N+2}\Theta _{k} T_{k}(z). \end{aligned}$$

The weighted inner product of each equation is taken with some $T_{k}$ and the orthogonality of the Chebyshev polynomial is utilised to form the generalised eigenvalue problem

$$\begin{aligned} \left( \begin{array}{ccc} 4D_{2}-a^{2}I &{} -I &{} O \\ O &{} 4D_{2}-a^{2}I &{} O \\ O &{} O &{} 4D_{2}-a^{2}I \\ \end{array} \right) \mathbf Q =R \left( \begin{array}{ccc} O &{} O &{} O \\ O &{} O &{} a^{2}I\\ -I &{} O &{} O \\ \end{array} \right) \mathbf Q , \end{aligned}$$

where $\mathbf Q =\left( \begin{array}{ccc} \hat{W}, &{} \hat{A}, &{} \hat{\Phi }, \\ \end{array} \right) ^{T}, \hat{W}=(W_{0},\ldots ,W_{N+2})^{T}, \hat{A}=(A_{0},\ldots ,A_{N+2})^{T}, \hat{\Theta }=(\Theta _{0},\ldots ,\Theta _{N+2})^{T}$, $D_{2}$ is the Chebyshev representation of $D^{2}$, and $I$ is the identity matrix. Using the boundary conditions and the fact that $T_{n}(\pm 1)=(\pm 1)^{n}$ we remove the last two rows of each $(N+2)\times (N+2)$ block and replace these rows by the discrete form of the boundary conditions (34) which have the form

$$\begin{aligned} BC_{1}:\,\,\,\,\, \sum _{n=0}^{N} W_{n}= & {} 0, \end{aligned}$$

(35)

$$\begin{aligned} BC_{2}:\,\,\,\,\,\sum _{n=0}^{N} (-1)^n W_{n}= & {} 0, \end{aligned}$$

(36)

$$\begin{aligned} BC_{3}:\,\,\,\,\, \sum _{n=0}^{N} A_{n}= & {} 0, \end{aligned}$$

(37)

$$\begin{aligned} BC_{4}:\,\,\,\,\,\sum _{n=0}^{N} (-1)^n A_{n}= & {} 0, \end{aligned}$$

(38)

$$\begin{aligned} BC_{5}:\,\,\,\,\, \sum _{n=0}^{N} \phi _{n}= & {} 0, \end{aligned}$$

(39)

$$\begin{aligned} BC_{6}:\,\,\,\,\,\sum _{n=0}^{N} (-1)^n \phi _{n}= & {} 0. \end{aligned}$$

(40)

Then, the final eigenvalue system has the following form

$$\begin{aligned} \left( \begin{array}{ccc} 4D^{2}-a^{2}I &{} -I &{} O \\ BC_{1} &{} 0\ldots 0 &{} 0\ldots 0\\ BC_{2} &{} 0\ldots 0&{} 0\ldots 0\\ O &{} 4D^{2}-a^{2}I &{} O\\ 0\ldots 0 &{} BC_{3}&{} 0\ldots 0\\ 0\ldots 0 &{} BC_{4}&{} 0\ldots 0\\ O &{} O &{} 4D^{2}-a^{2}I \\ 0\ldots 0 &{} 0\ldots 0&{} BC_{5}\\ 0\ldots 0 &{} 0\ldots 0&{} BC_{6}\\ \end{array} \right) \mathbf Q =R \left( \begin{array}{ccc} O &{} O &{} O\\ 0\ldots 0 &{} 0\ldots 0&{} 0\ldots 0\\ 0\ldots 0 &{} 0\ldots 0&{} 0\ldots 0\\ O &{} O &{} a^{2} I \\ 0\ldots 0 &{} 0\ldots 0&{} 0\ldots 0\\ 0\ldots 0 &{} 0\ldots 0&{} 0\ldots 0\\ -I &{} O &{} O \\ 0\ldots 0 &{} 0\ldots 0&{} 0\ldots 0\\ 0\ldots 0 &{} 0\ldots 0&{} 0\ldots 0\\ \end{array} \right) \mathbf Q . \end{aligned}$$

(41)

We solve this eigenvalue system by passing the above matrices to the Matlab’s QZ subroutine. For an initial value of the wave number $a,$ we obtain a spectrum of eigenvalues then we select the smallest eigenvalue $R_{m}$ such that the fluid is stable for all $R < R_{m}$. We then repeat for other values of $a$, so that we may build the neutral curve of $Ra = R^{2}_{m}$ against $a$, along which the fluid is neutrally stable. By iterating over $a$ we are thus able to obtain the critical Rayleigh number $Ra_{crit}$ as the minimum value on this curve, which occurs at wavenumber $a_{crit}$.

Also, we can adopt another technique to apply the boundary conditions where we can find the values of $W_{N+1}, W_{N+2}, A_{N+1}, A_{N+2} ,\Phi _{N+1}$ and $\Phi _{N+2}$ from the boundary conditions then we substitute these values in the system and thus we remove the last two rows. The new system have the order $(N+1) \times (N+1)$ instead of $(N+3) \times (N+3)$.

However, if $\lambda _{U}$ and $\lambda _{L}$ are fixed finite numbers, and thus our boundary conditions are slip-slip. The solution of (31)–(33) with slip-slip boundary conditions produces the same eigenvalue system to (41). The boundary conditions $BC1, BC2, BC5, BC6$ still the same as in (41), while the changes will be in the $BC3, BC4$ which have the following new forms:

$$\begin{aligned} BC_{3}:\,\,\,\,\,\sum _{n=0}^{N} 2N_{0} A_{n}+\sum _{n=0}^{N}n^{2} W_{n}= & {} 0, \end{aligned}$$

(42)

$$\begin{aligned} BC_{4}:\,\,\,\,\,\sum _{n=0}^{N} 2N_{0}(-1)^{n} A_{n}+\sum _{n=0}^{N}(-1)^{n+1}n^{2} W_{n}= & {} 0, \end{aligned}$$

(43)

Finally, letting $\lambda _{U} \rightarrow 0$ and $\lambda _{L}$ a fixed finite numbers, i.e. we have fixed-slip boundary conditions. We introduce new function $\Psi =DW$, then system (27) and (28) can written in the following form:

$$\begin{aligned} DW-\Psi= & {} 0, \end{aligned}$$

(44)

$$\begin{aligned} D^{3}\Psi -2a^{2}DA+a^{4}W-a^{2} R \, \Theta= & {} 0 \end{aligned}$$

(45)

$$\begin{aligned} (D^{2}-a^{2})\Phi + R W= & {} 0, \end{aligned}$$

(46)

hence, the boundary conditions have the form

$$\begin{aligned} W= & {} \phi =0,\,\,\,\,\,\,\,z=0, 1, \end{aligned}$$

(47)

$$\begin{aligned} N_{0}D\Psi +\Psi= & {} 0,\,\,\,\,\,z=1, \end{aligned}$$

(48)

$$\begin{aligned} \Psi= & {} 0,\,\,\,\,\,z=0, \end{aligned}$$

(49)

then, we can apply the same procedure which are used for free-free and slip-slip boundary conditions to produce the eigenvalue system.

Finite Difference Scheme

Standard Finite Difference

The standard second and fourth order central difference operators at grid point $i$ can be written as:

$$\begin{aligned} \left. \begin{array}{l} \delta ^2 u _{i} = \frac{{u _{i + 1} - 2 u _{i} + u _{i - 1} }}{{h^{2} }}, \\ \delta ^4 u _{i} = \frac{{u _{i + 2}-4u _{i + 1} +6 u _{i} -4 u _{i - 1}+u _{i - 2} }}{{h^{4} }}. \\ \end{array} \right\} \end{aligned}$$

(50)

The second and the fourth order derivatives for the function $u$ at grid point $i$ can be approximated by a second order accuracy as

$$\begin{aligned} \left. \begin{array}{l} \left. {\frac{{d ^2 u }}{{d z^2 }}} \right| _{i} = \delta ^{2} u _{i} - \frac{{h^{2} }}{{12}}\frac{{d ^{4} u }}{{d z^{4} }} + O(h^{4} ), \\ \left. {\frac{{d ^4 u }}{{d z^4 }}} \right| _{i} = \delta ^{4} u _{i} - \frac{{h^2 }}{{6}}\frac{{d ^{6} u }}{{d z^{6} }} + O(h^{4} ). \\ \end{array} \right\} \end{aligned}$$

(51)

By using these finite difference approximations, (27) and (28) can be discretized at a given grid point $i$ respectively as,

$$\begin{aligned} \delta _{z}^{4}W_{i}- 2a^{2}\delta _{z}^{2}W_{i}+ a^{4}W_{i}-a^{2} R \Theta _{i}= & {} 0, \end{aligned}$$

(52)

$$\begin{aligned} \delta _{z}^{2} \Theta _{i}-a^{2} \Theta _{i}+R W_{i}= & {} 0, \end{aligned}$$

(53)

The boundary conditions $D_{z}^{2}W=0$ at $z=0, 1$ are approximated using finite difference technique as $W_{-1}=-W_{1}$ and $W_{N+1}=-W_{N-1}$. In this manner, Eqs. (52) and (53) and the fixed boundary conditions lead to the finite difference equations

$$\begin{aligned}&\frac{W_{i+2}}{h^{4}}-\left( \frac{4}{h^{4}}+\frac{2 a^{2}}{h^{2}}\right) W_{i+1}+\left( \frac{6}{h^{4}}+\frac{4 a^{2}}{h^{2}}+a^{4}\right) W_{i}-\left( \frac{4}{h^{4}} +\frac{2 a^{2} }{h^{2}}\right) W_{i-1}\nonumber \\&\qquad \quad +\frac{W_{i-2}}{h^{4}}-R a^{2}\, \Theta _{i}=0, \quad i=2,\ldots ,N-2,\nonumber \\&\qquad \quad \frac{\Theta _{i+1}}{h^{2}}- \left( \frac{2}{h^{2}}+a^{2}\right) \Theta _{i}+\frac{\Theta _{i-1}}{h^{2}}+R W_{i}=0, \end{aligned}$$

(54)

$$\begin{aligned}&\qquad \quad i=1,\ldots ,N-1,\nonumber \\&\quad \frac{W_{3}}{h^{4}}-\left( \frac{4}{h^{4}}+\frac{2 a^{2}}{h^{2}}\right) W_{2}+\left( \frac{5}{h^{4}}+\frac{4 a^{2}}{h^{2}}+a^{4}\right) W_{1}+R a^{2} \Phi _{1} =0, \end{aligned}$$

(55)

which is the equation obtained from (52) with $i=1,$ and

$$\begin{aligned} \left( \frac{5}{h^{4}}+\frac{4 a^{2}}{h^{2}}+a^{4}\right) W_{N-1}-\left( \frac{4}{h^{4}}+\frac{2 a^{2}}{h^{2}}\right) W_{N-2}+\frac{W_{N-3}}{h^{4}}+R a^{2} \Phi _{N-1}=0, \end{aligned}$$

(56)

which arises from (52) with $i=N-1$.

For slip-slip boundary conditions, we adapt Eqs. (55) and (56). At $i=1$, we substitute $W_{0}=0$ and $W_{-1}=-(N0/h^{2}-1/2h)/(N0/h^{2}+1/2h)W_{1}$ in (54) to get the equivalent equation to (55). At $i=N$, we substitute $W_{N}=0$ $W_{N+1}=-(N0/h^{2}-1/2h)/(N0/h^{2}+1/2h)W_{N-1}$ in (54) to get the equivalent equation to (56). For fixed-slip boundary conditions, At $i=N$, we substitute $W_{N}=0$ $W_{N+1}=W_{N-1}$ in (54) to get the equivalent equation to (56). At $i=1$ the equation still the same with slip-slip boundary conditions.

High Order Finite Difference

The main idea of the high order finite difference scheme is to find the values of truncation errors from the original differential equation and substitute these values in the finite difference formula. In this scheme we can reduce the order of truncation errors.

The compact finite difference method is one of the methods which have been used to increase the accuracy of the numerical solutions of the partial differential equations. Therefore, we can use many techniques to derive the compact finite difference scheme for any problem. Spotz [37] used the Taylor series to drive a fourth and sixth compact finite difference for 2D and 3D Poisson equation.

Utilising the approach of [37], we can easily find the value of $D^{6}W$ and $D^{4}\Theta$ as follows

$$\begin{aligned} D^{6}W= & {} (2 a^{2}-M^{2})D^{4}W-a^{4} D^{2}W+a^{2} R D^{2} \Theta , \end{aligned}$$

(57)

$$\begin{aligned} D^{4}\Theta= & {} a^{2}D^{2}\Phi - R D^{2} W. \end{aligned}$$

(58)

Substituting the values of $D^{6}W$ and $D^{4}\Theta$ in (51), and then approximate $D^{2}W, D^{4}W$ and $D^{2}\Theta$ by using standard second order finite difference we have the following fourth order finite difference formula

$$\begin{aligned} \left. {\frac{{d ^{2} W }}{{d z^{2} }}} \right| _{i}= & {} \delta ^{2} W _{i} - \frac{{h^{2} }}{{12}}\delta _z^{4} (W _{i}+O(h^{2})) + O(h^{4} ), \end{aligned}$$

(59)

$$\begin{aligned} \left. {\frac{{d ^{2} \Theta }}{{d z^{2} }}} \right| _{i}= & {} \delta ^{2} \Theta _{i} - \frac{{h^{2} }}{{12}} \left\{ a^{2}(\delta ^{2}\Phi _{i}+O(h^{2}))- R (\delta ^{2} W_{i}+O(h^{2})) \right\} + O(h^{4} ), \end{aligned}$$

(60)

$$\begin{aligned} \left. {\frac{{d ^{4} W }}{{d z^{4} }}} \right| _{i}= & {} \delta ^{4} W _{i} - \frac{{h^{2} }}{{6}} \left\{ 2 a^{2}(\delta ^{4}W_{i}+O(h^{2}))-a^{4} (\delta ^{2}W_{i}+O(h^{2}))+a^{2} R (\delta ^{2} \Theta _{i}+O(h^{2}))\right. \nonumber \\&+ O(h^{4}). \end{aligned}$$

(61)

It clear that the overall truncation error will be of $O(h^{4})$. Using (59)–(61) finite difference approximations, (27) and (28) can be approximated at a grid point $i$ respectively as,

$$\begin{aligned} r_{2}\, \delta ^{4}W_{i}+ r_{1}\, \delta ^{2}W_{i}+ a^{4}W_{i}-a^{2} R \left( 1+\frac{h^{2}}{6}\delta ^{2}\right) \Theta _{i}= & {} 0, \end{aligned}$$

(62)

$$\begin{aligned} R\left( 1+\frac{\Delta z^{2}}{12} \delta ^{2}\right) W_{i}+r_{3}\, \delta _{z}^{2} \Theta _{i}-a^{2} \Phi _{i}= & {} 0, \end{aligned}$$

(63)

where $r_{1}=(h^{2}/6)\,\,a^{4}-2a^{2},\,\,r_{2}=1-(h^{2}/6)\,a^{2}$, $r_{3}=1-(h^{2}/12) \,a^{2}$ and $r_{4}=1-(h^{2}/6) \,a^{2}$. Thus, for free-free boundary conditions, (62) and (63) produce the following high order finite difference equations

$$\begin{aligned}&\frac{r_{2}}{h^{4}}W_{i+2}+\left( -\frac{4r_{2}}{h^{4}}+\frac{r_{1}}{h^{2}}\right) W_{i+1}+\left( \frac{6r_{2}}{h^{4}}-\frac{2r_{1}}{h^{2}}+a^{4}\right) W_{i}+\left( -\frac{4r_{2}}{h^{4}} +\frac{r_{1}}{h^{2}}\right) W_{i-1}\nonumber \\&\quad +\frac{r_{2}}{h^{4}}W_{i-2}-R a^{2}\,\left[ \frac{1}{6}\Theta _{i-1}+\frac{2}{3}\Theta _{i}+\frac{1}{6}\Theta _{i+1}\right] =0, \end{aligned}$$

(64)

$$\begin{aligned}&\quad i=2,\ldots ,N-2,\nonumber \\&-R\, \left[ \frac{-1}{12}W_{i-1}-\frac{5}{6}W_{i}-\frac{1}{12}W_{i+1}\right] +r_{3}\frac{\Theta _{i-1}}{h^{2}}-\left( \frac{2r_{3}}{h^{2}}+a^{2}\right) \Theta _{i}+r_{3}\frac{\Theta _{i+1}}{h^{2}} =0, \quad \quad \end{aligned}$$

(65)

$$\begin{aligned}&\quad i=1,\ldots ,N-1,\nonumber \\&\frac{r_{2}}{h^{4}}W_{3}+\left( -\frac{4r_{2}}{h^{4}}+\frac{r_{1}}{h^{2}}\right) W_{2}+\left( \frac{5r_{2}}{h^{4}}-\frac{2r_{1}}{h^{2}}+a^{4}\right) W_{1} -R a^{2}\,\left[ \frac{2}{3}\Theta _{1}+\frac{1}{6}\Theta _{2}\right] =0,\nonumber \\ \end{aligned}$$

(66)

which is the equation obtained from (62) with $i=1,$ and

$$\begin{aligned}&\left( \frac{5r_{2}}{h^{4}}-\frac{2r_{1}}{h^{2}}+a^{4}\right) W_{N-1}+\left( -\frac{4r_{2}}{h^{4}}+\frac{r_{1}}{h^{2}}\right) W_{N-2}+\frac{r_{2}}{h^{4}}W_{N-3}\nonumber \\&\quad -R a^{2}\,\left[ \frac{1}{6}\Theta _{N-2}+\frac{2}{3}\Theta _{N-1}\right] =0, \end{aligned}$$

(67)

which arises from (62) with $i=N-1$.

For slip-slip boundary conditions, we change the Eqs. (66) and (67) by new equations. At $i=1$, we substitute $W_{0}=0$ and $W_{-1}=-(N0/h^{2}-1/2h)/(N0/h^{2}+1/2h)W_{1}$ in (64) to get the equivalent equation to (66). At $i=N$, we substitute $W_{N}=0$ $W_{N+1}=-(N0/h^{2}-1/2h)/(N0/h^{2}+1/2h)W_{N-1}$ in (64) to get the equivalent equation to (67). For fixed-slip boundary conditions, At $i=N$, we substitute $W_{N}=0$ $W_{N+1}=W_{N-1}$ in (64) to get the equivalent equation to (67) . At $i=1$ the equation still the same with slip-slip boundary conditions.

Generally, the finite difference and high order finite difference schemes produce a generalized matrix eigenvalue problem of form

$$\begin{aligned} \Lambda \, \Sigma = R\, \Xi \, \Sigma , \end{aligned}$$

(68)

where $\Sigma$ is the eigenfunction vector, the Matrices $\Lambda$ and $\Xi$ have different values according to each case and $R$ represent the eigenvalues of our problem.

Finite Element Method

In this section we will discuss the finite element method which is used to solve systems (27) and (28). Firstly, we introduce a new variable $\Upsilon =D^{2}W$, Therefore, system (27) and (28) become as follows

$$\begin{aligned} D^{2}W= & {} \Upsilon ,\nonumber \\ D^{2}\Upsilon -2a^{2}\Upsilon +a^{4}W-a^{2} R \Theta= & {} 0,\nonumber \\ (D^{2}-a^{2})\Theta + R W= & {} 0, \end{aligned}$$

(69)

with free-free boundary conditions

$$\begin{aligned} W=\Upsilon =\Theta =0, \,\,\,\, \mathrm {at}\,\,\,\,\, z=0,1, \end{aligned}$$

(70)

or with slip-slip boundary conditions

$$\begin{aligned} \begin{array}{l@{\quad }l} W=\Theta =0, &{}\mathrm {at}\quad z=0,1,\\ \lambda _{L}\Upsilon -DW=0,&{} \mathrm {at}\quad z=0,\\ \lambda _{U}\Upsilon +DW=0,&{} \mathrm {at}\quad z=1, \end{array} \end{aligned}$$

(71)

or with fixed-slip boundary conditions

$$\begin{aligned} \begin{array}{l@{\quad }l} W=\Theta =0,&{} \mathrm {at}\quad z=0,1,\\ DW=0,&{} \mathrm {at}\quad z=0,\\ \lambda _{U}\Upsilon +DW=0,&{}\mathrm {at}\quad z=1, \end{array} \end{aligned}$$

(72)

Firstly, we divided the period $0\le x \le 1$ into $n$ elements, each element $e$ having $p+1$ nodes $n_{1},n_{2},\ldots ,n_{p+1}$. In term of its $p+1$ nodal values, the variables $W, \Upsilon , \Theta$ may be uniquely interpolated as a polynomial of $p+1$ order, the interpolation being given by

$$\begin{aligned} W^{e}=N^{e}\delta _{W}^{e},\,\,\,\,\,\Upsilon ^{e}=N^{e}\delta _{\Upsilon }^{e},\,\,\,\,\,\Theta ^{e}=N^{e}\delta _{\Theta }^{e}, \end{aligned}$$

(73)

where $\delta _{W}^{e}=\{W_{n_{1}}, W_{n_{2}}\ldots ,W_{n_{p+1}}\}$, $\delta _{\Upsilon }^{e}=\{\Upsilon _{n_{1}}, \Upsilon _{n_{2}},\ldots ,\Upsilon _{n_{p+1}}\}$ and $\delta _{\Theta }^{e}=\{\Theta _{n_{1}}, \Theta _{n_{2}},\ldots ,\Theta _{n_{p+1}}\}$, and the shape function matrix is

$$\begin{aligned} N^{e}=[N_{n_{1}}^{p},N_{n_{2}}^{p},\ldots ,N_{n_{p+1}}^{p}]. \end{aligned}$$

Therefore the over-all finite element approximation is given by

$$\begin{aligned} W=\sum _{e=1}^{n} W^{e},\,\,\,\,\,\Upsilon =\sum _{e=1}^{n} \Upsilon ^{e},\,\,\,\,\,\Theta =\sum _{e=1}^{n} \Theta ^{e}, \end{aligned}$$

(74)

The variational formulations of (69) is

$$\begin{aligned} \mathrm {Minimize}\,\,\,\, I_{W}[W]= & {} \int _{0}^{1}(-(DW)^{2}-2 \Upsilon W)dz,\end{aligned}$$

(75)

$$\begin{aligned} \mathrm {Minimize}\,\,\,\, I_{\Upsilon }[\Upsilon ]= & {} \int _{0}^{1}(-(D \Upsilon )^{2}-2a^{2} \Upsilon ^{2}+2a^{4}W \Upsilon -2a^{2}R \Theta \Upsilon ))dz\nonumber \\&+ \Upsilon (1) D \Upsilon (1)-\Upsilon (0) D \Upsilon (0), \end{aligned}$$

(76)

$$\begin{aligned} \mathrm {Minimize}\,\,\,\, I_{\Theta }[\Theta ]= & {} \int _{0}^{1}(-(D \Theta )^{2}-a^{2} \Theta ^{2}+2 R W \Theta )dz.\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \end{aligned}$$

(77)

Substitution of (74) into variational formulations (75)–(77) gives

$$\begin{aligned} I_{W}[W]= & {} \sum _{e=1}^{n} \int _{e}(-(DW^{e})^{2}-2\Upsilon ^{e}W^{e})dx=\sum _{e=1}^{n}I_{W}^{e}, \end{aligned}$$

(78)

$$\begin{aligned} I_{\Upsilon }[\Upsilon ]= & {} \int _{1}(-(D \Upsilon ^{1})^{2}-2a^{2}(\Upsilon ^{1})^{2}+2a^{4}W^{1} \Upsilon ^{1}-2a^{2}R \Theta ^{1} \Upsilon ^{1})dx-\Upsilon ^{1}(0) D \Upsilon ^{1}(0)\nonumber \\&+\sum _{e=2}^{n-1} \int _{e}(-(D \Upsilon ^{e})^{2}-2a^{2}(\Upsilon ^{e})^{2}+2a^{4}W^{e} \Upsilon ^{e}-2a^{2}R \Theta ^{e} \Upsilon ^{e})dx \nonumber \\&+\int _{n}(-(D \Upsilon ^{n})^{2}-2a^{2}(\Upsilon ^{n})^{2}+2a^{4}W^{n} \Upsilon ^{n}-2a^{2}R \Theta ^{n} \Upsilon ^{n})dx+\Upsilon ^{n}(1) D \Upsilon ^{n}(1) \nonumber \\= & {} \sum _{e=1}^{n}I_{\Upsilon }^{e}, \end{aligned}$$

(79)

$$\begin{aligned} I_{\Theta }[\Theta ]= & {} \sum _{e=1}^{n} \int _{e}(-(D \Theta ^{e})^{2}-a^{2} (\Theta ^{e})^{2}+2 R f W^{e} \Theta ^{e})dx=\sum _{e=1}^{n}I_{\Theta }^{e}, \end{aligned}$$

(80)

here, we use the fact that the functions $W^{e}, \Upsilon ^{e}, \Theta ^{e}$ are equivalent to zero outside the element $e$. Using (78)–(80), we obtain

$$\begin{aligned} I_{W}[W]= & {} I_{W}(W_{1}, W_{2},\ldots , W_{m}),\\ I_{\Upsilon }[\Upsilon ]= & {} I_{\Upsilon }(\Upsilon _{1}, \Upsilon _{2},\ldots , \Upsilon _{m}),\\ I_{\Theta }[\Theta ]= & {} I_{\Theta } (\Theta _{1}, \Theta _{2},\ldots , \Theta _{m}), \end{aligned}$$

where $m$ is the number of all nodes in all elements. Then using the Rayleigh-Ritz procedure to minimize $I_{W}[W], I_{\Upsilon }[\Upsilon ], I_{\Theta }[\Theta ]$, with respect to the variational parameters $W_{i}, \Upsilon _{i}, \Theta _{i}$ respectively, gives

$$\begin{aligned}&\frac{\partial I_{W}}{\partial W_{i}}=\sum _{e=1}^{n}\frac{\partial I_{W}^{e}}{\partial W_{i}}=0,\ \,\,\,\,\,\,\,\,\,\,i=1,\ldots ,m, \end{aligned}$$

(81)

$$\begin{aligned}&\frac{\partial I_{\Upsilon }}{\partial \Upsilon _{i}}=\sum _{e=1}^{n}\frac{\partial I_{\Upsilon }^{e}}{\partial \Upsilon _{i}}=0,\ \,\,\,\,\, i=1,\ldots ,m, \end{aligned}$$

(82)

$$\begin{aligned}&\frac{\partial I_{\Theta }}{\partial \Theta _{i}}=\sum _{e=1}^{n}\frac{\partial I_{\Theta }^{e}}{\partial \Theta _{i}}=0,\ \,\,\,\,\, i=1,\ldots ,m. \end{aligned}$$

(83)

Hence, we arrive to following formula

$$\begin{aligned} \frac{\partial I_{W}^{e}}{\partial W_{n_{i}}}= & {} -2\int _{e}\frac{d N_{n_{i}}^{p}}{d z} \frac{d N^{e}}{d z} \delta _{W}^{e}dz-2\int _{e} N_{n_{i}}^{p} N^{e} \delta _{\Upsilon }^{e}dz=0, \end{aligned}$$

(84)

$$\begin{aligned} \frac{\partial I_{\Upsilon }^{e}}{\partial \Upsilon _{n_{i}}}= & {} 2a^{4}\int _{e} N_{n_{i}}^{p} N^{e} \delta _{W}^{e}dz-2\int _{e}\left( \frac{d N_{n_{i}}^{p}}{d z} \frac{d N^{e}}{d z} +2a^{2} N_{n_{i}}^{p} N^{e} \right) \delta _{\Upsilon }^{e}dz\nonumber \\&-2a^{2}R\int _{e}N_{n_{i}}^{p} N^{e} \delta _{\Theta }^{e}dz=0, \end{aligned}$$

(85)

$$\begin{aligned} \frac{\partial I_{\Theta }^{e}}{\partial \Theta _{n_{i}}}= & {} -2\int _{e}\left[ \frac{d N_{n_{i}}^{p}}{d z} \frac{d N^{e}}{d z} +a^{2} N_{n_{i}}^{p} N^{e} \right] \delta _{\Theta }^{e}dz +2R \int _{e} N_{n_{i}}^{p} N^{e} \delta _{W}^{e}dz=0. \end{aligned}$$

(86)

Then, the matrix representation for the system of equation of element $e$ take the form

$$\begin{aligned} \left( \begin{array}{ccc} -D_{2}^{e} &{} -F_{1}^{e} &{} O \\ a^{4}F_{1}^{e} &{} -D_{2}^{e}-2a^{2}F_{1}^{e} &{} O \\ O &{} O &{} -D_{2}^{e}-a^{2}F_{1}^{e} \\ \end{array} \right) \left( \begin{array}{c} \delta _{W}^{e} \\ \delta _{\Upsilon }^{e} \\ \delta _{\Theta }^{e} \\ \end{array} \right) =R \left( \begin{array}{cccc} O &{} O &{} O \\ O &{} O &{} a^{2}\,F_{1}^{e} \\ -F_{1}^{e} &{} O &{} O \\ \end{array} \right) \left( \begin{array}{c} \delta _{W}^{e} \\ \delta _{\Upsilon }^{e} \\ \delta _{\Theta }^{e} \\ \end{array} \right) ,\quad \quad \quad \end{aligned}$$

(87)

where

$$\begin{aligned} O_{ij}= & {} 0,\,\,\,D_{2\,ij}^{e}=\int _{-1}^{1}\frac{d N_{n_{i}}^{p}}{d z} \frac{d N_{n_{j}}^{p}}{dz} dz, \\ F_{1\,ij}^{e}= & {} \int _{-1}^{1}N_{n_{i}}^{p} N_{n_{j}}^{p} dz, \,\,\,\,\,\, i=1,\ldots ,p+1,\,j=1,\ldots ,p+1. \end{aligned}$$

The above integral can be evaluated by applying the classical finite element nodal basis functions in one dimension on the standard element $\Omega _{st}=(-1,1)$. The standard shape functions are defined by the set of Lagrange polynomials

$$\begin{aligned} N_{n_{i}}^{p}(\zeta )=\prod _{j=n_{1},\, j\ne i}^{n_{p+1}} \frac{\zeta -\zeta _{j}}{\zeta _{i}-\zeta _{j}}, \end{aligned}$$

(88)

where $\zeta =2/h(z-z_{m})$, h is the length of element and $z_{m}$ is the mid-point. Thus these integrations take the form:

$$\begin{aligned} D_{2\,ij}^{e}= & {} \frac{2}{h}\int _{-1}^{1}\frac{d N_{n_{i}}^{p}}{d \zeta } \frac{d N_{n_{j}}^{p}}{d \zeta } d\zeta ,\\ F_{1\,ij}^{e}= & {} \frac{h}{2}\int _{-1}^{1}N_{n_{i}}^{p} N_{n_{j}}^{p} d\zeta \\&i=1,\ldots ,p+1,\,\,\,\,\,\,\,\,\,\,\,\,\,\,j=1,\ldots ,p+1. \end{aligned}$$

The integral were calculated analytically using Matlab routines. Finally, we assemble the systems of all elements $e=1,\ldots ,n$ to get the main system which has the form

$$\begin{aligned} \Lambda \,\Sigma =R \,\Xi \,\Sigma , \end{aligned}$$

(89)

where $\Sigma =(W_{1},\ldots , W_{m}, \Upsilon _{1},\ldots , \Upsilon _{m}, \Theta _{1}\ldots , \Theta _{m})^{T}$.

For free surface boundary conditions, we imposing the boundary conditions easily by removing $W_{1}, W_{m}, \Upsilon _{1}, \Upsilon _{m}, \Theta _{1}, \Theta _{m}$ from the system and thus we remove the rows and columns of order $1, m, m+1, 2m, 2m+1$ and $3m$.

For slip-slip boundary conditions, firstly, we transform these condition to the local coordinate $\zeta$ and thus these boundary condition take the form

$$\begin{aligned} \begin{array}{l@{\quad }l} W=\Theta =0,&{} \mathrm {at}\quad \zeta =-1,1,\\ \lambda _{L}\Upsilon -\frac{2}{h}DW=0,&{} \mathrm {at}\quad \zeta =-1,\\ \lambda _{U}\Upsilon +\frac{2}{h}DW=0,&{} \mathrm {at}\quad \zeta =1. \end{array} \end{aligned}$$

(90)

Now we can impose the boundary conditions (90)$_{1}$ by removing the rows and columns of order $1, m, 3m+1$ and $3m$, while we can apply the boundary conditions (90)$_{2,3}$ by substituting $\zeta =-1$ in (90)$_{2}$ and $\zeta =1$ in (90)$_{3}$ and then use (69) to arrive to the following final equations

$$\begin{aligned} \Upsilon _{1}= & {} \frac{2}{\lambda _{L}he_{1}}\left[ d_{2}W_{2}+\cdots +d_{p+1}W_{p+1}]-\frac{1}{e_{1}}[e_{2}\Upsilon _{2}+\cdots +e_{p+1}\Upsilon _{p+1}\right] ,\quad \end{aligned}$$

(91)

$$\begin{aligned} \Upsilon _{m}= & {} -\frac{2}{\lambda _{U}hg_{p+1}}\left[ f_{1}W_{m-p}+f_{2}W_{m-p+1}+\cdots +f_{p}W_{m-1}\right] \nonumber \\&-\frac{1}{g_{p+1}}\left[ g_{1}\Upsilon _{m-p}+\cdots +g_{p}\Upsilon _{m-1}\right] , \end{aligned}$$

(92)

where $d_{i}=dN_{i}(-1)/d\zeta$,  $e_{i}=N_{i}(-1)$, $g_{i}=N_{i}(1)$ and $f_{i}=dN_{i}(1)/d\zeta$.

Now, we substitute the values of $A_{1}$ and $A_{m}$ in (89) thus the columns of order $m+1$ and $2m$ will be zeros and thus we can remove the $m+1$ and $2m$ rows and columns.

However, for fixed-slip boundary conditions, we change the conditions $D W=0$ at $z=0,$ to another conditions related with function $\Upsilon$. To do this, let the first element is $[0, a]$, then, we integrate (69)$_{1}$ for first element and use the boundary condition $W=0$ at $z=0$, to arrive to the following conditions

$$\begin{aligned} DW^{1}(a)=\int _{0}^{a}\Upsilon ^{1}dz, \end{aligned}$$

(93)

where $1$ refer to the first element. Next, to make (93) are useful in our computational, we transform this condition to the local coordinate $\zeta$ and use (73) to obtain

$$\begin{aligned} \frac{2}{h}\frac{dN(1)}{d\zeta }\delta _{W}^{1}=\frac{h}{2}\int _{-1}^{1} N(\zeta )\delta _{A}^{1}d\zeta , \end{aligned}$$

(94)

Now, let $\tau _{i}=dN_{i}(1)/d\zeta$ and $\xi _{i}=\int _{-1}^{1}N_{i}(\zeta )d\zeta$. In addition, suppose that the nodes for the first element and for the last element have the order $1,2,\ldots ,p+1$ and $m-p,m-p+1,\ldots ,m$, respectively. Then, (94) leads to the following computational conditions

$$\begin{aligned} \Upsilon _{1}=\frac{4}{h^{2}\xi _{1}}[\tau _{2}W_{2}+\cdots +\tau _{p+1}W_{p+1}]-\frac{1}{\xi _{1}}[\xi _{2}\Upsilon _{2}+\cdots +\xi _{p+1}\Upsilon _{p+1}]. \end{aligned}$$

(95)

Now, we substitute the value of $\Upsilon _{1}$ from (95) and the value of $\Upsilon _{m}$ from (92) in (89) thus the columns of order $m+1$ and $2m$ will be zeros. Now, we can remove the rows and columns of order $1, m, m+1, 2m, 2m+1$ and $3m$.

Chebyshev Collocation Methods

In this section, we use the Chebyshev collocation method to solve the eigenvalue system (27) and (28). We apply two techniques of Chebyshev collocation method to impose the boundary conditions.

Method 1

For free-free boundary conditions, we use the same transformation which have been used in the finite element method, to arrive

$$\begin{aligned} \begin{array}{l} D^{2}W=\Upsilon ,\\ D^{2}\Upsilon -2a^{2}\Upsilon +a^{4}W-a^{2} R \Theta =0,\\ (D^{2}-a^{2})\Theta + R W=0, \end{array} \end{aligned}$$

(96)

with the same boundary conditions (70). In the first technique, we choose the trial functions such that these functions satisfy the boundary conditions as follows:

$$\begin{aligned} W= & {} \sum _{n=1}^{N} W_{n} \theta _{n}(z), \end{aligned}$$

(97)

$$\begin{aligned} \Upsilon= & {} \sum _{n=1}^{N} \Upsilon _{n} \theta _{n}(z), \end{aligned}$$

(98)

$$\begin{aligned} \Theta= & {} \sum _{n=1}^{N} \Theta _{n} \theta _{n}(z), \end{aligned}$$

(99)

where

$$\begin{aligned} \theta _{n}(z)=(1-z^2)T_{2n-2}(z). \end{aligned}$$

(100)

Here $T_{n}(z)$ is the $n$th-degree of Chebyshev polynomial of the first kind, which is defined by

$$\begin{aligned}&T_{0}(z)=1,\,\,\,\,\,\,\,T_{1}(z)=z,\\&\quad T_{n+1}(z)-2zT_{n}(z)+T_{n-1}(z)=0,\,\,\,\,\,-1 \le z \le 1, \end{aligned}$$

or

$$\begin{aligned} T_{n}(z)=\cos (n \arccos (z)),\,\,\,\,\,-1 \le z \le 1, \end{aligned}$$

It is clear that $w, \Phi$ and $\Theta$ satisfies the boundary conditions (70). Now, Substituting (97)–(99) into (96), and requiring that (96) be satisfied at $N$ collocation points $z_{1}, z_{2},\ldots ,z_{N}$, where

$$\begin{aligned} z_{i}=\cos \left( \frac{i-1}{2N-1}\pi \right) ,\quad i=1,\ldots ,N, \end{aligned}$$

(101)

we obtain $2N$ algebraic equations for $2N$ unknowns $w_{1},\ldots , w_{N}, A_{1},\ldots , A_{N}, \Upsilon _{1},\ldots , \Upsilon _{N}$:

$$\begin{aligned} \left( \begin{array}{ccc} 4D^{2} &{} -I &{} O\\ a^{4} I &{} 4D^{2}-2a^{2}I &{} O \\ O &{} O &{} 4D^{2}-a^{2}I\\ \end{array} \right) X =R \left( \begin{array}{ccc} O &{} O &{} O\\ O &{} O &{} a^{2} R I \\ - R I &{} O &{} O\\ \end{array} \right) X, \end{aligned}$$

where $X=(w_{1},\ldots , w_{N}, \Upsilon _{1},\ldots , \Upsilon _{N}, \Theta _{1},\ldots , \Theta _{N}), O$ is the zeros matrix, $I(n_{1},n_{2})= \theta _{n_{2}}(z_{n_{1}}),\,\,D^{2}(n_{1},n_{2})=\theta {^{\prime \prime }}_{n_{2}}(z_{n_{1}}),\,\,n_{1}=1,\ldots ,N,\,\,n_{2}=1,\ldots ,N$. The above differentiation matrices, which are corresponded to the trail functions (97)–(100), were computed analytically using Matlab routines.

Now, to deal with the solution of (27) and (28) with slip-slip boundary conditions and $\lambda _{L}=\lambda _{U}=\lambda$. Firstly, we introduce a new function $\vartheta =Dw$, then our system can written in the following form:

$$\begin{aligned}&Dw-\vartheta =0, \end{aligned}$$

(102)

$$\begin{aligned}&D^{3}\vartheta -2a^{2}D \vartheta +a^{4}W-a^{2} R \, \Theta =0, \end{aligned}$$

(103)

$$\begin{aligned}&(D^{2}-a^{2})\Theta + R W=0. \end{aligned}$$

(104)

Then, according to this transform, the boundary condition have the form

$$\begin{aligned} W=\Theta= & {} 0,\,\,\,\,\,\,\,\mathrm {at}\,\,\,\,\,\,\, z=0, 1, \end{aligned}$$

(105)

$$\begin{aligned} \lambda D\vartheta +\vartheta= & {} 0,\,\,\,\,\,z=1, \end{aligned}$$

(106)

$$\begin{aligned} \lambda D \vartheta -\vartheta= & {} 0,\,\,\,\,\,z=0. \end{aligned}$$

(107)

We choose the trial functions such that these functions satisfy the boundary conditions as follows:

$$\begin{aligned} W= & {} \sum _{n=1}^{N} W_{n} \theta _{n}(z), \end{aligned}$$

(108)

$$\begin{aligned} \vartheta= & {} \sum _{n=1}^{N} \vartheta _{n} \psi _{n}(z), \end{aligned}$$

(109)

where

$$\begin{aligned} \theta _{n}(z)= & {} (1-z^{2})T_{2n-2}(z), \end{aligned}$$

(110)

$$\begin{aligned} \psi _{n}(z)= & {} \left( 1-z^2+\frac{2\lambda }{1+\lambda (2n-1)^{2}}\right) T_{2n-1}(z). \end{aligned}$$

(111)

It is clear that $W$ and $\vartheta$ satisfy the boundary conditions (105)–(107). Now we can apply the same procedure which are used for free-free boundary conditions. For fixed-slip boundary conditions, we are unable to solve the problem because it is not easy to suggest a functions which satisfy two different boundary conditions. However, if $\lambda _{L} \ne \lambda _{U}$, it is not easy to find a trial functions such that these functions satisfy the boundary conditions.

Method 2

For free-free boundary conditions we adopt the modified system (96). We expand the solutions of the governing stability Eq. (96) as truncated series of Chebyshev polynomials

$$\begin{aligned} W=\sum _{n=0}^{N} w_{n} T_{n}(z),\,\,\,\,\,\,\,\,\,\Upsilon =\sum _{n=0}^{N} \Upsilon _{n} T_{n}(z),\,\,\,\,\,\,\,\,\,\Theta =\sum _{n=0}^{N} \Theta _{n} T_{n}(z), \end{aligned}$$

(112)

then, we insert (112) into the Eq. (96), and then substitute the Gauss–Labatto points which are defined by

$$\begin{aligned} y_{i}=\cos \left( \frac{\pi i}{N-3}\right) ,\,\,\,\,\,\,\, i=0,\ldots ,N-2. \end{aligned}$$

(113)

Thus, we obtain $3N-3$ algebraic equations for $3N+3$ unknowns $W_{0},\ldots , W_{N}, \Upsilon _{0},\ldots , \Upsilon _{N}, \Theta _{0},\ldots , \Theta _{N}$. Now, we can add six rows using the boundary conditions (70) as follows

$$\begin{aligned} BC_{1}:\,\,\,\,\, \sum _{n=0}^{N} W_{n}= & {} 0, \end{aligned}$$

(114)

$$\begin{aligned} BC_{2}:\,\,\,\,\,\sum _{n=0}^{N} (-1)^n W_{n}= & {} 0, \end{aligned}$$

(115)

$$\begin{aligned} BC_{3}:\,\,\,\,\, \sum _{n=0}^{N} \Upsilon _{n}= & {} 0, \end{aligned}$$

(116)

$$\begin{aligned} BC_{4}:\,\,\,\,\,\sum _{n=0}^{N} (-1)^n \Upsilon _{n}= & {} 0, \end{aligned}$$

(117)

$$\begin{aligned} BC_{5}:\,\,\,\,\, \sum _{n=0}^{N} \Theta _{n}= & {} 0, \end{aligned}$$

(118)

$$\begin{aligned} BC_{6}:\,\,\,\,\,\sum _{n=0}^{N} (-1)^n \Theta _{n}= & {} 0. \end{aligned}$$

(119)

Then, we obtain the generalised eigenvalue problem

$$\begin{aligned} \left( \begin{array}{ccc} 4D^{2} &{} -I &{} O\\ BC_{1} &{} 0\ldots 0 &{} 0\ldots 0\\ BC_{2} &{} 0\ldots 0&{} 0\ldots 0\\ a^{4} I &{} 4D^{2}-2a^{2}I &{} O \\ 0\ldots 0 &{} BC_{3}&{} 0\ldots 0\\ 0\ldots 0 &{} BC_{4}&{} 0\ldots 0\\ O &{} O &{} 4D^{2}-a^{2}I\\ 0\ldots 0 &{} 0\ldots 0&{} BC_{5}\\ 0\ldots 0 &{} 0\ldots 0&{} BC_{6}\\ \end{array} \right) X =R \left( \begin{array}{ccc} O &{} O &{} O\\ 0\ldots 0 &{} 0\ldots 0&{} 0\ldots 0\\ 0\ldots 0 &{} 0\ldots 0&{} 0\ldots 0\\ O &{} O &{} a^{2} R I \\ 0\ldots 0 &{} 0\ldots 0&{} 0\ldots 0\\ 0\ldots 0 &{} 0\ldots 0&{} 0\ldots 0\\ - R I &{} O &{} O\\ 0\ldots 0 &{} 0\ldots 0&{} 0\ldots 0\\ 0\ldots 0 &{} 0\ldots 0&{} 0\ldots 0\\ \end{array} \right) X, \end{aligned}$$

where $X=(W_{0},\ldots , W_{N}, \Upsilon _{0},\ldots , \Upsilon _{N}, \Theta _{0},\ldots , \Theta _{N}), O$ is the zeros matrix, $I(n_{1},n_{2})= T_{n_{2}}(z_{n_{1}}),\,\, D^{2}(n_{1},n_{2})=T''_{n_{2}}(z_{n_{1}}),\,\,n_{1}=0,\ldots ,N-2,\,\,n_{2}=0,\ldots ,N$. Again, we computed the differentiation matrices, which are corresponded to the trail functions (112) analytically using Matlab routines.

However, for numerical solutions of (27) and (28) with slip-slip boundary conditions, we adopt system (102)–(104) with boundary conditions (105)–(107) but with general values of $\lambda _{U}, \lambda _{L}$ i.e. $\lambda _{U} \ne \lambda _{L}$. Then, we apply the same procedure which are used for free-free boundary conditions to arrive to the following eigenvalue problem system

$$\begin{aligned} \left( \begin{array}{ccc} 2D &{} -I &{} O\\ BC_{1} &{} 0\ldots 0 &{} 0\ldots 0\\ BC_{2} &{} 0\ldots 0&{} 0\ldots 0\\ a^{4} I &{} 8D^{3}-4a^{2}D &{} O \\ 0\ldots 0 &{} BC_{3}&{} 0\ldots 0\\ 0\ldots 0 &{} BC_{4}&{} 0\ldots 0\\ O &{} O &{} 4D^{2}-a^{2}I\\ 0\ldots 0 &{} 0\ldots 0&{} BC_{5}\\ 0\ldots 0 &{} 0\ldots 0&{} BC_{6}\\ \end{array} \right) X =R \left( \begin{array}{ccc} O &{} O &{} O\\ 0\ldots 0 &{} 0\ldots 0&{} 0\ldots 0\\ 0\ldots 0 &{} 0\ldots 0&{} 0\ldots 0\\ O &{} O &{} a^{2} R I \\ 0\ldots 0 &{} 0\ldots 0&{} 0\ldots 0\\ 0\ldots 0 &{} 0\ldots 0&{} 0\ldots 0\\ - R I &{} O &{} O\\ 0\ldots 0 &{} 0\ldots 0&{} 0\ldots 0\\ 0\ldots 0 &{} 0\ldots 0&{} 0\ldots 0\\ \end{array} \right) X, \end{aligned}$$

where $X=(W_{0},\ldots , W_{N}, \vartheta _{0},\ldots , \vartheta _{N}, \Theta _{0},\ldots , \Theta _{N}), O$ is the zeros matrix, $I(n_{1},n_{2})= T_{n_{2}}(z_{n_{1}}),\,\,D(n_{1},n_{2})=T_{n_{2}}(z_{n_{1}}),\,\, D^{3}(n_{1},n_{2})=T{^{\prime \prime \prime }|}_{n_{2}}(z_{n_{1}}) ,\,\,n_{1}=0,\ldots ,N-2,\,\,n_{2}=0,\ldots ,N$. The boundary conditions $BC1, BC2, BC5, BC6$ have the same form of the free-free boundary conditions, while $BC3$ and $BC4$ have different form which is related with slip boundary conditions and these written in the following form:

$$\begin{aligned} BC_{3}:\,\,\,\,\,\sum _{n=0}^{N} (2\lambda _{U}n^{2}+1) \vartheta _{n}= & {} 0, \end{aligned}$$

(120)

$$\begin{aligned} BC_{4}:\,\,\,\,\,\sum _{n=0}^{N} (2(-1)^{n+1}\lambda _{L}n^{2}-1) \vartheta _{n}= & {} 0. \end{aligned}$$

(121)

For fixed-slip boundary conditions, the system have the same form of slip-slip with $\lambda _{L}=0$. This show the flexibility of this method, and thus, we believe that this method is more flexible than the other numerical methods. Moreover, this method is very accurate and not need to difficult numerical treatments.

Remark 1

The finite difference, High order finite difference and order finite element methods lead to non-singular matrices in the eigenvalue systems and so we may employ $LU$ decomposition, unlike the other numerical methods which necessarily have singular matrices in the eigenvalue systems and so necessitate use of the $QZ$ algorithm. The banded nature of matrices in the eigenvalue systems would also lead naturally to solution by an Arnoldi technique. In addition, we found no occurrence of spurious eigenvalues as frequently arise with the Chebyshev methods, (cf. [4, 8]).

Numerical Results of Thermal Convection and Conclusions

In this section we report our numerical results of (27) and (28) with free-free, slip-slip and fixed-slip boundary conditions. It is very important to make a comparison between the accuracy of the methods to the crucial measurement being the exact solution. As we mention in the introduction [36], showed that, in the case of free-free boundary conditions, we may obtain the analytical result of the Rayleigh number $Ra_{crit} = 27 \pi ^{4}/4$ and $a^{2}_{crit} = \pi ^{2}/2$. Thus, we compute the absolute error of the critical Rayleigh and the wave numbers from the relation:

$$\begin{aligned} e_{a}=|a_{exact}-a_{numer}|,\,\,\,\,\,\,\,\,\,e_{R}=|R_{exact}-R_{numer}| \end{aligned}$$

where $e_{a}, e_{R}, a_{exact}, a_{numer}, R_{exact}$ and $R_{numer}$ are the absolute error of wavenumber, the absolute error of critical Rayleigh number, the exact value of wavenumber, the numerical value of wavenumber, the exact value of the critical Rayleigh number and the numerical value of the critical Rayleigh number, respectively. We report the absolute error of the critical Rayleigh and the wave numbers for the numerical methods with free-free boundary conditions in Figs. 1, 2, 3, 4, 5 and 6.

Fig. 1

The absolute error of wavenumbers of finite element method plotted against the number of elements

Fig. 2

The absolute error of critical Rayleigh numbers of FE method plotted against the number of elements

Fig. 3

The absolute error of wavenumbers of Chebyshev tau (solid line), Chebyshev collection method-1 (dashed line) and Chebyshev collection method-2 (dotted line) plotted against the number of Chebyshev polynomials

Fig. 4

The absolute error of critical Rayleigh numbers of Chebyshev tau (solid line), Chebyshev collection method-1 (dashed line) and Chebyshev collection method-2 (dotted line) plotted against the number of Chebyshev polynomials

Fig. 5

The absolute error of wavenumbers of finite difference method (solid line) and high order finite difference method (dashed line) plotted against $h$

Fig. 6

The absolute error of critical Rayleigh numbers of finite difference method (solid line) and high order finite difference method (dashed line) plotted against $h$

In Fig. 1, the absolute errors of the wavenumber which are produce from the finite element method are introduced. The table demonstrates that the absolute error of the wavenumbers does not change when increasing the number of nodes and number elements, where, in general, the values of absolute error $9.15 \times 10^{-6}$ or very close to this value can be seen. When the number of elements are 2, 3, 4 and the number of nodes is 4, the accuracy is less than normal i.e. the accuracy is less than $9.15 \times 10^{-6}$ which is very much to be expected as the approximation in these situations is of a low order. However, with a high number of elements and nodes, the accuracy values oscillate where the absolute error is less or greater than the normal value. Furthermore, this behaviour is to be expected in studying hydrodynamic stability problems. As the number of nodes and element increase, theoretically, the accuracy of the finite element method should increase. However, as the number of nodes and elements increase, the computer calculations increase and thus computer error will be greater. Since hydrodynamic stability problems require solving the eigenvalue system many times in order to achieve the critical Rayleigh number which corresponds to the critical wavenumber, computer error is expected to make the absolute error greater than the theoretical one. This behaviour is very clear in Fig. 6, where the absolute error is higher than the truncation error of FD and HFD methods. For both the absolute error of the critical Rayleigh and wave numbers, the accuracy increases with a higher number of nodes and elements until it reaches a peak, at which point the behaviour of the accuracy oscillates.

Similar behaviour of the absolute error of the wavenumber can be found in the Chebyshev Tau, Chebyshev collection 1 and 2 schemes (see Figs. 3, 5). However, the accuracy of the wave number for standard and high order finite difference methods is less than in the Chebyshev methods and finite element methods. Generally, the accuracy of the finite difference schemes corresponds with the value of $h$, where accuracy increases as the value of $h$ decreases. However, the value of $h$ cannot be taken as less than $0.01$ for two reasons. First, with $h=0.01$ and after imposing the boundary conditions, the order of the eigenvalue matrices will be $199 \times 199$, and thus, according to the computer’s ability, this is an optimal choice. Secondly, it is difficult to believe that low $h$ values can provide more accurate results as the computations have to increase rapidly, which is very important spatially when the nonlinear stability problems are solved.

Figure 2 presents the absolute errors of the Rayleigh number generated from the finite element method. Unlike the behaviour of the wavenumber, whereby the required accuracy with a low number of elements and nodes can be achieved, the accuracy of the Rayleigh numbers require a greater number of elements and nodes. Additionally, the absolute errors of the Rayleigh number are less than the absolute errors of the wavenumber. Similar to the wave numbers accuracy, the accuracy of the Rayleigh number increases with an increased number of nodes and elements up to the normal absolute error which is $2 \times 10^{-10}$. At this point the accuracy behaviour oscillates more. Figure 2 shows that the required accuracy can be achieved with $2$ elements and $8$ nodes and this choice is the best according to the computer run time. In other words, using this choice the best accuracy with the smallest eigenvalue matrices can be achieved. This means that the order of eigenvalue matrices will be $30 \times 30$. Now, from the results of the other numerical methods in Figs. 3 and 4, it can be seen the following observations

• The Chebyshev Tau method can achieve the required accuracy using at least $11$ Chebyshev polynomials, thus the order of eigenvalue matrices will be $26 \times 26$ after adding the two rows of boundary conditions.

• The Chebyshev collection method 1 can achieve the required accuracy using $6$ Chebyshev polynomials. Thus the order of the eigenvalue matrices will be $16 \times 16$ after adding the two rows of boundary conditions. Here we should mention that the Chebyshev collection method 1 achieves the required accuracy with a smaller number of Chebyshev polynomials because it uses only the even polynomials, i.e. this method uses the polynomials $T_{0}, T_{2},\ldots$. So when the required accuracy with $6$ Chebyshev polynomials is reached, this signifies that it uses the polynomials of order $0, 2, 4, 6, 8, 10$. This is a very significant advantage of this method, achieving a high level of accuracy with low numbers of polynomials.

• The Chebyshev collection method 2 can achieve the required accuracy using $10$ Chebyshev polynomials. Thus the order of eigenvalue matrices will be $24 \times 24$ after adding the two rows of boundary conditions.

To place the efficiency of these numerical methods in context Figs. 7, 8 and 9 provide a visual representation of the computational time required to converge to the first eigenvalue of the spectrum as the number of polynomials is increased. Figure 7, shows that the computational time of the finite element method increases with increasing the number of nodes and elements. From Fig. 8, in a direct comparison between both the Chebyshev methods, the Chebyshev collection method-1 variant was found to be marginally more computationally efficient. The difference is to be expected, as the Chebyshev collection method-1 has the highest accuracy with the same number of polynomials, although the results are sufficiently close to make any visual representation unnecessary. Figure 9 clearly demonstrate that as the finite difference method is substantially more computationally efficient than the high order finite difference method.

Fig. 7

The computational time (seconds) of finite element method plotted against the number of elements

Fig. 8

The computational time (seconds) of Chebyshev tau (solid line), Chebyshev collection method-1 (dashed line) and Chebyshev collection method-2 (dotted line) plotted against the number of Chebyshev polynomials

Fig. 9

The computational time (seconds) of finite difference method (solid line) and high order finite difference method (dashed line) plotted against $h$

Here, it is very important to mention that these numerical methods reach a high level of accuracy with these numbers of polynomials, nodes and elements for free-free boundary conditions. However, there is a need to increase these numbers when dealing with other kinds of boundary conditions. Accordingly, the results of numerical methods are reported in Tables 1 and 2 for slip-slip and fixed-slip boundary conditions with selected numbers of polynomials, nodes and elements and $h$. The convergence has been checked and with the result that $8$ decimal places for slip-slip and fixed-slip boundary conditions is achieved using $20$, $10$ and $20$ Chebyshev polynomials for Chebyshev Tau, Chebyshev collection method 1 and Chebyshev collection method 2, respectively, while for the finite element method the convergence to $8$ decimal places required at least three elements and ten nodes. Additionally, the required convergence was satisfied with $h=0.001$, $h=0.01$ for FD and HFD, respectively. Therefore, the numerical results in Tables 1 and 2 are reported according to these choices.

Table 1 The critical Rayleigh and wave numbers for symmetric-slip case for a selection of $\lambda$ values

Table 2 The critical Rayleigh for fixed-slip case for a selection of $\lambda _{U}$ values

One of the important reasons for applying different numerical methods for solving the system (27) and (28) is to make a comparison between these methods so that a conclusion can be reached regarding which of these methods is the best in solving hydrodynamic stability problems. The advantage in using the Chebyshev Tau method is that it can achieve the required accuracy using a small number of polynomials, allowing the achievement of highly accurate results with a short run time. On the other hand, the Chebyshev Tau method is not easy to apply, as a considerable effort to solve any system of equations is required. With regard to problems of variable coefficients, the Chebyshev Tau method presents complications in finding the solution because this method depends on the writing of all functions in the system of the equation in the form of Chebyshev polynomial series, which proves difficult when functions such as triangular and hyperbolic functions are present.

The finite difference method needs a large number of divisions to reach the required accuracy, while the high order finite difference method can reach to the desired accuracy by using less number of divisions. We see from the results that the finite difference method requires $h=0.001$ to achieve a very good accuracy and convergence results, while we use $h=0.01$ for High order finite difference method. Now, let we define the rate of convergence between two norms, $e_{1}$ and $e_{2}$, with two different grid spacings $h_{1}$ and $h_{2}$, with $h_{1} < h_{2}$, as

$$\begin{aligned} P=\log (e_{1}/e_{2})/\log (h_{1}/h_{2}) \end{aligned}$$

For standard and high order finite difference methods, we computed the rate of convergence with $h_{1}=0.01$ and $h_{2}=0.02$. The value of the rates of convergence are $1.93$ and $2.7$ for standard and high order finite difference methods, respectively, which indicates that the standard finite difference has less than second order of accuracy and the high order finite difference has less than third order of accuracy.

We believe that the finite element method is one of the best choices in solving the hydrodynamic stability problems as it is a flexible method and it can give very accurate results. However, the disadvantage of this method is that the required run time is longer than the Chebyshev numerical methods. Thus, using this method for linear theory problem is recommended. Moreover, it could be a suitable choice for nonlinear problems but with no more than two maximised parameters. The Chebyshev collection method 1 is an inflexible method where if the boundary conditions are not symmetric it becomes impossible to find trial functions satisfying the boundary conditions. However, this method has the highest accuracy between the numerical methods and requires a smaller number of polynomials to achieve excellent accuracy and convergence. This point is of great value because our numerical calculations require great computational time, with a run time in some cases of more than ten hours.

Finally, the Chebyshev collection method 2 is the most flexible method among the other numerical methods and can achieve high accuracy using a reasonable number of Chebyshev polynomials. Usually, we used this method in my work because of the highly accurate results within a short time. The author strongly recommends this method for use in solving hydrodynamic stability problems.

Instability in Poiseuille Flow in a Porous Medium with no Slip Boundary Conditions

The equations for Poiseuille flow in a Brinkman porous material are derived in [33] and in [39], p. 234. With $v_i$ being the velocity field, $p$ the pressure, $R$ being the Reynolds number, and $M^2$ a non-dimensional Darcy (friction) coefficient the governing equations are

$$\begin{aligned} \begin{aligned}&R\biggl (\frac{\partial \mathbf {v}}{\partial t}+ (\mathbf {v} \cdot \nabla ) \mathbf {v} \biggr )=- \nabla p + \nabla ^{2} \mathbf {v} -M^{2} \mathbf {v},\\&\nabla \cdot \mathbf {v}=0, \end{aligned} \end{aligned}$$

(122)

where standard notation is employed, and Eq. (122) hold in the domain $\{(x,y)\in \mathbb {R}^2\}\times \{z\in (-1,1)\}\times \{t>0\}$.

The basic solution whose stability we are interested in is one where the fluid is driven along the channel in the $x$-direction by a constant pressure gradient of form

$$\begin{aligned} -\frac{\partial {\bar{p}}}{\partial x}=K^2>0, \end{aligned}$$

where an overbar denotes the basic state and $K^2$ is the velocity at the center line of the channel. The basic velocity field corresponding to this pressure gradient has form ${\bar{\mathbf{v}}}=(U(z),0,0)$. Then one finds $U$ must satisfy the boundary value problem

$$\begin{aligned} -K^2=U^{\prime \prime }-M^2U,\qquad -1<z<1, \end{aligned}$$

(123)

subject to boundary conditions

$$\begin{aligned} U=0,\qquad z=\pm 1. \end{aligned}$$

(124)

If in our non-dimensionalization we pick $K^2$ so that $K^2=M^2\,\cosh M/(\cosh M-1)$ then $U$ is found to be

$$\begin{aligned} U=\frac{\cosh M-\cosh Mz}{\cosh M-1}. \end{aligned}$$

(125)

We now wish to investigate the stability of solution (125) and so let $\mathbf{u}=(u,v,w)$ be a perturbation to ${\bar{\mathbf{v}}}$ with corresponding pressure perturbation $\pi$. The perturbation equations are derived in detail in [39], p. 235, and he shows that after linearization and assuming spatial and time dependence like $\exp (i\alpha x+i\beta y-ict)$ then one may show that $w(z)$ satisfies the equation

$$\begin{aligned} (D^2-a^2)^2w-M^2(D^2-a^2)w=iaR(U-c)(D^2-a^2)w-iaRU^{\prime \prime }w, \end{aligned}$$

(126)

where $D=d/dz, a^2=\alpha ^2+\beta ^2,$ and $z\in (-1,1)$. The boundary conditions which $w$ must satisfy are

$$\begin{aligned} w=Dw=0,\qquad z=\pm 1. \end{aligned}$$

(127)

In the next section we briefly describe the numerical methods employed to solve Eq. (126) together with (127).

Numerical Techniques

To the best of our knowledge, there are three main texts which deal with the numerical solution of the eigenvalue systems of the equations of the Poiseuille flow stability problem. The first article was introduced by [33] who developed a linearised instability analysis. Then, [39], section 5.8.1 shows that the analysis of [33] is incomplete in that he omits a term from one equation. [24] presented a numerical solution of the corrected system of equations and they stress that the fundamental model is due to [33]. Here we apply three numerical methods and discuss the advantage and disadvantage of each one.

Chebyshev Collocation Methods

Since solving (126) and (127) is a difficult numerical problem, we adopt the Chebyshev collocation method to solve the eigenvalue systems. We apply two variations of the Chebyshev collocation method in order to incorporate the boundary conditions, and also two methods provide an independent check. Firstly we introduce the function $\varphi =Dw$, and then Eq. (126) may be written as the system

$$\begin{aligned} \begin{aligned}&Dw-\varphi =0,\\&D^3 \varphi -(2a^2+M^2)D \varphi +(a^4+a^2M^2)w+iaRU^{\prime \prime }w\\&\qquad =iaR(U-c)(D \varphi -a^2w). \end{aligned} \end{aligned}$$

(128)

The boundary conditions (127) now become

$$\begin{aligned} w=\varphi =0,\quad z=\pm 1, \end{aligned}$$

(129)

Method 1

Here we expand $w$ and $\varphi$ as (truncated) series in trial function $\theta _n$, so that

$$\begin{aligned} w=\sum _{n=1}^Nw_n\theta _n(z),\qquad \varphi =\sum _{n=1}^N \varphi _n \theta _n(z), \end{aligned}$$

(130)

where $\theta _n$ is defined by

$$\begin{aligned} \begin{aligned}&\theta _n(z)=(1-z^2)T_{2n-2}(z),\\ \end{aligned} \end{aligned}$$

(131)

and $T_n(z)$ are the Chebyshev polynomials of the first kind, cf. [4]. The reason for the choice of the basis function $\theta _n$ is that in this way the functions $w$ and $\varphi$ satisfy the boundary conditions (129). The next step is to substitute expressions (130) into Eqs. (128) and (129), and then require that Eqs. (128) and (129) be satisfied at $N$ collocation points $z_1,\ldots ,z_N,$ where $z_i$ are defined by $z_i=\cos \,([i-1]/[2N-1]\pi ),\,\,i=1,\ldots ,N.$ This results in a $2N\times 2N$ system of algebraic equations in the coefficients $w_1,\ldots ,w_N, \varphi _1,\ldots ,\varphi _N$:

$$\begin{aligned} \left( \begin{array}{cc} D &{} -I\\ \ A_{1} &{} A_{2} \\ \end{array} \right) X =c \left( \begin{array}{cc} O &{} O\\ i a^{3} R I &{} -i a R D\\ \end{array} \right) X, \end{aligned}$$

where $X=(w_{1},\ldots , w_{N}, \varphi _{1},\ldots , \varphi _{N})$, $O$ is the zeros matrix,

$$\begin{aligned} A_{1}(n_{1},n_{2})= & {} (a^4+a^2 M^2+i a RU{^{\prime \prime }}(z_{n_{1}})+i a^{3} RU(z_{n_{1}}))I(n_{1},n_{2}),\\ A_{2}(n_{1},n_{2})= & {} D^{3}(n_{1},n_{2})-(2a^2+M^2+i a RU(z_{n_{1}}))D(n_{1},n_{2})\\ I(n_{1},n_{2})= & {} \theta _{n_{2}}(z_{n_{1}}),\, D(n_{1},n_{2})= \theta {^\prime }_{n_{2}}(z_{n_{1}}),\, D^{3}(n_{1},n_{2})= \phi {^{\prime \prime \prime }}_{n_{2}}(z_{n_{1}}),\,\\ n_{1}= & {} 1,\ldots ,N,\,n_{2}=1,\ldots ,N. \end{aligned}$$

For a given wave number $a$ and Reynolds number $R$ we solve the above generalized eigenvalue problem using he Matlab QZ algorithm, obtaining the spectrum $\{c_{1}, c_{2}, \ldots \}$ ordered such that $\mathrm {Imag}\{c_{1}\} \ge \mathrm {Imag} \{c_{2}\}$. We iterate over $R$ using secant method until we have $\mathrm {Imag} \{c_{1}\}=0$ (and therefore neutral stability), and then repeat, iterating over a with golden section search until we have the minimum point $(a_{L},R_{L})$ such that the system is linearly stable for all $R < R_{L}$.

Method 2

In order to implement the second technique we approximate the solutions to Eqs. (128) and (129) as truncated series of Chebyshev polynomials as follows,

$$\begin{aligned} w=\sum _{n=0}^Nw_nT_n(z),\qquad \varphi =\sum _{n=0}^N \varphi _nT_n(z). \end{aligned}$$

(132)

These expressions are employed in Eqs. (128) and (129) and then the resulting equations are evaluated at Gauss–Lobatto points $y_i$ defined by $y_i=\cos (\pi i/[N-3]),\,\,i=0,\ldots ,N-2$. This leads to $2N-2$ algebraic equations for $2N+2$ unknowns $w_0,\ldots ,w_N, \varphi _0,\ldots , \varphi _N.$ The remaining four equations are furnished by the boundary conditions (129) which become

$$\begin{aligned} \begin{aligned}&\sum _{n=0}^N \varphi _n=0,\qquad \sum _{n=0}^N(-1)^n \varphi _n=0,\\&\sum _{n=0}^N \varphi _n=0,\qquad \sum _{n=0}^N(-1)^n \varphi _n=0, \end{aligned} \end{aligned}$$

(133)

and these equations are added as rows to the matrices generated above to yield a $(2N+2)\times (2N+2)$ matrix eigenvalue equation. Then, we obtain the generalised eigenvalue problem:

$$\begin{aligned} \left( \begin{array}{cc} D &{} -I\\ BC_{1} &{} 0\ldots 0\\ BC_{2} &{} 0\ldots 0\\ A_{1} &{} A_{2} \\ 0\ldots 0 &{} BC_{3}\\ 0\ldots 0 &{} BC_{4}\\ \end{array} \right) X =c \left( \begin{array}{cc} O &{} O\\ 0\ldots 0 &{} 0\ldots 0\\ 0\ldots 0 &{} 0\ldots 0\\ i a^{3} R I &{} -i a R D\\ 0\ldots 0 &{} 0\ldots 0 \\ 0\ldots 0 &{} 0\ldots 0 \\ \end{array} \right) X, \end{aligned}$$

where $X=(w_{0},\ldots , w_{N}, \varphi _{0},\ldots , \varphi _{N})$, $O$ is the zeros matrix,

$$\begin{aligned} A_{1}(n_{1},n_{2})= & {} (a^4+a^2m^2+i a RU{^{\prime \prime }}(y_{n_{1}})+i a^{3} RU(y_{n_{1}}))I(n_{1},n_{2}),\\ A_{2}(n_{1},n_{2})= & {} D^{3}(n_{1},n_{2})-(2a^2+m^2+i a RU(y_{n_{1}}))D(n_{1},n_{2})\\ I(n_{1},n_{2})= & {} T_{n_{2}}(y_{n_{1}}),\,\,\,\,\,\,\,\,\,\, \,\,\,\,\,D(n_{1},n_{2})= T{^\prime }_{n_{2}}(y_{n_{1}}),\\ D^{3}(n_{1},n_{2})= & {} T{^{\prime \prime \prime }}_{n_{2}}(y_{n_{1}}),\\ n_{1}= & {} 0,\ldots ,N-2,\,\,\,\,\,\,\,\,n_{2}=0,\ldots ,N. \end{aligned}$$

This is solved by the QZ algorithm.

Finite Element Method

In this section we will discuss the finite element method which is used to solve systems (126) and (127). Firstly, we introduce a new variable $\chi =D^{2}w$, then, system (126) has the following form

$$\begin{aligned}&D^{2}w=\chi , \end{aligned}$$

(134)

$$\begin{aligned}&\quad D^{2}\chi -(2a^2+m^2)\chi +(a^4+a^2m^2)w+i a RU{^{\prime \prime }}w-i a RU(\chi -a^2w)\nonumber \\&\qquad =-i a Rc(\chi -a^2w), \end{aligned}$$

(135)

with the boundary conditions

$$\begin{aligned} w=0,\,\,\,\,\,\,\,\,Dw=0,\,\,\,\,\,\,\,\, \mathrm {at} \,\,\,\,\,\,\,\,z=\pm 1. \end{aligned}$$

(136)

Using the finite element procedure which is shown in “Finite Element Method” section, the matrix representation for the system of equations of element $e$ take the form

$$\begin{aligned}&\left( \begin{array}{cc} -D_{2}^{e} &{} -F_{1}^{e} \\ (a^4+a^2m^2)F_{1}^{e}+i a RF_{3}^{e}+i a^{3} RF_{2}^{e} &{} -D_{2}^{e}-(2a^2+m^2)F_{1}^{e}-i a RF_{2}^{e}\\ \end{array} \right) \left( \begin{array}{c} \delta _{W}^{e} \\ \delta _{\chi }^{e} \\ \end{array} \right) \nonumber \\&\quad =c\left( \begin{array}{cc} O &{} O \\ i R a^{3} F_{1}^{e} &{}-i a R F_{1}^{e}\\ \end{array} \right) \left( \begin{array}{c} \delta _{W}^{e} \\ \delta _{\chi }^{e} \\ \end{array} \right) , \end{aligned}$$

(137)

where

$$\begin{aligned} O_{ij}= & {} 0,\,\,\,D_{2\,ij}^{e}=\int _{-1}^{1}\frac{d N_{n_{i}}^{p}}{d z} \frac{d N_{n_{j}}^{p}}{dz} dz, \,\,\,\,\,\,F_{1\,ij}^{e} =\int _{-1}^{1}N_{n_{i}}^{p} N_{n_{j}}^{p} dz,\\ F_{2\,ij}^{e}= & {} \int _{-1}^{1}U(z)\,N_{n_{i}}^{p} N_{n_{j}}^{p}dz, F_{3\,ij}^{e} =\int _{-1}^{1}U{^{\prime \prime }}(z)\,N_{n_{i}}^{p} N_{n_{j}}^{p}dz,\,\,\,\,\, i=1,\ldots ,p+1,\\&j=1,\ldots ,p+1. \end{aligned}$$

The above integral can be evaluated by applying the classical finite element nodal basis functions in one dimension on the standard element $\Omega _{st}=(-1,1)$. The standard shape functions are defined by the set of Lagrange polynomials

$$\begin{aligned} N_{n_{i}}^{p}(\zeta )=\prod _{j=n_{1},\, j\ne i}^{n_{p+1}} \frac{\zeta -\zeta _{j}}{\zeta _{i}-\zeta _{j}}, \end{aligned}$$

(138)

where $\zeta =2/h(z-z_{m})$, $h$ is the length of element and $z_{m}$ is the mid-point. Thus these integrations take the form:

$$\begin{aligned} D_{2\,ij}^{e}= & {} \frac{2}{h}\int _{-1}^{1}\frac{d N_{n_{i}}^{p}}{d \zeta } \frac{d N_{n_{j}}^{p}}{d \zeta } d\zeta ,\\ F_{1\,ij}^{e}= & {} \frac{h}{2}\int _{-1}^{1}N_{n_{i}}^{p} N_{n_{j}}^{p} d\zeta \\ F_{2\,ij}^{e}= & {} \frac{h}{2}\int _{-1}^{1}U\left( \frac{h}{2}\zeta +z_{m}\right) \,N_{n_{i}}^{p} N_{n_{j}}^{p} d\zeta ,\\ F_{3\,ij}^{e}= & {} \frac{h}{2}\int _{-1}^{1}U{^{\prime \prime }}\left( \frac{h}{2}\zeta +z_{m}\right) \,N_{n_{i}}^{p} N_{n_{j}}^{p} d\zeta ,\\&i=1,\ldots ,p+1,\,\,\,\,\,\,\,\,\,\,\,\,\,\,j=1,\ldots ,p+1. \end{aligned}$$

All these integral was calculated analytically using Matlab routines. Finally, we assemble the systems of all elements $e=1,\ldots ,n$ to get the main eigenvalue system which have the form

$$\begin{aligned} \Lambda \,\Pi =c\, \Xi \,\Pi , \end{aligned}$$

(139)

where $\Pi =(w_{1},\ldots , w_{m}, \chi _{1},\ldots , \chi _{m})^{T}$.

Now, we deal with the boundary conditions using the same techniques which have been used previously in the finite element section.

Table 3 The critical Rayleigh and wave numbers which are evaluated by using Chebyshev collocation method-1, Chebyshev collocation method-2 and finite element method for a selection of $M$ values

Numerical Results of Poiseuille Flow and Conclusions

The numerical results reported are based on the leading eigenvalue of the system (128) and (129). By this we mean that when one employs the time representation in $w$ and $\varphi$ of form $e^{-ict}$ with $c=c_r+ic_i$ then this results in $w$ and $\varphi$ having terms of form $\exp (-ic_rt).\exp (c_it).$ The eigenvalues are found such that the largest value of $c_i$ is $c_i=0$ and then the result is minimized over the wavenumber $a$. The resulting $R$ value is then the critical Reynolds number with corresponding wavenumber. The value $c_i=0$ is chosen because this is the threshold at which the solution becomes unstable according to linearized theory. For, if $c_i>0$ then $w$ and $\varphi$ grow rapidly like $\exp {(c_it)}$ and the solution is unstable.

We have found the finite element method to be more unstable numerically than either of the two collocation methods reported above. However, we include in Table 3 numerical results to compare the performance of the collocation and finite element methods. All methods require more care as $M$ increases and value of corresponding critical Rayleigh number increase with $M$. The values reported in Table 3 give the number of polynomials and elements required to achieve the same level of accuracy for all three methods. From Table 3, we see that as $M$ increases the number of Chebyshev polynomials increases for both collocation methods-1 and -2. For the finite element method both the polynomial order of the element and the number of elements must strongly increase as $M$ increases and, indeed, for $M$ greater than 5 we found it impossible to obtain satisfactory results, whereas the collocation methods still worked.

The spectrum which is plotted in Figs. 10 and 11, is similar to that found in Poiseuille flow in a porous medium with no-slip boundary conditions by [24], the eigenvalues displaying a $Y$ shape in the $(c_r,c_i)$ diagram. As $M$ increase, the eigenvalues at the intersection of the three lines in the $Y$ become more numerically unstable, and as $M$ increases this effect is very pronounced. The spectrum of (128) and (129) behaves very like that of the Orr–Sommerfeld problem for classical Poiseuille flow. For example, for higher Reynolds numbers we witnessed mode crossing of eigenvalues. For example, for $M = 0,0.5,$ and $1.0,$ the first and second eigenvalues interchange places for $R$ between 80822 and 80828, 86852 and 86854, and 106618 and 106620, respectively, with the previous first eigenvalue moving down the list as R increases. This behaviour is very similar to that observed by [4]. Moreover, the spectrum is very sensitive and care must be taken with the number of polynomials used in the numerical approximation, and in the arithmetical precision used in the calculation (those presented here are all in 64 bit arithmetic). The spectrum for $M = 2$ at the critical values are shown in Figs. 10 and 11 for Chebyshev collection method $1$ and $2$, respectively. Since $U$ is an even function of $z$, the proper solution of eigenvalue system (128) and (129) falls into two non-combining groups of even and odd solutions. However, Chebyshev collection method $1$ produce the approximate eigenvalues with even modes for plane Poiseuille flow, while the Chebyshev collection method $2$ give the approximate eigenvalues with even and odd modes, and this behaviour is clear in 10 and 11.

Fig. 10

Spectral of growth rate $c=c_{r}+ic_{i}$ at critical values and $M=2,$ by using Chebyshev collocation method-1. Here, the approximate eigenvalues associated with even modes for plane Poiseuille flow. The number of polynomials used in the numerical approximation correspond to a 40, b 80, c 120, d 200, e 300, f 400

Fig. 11

Spectral of growth rate $c=c_{r}+ic_{i}$ at critical values and $M=2,$ by using Chebyshev collocation method-2. Here, the approximate eigenvalues associated with even and odd modes for plane Poiseuille flow. The number of polynomials used in the numerical approximation correspond to a 100, b 200, c 300, d 400, e 500, f 600