Stability of a Buoyant Oldroyd-B Flow Saturating a Vertical Porous Layer with Open Boundaries

The performance of several engineering applications are strictly connected to the rheology of the working fluids and the Oldroyd-B model is widely employed to describe a linear viscoelastic behaviour. In the present paper, a buoyant Oldroyd-B flow in a vertical porous layer with permeable and isothermal boundaries is investigated. Seepage flow is modelled through an extended version of Darcy’s law which accounts for the Oldroyd-B rheology. The basic stationary flow is parallel to the vertical axis and describes a single-cell pattern where the cell has an infinite height. A linear stability analysis of such a basic flow is carried out to determine the onset conditions for a multicellular pattern. This analysis is performed numerically by employing the shooting method. The neutral stability curves and the values of the critical Rayleigh number are evaluated for different retardation time and relaxation time characteristics of the fluid. The study highlights the extent to which the viscoelasticity has a destabilising effect on the buoyant flow. For the limiting case of a Newtonian fluid, the known results available in the literature are recovered, namely a critical value of the Darcy–Rayleigh number equal to 197.081 and a corresponding critical wavenumber of 1.05950.


Introduction
The convective instability in a plane vertical porous layer saturated by a fluid is a topic of great interest for its applications spanning from geophysical systems to building insulation. Although the analysis of geophysical systems typically involves groundwater or hydrocarbons, the study of thermal insulation panels features air or another gas as a working fluid. Hydrocarbons saturating a porous medium may occasionally display a non-Newtonian viscoelastic behaviour. This circumstance will be that envisaged in the present study.
The classical paper by Gill [1] offered a rigorous mathematical proof that a vertical porous layer, saturated by a Newtonian fluid endowed with impermeable boundaries, kept at different uniform temperatures, displays a stationary conduction regime with a stable parallel vertical buoyant flow, regardless of the temperature gap between the boundaries. In other words, such basic buoyant flow is always stable for every value of the Rayleigh number. The linear stability analysis carried out by Gill [1] is relative to a Newtonian fluid subject to Darcy's law. In the context of Newtonian fluids, other authors further developed Gill's important result by including other features such as the non-linearity of the perturbation dynamics or other momentum transfer models of the seepage flow in the porous medium [2][3][4][5][6][7]. Further investigations provided evidence that the stability theorem proved by Gill does not hold when the hypothesis of impermeable boundaries for the porous layer is released [8][9][10][11][12]. As proved in these studies, with perfectly or partly permeable boundaries, the onset of convective instability is possible for sufficiently large values of the Rayleigh number.
Stability analyses of the buoyant flow in a vertical porous layer envisage a non-Newtonian saturating fluid only in a few instances [13][14][15]. Barletta and Alves [13] considered a power-law fluid, while Shankar and Shivakumara [14,15] considered a viscoelastic fluid described according to the Oldroyd-B model as developed, for saturated porous media, by Khuzhayorov et al. [16] and employed by other authors [17,18]. However, the non-Newtonian versions of the stability analysis of Gill's system developed by Barletta and Alves [13] and by Shankar and Shivakumara [14,15] are relative to a layer with impermeable boundaries.
Following the path devised by Barletta [8], the aim of this paper is to relax the assumption of impermeable boundaries and reconsider the linear stability analysis of the vertical buoyant flow in a vertical layer saturated by a viscoelastic Oldroyd-B fluid. The special case analysed by Barletta [8] is retrieved when the relaxation time of the viscoelastic fluid coincides with the retardation time. On the other hand, the viscoelastic behaviour shows up when such times are different. In this case, the threshold to convective instability is affected by the viscoelastic rheology of the fluid. The neutral stability condition and the critical Rayleigh number will be obtained for different values of the characteristic governing parameters of viscoelasticity.

Mathematical Model
A vertical porous layer saturated by a viscoelastic fluid is studied. The layer is infinitely wide in the y and z directions and has thickness L in the horizontal direction x. A sketch of the layer is presented in Figure 1. The Oldroyd-B model is employed to describe the viscoelasticity. Darcy's law for Oldroyd-B type of fluids [16] is written where τ 1 and τ 2 are two characteristic relaxation/retardation times, respectively, with τ 2 τ 1 . The buoyancy force is modelled by the Oberbeck-Boussinesq approximation, x = (x, y, z) is the Cartesian position vector, u = (u, v, w) is the filtration velocity vector, T is the temperature and g is the gravitational acceleration vector (opposite to the z direction). The energy balance equation employed to model the heat transfer is the convection/conduction equation where no source or sink term is considered.
The vertical boundaries of the layer are permeable and the external environments in the regions x < −1/2 and x > 1/2 are considered as isothermal fluid reservoirs in a motionless state, so that the pressure distribution along the boundaries is purely hydrostatic. Since the external fluid reservoirs are kept at different uniform temperatures T 1 and T 2 , where T 2 > T 1 , the layer is subject to side heating, as depicted in Figure 1.
Thus, the governing equations are given by where g and e z are the modulus of g and the unit vector along the z direction, respectively, while p 0 is the reference hydrostatic pressure of the external environments. In Equation (2) µ is the dynamic viscosity of the fluid, K is the permeability of the porous medium, ρ f is the fluid density evaluated at the reference temperature T 0 , β is the thermal expansion coefficient of the fluid, σ is the heat capacity ratio, χ is the ratio between the effective thermal conductivity k e f f of the fluid saturated porous medium and the volumetric heat capacity of the fluid, T 0 = (T 1 + T 2 )/2. The parameters σ, k e f f and χ are defined as [19] where ϕ is the porosity, c f is the specific heat capacity of the fluid, c s is the specific heat capacity of the porous medium, ρ s is the density of the solid phase, k f is the thermal conductivity of the fluid and k s is the thermal conductivity of the solid phase.  The following scaling, where ∆T = T 2 − T 1 , is employed to obtain the dimensionless formulation of the problem: By substituting Equation (4) in Equation (2), the following dimensionless governing equations are obtained where Basic State Let us consider as basic state the steady and fully developed buoyant flow in the vertical direction z and adopt the subscript b to denote the quantities relative to the basic state. Thus, the only non-vanishing component of the filtration velocity vector u b is w b , along the z-direction. According to Equation (5), u b results in a solenoidal field. Thus, the basic stationary solution is characterised by a zero net mass flow rate throughout the layer, namely The solution for the basic state describes the behaviour of a stationary flow which is parallel to the vertical axis, thus yielding a single-cell vertical pattern where the cell has an infinite vertical height.

Linear Stability Analysis
To perform a linear stability analysis, we proceed by perturbing the basic state just defined by applying small-amplitude disturbances, namely where ε 1. By substituting Equation (8) into Equation (5) and by neglecting O(ε 2 ) terms, one obtains By applying the divergence operator to the perturbed momentum balance equation and by applying the linear differential operator (1 + λ 2 ∂/∂t) to the perturbed energy balance equation, one obtains the following pressure-temperature formulation, which is more convenient than Equation (9) since the boundary conditions are relative to P and Θ We assume the pressure and temperature perturbations to have the form of normal modes, where η is the growth or decay rate, ω is the angular frequency and k = (0, k y , k z ) is the wave vector. The parameters (k y , k z , η, ω) are real while ( f , h) are, in general, complex functions. The growth rate η marks the difference between stability (η < 0) and instability (η > 0). The neutrally stable configuration is identified by η = 0. The condition of minimum R among the neutrally stable modes defines the critical values k c and R c . By substituting definitions (11) into Equation (10), one obtains the following eigenvalue problem for neutrally stable modes, where We mention that 0 < γ 1 where, according to Equation (1), the limiting case γ = 1 represents the case of a Newtonian fluid.
Moreover, the rescaled Darcy-Rayleigh number S accounts for the inclination of the wave vector to the vertical z axis. When k y = 0, one has a transverse mode propagating along the z axis and S = R. On the other hand, when k z = 0, one has a longitudinal mode propagating along the y axis and S = 0. Since the latter result can be obtained also when R = 0, horizontally propagating modes are equivalent to transverse modes having a vanishing Darcy-Rayleigh number. Thus, modes having k z < k and a given value of S correspond to transverse modes (k z = k) with the same S and a larger value of R. As a consequence, the transverse modes necessarily are the most unstable.

Numerical Solution and Discussion of the Results
The eigenvalue problem formulated through the dimensionless governing Equation (12) can be solved numerically by means of the shooting method. In the considered approach, the associated initial value problem is solved with the Runge-Kutta method, while Brent's method is adopted for the root searching procedure. The numerical solution of the problem is here found through the built-in functions NDSolve and FindRoot available within the Mathematica 12.0 environment.
In detail, the numerical approach used to solve the eigenvalue problem given by Equation (12) requires the input values (λ 1 , γ) and determines both S(k) and ω(k) as eigenvalue data. The minimum of S(k) yields the critical values (k c , S c , ω c ) that define the onset of the convective instability.
The main requirement for an efficient solution of Equation (12) through the shooting method is an accurate initial guess for the eigenvalue quantities S and ω, for a given prescribed input dataset (λ 1 , γ).

Limiting Case of a Newtonian Fluid
As already emphasised, the limiting case of a Newtonian fluid occurs when γ = 1. In this case, regardless of the value of λ 1 , the neutral stability curve, displayed in Figure 2, shows up a minimum in the plane (k, S) for the critical values which hold for every choice of λ 1 = λ 2 . It is worth noticing that, as it can be inferred from Equation (

Results for Viscoelastic Fluids
The neutral stability curves for the considered data (λ 1 , γ) are reported in Figure 3. Each frame refers to a different value of γ, while the solid lines are drawn for given values of λ 1 . For decreasing values of λ 1 , the curves move upward, as well as for increasing values of γ. It is worth noticing that the vertical axis has the range 0 S 197.081, since this value corresponds to the onset of instability for a Newtonian fluid, namely the maximum possible critical value of S. The threshold values of k, S, ω that have been found numerically are reported in Table 1. The obtained results (k c , S c ) are insensitive to the sign of ω c . Thus, in Table 1, the values of ω c are to be intended as |ω c |. As can be inferred from Table 1, for a given value of γ, an increase in the dimensionless relaxation time λ 1 leads to a decrease in the threshold values k c , S c , ω c (at least, whenever ω c = 0). The minimum (k c , S c ) of the neutral stability curve moves to the left and becomes smaller. Table 1. Critical values of (k, S, ω) versus λ 1 , γ.

Conclusions
The effects of viscoelasticity on the onset of convection in a fluid saturated vertical porous layer have been analysed. An extended version of Darcy's law has been utilised in order to implement the Oldroyd-B rheological model for the seepage flow in the porous medium. The boundaries of the vertical porous layer have been considered as permeable to external fluid reservoirs kept at different temperatures. The boundary temperature difference causes a side-heating mechanism and, hence, a stationary buoyant flow driven by the uniform horizontal temperature gradient. Such conditions define the basic conduction state. The linear stability of this basic state has been studied by employing a three-dimensional normal mode scheme, with a suitable Squire transformation mapping all possible modes to equivalent transverse rolls. Finally, the stability eigenvalue problem has been formulated and solved numerically by employing the shooting method.
The main highlights of our study are the following: • The dimensionless governing parameters identifying the viscoelastic behaviour are the relaxation parameter, λ 1 , and the ratio between the retardation time and the relaxation time, γ. The physically significant domain is one where 0 < γ 1. When γ → 1, the Newtonian fluid behaviour is recovered: the critical value of the modified Darcy-Rayleigh number is 197.081 and the corresponding critical value of the wavenumber is 1.05950; • The neutral stability threshold to the convective instability obtained for the Newtonian case always exists, whatever is the choice of λ 1 and γ. However, in most cases, the Newtonian branch of neutral stability is not the lowest one. In these cases, the neutral stability condition is characterised by travelling modes, i.e., modes with a non-zero angular frequency. Thus, the effect of viscoelasticity is generally destabilising with respect to the Newtonian fluid case; • There exist input values of (λ 1 , γ), with γ < 1, such that the critical conditions for the onset of the convective instability coincide with those for a Newtonian fluid or, equivalently, the Newtonian branch of neutral stability is the lowest one.
There are several possible directions where the study presented in this paper can be extended. In Authors' opinion, one of the most interesting developments of this study is the introduction of a model for partially permeable boundary conditions. With such a model one can investigate the gradual transition from permeable to perfectly impermeable boundaries, and its effects on the onset conditions for the instability.