Natural convection in a fluid saturating an anisotropic porous medium in LTNE: effect of depth-dependent viscosity

Thermal convection in a fluid saturating an anisotropic porous medium in local thermal nonequilibrium (LTNE) is investigated, with specific attention to the effect of variable viscosity on the onset of convection. Many fluids show a remarkable dependence of viscosity on temperature that cannot be neglected. For this reason, we take into account a fluid whose viscosity decreases exponentially with depth, according to Straughan (Acta Mech. 61:59–72, 1986), Torrance and Turcotte (J. Fluid Mech. 47(1):113–125, 1971). The novelty of this paper is to highlight how variable viscosity coupled with the LTNE assumption affects the onset of convection. A numerical procedure shows the destabilising effect of depth-dependent viscosity. Moreover, it comes out that the LTNE hypothesis makes the influence of viscosity more intense. Linear instability analysis of the conduction solution is carried out by means of the Chebyshev-tau method coupled to the QZ algorithm, which provides the critical Rayleigh number for the onset of convection in a straightforward way. The energy method is employed in order to study the nonlinear stability. The optimal result of coincidence between the linear instability threshold and the global nonlinear stability threshold is obtained. The influence of anisotropic permeability and conductivity, weighted conductivity ratio, and interaction coefficient on the onset of convection is highlighted.

rise in temperature [8], while viscosity of water changes more than 1 order of magnitude between 25 and 350 • C [9]. Over the same temperature interval, water thermal conductivity exhibits only a 1% change [10]. Actually, as stated in [11], for most Newtonian liquids thermal conductivity is essentially constant over the temperature range in which viscosity shows remarkable variations. For this reason, in this paper, we take into account only viscosity variations.
Over the years, thermal convection in variable viscosity fluids has represented a topic of great interest for many researchers. Early studies on thermal convection in clear fluids with temperature-dependent viscosity are those of [2,12]. Since then, such a problem has been extensively investigated [4,[13][14][15][16][17][18]. In [2], the influence of strongly temperature-and/or depth-dependent viscosity has been studied, with particular attention to convection in Earth's mantle. In [4], the conditional nonlinear stability is investigated for a clear fluid whose viscosity decreases linearly with temperature. In [15], a global nonlinear stability result is obtained when viscosity shows a maximum. Nonlinear stability results have been obtained also in [16], where viscosity shows exponential dependence on temperature, and in [17] where viscosity is a convex function of temperature.
The problem of thermal convection in variable viscosity fluids saturating a porous medium has been widely investigated, as well [1,9]. Nonlinear stability results have been determined in [19] and [20]. In [5], linear stability analysis of convection in water-saturated porous medium is performed. The above-cited papers deal with the problem of thermal convection in porous media under the hypothesis of local thermal equilibrium (LTE) between solid and fluid phase, i.e. by employing a single temperature model. Nevertheless, once the fluid velocity is sufficiently high, or as long as the solid thermal conductivity is much different from the fluid one, the assumption of LTE is no longer suitable to describe the physical phenomenon [21]. Applications are numerous and can be found in processes involving quick heat transfer, metal foams, and in everyday technology such as microwave heating [8], tube refrigerators and heat exchangers (see [22] for more details). As a result, the hypothesis of Local Thermal Non Equilibrium (LTNE) is employed to obtain more accurate results. This assumption involves two different temperatures (one for the solid, one for the fluid) so as to take into account heat exchanges between the phases. The problem of thermal convection in porous media in LTNE has been widely studied over the last years, starting from [23,24], in which both a linear and nonlinear analysis is carried out. Many papers analysed the effect of LTNE on the onset of convection, coupled to rotation and/or anisotropy [25][26][27][28][29], magnetic field [30], vertical throughflow [31,32], and also in non-horizontal porous layer [33][34][35]. Interesting results are obtained once the assumption of LTNE is combined with the Cattaneo's law for heat conduction in the solid matrix, since in such a case convection by means of oscillatory motions can occur [36][37][38][39].
At the best of our knowledge, natural convection in depth-dependent viscosity fluids in porous media in LTNE has not received a proper attention so far, that is why we believe this paper may be of great interest. The LTNE model has been employed in [8], where the onset of Darcy-Benard convection is investigated for a temperature-dependent viscosity fluid with quadratic density constitutive equation.
The plan of the paper is the following. Section 2 is devoted to the formulation of the mathematical model describing the motion of a depth-dependent viscosity fluid in presence of an anisotropic porous medium, in LTNE regime. In this Section, the basic steady solution and the dimensionless model governing the evolution of perturbation fields are determined. In Sect. 3, the instability analysis is performed in order to determine the generalised eigenvalue problem which can be solved by means of the Chebyshev-tau method coupled to the QZ algorithm. In Sect. 4, the global nonlinear stability analysis of the conduction solution is carried out by employing the energy method. Finally, in Sect. 5, numerical results are reported, with particular attention devoted to the influence of variable viscosity, mechanical and thermal anisotropy, weighted conductivity ratio, and interaction coefficient on the onset of instability.

Mathematical model
Let us take into account a horizontal porous layer, whose depth is d, saturated by an incompressible fluid F initially at rest. We choose a Cartesian frame of reference (x, y, z), being the z-axis vertically upward. A uniform temperature gradient is imposed and maintained across the medium. Let T L be the temperature on the lower plane z = 0 delimiting the layer, and let T U be the temperature on the upper plane z = d. We assume the layer to be heated from below, i.e. T L > T U . Moreover, the hypothesis of local thermal non-equilibrium holds for the porous medium at stake, i.e. the heat exchange between solid matrix and fluid is allowed. Based on this assumption, we need to define two different temperatures, one for the fluid phase, another for the solid one.
Let us denote by T f the temperature of the fluid phase and by T s the temperature of the solid phase. Then, by virtue of previous considerations, with T L > T U . The porous medium exhibits some anisotropic properties. Specifically, permeability and solid thermal conductivity are horizontally isotropic. Let K be the permeability tensor and let D s be the solid thermal conductivity one. By assuming that both tensors share the same principal axis (x, y, z), they can be written in diagonal form, i.e., In addition, we allow viscosity to be strongly dependent on temperature, at least at the steady state where temperature linearly increases with depth [10]. As a consequence, we assume the following behaviour of viscosity with respect to z [1,2]: being μ 0 and c positive constants. Under the previous assumptions, by employing the Oberbeck-Boussinesq approximation, the Vadasz-Darcy model is [8] ⎧ where v, p , T f , and T s are (seepage) velocity, pressure, fluid phase temperature, and solid phase temperature, respectively; while μ(z) is the fluid viscosity according to (3); κ f , ρ f , ρ s , c f , c s , g, α, c a , ε, h are fluid thermal conductivity, fluid and solid density, specific heat of fluid and solid phase, gravity acceleration, thermal expansion coefficient, acceleration coefficient, porosity, angular velocity and interaction coefficient, respectively.
The following boundary conditions are coupled to (4): where T L > T U and n is the unit outward normal to planes z = 0, d.
The basic steady solution of (4) is: is the adverse temperature gradient. In order to study the stability of m 0 , let us introduce the perturbation fields (u, θ, φ, π) to velocity field, fluid temperature field, solid temperature field, and pressure field, respectively.
After introducing perturbations to initial data, the following solution of (4) arises: Once the following dimensionless quantities are defined, omitting all asterisks, we obtain the following dimensionless model governing the evolution of perturbation fields: 3) (11.4, 5) are diffusivity ratio, weighted conductivity ratio, inter-phase heat transfer coefficient, Rayleigh number, and Vadasz number, respectively. To (10.1), we add the initial conditions: where ∇ · u 0 = 0, while the boundary conditions are Let be the periodicity cell and let us define (·, ·) and · inner product and norm on L 2 (V ), respectively. Perturbations are assumed to be periodic in x and y directions with periods 2π a x and 2π a y , respectively. Moreover, ∀t ∈ R + , and g can be expanded in a Fourier series uniformly convergent in V .

Linear instability analysis
In order to perform the linear instability analysis of m 0 , let us consider the linear version of (10 that can be written as being F = (u, v, w, π, θ, φ), N with all zero entries and such that diag It is straightforward to notice that the spatial operator M related to (15) is symmetric with respect to the L 2 -scalar product. As a consequence, its spectrum involves only real numbers, and the strong version of the principle of exchange of stabilities holds. Moreover, system (15) is autonomous. Hence, we can look for solutions such that Substituting (18) into (15), we get In order to determine the critical threshold of the Rayleigh number for the onset of instability, we focus on the marginal state in (19). Accounting for the principle of exchange of stabilities, we set σ = 0 in (19), so that By applying the double curl to (20.1) and retaining only the third component, from (20.1), one gets: where the prime denotes ordinary derivative with respect to z. Because of periodicity of the perturbations, unknown fields in (21) can be written as: with a 2 = a 2 x +a 2 y . Now, let us substitute (22) in (21) and retain only the n-th component. Hence, by denoting D ≡ d dz and which will be solved subject to the boundary conditions System (24) represents an eigenvalue problem of ordinary differential equations, where R is the eigenvalue. Starting from (24), one can determine the function λ(a 2 ) which describes how R depends on a 2 . As a consequence, the critical Rayleigh number for the onset of stationary instability is given by In order to solve the eigenvalue problem (24), we employ the Chebyshev-tau method coupled with the QZ algorithm. This numerical procedure takes advantage of the orthogonality of Chebyshev polynomials with respect to the scalar product, and many of their properties such as a recursive formula for multiplication between two polynomials, i.e.
2T r T q = T r +q + T |r −q| .
We will not go into further details of this numerical method, but we are going to give the main ideas to carry out this procedure. For more details, interested readers can refer to [10,[40][41][42][43].
Let T k , k ∈ N, be the k-th Chebyshev polynomial. The solution of (24) can be expanded as finite series of Chebyshev polynomials, i.e.
We substitute (29) in (24), and then we take the inner product with T k , k = 0, . . . , N , on the weighted Chebyshev polynomial space. In such a way, we manage to write system (24) as a generalised eigenvalue problem and being D 2 , D the Chebyshev representation of d 2 dz 2 and d dz , respectively, I the identity matrix, and F * the matrix representation of function f −1 (z).
For coding purpose, we work with k = 1, . . . , N . Then, the above matrices are 3N × 3N , and F * is N × N . When the exponential function f −1 (z) can be well approximated by a second-order polynomial, matrix F * is computed by means of the following recursive formulas: The eigenvalue problem (30)-(32) is then solved by the QZ algorithm which provides eigenvalues with no trouble. In Sect. 5, we show numerical results.

Nonlinear stability
This Section is devoted to develop a nonlinear stability analysis for the conduction solution m 0 , in order to determine a sufficient condition for its stability. We perform the scalar multiplication of (10.1) by u, (10.3) by θ and (10.4) by φ, and then we integrate over the periodicity cell V . The sum of the resulting equations can be written as where Lyapunov Function, (35.1) Starting with (34), by defining 1 where H = (u, θ, φ) : w = θ = φ = 0 on z = 0, 1; periodic in x and ydirections, with period 2π a x , 2π a y respectively; ∇ · u = 0; D < ∞ , we obtain dE dt By employing the Poincaré inequality in (35.3), we get being ξ * = min{1, ξ −1 } and ζ * = min{1, ζ }. Hence, (38), (39) yield the exponential decay of the energy function E(t). In fact, by virtue of (39), inequality (38) becomes where k = min 2π 2 , By virtue of (41), as long as R < R E , perturbation fields on seepage velocity, fluid and solid temperature decay exponentially in time. Thus, we have proved that R < R E is a sufficient condition for the global nonlinear and exponential stability of the conduction solution m 0 . Now, let us proceed to determine the critical threshold R E . With this aim, we write the Euler-Lagrange equations by solving the variational problem (36) where ω is a Lagrange multiplier arising from the incompressibility of u. System (42) coincides with (20.1). As a consequence, we manage to obtain the coincidence between the global nonlinear stability threshold R E and the linear instability threshold R L . As a result, condition R < R E is not only sufficient, but also necessary for the stability of m 0 .

Numerical results
In the present Section, we report results obtained once the generalised eigenvalue problem (30) and the subsequent minimum problem (26) have been solved. Given the coincidence between linear and nonlinear threshold, we denote by R 2 c the critical Rayleigh number beyond which thermal convection occurs. The aim of this Section is to highlight the influence of variable viscosity, thermal and mechanical anisotropy, weighted conductivity ratio, and interaction coefficient on the onset of convection. Let us remark that, since the interaction heat transfer coefficient H cannot be easily measured, we set a range in which H can vary. According to the choice done in [29,44], H ∈ (10 −2 , 10 6 ). As a consequence, the critical Rayleigh number is plotted as function of H in every picture throughout the Section.   Table 3 in [1] and the critical wavenumber a 2 c and the rescaled critical threshold It may be observed that the effect of variable viscosity on the critical Rayleigh number is substantial. As shown in Fig. 1, as c increases, R 2 c decreases for every choice of the interaction coefficient H , namely the onset of convection gets easier. This result is physically reasonable and in agreement with what stated in [9], i.e. fluids whose viscosity near the hot boundary is much lower than everywhere else in the layer are more likely to become unstable sooner.
Let us remark that, as H → ∞, fluid and solid phases can be treated as a single phase since they exchange heat so rapidly that they have nearly identical temperature. As a consequence, we recover the regime of LTE. Whereas, when H → 0, fluid instability is not affected by the properties of the solid matrix. As pointed out in [23], the respective mathematical problems are identical, except for a rescaling of the Rayleigh number. Table  1 shows the coincidence between results delivered in [1] under the assumption of LTE and results obtained in  Table 1, where c = 0, we recovered the classical Rayleigh number 4π 2 for the Darcy-Benard problem. The decreasing trend of R 2 c with respect to increasing c, i.e. a higher gradient of viscosity across the layer, depends strongly on H . So, we analyse the behaviour of R 2 c both for high values of H (for which the LTE regime holds [23]) and for H = 100 (for which the LTNE regime holds). In Table 2, one can notice that the decreasing trend of R 2 c is less remarkable for H → ∞ rather than for H = 100, even though the critical threshold for the LTE case is always lower than the one in LTNE. Such a result suggests that the effect of variable viscosity on the onset of convection is less intense under the assumption of LTE rather than in the case of LTNE. A similar result has been pointed out in [8].
When c = 0, function f (z) = 1, i.e. viscosity is constantly equal to μ 0 , the value assumed in the middle of the layer. In case of isotropic porous medium, by assuming c = 0 we would expect to recover the same results found in [45]. In Fig. 2 we compare results delivered in the present paper by the Chebyshev-tau method with the analytical expression of the critical Rayleigh number found in [45]. The two neutral curves are perfectly overlapped. In fact, one can compute the absolute error, which is 1.85 × 10 −10 at most.
In Fig. 3, the influence of mechanical anisotropy on the critical Rayleigh number is shown. The destabilising effect of increasing permeability on the conduction solution is clear. The Rayleigh number decreases with ξ for any H , in agreement with findings in [29,44,46,47]. This result is physically admissible. Based on definition (2.1), increasing ξ is due to an increase in the horizontal permeability K H , which means that the fluid is allowed to move easier and easier in the horizontal direction. Thus, resistance to fluid motion reduces, and the advective transport enhances [48]. On the other hand, solid thermal conductivity has a stabilising effect on m 0 , delaying the onset of convection. In Fig. 4, such a behaviour is evident, specifically for large values of H . If κ s H grows, the onset of convection is delayed because the solid matrix easily absorbs heat from the fluid. Actually, when H goes to zero, solid thermal conductivity affects weakly the onset of convection. This result is in agreement with what is found in [29], and it is physically reasonable since small values of H imply that the heat exchange between fluid and solid phases occurs so slowly that the solid matrix does not manage to absorb enough heat to cool down the fluid.
The onset of convection is encouraged by increasing values of the weighted conductivity ratio γ , as shown in Fig. 5. Increasing values of κ f lead to an increasing γ , which imply a reduction of R 2 c , i.e. a destabilising effect. This behaviour is expected from a mathematical point of view, since looking at the Rayleigh number definition in (11.4) we can notice that R 2 is inversely proportional to κ f .

Conclusions
A linear and nonlinear stability analysis of the conduction solution of the Vadasz-Darcy model describing the motion of a depth-dependent viscosity fluid in an anisotropic porous medium in LTNE is performed. Linear instability analysis leads to a generalised eigenvalue problem which can be solved by means of the Chebyshevtau method coupled with the QZ algorithm. Nonlinear stability analysis is carried out by employing the energy method. The optimal result of coincidence between the linear instability threshold and the global nonlinear stability threshold is obtained.
Viscosity decreasing exponentially with depth has a destabilising effect on the onset of convection, leading to lower values of the critical Rayleigh number. In fact, viscosity drag reduces with depth, namely the fluid meets fewer obstacles to its motion. Moreover, comparing results with those ones found under LTE assumption in [1], we can remark that the presence of LTNE makes the influence of variable viscosity more intense. Applications can be found in oil cooling systems in hydraulic units that involve heat exchangers. Within the heat exchanger, maintaining the oil cold is important in order to preserve its characteristics and proper operating conditions. Findings in the paper show that the effect on the onset of convection of replacing a kind of oil with another with a different law of viscosity is larger under the LTNE hypothesis than under the LTE one. A practical example is given by two oils, i.e. SAE 10W-60 and SAE 0W-30, whose values of kinematic viscosity ν are shown in the Furthermore, we proved the destabilising effect of both mechanical anisotropy and the weighted conductivity ratio, other than the stabilising effect of thermal anisotropy.