On the stability of parallel flow in a vertical porous layer with annular cross-section

The linear stability of buoyant parallel flow in a vertical porous layer with an annular cross-section is investigated. The vertical cylindrical boundaries are kept at different uniform temperatures and they are assumed to be impermeable. The emergence of linear instability by convection cells is excluded on the basis of a numerical solution of the linearised governing equations. This result extends to the annular geometry the well-known Gill's theorem regarding the impossibility of convective instability in a vertical porous plane slab whose boundaries are impermeable and isothermal with different temperatures. The extension of Gill's theorem to the annular domain is approached numerically by evaluating the growth rate of normal mode perturbations and showing that its sign is negative, which means asymptotic stability of the basic flow. A concurring argument supporting the absence of linear instability arises from the investigation of cases where the impermeability condition at the vertical boundaries is relaxed and a partial permeability is modelled through Robin boundary conditions for the pressure. With partially permeable boundaries, an instability emerges which takes the form of axisymmetric normal modes.


Introduction
In a short paper, Gill (1969) captured the core thermal property of porous insulating slabs employed for the thermal insulation of buildings. Heat transfer in a vertical plane layer of fluid-saturated porous material with impermeable boundaries having different uniform temperatures is always in a conduction regime, no matter how large is the imposed temperature difference. This result is far from being obvious as one assumes that a conduction regime forms when the temperature difference is sufficiently small, while cellular convection flow is expected to arise for larger temperature differences. Incidentally, this is precisely what happens if we have a vertical fluid layer instead of a saturated porous medium (Vest and Arpaci, 1969).
The impossibility of a convective regime, with an enhanced heat transfer rate compared to the conduction regime, means that a vertical porous slab saturated by air provides a much more efficient insulation than does a vertical air gap free of porous material. The core of Gill's theorem (Gill, 1969) is the linear stability analysis of the basic conduction regime in the vertical porous slab in which he used an integral analysis to show that the exponential growth rate of small disturbances always remains negative. This work has offered a fertile ground for further developments; examples include the detailed analysis of the growth rate of perturbations (Rees, 1988;Lewis et al., 1995), the inclusions of other effects (Kwok and Chen, 1987;Rees, 2011) and the extension to the nonlinear regime (Straughan, 1988). An important feature of Gill's theorem for the absence of thermoconvective instability in a vertical porous slab is that it relies on the hypothesis that the bounding planes are impermeable. If the boundaries are modelled as permeable, then an instability occurs (Barletta, 2015).
The aim of this paper is to investigate the validity of Gill's theorem when its formulation is adapted to a vertical annular layer of saturated porous material. With reference to thermal insulation techniques, this result may be interesting when we focus, say, on heat transfer from a hot fluid flowing in a vertical round pipe. In fact, insulation of the hot pipe can be done by cladding a low conductivity porous layer around the pipe wall.
The main difference with respect to the plane slab examined by Gill (1969) is that the case of an annular porous layer does not allow for a simple rigorous proof of stability. In the present paper we have adopted a strategy based on the numerical evaluation of the perturbation growth rate in order to test the stability of the flow, and it is concluded that the conduction regime is always stable. This conclusion is further validated by considering cases when the impermeability of the boundaries is made imperfect. This imperfection is monitored by means of a dimensionless parameter τ where τ = 0 corresponds to impermeable boundaries. It was found that instability in the form of an axisymmetric mode arises whenever τ > 0, but that the critical Rayleigh number becomes infinite as τ → 0.

Mathematical model
We aim to model a vertical porous layer with cylindrical shape and annular cross-section, bounded by an internal radius r 1 and an external radius r 2 . The radial boundaries r = r 1 and r = r 2 are impermeable and isothermal with temperatures T 1 and T 2 , respectively.

Governing equations
By adopting the Oberbeck-Boussinesq approximation and Darcy's law (Nield and Bejan, 2017), the local balance equations for mass, momentum and energy can be expressed in a dimensionless form as with the boundary conditions r = 1 : Here, u = (u, v, w) is the seepage velocity with u, v and w denoting the radial, angular and axial velocity components, p is the local difference between the pressure and the hydrostatic pressure, T is the temperature, t is time and e z is the unit vector along the axial z axis. The cylindrical coordinates (r, φ, z) are employed. The dimensionless parameters R, s and ζ are defined as where g is the modulus of the gravitational acceleration g = −g e z , β is the coefficient of thermal expansion, K is the permeability, ν is the kinematic viscosity, α is the average thermal diffusivity and T 0 is the average temperature employed as the reference value within the Oberbeck-Boussinesq approximation and the definition of the buoyancy force. The scaling employed to exploit the dimensionless equations (1) and the boundary conditions (2) are given by where µ is the dynamic viscosity and we denoted with σ the ratio between the volumetric heat capacity of the saturated porous medium and that of the fluid.

Buoyant parallel flow
There exists a stationary basic solution of Eqs.
(1) and (2) characterised by a parallel velocity field, a vanishing flow rate and a purely radial temperature gradient. This is given by where the subscript "b" means basic state.

Stability of the basic state
In order to examine whether the basic flow state defined by Eq. (5) is stable or not, we first rewrite Eqs. (1) and (2) according to a pressure-temperature formulation, Then, we perturb the basic state (5) by small amplitude normal modes, where "c.c." is a shorthand for complex conjugate, ε is a perturbation parameter such that |ε| ≪ 1, m is a non-negative integer, k is the wavenumber, ω is the angular frequency, η is the time growth rate, while f and h are radial amplitude functions. By substituting Eq. (7) into Eqs. (6), by taking into account Eq. (5), and by neglecting all terms O(ε 2 ), we obtain Equations (8) define an eigenvalue problem, where the eigenvalue is the complex quantity η − i ω, while the pair (f, h) is the eigenfunction. The solution is sought for fixed input values of (m, k, R, s). The real part of the eigenvalue, i.e. the growth rate η, is an important parameter as it allows one to detect the linear stability (η 0) or instability (η > 0) of the basic solution. We could not find a rigorous proof that the solution of Eqs. (8) yields η ≤ 0, as for the stability theorem proved by Gill (1969) for a plane slab. However, there are some results that can be proved and which are the starting point for a numerical solution of Eqs. (8).
3.1 Asymptotic case R → 0 The limit of a vanishing Rayleigh number is a case where no buoyant flow exists and the fluid is isothermal in the basic state (5). Then, we expect η ≤ 0 as the basic state will be stable. However, it is interesting to determine the value of η independently of the information that its sign cannot be positive. By assuming R = 0, we multiply Eq. (8a) by rf , where the bar over the symbol denotes complex conjugation. Then, we integrate by parts over the interval 1 r s by taking into account the boundary conditions (8c). The result is Equation (9) can be satisfied only with df /dr = 0, if m = 0 and k = 0, or with f = 0 in every other conditions. In either cases, the last term on the left hand side of Eq. (8b) is zero. Thus, if we multiply Eq. (8b) by rh and we integrate by parts over the interval 1 r s by taking into account the boundary conditions (8c), we obtain Given that h cannot be identically zero, Eq. (10) implies that ω = 0. In addition, we can conclude that Eq. (10) can be satisfied only if η < 0. In order to determine η for given s and k, we rewrite Eqs. (8b) and (8c) as The solution of Eq. (11) can be expressed in terms of modified Bessel functions of order m, namely  15) provided that the dispersion relation, is satisfied. Here, C is an arbitrary constant, while I m and K m are the modified Bessel functions of first and second kind, respectively. Thus, by employing the properties of Bessel functions (see, for instance, chapter 10 of Olver et al., 2010), we can evaluate the growth rate η as where γ is a root of For each choice of (m, s), we are interested in detecting through Eqs. (14) and (15) the smallest γ as it yields the less stable condition, i.e. that where η is at its largest. Plots of the value of γ versus s for m = 0 to 3 are reported in Fig. 2. We can draw some afterthoughts: the growth rate η for R = 0 is always strictly negative; the value of η decreases with m and k while it increases with s. The former feature means that the basic state is asymptotically stable, according to our linear analysis. However, the most important fact is that Eq. (14) provides an accurate starting point for the computation of the eigenvalue η − i ω when the Rayleigh number gradually increases from zero, for fixed (m, k, s).

Asymptotic case k → 0
The case of an infinitely large wavelength, or an infinitely small wavenumber, can be tackled by a treatment similar to that presented in Section 3.1. In fact, since R appears in Eqs. (8) only through the expression k R, the limit k → 0 is just a sub-case of that analysed in Section 3.1. Then, we can just use the results drawn in that section and specialise them by setting a zero wavenumber.
In physical terms, we infer that axially invariant, or z independent, modes cannot activate the instability for every value of R.

Computation of the eigenvalue
The solution of the differential eigenvalue problem (8) can be approached numerically by employing the shooting method (see, for instance, chapter 9 of Straughan, 2008, or chapter 10 of Barletta, 2019). The main stages of the numerical procedure are the following: 1. We solve Eqs. (8a) and (8b) as an initial value problem where the initial conditions are those imposed at r = 1. They are given by Eq. (8c). However, as they stand, they are insufficient to match the differential order of Eqs. (8a) and (8b). Thus, we need two additional initial conditions. One of them relies on the scale invariance of the homogeneous problem (8), which can be broken by imposing dh/dr = 1 at r = 1. The second one is the statement f (1) = ξ 1 + i ξ 2 which does not imply any loss of generality inasmuch as ξ 1 and ξ 2 are general real parameters.
2. We determine the eigenvalue η − i ω together with the parameters ξ 1 and ξ 2 by employing a root finding algorithm, such as the Newton-Raphson method, applied to the target conditions expressed by Eq. (8c), at r = s, namely df /dr = 0 and h = 0. In fact, since the eigenfunctions f and h are complex-valued, such target conditions entail four different real constraints.
Both stages are implemented by employing the Mathematica software ( c Wolfram Research) and, in particular, we use the built-in functions NDSolve and FindRoot.
The shooting method can be utilised for every assignment of the input data (m, k, R, s). The analysis carried out in Section 3.1 establishes the values of η, ω, ξ 1 and ξ 2 for R = 0, For every assigned (m, k, s), we gradually increase the value of R above zero step-by-step. At each step, we call the root finding algorithm by initialising the search routine with the solution data found at the previous step. The step size for the increment of R is dynamically adapted according to the extent of the eigenvalue change.

Analysis of the growth rate
The numerical computation of the complex eigenvalue η − i ω, for assigned s, k, m and R, allows one to track the growth rate of the normal modes. The sign of the parameter η yields the most important information regarding the stable/unstable behaviour of the system. The numerical data collected by varying the governing parameters support the conclusion that η < 0 in every case. Thus, we reach the conclusion that no instability is possible, exactly as in the case of the plane slab examined by Gill (1969). Figure 3 displays the growth rate η versus R for three different aspect ratios, s = 1.5, 2 and 3. The numerical data are relative to m = 0, 1 and 2. These three types of normal mode are identified by solid lines, black dotted lines and grey dotted lines, respectively. There are some distinctive features that may be pinpointed at a glance. The growth rate η is a negative, monotonically decreasing, function of R for all the values of s, k and m considered in Fig. 3. The axisymmetric modes (m = 0) yield the less stable conditions, i.e. those leading to the largest growth rate, only when R is sufficiently small. At larger values of R, the m = 1 or the m = 2 modes prevail. The selection of the less stable modes largely depends on the value of k. We note that, with k = 0.1, η displays a weak dependence on k. This is expected as in Section 3.2 we pointed out that η is independent of R when k → 0. On the other hand, η is poorly influenced by the value of m when k is larger as it becomes apparent in Fig. 3 with k = 2.
All the data reported in Fig. 3 show that η < 0 and that the growth rate is a decreasing function of R and of k, while the influence of the angular number changes when k increases. As anticipated, the physical information gathered from this analysis is that no linear perturbation mode grows in time, so that the basic flow is always linearly stable.

Changing the impermeability condition at the boundaries
One method of validating the conclusion that the basic flow (5) is stable is by relaxing the impermeability boundary conditions at the inner and outer radii of the cylinder and then investigating the limit of zero permeability. This may be achieved by replacing Eq. (8c) by The parameter τ , which is assumed to be non-negative, marks the departure from impermeability. Obviously, the boundaries return to being perfectly impermeable only when τ → 0, a limit where Eqs. (8c) and (17) coincide. Equation (17) follows with Robin boundary conditions for the pressure. Such conditions physically mean that the normal component of the seepage velocity at the boundary is proportional to the pressure difference between the boundary and the external environment. When τ = 0 our computations show that the growth rate η may undergo a transition from negative to positive when R becomes sufficiently large, with a neutral stability curve (η = 0) delimiting the boundary between the regions of parametric stability and instability.
This behaviour is displayed in Fig. 4, where neutral stability curves are drawn for the case s = 2. It is evident that decreasing the value of τ starting from τ = 10 turns into a stabilization of the flow as the neutral stability curves move up and left in the (k, R) plane. The wide gap, for m = 0, between the minimum of the τ = 1 curve and that of the τ = 0.5 curve, if compared to that between the curves for τ = 10 and τ = 4, suggests a very steep variation of the neutral stability curve as τ becomes smaller and smaller. Just the same conclusion is drawn by comparing the curves with τ = 4 and τ = 2, relative to the m = 1 modes. Figure 4 also indicates that, for a given τ , the m = 0 modes are those activating the instability first or, equivalently, that such modes are the most unstable. Figure 5 allows one to trace the trend of the critical value R c versus τ for two sample cases, corresponding to the ratios s = 2 and s = 3. First of all, it must be noted that, in both cases, the most unstable branch is that relative to the axisymmetric normal modes, i.e. those with m = 0. The modes with m = 1 provide a higher branch of instability, while those with m = 2 yield even larger values of R c and are not even visible within the vertical range of Fig. 5. The behaviour reported in Fig. 4, suggesting a steep increase in the critical value of R, is confirmed quite clearly by the curves displayed in Fig. 5. In fact, it is quite evident that the critical value of R tends to infinity when τ → 0. In particular, the numerical values for R c for the case, s = 2 and m = 0, satisfy the relation R c ∼ 267.671τ −1/2 + 52.55τ 1/2 for τ ≪ 1.
When τ = 0.05 this formula yields R c = 1208.81 whereas the present numerical calculation gives R c = 1208.83. Such a result shows clearly that R c → ∞ as τ → 0, and it yields an indirect supporting argument for the conclusions drawn in Section 3.4. The basic buoyant flow described in Section 2.2 is always linearly stable. We have in fact some further information, as the analysis illustrated in Figs. 4 and 5 conveys the awareness that breaking the impermeability condition as it is given by Eq. (8c) induces the instability of the basic parallel flow.

Conclusions
The validity of Gill's theorem (Gill, 1969), originally formulated for a vertical plane porous slab subjected to a boundary temperature difference, is examined for the case of an annular porous layer. More precisely, the stability of the basic stationary conduction regime, implying a vertical buoyant flow through the porous layer, is tested by the introduction of linear normal mode perturbations. Analytical results have been discussed for the special cases where either the Rayleigh number tends to zero, or the wavenumber of the normal modes tends to zero.
In the absence of a general rigorous proof of stability, a numerical solution method has been employed to compute the growth rate of the perturbation modes. In all the tested cases the growth rate turned out to be negative, meaning that the basic flow state is always stable. An additional argument supporting this conclusion has been provided by examining the effects of relaxing the impermeability boundary conditions. Robin pressure conditions, as mediated by the dimensionless parameter τ , are introduced instead of Neumann conditions. This change turned out to activate an instability which, however, disappears in the limiting case where the Robin conditions for the pressure tend towards the Neumann boundary conditions that correspond to impermeable boundaries.