The energy method analysis of the Darcy–Bénard problem with viscous dissipation

A nonlinear analysis of the effect of viscous dissipation on the Rayleigh–Bénard instability in a fluid saturated porous layer is performed. The saturated medium is modelled through Darcy’s law, with the layer bounded by two parallel impermeable walls kept at different uniform temperatures, so that heating from below is supplied. While it is well known that viscous dissipation does not influence the linear threshold to instability, a rigorous nonlinear analysis of the instability when viscous dissipation is taken into account is still lacking. This paper aims to fill this gap. The energy method is employed to prove the nonlinear conditional stability of the basic conduction state. In other words, it is shown that a finite initial perturbation exponentially decays in time provided that its initial amplitude is smaller than a given finite value.


Introduction
The effects of viscous dissipation on fluid flow have been widely investigated from different perspectives. It is well known that viscous dissipation can be a significant heat source when a forced flow is imposed. Its strength can be sufficient to yield important changes for the temperature field and for the heat transfer rate. In the case of natural convection, when the flow is induced only by the buoyancy force, the internal heating caused by viscous dissipation is less intense and, in most cases, it can be neglected [1]. In this framework, an important question arises on whether the viscous dissipation effect can cause, or influence, the onset of convective instability when a fluid layer, or a fluid-saturated porous layer is subject to a downward-oriented temperature gradient. In this regard, an important difference occurs between basic states where the fluid is at rest and situations where a basic steady flow exists [2]. A basic state with a fluid at rest and heated from below is typical of the Rayleigh-Bénard system [3], or of its saturated porous medium version, well known as the Darcy-Bénard system [4][5][6][7].
The effect of viscous dissipation cannot modify the basic downward temperature gradient which triggers the Darcy-Bénard instability as the fluid is at rest in the basic state. However, when a disturbance acts on the basic state inducing a transient flow, then viscous dissipation may play a role. This is not the case if Communicated by Andreas Öchsner.
A. Barletta  a linearised scheme for the analysis of the instability is carried out. In fact, the viscous dissipation in the perturbed flow would be quadratic with respect to the perturbation amplitude and, hence, negligible in this scheme. The scenario may markedly change when a fully nonlinear approach to the dissipation Darcy-Bénard system is pursued. This approach can be based on the energy method [3,4], where the growth or decay of a disturbance acting on the basic state is established by employing norm inequalities based on functional analysis results relative to disturbance fields having compact support. The compact support is usually due to the finite periodicity cell characterising the disturbances. In fact, this analysis scheme presumes a periodic cellular pattern tiling the horizontal plane.
The aim of this paper is the application of the energy method to investigate the nonlinear instability threshold when the Darcy-Bénard system is modelled by including the effect of viscous dissipation. The result proved through this study is that there is a nonlinear conditional stability in the subcritical domain. This means that finite disturbances decay in time whenever the energy functional expressing their amplitude is smaller than a threshold value. We mention a previous study of the effect of viscous dissipation in the Darcy-Bénard system instability [8]. These authors performed a two-dimensional numerical study where the existence of subcritical instability caused by the viscous dissipation effect was ruled out. Another recent numerical study regarding the nonlinear instability of the Darcy-Bénard system including the effect of viscous dissipation was carried out in [9]. This study is focussed on the weakly nonlinear approach to the investigation of the convection patterns and, in particular, to the formation of hexagonal cells. We believe that the present paper provides both a validation and a consolidation of the conclusions drawn in [8], as well as in [9]. In fact, the analysis carried out here is not based on numerical results relative to the special two-dimensional case, but on an entirely analytical procedure relative to general three-dimensional disturbances.

Mathematical model
Let us consider a horizontal porous layer with infinite horizontal width and vertical thickness H , bounded by two impermeable isothermal planes at temperatures T 0 + T (lower boundary), T > 0, and T 0 (upper boundary). The z axis is vertical with the origin on the lower boundary plane, while the x and y axes are horizontal. Here, the stars denote the dimensional time, coordinates and fields. We write the dimensionless governing equations for a fluid saturated porous medium, according to Darcy's law and to the Oberbeck-Boussinesq approximation, and including the viscous dissipation contribution in the local energy balance equation, where u i is the velocity, T is the temperature, p is the pressure, (x i ,ŷ i ,ẑ i ) are the unit vectors along the (x, y, z) axes, R is the Darcy-Rayleigh number, and B is the Brinkman number. When B 1, the viscous dissipation contribution, Bu i u i , becomes negligible. The dimensionless quantities are defined through the scaling where ρ is the reference fluid density, α is the thermal diffusivity, μ is the dynamic viscosity, β is the coefficient of thermal expansion, g is the modulus of the gravitational acceleration, K is the permeability of the porous medium, c is the specific heat of the fluid, and σ is the ratio between the average volumetric heat capacity of the saturated porous medium, (ρc) m , and the volumetric heat capacity of the fluid, ρc. Hereafter, the components of u i are also denoted as (u, v, w). The boundary conditions are In (1) and in the following, we use Einstein's notation for the sums over repeated indices.

The basic solution and the disturbances
A stationary solution of (1) subject to the boundary conditions (3) is We now express the disturbances relative to the basic solution (4) as On substituting (5) into (1) and (4), we obtain the governing equations for the disturbances where W is the z-component of U i . Equation (6) are the starting point for a nonlinear stability analysis of the basic solution (4). If we neglect the nonlinear terms, U i θ ,i and BU i U i , in (6) 3 we obtain the linearised eigenvalue problem The solution of (7) yields the linear stability bound R = R , such that the instability is defined by the values of R greater than R . We note that the viscous dissipation contribution is absent in (7), and hence it has no effect in the linear stability analysis. Therefore, we conclude that the well-known result of the linear stability analysis of the Darcy-Bénard problem [4-6,10-13], namely R = 4π 2 , is not altered if the effect of viscous dissipation is taken into account. Due to the possible existence of the subcritical instability, R < R generally does not ensure the stability of the basic solution. This is the motivation of the nonlinear stability analysis, based on (6). The study of (6) through the energy method allows one to ascertain if there exists a nonlinear stability bound R E ≤ R , implying the existence of a subcritical instability [3,4]. The condition R < R E ensures the stability of the basic solution.
We know that, if the effect of viscous dissipation is neglected, R E = R and hence no subcritical instability occurs [4].
Because of the presence of the quadratic term BU i U i , here we investigate if the inequality R < R will imply conditional nonlinear stability. If this is the case, then R E = R .

The energy method
We assume that the convective instability implies the onset of periodic patterns tiling the horizontal plane (x, y). Let V be the spatial periodicity cell, V = × [0, 1], where is the general tile in the plane (x, y). At the vertical boundaries of V we assume the symmetry conditions U ini = 0 and θ ,ini = 0, wheren i is the unit normal to these boundaries. We will use henceforth the notations where F ∈ L 2 (V ) and S i ∈ L 2 (V ) are any square-integrable scalar and vector functions over V .
We will also assume that U i ∈ L 2 (V ) (U i is square-integrable over V ), and that θ ∈ H 1 (V ) (both θ and θ ,i are square-integrable over V ). Both U i and θ satisfy the symmetry conditions on the boundary ∂ , as well as the conditions (7) 4 at z = 0 and z = 1. Use will be made of the Poincaré inequality [14][15][16] and of the arithmetic-geometric mean inequality, We now multiply (6) 2 by U i and integrate over V . Thus, we write After an integration by parts and use of (6) 1 , (6) 4 , as well as of the symmetry conditions U ini = 0 at the vertical boundaries of V , we obtain We now multiply (6) 3 by θ and integrate over V . We write After an integration by parts and use of (6) 1 , (6) 4 , as well as of the symmetry conditions U ini = 0 and θ ,ini = 0 at the vertical boundaries of V , we conclude that where ∂ V is the boundary of V . As a consequence, (13) can be rewritten as It is now convenient to define the scaling so that (12) and (16) now read We sum (17) and λ times (18), where λ is a real positive coupling parameter [4]. Thus, we obtain where The nonlinear stability bound, R E , can be determined so that where the maximum is evaluated over all the possible U i ∈ L 2 (V ), ∈ H 1 (V ) fulfilling the constraint U i,i = 0, the boundary conditions = 0 = W at z = 0, 1, and the symmetry conditions U ini = 0 and ,ini = 0 at the vertical boundaries of V . In order to justify (23), we mention that the optimal definition of the stability bound is obtained when R E is at its maximum. In fact, the condition (23) means that R E is the maximum positive real number such that I /D < 1/R E .

The Euler-Lagrange equations for the nonlinear stability analysis
In order to determine the nonlinear stability bound R E = R 2 E , we must obtain the maximum of I /D, as stated in (23). This has been done by Straughan [4]. For completeness, we here report his calculations. We evaluate the variation of I /D, The maximum is determined by setting δ{I /D} = 0. By taking into account (23), we obtain It is now convenient to rescale by replacing it with / √ λ. Then, (21) and (22) yield Therefore, where (h i , η) are arbitrary fields that satisfy the same boundary conditions as (U i , ) on ∂ V . The constraint U i,i = 0 can be taken into account by incorporating an additional term in the definition of I , so that Here is an arbitrary scalar field playing the role of a Lagrange multiplier. In (27) 1 , an integration by parts has been done, and use has been made of the boundary condition U ini = 0 on ∂ V . From (26) and (27) we obtain In (27), an integration by parts has been done, and use has been made of the condition ,ini = 0 on the vertical boundaries of V , and of the condition η = 0 on z = 0, 1. We now substitute (??) and (27) Thus, we may write the system of differential equations It is easy to see that the optimal value of λ is 1. With this value, provided that one identifies θ with /R E , the stationary linear perturbation equation (7) coincide with Eq. (30) for R E = R , where R = 4π 2 is the critical Rayleigh number for linear instability.

Conditional nonlinear stability in the energy norm for R < R E
From (19) and (23), we obtain the inequality where We now use first the Schwarz inequality and, then, the Hölder inequality [3] to write where · 4 and · 6 are the usual norms in L 4 (V ) and in L 6 (V ), and C 1 = |V | 1/6 (|V | is the measure of V ). Then, we apply the Sobolev inequality and the generalised Poincaré inequality (see [17] or Lemma 2.1 in [18]) where C 2 is a positive constant and C 3 is a positive constant depending on V . Our next step is to prove the inequality To obtain (34), we take the Laplacian of (6) 2 : Multiplying by −U i and integrating over V , we have We note that U = (U 1 , U 2 , W ) is a divergence free vector and that the boundary conditions on z = 0, z = 1, as well as the symmetry conditions at the vertical boundaries, imply that the normal component of the velocity vanishes on the whole boundary of V . Then, the first term on the right hand side of (36) vanishes. Moreover, the second term integrates to obtain R θ , j W , j . It is more difficult to employ the integration by parts for the term − U i, j j U i as the boundary conditions yield W = 0 at z = 0 and z = 1, while no boundary conditions on U 1 and U 2 are imposed there. Integrating the left hand side of the previous identity, because of the boundary conditions and the symmetry conditions, we have From (6_1,2) we have U 1 = −P ,1 , U 2 = −P ,2 , W = −P ,3 + Rθ. From this we deduce − U 1 U 1,33 − U 2 U 2,33 = − P ,1 P ,133 + P ,2 P ,233 = P ,11 P ,33 + P ,22 P ,33 We now compute By employing (37)-(39) we finally infer Then, by using the Schwarz inequality and (16), we obtain (34). From (32)-(34), we have where C 4 = C 1 C 2 3 R 2 . From the last inequality, by recalling that the optimal value of λ is 1, one obtains where C 5 = √ 2C 4 BR. By employing (31) and (42), one may write the energy differential inequality Now, assuming that by a recursive argument (see [19]), we prove that dE/dt < 0, ∀t ≥ 0, and that We denote with α the quantity −β + C 5 E 1/2 (0), and we observe that α < 0 and dE/dt α D. Since ∇ 2 D, on using the Poincaré inequality (9), namely π 2 2 ∇ 2 , one has From (43) and (46) we infer dE dt Finally, an integration yields If R < R E , so that β > 0, and if (44) holds, then α < 0 and the exponential time decay of E(t) means the time decay of 2 . This implies the nonlinear conditional stability, at least for the temperature field. The conditional stability of the velocity field holds as well. In fact, we can start from (17) and use the arithmetic-geometric mean inequality (10), with γ = 1/R, so that we obtain As a consequence, the time decay of 2 yields the time decay of U 2 .

Conclusions
A nonlinear analysis of the onset of the Rayleigh-Bénard instability in a porous layer has been carried out. The porous medium is described through the Oberbeck-Boussinesq approximation, by assuming a non-negligible effect of viscous dissipation. Momentum transfer is modelled by employing Darcy's law. The layer is assumed to be bounded by impermeable planes where uniform unequal temperatures cause heating of the layer from below. The energy method is adopted as a tool for the investigation of the growth or decay of perturbations superposed to the basic state. The latter is one of zero velocity and a uniform downward-oriented temperature gradient. A suitable energy functional is defined whose time evolution captures the stable or unstable nature of the perturbations. Euler-Lagrange equations are employed to build an inequality for the time-derivative of the energy functional. Further information has been gained by employing a chain of inequalities based on functional analysis. The final conclusion is that, in the subcritical regime, there is a nonlinear conditional stability of the basic state.
The nonlinear stability condition obtained here is to be expected because of the cubic term N = BR U i U i that could increase when U i and are sufficiently large (the same happens for the cubic terms in other Bénard problems [19]). Nevertheless, conditions (44) and (45) imply conditional stability.
There are opportunities for further developments of the study performed in this paper. In fact, a significant feature not investigated in our study is the presence of a basic horizontal throughflow in the porous layer. The throughflow triggers the effect of viscous dissipation in such a way as to alter the basic state. These conditions have been envisaged by [20,21] and reconsidered, within a numerical nonlinear scheme, by [8]. The energy method can add new important insights into the Rayleigh-Bénard instability with throughflow, where the effect of viscous dissipation may affect significantly the critical value of the Darcy-Rayleigh number at the onset of the linear instability.