Royal Society Publishing

Tipping points in Cattaneo–Christov thermohaline convection

B. Straughan


A model is proposed for thermohaline convection when the heat flux satisfies a relaxation time law rather than the classical Fourier one. Here, we adopt the recent law due to Christov. That a Cattaneo-like law would be relevant in thermohaline convection in star evolution was proposed in 1995 by Herrera and Falcón. They do not develop a detailed model which we do here. We find that with the Cattaneo–Christov law, there is a transition curve (depending on the salt concentration) such that for Cattaneo numbers greater than this transition, the nature of convection changes dramatically.

1. Introduction

The propagation of thermal waves (heat waves) is a topic of great current interest. It is now realized that this is not simply a low-temperature phenomenon, but is something that has application in everyday life. For example, thermal waves are important in skin burns (Dai et al. 2008), in laser heating in the cornea or in radiofrequency heating for the elimination of cardiac arrhythmias, destruction of tumours, treatment of gastro-oesophageal reflux disease, corneal heating (Berjano et al. 2005; López Molina et al. 2008; Tung et al. 2009), in phase changes (Liu et al. 2009; Miranville & Quintanilla 2009, 2010), in biological materials (Saidane et al. 2005), in planetary and stellar evolution (Herrera & Falcón 1995; Falcón 2001; Bargmann et al. 2008) and in nanofluids (Vadasz 2005, 2006; Vadasz et al. 2005; Straughan 2008, 2010ac). In particular, Vadasz et al. (2005) question whether traditional treatments of thermal conductivity measurement can account for the apparent increase in thermal conductivity of a nanofluid suspension, because the experiments are interpreted using classical Fourier heat conduction theory. They suggest that a thermal relaxation theory like that of Cattaneo (1948) might play a major role since heat is then transported as a wave.

Recently, Christov (2009) has suggested a Lie derivative invariant form of time derivative in Cattaneo’s heat flux equation. Straughan (2010a,b) has studied thermal convection in a fluid and in a fluid-saturated porous medium using the theory which follows from the work of Cattaneo (1948) and Christov (2009), herein called Cattaneo–Christov theory. Cattaneo’s (1948) theory has been well studied in other areas of fluid and continuum mechanics (cf. Puri & Jordan 1999a,b; Christov & Jordan 2005, in press), although they did not employ the objective derivative of Christov (2009). The work of Straughan (2010a,b), shows that Cattaneo–Christov theory has the potential to yield interesting novel results in thermal convection.

The motivation for the present article is the work of Herrera & Falcón (1995) and of Falcón (2001). Herrera & Falcón (1995) discussed employing Cattaneo’s (1948) theory to study thermohaline convection and concluded that it could have important implications on stellar atmospheres, for example in a close binary system where helium-rich material may be transferred to a main sequence star, or in a neutron star, or in the core of a collapsed supernovae. Falcón (2001) discussed the cooling of white dwarfs and neutron stars by employing a Cattaneo theory. Herrera & Falcón (1995) derive physical balances and exploit Cattaneo’s (1948) law, but they do not produce a precise model for thermohaline convection employing Cattaneo’s (1948) heat flux law. The object of this work is to do precisely this. It is worth pointing out that in the stellar context, the phenomenon is referred to as thermohaline convection but the salt may not be sodium chloride. The idea is that a dissolved salt may render a part of the fluid heavier, thereby wanting to fall under gravity and generate a convective motion. For example, Herrera & Falcón (1995) write that, ‘… helium-rich material may be transferred to a main sequence star’ and ‘… a helium-rich outer layer is formed, leading thereby to the occurrence of a thermohaline instability at the interface between that layer and the original stellar material’. We derive a model for thermohaline convection in a porous medium saturated with a linear viscous fluid when the heat flux law is one of Cattaneo–Christov type. Given that the stellar material may be of a porous type, this is one area of potential application. We concentrate on the case where the saturated porous layer is heated and saturated below. This is a difficult mathematical case even when the classical Fourier heat law is employed (cf. Nield 1968; Nield & Bejan 2006; Straughan 2008, pp. 168–172). The problem is that the salt field stabilizes while the temperature field is destabilizing and this competition leads to a non-symmetric operator in the governing partial differential equations. For example, I do not think a nonlinear analysis of energy stability yielding a global (or even finite radius of attraction) threshold close to that of linear theory has been achieved (e.g. Straughan 2008, pp. 168–172). In the Cattaneo–Christov theory of thermohaline convection in a porous material, we find that the thermal relaxation effect plays a key role. As far as we know this is a very novel effect and we find that an instability threshold (tipping point) is established, which depends strongly on the Cattaneo number (relaxation time), but which is influenced by the porous properties and by the salt concentration.

2. Equations

To derive the model we present equations for fluid momentum, balance of mass, salt concentration, energy balance and the Cattaneo–Christov heat flux law. Standard indicial notation is adopted throughout.

For the fluid momentum, we adopt Darcy’s law (cf. Nield & Bejan 2006; Straughan 2008, pp. 10–14),Embedded Image 2.1 where Embedded Image are pressure, dynamic viscosity, permeability, density and gravity, k=(0,0,1) and vi is the pore averaged velocity. If Vi is the actual velocity in the fluid and ϕ is the porosity, then vi=ϕVi (cf. Nield & Bejan 2006; Straughan 2008, pp. 14–15). We adopt a Boussinesq approximation where density variations are present only in the buoyancy term and thenEmbedded Image 2.2 where T,C are the temperature and salt concentration, α is the thermal expansion coefficient of the fluid, αs is the salt expansion coefficient, T0,C0 are reference values and ρ0 is the density when T=T0, C=C0. With equation (2.2), Darcy’s law (2.1) becomesEmbedded Image 2.3 where Embedded Image is a modified pressure, z being the vertical coordinate. Since the fluid is incompressible,Embedded Image 2.4

The equation for the salt field is (cf. Straughan 2008, p. 15),Embedded Image 2.5 where kC is the salt diffusion coefficient.

We now write energy balances for the solid and fluid parts of the porous medium separately, and likewise write Cattaneo–Christov heat flux laws for each of the solid and fluid components (cf. Straughan 2008, pp. 14–15), which only deals with the classical Fourier theory. The solid and fluid parts are denoted by subscript s,f. Then we have the energy balance and Cattaneo law for the solid,Embedded Image 2.6 andEmbedded Image 2.7 where c is the specific heat, τs is the relaxation time, ks is thermal conductivity and Embedded Image the heat flux. For the fluid, the energy balance and Christov (2009) heat flux laws take the formEmbedded Image 2.8 andEmbedded Image 2.9 where cp is the specific heat at constant pressure, kf is the thermal conductivity and τf is the relaxation time. (The classical Cattaneo law is found from equation (2.9) by putting V0. If then τf=0, one recovers Fourier’s law Embedded Image. Equation (2.9) is an invariant form of Cattaneo law valid in a moving body and proposed by Christov (2009)).

To derive an averaged equation governing the (averaged) porous media properties, we multiply equation (2.6) by (1−ϕ), equation (2.8) by ϕ and add, and likewise multiply equation (2.7) by (1−ϕ) and equation (2.9) by ϕ and add. We now define the quantities Qi,(ρ0c)m,km,κ,M and τ byEmbedded Image Then, one shows that equations (2.6)–(2.9) lead toEmbedded Image 2.10 andEmbedded Image 2.11 The complete system of equations for our model is then (2.3), (2.4), (2.5), (2.10) and (2.11), which we denote by system Cattaneo–Christov–Thermal (CCT).

3. Conduction solution, perturbation equations

We now suppose the fluid saturated porous layer occupies the region Embedded Image. On the boundaries z=0,d, we have the boundary conditionsEmbedded Image 3.1 where TL,TU,CL,CU are constants with TL>TU,CL>CU, i.e. consistent with heating and salting below. Under these boundary conditions, system CCT admits the stationary solution whose stability we investigate, namely,Embedded Image 3.2 The steady pressure Embedded Image may then be found from equation (2.3).

To study the instability, we introduce perturbations (ui,π,θ,qi,φ) to the steady solution Embedded Image byEmbedded Image The nonlinear perturbation equations have the formEmbedded Image 3.3 where w=u3.

We now introduce non-dimensional scales,Embedded Image together with temperature, concentration, velocity U, pressure P and heat flux Q* scales asEmbedded Image andEmbedded Image We also need non-dimensional Cattaneo numbers, Embedded Image and the Lewis number, Le, given byEmbedded Image Important non-dimensional parameters are the Rayleigh number, R2, and salt Rayleigh number, R2s, defined byEmbedded Image Put ε=ϕM and then equations (3.3) in non-dimensional form (dropping stars) becomeEmbedded Image 3.4 These equations are to be solved together with the boundary conditionsEmbedded Image 3.5 together with the fact that the (x,y) behaviour of the solution satisfies a plane tiling periodic pattern.

4. Linear instability

To study the instability of the steady solution (3.2), we linearize equations (3.4) to obtainEmbedded Image 4.1 Now remove π by taking curlcurl of equation (4.1)1 and retain the third component. We take div of equation (4.1)5 and set Q=qi,i. Since the system is linear, we write w=eσtw(x), θ=eσtθ(x), φ=eσtφ(x), Q=eσtQ(x). (Actually, we write w,θ,φ,Q as Fourier series of terms like this, but since we are seeking an instability threshold, we can concentrate on one typical Fourier element.) We now define Embedded Image and ϵ=εLe and from equation (4.1) we must solveEmbedded Image 4.2 where Δ* is the horizontal Laplacian, Δ*=∂2/∂x2+∂2/∂y2. This is an eigenvalue problem for σ to be solved subject to the boundary conditionsEmbedded Image Owing to periodicity in (x,y), we introduce a plane tiling function f, and a wavenumber, a, with Δ*f=−a2f. (Details of the form of f for particular cell shapes are discussed in Chandrasekhar (1981, pp. 43–52), and in Straughan (2008, p. 152.) With the representationsEmbedded Image andEmbedded Image equations (4.2) reduce to solvingEmbedded Image 4.3 on z∈(0,1), where D=d/dz. Equations (4.3) are to be solved subject to the boundary conditionsEmbedded Image 4.4

We solve equations (4.3) and (4.4) exactly by analytical means, but also numerically using a D2-Chebyshev tau numerical method as a check (cf. Dongarra et al. 1996; Straughan 2010a). Since the instability surface is non-trivial, the numerical check proves very useful.

By writing W,Θ,Φ,Q as a series of terms of form Embedded Image (consistent with boundary conditions (4.4)), one may derive linear instability thresholds. The stationary convection boundary (σ=0) is given byEmbedded Image where Λ=n2π2+a2. Minimizing firstly over n, one finds n=1 and then a further minimization over a2 yieldsEmbedded Image 4.5 ac being the critical value of a which yields a minimum in R2. To find the oscillatory convection boundary, we follow the method of Chandrasekhar (1981) and put σ=1 in equation (4.3). After some calculation one may show this results in having to solve the equationEmbedded Image 4.6 where X=R2a2, with the coefficients A,B,C1 given byEmbedded Image Thus, the critical Rayleigh number for oscillatory convection is found by minimizingEmbedded Image 4.7 in a,n with Embedded Image. We find numerically the minimum is given with the minus sign and n=1.

We analyse the numerical results in §5. However, we note that Embedded Image. Also, taking the minus sign in equation (4.7) and seeking the asymptotic solution Embedded Image using the binomial theorem, we see that as Embedded Image, Embedded Image. This fact proves very useful in the next section.

5. Numerical results and conclusions

The surface of transition to instability one obtains is very interesting. This is a surface in Embedded Image space which we describe below. To understand this surface, it is constructive to first consider the two cases Embedded Image and R2s=0. When Embedded Image, one finds the stationary convection curve is given by equation (4.5), i.e.Embedded Image This yields the transition threshold for Embedded Image small. However, employing the Chandrasekhar (1981) method as in §4, we find the oscillatory convection curve is given by the straight lineEmbedded Image 5.1 (The coefficient A in equation (4.6) is zero when Embedded Image and equation (5.1) follows directly.) Thus, we have a situation where the instability curve is equation (4.5) until R2s reaches the value 4π2/(ϵ−1) and then equation (5.1) becomes the instability threshold. (For practical values of the parameters ϵ>1, e.g. if one employs values for salt as given in Straughan (2008, p. 254), then ϵ≈30.) On the other hand, when R2s=0, Straughan (2010b) shows that the instability threshold is R2=4π2 until Embedded Image, when the instability curve is Embedded Image, with the wavenumber changing from π to be infinite. Thus, the situation is as shown in figure 1.

Figure 1.

Instability curves for Embedded Image.

We now describe the instability surface in Embedded Image space. As we noted above, for practical values of a typical porous material saturated with brine, ϵ≈30. We have computed the instability surface for ϵ=10,100, and the behaviour is similar. However, ϵ=10 allows us to illustrate the effect better, the R2 values are larger and so we concentrate on this case.

When R2s is small, we find that the behaviour of the R2s=0 case continues. Thus, for ϵ=10, we find that the instability surface is given by R2=R2s+4π2 for Embedded Image, after which an increase in Embedded Image switches the surface to Embedded Image with the wavenumber changing from π to be infinite. Thus, we see stationary convection for Embedded Image up to the transition when oscillatory convection prevails but with Embedded Image. The Embedded Image value where we observe this behaviour numerically to change lies in the interval Embedded Image. Then, there is a range of Embedded Image, for ϵ=10, Embedded Image, where stationary convection is witnessed for Embedded Image small, then oscillatory convection holds with a finite wavenumber minimum until a second transition in Embedded Image is reached, when Embedded Image becomes the convection threshold. For ϵ=10, Embedded Image. For Embedded Image, we only witness oscillatory convection with finite wavenumber until a Embedded Image transition value is reached, whereupon Embedded Image is the instability value. Thus, the instability surface that emanates from the curves drawn in figure 1 is a complicated one in which we have three regimes in R2s. In one, we go from stationary convection to Embedded Image for Embedded Image; then in the second, for Embedded Image we have stationary convection → oscillatory convection →Embedded Image and then finally oscillatory convection →Embedded Image for Embedded Image. The transitions are observed as Embedded Image increases, for R2s fixed. For all values though the instability switches to Embedded Image once Embedded Image reaches a threshold. We detail this threshold in tables 1 and 2 for R2s taking values 0–20. In tables 1 and 2 the Embedded Image and aosc values refer to where the minimum given by equation (4.7) is achieved. Table 3 provides details of the various crossings of eigenvalues in different regimes; in particular, it presents some values of the transition from stationary convection to where oscillatory convection occurs with a finite wavenumber a. We note that the threshold represents for each R2s value a ‘tipping point’, where once Embedded Image is larger than the threshold, convection occurs once R2 is greater than Embedded Image. This threshold is dependent upon ϵ and also upon R2s. In particular, we note that as R2s increases, the critical value of Embedded Image decreases. Thus, increasing the basic concentration gradient of salt lowers the Embedded Image threshold beyond which convection occurs. With porous materials now being discovered which have a non-small relaxation time, effectively Embedded Image, this could have important consequences on thermohaline convection.

View this table:
Table 1.

Embedded Image values where transition to Embedded Image occurs. The transition occurs from stationary to oscillatory convection and a becomes infinite. Boldface in tables 13 indicates instability value.

View this table:
Table 2.

Embedded Image values where transition to Embedded Image occurs. NR indicates not relevant.

View this table:
Table 3.

Transition to oscillatory convection. For Embedded Image, a=π. The a value presented represents that value for oscillatory convection.


This work was supported in part by a grant from the Leverhulme Trust, ‘Tipping points: mathematics, metaphors and meanings’. I am indebted to three anonymous referees for their comments that have led to improvements in the manuscript.

  • Received February 21, 2010.
  • Accepted April 26, 2010.


View Abstract