Three-Dimensional Convective Planforms for Inclined Darcy-Bénard Convection

: We investigate the onset of convection in an inclined Darcy-Bénard layer. When such a layer is unbounded in the spanwise direction it is generally known that longitudinal rolls comprise the most unstable planform. On the other hand, when a layer has a sufﬁciently small spanwise width, then transverse rolls form the most unstable planform. However, the layer remains stable to transverse roll disturbances when the inclination is above roughly 31 degrees from the horizontal. This paper considers the transition between these two extreme cases where the spanwise width takes moderate values and where rectangular cells are considered. It is found that the most unstable planform is quite strongly sensitive to the magnitude of the spanwise width and that there are large regions of parameter space within which three-dimensional convection patterns have the smallest critical Darcy-Rayleigh number.


Introduction
The Rayleigh-Bénard instability in a porous layer saturated by a fluid and heated from below has been the subject of many investigations over the last decades. A thorough and updated survey of the most important results in this field can be found in the book by Nield and Bejan [1]. The alternative name used for this instability when a porous medium is involved is Darcy-Bénard as the momentum transfer for seepage flows is modelled through Darcy's law. Usually, the porous layer where the thermoconvective instability arises in the form of cellular patterns is assumed to be horizontal and bounded by impermeable walls kept at different uniform temperatures. This situation has a stationary basic state where no flow occurs and a vertical temperature gradient directed downward causes the instability. There is a different scenario when the porous layer is inclined at an angle α to the horizontal. In this case, the basic seepage velocity is a parallel vector field with a non-uniform linear profile yielding a longitudinal flow with zero net rate.
Pioneering studies of the Darcy-Bénard problem in an inclined porous layer were carried out by Bories and Combarnous [2] and by Weber [3]. The main consequence of the inclination is that the critical Rayleigh number at the onset of the instability depends on the inclination α through a factor 1/ cos α. The ultimate effect of this factor is that no convective instability can arise when the porous layer becomes vertical (α = 90 • ). A rigorous analysis of this effect is reported in a short paper by Gill [4]. The consequences of Gill's no-go theorem for the instability in a vertical layer were studied further by Lewis, Bassom and Rees [5] using asymptotic theory and by Straughan [6] using nonlinear stability theory. It must be mentioned that the scaling of the critical Rayleigh number through the factor 1/ cos α holds whenever the most unstable modes at onset are longitudinal rolls. This is the typical situation when the inclined porous layer has no lateral confinement. On the other hand, the effect of

Governing Equations and Stability Analysis
The aim is to investigate the effect of the presence of sidewalls on the onset of convective instability for the inclined Darcy-Bénard problem. Specifically the intention is to determine the identity of the most unstable mode without being restricted solely to transverse or longitudinal rolls. The nondimensional governing equations are where z is a horizontal coordinate and where α is the angle of rotation of the layer about the z-axis.
In addition, y is the coordinate which acts in the normal direction to the heated surfaces, and x is the coordinate up the layer. The above equations form the standard set for problems of this type, and the Darcy-Rayleigh number, Ra, is given by where d is the thickness of the uniform layer where ρ is the reference density, g the strength of gravity, β the coefficient of cubical expansion, ∆T the temperture drop across the layer, µ the dynamic viscosity, and κ the thermal diffusivity of the porous medium. The boundary conditions which are required to complete the specification of the problem are, We note that the conditions on the spanwise sidewalls at z = 0, L also imply that ∂p/∂z = 0 there. From these equations we see that the basic state, whose stability forms the subject of the paper, is given by, These profiles may be subtracted out of the governing equations and the resulting system is then linearised. Equations (1) to (4) then remain unchanged, but Equation (5) takes the form, The velocity components may now be removed and the resulting pair of equations in pressure/temperature form are, The boundary conditions are now that, We shall consider three-dimensional disturbances which are proportional to cos(nπz/L), i.e., where n rolls have been fitted into the spanwise width of the layer. This form may be used for both the pressure and the temperature because each satisfies a zero Neumann condition at the sidewalls. However, there also needs to be a component in the x-direction in order to allow for the presence of three-dimensional convection, and therefore we may take the disturbance temperature field to be proportional to cos(nπz/L) cos(mx) where m is a wavenumber. This product of cosines, which represents convection with a rectangular planform, may also be split into the sum of two independent oblique rolls, cos(nπz/L) cos(mx) = 1 4 e i(nπz/L−mx) + e i(nπz/L+mx) + complex conjugate, which are equally inclined away from the direction of the x-axis but in opposite directions. We note that it is sufficient to consider just one of these rolls since the other has identical stability properties.
Therefore in the rest of the paper we shall refer to an oblique roll disturbance as a proxy for the 3D cell with the rectangular planform given that it may be visualised easily. The presentation of the stability analysis is made easier by defining an overall wavenumber, k, where k is given by, and where φ is its orientation. We note that φ = 0 corresponds to longitudinal rolls, and φ = π/2 to transverse rolls (although such rolls will also require n = 0). Although our computations will use k, it is the spanwise width, L, and the nuber of cells, n, which are imposed on the configuration. Therefore we have, this dependence of k on φ will be relied upon later. Therefore we shall take the following form for the disturbances, where λ = λ r + iλ i is the complex exponential growth rate. The functions f (y) and g(y) satisfy, This will be treated as an eigenvalue problem for both Ra and λ i where λ r = 0 is imposed and which corresponds to marginal stability. The boundary conditions are taken to be, where the last condition is a normalisation condition. This complex eigenvalue problem was solved using a fourth order Runge Kutta method as part of a shooting method code. Computations made with this code, which uses a pressure/temperature formulation, compare perfectly with those of Rees and Bassom [7] who used a streamfunction/temperature formulation. Table 1 shows how the computed values of the critical Darcy-Rayleigh number vary with successive interval-halving where N is the number of uniform intervals. Two different cases are presented and they demonstrate that the numerical errors reduce by a factor of roughly 2 4 as the intervals are halved, which also confirms the 4th order accuracy of the method. We have therefore decided to use a uniform grid of 100 intervals in all of our computations; this may be seen to provide eight significant figures of accuracy. It is important to note that while some of the computations presented in Rees and Bassom [7] correspond to travelling modes, i.e., they correspond to nonzero values of λ i , these cases arise in regions of parameter space which are quite distant from where the minima in the neutral curves are located. Therefore the present computations, which are concerned with such minima, always give rise to stationary modes for which λ i = 0.

Results and Discussion
The ultimate aim is to provide an informative modal map in the spirit of Beck's well-known analysis of Darcy-Bénard convection in a cuboid [21]. This will show the preferred planform for convection at onset as a function of both the spanwise width of the layer and the inclination angle. We shall begin with an attempt to understand how the neutral curves vary with the governing parameters.

Sample Neutral Curves
The four subfigures which comprise Figure 1 display how the critical Darcy-Rayleigh number varies with the roll orientation, φ for the inclination angle, α = 5 • . The subfigures are distinguished by having different numbers of rolls occupying the layer. In all cases the equivalent critical Darcy-Rayleigh number for the transverse mode is given by the blue line. A rather large number of curves are displayed here and it is impractical to attempt to label them. Therefore the manner in which the curves have been plotted serves to identify the values of L. Thus red curves correspond to values of L which are less than n, and dashed curves denote cases where L is an integer. The right-most curve corresponds to L = 5. Therefore, in all four cases shown, the curve with the lowest minimum value may be seen to correspond to when L = n, given that it is dashed and is the first instance of a black curve as one moves to the right. In fact, each of the L = n curves shown have identical shapes because of the form of Equation (15) which means that the linear stability equations are solved using precisely the same value of k. This observation may also be generalised to saying that any two cases with identical values of n/L, φ and α are identical mathematically and have the same critical Darcy-Rayleigh number.
We note that the curves corresponding to different widths, L, always have their minimum below the critical value of Ra for transverse rolls, and therefore we do not expect transverse rolls to form the most unstable mode unless the layer is sufficiently narrow, i.e., unless L is sufficiently small. In addition, the roll orientation is not always that of the longitudinal roll (φ = 0), but φ often takes nonzero values and corresponds to an oblique mode.
We mention that the locus of the minima as L increases continuously traces the same curve independently of the value of n, and this will also be true for a similar set of curves for different inclinations. Therefore we shall not present such a large amount of data for other inclination. Figure 2 shows how such neutral curves vary with the inclination angle, and we have chosen the 2-roll case, n = 2, for these. The main general observation is that the Darcy-Rayleigh number increases with increasing inclination, should n, L and φ be held fixed, and this also includes the transverse roll. Yet again, for any chosen value of L and α the minimum value of the Darcy-Rayleigh number lies below that of the transverse roll, and it may or may not correspond to φ = 0.
Let us consider, for example, the α = 10 • case, and specifically at the two black curves which intersect the vertical axis one moves upwards from Ra = 40. The dashed line, which corresponds to L = 2, has its minimum at φ = 0 and therefore it is a longitudinal roll. The second curve, a continuous curve, has its minimum at φ 22 • and is an oblique roll. Therefore this is a case where neither the transverse nor the longintudinal roll forms the preferred mode, a situation which, as will be seen, is quite common for low inclinations.
Curves for L > 5 follow the same trend, namely that their minima correspond to increasing values of both Ra and φ. In fact Ra will never exceed the value for the tranverse roll while φ tends towards 90 • . 2) · · · 5 when the inclination from the horizontal is α = 5 • and for the n = 1, n = 2, n = 3 and n = 4 modes. The red curves denote those values of L which are less than n. The dotted lines correspond to integer values of L while the rightmost curve is for L = 5. The blue horizontal line corresponds to transverse rolls.
In Figure 3 neutral curves are shown again but from a different point of view. We have taken the inclination α = 20 • to represent a typical case, and have concentrated on aspect ratios varying in steps of 0.2 from L = 3 to L = 4. The different curves correspond to different values of n and these are colour-coded in the manner defined in the caption.
When L = 3, the preferred mode is the n = 3 mode which is displayed in red. When L increases to 3.2, the various minima in the curves tend to move either to the right (when φ is nonzero) or downwards (when φ is zero). This process continues in the same way except that when L = 3.4, the red n = 3 curve now has a local maximum at φ = 0 and the global minimum, which is difficult to see in the figure, is at φ 25.8 • . Therefore the preferred mode is now an oblique roll. At L = 3.6 the critical Darcy-Rayleigh number for the n = 3 longitudinal is now larger than that for the n = 4 roll, shown in purple. However, the oblique n = 3 mode remains the preferred pattern. When L = 3.8 and L = 4, the n = 4 longitudinal roll assumes the role of being the most dangerous.
Although we have chosen to consider this transition using the case α = 20 • and when L varies between 3 and 4, this scenario is ubiquitous. There are two exceptions, however: (i) when L is sufficiently far below 1 that transverse rolls are favoured, and (ii) when the inclination is sufficently large. So Figures 4-6 begin the the process of understanding the global picture of the linear stability properties of this system. These figures, which correspond respectively to α = 10 • , α = 20 • and α = 25 • , show how the critical Darcy-Rayleigh number and the roll orientation varies with aspect ratio, L. In the Darcy-Rayleigh number plots we show curves for transverse, longitudinal and oblique rolls. The blue curves denote longitudinal rolls, and in Figure 4, for the longitudinal mode which has its minimum at L = 1, it needs only a very slight increase in L from L = 1 to L = 1.023303 for the longitudinal mode's variation in φ to make the above-mentioned transition from φ = 0 being a minimum to being a local maximum, and then oblique rolls form the preferred planform. The transitional value, L = 1.023303, may be referred to as a quartic point.
The detailed variation of φ with L is shown in the lower part of the figure and it resembles the separate buildings of the Sydney Opera House in its disjointed shape. The sudden transitions in φ back to zero correspond to the blue disks in the upper plot where, for example, the n = 1 red curve crosses the n = 2 blue curve and the preferred mode changes its identity. Figures 5 and 6 have similar features except that it becomes clear that there are only a finite number of 'buildings' in the Sydney Opera House shape. The reason is that the quartic point for each inclination, the red disks in Figures 4-6, all occur at the same value of the Darcy-Rayleigh number, while the intersections between the curves corresponding to neighbouring longitudinal mode curves occur at decreasing values of the Darcy-Rayleigh number. So for α = 20 • there are four disconnected ranges of values of L for which oblique modes arise, while for α = 25 • it is only two, and for α = 30 • , which is not shown here, it is just the one. It is possible to determine quite easily what happens when α = 0. The neutral curve for an unbounded horizontal plane layer is given by, where, for the present bounded problem, k is given by Equation (15). Hence we have, Ra = π 2 (L 2 cos 2 φ + n 2 ) 2 n 2 L 2 cos 2 φ .
For a given number of rolls, n, this expression is maximised when L cos φ = n or cos φ = n/L, at which point Ra = 4π 2 . When L is an integer, the case n = L yields an n-cell longitudinal roll pattern (φ = 0) as the most unstable disturbance; other values of n give larger values of the Darcy-Rayleigh number irrespective of the roll orientation. When L is not an integer, then the most unstable disturbance takes the form of oblique rolls with n cells where n is any integer which is smaller than L. For example, when L = 3.5, then there are three possible planforms where Ra c = 4π 2 , namely n = 1 (cos φ = 1/3.5), n = 2 (cos φ = 2/3.5) and n = 3 (cos φ = 3/3.5). These curves are shown as the black dashed lines in Figure 7 when we see that there will always be a degeneracy whenever L > 2 by having more than one choice of pattern with the smallest possible critical Darcy-Rayleigh number. Although we shall not prove it here once α rises about zero be it ever so slightly, this degeneracy in the identity of the preferred mode resolves. The unique preferred mode then becomes the one with the largest value of n which is less than L, and this is shown as the red line. The variation in the orientations of the oblique mode with L for this very slightly-inclined case is shown in Figure 7, and these too resemble the Sydney Opera House shape. Moreover, if one compares the dependence of Φ on L as shown in Figures 4-6 in that order, then the dependence shown in Figure 7 is clearly seen to be reasonable. One implication of Figure 7 is that longitudinal rolls can only form the preferred onset pattern in a horizontal layer when L is an integer.

The Modal Map
Finally, we are in a position to provide the modal map which summarises the stability properties of the inclined layer with a finite spanwise width, and this is shown in Figure 8. There are three main regions and these are colour-coded for easy reference. The blue region corresponds to when we expect transverse rolls to appear first as Ra increases. The region is bounded to the right at L = 1 and above by α = 31.49032 • , which is the maximum inclination for linear instability (Rees and Bassom 2000).
The white regions are when longitudinal rolls are preferred. For a randomly chosen value of L, the appearance of longituinal rolls becomes decreasingly likely as the layer tends towards being horizontal, and the ever-decreasing ranges of values of L tend to be focused near integer values of L.
The orange shark fin shapes correspond to when oblique rolls (i.e., rectangular cells) are favoured. The left hand boundaries correspond to when the neutral curves (in terms of Ra against φ) have a quartic point, and therefore Ra satisfies d 2 Ra/dφ 2 = 0; such locations are shown as the blue disks in Figures 4-6. Numerically, this curve was obtained by solving Equations (A6)-(A11) given in the Appendix A (although the Appendix A also provides the analytical solution).
The right hand boundaries correspond to the red disks in those Figures, i.e., where an oblique mode for n rolls has the same critical Darcy-Rayleigh number as a longitudinal mode with n + 1 rolls. This curve was obtained by solving Equations (5) and (7) together with (A8) and (A9) which guarantees that the minimum in the neutral curve has been obtained. The value of α as a function of L was obtained by adding the further constraint that the oblique mode with n cells has the same critical Darcy-Rayleigh number as the longitudinal mode with n + 1 cells.
When these two curves meet at the sharp end, both sides of the fin tip represent longitudinal rolls, the left with say n rolls and the right with n + 1 rolls. This is well-known to take place when L = n(n + 1), for this is what happens for 2D convection in a porous box of aspect ratio, L (cf. p313 of Nield and Bejan 2017).
The Appendix contains a detailed analysis for the left hand sides of the fins, and this provides the following formula relating α to n and L: which compares favourably to 8 significant figures with the numerical calculations shown in Figure 7. When L = n(n + 1) is substituted into Equation (23) we obtain tan 2 α = 96n(n + 1) 2 (2n + 1) 3 which gives the value of α at the tip of the fin. For n = 1 and n = 5, for example, Equation (24) yields, tan 2 α = 128 9(21 + π 2 ) and tan 2 α = 17820 1331(93 + π 2 ) . (25) Some numerical data which has been obtained using Equation (24) are given in Table 2, below.   Table 2 shows that, as L increases, the range of inclinations within which we might observe oblique rolls decreases as the aspect ratio increases. Indeed, it is quite easy to use Equation (24) to show that We have also presented lines of constant φ within the orange regions of Figure 8 to indicate concisely which oblique modes are to be expected. The form of these contours suggests that the maximum possible value of φ within any chosen shark fin shape is to be found at the bottom right hand corner. This value of φ may be found using Equation (22) using L = n + 1. Therefore we have, cos φ = n/(n + 1).
For the left-most shark fin we have φ = cos −1 (1/2) = 60 • while for the second we have φ = cos −1 (2/3) = 48.1897 • , and so on. As an extreme example, for n = 100, which is the final entry in Table 2, φ = cos −1 (100/101) = 8.0693 • . Therefore, when L is large, not only is there a small range of α within which which oblique rolls are preferred, but the deviation from the direction of the longitudinal roll is also small.

Conclusions
In this paper we have studied the effect of having a spanwise constraint on the onset of Darcy-Bénard convection in an inclined porous layer. When the layer is unconstrained in this way, it is well-known that longitudinal rolls always form the most unstable convection pattern and that the critical Darcy-Rayleigh number is 4π 2 / cos α. It is also well-known that if the layer is sufficiently well-constrained by having an exceptionally narrow layer then transverse rolls arise. Here we have relaxed the assumption that convection takes either of these two forms by considering cells with a rectangular planform, for which the oblique rolls are a proxy. It has been found that, for moderate values of the aspect ratio, L, and for inclinations which are less than, say, 30 • , quite a large amount of the (L, α) parameter space corresponds to oblique rolls forming the preferred pattern of convection.
We have been able to provide an analytical expression for when a longitudinal roll loses its dominance to oblique modes as L increases. We have also found a precise expression for when oblique rolls give way to longitudinal rolls as the inclination increases (viz. the apex of the shark fin shapes). Although transverse rolls are well known to be unstable only when α < 31.49032 • , we have found that it is possible to obtain oblique rolls at slightly higher inclinations than that, although α = 34.167212 is the upper limit.
If one were to wish to extend this analysis still further by considering an inclined cuboid, then the task will be very substantially more difficult. In the present case the basic state may be written down analytically, while this cuboidal extension will need its basic state computed numerically. But even then, there are complications because the concept of the onset of convection doesn't always apply (Guerrero-Martínez, Karimi and Ramos [22]). Even when confined to two-dimensional convection, some aspect ratios (typically those where L is close to being an even integer) will exhibit a classical bifurcation from the basic state to a strongly convecting solution, but others (where L is close to being an odd integer) exhibit a smooth transition to strong convection. Another complicating factor is that the work of Wen and Chini [19] has uncovered subcritical onset of convection. Solutions are sought in the form of a power series in φ. Therefore we begin with the expansions,

Abbreviations
and then it is necessary to solve the following systems, The terms typeset in blue in Equations (A10) and (A11) require some explanation. Equation (15) may be expanded according to, This means that the k 2 terms in Equations (A2) and (A3) are, formally, where k 0 = nπ/L. In Equations (A6) to (A11) the 0-subscripts have been omitted for clarity of presentation.
In the above systems we have also set Ra 1 = 0 for reasons of symmetry, and therefore this term has been omitted. Although we have retained Ra 2 in the above, it is our intention to set Ra 2 = 0 in order to determine the relationship between L and α which defines the left hand side of each of the sharks fins in Figure 8.
At leading order, we shall use the following as the solutions, At O(φ) we have the equations, which defines the quantity, A = ikRa 0 sin α, for the sake of brevity. First, we note the presence of the following Complementary Function, where and where B is presently arbitrary. The Particular Integral may be found making the following substitutions into Equations (A15) and (A16), When considering those terms which are proportional to (y − 1 2 ) 2 , Equations (A15) and (A16) yield, These equations are multiples of one another and they both yield, When considering those terms which are proportional to (y − 1 2 ), we obtain, If we multiply Equation (A25) by π, Equation (A26) by (k 2 + π 2 ) and add, then the left hand sides cancel and we are left with the solvability condition, This, together with Equation (A24) yields the following expressions for C and F: When these expressions are substituted into Equations (A25) and (A26), both give the following relation betwen B and E: Now we consider those terms which are not proportional to a power of (y − 1 2 ). These give, The left hand sides of these equations may be eliminated by multiplying Equation (A30) by π and subtracting Equation (A31) multiplied by (k 2 + π 2 ). This eventually simplifies to, Equations (A32) and (A29) forms a pair of equations for B and E. These values are found to be When these expressions are substituted into Equations (A30) and (A31), both give the following relation between A and D: Thus far, four of the constants which were introduced in Equations (A20) and (A21) have been found together with a relationship between the other two, namely Equation (A34). When this Particular Integral is added to the Complementary Function given in Equations (A17) and (A18) there remains two constants to determine. This is undertaken by ensuring that the appropriate boundary conditions are satisfied by f 1 and g 1 : Although there are four boundary conditions the symmetry of the solutions shown in Equations (A20) and (A21) means that they represent only two in practice and therefore the full set of constants may now be determined.
The condition, f 1 (1) = 0, yields, − πD − E − 1 4 πF + (2k 2 + π 2 )(k 2 + π 2 ) k 2 B sinh 1 2 σ = 0, while the condition, g 1 (1) = 0, yields, We may combine these two equations together so that Equation (A34) may be used to eliminate A and D. Thus we subtract Equation (A36) multiplied by k 2 from Equation (A37) in order to obtain, (k 2 + π 2 ) π 2 − (2k 2 + π 2 ) B sinh 1 2 σ = π 2 (k 2 + π 2 )A − πk 2 D + 1 4 π 2 (k 2 + π 2 )C − k 2 (E + 1 2 πF). (A38) Using the above solutions for C, E and F together with Equation (A37) we obtain, and therefore this complementary function plays no role in the final solution. Now we use Equations (A37) and (A34) separately to show that, At O(φ 2 ) it is necessary to guarantee that Equations (A10) and (A11) have a solution. At this stage in such an analysis it would be usual to be computing the value of Ra 2 , but here we are attempting to find an analytical expression which relates L and α for the specific case where φ = 0 corresponds to where ∂ 2 Ra/∂φ 2 = 0 -such a locus delineates the region where φ = 0 is a local minimum in the neutral curve from the region where nonzero values of φ yield the minimum.