Onset of Convection in an Inclined Anisotropic Porous Layer with Internal Heat Generation

Abstract

The onset of convection in an inclined porous layer which is heated internally by a uniform distribution of heat sources is considered. We investigate the combined effects of inclination, anisotropy and internal heat generation on the linear instability of the basic parallel flow. When the Rayleigh number is sufficiently large, instability occurs and a convective motion is set up. It turns out that the preferred motion at convection onset depends quite strongly on the anisotropy ratio, $\xi$, and the inclination angle. When $\xi <1$ the preferred motion is in the form of longitudinal rolls for all inclinations. When $\xi >1$ transverse rolls are preferred for small inclinations but, at high inclinations, longitudinal rolls are preferred. At intermediate inclinations the preferred roll orientation varies smoothly between these two extremes.

1. Introduction

The onset of convection in porous layers forms one of the classical problems of stability theory and the most well-known of these is the the Darcy-Bénard (or Horton-Rogers-Lapwood) problem [1,2]. In its simplest form a horizontal porous layer for which Darcy’s law applies has a lower boundary which is impermeable and uniformly hot while the impermeable upper boundary is uniformly cold. It is well-known that the critical Darcy-Rayleigh number for the onset of convection is $4{\pi }^2$ when the layer is unbounded horizontally. The corresponding wavenumber is $\pi$. A second well-known problem of the same type comprises a porous layer with uniformly cold boundaries but which is heated via a uniform internal heat generation mechanism. Practical examples of such systems include heat removal from fuel debris in nuclear reactors, the underground disposal of radioactive waste materials, agricultural applications involving the storage of foodstuff, and exothermic chemical reactions in packed-bed reactors. It is this configuration which we study here and we shall determine the combined effects of layer inclination and anisotropy on the onset of convection.

The first experimental and theoretical studies on the topic of an internally-heated layer were undertaken in the late 1970s by Kulacki and Ramchandani [3], Gasser and Kazimi [4], Buretta and Berman [5], Hardee and Nilson [6], Hwang and Marr [7], Tveitereid [8], Rhee et al. [9] and Kulacki and Freeman [10]. These papers represent quite a wide range of different types of boundary condition such as a free upper surface in [4] or the presence of an insulated bounding surface [8,10]. The case where both surfaces are impermeable and are at identical constant temperatures has a critical Darcy-Rayleigh number of $471.38$ and a critical wavenumber of 4.6752, values which have been obtained by various authors such as Barletta et al. [11], Nandal and Mahajan [12] and Nouri-Borujerdi et al. [13,14]. This ten-fold rise in the critical Darcy-Rayleigh number may be explained at least qualitatively by the facts that (i) the internally-heated case admits a destabilising temperature gradient only in the upper half of the layer and (ii) the maximum temperature difference is only $\frac{1}{8}$ of that of that of the classical Darcy-Bénard layer. The increase in the wavenumber may be traced to the fact that convection cells tend to remain of roughly a unit aspect ratio, and given that most of the disturbance occuies the upper half of the layer its horizontal period is reduced from that of the Darcy-Bénard problem.

Many authors have considered different extensions to this ‘pure’ internally-heated problem as described above. Nield and Kuznetsov [15] considered different types of nonuniform internal heating profiles, while Nield and Kuznetsov [16] and Kuznetsov and Nield [17] studied the effects of inhomogeneity in the structure of the porous medium itself. Matta and Hill [18], on the other hand, considered an equivalent double-diffusive configuration where the heat source depends on the solutal concentration. Nouri-Borujerdi et al. [13] considered the equivalent Darcy-Brinkman case and the critical Darcy-Rayleigh number was found to rise as the Darcy number increases. Nouri-Borujerdi et al. [14] relaxed the assumption of local thermal equilibrium between the phases. In addition Mahajan and Nandal [19] considered how anisotropy and Brinkman effects alter the stability criteria, while Yadav et al. [20] consider the combined effects of rotation and the Brinkman terms. We also mention some works which include anisotropy but with heating from below: Malashetty and Swamy [21] and Yavogan et al. [22] both analysed the effect of a viscoelastic fluid and obtained an analytical expression for the critical Darcy-Rayleigh number, Raghunatha et al. [23] considered the case when the porous medium is saturated by an Oldroyd-B fluid, Abdelhafez and Tsybulin [24] who considered the presence of a binary fluid, Yadav and Kim [25] who studied the instability of a time-evolving boundary layer flow and Storesletten and Rees [26] on the instability of a boundary layer induced by uniform surface suction.

The present paper consists of a linear stability analysis of the above-described uniformly-heated porous layer, but we also consider the extra effects of anisotropy and inclination. It is well-known from studies of the Darcy-Bénard problem that the introduction of anisotropy by the modification of one of the principal permeabilities in the horizontal direction removes the directional invariance of the roll pattern to one where there is a preferred direction. Thus if the y-direction is vertical and the altered permeability is in the x-direction, then the axis of resulting rolls will be in the x-direction if the permeability in the x-direction is larger than the other principal permeabilities, and in the z-direction otherwise. The effect of inclination on the Darcy-Bénard problem is to cause longitudinal rolls, i.e., those with axes up the layer, to be preferred. The combination of these two effects, especially when they conflict with regard to their prediction of the preferred roll direction, is of some great interest. In this regard the works by Rees and Postelnicu [27], Postelnicu and Rees [28] and Rees et al. [29] are important for they display both sudden and smooth transitions in the preferred roll direction as the inclination is increased.

This work, then, represents the combined effects of internal heating, anisotropy and inclination. It may be regarded as being an extension (i) of the works of Rees and Postelnicu [27] and Storesletten and Tveitereid [29] by having the presence of internal heating, (ii) of the work of Barletta et al. [11] by the inclusion of the effect of anisotropy and (iii) of the works of Mahajan and Nandal [19] and Capone et al. [30] by inclining the layer. We concentrate first on two-dimensional convection (transverse rolls) which will always arise whenever the porous layer is sufficiently narrow in the spanwise direction. Then we consider rolls of all orientations, the ultimate aim being to determine how the critical orientation varies with the governing parameters, namely the anisotropy ratio, $\xi$, and the inclination angle, $\varphi$.

2. Governing Equations and the Basic State

We consider the onset of convection in an inclined porous layer of infinite lateral extent which is heated internally by a uniform distribution of heat sources. The layer is bounded by two parallel solid plates which are separated by a distance h and inclined at an angle $\alpha$ to the horizontal. Both boundaries are held at a constant temperature $T_0$. The horizontal z-axis forms the direction about which the layer has been rotated, the y-axis is perpendicular to the bounding plates, and the x-axis lies in the lower surface pointing up the plane. The detailed configuration is shown in Figure 1. Moreover, the anisotropic porous layer is taken to be transversely isotropic in the permeability but homogenously isotropic in the thermal diffusivity. The saturating fluid is assumed to be such that Darcy’s law holds and the Oberbeck-Boussinesq approximation is valid. The governing equations are then,

$${\displaystyle \frac{\partial u^{\ast }}{\partial x^{\ast }}}+{\displaystyle \frac{\partial v^{\ast }}{\partial y^{\ast }}}+{\displaystyle \frac{\partial w^{\ast }}{\partial z^{\ast }}}=0,$$

(1)

$$u^{\ast }={\displaystyle \frac{K_L}{\mu }}\left(-{\displaystyle \frac{\partial p^{\ast }}{\partial x^{\ast }}}+{\rho }_{0}g\beta (T-T_0)sin\alpha \right),$$

(2)

$$v^{\ast }={\displaystyle \frac{K_T}{\mu }}\left(-{\displaystyle \frac{\partial p^{\ast }}{\partial y^{\ast }}}+{\rho }_{0}g\beta (T-T_0)cos\alpha \right),$$

(3)

$$w^{\ast }=-{\displaystyle \frac{K_T}{\mu }}{\displaystyle \frac{\partial p^{\ast }}{\partial z^{\ast }}},$$

(4)

$$\sigma {\displaystyle \frac{\partial T}{\partial t^{\ast }}}+u^{\ast }{\displaystyle \frac{\partial T}{\partial x^{\ast }}}+v^{\ast }{\displaystyle \frac{\partial T}{\partial y^{\ast }}}+w^{\ast }{\displaystyle \frac{\partial T}{\partial z^{\ast }}}=\kappa \left({\nabla }^{2}T+{\displaystyle \frac{q}{k_m}}\right),$$

(5)

where q is the rate of heat generated within the porous medium per unit volume, $\beta$ is the coefficient of cubical expansion, $k_m$ is the thermal conductivity of the porous medium, and p is the pressure, which includes the hydro-static pressure. Terms in Equations (1)–(5) have their usual meaning in the porous medium context, but we mentioned that $K_L$ is the permeability in the x-direction, $K_T$ is the permeability in both the y- and z-directions, that $\sigma$ is the heat capacity ratio, and $\kappa$ is the thermal diffusivity. A zero subscript denotes reference quantities. Further details are in the Abbreviations at the end of the paper. The variables are made dimensionless by setting

$$(x^{\ast },y^{\ast },z^{\ast })=h(x,y,z),(u^{\ast },v^{\ast },w^{\ast })={\displaystyle \frac{\kappa }{h}}(u,v,w),$$

(6)

$$T-T_0={\displaystyle \frac{q{h}^2}{k_m}}\theta ,p^{\ast }={\displaystyle \frac{\kappa \mu }{K_T}}p,t^{\ast }={\displaystyle \frac{\sigma h^2}{\kappa }}t,$$

(7)

so that the governing equations reduce to,

$${\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{\partial v}{\partial y}}+{\displaystyle \frac{\partial w}{\partial z}}=0,$$

(8)

$$u=-\xi {\displaystyle \frac{\partial p}{\partial x}}+\xi \mathrm{Ra}\theta sin\alpha ,v=-{\displaystyle \frac{\partial p}{\partial y}}+\mathrm{Ra}\theta cos\alpha ,w=-{\displaystyle \frac{\partial p}{\partial z}},$$

(9)

$${\displaystyle \frac{\partial \theta }{\partial t}}+u{\displaystyle \frac{\partial \theta }{\partial x}}+v{\displaystyle \frac{\partial \theta }{\partial y}}+w{\displaystyle \frac{\partial \theta }{\partial z}}={\nabla }^2\theta +1.$$

(10)

The dimensionless parameters which appear in these equations are the Darcy-Rayleigh number, $\mathrm{Ra}$, and the anisotropy parameter, $\xi$, which are given respectively by,

$$\mathrm{Ra}={\displaystyle \frac{{\rho }_{0}g\beta K_{T}h}{\kappa \mu }}\left({\displaystyle \frac{q{h}^2}{k_m}}\right),\xi ={\displaystyle \frac{K_L}{K_T}}.$$

(11)

In this paper we shall use two different formulations of the governing equations for the numerical work. In the first instance two-dimensional cases will be considered using the streamfunction/temperature approach whereas for three-dimensional flows a pressure/temperature approach will be used. In the former case the streamfunction is defined using

$$u=-{\displaystyle \frac{\partial \psi }{\partial y}},v={\displaystyle \frac{\partial \psi }{\partial x}},w=0,$$

(12)

and therefore Equation (8) is satisfied while Equations (9) and (10) transform into,

$$\xi {\displaystyle \frac{{\partial }^2\psi }{\partial x^2}}+{\displaystyle \frac{{\partial }^2\psi }{\partial y^2}}=\xi \mathrm{Ra}\left({\displaystyle \frac{\partial \theta }{\partial x}}cos\alpha -{\displaystyle \frac{\partial \theta }{\partial y}}sin\alpha \right),$$

(13)

$${\displaystyle \frac{\partial \theta }{\partial t}}={\displaystyle \frac{{\partial }^2\theta }{\partial x^2}}+{\displaystyle \frac{{\partial }^2\theta }{\partial y^2}}+{\displaystyle \frac{\partial \psi }{\partial y}}{\displaystyle \frac{\partial \theta }{\partial x}}-{\displaystyle \frac{\partial \psi }{\partial x}}{\displaystyle \frac{\partial \theta }{\partial y}}+1.$$

(14)

In the latter case $(u,v,w)$ may be eliminated from Equations (8)–(10) and we obtain the following governing equations for the unsteady flow:

$$\xi {\displaystyle \frac{{\partial }^{2}p}{\partial x^2}}+{\displaystyle \frac{{\partial }^{2}p}{\partial y^2}}+{\displaystyle \frac{{\partial }^{2}p}{\partial z^2}}=\mathrm{Ra}\left(\xi {\displaystyle \frac{\partial \theta }{\partial x}}sin\alpha +{\displaystyle \frac{\partial \theta }{\partial y}}cos\alpha \right),$$

(15)

$${\displaystyle \frac{\partial \theta }{\partial t}}+\mathrm{Ra}\theta \left(\xi {\displaystyle \frac{\partial \theta }{\partial x}}sin\alpha +{\displaystyle \frac{\partial \theta }{\partial y}}cos\alpha \right)-\left(\xi {\displaystyle \frac{\partial p}{\partial x}}{\displaystyle \frac{\partial \theta }{\partial x}}+{\displaystyle \frac{\partial p}{\partial y}}{\displaystyle \frac{\partial \theta }{\partial y}}+{\displaystyle \frac{\partial p}{\partial z}}{\displaystyle \frac{\partial \theta }{\partial z}}\right)={\displaystyle \frac{{\partial }^2\theta }{\partial x^2}}+{\displaystyle \frac{{\partial }^2\theta }{\partial y^2}}+{\displaystyle \frac{{\partial }^2\theta }{\partial z^2}}+1.$$

(16)

For these two formulations we may write the boundary conditions in the form,

$$\psi ={\displaystyle \frac{\partial p}{\partial y}}=\theta =0\mathrm{on}\mathrm{both}y=0\mathrm{and}y=1.$$

(17)

There exists a unique steady basic flow, $\mathbf{v}=(u_b,0,0)$, $\psi ={\psi }_b\left(y\right)$, $\theta ={\theta }_b\left(y\right)$ and $p=p_b\left(y\right)$ where the net mass flux through a cross section $x=\mathrm{constant}$ is zero. This steady flow with nonlinear velocity and temperature profiles is given by,

$$u_b=\xi \mathrm{Ra}(-\frac{1}{2}{y}^2+\frac{1}{2}y-\frac{1}{12})sin\alpha ,{\psi }_b=\xi \mathrm{Ra}(\frac{1}{6}{y}^3-\frac{1}{4}{y}^2+\frac{1}{12}y)sin\alpha ,$$

(18)

$${\theta }_b=\frac{1}{2}(y-y^2),p_b=\mathrm{Ra}[\frac{1}{12}xsin\alpha +(\frac{1}{4}{y}^2-\frac{1}{6}{y}^3)cos\alpha ]+p_0,$$

(19)

where $p_0$ is the pressure at $(x,y,z)=(0,0,0)$.

3. Linear Stability Analysis for Two-Dimensional Disturbances

3.1. Reduction to ODE Eigenvalue Form

The linear stability characteristics of the basic steady state given by Equations (18) and (19) may be investigated by setting,

$$\psi ={\psi }_b+{\psi }_1,\theta ={\theta }_b+{\theta }_1,$$

(20)

where ${\psi }_1$ and ${\theta }_1$ are asymptotically small perturbations in the streamfunction and temperature, respectively. The substitution of the expressions (20) into the governing Equations (13) and (14) followed by linearisation yields the following set of perturbation equations,

$$\xi {\displaystyle \frac{{\partial }^2{\psi }_1}{\partial x^2}}+{\displaystyle \frac{{\partial }^2{\psi }_1}{\partial y^2}}=\xi \mathrm{Ra}\left[{\displaystyle \frac{\partial {\theta }_1}{\partial x}}cos\alpha -{\displaystyle \frac{\partial {\theta }_1}{\partial y}}sin\alpha \right],$$

(21)

$${\displaystyle \frac{{\partial }^2{\theta }_1}{\partial x^2}}+{\displaystyle \frac{{\partial }^2{\theta }_1}{\partial y^2}}=\left[\frac{1}{2}-y\right]{\displaystyle \frac{\partial {\psi }_1}{\partial x}}-\xi \mathrm{Ra}sin\alpha \left[\frac{1}{12}-\frac{1}{2}y+\frac{1}{2}{y}^2\right]{\displaystyle \frac{\partial {\theta }_1}{\partial x}}+{\displaystyle \frac{\partial {\theta }_1}{\partial t}},$$

(22)

which are to be solved subject to homogeneous Dirichlet conditions. The disturbances may be reduced to ordinary differential form by factoring out a monochromatic horizontal component: we let,

$$\left(\begin{array}{c}{\psi }_1\\ {\theta }_1\end{array}\right)=e^{ikx+\lambda t}\left(\begin{array}{c}i\Psi \left(y\right)\\ \Theta \left(y\right)\end{array}\right)+\mathrm{complex}\mathrm{conjugate}.$$

(23)

Hence the disturbance equations reduce to,

$${\Psi }^{″}-\xi k^2\Psi =\xi \mathrm{Ra}\left[k\Theta cos\alpha +i{\Theta }^{\prime }sin\alpha \right],$$

(24)

$${\Theta }^{″}-k^2\Theta =-(\frac{1}{2}-y)k\Psi -ik\xi \mathrm{Ra}sin\alpha \left[\frac{1}{12}-\frac{1}{2}y+\frac{1}{2}{y}^2\right]\Theta +\lambda \Theta ,$$

(25)

and the boundary conditions are that,

$$\Psi =\Theta =0\mathrm{on}\mathrm{both}y=0\mathrm{and}y=1.$$

(26)

This system may be regarded as an eigenvalue problem for the complex growth rate, $\lambda$, as a function of k, $\xi$, $\mathrm{Ra}$ and $\alpha$. Alternatively, if we let $\lambda ={\lambda }_r+i{\lambda }_i$ and set ${\lambda }_r=0$ to model marginal stability, then this system may instead be regarded as an eigenvalue problem for both the Darcy-Rayleigh number, $\mathrm{Ra}$, and for ${\lambda }_i$ as functions of k, $\xi$ and $\alpha$.

We may now define the phase velocity of the disturbances by modifying the exponential in Equation (23). At marginal stability $\lambda =i{\lambda }_i$ and therefore,

$$e^{ikx+\lambda t}=e^{ik(x-ct)}\mathrm{where}c=-{\lambda }_i/k.$$

(27)

3.2. Numerical Results

In the first instance Equations (24)–(26) were discretised using second order accurate central difference approximations and the resulting difference equations were rearranged into a complex matrix eigenvalue problem for $\lambda$ precisely in the manner described in Rees and Bassom [31]. For a chosen pair of values of $\xi$ and inclination, $\alpha$, the first few eigenvalues for $\lambda$ were computed on a fine grid of values of $\mathrm{Ra}$ and k, and contours of ${\lambda }_r=0$ were plotted to provide the neutral curves. The use of $N=10$ intervals in the range $0\le y\le 1$ gave poor results particularly when $\mathrm{Ra}$ or k was too large. The number of intervals was then increased first to $N=20$, then to $N=40$ and ultimately to $N=80$, we found that, that the neutral curves presented in Figure 2 and which correspond to $N=40$ exhibit no visual changes when comparing them with $N=80$ solutions.

In more detail, Figure 2 displays four sets of neutral curves for $\xi =5$, 2, $0.5$ and $0.2$. For this set of values of $\xi$ the various frames of Figure 2 display the manner in which the neutral curves evolve as $\alpha$ increases in $1^{\circ }$ increments from the horizontal case, $\alpha =0^{\circ }$. For ease of comparison each subfigure corresponds to k and $\mathrm{Ra}$ lying in the ranges, $0\le k\lesssim 3k_c$ and $0\le \mathrm{Ra}\lesssim 12{\mathrm{Ra}}_c$ where $k_c$ and ${\mathrm{Ra}}_c$ are the critical wavenumber and Darcy-Rayleigh number when the layer is horizontal ($\alpha =0$). We have plotted the neutral curves corresponding to the first two modes.

It is clear from Figure 2 that the general form taken by the neutral curves remains almost unchanged as $\xi$ varies. Thus the critical values of $\mathrm{Ra}$ (which correspond to the minimum of each neutral curve and which are labelled as ${\mathrm{Ra}}_c$ hereinafter) always increase as the inclination, $\alpha$, increases and this may be expected given that buoyancy decreases with increasing inclination. For three of the cases shown in Figure 2 there is a local maximum inclination angle which is at an isola point (shown as black disks). These points are characterised by being surrounded by elliptic-shaped neutral curves where the $\xi =2$ case provides a good example since the value of $\alpha$ at the isola itself lies just above ${27}^{\circ }$. In these three cases ($\xi =5$, 2 and $0.5$) there also exists a saddle point; for $\xi =0.5$ this is depicted as a red disk but the saddle points for $\xi =5$ and $\xi =2$ lie above the maximum value of $\mathrm{Ra}$ used for plotting.

Figure 3 gives an alternative view which both summarises and extends the above discussion and shows, for a wide selection of values of $\xi$, how the local minimum and maximum of $\mathrm{Ra}$ on each neutral curve varies with inclination. The lower set of black disks and the associated dotted curve correspond to how the value of $\mathrm{Ra}$ for the isolae varies with inclination, and the upper set correspond to the saddle points. These two curves eventually merge when $\xi$ is roughly $0.27843$, which is in accord with the fact that no isola is displayed in Figure 2 for $\xi =0.2$. We presume that each of the solid curves shown in Figure 3 tend asymptotically toward a constant value of $\alpha$ as ${\mathrm{Ra}}_c\to \infty$, as was shown in Rees and Bassom [31] for the inclined Darcy-Bénard problem.

Some indication of how the critical wavenumber varies with $\alpha$ may be gleaned from Figure 2, but a more informative depiction is shown in Figure 4. Thus for relatively large values of $\xi$ the critical wavenumber decreases at first as $\alpha$ increases from zero whereas $k_c$ increases when $\xi$ takes smaller values. The presence of isolae is marked with black disks, and now it is clear that they correspond to a local maximum value of $\alpha$ in those curves. On the other hand, the saddle points, which are also marked with black disks, correspond to local minima. As $\xi$ decreases the curves shown in Figure 4 gradually unkink themselves as the isola and the saddle point approach one another. Eventually, $k_c$ becomes a single-valued function of $\alpha$.

In Figure 5 we show how the corresponding phase velocity, c, of the disturbances vary with $\alpha$ for the same set of values of $\xi$. Postive values mean that the disturbance ascends the layer when $\alpha$ is positive while negative values mean that the disturbance descends. Figure 5 shows that disturbances tend to ascend the layer when the inclination is small, but eventually descend when the inclination is sufficiently large. The reason for this behaviour is related to the overall basic flow pattern and specifically to the location of the disturbance within the layer. Using Equation (18) we see that

$$u_b\left(0\right)=u_b\left(1\right)=-\frac{1}{12}\xi \mathrm{Ra}sin\alpha ,u_b\left(\frac{1}{2}\right)=\frac{1}{24}\xi \mathrm{Ra}sin\alpha ,$$

(28)

and therefore the basic state corresponds to descending flow near the boundaries and ascending flow in the centre.

Figure 6 clarifies the reason why the phase velocity varies as it does in Figure 5. This figure shows instantaneous snapshots of the disturbance isotherms for six different inclinations when $\xi =5$. The first is the horizontal case where $c=0$, while the second, third and fourth correspond to positive values of c where the fourth also corresponds to the maximum phase velocity. For larger inclinations c then decreases. The fifth frame corresponds to $c=0$ while the sixth is at the isola point. Disturbances tend to be in the upper half of the layer because this is where the basic temperature gradient is destabilising. When inclinations are sufficiently small, the disturbances lie mostly in the region where the flow is up the layer. Moreover, the expression for $u_b\left(\frac{1}{2}\right)$ in Equation (28) shows that this value increases as $\alpha$ increases, and therefore the phase velocity of the disturbances also tends to increase with $\alpha$ at first. But as $\alpha$ increases further, disturbances become increasingly concentrated near the upper boundary where the basic flow is negative.

When other values of $\xi$ are considered then figures like those shown in Figure 6 are obtained. The main change in the appearance of the analogous isotherm contours is the period of the cellular pattern because the critical wavenumbers vary quite strongly with $\xi$, as may be seen in Figure 4.

4. Linear Stability Analysis for Three-Dimensional Disturbances

4.1. Reduction to ODE Eigenvalue Form

The linear stability characteristics of the basic steady state to generally oblique roll disturbances may be investigated by setting,

$$p=p_b+p_1,\theta ={\theta }_b+{\theta }_1,$$

(29)

where $p_1$ and ${\theta }_1$ are asymptotically small perturbations in the pressure and temperature, respectively. The substitution of the expressions (29) into the governing Equations (15) and (16) gives the following set of linearised perturbation equations

$$\xi {\displaystyle \frac{{\partial }^2{p}_1}{\partial x^2}}+{\displaystyle \frac{{\partial }^2{p}_1}{\partial y^2}}+{\displaystyle \frac{{\partial }^2{p}_1}{\partial z^2}}=\mathrm{Ra}\left(\xi {\displaystyle \frac{\partial {\theta }_1}{\partial x}}sin\alpha +{\displaystyle \frac{\partial {\theta }_1}{\partial y}}cos\alpha \right),$$

(30)

$$\begin{array}{c}{\displaystyle \frac{\partial {\theta }_1}{\partial t}}+\mathrm{Ra}{\theta }_b\left(\xi {\displaystyle \frac{\partial {\theta }_1}{\partial x}}sin\alpha +{\displaystyle \frac{\partial {\theta }_1}{\partial y}}cos\alpha \right)+\mathrm{Ra}{\theta }_1{\displaystyle \frac{\partial {\theta }_b}{\partial y}}cos\alpha \\ -\left(\xi {\displaystyle \frac{\partial p_b}{\partial x}}{\displaystyle \frac{\partial {\theta }_1}{\partial x}}+{\displaystyle \frac{\partial p_b}{\partial y}}{\displaystyle \frac{\partial {\theta }_1}{\partial y}}+{\displaystyle \frac{\partial p_1}{\partial y}}{\displaystyle \frac{\partial {\theta }_b}{\partial y}}\right)={\nabla }^2{\theta }_1\end{array}$$

(31)

We Fourier-decompose the perturbations in the x- and z-directions by setting

$$\left(\begin{array}{c}{p}_1\\ {\theta }_1\end{array}\right)=\left(\begin{array}{c}P\left(y\right)\\ \Theta \left(y\right)\end{array}\right)exp[ik(xsin\varphi +zcos\varphi )+\lambda t],$$

(32)

where $\varphi$ is the orientation of the roll relative to the x-axis, and where k and $\lambda$ are once more the wave number of the vortex/cell amd the exponential growth rate, respectively. Expressions (32), when substituted into Equations (30) and (31), reduce the perturbation equations to a set of ordinary differential equations for the reduced disturbance amplitudes, P and $\Theta$:

$$P^{″}-k^2\left(\xi {sin}^2\varphi +{cos}^2\varphi \right)P=\mathrm{Ra}\left(ik\xi \Theta sin\varphi sin\alpha +{\Theta }^{\prime }cos\alpha \right),$$

(33)

$${\Theta }^{″}-k^2\Theta =\frac{1}{2}\mathrm{Ra}(-y^2+y-\frac{1}{6})ik\xi \Theta sin\varphi sin\alpha +\frac{1}{2}\mathrm{Ra}(1-2y)\Theta cos\alpha -\frac{1}{2}(1-2y)P^{\prime }+\lambda \Theta .$$

(34)

The boundary conditions are,

$$P^{\prime }\left(0\right)=\Theta \left(0\right)=0,P^{\prime }\left(1\right)=\Theta \left(1\right)=0.$$

(35)

4.2. Numerical Results

Equations (33) and (34) subject to boundary conditions (35) were solved using precisely the same methodology as was used in [27,29], which is a variant of the well-known Keller box method. First, Equation (33) was differentiated once to obtain a second order equation in $P^{\prime }$ so that the boundary conditions for both variables, $P^{\prime }$ and $\Theta$ are Dirichlet and thereby absolute accuracy is increased. Second, the resulting equations are supplemented by a system obtained by differentiating these equations with respect to k where $\partial \mathrm{Ra}/\partial k=0$ is insisted upon in order to minimise the critical Rayleigh number with respect to the wavenumber. Finally the full system was discretised using second order accurate central differences and solved using the modified Keller box method which is described in [27]. This method allows for the computation of eigenvalues which, in this case, are ${\mathrm{Ra}}_c$, $k_c$, ${\lambda }_i$ and $\partial {\lambda }_i/\partial k$. The usual block tridiagonal structure of the iteration matrix is replaced by with an addition row and column of blocks because of the presence of the eigenvalues. A curve-following methodology was also implemented which allows the neutral curve to be followed in $\mathrm{Ra}$-$\alpha$ space for each chosen value of $\xi$ and $\varphi$. Thus a fifth ‘eigenvalue’ is required which measures the length of the curve being followed. The code used for this work was cross-calibrated by solving the same equations using a standard fourth order Runge-Kutta scheme embedded in a boundary value problem solver.

Figure 7 displays the variation of the critical value of ${\mathrm{Ra}}_{c}cos\alpha$ with $\alpha$ for a range of roll orientations for $\xi =0.5$. For these three-dimensional cases we have chosen to use ${\mathrm{Ra}}_{c}cos\alpha$ rather than ${\mathrm{Ra}}_c$ in order to visualize more easily the behaviour of the system and also to avoid ${\mathrm{Ra}}_c$ tending towards infinity as the layer tends towards being vertical. In particular we obtain that ${\mathrm{Ra}}_c=471.17/cos\alpha$ with $k_c=4.675$. The message which arises from Figure 7 is that longitudinal rolls are preferred over all other orientations when $\xi <1$, a general result which accords with the analysis of [27]. In these circumstances there is more resistance to flow in the longitudinal or x-direction than in the other two directions, and it is natural therefore to expect the favoured rolls to have axes in the x-direction.

When $\xi >1$ the situation becomes more complicated. At small inclinations, i.e., when $\alpha$ is small, then transverse rolls ($\varphi ={90}^{\circ }$) form the preferred orientation, but longitudinal rolls are preferred at much higher inclinations. In the inclined Darcy-Bénard problem studied by Storesletten and Tveitereid [32] it was found that there exists a critical inclination below which transverse rolls are preferred and above which longitudinal rolls are preferred. In the more general problem studied by Rees and Postelnicu [27] the transition was found not to be sudden, but, when $\xi$ takes values which are not too large, there exists a smooth transition where the preferred roll orientation varies continuously between $\varphi ={90}^{\circ }$ and $\varphi =0^{\circ }$ as $\alpha$ increases. It is this latter behaviour which we have found to occur in the present problem.

In Figure 8 and Figure 9 we display the variation of $\mathrm{Ra}cos\alpha$ with $\alpha$ for a range of values of $\varphi$ ($0^{\circ }$, $5^{\circ }$, ${10}^{\circ }$, ⋯, ${90}^{\circ }$) for $\xi =2$ and 5, respectively. In Figure 8 we see that the transition between transverse rolls at lower inclinations and longitudinal rolls at higher inclinations takes place close to $\alpha =\mathrm{}{19.6}^{\circ }$. This transition is smooth, as shown in the second frame of Figure 8, and it takes place very rapidly over an inclination range of roughly ${0.4}^{\circ }$.

Figure 9 shows the general behaviour of the neutral curves for the more anisotropic case, $\xi =5$. Again there exists a smooth transition between transverse and longitudinal rolls, but it now extends over a slightly larger inclination range of ${0.54}^{\circ }$ and is centred at roughly $\alpha ={15.75}^{\circ }$.

In both the cases considered above and for general values of $\xi$, it is possible to minimise the value of ${\mathrm{Ra}}_{c}cos\alpha$ with respect to $\varphi$ in order to determine the overall minimum critical Rayleigh number as a function of both the inclination and the anisotropy parameter. The results of this process are displayed in Figure 10 which shows how the critical Rayleigh number, and the corresponding roll orientations and wavenumbers vary, and we have selected six different values of $\xi >1$ to illustrate the general behaviour. The curves for $\xi =2$ and 5 follow the lower envelope of those shown in Figure 8 and Figure 9, respectively.

It is clear that the transition from transverse rolls ($\varphi ={90}^{\circ }$) to longitudinal rolls ($\varphi =0^{\circ }$) increases from $\alpha =0^{\circ }$ initially as $\xi$ increases from 1 upwards. When $\xi$ is close to 2 this increase in the transition inclination reverses direction and begins to decrease once more as $\xi$ increases further. The rapid variation in the transition process is captured well in the frames displaying $k_c$ and the roll orientation, $\varphi$.

Finally, it is of interest to determine an overall summary of the properties of the present layer in terms of the identity of the convection pattern which is to be expected in practice. This is given in Figure 11, which displays the regions in $(\alpha /\xi )$-space wherein longitudinal, oblique or transverse rolls form the preferred pattern at onset. The rapid transition between transverse and longitudinal rolls as $\alpha$ increases is now very evident due to the closeness of the two curves, and the above-mentioned maximum inclination angle corresponding to $\xi \simeq 2$ is also clear. The curves become extremely close when $\xi \to 1^{+}$ and our numerical data suggests that,

$${\displaystyle \frac{\alpha }{{(\xi -1)}^{1/2}}}\sim \left\{\begin{array}{cc}41.15-48.2(\xi -1)\hfill & \varphi =0^{\circ },\hfill \\ 41.15-49.2(\xi -1)\hfill & \varphi ={90}^{\circ }.\hfill \end{array}\right.$$

(36)

Therefore the transition range between these extreme orientations is roughly ${(\xi -1)}^{3/2}$ as $\xi \to 1$. On the other hand, when $\xi$ takes values which are larger than those shown in Figure 11 then the inclinations at which transitions happen continue to decrease. Some more extreme examples are given in

Table 1. Transitional values of $\alpha$ when $\xi$ is large; see Figure 11.

$\mathit{\xi }$	for $\mathit{\varphi }={90}^{\circ }$	for $\mathit{\varphi }=0^{\circ }$
125	$3.542304$	$3.660913$
250	$2.512439$	$2.595616$
500	$1.779291$	$1.837828$
1000	$1.259118$	$1.300388$

. These values were obtained by insisting that ${\partial }^2{\mathrm{Ra}}_c/\partial {\varphi }^2=0$, and from this data it is possible to shown (using some error analysis and Richardson extrapolation) that the transition zone (i.e., the range of values of $\alpha$ for which oblique rolls are preferred) occupies the range,

$$39.85{\xi }^{-1/2}\le \alpha \le 41.15{\xi }^{-1/2}.$$

(37)

when $\xi \gg 1$.

5. Conclusions

The onset of convection has been considered for a porous layer which is inclined to the horizontal, is anisotropic and is subject to a uniform internal heat generation. When the Rayleigh number is sufficiently large, instability occurs, and a convective motion is set up. The critical Rayleigh number at the onset of convection then depends very strongly on the anisotropy ratio as well as the inclination angle. Indeed the preferred orientation of convection rolls at onset can depend quite strongly on these parameters.

We have shown that, when $\xi <1$ then convection rolls are longitudinal (i.e., with axes in the x-direction) for all inclinations. Numerically we found that ${\mathrm{Ra}}_{c}cos\alpha =471.17$ and critical wavenumber $k_c=4.675$. However, when the flow is confined to being two dimensional, it has been shown that convection happens for only a finite range of inclinations which is dependent on the anisotropy parameter, $\xi$. Moreover, when $\xi >0.27843$ the neutral curves display an isola point meaning that there a local maximum inclination above which the basic flow is linearly stable.

When $\xi >1$ the general situation is more complicated. At small inclinations transverse rolls ($\varphi ={90}^{\circ }$) form the preferred orientation, but longitudinal rolls ($\varphi =0^{\circ }$) are preferred at much higher inclinations. In general, the critical Rayleigh number increases as the inclination of the layer increases, but the identity of the favoured roll orientation changes from transverse to longitudinal. The transition between these extreme orientations takes place smoothly but rather rapidly as the inclination, $\alpha$, increases. The associated variation in the critical wavenumbers becomes more extreme as the anisotropy ratio $\xi$ increases. We have also ascertained that the narrow transition region is centred on values of $\alpha$ which are close to zero both as $\xi ⟶1^{+}$ and as $\xi ⟶\infty$, and that transition region itself tends towards a zero width for each of those limits.