## Abstract

This paper explores an operative technique for deriving nonlinear stability by studying double-diffusive porous convection with a concentration-based internal heat source. Previous stability analyses on this problem have yielded regions of potential subcritical instabilities where the linear instability and nonlinear stability thresholds do not coincide. It is shown in this paper that the operative technique yields sharp conditional nonlinear stability in regions where the instability is found to be monotonic. This is the first instance, in the present literature, where this technique has been shown to generate sharp thresholds for a system with spatially dependent coefficients, which strongly advocates its wider use.

## 1. Introduction

Two of the key techniques employed in a stability analysis are the nonlinear energy method and the more widely used method of linearized instability. A linear instability analysis, by definition, only provides boundaries for instability due to the presence of nonlinear terms (cf. 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 behaviour of an instability. Another key advantage of performing a nonlinear energy analysis is that it provides rigorous conditions for stability without exploiting linearization or weakly nonlinear approximations.

Nonlinear stability analyses have been performed successfully in a variety of applications, see Hill & Carr (2010) and Singh (2010). However, for a considerable range of problems, the standard energy approach does not generate sharp results (i.e. where the linear and nonlinear thresholds coincide), see Carr (2003), Sunil & Mahajan (2008) and Saravanan & Brindha (2010). This may, of course, be due to the significant possibility that the linear theory has not predicted the onset of instability accurately; however, there are several new nonlinear energy techniques in development which show that sharp thresholds can still be achieved.

In the field of porous media convection, a differential constraint approach (cf. van Duijn *et al.* 2002) has been introduced, where the fluid-governing Darcy equation is retained as a constraint. This has been shown to reduce the number of required coupling parameters (cf. Hill 2009) and to provide sharp nonlinear thresholds in contrast to the standard approach, see Capone *et al.* (2010).

In a separate recent development, an operative technique (for systems with constant coefficients) was introduced by Mulone (2004) that involves a change of variable in the energy method, based on the eigenvalues of the linearized system. In this paper, which was concerned with porous media convection, the coincidence of the critical linear and nonlinear stability parameters was proved. Further applications of this technique have seen sharp thresholds being produced for biological systems with diffusion (Mulone & Straughan 2009).

To assess the feasibility of an operative approach that could include systems with spatially dependent coefficients (which extends the use of the technique considerably), Hill (2008) further developed the method in the context of penetrative double-diffusive convection. This is an interesting test case, as the standard energy method approach yields a nonlinear threshold that is independent of the salt field. By adopting this new approach, a nonlinear threshold that was dependent on the salt field was shown, although coincidence between the linear and nonlinear thresholds was not achieved.

It is important to note that the operative technique yields regions of conditional stability, for which the initial data must be bounded by a threshold. Conditional stability does not necessarily preclude subcritical instabilities if this restriction on the initial data is not considered.

The purpose of this paper is twofold. Firstly, to improve on the design of the method from Hill (2008), now that its feasibility has been established, by basing the transformation solely on the linear eigenvalues, making the method much easier to solve numerically and, therefore, more widely useable. Secondly, to apply it to a physical system where all previous attempts to generate sharp nonlinear thresholds have not been successful in regions where the instability is monotonic. The technique, at present, has not been developed to be applicable in regions of oscillatory instability. The energy method itself (which underpins all of the nonlinear techniques) does not tend to generate sharp thresholds in such regions (cf. Straughan 2004). The operative technique will be used in conjunction with the differential constraint approach, which is a new development in the current literature.

Due to the competing temperature and salt gradients, double-diffusive convection poses a challenging nonlinear stability problem, and is therefore a good class of problems to test the method on. The introduction of a concentration-based internal heat source yields a system that forms an accurate model of cumulus convection as it occurs in the atmosphere (Krishnamurti 1997), and can also be interpreted as a simplified model of the salt gradient layer in a solar pond (Kudish & Wolf 1978). Although there has been considerable interest in this problem, see Krishnamurti (1997), Straughan (2002), Hill (2003, 2005) and Chang (2004), the coincidence of the linear and nonlinear thresholds for the double-diffusive case has not been achieved, even with the application of the differential constraint approach (Hill 2009). In this paper, we apply our operative technique to this problem (based on the results of Hill 2005), the purpose of which is to generate sharp nonlinear thresholds.

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.

## 2. Formation of the problem

Let us consider the same configuration as Hill (2005), where a fluid-saturated porous layer is contained between the horizontal parallel planes *z*=0 and *L*. O*xyz* is the standard Cartesian frame of reference with unit vectors **i**,**j** and **k**, respectively. The model is of double-diffusive convection induced by the selective absorption of radiation, where the convection mechanism is essentially a penetrative one effectively modelled via an internal heat source.

The Darcy equation is assumed to govern the fluid motion in the layer, such that
2.1
where **v**=(*u*,*v*,*w*) and *p* are velocity and pressure, *g* is the acceleration due to gravity, *μ* is the dynamic viscosity of the fluid and *K* is the permeability.

Denoting *T* to be the temperature and *C* to be the concentration of the dissolved species, the density *ρ*(*T*,*C*) is given by
where *ρ*_{0}, *T*_{0} and *C*_{0} are reference values of density, temperature and concentration, respectively, and *α*_{t} and *α*_{c} are the coefficients of thermal and solutal expansion.

Equation (2.1), together with the incompressibility condition and the equations of energy and solute balance, yield the following system of governing equations:
2.2
2.3
2.4
and
2.5
The internal heat source is modelled linearly with respect to concentration, which is represented by the introduction of the *βC* term in equation (2.4), where *β* is some constant of proportionality. In these equations, *ε* is the porosity, *κ*_{t} and *κ*_{c} are the thermal and solutal diffusivities, and *M*=(*ρ*_{0}*h*_{p})_{f}/(*ρ*_{0}*h*)_{m} where (*ρ*_{0}*h*)_{m}=(1−*ε*)(*ρ*_{0}*h*)_{s}+*ε*(*ρ*_{0}*h*_{p})_{f}, *h*_{p} is the specific heat of the fluid and *h* is the specific heat of the solid, with the subscripts ‘f’, ‘s’ and ‘m’ referring to the fluid, solid and porous components of the medium, respectively.

The boundary conditions for the problem are **v**=**0**, *T*=*T*_{L} and *C*=*C*_{L} at *z*=0 and **v**=**0**, *T*=*T*_{U}<*T*_{L}, *C*=0 at *z*=*L*.

Assuming a steady-state solution (, ) that corresponds to no fluid flow (i.e. ), and introducing a perturbation to the steady state of the form we may derive the governing non-dimensionalized perturbation equations, 2.6 2.7 2.8 and 2.9

In these equations, *Le*=*κ*_{t}/*κ*_{c} is the Lewis number, is the normalized porosity, *τ*=*βc*_{L}*L*^{2}/(*κ*_{t}(*T*_{L}−*T*_{U})) is a measure of the internal heat source generated by the radiation absorbing concentrate and *R*_{T}^{2}=*α*_{t}*gρ*_{0}(*T*_{L}−*T*_{U})*KL*/(*μκ*_{t}) and *R*^{2}_{S}=*α*_{c}*gρ*_{0}*C*_{L}*KL*/(*μκ*_{t}) are the thermal and solute Rayleigh numbers, respectively. The term *F*(*z*)=1−*τ*/3+*τ*(*z*−*z*^{2}/2) is the non-dimensionalized temperature gradient. Although a slightly different non-dimensionalization is used in this paper, a detailed derivation of equivalent governing equations and explicit parameter definitions may be found in Hill (2005).

The boundary conditions that follow from the perturbed quantities are *w*=0,*θ*=0 and *ϕ*=0 at *z*=0,1, where (**u**,*θ*,*ϕ*) have a periodic plan-form tiling the (*x*,*y*) plane, and *Ω* is the period cell for the perturbations.

## 3. Linear analysis

The linearized equations are derived from equations (2.6)–(2.9) by discarding the nonlinear terms. Introducing normal modes of the form
where is the growth rate, and the wavenumber is a measure of the width of the convection cells that form at the onset of instability, and taking the double *curl* of equation (2.6) to remove the pressure term, system (2.6)–(2.9) yields
3.1
3.2
and
3.3
with boundary conditions *w*=*θ*=*ϕ*=0 at *z*=0,1. If *Re*(*σ*)>0, then the perturbation will grow exponentially in time, clearly leading to linear instability. The sixth-order system (3.1)–(3.3) was solved using the Chebyshev-tau method (cf. Dongarra *et al.* 1996), which is a spectral technique coupled with the QZ algorithm. Numerical results for the linear theory are presented in §5.

## 4. An operative method approach to derive sharp nonlinear thresholds

The operative method approach that is used in this work is a development of the argument presented in Mulone (2004) (and later in Mulone & Straughan 2006) to address systems where the coefficients have spatial dependence. The presence of spatially dependent coefficients is highly common in double-diffusive porous convection problems, which makes this technique particularly applicable. A version of this technique was first applied to a simple penetrative double-diffusive convection problem in Hill (2008) to assess its feasibility.

### (a) Transforming the nonlinear system

Our aim is to construct a transformation associated, in a canonical way, to the linear system via the first eigenvalue of the Laplacian operator that we can apply to our full nonlinear system (2.6)–(2.9).

However, in systems with spatially dependent coefficients (in our case, the *F*(*z*) term in equation (3.2)), the eigenvalues of the Laplacian operator cannot be directly evaluated. To overcome this, we construct the transform based on the eigenvalues of a modified version of linearized equations. For equations (3.1)–(3.3), this is achieved by replacing *F*(*z*) with the average of its value over the range *z*∈[0,1] (namely 1).

Under this modified system, all the even derivatives of *w*, *θ* and *ϕ* now vanish at the boundaries, so we may take solutions of the form
It can easily be shown that the thermal Rayleigh number of this modified linearized system is minimized by the mode *p*=1. Since the sin functions are linearly independent, we may study only the most unstable mode *p*=1. Letting *s*=*a*^{2}_{L}+*π*^{2} and equation (3.1) may be written
4.1
where and *R*^{2}_{SL} are the wavenumber, and thermal and solute Rayleigh numbers associated with the modified linear system, respectively. Using equation (4.1), the remaining linearized equations (3.2) and (3.3) may be written in the matrix form
where **X**=(*θ*,*ϕ*)^{T} and
The eigenvalues of matrix *M* are given by
4.2
where
Under the condition
4.3
the matrix *Q* of eigenvectors associated with eigenvalues equation (4.2) can be constructed, such that
and
Introducing the new variables
we can now transform the full nonlinear system (2.6)–(2.9). Note that as all the values in the transformation are taken from the linear system (which are derived first), matrices *Q* and *Q*^{−1} contain only real numbers.

After making this substitution (and taking the third component of the double curl of equation (2.6) to remove the pressure term), we have the system
4.4
4.5
and
4.6
where Δ*=∂^{2}/∂*x*^{2}+∂^{2}/∂*y*^{2} and
with boundary conditions *w*=*γ*=*ψ*=0 at *z*=0, 1.

### (b) Deriving the sharp nonlinear stability bounds

To obtain sharp nonlinear stability bounds in the stability measure *L*^{2}(*Ω*), (following an analogous argument to Lombardo *et al.* (2001)), we multiply equation (4.5) by *γ* and Δ*γ*, respectively, and equation (4.6) by *ψ* and Δ*ψ*, respectively, and integrate over *Ω* to obtain
where ∥⋅∥ and 〈⋅ ,⋅〉 denote the norm and inner product on *L*^{2}(*Ω*). Adopting the differential constraint approach of van Duijn *et al.* (2002) and Hill (2009), equation (4.4) is retained as a constraint. This is the first instance, in the current literature, in which the differential constraint approach has been used in conjunction with an operative method.

Letting
where *λ*>0 is a coupling parameter, and *b* is a positive constant, we can write
4.7
where
To ensure that (i.e. 1+*B*_{1}>0 and 1+*B*_{2}>0), we have the condition
4.8
Note that if equation (4.8) holds, then condition (4.3) follows.

By defining
where is the space of admissible functions, it follows that and for a positive constant *p*_{1} (see appendix A), where *R*_{E}>1.

From equation (4.7), we see that
If *R*_{E}>1, then by the Poincaré inequality (e.g. *π*^{2}∥*γ*∥^{2}≤∥∇*γ*∥^{2}), for some positive constant *p*_{2}. Hence, it follows that
4.9
We will show that equation (4.9) ensures conditional nonlinear stability as long as
4.10a
and
4.10b
For equation (4.10b) to hold, there are two possibilities as follows:

*E*^{1/2}(*t*)<1/*p*_{1},*t*≥0,

or

there exists

*ζ*>0 such that*E*^{1/2}(*ζ*)=1/*p*_{1}, with*E*^{1/2}(*t*)<1/*p*_{1}, for*t*∈[0,*ζ*).

Suppose (ii) holds. From equation (4.9), it follows that d*E*/d*t*≤0 for *t*∈[0,*ζ*), hence, *E*^{1/2}(*t*)≤*E*^{1/2}(0)<1/*p*_{1} for any *t*∈[0,*ζ*). As *E*(*t*) is a continuous function of *t* on [0,*ζ*], then it is impossible that *E*^{1/2}(*ζ*)=1/*p*_{1}. This contradiction implies (i) and consequently d*E*/d*t*≤0 for *t*≥0 such that *E*^{1/2}(*t*)≤*E*^{1/2}(0). Thus,
Integrating, we have
which shows the decay of *γ* and *ψ*. Using equation (2.6),
4.11
The decay of **u** clearly follows.

It is important to note that we have derived conditional stability (under the restriction (4.10b)), such that the initial data *E*^{1/2}(0) must be bounded by 1/*p*_{1}. From the definition of *p*_{1} (see appendix A) as *R*_{E}→1, the bound on the initial data tends to 0. However, for , we may take *b*=0, which yields unconditional nonlinear stability, such that there is no bound on the initial data.

Following the differential constraint approach of van Duijn *et al.* (2002) and Hill (2009) (where we retained equation (4.4) as a constraint), we introduce an Euler–Lagrange multiplier *μ* such that
The Euler–Lagrange equations for the maximization problem 1/*R*_{E} = are now given by
4.12
with boundary conditions *w*=*μ*=*γ*=*ψ*=0 at *z*=0,1.

After adopting normal modes, the eighth-order eigenvalue problem (4.12) for *R*_{E} was solved using the Chebyshev-tau method (cf. Dongarra *et al.* 1996). Under the analysis of this section, the ranges of the remaining parameters (restricted by equation (4.8)) with eigenvalues *R*_{E}>1 correspond to regions of stability. An overview of the numerical approach and the results for the operative nonlinear energy approach are presented in §5.

## 5. Results and conclusions

For fixed values of *Le*, *τ*, and *R*^{2}_{S}, the region (in terms of *R*_{T}^{2}) that yielded conditional stability under the operative method was numerically derived by adopting the following approach:

solve the linear instability problem to yield the critical values of

*R*^{2}_{TL}and*a*^{2}_{L}(noting that );if condition (4.8) is satisfied, then the nonlinear operative approach may be applied;

use the derived

*R*^{2}_{TL}and*a*^{2}_{L}values to define the constants*A*_{j},*B*_{j},*C*_{j}and the functions*F*_{j}(*j*=1,2) in equation (4.12);fix the value of

*R*_{T}^{2};minimize

*R*_{E}over the wavenumber*a*^{2}and maximize*R*_{E}over the coupling parameter*λ*by repeatedly solving eigenvalue problem (4.12) for*R*_{E};if

*R*_{E}>1, then we have stability; andrepeating steps 5 and 6 for a variety of

*R*_{T}^{2}values yields the region of stability.

Eigenvalue problem (4.12) was solved using the Chebyshev tau-QZ method (Dongarra *et al.* 1996), which is a spectral method coupled with the QZ algorithm. The results were checked by varying the number of polynomials to verify convergence. To compare our operative technique with the existing energy methods, standard nonlinear stability results are also derived (see Hill 2005 for more details).

Figure 1 presents stable and unstable regions, which have been generated numerically, in a graphical form. As the linear instability thresholds (which give the unstable regions) have been presented before for a range of parameter values in Hill (2005), figure 1 simply gives a representative example of how the operative method performs. The conditionally stable and globally stable regions correspond to the results generated by the operative and standard nonlinear stability approaches, respectively.

It is immediately obvious from figure 1 that the operative technique produces sharp thresholds. In figure 1*a*,*b*, as *R*^{2}_{S} is increased, the region between the linear and standard nonlinear thresholds becomes very substantial. This region is entirely removed by the operative approach, clearly demonstrating the highly substantial improvement this new technique brings over the existing methods.

However, in contrast to the operative method, the standard approach yields unconditional stability (such that there is no bound on the initial data). The conditional stability generated by the operative approach does not necessarily preclude sub-critical instabilities if the restriction on the initial data is not considered.

From figure 1*a*, we can see that the sharpness of the threshold is maintained in nearly all the regions where the method is valid (under condition (4.8)). This behaviour was repeated in the parameter ranges explored. It is important to note, though, that condition (4.8) only appeared to fail near to regions where the instability becomes oscillatory. The energy method itself (which underpins all of the nonlinear techniques) does not tend to generate sharp thresholds in such regions (cf. Straughan 2004).

It has been demonstrated that for a system with spatially dependent coefficients, the operative technique presented in this paper can generate sharp thresholds, which is not achieved by the standard approach. This strongly advocates its wider use, with the potential to sharpen nonlinear stability thresholds in a substantial number of problems. The stability results generated by this approach are conditional, though, such that the initial data must be bounded by a defined threshold.

## Acknowledgements

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

## Appendix A

Recall that for *R*_{E}>1 and *b*>0,
where under condition (4.8), 1+*B*_{1}>0 and 1+*B*_{2}>0.

Letting and using Young's inequality (for constants *c*_{1}, *c*_{2}>0) and noting that *C*_{2}<0, it follows that
Choosing *c*_{1} and *c*_{2} such that
which holds for 1<*ϕ* *Le*<3, we have
Using equation (4.11) yields
where
Using the Poincaré inequality and letting
we have
as required.

Concerning using the inequality sup_{Ω}|*f*|≤*c*_{6}∥Δ*f*∥, where
and *f*=0 on *z*=0,1, (cf. Straughan 2004), the Cauchy–Schwarz inequality and equation (4.11), we have
where
It is important to note that as *R*_{E}→1.

Again, using the inequality the Cauchy–Schwarz inequality and equation (4.11), we obtain
where
Hence,
where as required, with as *R*_{E}→1.

- Received February 28, 2011.
- Accepted August 22, 2011.

- This journal is © 2011 The Royal Society

This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.