## Abstract

A thermal convection model is considered that consists of a layer of viscous incompressible fluid contained between two horizontal planes. Gravity is acting vertically downward, and the fluid has a density maximum in the active temperature range. A heat source/sink that varies with vertical height is imposed. It is shown that in this situation there are three possible (different) sub-layers that may induce convective overturning instability. The possibility of resonance between the motion in these layers is investigated. A region is discovered where a very sharp increase in Rayleigh number is observed. In addition to a linearized instability analysis, two global (unconditional) nonlinear stability thresholds are derived.

## 1. Introduction

Penetrative convection, where part of a fluid layer becomes thermally unstable and then the ensuing convective fluid motion ‘penetrates’ into the otherwise stable layer(s), is a subject with immense application. Straughan (2004*a*) mentions several applications in geophysics, Mharzi *et al.* (2000) dealt with applications in building design and heat transfer, Zhang & Schubert (2000, 2002) are concerned with applications in astrophysics, and a very interesting application is highlighted in a study by Kaminski *et al.* (2011), where they point out that penetrative convection above large lava flows may be responsible for enhancing the rise of volcanic plumes into the Earth's stratosphere.

Mathematical models to describe penetrative convection, based on either an internal heat source or sink, or using a nonlinear density–temperature relationship in the buoyancy term, have been developed and analysed intensely. In particular, recent studies of linear instability and of nonlinear stability have been developed by Bennacer *et al.* (1993), Capone *et al.* (2010, 2011*a*,*b*), Carr (2004), Carr & Straughan (2003), Gobin & Bennacer (1996), Gobin & Benoit (2012), Gobin *et al.* (2005), Hill (2004, 2005), Krishnamurti (1997), Larson (2000, 2001), Normand & Azouni (1992), Papanicolaou *et al.* (2011), Payne & Straughan (2000), Shivakumara *et al.* (2010), Straughan (2002*a*,*b*, 2004*a*), Straughan & Walker (1996) and Zhang & Schubert (2000, 2002).

In addition to penetrative convection, there have been some very interesting recent analyses of resonance in thermal convection. By this, we mean that instability in one part of a fluid layer may simultaneously occur with instability in another part of the layer. Such resonance can lead to unusually high Rayleigh numbers at the onset of thermal convection and so this is very much of interest to the heat transfer industry. Typically, a high Rayleigh number may be synonymous with delaying or prohibiting heat transfer and is thus important in insulation, while low Rayleigh numbers may be desirable when one requires rapid heat transfer such as in cooling pipes used in many modern devices such as computers. The implications of resonance in the energy sector may be important, especially with nano-devices (see e.g. Duan (2012) and Mharzi *et al.* (2000)); so we deem a mathematical analysis of such resonance to be useful.

Narasimhan & Reddy (2011) investigated the possibility of thermal resonance in a bidispersive porous medium where resonance is between small-scale perturbations in the porous blocks and large-scale behaviour in the open fluid. Normand & Azouni (1992) showed that if one used a constant heat source/sink together with a quadratic density–temperature relationship in the fluid, then two-layer resonance was possible, and greatly increased Rayleigh numbers were observed in an appropriate parameter range where convection switched from stationary convection on two separate real branches to a range where oscillatory convection dominated. Straughan (2004*b*) analysed an analogous situation to that of Normand & Azouni (1992) but in a Darcy porous medium. He found that the stationary–oscillatory–stationary structure of Normand & Azouni (1992) still prevailed, but the oscillatory branch was never the one where convection actually occured, thereby demonstrating that the porous case was very different from the pure fluid one.

The objective of this study is to consider thermal convection in a plane layer when the density–temperature relation in the buoyancy term is quadratic. However, we allow a heat source/sink that varies linearly with the vertical height *z*. This function is a heat sink in part of the layer and a heat source in the rest. This has some resemblance to the experimental situation of Krishnamurti (1997), although it is different, and appropriate experiments would be appreciated. We derive linear instability results and global nonlinear stability (for all initial data) and we show that there is a very interesting parameter range where the Rayleigh number for instability rises very rapidly. The basic steady-state configuration of interest is one where there are three potential sub-layers which may give rise to instability that can then penetrate into the other layer(s). With the modern need for heat transfer or insulation devices in industry, especially connected with nanotechnology, cf. Duan (2012) and Mharzi *et al.* (2000), we believe the study enclosed herein is of value.

## 2. Equations of motion

The balance of linear momentum equation has the form
2.1where **v** and *p* denote velocity and pressure, and standard indicial or vector notation is used throughout. The quantities *ν*,*g*,*ρ*_{0} are kinematic viscosity, gravity and the constant density, where only density variations in the buoyancy term are considered. The vector **k**=(0,0,1), Δ is the Laplacian and equation (2.1) holds in the layer for all positive time, *t*>0. The density function *ρ*(*T*) has the form
2.2where *α* is a thermal expansion coefficient, and *T*_{m} is that value of temperature where the density *ρ* achieves a maximum. (Several real fluids display a density maximum, perhaps the most common being water, where *T*_{m}≈4^{°}*C*.)

The balance of mass equation is
2.3and the balance of energy equation is
2.4with equations (2.3) and (2.4) holding on . In equation (2.4), *T* is temperature, *κ* is thermal diffusivity and
2.5where *a*,*b* are constants to be defined later.

We wish to consider a steady-state solution to equations (2.1)–(2.5), where the velocity and the temperature is a cubic in *z* as in figure 1. (We have considered *T*_{0}<*T*_{m}, but we could equally well analyse the analogous situation where *T*_{0}>*T*_{m}.) The points *d*_{1},*d*_{3} and *h* are where , and *d*_{2} and *d*_{4} are where . In this situation, there are three potentially unstable layers where thermal convection may arise, namely (0,*d*_{1}), (*d*_{2},*d*_{3}) and (*d*_{4},*h*), this being due to relation (2.2), which ensures the fluid density is maximum at *z*=*d*_{1},*d*_{3} and *h*. Thus, in principle, there is a possibility of resonant convection between each or all three of these layers. The study of this mathematical and physical problem is the goal of this work.

In the steady state, , and then equation (2.4) reduces to solving
with solution
Let us denote Δ*T*_{1}=*T*_{m}−*T*_{0}, then one finds the solution that matches the scenario of figure 1 is given when
2.6
2.7
and
2.8where *ζh*=*d*_{1} and *ξh*=*d*_{3}. This leads to the temperature profile
2.9The turning points *d*_{2} and *d*_{4} are found to be given by
2.10It is worth observing that with the values of *a*,*b* given by (2.6) and (2.7), the heat source/sink satisfies
while
thus we have a heat source in the lower part of the layer, whereas it is a heat sink in the upper.

To study nonlinear stability and linear instability of the basic steady-state solution above, we introduce perturbations *u*_{i},*θ*,*π* by , and , where is that function (up to an appropriately chosen constant) which arises from equation (2.1). Length, velocity and temperature scales are chosen as , and the Prandtl and Rayleigh numbers are defined by
2.11Let now and define the non-dimensional temperature function *F*(*z*) by
2.120≤*z*≤1. One then derives the non-dimensional perturbation equations as
2.13where equations (2.13) hold on , with *w*=*u*_{3}, and where *F*^{′}=d*F*/d*z*. The boundary conditions to be satisfied are
2.14with *u*_{i},*θ*,*π* satisfying a plane-tiling periodicity in the horizontal plane.

## 3. Linear instability

In order to study linear instability, we discard the quadratic terms in (2.13)_{1} and (2.13)_{3}. A time dependence such as *u*_{i}=e^{σt}*u*_{i}(**x**), *θ*=e^{σt}*θ*(**x**), *π*=e^{σt}*π*(**x**) is now assumed and then, after removing, the pressure perturbation the linearized instability equations that arise from (2.13) are found to be
3.1where Δ*=∂^{2}/∂*x*^{2}+∂^{2}/∂*y*^{2}, and (3.1) hold on .

When looking for resonance, one might expect such behaviour when the three potentially unstable layers have the same depths, i.e. when
3.2Additionally, similar resonance may be expected if the Rayleigh numbers pertaining to each layer have approximately the same values. From equation (2.11), we see that *Ra*∝(Δ*T*_{1})^{2}*h*^{3}, and so Rayleigh numbers for each of the layers may be defined by
3.3where *T*_{eu} and *T*_{el} indicate extreme upper and lower values of achieved when *z*=*d*_{2} and *z*=*d*_{4}, respectively. The relations *d*_{1}=*ζh*, *d*_{3}=*ξh* with *d*_{2},*d*_{4} given by (2.10) are used, and a numerical search was performed to determine whether (and when) the depths are approximately equal, and/or the Rayleigh numbers are equal as in (3.2) and (3.3). Our numerical search indicated that we should concentrate on values around *ζ*=0.128,*ξ*=0.532 and *ζ*=0.191,*ξ*=0.531. (A referee has kindly pointed out that equations (3.2) have the unique analytical solution . This agrees with the numerically computed pair.) This gives a very useful guide as to where we might expect complex growth rates and possible resonance. Such a guide is invaluable, because we find that the eigenvalues are very sensitive to parameter changes and one must compute with precise values of *ζ* and *ξ*. Thus, a mass search over all possible *ζ* and *ξ* values would be overly time-consuming and we compute very carefully in the neighbourhood of the (*ζ*,*ξ*) pairs given earlier.

To proceed from equations (3.1), a plane tiling form *f* is introduced, see e.g. Straughan (2004*a*, p. 51), and then we put *w*=*W*(*z*)*f*(*x*,*y*), *θ*=*Θ*(*z*)*f*(*x*,*y*) and introduce the wavenumber *a* by Δ**f*=−*a*^{2}*f*. Equations (3.1) then yield the eigenvalue problem
3.4In (3.4) *D*=d/d*z*, *z*∈(0,1), and the boundary conditions are
3.5The system (3.4), (3.5) is solved by a *D*^{2} Chebyshev tau method, cf. Dongarra *et al.* (1996). Detailed numerical results are presented in §5.

## 4. Global nonlinear stability

Linearized instability theory yields a threshold above which we can be sure that instability occurs. However, it does not yield any information on whether the solution is stable below this parameter threshold. One way of achieving a stability threshold is to use nonlinear energy stability theory, cf. Straughan (2004*a*). This theory is particularly useful if the stability so obtained is global, i.e. holds for all values of the initial perturbations, or at least holds for sufficiently large values. Such non-vanishing values are usually difficult to obtain for equations as complicated as (2.13), (2.14). In this study, we use a weighted energy method, and additionally use a generalized energy method involving *L*^{2} and *L*^{3} integrals in the energy (Lyapunov) function, to achieve global stability thresholds. We might point out that nonlinear energy stability methods have been increasingly applied in penetrative convection and in other thermal convection studies, cf. Capone *et al.* (2010, 2011*a*,*b*), Carr (2004), Falsaperla *et al.* (2012), Falsaperla & Mulone (2010), Hill (2004, 2005), Hill & Carr (2010), Hill & Malashetty (2012), Larson (2000, 2001), Lombardo *et al.* (2008), Papanicolaou *et al.* (2011), Payne & Straughan (2000), Rionero (2009, 2011, 2012), Shivakumara *et al.* (2010), Straughan (2002*a*,*b*, 2004*a*,*b*) and Straughan & Walker (1996).

### (a) Weighted energy method

With the standard energy technique, one would multiply equation (2.13)_{1} by *u*_{i} and integrate over a period cell *V* of the nonlinear convection disturbance. For equations (2.13), this presents a problem because one is left with a cubic term of form to estimate. To overcome this, we use a weighted energy functional. Let ∥·∥ and (·,·) denote the norm and inner product on *L*^{2}(*V*) and let 〈·〉 denote integration over *V* . The energy we now use is
4.1where , *μ*>2 being a constant we select optimally.

Upon differentiation of *E*, the use of (2.13) and (2.14) and some integrations by parts one finds
4.2where
4.3
4.4and
4.5To proceed from (4.2), put
4.6where is the space of admissible solutions, i.e. *u*_{i}∈*H*^{1}(*V*), *θ*∈*H*^{1}(*V*), with *u*_{i,i}=0, and with the boundary conditions satisfied. Then, from equation (4.2),
4.7If now *R*<*R*_{w}, then by using Poincaré's inequality, one may determine a constant *c*>0 such that *D*≥*cE* and then one deduces exponential decay of *E* for all initial data. Thus, *R*_{w} represents a global nonlinear stability threshold. The Euler–Lagrange equations that arise from (4.6) yield an eigenvalue problem for *R*_{w}, and these are found to be
4.8The Lagrange multiplier is removed from (4.8)_{1} by taking curlcurl and keeping the third component of the result. This leads to the system
4.9The functions *w* and *θ* are written as *w*=*W*(*z*)*f*(*x*,*y*), *θ*=*Θ*(*z*)*f*(*x*,*y*), and this results in having to solve the following eigenvalue problem with wavenumber *a*, for *R*_{w}:
4.10The boundary conditions that apply to *W*,*Θ*, are as in (3.5). This system is solved by a *D*^{2} Chebyshev tau method, and numerical output is discussed in §5.

### (b) A generalized energy method

One may argue whether the form of the Navier–Stokes equations given in (2.1) is actually adequate for a nonlinear analysis. Indeed, Straughan (2002*a*) adapts several models of Ladyzhenskaya to be applicable to thermal convection. In addition, a recent study of Bulíček *et al.* (2007) likewise argues for a nonlinear dependence of stress upon the symmetric part of the velocity gradient. Hence, we now consider a modification of (2.1), which is consistent with model III of Ladyzhenskaya, as described by Straughan (2002*a*), but adapted to the penetrative convection problem in hand. Thus, we consider equation (2.1) modified to
4.11where is a constant. Equations (2.2)–(2.5) still hold as does the steady solution in §2 and the linearized analysis in §3. By adopting equation (4.11), one now derives, instead of the nonlinear perturbation equations (2.13), the system
4.12where *ν*_{1} is a non-dimensional form of . The boundary conditions are as in (2.14).

To develop a nonlinear stability analysis from equations (4.12), we construct the identities
4.13
4.14
and
4.15where ∥·∥_{3} denotes the *L*^{3}(*V*) norm. We now let *λ*>0 and *a*_{1}>0 be constants to be selected and consider the energy function *E*,
4.16and the generalized energy function, ,
4.17Upon using equations (4.13)–(4.15), we may derive
4.18where now
4.19and
4.20

We now use the Cauchy–Schwarz and Poincaré inequalities to find
and use these in (4.18) after estimating the *wθ*^{2} terms with Young's inequality. In this way, we may see that
4.21where *α*>0 is a constant to be chosen and denotes the maximum of *F*^{′} on (0,1), i.e.
(provided *ζ*+*ξ*>1/2). From inequality (4.21), we may deduce
4.22

To progress beyond inequality (4.22), we require
4.23and
4.24Choose now *α* such that
Then from (4.23), we may obtain
4.25Optimize the right-hand side of (4.25) to obtain . This then leads to the restriction on *ν*_{1},
4.26and we note that . With the above choices, the coefficient of in (4.22) becomes
Then from inequality (4.22), we may find
4.27where
4.28Collecting together terms, one finds
and then with use of Poincaré's inequality
Put *c*=1−*R*/*R*_{E} and then from (4.27), if *R*<*R*_{E},
4.29Exponential decay of may be deduced from (4.29) for all initial data, provided *ν*_{1} satisfies condition (4.26), and thus global nonlinear stability holds when *R*<*R*_{E}.

The number *R*_{E} is found by maximizing
4.30where and
4.31The Euler–Lagrange equations for are found to be
4.32and by eliminating the Lagrange multiplier *ω*, we find
4.33With the representations *w*=*W*(*z*)*f*(*x*,*y*), *ϕ*=*Φ*(*z*)*f*(*x*,*y*), this leads to the eigenvalue problem
4.34to be solved for *R*_{E} subject to boundary conditions (3.5) with *Θ* replaced by *Φ*. This system is solved numerically by the *D*^{2} Chebyshev tau method, and numerical output is discussed in §5.

## 5. Numerical results and conclusions

As stated in §2, a calculation matching depths and Rayleigh numbers for the three potentially unstable layers suggested *ζ*=0.128,*ξ*=0.532 and *ζ*=0.191,*ξ*=0.531. We have computed critical Rayleigh numbers and wave numbers for many combinations of *ζ* and *ξ*, 0<*ζ*<*ξ*<1. Our computations are based on solving equations (3.4) and (3.5) for the linear instability thresholds; equations (4.10) and (3.5) for the weighted energy thresholds; and equations (4.34) and (3.5) for the generalized energy thresholds. For the linear instability thresholds, we concentrate mainly on reporting values for *ζ* near 0.16 and 0.17 with *ξ* near 0.51–0.53. This is the range of values in (*Ra*,*ζ*,*ξ*) space where we found a very rapid rise in Rayleigh number, and the possibility of oscillatory convection with the major component of convection switching between the three basic layers, and penetrative convection ensuing.

In tables 1–3, we report critical instability values based on linearized theory. For *ζ*,*ξ* in the ranges [0.1626,0.1628] and [0.51,0.53], we find that the *Ra* versus *a* curves display behaviour like that of Normand & Azouni (1992), figures 2*f*–*h*, or see also Straughan (2004*b*), figures 3–5; however, we stress that the situation here is more like that of Normand & Azouni (1992), where the oscillatory convection branch is lowest, leading to this kind of instability at onset. For example, in table 1, we display minimum *Ra* values for each of the three eigenvalue branches, corresponding to different wavenumbers. Each branch is U-shaped and we report the minimum on each branch. Thus, the actual instability in table 1 is achieved by stationary convection with *a*∈[5.89,5.91], and so with larger width cells. In table 2, we again have three eigenvalue branches, but now the minimum is on the middle one, which is a complex branch leading to oscillatory instability, with *a*∈[8.44,8.46], and a cell that is less wide than those pertaining to table 1. In table 3 when *ξ*=0.51, the real branch on the left appears to disappear, and we find only a complex branch, then a real branch with larger critical wavenumber. The onset of convection is again via oscillatory convection. However, as tables 2 and 3 show, as *ζ* increases the nature of convection changes to stationary convection, even though the topology of the *Ra*−*a* curves is still based on a complex and a real branch. Then the wavenumber is larger; so the cells are even narrower. When *ξ* decreases very slightly, there are still only two eigenvalue branches, a complex one and a real one to the right, and the critical Rayleigh number may be even higher. For example, we computed
For other combinations of *ζ*,*ξ*, we did not discover oscillatory convection, and the Rayleigh numbers we computed were much lower. For example, with *ξ*=0.8, *ζ*∈[0.4,0.49] *Ra* varies from 15914 to 10911, with *a* varying from 4.5 to 4.0.

The critical Rayleigh numbers for *ξ*=0.52 with *ζ* varying are displayed in figure 2, with greater detail in figure 3 of the region where instability changes from stationary to oscillatory and then to stationary again.

The behaviour of the eigenfunctions is displayed in figures 4–13. The *W* eigenfunction leads to streamlines, although these are not displayed explicitly. Figure 4 corresponds to the minimum on the real eigenvalue branch to the left of the complex branch when *ξ*=0.52,*ζ*=0.1627, with *a*=5.26. Thus, this *W* function is not the one realized, but, figures 5 and 6 are realized being on the same branch with *ζ* and/or *ξ* slightly shifted. In each case, the three instability zones are clearly displayed with the middle one being the dominant zone. We also clearly see penetrative convection with a two-cell structure. Figure 7 displays the real and imaginary parts of *W* on the complex branch when *ξ*=0.52,*ζ*=0.1627,*a*=8.45. The solution will oscillate between these two positions. The real part shows again penetrative convection with dominance in the second potentially unstable layer, while the imaginary part actually shows four cells with penetrative convection, with the dominant motion being in the layer nearest to the upper plate. The real part is, in some ways, a natural extension from figures 5 and 6. However, once we move to the minimum on the real branch with larger wavenumber, *ξ*=0.52, *ζ*=0.1627, *a*=14.12, the imaginary part as in figure 7 is not pursued, as seen in figure 8. While the eigenfunction in figure 8 will not be realized in practice, the same shape is realized, for nearby values of *ζ*,*ξ*, as seen in figures 9 and 10. Here, the dominant layer for *W* is the one nearest the lower plate. Again, penetrative convection is observed but now with a three-cell structure; the cells following the pattern of the potentially unstable layers outlined earlier.

Finally, figures 11–13 display the temperature eigenfunction. Figure 11 corresponds to that of figure 4, and we see a behaviour that is not dissimilar to that of *W* there. Figure 12 shows the real and imaginary parts of *Θ* corresponding to the *W* values in figure 7. The solution will oscillate between these values, and a triple layer structure is evident. Figure 13 shows *Θ* corresponding to *W* in figure 8. The temperature structure has a clear three-layer pattern and is in accordance with the corresponding value of *W*.

Table 4 displays values of the linearized critical values, *Ra*_{L},*a*_{L} and *σ*_{i}, together with the global nonlinear stability values *Ra*_{E}, *a*_{E} and *λ*, and *Ra*_{w}, *a*_{w} and *μ*. The nonlinear thresholds clearly display a potentially large region where subcritical instabilities may arise. However, we have computed these in the (*ζ*,*ξ*) region where a very rapid rise in the Rayleigh number is observed. A similar discrepancy in linear versus nonlinear thresholds was observed earlier by Straughan (2004*b*) in a porous convection model, although the situation there was not as complicated as the present. In future work, we intend solving numerically the actual convection equations in the region of possible subcritical instabilities, and also once *Ra*_{L} is exceeded to actually calculate heat transfer values.

## Acknowledgements

This work was supported by the Leverhulme research grant, ‘Tipping Points, Mathematics, Metaphors and Meanings’. I am indebted to three anonymous referees for their pointed remarks that have led to improvements in the manuscript. In particular, one referee pointed out the analytical solution to (3.2).

- Received April 8, 2012.
- Accepted July 24, 2012.

- This journal is © 2012 The Royal Society