Effect of Temperature Upon Double Diffusive Instability in Navier–Stokes–Voigt Models with Kazhikhov–Smagulov and Korteweg Terms

We present models for convection in a mixture of viscous fluids when the layer is heated from below and simultaneously the pointwise volume concentration of one of the fluids is heavier below. This configuration produces a problem of competitive double diffusion since heating from below promotes instability, but the greater density of fluid below is stabilizing. The fluids are of linear viscous type which may contain Kelvin–Voigt terms, but density gradients due to the mixture appear strongly in the governing equations. The density gradients give rise to Korteweg stresses, but may also be described by theory due to Kazhikhov and Smagulov. The systems of equations which appear are thus highly nonlinear. The instability surface threshold is calculated and this is found to have a complex nonlinear shape, very different from the linear ones found in classical thermohaline convection in a Navier–Stokes fluid. It is shown that the Kazhikhov–Smagulov terms, Korteweg terms and Kelvin–Voigt term play a key role in acting as stabilizing agents but the associated effect is very nonlinear. Quantitative values of the instability surface are displayed showing the effect Korteweg terms, Kazhikhov–Smagulov terms, and the Kelvin Voigt term have. The nonlinear stability problem is addressed by means of a generalized energy theory deriving different results depending on which underlying theory is employed.

studies is increasingly diverse. Mixtures of fluids are highly complex and while they may be described by a stess tensor which depends linearly on the velocity gradient it is known that there has to be additionally a dependence on density gradients due to the different constituents. Such complex fluids define a class which is often referred to as generalized Navier-Stokes fluids, and these generalized fluids also encompass those which include viscoelastic effects where the stress tensor depends on the history of the velocity gradient, and this area is vast, as is witnessed by the work of [1-3, 16, 17, 23, 33, 34, 75] In this article we study a mixture of miscible vicous fluids but we also allow the mixture to be comprised of a particular class of viscoelastic fluids associated with the names of Kelvin and of Voigt, cf. [3,5,57]. Analysis studies of Kelvin-Voigt fluids (which encompass Navier-Stokes-Voigt) have been presented in much detail by [49,50], and generalizations of these to incorporate temperature effects are given by [42,48,68] and [67], see also [40,47,53,69]. Kelvin-Voigt fluids of order zero are also known as Navier-Stokes-Voigt fluids and the equations for these have been analysed in great detail with respect to questions of existence, behaviour of attractors by e.g. [10,20,36,37]. In connection with classes of viscoelastic fluids we point out that a very useful account of Maxwell, Oldroyd and Kelvin-Voigt fluids of various orders is given by [51] who discuss the solution existence question at length. Double diffusion is caused by two effects such as a temperature gradient and a concentration gradient and these may often be in opposition. For example if a layer of water (or other solvent) has sodium chloride dissolved in it then instability may be caused by the salt concentration being heavier at the top or by heating the base of the layer and the resulting convective motion is witnessed in the fluid (Newtonian or viscoelastic). This phenomenon is widely studied in the literature, see e.g. [7,13,14,31,52,[62][63][64][65][66]73][pp. 238-268]. However, double or multi diffusive instability may also be witnessed in a mixture of miscible liquids subject to a temperature gradient, but where the density gradient due to the different fluids is also an ingredient for instability. Indeed, double or multi diffusive instability for a mixture of miscible fluids is an area which is increasingly occupying application areas. For example, diffusion in glycerol-water mixtures has many applications, see e.g. [12,21,74]. Double or multi diffusion convective instability is particularly important in the area of renewable energy. A solar pond is a device which employs the Sun's rays to produce electricity and recent research in this area is concentrating on using different solvents and dissolved salts, being a precise application of multi fluid mixture convection, see e.g. [4,59].
In a mixture of viscous fluids one has the usual dependence of the stress tensor upon the symmetric part of the velocity gradient, but there will in general also be a dependence upon gradients of density. The idea of incorporating density gradients into the stress is attributed to Korteweg, see e.g. [70, p. 514]. Further information on Korteweg stress formulations be may found in [41,46], and these writers analyse in detail thermodynamic aspects and different flow regimes.
Apart from the applications mentioned above, Korteweg stresses are important in real life and have application in vulcanology, see [71], and in other areas of geophysics, see [45]. Beginning with equations of [26,35] developed a model for a layer of fluid mixture where the stress incorporated Korteweg terms, and they particularly studied the stability of an isothermal liquid layer where convective motion may be driven by gravity if the upper part of the layer is denser than that below. [26] derived a set of Korteweg-Boussinesq equations which are highly useful in this process. The question of uniqueness of a solution to the isothermal Korteweg-Boussinesq equations is addressed by [27]. Temperature effects were incorporated to stability problems involving incompressible viscous fluids with Korteweg stresses by [54]. These writers investigate instability arising due to density gradients via Korteweg stresses, but they do not incorporate buoyancy effects as in [26,27]. The Korteweg effect will dominate in a zero gravity atmosphere, but even in low gravity conditions the joint effect of buoyancy and density gradients should be important. One of the aims of this work is to include the effect of gravity on instability in addition to the density gradient contribution.
The effect of Korteweg terms upon fluid mechanics is increasingly becoming popular from a continuum mechanics point of view, see e.g. [18,19,28], but particularly from a perspective of analysis in partial differential equations, see e.g. [11,15,32,58,60,72].
A separate development to deriving a theory for a mixture of isothermal fluids which incorporates density gradient terms is due to [38,39], and this theory was employed to analyse gravitational instability when the fluid layer is heavier at the top by [24]. The [38,39], theory has been investigated in some depth from an anlysis viewpoint of existence and solution behaviour by [55,56], and also by [8,9]. More recent papers on analysis of the Kazhikhov-Smagulov models is due to [29,30].
In this work we firstly extend the Korteweg-Boussinesq model of [26] to include temperature effects and thus be applicable to the double diffusive convection scenario. We are particularly interested in the competitive case where temperature is destabilizing whereas the concentration is stabilizing. However, we also develop the [38,39] model to include temperature effects. To do this we employ a suitable Boussinesq approximation, see [6], and derive the model from the completely nonlinear set of equations analysed in the isothermal case by [8]. This development is very interesting because we show that in addition to the Navier-Stokes terms present in standard double diffusion theory there arise two new sets of terms in the momentum equation. One set of terms we attribute to the Kazhikhov-Smagulov development. The second set of terms has the same form as those which arise from Korteweg theory. Hence, the Kazhikhov-Smagulov full theory analysed by [8] contains the Korteweg terms naturally. Thus, we effectively show that there are three models which arise. One corresponds to the Korteweg theory, the second corresponds to Kazhikhov-Smagulov theory, and the third is a combination of both.
We perform a linear instability analysis and a complementary nonlinear stability analysis for all three models for a mixture of fluids in the competitive double diffusion situation. For the Korteweg model the nonlinear stability obtained is global, in the sense that the stability is for all initial data. For the model derived by [8] from the full set of Kazhikhov-Smagulov equations the nonlinear stability analysis obtained is conditional upon the size of the initial data.

Equations
Let x denote a spatial point in a three-dimensional body and let t denote time. If ρ(x, t) denotes the density of a mixture of two miscible fluids then the Cauchy stress tensor, t i j , for a linearly viscous fluid with Korteweg stress terms may be written as (1) see [41, 70, p. 514]. Here d i j is the symmetric part of the velocity gradient, i.e.
where v i is the velocity field. We employ standard indicial notation throughout together with the Einstein summation convention. For example, the divergence of the velocity field is for a function C depending upon x, t. In (1) P is the pressure in the fluid and we here treat the coefficients, β 1 , β 2 , γ, H 1 , H 2 as constants. For an incompressible fluid the γ term is omitted from (1) and since the stress tensor appears in the momentum equation as t i j, j the term in H 2 can also be incorporated into a modified pressure of form The momentum equation for a fluid in the Korteweg-Boussinesq equations as derived by [26], employs a Boussinesq approximation, see [6]. The momentum equation of [26] in our notation may be written where ν(C) = μ(C)/ρ 0 is the kinematic viscosity of the fluid, ρ 0 is a constant reference value of density, g is gravity, k = (0, 0, 1), and we have written K 1 = −H 1 .
[26] mainly discuss the case H 1 < 0 but they leave the possibility that it could be positive. We choose to retain H 1 negative since experimental values, cf. [54], suggest this, and also this sign of coefficient arises naturally in the development from the full set of Kazhikhov-Smagulov equations. The function C is the concentration of one of the fluids (the volume fraction) and we observe that [35] defines this quantity for a mixture of water and glycerin to be C = V W /V T , where V W is the volume of water and V T is the total volume. (This must be interpreted in the continuum mixture sense where both constituents exist at each point and one considers a representative elementary volume around a point x containing both water and glycerin and V W is the limit as the volume tends to zero. Of course, V W in this way changes with time, in general.) In [26] the density is written as where C 0 is a constant reference value, and α c is the expansion coefficient associated with the concentration variation. Upon employing the density equation (3) in (2) we arrive at the complete momentum equation in the Korteweg-Boussinesq approximation of [26], namely where the constant term in the density, (3), has also been absorbed into the pressure. The fluid is (in a sense) incompressible, see [26], and then the continuity equation is The evolution equation for the concentration is The Korteweg-Boussinesq system of equations of [26], Eqs. (4.11)-(4.13), is comprised of (4), (5), (6), and for clarity we rewrite them together here, To extend the model of [26] to incorporate temperature we need to modify the density in (3) to allow for variable temperature and thus the density is now written as where ρ 0 , T 0 , C 0 are constant reference values, T (x, t) is the temperature and α is the expansion coefficient of the fluid. The evolution equation for the temperature field may now be written as cf. [54], where κ is the thermal diffusivity which is assumed constant.
There is a case for the viscosity to depend on both temperature and concentration, but as this is the first time we have seen work on the double diffusive convection problem we choose to keep the kinematic viscosity constant and so proceeding in this manner we may arrive at the full system of Korteweg-Boussinesq equations for double diffusive convection in the form The momentum equation of the [38,39] equations is given by [8] as where λ > 0 is a constant we call the Kazhikhov-Smagulov parameter. In Eq. (11), f i is the body force and P is a modified pressure given by wherep is the actual fluid pressure andλ is another constant, cf. [24], Eqs. (2.5), (2.6). [8] writes that Kazhikhov and Smagulov consider the simplified equation which follows from (11) by omitting the term involving λ 2 . This simplified system leads to what [24] call the Kazhikhov-Smagulov equations. [8] also notes that Kazhikhov and Smagulov require the condition where M and m are the supremum and infimum of the initial density of the mixture.
[8] removes the condition (12) and establishes existence of a solution in a precise sense working with the full equation (11). We divide equation (11) by ρ > 0 and we observe that the resulting λ 2 terms may be rewritten as Now defined a modified pressure p by and then we may write equation (11) as We next employ a Boussinesq approximation on (13) and in a sense we follow the development of [26] where they proceed from their Eqs. (4.3)-(4.7) to arrive at their Korteweg-Boussinesq Eqs. (4.11)-(4.13). We assume that ρ is a constant, ρ 0 , when it appears undifferentiated in that equation but we retain the spatial derivative terms. The density in the body force component is assumed given by and then we replace the derivatives of ρ in terms of C (the volume concentration) and appending equations (9), (6) for T and C we arrive at a Boussinesq-approximation set of equations arising from the full set of nonlinear equations addressed by [8].
The full system of equations we derive has form It should be observed that we have added a Kelvin-Voigt regularizing term, the term involvingδ, to the left hand side, [20]. Whenδ = 0 we obtain a Navier-Stokes type theory whereas whenδ > 0 there results a Navier-Stokes-Voigt type theory. We also note that if we omit the λ/ρ 0 term then this becomes analogous to the Korteweg theory discussed earlier whereas when we omit the λ 2 /ρ 2 0 term we have what might be termed a Kazhikhov-Smagulov type theory.

Double Diffusive Convection
Suppose Eq. (14) hold in the horizontal layer R 2 × {z ∈ (0, d)} with the velocity, temperature and concentration prescibed on the boundaries z = 0 and z = d, as for prescribed constant values T L , T U , C L , C U , with T L > T U > 0 and C L > C U where T L , T U are in • K. Thus temperature has a tendency to destabilize whereas concentration stabilizes.
The steady solution to (14) and (15) in whose stability we are interested is where the temperature and concentration gradients, β, β c , are given by The steady pressurep may then be derived up to a constant at ones disposal from (14) 1 .
To investigate stability of the steady solution (16) we introduce perturbations wherep is the steady state pressure arising from (14) 1 .
The equations for the perturbations are then derived from (14) and these equations are non-dimensionalized with the scalings, where Le, Sc and Pr are the Lewis, Schmidt and Prandtl numbers. The Rayleigh number, R, and concentration Rayleigh number, S, are introduced as and we introduce two further non-dimensional parameters, namely the Kazhikhov-Samgulov parameter, α 2 , and the Korteweg parameter, α 1 , by The fully nonlinear non-dimensional equations for (u i , θ, φ, π) then become (omitting the *s) where w = u 3 , and these equations hold on the domain The boundary conditions to be satisfied are together with the fact that u i , φ, θ and π satisfy a plane tiling periodicity in the x, y plane. We believe this is the first analysis of this model and so while α 1 and α 2 are connected we here treat them as independent constants in order to highlight the Kazhikhov-Smagulov and Korteweg effects on stability and instability. Let us observe that if α 2 = 0, α 1 = 0, then (19) represents a direct extension of the Korteweg-Boussinesq model of [26], to the double diffusive convection problem. If α 1 = 0, α 2 = 0, then (19) may be thought of as a Kazhikhov-Smagulov equation model for double diffusive convection. When α 1 = 0, α 2 = 0, equations (19) are a double diffusive convection model arising from the full set of Kazhikhov-Smagulov equations analysed by [8]. We point out that if we take δ = 0 then (19) is representative of Navier-Stokes theory whereas when the Kelvin-Voigt parameter δ > 0, (19) is representative of Navier-Stokes-Voigt theory.
We seek now solutions like (x, y), where h is a function compatible with tiling the plane, and h satisfies * h = −a 2 h, a being a wavenumber. To illustrate the novel effects of the Korteweg and the Kazhikhov-Smagulov terms, in Sect. 6 we concentrate on two stress free surfaces, cf. [22,31], and then the boundary conditions required are where D = d/dz. Care has to be taken with deriving the boundary conditions. When one solves for two fixed surfaces the condition D 2 W = 0 on z = 0, 1, is replaced by DW = 0 on z = 0, 1. In deriving the appropriate condition for two stress free surfaces we recall that the Cauchy stress tensor has the form and this yields the following conditions on the stress tensor on the planes z = 0, 1, where α = 1, 2. For the linearized case the K 1 term contribution leads to −K 1 φ ,αC,z − K 1C,α φ ,z , on z = 0, 1, and since φ = 0 on z = 0, 1, andC =C(z), these terms vanish, which leads to the condition on D 2 W .
We firstly study the problem with only Korteweg terms and no Kazhikhov-Smagulov contributions, so that α 2 = 0, α 1 = 0. In this case the eigenvalue problem for instability reduces to solving the determinant equation where = π 2 + a 2 .
One sees that the stationary convection boundary, σ = 0, is given by and the critical value for instability is obtained by minimizing R stat in a 2 .
To determine the oscillatory convection critical Rayleigh number one solves (22). This yields a cubic in σ which is then resolved into real and imaginary parts. After some calculation one may then show that the oscillatory Rayleigh number is given by The critical value for oscillatory convection is found by minimizing R osc in a 2 .
Numerical results for stationary and oscillatory convection are displayed in Sect.

6.
To determine the critical Rayleigh numbers in general (including for the cases α 1 = 0, α 2 = 0 or α 1 = 0, α 2 = 0), requires numerical solution of (21). This we do by employing a Chebyshev-tau method coupled with the QZ algorithm. In this case (21) are written in the form For fixed R, S, Sc, Le, δ, α 1 , α 2 , the eigenvalues σ are computed and the secant method is employed to locate where σ r = 0 (σ = σ r + iσ 1 ), and then this value of R is minimized in a to yield the critical values of R.

Nonlinear Stability
We now analyse nonlinear stability of the solution (16) for the fully nonlinear equations (19), (20). To facilitate the nonlinear analysis we introduce the variable ψ by the transformation ψ = Rθ, where R 2 = R. This changes the Rθ k i term in (19) 1 to Rψk i and it changes the w term in (19) 3 to Rw. Let V denote the period cell for the perturbation solution and denote by · and (·, ·) denote the norm and inner product on L 2 (V ). We have observed there are effectively three systems of equations of partial differential equations in (19), (20).
Define R E F I to be the solution to the maximization problem for the standard Bénard problem 1 where and H in (26) , with the constants in c being defined in the proof.

Proof.
Multiply (19) 1 by u i (with Rθ ≡ Rψ) and integrate over V . After some integration by parts and use of the boundary conditions one finds d dt and d dt We require a further identity and so we follow an idea in [26,27]. Hence, multiply (19) 4 by − φ and integrate over V to find d dt From (27)- (30) we may now arrive at a generalized energy equation of form where and I = 2R(ψ, w), For part A of the theorem, α 2 = 0, and then N = 0, and we set Since I only involves w and θ , R E may be replaced by R E F I and using Poincaré's inequality where For R < R E F I (≡ R E ) then put k = 1 − R/R E F I > 0, and using (32) and (33) we may arrive at and so E(t) ≤ E(0) exp(−kc 2 t).
Thus we have global nonlinear stability in the measure E(t), and part A is proved.
To prove part B we note that using integration by parts where · 4 denotes the norm on L 4 (V ) and c 1 is the constant in the Sobolev embedding of where the constant c 3 may be found in [25, pp. 280-283]. From the forms for E(t) and D(t) one may show and so employ this inequality together with the bounds (34) and (35) to obtain where Next, employ (36) in (31) together with the procedure leading to (32), (33) to find If now E 1/2 (0) < k/c 4 then one may use a continuity argument to show E(t) → 0, exponentially, see e.g. [62, pp. 14-16], and part B follows.

Remark.
When α 2 = 0, α 1 = 0 the above proof does not work. [24] also found difficulty numerically to proceed with the reduced equations of Kazhikhov-Smagulov in their analysis of convective overturning in the isothermal problem.

Numerical Results
To understand the novel effects of this work we firstly discuss the instability diagram for classical thermohaline convection in a linearly viscous fluid, see Fig. 1. This is when α 1 = 0, α 2 = 0, δ = 0. Note that for S less than the threshold value S = 51.9 the instability curve is a straight line of slope one which is intersected by another straight line which is the oscillatory convection threshold. Once S is larger than the transition value oscillatory convection is the mechanism leading to instability.
In Fig. 2 we show a typical analogous picture when Korteweg terms are present, i.e. α 2 = 0, δ = 0, but α 1 = 0. Observe that the Korteweg terms have a very strong effect on the instability curves. The curve for stationary convection is very nonlinear as S increases. Once again there is a transition to oscillatory convection but at a much smaller value of S = 5.6, indicating the strong stabilizing effect of the Korteweg terms. The oscillatory convection curve is actually also very nonlinear. For example, using the parameter values in Fig. 2, R increases from 709.424 at S = 0 to 713.043 when S = 10. However, when S = 100, R has value 1193.81. Figure 3 shows the stationary and oscillatory curves for the same parameter values as Fig. 2, but now the effect of the Kelvin-Voigt parameter, δ, is seen. This parameter only affects oscillatory convection, but it is seen that δ has a stabilizing effect. We have only presented details for Le = 20 and Sc = 40, but similar behaviour is observed for all the other values of Le and Sc we have employed. A referee has pointed out that oscillatory convection only occurs when Le > 1. This is true, although for most fluids encountered in real life this condition holds.
The effect of the Kazhikhov-Smagulov parameter, α 2 , upon the stationary convection threshold is observed in Tables 1 and 2. The parameter α 2 , like α 1 , has a strong stabilizing effect. For the scaling chosen here the effect of α 1 is greater than that of α 2 , as seen in Table 1, but both parameters have a pronounced effect on the threshold of convective motion. We also observe that as α 1 increases a 2 displays a strong decrease whereas when α 2 increases, a 2 increases. This means that as α 1 increases the aspect ratio of the convection cells at the onset of convection is increasing and so the cells become wider with increasing α 1 . Thus increasing internal stresses due to density gradients are leading to wider cells. The Kazhikhov-Smagulov parameter has the opposite effect. As α 2 increases the cells become narrower. Hence, while both effects stabilize the onset of stationary convection, each term has a very different effect on the behaviour of the cell shape. We should point out that the eigenvalue problem to calculate σ is a difficult one numerically, and spurious eigenvalues are definitely witnessed.

Conclusions
We have produced a model to describe the behaviour of a mixture of two miscible fluids where differences in concentration by volume of a given fluid lead to density gradients in addition to the expansion by heating. The double diffusive convection problem is analysed in detail by a linear instability technique, but also by a fully nonlinear energy argument. There are two sources of terms which are not present in classical thermosolutal convection, these being associated with Korteweg effects and also with effects due to work of Kazhikhov and Smagulov. Both effects are shown to have a strong influence on the convection threshold and also upon the size of the resulting convection cells.
Froma mathematical viewpoint the novelty of the Korteweg and of the Kazhikhov-Smagulov terms may be seen by examining the linearized theory. Omitting the nonlinear terms in (19) the equations governing linearized instability theory are It is convenient to cast equations (38) together with the boundary conditions (20) in the setting of an abstract system, defined in a Hilbert space, cf. [61, p. 156]. For our purpose A and L are symmetric, unbounded linear operators, M is a skew-symmetric linear operator, and N (u) is a nonlinear operator. We operate on (38) 4 by S − α 1 S 2 , and derive the following equation Then (38) and (40) may be written in the form (39) where u = (u, v, w, ψ, φ) T and essentially L and M have the forms [61, p. 156] points out that when = 0 (i.e. α 1 = α 2 = 0, or S = 0) then, under suitable conditions on N , the linear and nonlinear stability boundaries coincide. If M is bounded then σ will be real and the instability and nonlinear stability boundaries remain within O( ) of each other if is small enough. However, he also notes that  The second and third columns refer to the case where α 1 = 0, whereas columns five and six refer to where α 2 = 0. The parameter S = 10 there are many problems in applied mathematics in which is not small, citing the rotating Bénard problem and the magnetohydrodynamic Bénard problem as examples. The current work is also a novel example of this. It is typical that non-zero skew symmetric terms in (39) lead to stabilizing effects and the α 1 and α 2 terms are each involved with skew symmetric operators. We also draw attention to the fact that [43,44] have identified interesting attractors and behaviour for ordinary differential equation systems derived from double diffusive convection theory. We believe the new effects identified here could reveal yet further interesting attractor behaviour.
It is worth stressing the application area of solar ponds for the type of work investigated here. Use of sodium chloride in a solar pond can have an adverse effect on agriculture fields, and salts such as muriate of potash and urea are being investigated as alternatives. Other phase change materials and nanofluids may lead to better thermal efficiency and storage, cf. [59], and a solvent of Navier-Stokes-Voigt type may be preferable.