Modal Selection for Inclined Darcy-Bénard Convection in a Rectangular Cavity

Abstract

Nonlinear free convection in an inclined rectangular porous cavity heated from below has been studied using a two-dimensional spectral decomposition. The code uses pseudo-arclength continuation to follow solution curves around fold bifurcations. The evolution with inclination of the pattern of convection is complicated and it relies strongly on both the Darcy–Rayleigh number and the aspect ratio of the cavity. When the inclination is large it is generally true that only one cell appears, and that it has a circulation that is consistent with the direction of the buoyancy forces along the heated and cooled boundaries. However, as the inclination decreases back towards the horizontal, this unicellular pattern evolves, sometimes initially via fold bifurcations, into patterns with different numbers of cells. Such evolutions always conserve the parity of the number of cells (such as one cell becoming three and then five, or two cells becoming four), but bifurcations also arise between patterns with different parities. These phenomena are illustrated using a suitable selection of solution curves that show the dependence of the Nusselt number on the inclination.

1. Introduction

The instability of a uniform horizontal porous layer that is heated from below is the direct analogue of the Rayleigh-Bénard problem and is known as the Darcy-Bénard problem, the porous Bénard problem, or the Horton–Rogers–Lapwood problem. The first works on this linear instability problem were by Horton and Rogers [1], and Lapwood [2]. Since then researchers have extended these pioneering works, not only into the nonlinear regime, but also by adding further complications such as modified boundary conditions, internal heating, the presence of local thermal nonequilibrium, more than one diffusing species, non-Newtonian fluids, rotation, and the like. Given that the monograph by Nield and Bejan [3] devotes more than 100 pages to this, we shall refer the reader to their review of the state of the art.

The topic of the present work restricts the porous layer (which is often assumed to be unbounded horizontally) to that of a finite cavity, further restricts it to admitting solely two-dimensional flow, but allows the cavity to be inclined. Once more, Nield and Bejan [3] devote more than six pages to a review of the general topic of inclined layers, and once more, we refer the reader to their monograph. Rather, we shall review briefly those papers which are of very direct relevance to our work.

The first paper which considered the onset of two-dimensional convection in an inclined porous layer and which presented an exhaustive account of the linearized stability theory was by Rees and Bassom [4]. Using a finite difference method to convert the linear stability equations into a matrix eigenvalue, problem they obtained the neutral curves for various inclinations of the layer. They found that there is a maximum inclination above which the basic state is linearly stable. Thus, we may expect convection to arise in an unbounded inclined layer when the inclination satisfies $\left|\alpha \right|<31.49032°$, but the basic state is otherwise linearly stable.

A much more recent work by Wen and Chini [5] considered nonlinear convection at moderate values of the Darcy–Rayleigh number and presented solutions when the inclination is as high as $35°$, which is above the linear threshold given in [4]. Further work by Wen and Chini [6], which concentrated mainly on convection at very high values of the Darcy–Rayleigh number, used an energy method to show that nonlinear convection can occur at all inclinations as long as the Darcy–Rayleigh number exceeds $91.6$. Clearly, the source of these nonlinear solutions owes little to linearized stability theory, and this motivated Rees and Barletta [7] to investigate how nonlinear solutions may be generated as the inclination increases past $31.49032°$. It was suspected that it could be explained only by the presence of subcritical instabilities arising somewhere within the range of $0<\alpha <31.49032°$, and therefore, a weakly nonlinear theory was used to look for clues. It was found that the onset of convection is supercritical when $0\le \alpha <24.247627°$, but is otherwise subcritical. The analysis was confirmed by some strongly nonlinear finite difference computations for $24.247627°<\alpha <31.49032°$, where for any chosen value of $\alpha$ within that range, strongly nonlinear flows were found at subcritical values of the Darcy–Rayleigh number. Further work on this aspect of the layer analogue of the present cavity problem is in progress.

The works of Wen and Chini [5,6] were concerned with porous layers. At roughly the same time, Guerrero-Martínez et al. [8] used a finite difference method to determine the effect of inclination on the convective flows within a porous cavity as a function of the Darcy–Rayleigh number, the inclination, and the aspect ratio. They too discovered strongly nonlinear flows at $\alpha =45°$, and this fact alone forms the main motivation for the present study. So, one possible question to be answered is: is there a maximum inclination for which nonlinear convection exists, by which is meant something other than a buoyancy-induced single cell state? The paper by Guerrero-Martínez et al. [8] presents results in the form of the variation of the Nusselt number with the inclination, but the solution curves often just terminate. This phenomenon tends to be indicative of the unseen presence of unstable steady states that cannot be computed using unsteady methods. It is also indicative of fold bifurcations. Therefore, we have selected to use a spectral method, which reduces the governing steady-state equations to a system of nonlinear algebraic equations which are then solved using the Newton-Raphson method. The presence of fold bifurcations also dictated the need to use a pseudo-arclength continuation scheme.

It is important to note that a similar type of analysis was undertaken by Riley and Winters [9] using a continuation scheme, but where the basic solver involved a finite element discretization. They confined their analysis to a tilted unit cavity, and displayed their results in terms of increasing values of the Darcy–Rayleigh number.

In the following section, we present the governing equations and then we describe the numerical scheme in Section 3. Section 4 contains the result of our numerical validation exercise, together with a comparison with one case shown in [8] in order to demonstrate the utility of the present numerical scheme. Detailed remarks are also given about modal exchange mechanisms in terms of the direction of circulation of the convection cells. Section 5 forms our main results section where we show various examples of modal interaction for different cavity aspect ratios and values of Darcy–Rayleigh numbers. The general aim in that section is to acquire some intuitive understanding of what interactions are likely to happen in general. Finally, we present some conclusions in Section 6.

2. Governing Equations

We will consider nonlinear two-dimensional convection in a porous cavity that has aspect ratio, L, and is inclined at the angle $\alpha$ to the horizontal. It is assumed that the solid matrix and the fluid are in local thermal equilibrium, that Darcy’s law and the Boussinesq approximation apply, and that the fluid density is a linear function of the temperature. Subject to these assumptions, the dimensional governing equations are

(1)${\displaystyle \frac{\partial \overline{u}}{\partial \overline{x}}}+{\displaystyle \frac{\partial \overline{w}}{\partial \overline{z}}}=0,$

(2)$\overline{u}=-{\displaystyle \frac{K}{\mu }}{\displaystyle \frac{\partial \overline{p}}{\partial \overline{x}}}+{\displaystyle \frac{{\rho }_{f}g\beta K}{\mu }}(T-T_{\mathrm{c}})sin\alpha ,$

(3)$\overline{w}=-{\displaystyle \frac{K}{\mu }}{\displaystyle \frac{\partial \overline{p}}{\partial \overline{z}}}+{\displaystyle \frac{{\rho }_{f}g\beta K}{\mu }}(T-T_{\mathrm{c}})cos\alpha ,$

and

(4)$\sigma {\displaystyle \frac{\partial T}{\partial \overline{t}}}+\overline{u}{\displaystyle \frac{\partial T}{\partial \overline{x}}}+\overline{w}{\displaystyle \frac{\partial T}{\partial \overline{z}}}=\kappa \left({\displaystyle \frac{{\partial }^{2}T}{\partial {\overline{x}}^2}}+{\displaystyle \frac{{\partial }^{2}T}{\partial {\overline{z}}^2}}\right).$

The coordinates are $\overline{x}$ and $\overline{z}$, which are, respectively, the coordinate up the layer and the coordinate that is perpendicular to the bounding surfaces, while the corresponding Darcy velocities are $\overline{u}$ and $\overline{w}$. The temperature and the dynamic pressure are T and $\overline{p}$. Here, K is the permeability, $\mu$ the dynamic viscosity, $\beta$ the coefficient of cubical expansion, g the acceleration due to gravity, $\sigma$ the ratio of the heat capacity of the saturated porous medium to that of the saturating fluid, $T_{\mathrm{c}}$ the reference temperature, and ${\rho }_f$ the fluid density at the reference temperature. Finally, $\kappa$ is the ratio of the average thermal conductivity of the saturated medium and the product of the density and the specific heat of the fluid.

The boundary conditions are as follows,

(5)$\overline{x}=0,Lh:\overline{u}=0,{\displaystyle \frac{\partial T}{\partial x}}=0,$

and

(6)$\overline{z}=0:\overline{w}=0,T=T_h,\overline{z}=h:\overline{w}=0,T=T_c,$

where h is the thickness of the porous cavity, L is its aspect ratio, and $T_c<T_h$. We shall restrict our attention to two-dimensional steady flows, and therefore, we shall set $\overline{v}=0$ and eliminate all $\overline{y}$- and $\overline{t}$-derivatives.

The full system of equations and boundary conditions may now be rendered dimensionless by using the following scalings,

(7)$(x,z)={\displaystyle \frac{1}{h}}(\overline{x},\overline{z}),(u,w)={\displaystyle \frac{h}{\kappa }}(\overline{u},\overline{w}),p={\displaystyle \frac{K}{\kappa \mu }}\overline{p},\theta ={\displaystyle \frac{T-T_c}{T_h-T_c}}.$

The governing equations and boundary conditions become,

(8)${\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{\partial w}{\partial z}}=0,$

(9)$u=-{\displaystyle \frac{\partial p}{\partial x}}+\mathrm{Ra}\theta sin\alpha ,$

(10)$w=-{\displaystyle \frac{\partial p}{\partial z}}+\mathrm{Ra}\theta cos\alpha ,$

(11)$u{\displaystyle \frac{\partial \theta }{\partial x}}+w{\displaystyle \frac{\partial \theta }{\partial z}}={\displaystyle \frac{{\partial }^2\theta }{\partial x^2}}+{\displaystyle \frac{{\partial }^2\theta }{\partial z^2}},$

(12)$x=0,L:u=0,{\displaystyle \frac{\partial \theta }{\partial x}}=0,$

and

(13)$z=0:w=0,\theta =1,z=1:w=0,\theta =0.$

In the above, the Darcy–Rayleigh number, $\mathrm{Ra}$, is given by,

(14)$\mathrm{Ra}={\displaystyle \frac{{\rho }_{f}g\beta Kh(T_h-T_c)}{\mu \kappa }}.$

Given that the flow is now taken to be two-dimensional, we may define the streamfunction $\psi$, using

(15)$u=-{\displaystyle \frac{\partial \psi }{\partial z}},w={\displaystyle \frac{\partial \psi }{\partial x}},$

which automatically satisfies the equation of continuity, Equation (8), and hence, we obtain the system

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

(17)${\displaystyle \frac{{\partial }^2\theta }{\partial x^2}}+{\displaystyle \frac{{\partial }^2\theta }{\partial z^2}}={\displaystyle \frac{\partial \psi }{\partial x}}{\displaystyle \frac{\partial \theta }{\partial z}}-{\displaystyle \frac{\partial \psi }{\partial z}}{\displaystyle \frac{\partial \theta }{\partial x}}.$

Figure 1 shows the precise boundary conditions for $\psi$ and $\theta$ that need to be satisfied. In addition, the figure also depicts the isotherms for a weak free-convective circulation where the hot lower surface induces a buoyant flow up that surface, and the cold upper surface causes a corresponding flow down that surface.

3. Numerical Method

3.1. Spectral Method

We intend to use a spectral method in the form of a double half-range Fourier Sine Series to solve the governing equations. Given that the thermal boundary condition at $z=0$ is inhomogeneous, we may replace $\theta$ by $\theta +1-z$ in Equations (16) and (17) to render it homogeneous. This yields

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

(19)${\displaystyle \frac{{\partial }^2\theta }{\partial x^2}}+{\displaystyle \frac{{\partial }^2\theta }{\partial z^2}}={\displaystyle \frac{\partial \psi }{\partial x}}{\displaystyle \frac{\partial \theta }{\partial z}}-{\displaystyle \frac{\partial \psi }{\partial z}}{\displaystyle \frac{\partial \theta }{\partial x}}-{\displaystyle \frac{\partial \psi }{\partial x}}.$

We note that $(1-z)$ is the basic state when the layer is unbounded in the x-direction. The boundary conditions for $\theta$ are unchanged from what is displayed in Figure 1, except that $\theta =0$ on $z=0$.

Equations (18) and (19) may be solved using the following substitutions,

(20)${\displaystyle \psi =\sum _{i=1}^{\infty }\sum _{j=1}^{\infty }{A}_{i,j}sinikxsinj\pi z,}$

(21)$\theta =\sum _{j=1}^{\infty }{B}_{0,j}sinj\pi z+\sum _{i=1}^{\infty }\sum _{j=1}^{\infty }{B}_{i,j}cosikxsinj\pi z,$

where $k=\pi /L$ is the wavenumber that corresponds to the presence of one cell within the cavity, $0\le x\le L$.

The streamfunction Equation (18) becomes,

(22)$\begin{array}{ccc}\hfill {\displaystyle {\displaystyle \frac{{\partial }^2\psi }{\partial x^2}}+{\displaystyle \frac{{\partial }^2\psi }{\partial z^2}}}& =& {\displaystyle \mathrm{Ra}sin\alpha -\mathrm{Ra}cos\alpha \sum _{i=1}^{\infty }\sum _{j=1}^{\infty }ik{B}_{i,j}sinikxsinj\pi z}\hfill \\ & -& {\displaystyle \mathrm{Ra}sin\alpha \sum _{j=1}^{\infty }j\pi B_{0,j}cosj\pi z-\mathrm{Ra}sin\alpha \sum _{i=1}^{\infty }\sum _{j=1}^{\infty }j\pi B_{i,j}cosikxcosj\pi z.}\hfill \end{array}$

However, the first, third, and fourth terms on the right-hand side are not all zero on the four boundaries of the cavity, and therefore, they need to be expanded as a double half-range Fourier Sine Series. We obtain,

(23)$\begin{array}{ccc}\hfill \psi & =& {\displaystyle -\mathrm{Ra}sin\alpha \underset{i+j=\mathrm{odd}}{\sum _{i=1}^{\infty }\sum _{j=1}^{\infty }}\frac{16}{ij(i^2{k}^2+j^2{\pi }^2)k\pi L}sinikxsinj\pi z}\hfill \\ & & \\ & +& {\displaystyle \mathrm{Ra}cos\alpha \sum _{i=1}^{\infty }\sum _{j=1}^{\infty }\frac{ik{B}_{i,j}}{(i^2{k}^2+j^2{\pi }^2)}sinikxsinj\pi z}\hfill \\ & & \\ & +& {\displaystyle \mathrm{Ra}sin\alpha \underset{i_2=\mathrm{odd}}{\sum _{i_2=1}^{\infty }}\underset{j+j_2=\mathrm{odd}}{\sum _{j_2=1}^{\infty }\sum _{j=1}^{\infty }}\frac{16j_{2}j{B}_{0,j}}{i_2(j_2^2-j^2)(i_2^2{k}^2+j_2^2{\pi }^2)k\pi L}sin{i}_{2}kxsin{j}_2\pi z}\hfill \\ & & \\ & +& {\displaystyle \mathrm{Ra}sin\alpha \underset{i+i_2=\mathrm{odd}}{\sum _{i_2=1}^{\infty }\sum _{i=1}^{\infty }}\underset{j+j_2=\mathrm{odd}}{\sum _{j_2=1}^{\infty }\sum _{j=1}^{\infty }}\frac{16i_2{j}_{2}j{B}_{i,j}}{(i_2^2-i^2)(j_2^2-j^2)(i_2^2{k}^2+j_2^2{\pi }^2)k\pi L}sin{i}_{2}kxsin{j}_2\pi z.}\hfill \end{array}$

If we should know all of the $B_{i,j}$ coefficients for the temperature field, then this expression defines the values of $A_{i,j}$ in Equation (20).

The heat transport Equation (19), becomes

(24)$\begin{array}{ccc}\hfill 0& =& {\displaystyle -\sum _{j=1}^{\infty }{j}^2{\pi }^2{B}_{0,j}sinj\pi z+\sum _{i=1}^{\infty }\sum _{j=1}^{\infty }ik{A}_{i,j}cosikxsinj\pi z}\hfill \\ & & \\ & -& {\displaystyle \sum _{i=1}^{\infty }\sum _{j=1}^{\infty }(i^2{k}^2+j^2{\pi }^2)B_{i,j}cosikxsinj\pi z}\hfill \\ & & \\ & +& {\displaystyle \sum _{i_1=1}^{\infty }\sum _{j_1=1}^{\infty }\sum _{i_2=1}^{\infty }\sum _{j_2=1}^{\infty }\frac{1}{4}\left(i_2{j}_{1}k\pi \right)A_{i_1,j_1}{B}_{i_2,j_2}[cos(i_1+i_2)kx-cos(i_1-i_2)kx]}\hfill \\ & & \times \left[sin(j_1+j_2)\pi z-sin(j_1-j_2)\pi z\right]\hfill \\ & & \\ & -& {\displaystyle \sum _{i_1=1}^{\infty }\sum _{j_1=1}^{\infty }\sum _{i_2=1}^{\infty }\sum _{j_2=1}^{\infty }\frac{1}{4}\left(i_1{j}_{2}k\pi \right)A_{i_1,j_1}{B}_{i_2,j_2}\left[cos(i_1+i_2)kx+cos(i_1-i_2)kx\right]}\hfill \\ & & \times \left[sin(j_1+j_2)\pi z+sin(j_1-j_2)\pi z\right]\hfill \\ & & \\ & -& {\displaystyle \sum _{i_1=1}^{\infty }\sum _{j_1=1}^{\infty }\sum _{j_2=1}^{\infty }\frac{1}{2}\left(i_1{j}_{2}k\pi \right)A_{i_1,j_1}{B}_{0,j_2}cos{i}_{1}kx\left[sin(j_1+j_2)\pi z+sin(j_1-j_2)\pi z\right].}\hfill \end{array}$

This system may be regarded as being a set of nonlinear algebraic equations for the $B_{i,j}$ coefficents where, again, we note that each $A_{i,j}$ term that appears above may be written in terms of the values of the $B_{i,j}$ terms. All of the summations were truncated so that the summation counters satisfy $0\le i\le \mathrm{NX}$ and $1\le j\le \mathrm{NZ}$. This means that we will need to solve for $(\mathrm{NX}+1)\times \mathrm{NZ}$ values of the $B_{ij}$ coefficients. Finally, we note that the case where $\mathrm{NX}=1$ and $\mathrm{NZ}=2$ represents the minimal truncation that still retains some nonlinearity; generally, this will be quite inaccurate.

3.2. Implementation

We employed a user-written multi-dimensional Newton-Raphson code, the iteration scheme for which requires the use of the Jacobian matrix. Such matrices are usually encoded directly into the program, but the presence of the quadruple summations in Equations (23) and (24) means that the direct encoding of the Jacobian is impractical. Therefore, we used a numerical differentiation routine as a practical alternative. The Jacobian matrix was then created by varying in turn each $B_{i,j}$-value by ${10}^{-8}$ (using double precision Fortran90), and this allowed us to evaluate all of the required derivatives.

One of the chief reasons for selecting this direct method of solution was to enable us to determine solutions that are unstable in practice, and that are therefore unable to be computed using unsteady methods. It is also clear from the shapes of the solution curves given in [8] that fold bifurcations arise, and therefore, it was essential to include a curve-tracking facility.

We used a pseudo arc-length methodology that involved the Nusselt number, and the inclination where both $\mathrm{Ra}$ and L take fixed values. We define the Nusselt number as being the mean rate of heat transfer per unit length along the lower surface:

(25)$\mathrm{Nu}=-{\displaystyle \frac{1}{L}}{\int }_0^L{\displaystyle \frac{\partial \theta }{\partial z}}dx=-\sum _{j=1}^{\mathrm{NZ}}j\pi B_{0,j},$

where the final summation was obtained by integrating Equation (11). If the values of $\delta \mathrm{Nu}$ and $\delta \alpha$ are the maximum allowable absolute changes in $\mathrm{Nu}$ and $\alpha$ along a solution curve, and where the superscript, m, denotes the $m^{\mathrm{th}}$ point along a solution curve, then the additional equation that we need to solve is,

(26)${\left({\displaystyle \frac{{\mathrm{Nu}}^{m+1}-{\mathrm{Nu}}^m}{\delta \mathrm{Nu}}}\right)}^2+{\left({\displaystyle \frac{{\alpha }^{m+1}-{\alpha }^m}{\delta \alpha }}\right)}^2=1.$

In practice, if the solutions are already known at points m and $m-1$, then the linear extrapolation of the values of $\mathrm{Nu}$ and $\alpha$ form the initial iterate for the required values at point $m+1$.

While it is quite standard to devise numerical validation tests for one’s codes, it is of particular importance here because of the complexity of the present algorithm. Therefore, we selected to employ a finite difference scheme for the purpose. Second-order accurate central differences in space were used in a time-stepping code that was run to steady-state with more than six decimal places of accuracy. We chose to check the case of $L=1$, $\mathrm{Ra}=100$, and $\alpha =20°$ because this selection of parameters means that every term in the above summations was used. The value of the Darcy–Rayleigh number is roughly $2.5$ times the critical value for the onset of convection for a single cell occupying a horizontal cavity, and it is therefore quite strongly nonlinear. The steady-state obtained in this way was then processed via the computation of the equivalent values of $B_{1,1}$ and $B_{0,2}$ using the trapezium rule to evaluate the integrals,

(27)$B_{1,1}\simeq {\displaystyle \frac{4}{L}}{\int }_0^L{\int }_0^1\theta (x,z)sinkxsin\pi zdzdx,B_{0,2}\simeq {\displaystyle \frac{4}{L}}{\int }_0^L{\int }_0^1\theta (x,z)sin2\pi zdzdx,$

where $k=\pi /L$. A selection of grids was used, and Richardson’s extrapolation was employed recursively first to increase the order of accuracy from second to fourth order, and then secondly, to sixth order. The result of this process yielded the data given in Table 1.

The convergence of both numerical methods as the resolutions increase is quite evident. The nested application of Richardson’s extrapolation for the finite difference method allowed us to be able to obtain values of $B_{1,1}$ and $B_{0,2}$ that are accurate to six decimal places, and these values agree perfectly with the $17\times 18$ spectral results. Therefore, we conclude that both methods have been encoded correctly.

We also see that the $9\times 10$ case yields values of $B_{1,1}$ and $B_{0,2}$ that have relative errors of $2\times {10}^{-4}$ and $6\times {10}^{-5}$, respectively. Given the presence of the quadruple summations in Equations (23) and (24), the number of multiplications per Newton-Raphson iteration increases roughly as the fourth power of $\mathrm{NX}\times \mathrm{NZ}$, and therefore, we considered this resolution to be a good compromise between the accuracy and speed of computation. For other values of the aspect ratio L, we chose to use the spectral resolution $9L\times 10$ when L is an integer, and to use an $\mathrm{NX}$ value that is close to $9L$ otherwise.

4. Comparison with Guerrero-Martínez et al. [8] and Some Flow Patterns

It also seemed apposite to undertake a final check of the present encodings by comparing our results with those of Guerrero-Martínez et al. [8] in order to show the utility of the pseudo-arclength continuation in the present context. In what follows, the value n corresponds to the number of distinct cells within the cavity. Our mathematical definition of n is one more than the number of sign changes in $\psi$ at $z=\frac{1}{2}$ within the range of $0<x<L$.

Figure 2 displays the variation of the Nusselt number with inclination for the case of $\mathrm{Ra}=100$ and $L=3$. Figure 2a shows the finite difference computations of [8] where their definition of the Nusselt number, which is ${\mathrm{Nu}}_{\mathrm{GM}}=L\mathrm{Nu}$, varies with $\alpha$. The separate (and differently coloured) curves correspond to different numbers of cells within the cavity ($n=1$, $n=3$ and $n=4$), while the vertical dotted lines correspond to when their converged solution for one value of $\alpha$ has jumped suddenly to a different convection pattern as $\alpha$ has been increased or decreased incrementally.

Figure 2b shows the solution curves that were computed using the present spectral method with pseudo-arclength continuation. Those curves that are shown in Figure 2a are reproduced in Figure 2b with no discernible difference, and therefore, our computations agree with those in [8]. However, it is evident that Figure 2b also displays a large number of extra curves, all of which are fully converged steady-state solutions of the governing equations. These extra curves are either very difficult to be computed with an unsteady solver, or (which is much more likely to true) they represent solutions that are unstable.

The uppermost curve corresponds to a three-cell convection pattern which appears to approach a fold bifurcation as $\alpha$ increases because the slope of the Nusselt number curve becomes very large. In Figure 2a, the solution then drops down to the red curve, which is a one-cell pattern. In Figure 2b the corresponding curve turns back towards smaller inclinations. Finally, the curve undergoes a second fold bifurcation to the one-cell pattern that is characteristic of flows at high inclinations.

For this choice of parameters, Figure 3 shows the streamlines and isotherms corresponding to the four representative examples that correspond to the black disks in Figure 2b. Figure 3a,b represent three-cell flows for $\alpha =5°$ and $\alpha =25°$. In both cases, the outer cells are stronger than the middle one because their anticlockwise circulations are aided by the effect of buoyancy forces along the boundaries; we will call these natural circulations. By contrast, the inner cell is weaker because its circulation is inhibited by the action of buoyancy forces along the boundaries; this is an anti-natural circulation. This increasing inhibition of the middle cell as $\alpha$ increases not only leads to weaker circulations, but also to a decreased width.

As the $n=3$ curve is traversed, it passes a fold bifurcation when $\alpha \simeq 32.06°$ descends, transitions from representing a three-cell flow to representing a one-cell flow (see the change in colour), and encounters a second fold bifurcation at $\alpha \simeq 23.97$. Subsequently, the curve continues to represent a single cell until and beyond when the layer becomes vertical.

The process of transition from a three-cell state to a one-cell state may be seen by comparing Figure 3b–d. The middle cell seen in Figure 3b is then squeezed out of existence, leaving behind the last vestiges of the naturally circulating cells, as seen in Figure 3c. In this case, the streamfunction is single-signed within the computational domain, and therefore, this is regarded as a single cell. The precise point along the $\mathrm{Nu}$-curve where three cells become one is again marked by the change in the colour. Finally, at larger inclinations, all of the fluid circulates about the center of the cavity, as shown in Figure 3d.

This type of transition between patterns with odd values of n is quite ubiquitous, and it is also possible to have a similar transition between two even values of n.

Figure 4 shows the two $n=3$ cases, which are shown as yellow disks in Figure 2. The inclination is much shallower than those used in Figure 3, and therefore, the small difference between the streamlines may only be seen upon close inspection, but the underlying isotherm patterns are very different from one another. In Figure 4a, the outer cells display natural circulations, while in Figure 4b, they display anti-natural circulations. Therefore, it is to be expected that the former flow has a larger Nusselt number than the latter.

Figure 5 is an example of a transition between a pattern with an odd number of cells and a pattern with an even number. Again, the parameters are the same: $\mathrm{Ra}=100$, $L=3$, and $\alpha =5°$, but they correspond to the green disks in Figure 2. When one has four cells then one of the outer cells will be anti-natural; in Figure 5a, which corresponds to the uppermost green disk in Figure 2, it is the right-hand cell that plays this role. Indeed, it is quite clearly narrower and weaker than the remaining cells.

As the $n=4$ curve is traversed towards where the middle green disk is, the right-hand cell continues to weaken and it is now barely present in Figure 5b.

When one continues along the $n=4$ curve, one eventually passes through $\alpha =0$ and towards negative values. This will bring one to the location of $\alpha =-5°$, which is the mirror image of the lowest green disk in Figure 2c. The streamline and isotherm pattern for that point will be precisely the reflection of Figure 5c about $x=L/2$. However, the pattern we see represents $\alpha =5°$, and this point is very close to where the $n=4$ solution bifurcates from the middle branch of the $n=3$ (blue) curve. The bifurcation point itself marks when the former $n=4$ solution has finally removed the right-hand anti-natural circulation. This evolution provides a mechanism whereby one cell is either lost or gained at a bifurcation point.

5. Results Using the Spectral Method

It is clearly impossible to present a comprehensive account of how the pattern of convection changes as a function of $\mathrm{Ra}$, L, and $\alpha$, despite this being just a three-parameter set. Even with the results shown in Figure 2, it is also clear that there exist many different solutions corresponding to different numbers of cells, even when two of the three parameters have been fixed. That said, it is possible to be comprehensive when L is relatively small and $\mathrm{Ra}$ is not too large. Therefore, we shall present a selection of different cases with the overall aim of attempting to acquire a good intuition about how this cavity problem behaves in general. We shall split the presentation into subsections where the different aspect ratios are considered in each.

Before that, though, it is good to recall critical values of $\mathrm{Ra}$ as a function of L for when the cavity is horizontal. It is well-known that the critical Rayleigh number is

(28)$\mathrm{Ra}={\displaystyle \frac{{(k^2+{\pi }^2)}^2}{k^2}}$

for the horizontal porous layer [10], where k is the wave number. For an unbounded layer, k may take any real value, but for a cavity of width L, the wavenumber takes the discrete values $k=n\pi /L$, where n is the number of cells within the cavity. Therefore, Equation (28) translates into

(29)$\mathrm{Ra}={\displaystyle \frac{{(L^2+n^2)}^2}{n^2{L}^2}}{\pi }^2.$

In all of the subsequent figures, the solution curves will be coloured to indicate the parity of the flow, i.e., the number of cells that exist within the cavity:

(30)$\begin{array}{cc}n=1& n=2\\ n=3& n=4\\ n=5& n=6\\ n=7& n=8\end{array}$

So, modes with an odd parity have been assigned rustic colours and modes that have an even parity when $\alpha =0$ are formed as dashed lines. We note that Figure 2 in this paper has adopted the colouring scheme used in [8].

5.1. Solution Curves for $L=1$

This is the unit aspect ratio cavity. Given that the critical wave number for the unbounded horizontal layer is $k=\pi$, which corresponds to a length of 2 or, equivalently, a unit cell width, we would expect a single-cell flow to dominate the nonlinear dynamics for the unit cavity. If we define ${\mathrm{Ra}}^{\left(n\right)}$ to be the critical Darcy–Rayleigh number for an n-cell flow when the cavity is horizontal, then we may use Equation (29) to provide the following values:

(31)$\begin{array}{cc}& {\mathrm{Ra}}^{\left(1\right)}=4{\pi }^2=39.478,{\mathrm{Ra}}^{\left(2\right)}=6.25{\pi }^2=61.685,\\ & {\mathrm{Ra}}^{\left(3\right)}=\frac{100}{9}{\pi }^2=109.662,{\mathrm{Ra}}^{\left(4\right)}=\frac{289}{16}{\pi }^2=178.270,\end{array}$

where the numerical values are correct to three decimal places. These form the context for the following figures.

Figure 6 shows the evolution as $\mathrm{Ra}$ increases in the single-cell $n=1$ curves through a whole rotation of the cavity through 360 degrees. When $\alpha =\pm 180°$, the heating is from above and there is no flow at all, and hence, $\mathrm{Nu}=1$ which corresponds to pure conduction. At other inclinations, a flow is generated and this provides an increased Nusselt number.

When $\alpha$ increases from zero, a relatively weak buoyancy drives a slow flow and the Nusselt number increases. Eventually, the Nusselt number reaches a maximum and decreases back to $\mathrm{Nu}=1$ when $\alpha =180°$. All of the curves shown in Figure 6a correspond to natural circulations. However, once $\mathrm{Ra}$ exceeds ${\mathrm{Ra}}^{\left(1\right)}=4{\pi }^2$, there exists the possibility of having a thermoconvective instability and hence, a strong convection when $\alpha =0$. We see this in Figure 6b, where each of the solution curves exhibits a loop. It is perhaps easiest to focus on the $\mathrm{Ra}=100$ curve, where as $\alpha$ decreases from $180°$, we eventually obtain a non-unit value of $\mathrm{Nu}$ when $\alpha =0$, and this is a strongly nonlinear solution. Within this range of inclinations, the convection is a natural circulation, but once $\alpha$ becomes negative along the same curve, then the corresponding convection cell consists of an anti-natural circulation where buoyancy forces along the the boundaries again serve to inhibit the convective instability. Eventually, the curve passes through $\alpha =0$ where $\mathrm{Nu}=1$ and there is no flow. Everything that happens thereafter is the mirror image of what has already been described.

Figure 7 shows some solution curves where the values of $\mathrm{Ra}$ have been chosen to show the evolution of the system as $\mathrm{Ra}$ increases. When $\mathrm{Ra}=65$, we see immediately that there is a short solution branch corresponding to two cells, $n=2$. This is consistent with the fact that $\mathrm{Ra}>{\mathrm{Ra}}^{\left(2\right)}$; see Equation (31). This particular value of $\mathrm{Ra}$ is not strongly supercritical, and therefore, the $n=2$ flow is weak and $\mathrm{Nu}$ is only slightly larger than 1. The $n=2$ branch terminates at a bifurcation from the $n=1$ curve, which is of the same nature as was shown in Figure 5; in this case, the two-cell flow that exists at $\alpha =0$ gradually loses the cell with the anti-natural circulation, and hence, an $n=1$ state is eventually obtained as $\alpha$ increases. We also note that the red $n=2$ curve is symmetric about the vertical axis.

As $\mathrm{Ra}$ increases further, the $n=2$ curve develops a cusp at $\alpha =0$ (see $\mathrm{Ra}=80$) and then a loop (see $\mathrm{Ra}=100$ and beyond). In the meantime, a small section of curve close to $\alpha =0$ and $\mathrm{Nu}=1$ is being identified as being an $n=3$ solution and is coloured in green. At $\alpha =0$, a cusp may then be seen when $\mathrm{Ra}=110$, and this marks the onset of the convection of a three-cell mode, since $\mathrm{Ra}$ is now almost exactly the value of ${\mathrm{Ra}}^{\left(3\right)}$; see Equation (31). When $\mathrm{Ra}=120$, the $n=3$ cusp has now developed into a loop, reflecting the fact that we now have a three-cell solution of moderate amplitude.

At $\mathrm{Ra}=150$, the $n=3$ curve has become more complicated, while the Nusselt number for $n=2$ has increased to such an extent that it is only just smaller than that for the $n=1$ solution.

In all of the above, we have merely computed the solution curves, and the present code/method, despite its complexity, is incapable of providing information on the stability of the computed solutions. It is also worth mentioning that the solution curves become increasingly intricate and difficult to obtain as $\mathrm{Ra}$ increases further. Thus far, only the $n=1,2$ and 3 solutions have appeared, but we would need to consider larger values of $\mathrm{Ra}$ in order to see some $n=4$ solutions, since $R{a}^{\left(4\right)}=178.270$ when $L=1$.

5.2. Solution Curves for $L=2$

For $L=2$, we may use Equation (29) again to provide the following critical values for $\mathrm{Ra}$ for different values of n:

(32)$\begin{array}{ccc}{\mathrm{Ra}}^{\left(1\right)}=6.25{\pi }^2=61.685,\hfill & {\mathrm{Ra}}^{\left(2\right)}=4{\pi }^2=39.478,\hfill & {\mathrm{Ra}}^{\left(3\right)}=\frac{169}{36}{\pi }^2=46.332,\hfill \\ {\mathrm{Ra}}^{\left(4\right)}=6.25{\pi }^2=61.685,\hfill & {\mathrm{Ra}}^{\left(5\right)}=\frac{841}{100}{\pi }^2=83.003,\hfill & {\mathrm{Ra}}^{\left(6\right)}=\frac{100}{9}{\pi }^2=109.662.\hfill \end{array}$

Naturally, for a cavity of length, $L=2$, we would expect the primary mode of instability to be formed from two cells. Clearly, an $n=3$ pattern is expected to appear at a slightly higher Darcy–Rayleigh number, and then both the $n=1$ and $n=4$ patterns arise subsequently at the same value of $\mathrm{Ra}$.

Figure 8 shows some subcritical cases for $\mathrm{Ra}=10$, 20, 30, and 40. Qualitatively, these single-cell solutions have exactly the same variation with $\alpha$ as was seen in Figure 6a for $L=1$.

Figure 9 shows the solution curves for six chosen values of $\mathrm{Ra}$. When $\mathrm{Ra}=40$, which is just above the critical value for the $n=2$ pattern, we see an $n=2$ curve within a small range of inclinations centered on $\alpha =0$. This weak nonlinear flow bifurcates from the unicellular state when $\left|\alpha \right|\simeq 1.2°$, and the $n=1$ curve then forms the unique solution for larger inclinations.

At the slightly higher $\mathrm{Ra}=50$, which is above the critical value of $\mathrm{Ra}$ for $n=2$ and $n=3$, we see that the solution curves have already evolved quite substantially. The black $n=1$ curve that is present for large inclinations eventually transforms into the green $n=3$ labeling as $\alpha$ decreases, and a loop is formed near to $\alpha =0$. So, when the layer is horizontal, there are only two nonlinear solutions, one with two cells and one with three. We also note that the red $n=2$ curve now bifurcates from a green $n=3$ curve.

As $\mathrm{Ra}$ increases still further, the evolution of the solution curves becomes increasingly complicated. Thus, as $\mathrm{Ra}$ increases from 60 through ${\mathrm{Ra}}^{\left(1\right)}={\mathrm{Ra}}^{\left(4\right)}=61.685$, the black $n=1$ curve centered on $\alpha =0$ develops a cusp that then transforms into a loop with a nonlinear solution at $\alpha =0$. At the same values, $\mathrm{Ra}=61.685$, a blue $n=4$ solution first appears, and a weak nonlinear solution may be seen when $\mathrm{Ra}=70$. These curves continue to evolve and develop further loops as $\mathrm{Ra}$ continues to increase, and when $\mathrm{Ra}=100$, we also see the appearance of the orange $n=5$ solution, given that $\mathrm{Ra}$ is now above the value of ${\mathrm{Ra}}^{\left(5\right)}$ given in Equation (32).

5.3. Solution Curves for $L=3$

Figure 10 is the $L=3$ counterpart to Figure 6a and Figure 8, and it has been plotted using the same vertical scale for the sake of comparison. The shapes of these curves offer no surprises, except that the $\mathrm{Ra}=40$ curve is shaded in green when the layer is close to being horizontal. For such a cavity, the $n=3$ mode has $4{\pi }^2$ as its critical Darcy–Rayleigh number when $\alpha =0$, and therefore, the value of $\mathrm{Ra}=40$ is slightly supercritical. Clearly, the $n=1$ solution evolves into a three-cell solution as $\alpha$ tends towards zero because of the presence of thermoconvective instability.

For $L=3$ we may again use Equation (29) to determine the following critical values of $\mathrm{Ra}$ when the layer is horizontal:

(33)$\begin{array}{ccc}{\mathrm{Ra}}^{\left(1\right)}=\frac{100}{9}{\pi }^2=101.662,\hfill & {\mathrm{Ra}}^{\left(2\right)}=\frac{169}{36}{\pi }^2=46.332,\hfill & {\mathrm{Ra}}^{\left(3\right)}=4{\pi }^2=39.478,\hfill \\ {\mathrm{Ra}}^{\left(4\right)}=\frac{625}{144}{\pi }^2=42.837,\hfill & {\mathrm{Ra}}^{\left(5\right)}=\frac{1156}{225}{\pi }^2=50.7078,\hfill & {\mathrm{Ra}}^{\left(6\right)}=6.25{\pi }^2=61.685,\hfill \\ {\mathrm{Ra}}^{\left(7\right)}=\frac{3364}{441}{\pi }^2=76.281,\hfill & {\mathrm{Ra}}^{\left(8\right)}=\frac{5329}{576}{\pi }^2=91.311,\hfill & {\mathrm{Ra}}^{\left(9\right)}=\frac{100}{9}{\pi }^2=101.662.\hfill \end{array}$

Figure 11 shows the detailed solution curves for $\mathrm{Ra}=40$, 50, 70, and 100. The subfigure corresponding to $\mathrm{Ra}=40$ is a close-up of the uppermost curve in Figure 10, but there is a very narrow green loop centered on $\alpha =0$ that shows that $\mathrm{Ra}$ is supercritical and that a nonlinear $n=3$ solution exists, albeit a fairly weak one.

A small increase in $\mathrm{Ra}$ to 50 brings with it the potential for three patterns to exist: $n=2$, $n=3$, and $n=4$; see Equation (33). However, first, we note the $n=1$ curve, which is the sole solution for large inclinations, transforms into an $n=3$ curve as $\alpha$ decreases, and as $\alpha$ decreases still further towards zero, a final transformation to a five-cell state takes place. The orange colour is a herald of the fact that the critical value, ${\mathrm{Ra}}^{\left(5\right)}$, is only very slightly higher than 50.

We do indeed see the presence of the two-cell and four-cell solutions, as predicted, where the latter is stronger when $\alpha =0$, and this is possibly due to the fact that it has the lower critical Darcy–Rayleigh number. Both the $n=2$ and $n=4$ solutions bifurcate from the central loop, which transitions between being $n=3$ and an $n=5$ curve. A more complicated scenario may be seen for $\mathrm{Ra}=70$. The $n=2$ and $n=4$ solutions are now much stronger, and the $n=5$ solution, hinted at when $\mathrm{Ra}=50$, is also well into the nonlinear regime. It is of interest to note that the straightforward bifurcation of the $n=4$ (blue) solution from the $n=3$ (green) solution that we saw for $\mathrm{Ra}=50$ now assumes a more complicated form with the presence of a very tight loop. Indeed the presence of tight loops turns out to be quite ubiquitous when either $\mathrm{Ra}$ or L increases. We also see the relatively weak presence of an $n=6$ (purple) pattern which bifurcates from the $n=1$ branch at the bottom of the subfigure; its presence somewhere was expected, given the value of ${\mathrm{Ra}}^{\left(6\right)}$ presented in Equation (33).

When $\mathrm{Ra}=100$, we obtain the case that we considered in Figure 2. Given Equation (33), we could have predicted the presence of modes $n=2$ through to $n=8$, and the newly added modes $n=7$ (brown) and $n=8$ (grey) are quite weak. Mode $n=8$ bifurcates from an $n=1$ branch in a manner that is now quite familiar. On the other hand, the way in which the $n=7$ mode appears is quite unusual. That said, one of main characteristics of this subfigure is the number of loops that occur. All three solutions with even parity exhibit a loop centered about $\alpha =0$, while another arises close to $\alpha =6°$ and $\mathrm{Nu}=1.35$. There are also two very sharp changes in direction that are difficult to follow in the present code, despite it being designed to do so. These are at $(\alpha ,\mathrm{Nu})=(2.4°,2.1)$ and $(\alpha ,\mathrm{Nu})=(14°,1.95)$.

5.4. The Effect of Varying L When $\mathrm{Ra}=65$

Thus far, we have chosen to consider three different aspect ratios of the cavity, and to concentrate on how the solution curves vary as $\mathrm{Ra}$ increases. All of these aspect ratios correspond to integer values of L, and therefore, the onset of convection arises when $\mathrm{Ra}=4{\pi }^2$ when the cavity is horizontal, and $n=L$ forms the most unstable mode of convection. It is well-known, and easily derivable from Equation (29), that the most unstable mode has n cells when L satisfies

(34)$\sqrt{n(n-1)}<L<\sqrt{n(n+1)}.$

Therefore, this subsection is devoted to illustrating how inclination and nonlinearity affect this conclusion when L takes noninteger values. We shall allow L to increase from $L=1.2$ to $L=2.8$ in steps of $0.2$. This range spans the full transition from an odd parity pattern through to an even parity, and on to the next odd parity pattern. We have chosen to consider just the one case, $\mathrm{Ra}=65$, which represents a moderately supercritical Darcy–Rayleigh number. Nine different cases may be seen in Figure 12, and all of the subfigures use the same vertical scale.

When $L=1.2$, we are well within the range of values of L for which we expect one cell when $\alpha =0$, but an $n=2$ solution bifurcates from the $n=1$ branch and is quite weak. The large-$\alpha$ solution curves are exactly like those seen in Figure 6b, for example. When $L=1.4$, we are just below $L=\sqrt{2}$, which is when the $n=1$ and $n=2$ modes have the same critical Darcy–Rayleigh number. Nevertheless, the two-cell solution is slightly stronger, at least in terms of the Nusselt number, over a narrow range of inclinations. A green colouring at the bottom of the subfigure corresponds once more to the imminent appearance of an $n=3$ solution for which the critical Darcy–Rayleigh number is $\mathrm{Ra}=67.208$.

The green $n=3$ pattern may be seen in the $L=1.6$, $L=1.8$, and $L=2.0$ subfigures, and these flows increase in strength as L increases within that range. However, the $n=2$ pattern now dominates all other modes over an increasing range of inclinations centered on $\alpha =0$. When $L=2$ this range is roughly $-11°<\alpha <11°$. In all of the subfigures, up to and including $L=2$, the $n=2$ solution itself bifurcates from the $n=1$ solution.

We note that the onset of an $n=4$ pattern when $L=2$ is when $\mathrm{Ra}=61.685$, which is less than 65, and indeed, a close examination of the lowest part of the $n=1$ curve (i.e., under quite strong magnification) shows the initial appearance of the blue $n=4$ solution. We have examined this in detail, and it displays the same qualitative transition via the red $n=2$ state to a bifurcation from the black $n=1$ as the $L=2.2$ subfigure shows the red curve transitions via the green $n=3$ state to a bifurcation from the (continuous) $n=3$ curve.

As L increases further via $L=2.2$, we see a strengthening of the $n=3$ flow relative to the $n=2$ flow. The $n=4$ pattern also continues to grow, which is opposite to that of the $n=2$ flow. These blue curves are incomplete and we have been unable to compute them. However, Figure 11 ($\mathrm{Ra}=70$, a close case) shows what we believe to be the shape which will be adopted here.

Our conclusion based on this sample of computations is (i) that when L is close to being an odd integer, m, then the $n=m$ pattern dominates for quite a wide range of inclinations; and (ii) that when L is close to being an even integer, m, then the $n=m$ pattern dominates, but over a relatively narrow range of inclinations; that when L approaches an odd integer value, m, then the strength of the $n=m-1$ pattern weakens while that of the $n=m+1$ pattern increases (e.g., comparing the behaviours of the red and blue curves between $L=2.6$ and $L=2.8$). We presume that a similar observation may be made concerning when L approaches an even integer value.

5.5. Some Solutions for Yet Larger Aspect Ratios

It is perhaps no surprise that the computational difficulties increase when either $\mathrm{Ra}$ or L are increased. Larger values of the Darcy–Rayleigh number allows for further modal patterns to fit within the cavity, and for an accompanying increase in the intricacy for how those different patterns interact. While the same observations are also true when L increases, not only does the number of computations increase, but the length of each computation itself increases as the third power of $\mathrm{NX}$, and also the number of spectral modes used in the x-direction.

In this subsection, we shall confine ourselves to $\mathrm{Ra}=100$ and to the following cavity aspect ratios: $L=4$, 5, 6, and 7. The aim will be to determine a qualitative understanding of the modal exchanges that take place in a shallow cavity as the aspect ratio increases. Figure 13 depicts the corresponding solution curves, although none of the individual subfigures is complete. Our initial numerical experimentation showed that an astonishing degree of complexity is obtained in the region fairly closely centered on $\alpha =0$, and for smaller values of $\mathrm{Nu}$. This has already been seen in Figure for $\mathrm{Ra}=100$. Therefore, Figure 13 has retained only those solution curves that exist for the larger inclinations (usually those with odd parity) and a selection of even parity solutions where n and L are close to one another.

If attention is confined to the even-parity modes shown in Figure 13, then we see that there is a natural tendency for curves corresponding to a fixed value of n to rise as L increases towards the value of n, and to descend as L increases further. The presence of even parity cells is generally confined to within approximately $12°$ of the horizontal, although there is a tendency for that limit to increase slightly as L increases. It is also interesting to note that the $n=6$ mode transfers more heat than any other mode when $L=5$ and $\alpha$ is close to zero, and the same may also be said for the $n=8$ mode when $L=7$.

If we now confine attention to the larger inclinations where only those solutions with an odd parity exist, then we see the presence of many fold bifurcations. Generally, the right-most structure takes an S-shape and it connects a pair of fold bifurcations. Other structures at smaller inclinations take a distinctive downward-pointing tongue-like shape, which is illustrated well for the $L=5$ case.

Interactions and couplings between the solution curves as L increases are many and of different kinds, and it is difficult to generalize or even to describe them easily. For example, the green $n=3$ loop for $L=6$ has become a cusp when $L=7$, and although we do not display it here, it becomes a smooth fold when $L=8$. Further, the $n=1$ curve, which forms a unique nonlinear solution at high inclinations, eventually transforms into the green $n=3$ pattern when $L=5$, but it transforms into the orange $n=5$ pattern when $L=7$. In this regard, the transitional case, $L=6$, is not clear in Figure 13, although the numerical data show that the transition from $n=1$ is to $n=5$, as it is for $L=5$. We illustrate how this happens in Figure 14 where we have used the aspect ratios, $L=5.8$ and $L=6.2$. For both the cases presented, we note that each consists of just one continuous solution curve. The computation itself is displayed as a thick line, while its mirror image is displayed using a thin line in order to clarify what happens to the solution trajectory when L increases from $5.8$ to $6.2$. There is clearly more than one qualitatve difference between the two trajectories.

5.6. A Maximum Inclination for Convective Instability?

Finally, and given that one of the original motivations for this study was the observation by Guerrero-Martínez [8] that buoyant instability is manifested in a porous cavity at an inclination that is well above that given by the linearized stability theory for the unconfined layer, it is worth presenting some initial data on this aspect.

The preceding figures demonstrate that when $L=3$ or larger, then the right-most solution curve displaying the $n=1$ to $n=3$ modal transition is the one that provides the maximum inclination for the existence of buoyancy-induced instability. This inclination may be obtained by the simple expedient of fitting a suitable parabola to the ‘nose’ of such a curve and finding its extremum. Table 2 shows the result of this process, and also of the second ‘nose’, the one with the $n=1$ to $n=5$ transition; these data are also displayed in Figure 15.

A brief glance at Figure 15a shows that there is an increasing freedom for strongly nonlinear convection to exist as the aspect ratio increases, but it also suggests that there could be a limiting value of the maximum inclination when L is asymptotically large. When plotted against $L^{-1/2}$ in Figure 15b, there is a suspicion that $\alpha$ could be a linear function of $L^{-1/2}$, and that the limiting inclination when $\mathrm{Ra}=100$ could therefore be very close to $\alpha =54°$ for the green disks. Intriguingly, the same limiting value appears to apply for the orange disks. At present, we do not know whether this is a significant observation or not, but we shall deem it to be worthy of further investigation later.

An alternative viewpoint is provided by Figure 16, which shows solution curves for $L=7$ and for values of $\mathrm{Ra}$ in steps of 10 between $\mathrm{Ra}=80$ and $\mathrm{Ra}=150$. Here, we show only a section of the right-most solution curve in each case, i.e., that part which is close to the maximum inclination for cellular convection. The locus of the respective maximum inclinations are given by the dashed red line, and it appears that the variation of ${\alpha }_{\max }\left(\mathrm{Ra}\right)$ is not too far from being linear. Clearly, it is impossible to predict whether the maximum inclination will continue to increase in this way as $\mathrm{Ra}$ increases further, or whether there will be an maximum inclination that is independent of $\mathrm{Ra}$ for a fixed value of L, or whether there is a global maximum inclination that is independent of both $\mathrm{Ra}$ and L. Again, we intend that this aspect will be the subject of subsequent research, but we can certainly conclude that nonlinear convection is possible when $\alpha$ is as large as $60°$ from the horizontal.

6. Conclusions

This has been a three-parameter system, something for which a fairly comprehensive set of results may usually be obtained, and from which a definitive understanding of the problem is often acquired. It is clear that these desirable outcomes are not possible for the inclined Darcy-Bénard cavity. Part of the reason for this is the possibility of having more than one potential convection pattern. Another reason is the potential for different solution trajectories to interact and to change their morphologies as a governing parameter changes. We have uncovered certain principles for these changes and some features that are ubiquitous:

• Even parity solutions have a $\mathrm{Nu}\left(\alpha \right)$ relationship that is even about $\alpha =0$.

• Odd parity solutions consist of two different sets where one set consists of flows with one more cell having a natural circulation than having an anti-natural circulation. The other set has one more cell with anti-natural circulations. The former persists to larger inclinations and transfers more heat.

• Even parity solutions can become of odd parity by either losing or gaining a cell at a bifurcation.

• Odd parity solutions can transform into other patterns with odd parity by losing or gaining an even number of cells. The same mechanism applies when even parity solutions transform gradually into other even parity solutions.

• At sufficiently high inclinations, the flow is unique and consists of just one circulating cell that is driven by the natural buoyancy forces on the hot and cold boundaries.

• When the aspect ratio is sufficiently large, then S-shaped solution trajectories are obtained, as are tongue-like shapes. These have been found to interact as the aspect ratio changes.

• For a fixed value of $\mathrm{Ra}$, the evidence so far suggests that there is a maximum inclination angle above which no nonlinear convection arises, apart from the single-cell circulation.

• For a fixed aspect ratio, the evidence so far shows that it is possible to have strong convection when $\alpha$ is larger than $60°$, and this suggests strongly that larger inclinations may be achieved.

Whilst we have presented a carefully curated and very large set of results involving many hundreds of hours of computation, there remains much more that could be done. One possible desire would be the creation of an animation which marches slowly through the equivalents of Figure 12 and Figure 13 as L increases in small increments. This would give a sense of understanding that a few printed figures do not. The same would be true for a fixed value of L, but where $\mathrm{Ra}$ increases slowly, and we would then be able to witness how the flows might evolve as the physically variable parameter, $\mathrm{Ra}$, varies. We would anticipate that hysteretical effects will arise. Other possible questions include those of distinguishing between solutions that are stable (i.e., computable with a time-dependent solver) and those that are not. Another one would query how precisely does the maximum inclination for strongly nonlinear convections vary with both $\mathrm{Ra}$ and L, and can the large-L limit be subject to some type of combined analytical/numerical treatment? Additionally, whilst the present numerical scheme has allowed us to follow the twists and turns of the solution trajectories, it has not been sufficiently able to follow every curve easily, particularly when the change in direction is very rapid, or when there is a neighbouring but distinct trajectory that will sometimes cause the Newton–Raphson scheme to slow down considerably or not to converge at all.

Author Contributions

All the authors were invaluable to the creation of this paper. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:

Latin letters
$A_{i,j}$	Fourier coefficient for $\psi$
$B_{i,j}$	Fourier coefficient for $\theta$
g	gravity
h	height of the channel
k	wave number
K	permeability
L	cavity aspect ratio
n	number of cells
$\mathrm{Nu}$	Nusselt number
${\mathrm{Nu}}_{\mathrm{GM}}$	Nusselt number in [8]
$\mathrm{NX}$	number of Fourier modes in the x-direction
$\mathrm{NZ}$	number of Fourier modes in the z-direction
p	pressure
$\mathrm{Ra}$	Darcy–Rayleigh number
t	time
T	dimensional temperature
u	Darcy velocity along the layer
v	spanwise Darcy velocity
w	Darcy velocity across the layer
x	coordinate along the layer
y	spanwise coordinate
z	coordinate across the layer
Greek letters
$\alpha$	inclination angle
$\beta$	coefficient of cubical expansion
$\theta$	temperature
$\kappa$	thermal diffusivity
$\mu$	dynamic viscosity
${\rho }_f$	reference density
$\sigma$	heat capacity ratio
$\psi$	streamfunction
Other symbols
c	cold
h	hot
$\overline{}$	dimensional quantity

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Depicting an inclined cavity with aspect ratio, L, and inclination, [Forumla omitted. See PDF.]. The boundary conditions are also indicated, as is the isotherm pattern for a weak anticlockwise unicellular flow.

Figure 2. Comparing the solution curves of (a) Guerrero-Martínez et al. [8] with (b) the present computations when [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.]. Note that [Forumla omitted. See PDF.], which represents the total heat flux from the lower surface, while [Forumla omitted. See PDF.] represents the mean heat flux per unit length. The number of cells within the cavity is denoted by the value of n. The black disks indicate the cases used for the streamline and isotherm patterns displayed in Figure 3, while the yellow disks represent the two [Forumla omitted. See PDF.] cases shown in Figure 4, and the green disks the three [Forumla omitted. See PDF.] cases shown in Figure 5.

Figure 3. Four different flow patterns for the case of [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.]. The continuous lines are streamlines, while the dashed lines represent isotherms. These cases correspond to the black disks shown in Figure 2b.

Figure 4. Streamlines and isotherms for two different [Forumla omitted. See PDF.] flows when [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.]: (a) displays the case when the outer cells exhibit the natural circulation, while (b) displays the case when they exhibit anti-natural circulation. These cases correspond to the yellow disks shown in Figure 2b in descending order.

Figure 5. Streamlines and isotherms for three different [Forumla omitted. See PDF.] flows when [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.]. All three cases have an anti-natural cell next to the right-hand boundary. These cases correspond to the green disks shown in Figure 2b in descending order.

Figure 6. The variation of [Forumla omitted. See PDF.] with inclination for the unit cavity, and for the indicated selection of values of [Forumla omitted. See PDF.]. The solutions here are restricted to [Forumla omitted. See PDF.] cells.

Figure 7. The variation of [Forumla omitted. See PDF.] with inclination for the unit cavity, and for the indicated selection of values of [Forumla omitted. See PDF.]. The [Forumla omitted. See PDF.] curves are shown in black, [Forumla omitted. See PDF.] in red, and [Forumla omitted. See PDF.] in green. The dashed curves correspond to even values of n.

Figure 8. The variation of [Forumla omitted. See PDF.] with inclination for the [Forumla omitted. See PDF.] cavity, and for the indicated selection of values of [Forumla omitted. See PDF.]. The solutions here are restricted to [Forumla omitted. See PDF.] cells.

Figure 9. The variation of [Forumla omitted. See PDF.] with inclination for the [Forumla omitted. See PDF.] cavity, and for the indicated selection of values of [Forumla omitted. See PDF.]. The [Forumla omitted. See PDF.] curves are shown in black, [Forumla omitted. See PDF.] in red, [Forumla omitted. See PDF.] in green, [Forumla omitted. See PDF.] in blue, and [Forumla omitted. See PDF.] in orange. The dashed curves correspond to even values of n.

Figure 10. The variation of [Forumla omitted. See PDF.] with inclination for the [Forumla omitted. See PDF.] cavity, and for the indicated selection of values of [Forumla omitted. See PDF.]. The curve corresponding to [Forumla omitted. See PDF.] is identified as being an [Forumla omitted. See PDF.] solution for small values of [Forumla omitted. See PDF.].

Figure 11. The variation of [Forumla omitted. See PDF.] with inclination for the [Forumla omitted. See PDF.] cavity, and for the indicated selection of values of [Forumla omitted. See PDF.]. The same colour scheme for the modes has been adopted here, but with the addition of [Forumla omitted. See PDF.] in purple, [Forumla omitted. See PDF.] in brown, and [Forumla omitted. See PDF.] in grey.

Figure 11. The variation of [Forumla omitted. See PDF.] with inclination for the [Forumla omitted. See PDF.] cavity, and for the indicated selection of values of [Forumla omitted. See PDF.]. The same colour scheme for the modes has been adopted here, but with the addition of [Forumla omitted. See PDF.] in purple, [Forumla omitted. See PDF.] in brown, and [Forumla omitted. See PDF.] in grey.

Figure 12. The variation of [Forumla omitted. See PDF.] with inclination for when [Forumla omitted. See PDF.] cavity, and for the indicated selection of values of cavity aspect ratios, L.

Figure 13. The variation of [Forumla omitted. See PDF.] with inclination for the case of [Forumla omitted. See PDF.], and for the indicated aspect ratios, L. Not all solution curves have been displayed, see the text.

Figure 14. Comparison of the trajectories of the solution curves for [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.]. Both subfigures represent one continuous solution curve.

Figure 15. Showing the maximum inclination for the presence of nonlinear convection as a function of (a) the aspect ratio, L and of (b) [Forumla omitted. See PDF.]. The green disks represent the solution branch containing the transition from an [Forumla omitted. See PDF.] to an [Forumla omitted. See PDF.] pattern. The orange disks represent the solution branch containing the transition from an [Forumla omitted. See PDF.] to an [Forumla omitted. See PDF.] pattern.

Figure 16. The right-most solution curves for [Forumla omitted. See PDF.] and for the given values of [Forumla omitted. See PDF.]. The dashed red line connects the points of maximum inclination.

Table 1

Cross-validation between (i) the finite difference method with Richardson’s extrapolation and (ii) the spectral scheme. Values of $B_{1,1}$ and $B_{0,2}$ are shown for the case where $\mathrm{Ra}=100$, $L=1$, and $\alpha =20°$.

(i)	$\mathtt{N}$	${\mathit{B}}_{1,1}$	${\mathit{B}}_{0,2}$		(ii)	$\mathrm{NX}\times \mathrm{NZ}$	${\mathit{B}}_{1,1}$	${\mathit{B}}_{0,2}$
	16	$-0.404918$	$-0.228235$	raw data		1×2	$-0.295059$	$-0.218866$
	32	$-0.403361$	$-0.226810$			3×4	$-0.382166$	$-0.222802$
	64	$-0.402951$	$-0.226431$			5×6	$-0.399483$	$-0.225658$
	128	$-0.402847$	$-0.226334$			7×8	$-0.402290$	$-0.226203$
	32	$-0.402842$	$-0.226335$	1st		9×10	$-0.402719$	$-0.226288$
	64	$-0.402814$	$-0.226305$	extrapolate		11×12	$-0.402791$	$-0.226301$
	128	$-0.402812$	$-0.226302$			13×14	$-0.402806$	$-0.226302$
	64	$-0.402812$	$-0.226303$	2nd		15×16	$-0.402810$	$-0.226302$
	128	$-0.402812$	$-0.226302$	extrapolate		17×18	$-0.402812$	$-0.226302$

Table 2

Showing the maximum inclinations for the $n=1$ to $n=3$ and the $n=1$ to $n=5$ fold bifurcations.

L	$\mathit{n}=1$ to $\mathit{n}=5$	$\mathit{n}=1$ to $\mathit{n}=3$
3	$14.59°$	$32.06°$
4	$26.42°$	$41.07°$
5	$35.03°$	$45.75°$
6	$40.67°$	$48.20°$
7	$44.29°$	$49.46°$
8	$46.61°$	$50.09°$

References

1. Horton, C.W.; Rogers, F.T. Convection currents in a porous medium. J. Appl. Phys.; 1945; 16, pp. 367-370. [DOI: https://dx.doi.org/10.1063/1.1707601]

2. Lapwood, E.R. Convection of a fluid in a porous medium. Proc. Camb. Phil. Soc.; 1948; 44, pp. 508-521. [DOI: https://dx.doi.org/10.1017/S030500410002452X]

3. Nield, D.A.; Bejan, A. Convection in Porous Media; 5th ed. Springer: New York, NY, USA, 2017.

4. Rees, D.A.S.; Bassom, A.P. Onset of Darcy-Bénard convection in an inclined porous layer heated from below. Acta Mech.; 2000; 144, pp. 103-118. [DOI: https://dx.doi.org/10.1007/BF01181831]

5. Wen, B.; Chini, G.P. On moderate-Rayleigh-number convection in an inclined porous layer. Fluids; 2019; 4, 101. [DOI: https://dx.doi.org/10.3390/fluids4020101]

6. Wen, B.; Chini, G.P. Inclined porous medium convection at large Rayleigh number. J. Fluid Mech.; 2018; 837, pp. 670-702. [DOI: https://dx.doi.org/10.1017/jfm.2017.863]

7. Rees, D.A.S.; Barletta, A. When does the onset of convection in an inclined porous layer become subcritical?. Mech. Res. Comm.; 2022; 125, 103992. [DOI: https://dx.doi.org/10.1016/j.mechrescom.2022.103992]

8. Guerrero-Martínez, F.J.; Karimi, N.; Ramos, E. Numerical modelling of multiple steady-state convective modes in a tilted porous medium heated from below. Int. Com. Heat Mass Transf.; 2018; 92, pp. 64-72. [DOI: https://dx.doi.org/10.1016/j.icheatmasstransfer.2018.02.009]

9. Riley, D.S.; Winters, K.H. A numerical bifurcation study of convection in a tilted two-dimensional porous cavity. J. Fluid Mech.; 1990; 215, pp. 309-329. [DOI: https://dx.doi.org/10.1017/S002211209000266X]

10. Rees, D.A.S. The stability of Darcy-Bénard convection. Handbook of Porous Media; Vafai, K. Marcel Dekker: New York, NY, USA, 2000; pp. 521-558.