## Abstract

This paper investigates the instability thresholds and global nonlinear stability bounds for thermal convection in a fluid overlying a highly porous material. A two-layer approach is adopted, where the Darcy–Brinkman equation is employed to describe the fluid flow in the porous medium. An excellent agreement is found between the linear instability and unconditional nonlinear stability thresholds, demonstrating that the linear theory accurately emulates the physics of the onset of convection.

## 1. Introduction

Flow of a fluid overlying a porous medium saturated with the same liquid is encountered in a wide range of industrial and geophysical applications, and has consequently received much attention in the literature (see Nield 1977, 1991, 1998; Chen & Chen 1988, 1992; Chen 1990; Lu & Chen 1997; Das *et al.* 2002; Discacciati *et al.* 2002; Layton *et al.* 2003; Miglio *et al.* 2003; Chang 2004, 2005, 2006; Riviere 2005; Chang *et al.* 2006; Das & Lewis 2007; Hirata *et al.* 2007; Hoppe *et al.* 2007; Mu & Xu 2007; Hill & Straughan 2008; Straughan 2008).

The modelling of a fluid/porous system can be addressed using a two-layer approach, where the governing equations in the separate fluid and porous regions are coupled by appropriate interfacial boundary conditions. Modelling the fluid flow in the porous medium using Darcy's equation is extensive in previous studies (Nield 1977; Chen & Chen 1988; Carr & Straughan 2003; Carr 2004). However, when studying a highly porous material, the use of the higher order Darcy–Brinkman equation becomes more appropriate (Nield & Bejan 2006). In this paper, we allow our highly porous material to be a foametal, which is used extensively in industrial applications such as heat exchangers, chemical reactors and fluid filters. This material has well-documented values for its permeability and porosity (Straughan 2002; Goyeau *et al.* 2003), which we will use in our analyses.

An alternative modelling stance is to assume a one-domain approach, which couples the porous media flow with the viscous fluid flow in a single form of equations with highly discontinuous coefficients. The use of the Brinkman equation to solve the Darcy–Stokes problem in this one-domain approach is discussed in Xie *et al.* (2008). A comparison between the one- and two-domain approaches for the Navier–Stokes/Darcy–Brinkman problem may be found in Hirata *et al.* (2007).

This paper is concerned with constructing a stability analysis for thermal convection in a fluid overlying a highly porous material. There are at least two key techniques in a stability analysis, namely the energy method and the method of linearized instability. Owing to its inherent nature, a linear stability analysis only provides boundaries for instability due to the presence of nonlinear terms (Straughan 2004). The use of the energy method to construct stability thresholds is therefore crucial to assess whether the linear theory accurately encapsulates the physics of the onset and the behaviour of convection (see Payne *et al.* 1999; Flavin & Rionero 2003; Straughan 2006). Owing to the difficulties associated with handling the interfacial boundary conditions (Payne & Straughan 1998), nonlinear energy stability analyses of fluid/porous systems are not widespread in the current literature. The main difficulty arises owing to the nonlinear ** v**.∇

**term in the Navier–Stokes equations, which does not have a counterpart in the Darcy equations in a porous media. This is discussed at length in ch. 6 of Straughan (2008). In this paper, both linear instability and unconditional nonlinear stability results are developed, making this work particularly timely.**

*v*All numerical results were derived using the Chebyshev tau-QZ method (Dongarra *et al.* 1996), which is a spectral method coupled with the QZ algorithm, and were checked by varying the number of polynomials to verify convergence. Standard indicial notation is employed throughout.

## 2. Formation of the problem

Let us consider a fluid occupying the three-dimensional layer , overlying a homogeneous porous medium occupying the layer . The interface between the saturated porous medium and the fluid layer is at *z*=0 (figure 1).

Gravity acts in the negative *z*-direction, and we assume that the density *ρ* has a linear temperature dependence of the formwhere *T* is the temperature; *ρ*_{0} and *T*_{0} are the reference density and temperature values, respectively; and *α* is the coefficient of thermal expansion. The governing equations in the fluid are given by the Stokes flow equations(2.1)where (2.1)_{2} and (2.1)_{3} are the incompressibility condition and balance of energy, respectively. In these equations, *ν*_{i}, *p*, *t* and *μ* are the velocity, pressure, time and dynamic viscosity, respectively, and *κ*_{f}, *g* and *c*_{p} are the thermal conductivity, acceleration due to gravity and specific heat at a constant pressure, respectively. It is important to note that, although we adopt the Stokes equations for fluid flow, system (2.1) is still nonlinear due to equation (2.1)_{3}. Such an approach is referred to as the nonlinear Stokes problem (cf. Duka *et al.* 2007).

Standard indicial notation is employed throughout, with the symbol *Δ* representing the Laplace operator and ** k**=(0, 0, 1). The subscripts (or superscripts) f and m denote the fluid and porous layers, respectively.

Assuming that the porous medium is of high porosity *ϕ*>0.75, the governing equations are those of the Darcy–Brinkman model(2.2)In these equations, the variables , *p*^{m}, *T*^{m} and *K* are the velocity, pressure, temperature and permeability, respectively, in the porous medium. Adopting the approach of Whitaker (1986), the effective viscosity *μ*_{e} is defined as *μ*_{e}=*μ*/*ϕ*. The quantities and are defined in terms of the fluid and solid components of the porous medium, such that where or .

Assuming that the upper and lower planes are fixed, we have the boundary conditions(2.3)where *T*_{L}>*T*_{U}. At the interface *z*=0, the continuity of the velocity, temperature and heat flux yields(2.4)The quantity is the velocity of the fluid at a point ** x** averaged over a representative elementary volume containing

**, i.e. a small volume, but big enough to contain both fluid and solid components at the microscopic level. Thus, is really**

*x**ϕV*

_{i}, where

*ϕ*is the porosity and

*V*

_{i}is the velocity at

**averaged over that part of the representative elementary volume containing only the fluid (cf. Straughan 2008, p. 14). Thus, the equation of velocity continuity in (2.4) is really an equation of mass flux continuity.**

*x*The remaining boundary conditions at *z*=0 are the continuity of normal stresses,(2.5)and tangential stresses,(2.6)for *γ*=1, 2.

Using the boundary conditions (2.3) and (2.4), the governing equations (2.1) and (2.2) admit a steady-state solution, in which the velocity field is zero, andwhere *ϵ*_{T}=*κ*_{f}/*κ*_{m} with the overbar denoting the steady state.

To assess the stability of the steady-state solution, we introduce perturbations into the steady-state solution, such that

The governing equations (2.1) of the fluid layer are non-dimensionalized with the scalingswhere Ra=*R*^{2} is the fluid Rayleigh number and . The porous layer scalings follow analogously from those of the fluid layer by replacing *d* and *κ*_{f} by *d*_{m} and *κ*_{m}, respectively, where is the porous Rayleigh number.

Substituting the perturbations and non-dimensionalized variables into (2.1)–(2.2), we derive the systems(2.7)in and(2.8)in with Prandtl number *Pr*=*μ*/(*κ*_{f}*ρ*_{0}), Darcy number , , and .

## 3. Stability analysis

Although a linear instability analysis does provide an interesting insight into the onset of convection, there are potential regions of subcritical instabilities where the instability occurs prior to the thresholds predicted by the linear theory being reached. Quantifying the discrepancy between the linear instability and nonlinear stability thresholds makes it possible to provide an assessment of the suitability of linear theory to predict the onset of convection.

### (a) Linear instability analysis

The linearized equations are derived from (2.7) to (2.8) by discarding the nonlinear terms and assuming normal modes of the formAnalogous definitions apply to the porous layer (with superscript m). The growth rates *σ* and *σ*^{m} can now be used to assess whether the zero solution is unstable. If *Re*(*σ*,*σ*^{m})>0, then the perturbation will grow exponentially in time, clearly leading to linear instability.

Letting *D*=d/d*z* and taking the double *curl* of (2.7)_{1} and (2.8)_{1} to remove the pressure terms, systems (2.7) and (2.8) become(3.1)where and . The boundary conditions for the twelfth-order system (3.1) at *z*=1 areandat *z*=−1. On the interface *z*=0, we haveandThe numerical results are presented in §4.

### (b) Unconditional nonlinear stability analysis

Let us define *Ω*_{f} and *Ω*_{m} to represent the period cells in the fluid and porous layers, respectively, and introduce the notation of norm and inner product on the spaces *L*^{2}(*Ω*_{f}) and *L*^{2}(*Ω*_{m}), whereTo obtain global nonlinear stability bounds in the stability measure *L*^{2}(*Ω*_{f}), we multiply equations (2.7)_{1} and (2.7)_{3} by *u*_{i} and *θ*, and integrate over the period cell. An analogous process is applied to (2.8)_{1} and (2.8)_{3}. We may now define the functional *E*(*t*) byfor coupling parameters *λ*_{1}, *λ*_{2}, *λ*_{3}>0, such that(3.2)Two zero terms have been included in (3.2) using (2.7)_{2} and (2.8)_{2}. By integrating the first and third terms on the right-hand side of (3.2) by parts, and employing the divergence theorem, we see that(3.3)To address the boundary terms in (3.3), the boundary conditions (2.5) and (2.6) (at *z*=0) are non-dimensionalized to yieldThus, setting , and noting that at *z*=0, (3.3) shows that(3.4)Similarly, upon integration by parts and using the non-dimensionalized boundary conditionsat *z*=0, the second and fourth terms on the right-hand side of (3.2) yield(3.5)where *λ*_{1}=*λ* and

Combining (3.4) and (3.5) with (3.2), it follows that d*E*/d*t*=−, where

Introducing the maximization problem where is the space of admissible functions for solutions to equations (2.7)–(2.8), it follows that

Using Poincaré-like inequalities (cf. Payne & Straughan 1998), we can deduce that ≥*cE* for some constant *c*. Upon integration, we have as *t*→∞, where and convergence is at least exponential.

The corresponding Euler–Lagrange equations that arise at the sharpest threshold *R*_{E}=1 are(3.6)in the fluid layer and(3.7)in the porous layer, where *ω* and *ω*^{m} are the Lagrange multipliers.

By adopting normal mode representations and taking the double *curl* of equations (3.6)_{1} and (3.7)_{1} to remove the Lagrange multipliers, the twelfth-order eigenvalue problem (3.6)–(3.7) can be used to locate the critical nonlinear Rayleigh number Ra_{E}, which is given byNumerical results for the nonlinear energy approach are presented in §4.

## 4. Results and conclusions

Adopting the physical definitions of foametal C (cf. Straughan 2002; Goyeau *et al.* 2003), we define the permeability and porosity of the porous medium to be 8.19×10^{−8} m^{2} and 0.79, respectively. The remaining parameters, unless stated otherwise, are fixed at *δ*=5×10^{−6}, *G*_{m}=10, *Pr*=6 and *ϵ*_{T}=0.7. It is important to note that, for the parameter ranges explored, the spectra of the eigenvalue associated with the linear theory were found to be real (growth rates *σ* and *σ*_{m} in §3*a*). Consequently, this demonstrates that *G*_{m} and *Pr* have no impact on either the linear or nonlinear thresholds. Although the parameter ranges explored are different, the pattern of the linear instability results is consistent with previous work, which uses the Darcy–Brinkman equation (Hirata *et al.* 2007).

Figure 2 shows the neutral curves for a variation in depth ratio values, where the linear instability and nonlinear stability thresholds are represented by solid and dashed lines, respectively.

At , we see that the instability is dominated by the fluid mode. The change from dominance by the fluid layer to dominance by the porous layer is found in the range As we increase to 0.15, the porous mode entirely dominates the instability. The region of potential subcritical instabilities between the linear instability and nonlinear stability thresholds is very small.

Figure 3 shows the neutral curves for a variation in thermal conductivity ratio *ϵ*_{T} values, where the linear instability and nonlinear stability thresholds are represented by solid and dashed lines, respectively.

For *ϵ*_{T}=0.5, the porous mode dominates, although the fluid mode is still present. As *ϵ*_{T} increases from 0.5 to 0.9, the fluid mode moves down and eventually dominates the instability. The region of potential subcritical instabilities between the linear instability and nonlinear stability thresholds is again shown to be very small.

Since the linear instability and nonlinear stability results clearly show excellent agreement, we can conclude that the linear theory accurately encapsulates the physics of the onset of convection.

## Acknowledgments

This work was supported by a research project grant of the Leverhulme Trust (grant no. F/00128/AK).

## Footnotes

- Received July 21, 2008.
- Accepted September 10, 2008.

- © 2008 The Royal Society