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

Onset of Convection in an Inclined Anisotropic Porous Layer with Internal Heat Generation Leiv Storesletten 1, D. Andrew S. Rees 2,∗ 1 Department of Mathematics, University of Agder, Serviceboks 422, 4604 Kristiansand, Norway; leiv.storesletten@uia.no 2 Department of Mechanical Engineering, University of Bath, Bath BA2 7AY, U.K.; D.A.S.Rees@bath.ac.uk * Corresponding author Version April 10, 2019 submitted to Fluids Abstract: The onset of convection in an inclined porous layer which is heated internally by 1 a uniform distribution of heat sources is considered. We investigate the combined effects of 2 inclination, anisotropy and internal heat generation on the linear instability of the basic parallel flow. 3 When the Rayleigh number is sufficiently large, instability occurs and a convective motion is set up. 4 It turns out that the preferred motion at convection onset depends quite strongly on the anisotropy 5 ratio, ξ, and the inclination angle. When ξ < 1 the preferred motion is in the form of longitudinal 6 rolls for all inclinations. When ξ > 1 transverse rolls are preferred for small inclinations but, at 7 high inclinations, longitudinal rolls are preferred. At intermediate inclinations the preferred roll 8 orientation varies smoothly between these two extremes. 9


Introduction
The onset of convection in porous layers forms one of the classical problems of stability theory and the most well-known of these is the the Darcy-Bénard (or Horton-Rogers-Lapwood) problem [1,2]. In its simplest form a horizontal porous layer for which Darcy's law applies has a lower boundary which is impermeable and uniformly hot while the impermeable upper boundary is uniformly cold. It is well-known that the critical Darcy-Rayleigh number for the onset of convection is 4π 2 when the layer is unbounded horizontally. The corresponding wavenumber is π. A second well-known problem of the same type comprises a porous layer with uniformly cold boundaries but which is heated via a uniform internal heat generation mechanism. Practical examples of such systems include heat removal from fuel debris in nuclear reactors, the underground disposal of radioactive waste materials, agricultural applications involving the storage of foodstuff, and exothermic chemical reactions in packed-bed reactors. It is this configuration which we study here and we shall determine the combined effects of layer inclination and anisotropy on the onset of convection.
The first experimental and theoretical studies on the topic of an internally-heated layer were undertaken in the late 1970s by Kulacki and Ramchandani [3], Gasser and Kazimi [4], Buretta and Berman [5], Hardee and Nilson [6], Hwang and Marr [7], Tveitereid [8], Rhee et al. [9] and Kulacki and Freeman [10]. These papers represent quite a wide range of different types of boundary condition such as a free upper surface in [4] or the presence of an insulated bounding surface [8,10]. The case where both surfaces are impermeable and are at identical constant temperatures has a critical Darcy-Rayleigh number of 471.38 and a critical wavenumber of 4.6752, values which have been obtained by various authors such as Barletta et al. [11], Nandal and Mahajan [12] and Nouri-Borujerdi et al. [13,14]. This ten-fold rise in the critical Darcy-Rayleigh number may be explained at least qualitatively by the facts that (i) the internally-heated case admits a destabilising temperature gradient only in the upper

Governing Equations and the Basic State
We consider the onset of convection in an inclined porous layer of infinite lateral extent which is heated internally by a uniform distribution of heat sources. The layer is bounded by two parallel solid plates which are separated by a distance h and inclined at an angle α to the horizontal. Both boundaries are held at a constant temperature T 0 . The horizontal z-axis forms the direction about which the layer has been rotated, the y-axis is perpendicular to the bounding plates, and the x-axis lies in the lower surface pointing up the plane. The detailed configuration is shown in Figure 1. Moreover, the anisotropic porous layer is taken to be transversely isotropic in the permeability but homogenously isotropic in the thermal diffusivity. The saturating fluid is assumed to be such that Darcy's law holds and the Oberbeck-Boussinesq approximation is valid. The governing equations are then, where q is the rate of heat generated within the porous medium per unit volume, β is the coefficient of cubical expansion, k m is the thermal conductivity of the porous medium, and p is the pressure, which includes the hydro-static pressure. Terms in Equations (1)-(5) have their usual meaning in the porous medium context, but we mentioned that K L is the permeability in the x-direction, K T is the permeability in both the y-and z-directions, that σ is the heat capacity ratio, and κ is the thermal diffusivity. A zero subscript denotes reference quantities. Further details are in the Abbreviations at the end of the paper. The variables are made dimensionless by setting so that the governing equations reduce to, The dimensionless parameters which appear in these equations are the Darcy-Rayleigh number, Ra, and the anisotropy parameter, ξ, which are given respectively by, In this paper we shall use two different formulations of the governing equations for the numerical work. In the first instance two-dimensional cases will be considered using the streamfunction/temperature approach whereas for three-dimensional flows a pressure/temperature approach will be used. In the former case the streamfunction is defined using and therefore Equation (8) is satisfied while Equations (9) and (10) transform into, In the latter case (u, v, w) may be eliminated from Equations (8)-(10) and we obtain the following governing equations for the unsteady flow: For these two formulations we may write the boundary conditions in the form, ψ = ∂p ∂y = θ = 0 on both y = 0 and y = 1.
There exists a unique steady basic flow, v = (u b , 0, 0), ψ = ψ b (y), θ = θ b (y) and p = p b (y) where the net mass flux through a cross section x = constant is zero. This steady flow with nonlinear velocity and temperature profiles is given by, where p 0 is the pressure at (x, y, z) = (0, 0, 0).

Reduction to ODE Eigenvalue Form
The linear stability characteristics of the basic steady state given by Equations (18) and (19) may be investigated by setting, where ψ 1 and θ 1 are asymptotically small perturbations in the streamfunction and temperature, respectively. The substitution of the expressions (20) into the governing Equations (13) and (14) followed by linearisation yields the following set of perturbation equations, which are to be solved subject to homogeneous Dirichlet conditions. The disturbances may be reduced to ordinary differential form by factoring out a monochromatic horizontal component: we let, Hence the disturbance equations reduce to, and the boundary conditions are that, Ψ = Θ = 0 on both y = 0 and y = 1.
This system may be regarded as an eigenvalue problem for the complex growth rate, λ, as a function of k, ξ, Ra and α. Alternatively, if we let λ = λ r + iλ i and set λ r = 0 to model marginal stability, then this system may instead be regarded as an eigenvalue problem for both the Darcy-Rayleigh number, Ra, and for λ i as functions of k, ξ and α.
We may now define the phase velocity of the disturbances by modifying the exponential in Equation (23). At marginal stability λ = iλ i and therefore,

Numerical Results
In the first instance Equations (24)-(26) were discretised using second order accurate central difference approximations and the resulting difference equations were rearranged into a complex matrix eigenvalue problem for λ precisely in the manner described in Rees and Bassom [32]. For a chosen pair of values of ξ and inclination, α, the first few eigenvalues for λ were computed on a fine grid of values of Ra and k, and contours of λ r = 0 were plotted to provide the neutral curves. The use of N = 10 intervals in the range 0 ≤ y ≤ 1 gave poor results particularly when Ra or k was too large. The number of intervals was then increased first to N = 20, then to N = 40 and ultimately to N = 80, we found that, that the neutral curves presented in Figure 2 and which correspond to N = 40 exhibit no visual changes when comparing them with N = 80 solutions.  Figure 2 display the manner in which the neutral curves evolve as α increases in 1 • increments from the horizontal case, α = 0 • . For ease of comparison each subfigure corresponds to k and Ra lying in the ranges, 0 ≤ k 3k c and 0 ≤ Ra 12Ra c where k c and Ra c are the critical wavenumber and Darcy-Rayleigh number when the layer is horizontal (α = 0). We have plotted the neutral curves corresponding to the first two modes.
It is clear from Figure 2 that the general form taken by the neutral curves remains almost unchanged as ξ varies. Thus the critical values of Ra (which correspond to the minimum of each neutral curve and which are labelled as Ra c hereinafter) always increase as the inclination, α, increases and this may be expected given that buoyancy decreases with increasing inclination. For three of the cases shown in Figure 2 there is a local maximum inclination angle which is at an isola point (shown as black disks). These points are characterised by being surrounded by elliptic-shaped neutral curves where the ξ = 2 case provides a good example since the value of α at the isola itself lies just above 27 • . In these three cases (ξ = 5, 2 and 0.5) there also exists a saddle point; for ξ = 0.5 this is depicted as a red disk but the saddle points for ξ = 5 and ξ = 2 lie above the maximum value of Ra used for plotting. Figure 3 gives an alternative view which both summarises and extends the above discussion and shows, for a wide selection of values of ξ, how the local minimum and maximum of Ra on each neutral curve varies with inclination. The lower set of black disks and the associated dotted curve correspond to how the value of Ra for the isolae varies with inclination, and the upper set correspond to the saddle points. These two curves eventually merge when ξ is roughly 0.27843, which is in accord with the fact that no isola is displayed in Figure 2 for ξ = 0.2. We presume that each of the solid curves shown in Figure 3 tend asymptotically toward a constant value of α as Ra c → ∞, as was shown in Rees and Bassom [31] for the inclined Darcy-Bénard problem.
Some indication of how the critical wavenumber varies with α may be gleaned from Figure 2, but a more informative depiction is shown in Figure 4. Thus for relatively large values of ξ the critical wavenumber decreases at first as α increases from zero whereas k c increases when ξ takes smaller values. The presence of isolae is marked with black disks, and now it is clear that they correspond to a local maximum value of α in those curves. On the other hand, the saddle points, which are also marked with black disks, correspond to local minima. As ξ decreases the curves shown in Figure 4 gradually unkink themselves as the isola and the saddle point approach one another. Eventually, k c becomes a single-valued function of α.  is sufficiently large. The reason for this behaviour is related to the overall basic flow pattern and 175 specifically to the location of the disturbance within the layer. Using Eq. (18) we see that and therefore the basic state corresponds to descending flow near the boundaries and ascending flow 177 in the centre.
178 Figure 6 clarifies the reason why the phase velocity varies as it does in Figure 5. This figure   179 shows instantaneous snapshots of the disturbance isotherms for six different inclinations when ξ =    The linear stability characteristics of the basic steady state to generally oblique roll disturbances 195 may be investigated by setting,

Linear stability analysis for three-dimensional disturbances
where p 1 and θ 1 are asymptotically small perturbations in the pressure and temperature, respectively.

197
The substitution of the expressions (29) into the governing equations (15) and (16) gives the following 198 set of linearised perturbation equations where φ is the orientation of the roll relative to the x-axis, and where k and λ are once more the  In Figure 5 we show how the corresponding phase velocity, c, of the disturbances vary with α for the same set of values of ξ. Postive values mean that the disturbance ascends the layer when α is positive while negative values mean that the disturbance descends. Figure 5 shows that disturbances tend to ascend the layer when the inclination is small, but eventually descend when the inclination is sufficiently large. The reason for this behaviour is related to the overall basic flow pattern and specifically to the location of the disturbance within the layer. Using Equation (18) we see that and therefore the basic state corresponds to descending flow near the boundaries and ascending flow in the centre. Figure 6 clarifies the reason why the phase velocity varies as it does in Figure 5. This figure shows instantaneous snapshots of the disturbance isotherms for six different inclinations when ξ = 5. The first is the horizontal case where c = 0, while the second, third and fourth correspond to positive values of c where the fourth also corresponds to the maximum phase velocity. For larger inclinations c then decreases. The fifth frame corresponds to c = 0 while the sixth is at the isola point. Disturbances tend to be in the upper half of the layer because this is where the basic temperature gradient is destabilising. When inclinations are sufficiently small, the disturbances lie mostly in the region where the flow is up the layer. Moreover, the expression for u b ( 1 2 ) in Equation (28) shows that this value increases as α increases, and therefore the phase velocity of the disturbances also tends to increase with α at first. But as α increases further, disturbances become increasingly concentrated near the upper boundary where the basic flow is negative.
When other values of ξ are considered then figures like those shown in Figure 6 are obtained. The main change in the appearance of the analogous isotherm contours is the period of the cellular pattern because the critical wavenumbers vary quite strongly with ξ, as may be seen in Figure 4.

Reduction to ODE Eigenvalue Form
The linear stability characteristics of the basic steady state to generally oblique roll disturbances may be investigated by setting, where p 1 and θ 1 are asymptotically small perturbations in the pressure and temperature, respectively. The substitution of the expressions (29) into the governing Equations (15) and (16) gives the following set of linearised perturbation equations We Fourier-decompose the perturbations in the x-and z-directions by setting where φ is the orientation of the roll relative to the x-axis, and where k and λ are once more the wave number of the vortex/cell amd the exponential growth rate, respectively. Expressions (32), when substituted into Equations (30) and (31), reduce the perturbation equations to a set of ordinary differential equations for the reduced disturbance amplitudes, P and Θ: Θ − k 2 Θ = 1 2 Ra (−y 2 + y − 1 6 )ikξΘ sin φ sin α + 1 2 Ra (1 − 2y)Θ cos α − 1 2 (1 − 2y)P + λΘ. (34) The boundary conditions are,

Numerical Results
Equations (33) and (34) subject to boundary conditions (35) were solved using precisely the same methodology as was used in [27,29], which is a variant of the well-known Keller box method. First, Equation (33) was differentiated once to obtain a second order equation in P so that the boundary conditions for both variables, P and Θ are Dirichlet and thereby absolute accuracy is increased. Second, the resulting equations are supplemented by a system obtained by differentiating these equations with respect to k where ∂Ra/∂k = 0 is insisted upon in order to minimise the critical Rayleigh number with respect to the wavenumber. Finally the full system was discretised using second order accurate central differences and solved using the modified Keller box method which is described in [27]. This method allows for the computation of eigenvalues which, in this case, are Ra c , k c , λ i and ∂λ i /∂k. The usual block tridiagonal structure of the iteration matrix is replaced by with an addition row and column of blocks because of the presence of the eigenvalues. A curve-following methodology was also implemented which allows the neutral curve to be followed in Ra-α space for each chosen value of ξ and φ. Thus a fifth 'eigenvalue' is required which measures the length of the curve being followed. The code used for this work was cross-calibrated by solving the same equations using a standard fourth order Runge-Kutta scheme embedded in a boundary value problem solver. Figure 7 displays the variation of the critical value of Ra c cos α with α for a range of roll orientations for ξ = 0.5. For these three-dimensional cases we have chosen to use Ra c cos α rather than Ra c in order to visualize more easily the behaviour of the system and also to avoid Ra c tending towards infinity as the layer tends towards being vertical. In particular we obtain that Ra c = 471.17/ cos α with k c = 4.675. The message which arises from Figure 7 is that longitudinal rolls are preferred over all other orientations when ξ < 1, a general result which accords with the analysis of [27]. In these circumstances there is more resistance to flow in the longitudinal or x-direction than in the other two directions, and it is natural therefore to expect the favoured rolls to have axes in the x-direction.  When ξ > 1 the situation becomes more complicated. At small inclinations, i.e., when α is small, then transverse rolls (φ = 90 • ) form the preferred orientation, but longitudinal rolls are preferred at much higher inclinations. In the inclined Darcy-Bénard problem studied by Storesletten and Tveitereid [32] it was found that there exists a critical inclination below which transverse rolls are preferred and above which longitudinal rolls are preferred. In the more general problem studied by Rees and Postelnicu [27] the transition was found not to be sudden, but, when ξ takes values which are not too large, there exists a smooth transition where the preferred roll orientation varies continuously between φ = 90 • and φ = 0 • as α increases. It is this latter behaviour which we have found to occur in the present problem.
In Figures 8 and 9 we display the variation of Ra cos α with α for a range of values of φ (0 • , 5 • , 10 • , · · · 90 • ) for ξ = 2 and 5, respectively. In Figure 8 we see that the transition between transverse rolls at lower inclinations and longitudinal rolls at higher inclinations takes place close to α = 19.6 • . This transition is smooth, as shown in the second frame of Figure 8, and it takes place very rapidly over an inclination range of roughly 0.4 • . Figure 9 shows the general behaviour of the neutral curves for the more anisotropic case, ξ = 5. Again there exists a smooth transition between transverse and longitudinal rolls, but it now extends over a slightly larger inclination range of 0.54 • and is centred at roughly α = 15.75 • .
In both the cases considered above and for general values of ξ, it is possible to minimise the value of Ra c cos α with respect to φ in order to determine the overall minimum critical Rayleigh number as a function of both the inclination and the anisotropy parameter. The results of this process are displayed in Figure 10 which shows how the critical Rayleigh number, and the corresponding roll orientations and wavenumbers vary, and we have selected six different values of ξ > 1 to illustrate the general behaviour. The curves for ξ = 2 and 5 follow the lower envelope of those shown in Figures 8 and 9, respectively.    It is clear that the transition from transverse rolls (φ = 90 • ) to longitudinal rolls (φ = 0 • ) increases from α = 0 • initially as ξ increases from 1 upwards. When ξ is close to 2 this increase in the transition inclination reverses direction and begins to decrease once more as ξ increases further. The rapid variation in the transition process is captured well in the frames displaying k c and the roll orientation, φ.
Finally, it is of interest to determine an overall summary of the properties of the present layer in terms of the identity of the convection pattern which is to be expected in practice. This is given in Figure 11, which displays the regions in (α/ξ)-space wherein longitudinal, oblique or transverse rolls form the preferred pattern at onset. The rapid transition between transverse and longitudinal rolls as α increases is now very evident due to the closeness of the two curves, and the above-mentioned maximum inclination angle corresponding to ξ 2 is also clear. The curves become extremely close when ξ → 1 + and our numerical data suggests that, Therefore the transition range between these extreme orientations is roughly (ξ − 1) 3/2 as ξ → 1. On the other hand, when ξ takes values which are larger than those shown in Figure 11 then the inclinations at which transitions happen continue to decrease. Some more extreme examples are given in Table 1. These values were obtained by insisting that ∂ 2 Ra c /∂φ 2 = 0, and from this data it is possible to shown (using some error analysis and Richardson extrapolation) that the transition zone (i.e., the range of values of α for which oblique rolls are preferred) occupies the range, when ξ 1.  Figure 11. Showing the different regions in (ξ, α)-space wherein longitudinal rolls, oblique rolls (between the two curves) and transverse rolls form the most unstable planform of convection.

267
On the other hand, when ξ takes values which are larger than those shown in Figure 11 Table 1. These values were obtained by insisting that ∂ 2 Ra c /∂φ 2 = 0, and from this data it is when ξ ≫ 1.
273 Table 1. Transitional values of α when ξ is large; see Figure 11.

Conclusions
The onset of convection has been considered for a porous layer which is inclined to the horizontal, is anisotropic and is subject to a uniform internal heat generation. When the Rayleigh number is sufficiently large, instability occurs, and a convective motion is set up. The critical Rayleigh number at the onset of convection then depends very strongly on the anisotropy ratio as well as the inclination angle. Indeed the preferred orientation of convection rolls at onset can depend quite strongly on these parameters.
We have shown that, when ξ < 1 then convection rolls are longitudinal (i.e., with axes in the x-direction) for all inclinations. Numerically we found that Ra c cos α = 471.17 and critical wavenumber k c = 4.675. However, when the flow is confined to being two dimensional, it has been shown that convection happens for only a finite range of inclinations which is dependent on the anisotropy parameter, ξ. Moreover, when ξ > 0.27843 the neutral curves display an isola point meaning that there a local maximum inclination above which the basic flow is linearly stable.
When ξ > 1 the general situation is more complicated. At small inclinations transverse rolls (φ = 90 • ) form the preferred orientation, but longitudinal rolls (φ = 0 • ) are preferred at much higher inclinations. In general, the critical Rayleigh number increases as the inclination of the layer increases, but the identity of the favoured roll orientation changes from transverse to longitudinal. The transition between these extreme orientations takes place smoothly but rather rapidly as the inclination, α, increases. The associated variation in the critical wavenumbers becomes more extreme as the anisotropy ratio ξ increases. We have also ascertained that the narrow transition region is centred on values of α which are close to zero both as ξ −→ 1 + and as ξ −→ ∞, and that transition region itself tends towards a zero width for each of those limits.