Thermosolutal Convection with a Navier–Stokes–Voigt Fluid

We present a model for convection in a Navier–Stokes–Voigt fluid when the layer is heated from below and simultaneously salted from below, the thermosolutal convection problem. Instability thresholds are calculated for thermal convection with a dissolved salt field in a complex viscoelastic fluid of Navier–Stokes–Voigt type. The Kelvin–Voigt parameter is seen to play a very important role in acting as a stabilizing agent when the convection is of oscillatory type. The quantitative size of this effect is displayed. Nonlinear stability is also discussed, and it is briefly indicated how the global nonlinear stability limit may be increased, although there still remains a region of potential sub-critical instability, especially when the Kelvin–Voigt parameter increases.


Introduction
The Navier-Stokes equations for the flow of a linearly viscous, incompressible fluid are well known throughout the fields of applied mathematics and engineering. For a Navier-Stokes fluid the response of stress to change in the velocity gradient is instantaneous. This special response is not enjoyed by all real life fluids and, in particular, by many classes of viscoelastic and complex fluids. These fluids possess a stress which does not react instantaneously and instead they remember the velocity gradient history. Typically such fluids possess fading memory where the history dependence diminishes as the time advances into the past. Theoretical work on such fluids is vast, see e.g. Amendola and Fabrizio [1], Amendola et al. [2], Anand et al. [3], Anand and Christov [4], Christov and Christov [5], Fabrizio et al. [6], Franchi et al. [7][8][9], Gatti et al. [10], B Brian Straughan brian.straughan@durham.ac.uk 1 Department of Mathematical Sciences, University of Durham, Stockton Road, Durham DH1 3LE, UK Jordan and Saccomandi [11], Jordan et al. [12,13], Payne and Straughan [14], Yang et al. [15].
Recent work (arising from the Russian literature) is focussing on a very interesting class of complex (viscoelastic) materials associated with the names of Kelvin and of Voigt, cf. Chirita and Zampoli [16], Berselli and Bisconti [17], Layton and Rebholz [18]. Models for Kelvin-Voigt fluids have been presented by Oskolkov [19,20], and the references therein. Of particular relevance to the present article are the generalizations of the Kelvin-Voigt models to incorporate thermal effects as presented by Sukacheva and Matveeva [21], Kondyukov [22], Sukacheva and Kondyukov [23]. A special class of the thermal Kelvin-Voigt models are the so called Navier-Stokes-Voigt equations which are analysed here.
We are interested in the thermosolutal convection problem where the layer is heated from below and simultaneously salted from below resulting in a hard mathematical problem where the heat effect encourages the fluid to rise whereas salting below has the opposite effect and creates a competition between the physical attributes of heating and salting. This is the configuration of the promising method of generating renewable energy through electricity with a solar pond, see e.g. Abdullah et al. [35], and the present work may have a bearing on how to add an additive to stabilize the solar pond and ensure it functions in an optimal way.
The goal of this work is to develop a model and analyse thermosolutal convection in a layer of Navier-Stokes-Voigt fluid heated and salted from below, generalizing the models of Sukacheva and Matveeva [21] and Sukacheva and Kondyukov [23]. We believe this is the first analysis of this problem.

The Navier-Stokes-Voigt Equations for Thermosolutal Convection
Let v(x, t), T (x, t), C(x, t), p(x, t) be the velocity, temperature, concentration of a dissolved salt, and pressure at position x and time t of a body of fluid. The Navier-Stokes-Voigt equations for thermosolutal convection may be derived from Sukacheva and Matveeva [21], by incorporating the concentration effect, cf. Straughan [36, Sect. 14.1]. The appropriate system of equations may be written as Hereλ, ν, ρ 0 , α, γ, g, κ and κ s are the Kelvin-Voigt coefficient, the kinematic viscosity, a reference density, the coefficient of thermal expansion of the fluid, the coefficient of concentration in the density, gravity, thermal diffusivity, and salt diffusivity. The symbol denotes the Laplacian, k = (0, 0, 1), and standard indicial notation together with the Einstein summation convention is employed throughout. For example, , then we write the divergence of the velocity field in the forms For an example involving a nonlinearity, we write To derive Eq. (1) a Boussinesq approximation has been employed whereby the density in the buoyancy force is written as for a reference temperature T 0 and reference concentration C 0 . This gives rise to the αgT k i and γ gCk i terms. Details are similar to those in Sect. 14.1 of Straughan [36] where the Navier-Stokes analogue is analysed. As Oskolkov [19,20] notes, the Kelvin-Voigtλ term in (1) arises from a stress relation like where μ(= νρ 0 ) > 0 is a constant, d i j is the symmetric part of the velocity gradient, and σ i j is the Cauchy extra stress tensor, related to the stress tensor t i j by t i j = − pδ i j + σ i j . Then Eq. (1) 1 follows from the balance of linear momentum equation with f i being the body force term. Equation (2) is a linear relationship and as such is compatible with a linear instability analysis. One should perhaps consider an objective derivative in (2), as discussed in another context by Christov [37]. However, we here follow the procedure of Sukacheva and Matveeva [21], Sukacheva and Kondyukov [23]. Equations (1) are sometimes called the Navier-Stokes-Voigt equations, cf. Layton and Rebholz [18], or alternatively they can be called the Kelvin-Voigt equations of order zero, or the Oskolkov equations, see Oskolkov [20].
The boundary conditions on the planes z = 0 and z = d are given by The steady solution in whose stability we are interested is given bȳ where the temperature and concentration gradients, β, β s , are given by The steady pressurep is a quadratic in z which may then be derived from Eq. (1) 1 .
To investigate stability of the conduction solution (4) we introduce perturbations The equations for the perturbations are then derived and non -dimensionalized with the scalings, cf. details in Sect. 14.1 of Straughan [36], where d, T , U , C and T are the scales for length, time, velocity, concentration, and temperature, respectively. Define the Rayleigh number, Ra = R 2 , the salt Rayleigh number Rs = C 2 , the Lewis number Le, and the Prandtl number Pr, by The Kelvin-Voigt parameterλ is rescaled as Now drop the *s and treat x i and t as the non-dimensional variables. The nondimensional form of the perturbation equations to (1) is where w = u 3 .
Equations (5) are defined on the domain R 2 × {z ∈ (0, 1)} × {t > 0} and the boundary conditions are together with the fact that u i , φ, θ and π satsify a plane tiling periodicity in the x, y plane.
To analyse linear instability we linearize (5) and seek solutions of the form Next take curl curl of Eq. (7) 1 and retain the third component of the resulting equation to reduce Eq. (7) to where D = d/dz and Eq. (9) hold on z ∈ (0, 1).
We must specify further information on boundary conditions. We already have We believe this is the first analysis to address the thermosolutal convection problem for a Navier-Stokes-Voigt fluid. Thus we restrict attention to two stress free surfaces and so we additionally require The equations lead to even derivatives being zero on the surfaces z = 0, 1, and then we may write w, φ and θ as a sin series in z of form sin(nπ z), n = 1, 2, . . . We take n = 1 since computations show this yields the lowest value of the Rayleigh number.
With the sin representation of solutions, Eq. (9) then yield a determinant equation which leads to the following expression for R 2 , where = π 2 + a 2 .
The stationary convection threshold follows from (12) and takes the form Upon minimization in a 2 we find the critical wave number value is a 2 c = π 2 /2 and inserting this in (13) it follows that Further analysis requires us to take the real and imaginary parts of (12) and recollect that R 2 is real. To find the oscillatory convection curve we put σ = iω, ω ∈ R, and then the real part of (12) yields From the imaginary part of (12) we may show that The critical value of the oscillatory convection threshold R 2 is found by employing (16) in (15) and then minimising Ra = R 2 in a 2 for fixed values of λ, Le and Pr. This we do numerically and report the findings in Sect. 4.

Remark 1
It is possible to add a Rayleigh friction term to the right hand side of equation (5) 1 as do Di Plinio et al. [39] for the isothermal and non concentration case. In this situation the perturbation equations (5) are replaced by for a constant ξ > 0. By proceeding as outlined above one may perform a stationary convection and an oscillatory convection analysis for system (17).

Remark 2
One could also consider a Soret effect upon the convection threshold for a Navier-Stokes-Voigt fluid, cf. Straughan and Hutter [34]. In this case the perturbation equations (5) are replaced by for a constant Soret coefficient S > 0. Again, one may proceed as outlined above to perform a stationary convection and an oscillatory convection analysis for system (18).

Numerical Results
We report on numerical solutions to Eqs. (14), (15) and (16). To do this we require values for the Prandtl and Lewis numbers. The Prandtl number for water at 20 • C is 6.99 and for vegetable oils Rodenbush et al. [40] report values in the range 25-400. For the thermal and solute diffusivities we employ values from Caldwell [41], Ozbek and Phillips [42], Yanez Limon et al. [43] and We have computed values of the Rayleigh number and critical wave number versus the salt Rayleigh number for various combinations of Pr, Le and λ. For the many values we examined the qualitative behaviour is always the same and resembles that of Fig. 1, where Pr = 50, Le = 21.67, although the quantitative values at the transitions depend on the parameters Pr, Le and λ. We here report on values for Pr = 50, Le = 21.67, for Pr = 6.99, Le = 70.44, and for Pr = 6.99, Le = 111.2. The pattern of instability in all cases follows that of Fig. 1 in that the stationary convection curve begins at Ra = 27π 4 /4 for Rs = 0 and is then the straight line E-S as Rs increases. Depending on the value of the Kelvin-Voigt parameter λ there is a transition from stationary to oscillatory convection as shown in Fig. 1, where the 0-0 curve denotes the oscillatory convection threshold when λ = 0, 1-1 curve denotes the oscillatory convection threshold when λ = 1, and the 2-2 curve denotes the oscillatory convection threshold when λ = 2. The basic conduction solution (4) is unstable when the (Ra, Rs) values lie above the linear instability curve, and convective motion ensues. For example, when λ = 1, the solution is unstable above the straight line E-1, but for Rs = C 2 values greater than the transition the instability is oscillatory when Ra = R 2 lies above the curve 1-1.
Observe that λ increasing raises the instability threshold, and as we might expect λ is a stabilizing influence. Table 1 yields the transition values for the values of Pr and Le as indicated. The a 2 values reported in this table are the critical values for oscillatory convection. We observe that the critical oscillatory wave number values are smaller as λ increases indicating that the convection cells become wider for increasing λ. The Lewis number effect is to lower the critical Rayleigh number as Le increases.

Nonlinear Stability Considerations for Thermal Convection in a Navier-Stokes-Voigt fluid
Let (·, ·) and · denote the inner product and norm on the Hilbert space L 2 (V ), where V is a period cell for the solution.
To develop a nonlinear stability analysis a classic approach to nonlinear energy stability will multiply (5) 1 by u i , (5) 2 by θ and (5) 3 by φ, and integrate each in turn over the period cell V using the boundary conditions. In this manner we obtain the identities d dt and d dt together with d dt We may add these equations to derive the energy equation where the energy function is given by the production term I is defined by and the dissipation is From this one computes the value where H is the space of admissible solutions, cf. Straughan [36]. Since the C terms disappear from I the Euler-Lagrange equations which arise from (26) yield the global nonlinear stability threshold of R E = 27π 4 /4, i.e. the straight line E-E in Fig. 1. This means that values of Ra, Rs below this threshold yield a globally stable nonlinear solution. However, there is still the possibility of sub-critical convection in the region above E-E and below the linear instability threshold. If one is prepared to employ a generalized energy functional by introducing ψ = θ − δφ where δ = C/R, and work with the function then see Joseph [44], Mulone [30], it is possible with a lot of analysis and much effort, employing numerical optimization on the parameters ω 1 > 0, ω 2 > 0, while also minimizing in the wavenumber, see Straughan [45], to increase the global nonlinear stability boundary into the region in Fig. 1 denoted by E-0-0-E. However, a lot of analytical and computational effort is required. The threshold so obtained does not coincide with the linear instability one and so there remain regions of possible subcritical instability. However, for the Navier-Stokes-Voigt fluid the coupling parameter generalized energy method must involve λ in the maximization, Euler-Lagrange process, to be able to achieve a global nonlinear stability boundary which will exceed the 0-0 curve when λ > 0. We believe at present this problem is open.
In the present situation because the energy in (23) involves λ ∇u 2 , the rate of decay obtained is better than that observed in classical Bénard convection. This faster decay

Conclusions
We have developed an analysis for the thermosolutal convection problem for a Navier-Stokes-Voigt viscoelastic fluid. We concentrate on the problem where the layer of fluid is heated from below and simultaneously salted from below. The thermal convection problem for a Navier-Stokes-Voigt fluid without a salt field was suggested by Sukacheva and Matveeva [21], Sukacheva and Kondyukov [23]. The isothermal model for a Navier-Stokes-Voigt fluid is explained at length by Oskolkov [19,20]. We have shown that for small enough salt Rayleigh number the onset of thermally driven motion is by stationary convection. There is then a transition to oscillatory convection. The transition values are calculated for several realistic values of Prandtl and Lewis numbers. We have shown that the transition to oscillatory convection occurs at higher Rayleigh and salt Rayleigh numbers as the Kelvin -Voigt parameter λ increases. The fact that the convective motion instability threshold is increased when the Kelvin-Voigt parameter increases is very useful for solar pond design. A solar pond is a mechanism whereby solar energy is harnessed through a thermosolutal process and converted into renewable electric energy, cf. Abdullah et al. [35]. A solar pond consists of a layer of salt water approximately 1-2 m thick in a horizontal configuration with direct access to solar radiation. The salt field is arranged so that the salt concentration decreases approximately linearly from the base of the layer. The base of the layer is chosen to absorb solar radiation and this configuration can achieve temperatures close to 100 • C in the solution near the base. The naturally heated brine solution is drawn out and passed through a heat exchanger to generate renewable electricity. It is important that convective motion does not commence otherwise the solution mixes and so it is key that the conditions remain in the stable case for the problem studied herein. Our work suggest that adding a suitable additive to the salt solution which increases the Kelvin -Voigt parameter will ensure the solar pond does not commence overturning instability and will, therefore greatly improve the efficiency of the device.