## Abstract

Studies of the nonlinear stability of fluid/porous systems have been developed very recently. A two-domain modelling approach has been adopted in previous works, but was restricted to specific configurations. The extension to the more general case of a Navier–Stokes modelled fluid over a porous material was not achieved for the two-domain approach owing to the difficulties associated with handling the interfacial boundary conditions. This paper addresses this issue by adopting a one-domain approach, where the governing equations for both regions are combined into a unique set of equations that are valid for the entire domain. It is shown that the nonlinear stability bound, in the one-domain approach, is very sharp and hence excludes the possibility of subcritical instabilities. Moreover, the one-domain approach is compared with an equivalent two-domain approach, and excellent agreement is found between the two.

## 1. Introduction

The stability of thermal convection in a system composed of a fluid overlying a porous medium saturated with the same fluid remains a subject of particular attention (cf. Nield & Bejan 2006), owing to the wide range of geophysical and industrial applications.

The main approach adopted in the literature to study this configuration is that of a two-domain model, where the set of governing equations for each region are considered separately, with appropriate boundary conditions at the interface, cf. Ochoa-Tapia & Whitaker (1995), Alazmi & Vafai (2001), Goyeau *et al.* (2003), Chandesris & Jamet (2006, 2009), Nield & Bejan (2006) and Hirata *et al.* (2007*a*). These models generally use Darcy’s law (or Darcy–Brinkman for high porosities; cf. Nield & Bejan 2006) in the porous region, coupled with the Navier–Stokes equations in the fluid region. Recent contributions include Vafai (2005), Chang (2006), Mu & Xu (2007), Straughan (2008), Hill & Straughan (2009) and Hill & Carr (2010).

An alternative, less widely used approach, is that of a one-domain model, where the governing equations for both regions are combined into a unique set of equations, which are valid for the entire domain. The physical properties that change the governing equations for each region are then taken to be discontinuous functions, see e.g. Hirata *et al.* (2009*a*). Comparisons between the one- and two-domain approaches and the subsequent treatment of the interfacial region are discussed in Goyeau *et al.* (2003) and Hirata *et al.* (2007*b*). The equivalence of the one- and two-domain approaches for the Darcy–Brinkman/Navier–Stokes model is shown in Hirata *et al.* (2009*b*).

Although linear instability has been studied extensively for fluid/porous systems, the exploration of nonlinear stability has been a much more recent development. Payne & Straughan (1998) and Hill & Straughan (2009) developed the first nonlinear stability thresholds for a fluid–porous system. In both these papers, owing to the difficulties associated with incorporating the nonlinear **v**·∇**v** advection term in the Navier–Stokes equations into a stability analysis, the fluid was modelled using Stokes equations. Hill & Carr (2010) used a viscoelastic model proposed by Ladyzhenskaya (Ladyzhenskaya 1969; Straughan 2008) as an alternative to Navier–Stokes, and constructed nonlinear stability thresholds.

These previous papers on nonlinear stability adopt a two-domain approach, which restricts their models to specific systems (i.e. Stokes flow and a viscoelastic fluid) owing to the difficulty in incorporating the boundary conditions in a stability analysis. A one-domain approach overcomes this difficulty by avoiding the explicit formulation of the boundary conditions at the fluid/porous interface. However, as the physical parameters that change between the regions are modelled using discontinuous functions, this precludes the use of a nonlinear generalized energy stability analysis, as one cannot use the divergence theorem, cf. Straughan (2004).

In this paper, we approximate these discontinuous functions with exponentially converging continuous ones, essentially making the relevant physical parameters homogeneous throughout the layers, with a small heterogeneous transition zone at the interface. By adopting this approach, we show that it is possible to construct unconditional nonlinear energy stability thresholds, and compare this result with the linear instability of both the one- and two-domain approaches. A comparison of the linear/nonlinear thresholds allows for the assessment of the suitability of linear theory to predict the physics of the onset of convection, which is crucial in understanding and controlling the system under investigation. This is discussed in more detail in §3.

## 2. Formation of the problem

Consider a Newtonian fluid occupying the three-dimensional layer which saturates an underlying homogeneous porous layer The interface between the saturated porous medium and the fluid layer is at *z*=0, see figure 1.

The governing equations for both regions may be written in the form
2.1
where the variables *v*_{i}, *p*, *T*, *c*_{p}, *ρ*_{0} and *T*_{0} are the velocity, pressure, temperature, specific heat at a constant pressure and reference density and temperature values, respectively. The porosity *ϕ*, permeability *K* and thermal diffusivity *κ* take different values in the fluid and porous regions (table 1), where *m* and *f* refer to the value of the parameter in the porous and fluid region, respectively. The quantities *κ*_{m} and *G*_{m}=(*ρ*_{0}*c*_{p})_{m}/(*ρ*_{0}*c*_{p})_{f} are defined in terms of the fluid and solid components of the porous medium, such that *Q*_{m}=*ϕQ*_{f}+(1−*ϕ*)*Q*_{s}, where *Q*_{m}=*κ*_{m} or (*ρ*_{0}*c*_{p})_{m}.

In the two-domain approach, the set of governing equations for each region are considered separately and boundary conditions at the interface are used to couple the system. System (2.1) is equivalent to adopting the Navier–Stokes and the Darcy–Brinkman equations to govern the flow in the fluid domain and the porous region, respectively, cf. Hirata *et al.* (2007*b*, 2009*b*).

In this paper, we adopt the one-domain approach, where the porous layer is viewed as a pseudo-fluid, making the whole layer a continuum. System (2.1) then forms a unique set of equations, which are valid for the entire domain, avoiding the explicit formulation of boundary conditions at the interface. The physical properties that vary between the two regions are described in table 1. A comprehensive discussion of the variances and various physical attributes regarding the governing equations is given in Alazmi & Vafai (2000).

Letting discontinuous functions describe the changing parameters in table 1, linear instability thresholds can be derived, cf. Goyeau *et al.* (2003) and Hirata *et al.* (2007*b*, 2009*b*). However, this approach precludes the direct construction of a nonlinear stability threshold, as it prevents the use of the divergence theorem (which will be discussed in §3*b*). To address this issue, let us consider functions of the form
Clearly, as , the function becomes
Under this definition, we approximate the discrete variables by
Owing to the exponential nature of *F*(*F*_{1},*F*_{2}), these parameters are essentially homogeneous throughout the layers (for large *n*, which we can choose at our discretion), with a small heterogeneous transition zone at the interface (figure 2).

Fixing the temperatures at the upper and lower boundaries to be *T*_{U} and *T*_{L}, respectively, and assuming no slip at these surfaces, such that *v*_{i}=0 at *z*=−*d*_{m},*d*, governing equations (2.1) admit a steady-state solution in which the velocity field is zero and
In these equations, *ϵ*_{T}=*κ*_{f}/*κ*_{m} and with the overbar denoting the steady state and

To assess the instability of the steady state, we introduce perturbations to the steady solution of the form
and non-dimensionalize with scalings of
where Ra=*R*^{2} is the Rayleigh number. Substituting the perturbations and non-dimensionalized variables into equation (2.1) (and dropping the stars), we derive the non-dimensional perturbation equations
2.2
where *w*=*u*_{3}, is the Darcy number, *f*_{1}=*F*(1/*δ*,0), and *F*(*F*_{1},*F*_{2})=(*F*_{2}+*F*_{1}*e*^{−2nz})/(1+*e*^{−2nz}).

## 3. Stability analysis

The main impetus behind the use of a linear instability analysis on the thermal convection problem is to predict the onset of convection. To form such an analysis, it is assumed that the perturbation is small and so neglects terms of quadratic and higher order, which leads to this approach providing limited information on the behaviour of the nonlinear system. There is, therefore, the potential for regions of subcritical instabilities where the instability occurs prior to the thresholds predicted by the linear theory being reached.

To address this issue, one can establish stability results through the use of generalized nonlinear energy techniques (cf. Straughan 2004). By comparing the instability thresholds generated by the linear theory with the stability thresholds generated by the energy method, an assessment of the suitability of the linear theory to predict the onset of convection can then be made. The energy method has been shown to be highly successful in thermal convection problems, when the symmetric part of the linear operator associated with the governing equations dominates (cf. Straughan 2004).

### (a) Linear instability analysis

The linearized equations are derived from equations (2.2)_{1} and (2.2)_{3} by discarding the nonlinear terms and assuming normal modes of the form
where *a*_{1}, are the horizontal wavenumbers and Letting *D*=*d*/*dz* and taking the double curl of equation (2.2)_{1} to remove the pressure term, where the third component is chosen, leads to the linearized equations
and
where . The boundary conditions for the sixth-order system at *z*=−1 and are

Numerical results for the linear instability of this one-domain approach are presented in §4. Results for the equivalent two-domain approach to linear instability (which are taken from Hill & Straughan 2009) are also presented in §4 for comparison. As the continuous parameter approach is, essentially, an approximation to the discontinuous model, a comparison with the standard approach is crucial to evaluate if our method captures the correct stability behaviour.

### (b) Nonlinear stability analysis

To obtain global nonlinear stability bounds in the stability measure *L*^{2}(*Ω*), where *Ω* is the period cell for the perturbations, we multiply equations (2.2)_{1} and (2.2)_{3} by *u*_{i} and *θ*, respectively, and integrate over *Ω* to obtain
and

Owing to the exponential formulation taken in §2, *ϕ*, *f*_{1}, *f*_{2} and *G* are continuous functions, allowing for the manipulation of the terms containing these functions using the divergence theorem, cf. Straughan (2004).

Letting *λ* be a positive coupling parameter to be selected at our discretion, and defining an energy
we have the following identity:
where
and
Using Poincaré-like inequalities (cf. Payne & Straughan 1998; Hill & Carr 2010), it follows that for some positive constant *c*. Letting where is the space of admissible functions, yields
if *R*_{E}>*R*. Integrating, we have which tends to 0 as if *R*_{E}−*R*>0, so that convergence is at least exponential. Therefore, unconditional nonlinear stability can be established for *R*_{E}>*R*.

At the sharpest threshold *R*_{E}=*R*, the maximization problem is solved by the Euler–Lagrange equations
3.1
where *ω* is a Lagrange multiplier. Taking the double curl of equation (3.1)_{1} and adopting normal mode representations yield the sixth-order governing system
and
This system can now be used to find the critical nonlinear Rayleigh number *Ra*_{E}, such that
Numerical results for the nonlinear energy approach are presented in §4.

## 4. Numerical results

The numerical results were derived by using the compound matrix method, a technique belonging to the family of shooting methods, as described in Straughan & Walker (1996) and Straughan (2004). These results were then checked using the Chebyshev-tau technique (cf. Dongarra *et al.* 1996), which is a spectral technique coupled with the QZ algorithm. The parameters, unless stated otherwise, are fixed at *ϕ*=0.78, *δ*=5×10^{−5}, *G*_{m}=10, Pr=6 and *ϵ*_{T}=0.7.

Figure 3 shows the linear instability thresholds for a variation of *n* values, where *n* is the parameter we choose to decide the shape of the function *F*. From figure 2, we recall that increasing *n* improves our approximation of the discontinuous physical parameters in table 1.

Bimodal behaviour, which is expected for certain parameter ranges in fluid/ porous systems (e.g. Hirata *et al.* 2007*b*, 2009*b*; Hill & Straughan 2009), is captured in figure 3. In this bimodal behaviour, the fluid mode (corresponding to perturbations of large wavenumbers) relates to the case where the convective flow is mainly confined in the fluid layer, with the porous mode (corresponding to perturbations of small wavenumbers) relating to the case where the convective flow occurs in the entire porous region.

It is clear from figure 3 that for , the fluid mode is significantly affected by changes in *n*, whereas the effect on the porous mode is negligible. As the porous mode is dominant for the effect on the neutral curve caused by an increase in *n* beyond 100 is negligible. In all cases, *n* must be increased until the neutral curve converges. In addition, figure 3 suggests that *n* must be increased more substantially when studying systems for which the fluid mode may be dominant. In the simulations presented here, the neutral curves were found to converge for values of *n*≥350.

Figure 4 shows the neutral curves for a variation of depth ratio values, where the linear instability and nonlinear stability thresholds for the one-layer model are represented by solid lines, and the linear instability thresholds for the two-layer model (taken from Hill & Straughan 2009) are given by dashed lines.

In figure 4, the linear and nonlinear thresholds for the one-domain approach are so close that they are indistinguishable graphically, making the region of potential subcritical instabilities very small, demonstrating the suitability of the linear theory to predict the physics of the onset of convection. Thus, by adopting the continuous parameter approach, we can construct nonlinear stability thresholds that show excellent agreement with those of the linear instability theory.

Comparing the one-layer approach with that of the two domain in figure 4, we can see that the thresholds show good agreement, demonstrating that the continuous parameter approach captures the behaviour as modelled by the more widely used two-domain model.

Since the nonlinear thresholds are essentially the same as the linear ones (for the parameter ranges explored), an exploration of the various physical parameter effects on the stability bound may be found in those published for the linear theory, see, e.g. Hirata *et al.* (2007*b*, 2009*a*) and Hill & Straughan (2009).

## Acknowledgements

The authors would like to thank three anonymous referees whose comments have led to improvements in the paper.

- Received January 11, 2010.
- Accepted February 25, 2010.

- © 2010 The Royal Society