Competitive Double Diffusive Convection in a Kelvin–Voigt Fluid of Order One

We present a model for convection in a Kelvin–Voigt fluid of order one when the layer is heated from below and simultaneously salted from below, a problem of competitive double diffusion since heating from below promotes instability, but salting from below is stabilizing. The instability surface threshold is calculated and this has a complex shape. The Kelvin–Voigt parameters play an important role in acting as stabilizing agents when the convection is of oscillatory type. Quantitative values of the instability surface are displayed. The nonlinear stability problem is briefly addressed.


Introduction
The classes of fluids being employed or discovered in geophysical and industrial engineering applications is increasingly diverse. Many fluids are now highly complex and are not adequately described by a stess tensor which depends linearly on the velocity gradient. Complex fluids are those which include viscoelastic fluids where the stress tensor depends on the history of the velocity gradient, and nanofluid suspensions, cf. Haavisto et al. [1], which involve a suspension of extremely fine particles in a carrier fluid. Such a suspension may require dependence on second gradients of the velocity field to incorporate physical effects such as the flattening of the parabolic profile of Poiseuille flow, cf. Straughan [2]. The field of viscoelastic fluids which includes fading memory fluids is vast, as is witnessed by the work of Amendola and Fabrizio [3], Amendola et al. [4], Anand et al. [5], Anand and Christov [6,7], Anh and Nguyet [8], Christov and Christov [9], Fabrizio et al. [10], Franchi et al. [11][12][13][14], Gatti et al. [15], Jordan et al. [16], Jordan and Puri [17], Payne and Straughan [18], Yang et al. [19].
In this article we concentrate on a particular class of viscoelastic fluids associated with the names of Kelvin and of Voigt, cf. Avalos et al. [20], Anh and Nguyet [8], Chirita and Zampoli [21], El Arwadi and Youssef [22], Layton and Rebholz [23], Rivera and Racke [24]. Analytical studies of Kelvin-Voigt fluids have been presented by Oskolkov [25,26], and generalizations of these to to incorporate temperature effects are given by Sukacheva and Matveeva [27], Matveeva [28] and Sukacheva and Kondyukov [29]. The general complex relationship between the stress and the history of the velocity gradient in a viscolelastic fluid is often approximated by including time derivatives of the stress and/or velocity gradient of various orders. These typically result in what are known as Maxwell fluids, Oldroyd fluids, or Kelvin-Voigt fluids. A very useful account of Maxwell, Oldroyd and Kelvin-Voigt fluids of various orders is given by Oskolkov and Shadiev [30] who discuss the solution existence question at length, cf. also Christov and Jordan [31].
Double diffusive convection involves the motion of a fluid (Newtonian or viscoelastic) wherein a salt is dissolved in the fluid and temperature effects are considered. This phenomenon is widely studied in the literature, see e.g. Barletta and Nield [32], Capone et al. [33], Galdi et al. [34], Gentile and Straughan [35], Harfash and Hill [36], Nield [37], Matta et al. [38], Mulone [39], Payne et al. [40], Straughan [41][42][43], Straughan [44], chapter 12, Straughan and Hutter [45]. We here present an analysis for a double diffusion problem in a Kelvin-Voigt fluid of order one. We believe this is the first presentation and analysis of this problem. We refer to the phenomenon of competitive double diffusion because we treat instability in a fluid layer which is heated from below while simultaneously subject to a heavier salt concentration below. These two physical effects are opposing and lead to interesting mathematics and physics even for a Newtonian fluid. The double diffusion problem just described is called the thermosolutal convection or solar pond problem. When the convection Rayleigh number is increased the instability Rayleigh number boundary also increases, see e.g. Joseph [46], Mulone [39], Straughan [47]. However, for a Kelvin-Voigt fluid of order one the viscoelastic term in the momentum equation also has the effect of increasing the critical instability Rayleigh number as the appropriate viscoelastic coefficient is increased. Thus, we have two physical effects which each separately increase the Rayleigh number threshold for the onset of convective motion. When both effects increase the critical Rayleigh number, one may naively expect the additive effect to increase it further. For a Newtonian fluid the effects of rotation and a vertical magnetic field each separately increase the critical Rayleigh number, Chandrasekhar [48], figure 29, p. 121, and figures 39 and 43, pages 171 and 191. However, when both effects are combined the result is not one of each being additive and having an enhanced increasing effect, see Chandrasekhar [48], figure 47, p. 203. Therefore, one has to be careful and the problem of thermosolutal convection in a Kelvin-Voigt fluid of order one as analysed here leads to a complex array of behaviours.
We now present the equations of thermal and salt convection in a Kelvin-Voigt fluid of order one. We have not seen these presented before. After that we derive a detailed linear instability analysis and global nonlinear stability analysis for these equations.

Kelvin-Voigt Order One Equations
We now denote by v(x, t), T (x, t), C(x, t), p(x, t) the velocity, temperature, concentration of a dissolved solute, and pressure at position x and time t of a body of fluid. When this fluid is of Kelvin-Voigt order one then the governing equations may be taken to be, Sukacheva and Matveeva [27], Straughan [47], where W i is a viscoelastic variable defined by (1) 5 . In additionλ, ρ 0 , ν, α, ζ, g,β, κ, κ s andγ are, respectively, the Kelvin-Voigt coefficient, a reference density, the kinematic viscosity, the coefficient of thermal expansion of the fluid, the expansion coefficient due to the solute, gravity, the viscoelastic coefficient in the momentum equation, thermal diffusivity, solute diffusivity, and the strength of viscoelasticity coefficient in the definition of W i . The symbol Δ denotes the Laplacian, k = (0, 0, 1), gravity acts downward, and the body force term in (1) 1 arises from a Boussinesq approximation employing the density for reference temperature and concentration, T 0 , C 0 , cf. the procedure in Straughan [49], section 14.1, Straughan [44], chapter 12. We employ standard indicial notation throughout together with the Einstein summation convention. For example, the divergence of the velocity field is The Kelvin-Voigt terms, involvingλ andβ arise from a constitutive theory of form where d i j = (v i, j + v j,i )/2, μ = ρ 0 ν > 0 is the dynamic viscosity, and σ i j is the Cauchy extra stress tensor, as shown by Oskolkov [25,26]. In fact, σ i j is related to the stress tensor t i j by t i j = −pδ i j + σ i j and equation (1) 1 arises from the linear momentum equation where f i is the body force. To see this observe equation (1) 5 may be integrated with an integrating factor to deduce, recalling the fading memory effect, The coefficients κ 3 and ξ 1 have form κ 3 = 2λρ 0 and ξ 1 = 2βρ 0 . It is noteworthy that the Kelvin-Voigt order one theory contains two extra parameters ξ 1 andγ which are not present in the Kelvin-Voigt order zero theory (also known as Navier-Stokes-Voigt theory), cf. Straughan [47], and these extra parameters yield a more accurate fit to experiments and provide a more refined fit to the fading memory behaviour, see e.g. Greco and Marano [50]. One should perhaps consider an objective derivative in (2), as is discussed in another context by Christov [51], cf. also Jordan et al. [52], equation (5b). However, in this article we follow the procedure of Sukacheva and Matveeva [27] and Sukacheva and Kondyukov [29].
It is important to observe that Kelvin-Voigt theory is being employed in various industrial and engineering applications to describe real materials. Gidde and Pawar [53] employ this theory to describe polydimethylsiloxane in a micropump, Jayabal et al. [54] use it to model skin in the context of the cosmetics industry, and Jozwiak et al. [55] employ Kelvin-Voigt fluids in their work on the dynamic behaviour of biopolymer materials. Erdel et al. [56] employ the complex shear moduli of a Kelvin-Voigt fluid model to calculate time-dependent coefficients for anomalous diffusion in a living cell nucleus. Askarian et al. [57] employ Kelvin-Voigt fluid models for the foundation for pipes conveying fluid. An important use of Kelvin-Voigt fluids is in the field of viscous dampers which are employed to reduce the effects of vibrations in large civil engineering structures, see e.g. Greco and Marano [50], Lewandowski and Chorazyczewski [58], Xu et al. [59]. In particular, high structures require viscous dampers, such as in the tower Taipei 101 in the city of Taipei. This tower is 1667 feet high and is very close to a fault line in the Earth's crust and it has been constructed to withstand typhoons and earthquakes. To make this possible Taipei 101 employs a 730 ton mass damper which is connected to eight viscous fluid dampers which act like shock absorbers when the mass damper moves.
The boundary conditions on the planes z = 0 and z = d are given by Thus, we analyse the interesting problem where the layer is hotter below which physically causes the fluid to expand and rise upward, whereas the layer is saltier below which causes the layer to be stable. The competition between these two effects yields an interesting physical and mathematical problem. The steady solution of interest is given bȳ where the temperature and concentration gradients, β, β s , are given by The steady pressurep may then be derived up to a constant at ones disposal from (1) 1 .
To investigate stability of the steady solution (4) we introduce perturbations The equations for the perturbations are then derived and non-dimensionalized with the scalings, cf. similar details in section 14.1 of Straughan [49], The Prandtl number, Pr, salt Prandtl number, Pc, the Lewis number, Le, and the Rayleigh and salt Rayleigh numbers, Ra = R 2 , and Rs = C 2 , are introduced as We next drop the *s and treat x i and t as the non-dimensional variables. In this manner, we arrive at the following system of non -dimensional perturbation equations arising from (1), 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 , q i , φ, θ and π satsify a plane tiling periodicity in the x, y plane.

Instability Analysis
We commence with an analysis of linear instability for the conduction solution (4). This analysis guarantees a threshold for instability. To initiate the procedure we linearize (5) and seek solutions of the form We then remove the pressure term from the equation which arises from (5) 1 by taking curl curl and retaining the third component. In this way we reduce the system to solving the equations where is a planform which reflects the shape of the instability cell, cf. Chandrasekhar [48], pp. 43-52, and Δ * h = −a 2 h, for a wavenumber a. The resulting equations hold on the domain z ∈ (0, 1), and have form where D = d/dz. The boundary conditions are We now consider two stress free surfaces and so we additionally require This allows us to write w, q 3 , φ and θ as a sin series in z of form sin(nπ z), n = 1, 2, . . . We have performed extensive computations and we found n = 1 always yields the lowest value of the Rayleigh number. Therefore, we henceforth take n = 1 and equations (8) reduce to a 4 × 4 determinant in w, q 3 , θ, φ and Λ, where Λ = π 2 + a 2 .
After some calculation one shows the Rayleigh number Ra = R 2 is given by the expression, The stationary convection threshold follows from (11) by taking σ = 0 and so we find Upon minimization in a 2 we find the critical wave number value is a 2 c = π 2 /2 and then the critical stationary convection Rayleigh number is To progress we take the real and imaginary parts of (11) and recall that R 2 must be real. To find the oscillatory convection curve we put σ = iω, ω ∈ R, and then the real part of (11) yields The imaginary part of (11) leads to the equation For equation (15) to hold we require, at least, or both conditions to hold simultaneously. Thus, (16) represent necessary conditions for oscillatory convection. Equation (15) is rearranged to yield a quadratic equation for ω 2 , namely, where the coefficients A and B are given by The critical values of the oscillatory convection threshold, for R 2 , are then found by employing the two solutions, ω 2 + and ω 2 − , of (17), i.e. 2ω 2 = −A ± √ A 2 − 4B, in (14), and minimizing Ra = R 2 in a 2 for fixed values of λ, , γ , Pr and Pc, while simultaneously ensuring ω 2 + and/or ω 2 − are positive. In fact, for all the values we tried ω 2 − is always negative. The fact that ω 2 − < 0 means that we have not found the interesting neutral curve behaviour where "instability islands" may occur in the Ra, C, space (for fixed γ ). The stationary convection curve (13) is a plane in this space (for fixed γ ) and it means we do not find closed three dimensional structures underneath the stationary convection plane given by (13), where ω 2 − > 0 inside, as is found in other areas of multi-component convection, such as double diffusive convection in a rotating fluid, first observed by Pearlstein [60], convection with temperature and two salt fields, Pearlstein et al. [61], penetrative convection with temperature and two salt fields, Straughan and Walker [62], Straughan [49], section 14.2, or inclined convection in a bidisperse porous material, Falsaperla et al. [63].
The neutral curve behaviour obtained from (13), (17) and (14) is a complex relation between Ra, C, , γ, λ, Pr and Pc, and we report numerical findings in Sect. 6.

Global Nonlinear Stability
If one writes the perturbation equations (5) in the form of an abstract equation where A, L, N are operators mapping into a suitable Hilbert space, where L is the linear operator and N (u) represent the nonlinearities, with u ≡ (u, v, w, q 1 , q 2 , q 3 , θ, φ) T , then concentrating on the part of the linear operator pertaining to the variable (w, q 3 , θ, φ) T the essential part of the linear operator has form where we have operated first on equation (5) 5 by −Δ and L A represents the skew symmetric part of L. Upon inspection of L A there are two skew symmetric effects. One is manifest in the 12 and 21 terms and reflects the presence of the viscoelastic term in (5) 1 involving . In the absence of the salt field this term yields an increasing stability threshold Ra as increases, and may lead to oscillatory convection, Hopf bifurcation, as shown in figures 2 and 4 of Straughan [64]. The second skew symmetric term in L A is manifest in the 14 and 41 terms and is due to the fact that the basic layer is heavier in salt at the base. Likewise, without the term, this salt effect leads to an increase in the Rayleigh number stability threshold as C increases. It can also lead to oscillatory convection, see Straughan [47], figure 1. Oscillatory convection is increasingly important in multicomponent hydrodynamic systems as is seen in the recent work of Rionero [65], see also Pearlstein [60], Pearlstein et al. [61], Straughan and Walker [62], Falsaperla et al. [63].
In the current problem both skew symmetric effects are present simultaneously and the joint effect is analysed numerically in Sect. 6. However, we may still obtain a global nonlinear stability result by employing the energy functional where · denotes the norm on the Hilbert space L 2 (V ), V being a period cell for the solution. From (5) and (18) one may obtain the energy equation where the production term I , and the dissipation D, have form From equation (19) one may show that a global nonlinear energy stability threshold for Ra is given by Ra ≤ 27π 4 /4, for two surfaces free of stress. Details are similar to the procedure in Straughan [49], section 14.1. This global stability threshold is depicted in Fig. 5 of this work. We believe it is possible to increase the global stability threshold by employing a generalized energy functional which contains, for example, a term of form ψ = θ − δθ for a suitable constant δ > 0, cf. Joseph [46], Mulone [39], and possibly a combination of terms involving u i , q i , ∇u i , ∇q i . This will undoubtedly involve much analysis and effort employing numerical optimization of coupling parameters with simultaneous minimization in the wavenumber for a generalized energy, cf. Straughan [66], Straughan [41]. This will also require one to involve the Kelvin-Voigt parameter λ > 0 in a non-trivial manner. This problem of determining the region of possible sub-critical instabilities is currently open.
One thing which is evident from decay of the energy (18) is that for Ra less than the global stability threshold one has decay also of ∇u and ∇q , cf. the work of Layton and Rebholz [23], on Kelvin-Voigt vortex solutions, and the numerical computations of Matveeva [28], on a Navier-Stokes-Voigt fluid.

Numerical Results
In this section we report numerical results based on the stationary convection curve (13) and the results of minimizing (14) with ω 2 given by (17). For the parameters we explore, the values of ω 2 +,− we find are always such that ω 2 − < 0 and we concentrate only on the minimization when ω 2 + > 0. Figures 1, 2, 3, 4, and 5 and Tables 1, 2, 3, and 4 are based on this strategy.
For fixed γ > 0, (13) shows the stationary convection surface is a plane in (Ra, Rs, ) space, increasing in Ra as and Rs increase. For appropriate parameter values oscillatory convection may occur and interest is when the value of Ra, is less than the stationary convection one. The two-dimensional surface of instability for oscillatory convection is found to not be a plane in (Ra, Rs, ) space, for fixed γ . When = 0 oscillatory convection may occur as shown in Straughan [47] and then the oscillatory convection branches appear to be straight lines, see Straughan [47], Fig. 1, with positions depending on λ. When C = 0 oscillatory convection may occur as shown in Straughan [64], Figs. 2 and 4. The oscillatory convection branches again appear to be straight lines which increase in , and their position depends on λ. However, when C = 0 and = 0 we find the situation is more complex. We select suitable values of parameters for realistic fluids as in Straughan [47]. To do this recall Pr = ν/κ, Pc = ν/κ c and we note the Lewis number Le = κ/κ c . In Straughan [47] Fig. 2. We see that before the transition to oscillatory convection the stationary convection value of a 2 is always π 2 /2. However, at transition this value jumps discontinuously and thereafter increases with increasing in a nonlinear way.
Since the wavenumber, a, measures the (aspect) ratio of width/depth of a convection cell, larger wavenumber values mean the aspect ratio decreases and one witnesses narrower convection cells. Thus, at transition the cell jumps to a narrower one which then becomes more narrow with increasing in a nonlinear manner as indicated in Fig. 2. Figures 3 and 4 show similar transition curves for Ra vs. and a 2 vs. when γ = 0.05, Pr = 25, Pc = 1750, C 2 = 5.
For fixed but varying C 2 the oscillatory convection curve still has a transition but on the oscillatory convection curve the Ra values increase only very slowly.  viscoelastic effect of /γ . The a 2 value over the range just reported stays the same at a 2 = 4.935, to 3 d.p. Table 1 shows the transition to oscillatory convection values when γ = 0.2, Pr = 25, Pc = 541.75, C 2 = 10, emphasizing the effect of changing λ. It is seen that increasing λ leads to an increase in the equivalent value, and an increase in a 2 , at the transition point. Table 2 fixes γ = 0.5, Pr = 25, Pc = 1750 and λ = 1. The transition to oscillatory convection values are reported as is increased from 0.6 to 2. The C 2 values decrease as is increased, with a substantial increase in the values of Ra and a 2 . The increase in the Ra values is particularly noticeable and this shows that a large value of would probably help stabilize the layer.   Tables 3 and 4 fix Pr = 25, Pc = 1750, λ = 10 −4 and they employ γ = 0.5 and γ = 0.05, respectively. C 2 is increased from 1 to 8 in both tables and the transition values of , Ra and a 2 are observed. As C 2 increases decreases, Ra decreases and a 2 decreases in both cases. The instability values of and Ra are smaller when γ = 0.05 as opposed to γ = 0.5. The a 2 values decrease for increasing C 2 but the smaller value of γ = 0.05 produces a larger range of a 2 .
It is clear that the oscillatory convection surface in (Ra, Rs, ) space is dependent on a complicated interaction between the parameters Pr, Pc, γ, C 2 and λ. Once specific values are known for a particular fluid the exact oscillatory instability surface may be calculated.
An important point to observe is that Fig. 5 shows the linear instability transition curve together with the global nonlinear stability one. In Fig. 5, C 2 = 10, but an analogous picture holds for any value of C 2 , mutatis mutandis. Given the complicated nature of the linear instability surface, it would be interesting to see nonlinear stability results obtained from a weakly nonlinear analysis, or from three-dimensional numerical simulations.

Conclusions
In this article we develop an analysis for a double diffusive convection problem in a viscoelastic fluid of Kelvin-Voigt order one type. Attention is focussed on the physically interesting and mathematically challenging 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 or for a Kelvin-Voigt fluid of order one without a salt field was suggested by Sukacheva and Matveeva [27], Sukacheva and Kondyukov [29]. An analysis of the thermal convection problem for a Kelvin-Voigt fluid of order one without a salt field is given by Straughan [64], whereas an equivalent analysis of the thermal convection problem for a Navier-Stokes-Voigt fluid is detailed by Straughan [47]. The isothermal models for a Kelvin-Voigt fluid of orders zero, one and higher are explained at length by Oskolkov [25,26], see also Oskolkov and Shadiev [30]. The Kelvin-Voigt fluid of order one contains two extra parameters than the Navier-Stokes-Voigt fluid and these provide a richer structure to the viscoelasticity of the model. From a thermal convection point of view these new effects act like a skew symmetric operator in the instability process, just as the competing effects of heating and salting below act like another, but different, skew symmetric operator. The combination of these two competing effects yields an interesting instability problem which is investigated in some detail in this article.
For small enough Rayleigh number, or small enough viscoelastic parameter (with γ fixed), the onset of thermally driven motion is by stationary convection. There is then a transition to oscillatory convection as either parameter is increased, or as both are increased. The transition surface is calculated for several realistic values of Prandtl and salt Prandtl numbers. What we see here is that the transition to oscillatory convection occurs at higher Rayleigh numbers as the Kelvin-Voigt parameter λ increases, although the transition values are strongly dependent on the salt Rayleigh number C 2 = Rs.
The increase observed as λ or increases is very important in real applications. For example, in solar pond design, or in semiconductor crystal growth from a molten liquid. 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. [67]. 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. In crystal growth, oscillatory convection may lead to striations in the semiconductor crystal which may have a detrimental effect on the ability of the crystal to function in a technical device, see e.g. Jakeman and Hu [68]. Thus an accurate description of the instability surface may help in this field also.