Tipping points in Cattaneo–Christov thermohaline convection

B. Straughan

Abstract

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), 2.1 where 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 then 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) becomes 2.3 where is a modified pressure, z being the vertical coordinate. Since the fluid is incompressible, 2.4

The equation for the salt field is (cf. Straughan 2008, p. 15), 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, 2.6 and 2.7 where c is the specific heat, τs is the relaxation time, ks is thermal conductivity and the heat flux. For the fluid, the energy balance and Christov (2009) heat flux laws take the form 2.8 and 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 . 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 τ by Then, one shows that equations (2.6)–(2.9) lead to 2.10 and 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 . On the boundaries z=0,d, we have the boundary conditions 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, 3.2 The steady pressure may then be found from equation (2.3).

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

We now introduce non-dimensional scales, together with temperature, concentration, velocity U, pressure P and heat flux Q* scales as and We also need non-dimensional Cattaneo numbers, and the Lewis number, Le, given by Important non-dimensional parameters are the Rayleigh number, R2, and salt Rayleigh number, R2s, defined by Put ε=ϕM and then equations (3.3) in non-dimensional form (dropping stars) become 3.4 These equations are to be solved together with the boundary conditions 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 obtain 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 and ϵ=εLe and from equation (4.1) we must solve 4.2 where Δ* is the horizontal Laplacian, Δ*=∂2/∂x2+∂2/∂y2. This is an eigenvalue problem for σ to be solved subject to the boundary conditions 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 representations and equations (4.2) reduce to solving 4.3 on z∈(0,1), where D=d/dz. Equations (4.3) are to be solved subject to the boundary conditions 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 (consistent with boundary conditions (4.4)), one may derive linear instability thresholds. The stationary convection boundary (σ=0) is given by where Λ=n2π2+a2. Minimizing firstly over n, one finds n=1 and then a further minimization over a2 yields 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 equation 4.6 where X=R2a2, with the coefficients A,B,C1 given by Thus, the critical Rayleigh number for oscillatory convection is found by minimizing 4.7 in a,n with . 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 . Also, taking the minus sign in equation (4.7) and seeking the asymptotic solution using the binomial theorem, we see that as , . 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 space which we describe below. To understand this surface, it is constructive to first consider the two cases and R2s=0. When , one finds the stationary convection curve is given by equation (4.5), i.e. This yields the transition threshold for small. However, employing the Chandrasekhar (1981) method as in §4, we find the oscillatory convection curve is given by the straight line 5.1 (The coefficient A in equation (4.6) is zero when 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 , when the instability curve is , with the wavenumber changing from π to be infinite. Thus, the situation is as shown in figure 1.

Figure 1.

Instability curves for .

We now describe the instability surface in 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 , after which an increase in switches the surface to with the wavenumber changing from π to be infinite. Thus, we see stationary convection for up to the transition when oscillatory convection prevails but with . The value where we observe this behaviour numerically to change lies in the interval . Then, there is a range of , for ϵ=10, , where stationary convection is witnessed for small, then oscillatory convection holds with a finite wavenumber minimum until a second transition in is reached, when becomes the convection threshold. For ϵ=10, . For , we only witness oscillatory convection with finite wavenumber until a transition value is reached, whereupon 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 for ; then in the second, for we have stationary convection → oscillatory convection → and then finally oscillatory convection → for . The transitions are observed as increases, for R2s fixed. For all values though the instability switches to once reaches a threshold. We detail this threshold in tables 1 and 2 for R2s taking values 0–20. In tables 1 and 2 the 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 is larger than the threshold, convection occurs once R2 is greater than . This threshold is dependent upon ϵ and also upon R2s. In particular, we note that as R2s increases, the critical value of decreases. Thus, increasing the basic concentration gradient of salt lowers the threshold beyond which convection occurs. With porous materials now being discovered which have a non-small relaxation time, effectively , this could have important consequences on thermohaline convection.

View this table:
Table 1.

values where transition to 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.

values where transition to occurs. NR indicates not relevant.

View this table:
Table 3.

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

Acknowledgements

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.