The Effect of Conducting Boundaries on the Onset of Convection in a Porous Layer Which is Heated from Below by Internal Heating

The onset of convection in a porous layer heated from below is considered, and we determine how the presence of two solid but heat-conducting bounding plates of finite thickness alters the manner in which convection ensues. Heat is generated by the lower plate (with an insulating lower boundary), but the upper one is passive with a fixed upper boundary temperature. It is shown that this composite layer may mimic in turn one of the three different types of classical single-layer onset problems which are well-known in the literature. The type which is selected (or indeed whether it corresponds to a transitional case) depends quite critically on the precise values of the relative thickness of the solid layers and their conductivity ratio. It is also shown that care needs to be taken over declaring that the solid plates are thin: extreme values of the conductivity ratio can yield a stability criterion which appears to be different from that suggested by the imposed boundary conditions.


Introduction
The modelling of convective instability in porous layers now has a very strong heritage. The foundational works are by Lapwood (1948) and Horton and Rogers (1945), and these determined the criterion for the onset of convective instability in a horizontal porous layer which is heated from below. In addition, the upper and lower boundaries were considered to be impermeable and subject to constant fixed temperatures, and the porous medium was isotropic and homogeneous. Other types of boundary condition, such as a constant heat flux or constant pressure, may be applied in a variety of different ways; Table 6.1 in Nield and Bejan (2013) gives a comprehensive account of the critical values of the Rayleigh number and wavenumber for these different combinations. The effects of confinement, such as convection in a cuboidal cavity and in a cylindrical cavity (with a circular horizontal planform), were considered by Beck (1972) and Bau and Torrance (1981), respectively. Unlike the corresponding pure fluid configurations, there exist aspect ratios of the cuboid and radii of the cylinder for which the critical Darcy-Rayleigh number is identical with that of the infinite unconfined layer.
The literature is now replete with the study of how the onset and subsequent development of convection is altered by the presence of different effects, such as Brinkman effects, non-Boussinesq effects, the presence of imperfections such as sinusoidally varying boundary temperatures, throughflow effects, rotation, layering, anisotropy, internal heating, local thermal nonequilibrium, more than one diffusing species and unsteady basic states; note that this list, while long, is not comprehensive. Although some work has been reported on configurations where two porous sublayers are separated by a solid intermediate sublayer (see Postelnicu 1999;Jang and Tsai 1988;Patil and Rees 2014;Straughan 2014 for systems of this general type), we are presently interested three-layer systems where a porous sublayer is bounded by two heat-conducting but impermeable solid sublayers. Mojtabi and Rees (2011) considered such a three-layer system subject to constant heat flux boundary conditions, and this was motivated by the need to model an experimental system where solid boundaries do not have zero thickness and idealised boundary conditions cannot be applied in practice. It was found that the critical parameters (namely the Darcy-Rayleigh number and wavenumber) and the shape of the neutral curve depend quite strongly on the thickness of the solid boundaries and their conductivity relative to that of the porous sublayer. For example, when the solid sublayers are sufficiently highly conducting, the full system may be made to resemble a single porous layer system with constant temperature boundary conditions. In addition, there is a locus in parameter space where the neutral curve has a quartic minimum at zero wavenumber, and this locus separates the regime within which the more usual parabolic minimum at nonzero wavenumbers exists from that where the parabolic minimum is at zero wavenumber. A second study was undertaken by Rees and Mojtabi (2011) where constant temperature boundary conditions were applied. There are many qualitative similarities between this study and the previous one by the same authors except that there is no quartic minimum, and the zero critical wavenumber limit is obtained as the solid conductivities tend to zero. Rees and Mojtabi (2011) also provided a weakly nonlinear analysis to determine regimes within which postcritical convection takes the form of rolls or patterns with square planform; this was an extension of the work of Riahi (1983) who considered solid layers of infinite thickness. A paper by Rees and Mojtabi (2013) also considered the effect on the onset of convection of a horizontal pressure gradient. If was found that there is a thermal drag caused by the stationary solid boundaries which reduces the velocity of the cells along the layer.
In the present paper, we consider the linear stability of a configuration which has some similarities to those of Mojtabi and Rees (2011) and Rees and Mojtabi (2011). Here, the lower solid sublayer is insulated on its lower surface, but its interior is subject to a uniform rate of heat generation, while the upper sublayer is subject to a constant temperature. It is shown in the Appendix that this may be mapped onto a system where the lower solid sublayer is subject to a constant heat flux but without heat generation within the sublayer. While this may seem to be a simple intermediate case between those of Mojtabi and Rees (2011) and Rees and Mojtabi (2011), three different limiting cases for the critical parameters are found in different asymptotic regions of parameter space, and this is a novel finding. Fig. 1 Definition sketch of the configuration being studied. The porous layer is sandwiched between two impermeable but thermally conducting layers

Governing Equations and Basic Solution
We consider a three-layer system as shown in Fig. 1. It consists of a central uniform porous layer (layer 2) within which convective flow may occur, and this is embedded between two thermally conducting but impermeable layers. These solid layers have identical properties, namely thickness, conductivity and heat capacity. A uniformly cold temperature is imposed on the upper surface of the upper solid layer (layer 3), while the lower layer (layer 1) is subject to uniform Joule heating but its lower surface is insulated. Thus, all the heat is generated within the lowest sublayer and passes out through the upper surface, and a potentially unstable temperature gradient is formed. The origin of the coordinate system is located at the bottom of the porous layer withx being the horizontal coordinate andẑ the vertical coordinate. The full set of governing equations are as follows. For layer 1 (−h ≤ẑ ≤ 0): for the porous layer (0 ≤ z ≤ H ): and for the upper layer (H ≤ z ≤ H + h): In these equations, the numerical subscripts refer to the appropriate sublayer, and T ref is the basic temperature of the interface between the porous sublayer and the lower solid sublayer for the purely conducting situation. The boundary and interface conditions are given by, The basic state consists of no flow (û =ŵ = 0) and a steadyẑ-dependent temperature. It is a straightforward if tedious matter to determine the conduction profiles in each layer: The temperature drop across the porous layer will be used as the temperature scale for the purposes of nondimensionalisation. This value is and we see that the temperature, T ref , of the lower interface temperature is given by, We may now nondimensionalise by introducing the following scalings, Thus, Eqs. (2) and (3) become, Here, the Darcy-Rayleigh number, Ra, is defined as, where T is given in terms of q in Eq. (8). The three heat transport equations become, where are the respective conductivity and thickness ratios between the solid and porous sublayers. The boundary and interface conditions given in Eq. (6) become, The apparently unusual boundary condition for θ 3 at z = 1 + δ has been defined in such a way that the temperature field varies between 0 and −1 within the porous layer. Thus, the basic solution is given by While there is a unit change in temperature across the porous layer, the overall temperature drop over the whole system is 1+3δ/2d, and this increases as either the thickness of the solid sublayers increase or their conductivity decreases. Five representative profiles of the basic temperature field are given in Fig. 2 for the case δ = 0.5; profiles for other values of δ display the same qualititative features. The insulating condition at z = −1 − δ is evident because of the zero gradient. By design, the temperature profile in the porous layer is independent of the values of both d and δ. Figure 2 shows that there is little variation in the temperature within the solid sublayers when d is large. This happens quite naturally because of the need to conserve the heat flux across the interface. And then, when d takes small values, there is a very large variation in temperature within the solid sublayers.

Perturbation Analysis
The basic state consists of the temperature profiles given in Eq. (19) with no flow. Since Squire's theorem may be shown easily to apply, i.e. that all three-dimensional disturbances may be decomposed into a sum or integral of two-dimensional ones, we need to consider only two-dimensional disturbances. Therefore, we adopt a streamfunction/temperature formulation using, u = −ψ z and w = ψ x , which guarantees that the equation of continuity, Eq. (11), is satisfied. Although the heat transport equations in the two solid sublayers remain unchanged, Eqs. (12) and (15) take the following forms, We may now introduce a disturbance of infinitesimally small amplitude using, where ψ basic = 0 and the basic temperature fields are given in Eq. (19). On neglecting powers of the disturbances, the linearised stability equations take the form, These equations are to be solved subject to the homogeneous form of the boundary conditions given in (18), with the addition of = 0 on both z = 0 and z = 1.
It may also be shown easily that the principle of exchange of stabilities applies, and therefore, the onset of convection is stationary and the time derivatives in Eqs. (23)-(26) may be neglected. If the domain is of infinite horizontal extent, then we may Fourier-decompose disturbances into independent monochromatic cells. Thus, we set, where k is the horizontal wavenumber. Equations (23)-(26) reduce to where primes denote derivatives with respect to z. The equations for 1 and 3 may be solved analytically and converted into suitable boundary conditions for 2 upon application of the interface conditions, as follows. For Layer 1, Eq. (29) has the solution, The application of the boundary condition, 1 (−δ) = 0, means that this solution may now be written in the form, Application of the interface conditions at z = 0, namely that 1 = 2 and d 1 = 2 allows us to eliminate the arbitrary constant, E, to yield which is regarded as a boundary condition for 2 that models perfectly the variation of 1 in the lower sublayer. A similar analysis for 3 yields The linear stability problem has now been reduced to solving Eqs. (30) and (31) subject to the two thermal boundary conditions given in (35) and (36), together with those in Eq. (27).

Dispersion Relation
Equations (30) and (31) form a fourth-order homogeneous system, and there are four independent solutions. These may be written in the following form, where A, B, C and D are arbitary constants, and where The Darcy-Rayleigh number is found always to lie above the curve, Ra = k 2 , and therefore, the value of σ is always purely real.
If we now apply the boundary conditions, (27), (35) and (36), then we obtain the following matrix/vector system for the four arbitrary constants: ⎛ This system has nonzero solutions only when the determinant of the matrix is zero. After some manipulation, we find that this determinantal condition is equivalent to It is impossible to solve this equation analytically to give Ra in terms of k, d and δ, and therefore, we resort to numerical means. A simple Newton-Raphson may be used easily to determine neutral curves as k varies for fixed values of d and δ. A slightly more complicated version was also written in order to find critical values, i.e. those values of k and Ra at the lowest point of the neutral curve. The manner in which this was undertaken is precisely as described in Rees and Mojtabi (2011) and Rees and Genç (2011) and is therefore omitted here. The codes were written in quadruple precision in Fortran90 because it was found that this amount of precision was required to obtain accurate computations when d and δ take values which are as small as 10 −10 and as large as 10 10 .

Results and Discussion
In what follows there are three sets of critical data which are useful references. For the classical Darcy-Bénard layer, the critical values of the Darcy-Rayleigh number and wavenumber are, where the appelation, CT/CT, refers to the fact that both the upper and lower boundary conditions correspond to a constant temperature. When both surfaces are subject to a constant heat flux, we have CHF/CHF: Ra c = 12, k c = 0.

Onset Profiles
We begin by showing, in Fig. 3, a selection of onset profiles ( , ) so that the effect of variations in the two nondimensional parameters, d and δ, on those profiles may be seen. In both frames the temperature profile, , takes positive values, while the streamfunction disturbance profile, , takes negative values. The profiles have been scaled so that the maximum magnitude is precisely unity for easy comparison. Given our definition of the streamfunction, the values of turn out conveniently to be negative when is positive. Although we have set k = 2 in all cases, the qualitative trends shown in Fig. 3 do not change for different wavenumbers. Figure 3a shows the effect of variations in d when δ = 0.4. When d takes large values, the solid sublayers are highly conducting and we see therefore that the disturbance temperature is almost perfectly constant; this situation is close to the CT/CT single-layer case mentioned above. At the opposite extreme, d 1, the solid layers are poor conductors, and this gives rise to relatively large variations in within the solid regions, as was seen for the basic conduction profiles, and only small variations within the porous sublayer. This is approximately the CHF/CHF case mentioned above. It is also interesting to note that the streamfunction profile hardly varies at all over the whole range of values of d. Figure 3b shows the effect of changing the depth, δ, of the solid sublayers. We have chosen a fairly representative value for the conductivity ratio, d = 5, and this case shows an interesting transition in the profiles from one which is close to being of CT/CT type when δ is large (i.e. thick solid sublayers) to one which is almost indistinguishable from a CT/CHF type when δ is small. In such cases, the very thin solid sublayers appear to be yield stability properties and profiles which are identical to cases where the layers are entirely absent. We will see later that this initial hypothesis turns out not to be that straightforward in general, and attention must be paid to the magnitude of δ in comparison with that of d. Figure 4 shows a small selection of neutral stability curves for some typical values of the parameters, d and δ, where each lies in the range, 10 −3 -10 3 . These curves show the standard unimodal form with a clearly identifiable minimum value; numerical investigations indicate that unimodal curves are always obtained for the present configuration. The uppermost curve (d = 10 3 ) in all four cases is almost exactly identical to that for the classical Darcy-Bénard problem (CT/CT) for which the critical values are Ra c = 4π 2 and k c = π. This is because the solid layers are highly conducting which means that temperature variations are rapidly damped, and therefore the interfaces act as perfectly conducting boundary conditions. That said, the critical values for the d = 10 3 curve when δ = 10 −3 are very slightly smaller than the classical values, but this is difficult to see graphically. On the other hand, the lowest curves (d = 10 −3 ) for all four cases change quite substantially as δ increases, and the critical values decrease quite markedly. Indeed, when δ is large, the highly insulating bounding layers mimic the CHF/CHF porous layer for which Ra c = 12 and k = 0.

Neutral Curves
On the other hand, from these graphs, it is possible to find out how quickly the neutral curves change as d varies. It is clear that there is almost no change in the location of the neutral curve when d increases from 10 −2 to 10 2 when δ = 10 −3 . But when δ = 10 3 , there is a rapid change as d increases from 10 −1 to just below 10 1 ; outside of these ranges there is little variation. While we have gleaned this information from Fig. 4, it is clear that we need to consider the variations in the critical parameters themselves, rather than by observing further examples of neutral curves.
An extensive set of computations has therefore been undertaken to determine the critical values of Ra and k in the ranges 10 −6 ≤ δ ≤ 10 6 and 10 −10 ≤ d ≤ 10 10 , and these are shown in Fig. 5. Referring to both frames in this figure, we see that the general behaviour of the critical values is that they pass from the CHF/CHF case via the CT/CHF case to the CT/CT case as d increases from exceptionally small values to exceptionally large values. When the solid sublayers are comparable to the thickness of the porous layer, or larger, then there is a direct transition from CHF/CHF to CT/CT which is quite quick near d = 10 0 when d is increasing. But when δ < 10 −3 , there is range centred on d = 10 0 where the composite layer resembles the CT/CHF layer with very high accuracy. Figure 6 shows a close-up image of the central region of the wavenumber curves in Fig. 5. This shows that, when δ = 10 −5 , the composite layer mimics a CT/CHF single layer when 10 −2 < δ < 10 2 , and when δ = 10 −6 , this happens when 10 −3 < δ < 10 3 .
In the close-up view, we also see that there is an undershoot in the critical wavenumber as d decreases whenever δ is sufficiently small. This presages the passage towards the CT/CHF equivalence. It is unclear what the physical reason for this might be, but a similar behaviour does not occur for the critical Darcy-Rayleigh number.

Summary
An alternative view of the curves shown in Fig. 5 is provided by Fig. 7, which shows contour representations of both Ra c and k c in the ranges, 10 −10 ≤ d, δ ≤ 10 10 . A grid with 100 inter- Critical values of Ra and k as a function of log 10 d for δ = 10 −6 , 10 −5 , 10 −4 , 10 −3 , 10 −2 , 10 −1 , 10 0 (dotted) and 10 1 vals in both directions was used for these figures. We now see clearly how the (d, δ) plane is divided into three very distinct regions with smooth transitions between these regions. Generally, the right-hand side corresponds to being equivalent to a CT/CT single layer, the left-hand Close-up view of the critical values of k as a function of log 10 d for 10 −5 , 10 −4 , 10 −3 , 10 −2 , 10 −1 , 10 0 (dotted) and 10 1 side to a CHF/CHF layer, while the triangular wedge is a CT/CHF configuration. Therefore, our composite layer is clearly capable of supporting three different types of behaviour.
Earlier, a statement was made that what constitutes a narrow solid layer is not a straightforward concept. In Fig. 6, if we choose δ = 10 −5 to represent a very thin pair of solid sublayers, then we should expect this to represent the single-layer CT/CHF case simply because 10 −5 is exceptionally small. However, the critical values correspond to the single-layer CT/CHF critical values only when d lies in the range 10 −4 < d < 10 4 . Outside of this range, we obtain the CHF/CHF case when d < 10 −9 and the CT/CT case when d > 10 5 . These extreme values for d overrule our standard intuition of what 'thin' means in the context of the stability. Of course, if one takes the alternative viewpoint that d is fixed, but δ decreases, then one will always eventually obtain a CT/CHF configuration.
It is clear from the two subfigures in Fig. 7 that the transitions between the three main regions take place at d ∼ 10 0 10 0 < δ, d > 10 0 : dδ ∼ 1, d < 10 0 : d ∼ δ. 123 triangular region at the bottom of that contour plot and which are close to d = 10 0 are also within 10 −6 of the CT/CHF value in Eq. (44). This corresponds almost perfectly with the curves given in Fig. 6.

Conclusions
We have conducted a linear stability analysis to determine the effect on stability criterion of having two solid but identical bounding sublayers above and below a saturated porous layer where the lower sublayer generates heat internally. The equations were nondimensionlised in such a way that the Darcy-Rayleigh number is based on the temperature drop across the porous layer. There are two parameters, namely the conductivity ratio between the solid and porous sublayers, d and the relative thickness of the bounding plates, δ. Therefore, we have been able to provide a comprehensive account of the stability properties over a very wide range of values of the parameters. The stability criteria were determined by finding the zeros of a dispersion relation and by minimisation over the wavenumber.
It has been found that the addition of conducting layers above and below the porous layer allows the stability criteria to vary quite substantially and, in particular, to allow this composite layer to mimic one of three standard single-layer configurations, namely the CT/CT, CT/CHF and the CHF/CHF cases. This type of behaviour is quite different from that of Mojtabi and Rees (2011) where the composite layer which they studied only displayed two types of limiting behaviour.
From the practical point of view of constructing an experimental configuration, we have shown that the presence of a physically thin pair of solid bounding plates does not necessarily guarantee that the stability properties correspond to the type (namely CT/CHF) corresponding to the boundary conditions which are applied. Rather, when δ is very small, only a certain range of values of d will yield the CT/CHF onset criterion.
We do not know at this stage whether the rolls which we have considered will form the stable mode of convection at values of the Darcy-Rayleigh number which are just above the critical values. Given the conclusions of Riahi (1983) and Rees and Mojtabi (2011), we feel that it is quite likely that square cells are likely to be preferred when the conductivity of the bounding sublayers is sufficiently small, but this aspect would need to be considered quantitatively elsewhere. However, it is very likely indeed that both the roll and square cell patterns will arise from supercritical bifurcations, and therefore, all the results presented here will also be valid should square cells be the preferred pattern post-onset.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.