1 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.

2 Kelvin–Voigt Order One Equations

We now denote by \(\mathbf{v}(\mathbf{x},t),\) \(T(\mathbf{x},t),\) \(C(\mathbf{x},t)\), \(p(\mathbf{x},t)\) the velocity, temperature, concentration of a dissolved solute, and pressure at position \(\mathbf{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],

$$\begin{aligned} \begin{aligned}&(1-{{\hat{\lambda }}}\varDelta )v_{i,t}+v_jv_{i,j}=-\frac{1}{\rho _0}p_{,i}+\nu \varDelta v_i+\alpha Tgk_i-\zeta Cgk_i+{{\hat{\beta }}}\varDelta W_i\,,\\&v_{i,i}=0,\\&T_{,t}+v_iT_{,i}=\kappa \varDelta T\,,\\&C_{,t}+v_iC_{,i}=\kappa _s\varDelta C\,,\\&W_{i,t}+{{\hat{\gamma }}}W_i=v_i\,, \end{aligned} \end{aligned}$$
(1)

where \(W_i\) is a viscoelastic variable defined by (1)\(_5\). In addition \({{\hat{\lambda }}},\rho _0,\nu ,\alpha ,\zeta ,g,{{\hat{\beta }}},\kappa ,\kappa _s\) and \({{\hat{\gamma }}}\) 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 \(\varDelta \) denotes the Laplacian, \(\mathbf{k}=(0,0,1)\), gravity acts downward, and the body force term in (1)\(_1\) arises from a Boussinesq approximation employing the density

$$\begin{aligned} \rho =\rho _0[1-\alpha (T-T_0)+\zeta (C-C_0)] \end{aligned}$$

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

$$\begin{aligned} v_{i,i}\equiv \sum _{i=1}^3v_{i,i}&= \frac{\partial v_1}{\partial x_1}+\frac{\partial v_2}{\partial x_2}+\frac{\partial v_3}{\partial x_3}\\&=\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z} \end{aligned}$$

where \(\mathbf{v}=(v_1,v_2,v_3)\equiv (u,v,w)\) and \(\mathbf{x}=(x_1,x_2,x_3)\equiv (x,y,z)\). A further example is

$$\begin{aligned} v_iC_{,i}\equiv \sum ^3_{i=1}v_iC_{,i}= u\frac{\partial C}{\partial x} +v\frac{\partial C}{\partial y} +w\frac{\partial C}{\partial z}. \end{aligned}$$

The Kelvin–Voigt terms, involving \({{\hat{\lambda }}}\) and \({{\hat{\beta }}}\) arise from a constitutive theory of form

$$\begin{aligned} \sigma _{ij}=\kappa _3\frac{\partial d_{ij}}{\partial t}+2\mu d_{ij} +\xi _1\int _{-\infty }^t\exp \{-{{\hat{\gamma }}}(t-s)\}\,d_{ij}\,ds, \end{aligned}$$
(2)

where \(d_{ij}=(v_{i,j}+v_{j,i})/2\), \(\mu =\rho _0\nu >0\) is the dynamic viscosity, and \(\sigma _{ij}\) is the Cauchy extra stress tensor, as shown by Oskolkov [25, 26]. In fact, \(\sigma _{ij}\) is related to the stress tensor \(t_{ij}\) by \(t_{ij}=-p\delta _{ij}+\sigma _{ij}\) and equation (1)\(_1\) arises from the linear momentum equation

$$\begin{aligned} \rho _0(v_{i,t}+v_jv_{i,j})=t_{ji,j}+\rho _0f_i\,, \end{aligned}$$

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, \(W_i=\int _{-\infty }^t\exp \{-{{\hat{\gamma }}}(t-s)\}v_i\,ds\). The coefficients \(\kappa _3\) and \(\xi _1\) have form \(\kappa _3=2{{\hat{\lambda }}}\rho _0\) and \(\xi _1=2{{\hat{\beta }}}\rho _0\). It is noteworthy that the Kelvin–Voigt order one theory contains two extra parameters \(\xi _1\) and \({{\hat{\gamma }}}\) 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.

3 Double Diffusive Convection

We shall suppose the Kelvin–Voigt order one fluid occupies the horizontal layer \(0<z<d\) where gravity is acting downward. Equations (1) are defined on the spatial region \({\mathbb {R}}^2\times \{z\in (0,d)\}\), for \(t>0\).

The boundary conditions on the planes \(z=0\) and \(z=d\) are given by

$$\begin{aligned} \begin{aligned}&v_i=0,\quad W_i=0,\quad z=0,d;\\&T=T_L,\quad z=0,\qquad T=T_U,\quad z=d;\\&C=C_L,\quad z=0,\qquad C=C_U,\quad z=d; \end{aligned} \end{aligned}$$
(3)

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 \(^{\circ }\)K. 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 by

$$\begin{aligned} {{\bar{v}}}_i\equiv 0,\quad {{\bar{W}}}_i\equiv 0,\quad {{\bar{T}}}=-\beta z+T_L\,,\quad {{\bar{C}}}=-\beta _s z+C_L\,, \end{aligned}$$
(4)

where the temperature and concentration gradients, \(\beta ,\beta _s\), are given by

$$\begin{aligned} \beta =\frac{T_L-T_U}{d}\,, \qquad \beta _s=\frac{C_L-C_U}{d}\,. \end{aligned}$$

The steady pressure \({\bar{p}}\) 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 \((u_i,q_i,\theta ,\phi ,\pi )\) by

$$\begin{aligned} v_i={{\bar{v}}}_i+u_i,\quad W_i={{\bar{W}}}_i+q_i,\quad T={{\bar{T}}}+\theta ,\quad C={{\bar{C}}}+\phi ,\quad p={{\bar{p}}}+\pi \,. \end{aligned}$$

The equations for the perturbations are then derived and non-dimensionalized with the scalings, cf. similar details in section 14.1 of Straughan [49],

$$\begin{aligned}&x_i=x_i^*d,\qquad t=t^*{{\mathcal {T}}}, \qquad U=\frac{\nu }{d}\,,\qquad {{\mathcal {T}}}=\frac{d^2}{\nu }\,,\\&P=\frac{\rho _0\nu U}{d}\,,\qquad \epsilon =\frac{{{\hat{\beta }}}{{\mathcal {T}}}}{\nu }=\frac{{{\hat{\beta }}}d^2}{\nu ^2}=\frac{{{\hat{\beta }}}d}{U\nu }\,,\qquad \lambda =\frac{{{\hat{\lambda }}}}{{{\mathcal {T}}}\nu }=\frac{{\hat{\lambda }}}{d^2}\,,\\&\gamma ={{\hat{\gamma }}}{{\mathcal {T}}}=\frac{{{\hat{\gamma }}}d^2}{\nu }\,, \qquad T^{\sharp }=U\sqrt{\frac{\beta \nu }{\kappa g\alpha }}\,, \qquad C^{\sharp }=U\sqrt{\frac{\beta _s\nu }{\gamma g\kappa _s}}\,. \end{aligned}$$

The Prandtl number, Pr, salt Prandtl number, Pc, the Lewis number, Le, and the Rayleigh and salt Rayleigh numbers, \(Ra=R^2\), and \(Rs={{\mathcal {C}}}^2\), are introduced as

$$\begin{aligned} Ra=\frac{\alpha g\beta d^4}{\nu \kappa }\,,\quad Rs=\frac{\beta _s g\gamma d^4}{\kappa _s\nu }\,, \quad Le=\frac{\kappa }{\kappa _s}\,,\quad Pr=\frac{\nu }{\kappa }\,, \quad Pc=\frac{\nu }{\kappa _s}\,. \end{aligned}$$

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),

$$\begin{aligned} \begin{aligned}&u_{i,t}-\lambda \varDelta u_{i,t}+u_ju_{i,j}=-\pi _{,i}+\varDelta u_i+\epsilon \varDelta q_i+Rk_i\theta -{{\mathcal {C}}}\phi k_i\,,\\&u_{i,i}=0,\\&Pr(\theta _{,t}+u_i\theta _{,i})=Rw+\varDelta \theta \,,\\&Pc(\phi _{,t}+u_i\phi _{,i})={{\mathcal {C}}}w+\varDelta \phi \,,\\&q_{i,t}+\gamma q_i=u_i\,, \end{aligned} \end{aligned}$$
(5)

where \(w=u_3\).

Equations (5) are defined on the domain \({\mathbb {R}}^2\times \{z\in (0,1)\}\times \{t>0\}\) and the boundary conditions are

$$\begin{aligned} u_i=0,\quad q_i=0, \quad \phi =0,\quad \theta =0,\qquad z=0,1, \end{aligned}$$
(6)

together with the fact that \(u_i,q_i,\phi ,\theta \) and \(\pi \) satsify a plane tiling periodicity in the xy plane.

4 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 \(u_i=u_i(\mathbf{x})e^{\sigma t},\) \(q_i=q_i(\mathbf{x})e^{\sigma t},\) \(\phi =\phi (\mathbf{x})e^{\sigma t},\) \(\theta =\theta (\mathbf{x})e^{\sigma t},\) \(\pi =\pi (\mathbf{x})e^{\sigma t}.\) 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

$$\begin{aligned} \begin{aligned}&\sigma (1-\lambda \varDelta )\varDelta w=\varDelta ^2w+\epsilon \varDelta ^2q_3++R\varDelta ^*\theta -{{\mathcal {C}}}\varDelta ^*\phi ,\\&Pr\sigma \theta =Rw+\varDelta \theta \,,\\&Pc\sigma \phi ={{\mathcal {C}}}w+\varDelta \phi \,,\\&(\sigma +\gamma )q_3=w, \end{aligned} \end{aligned}$$
(7)

where \(\varDelta ^*= \partial ^2/\partial x^2 +\partial ^2/\partial y^2\) is the horizontal Laplacian. Introduce now the forms \(q_3=q_3(z)h(x,y),\) \(w=w(z)h(x,y),\) \(\theta =\theta (z)h(x,y)\), \(\phi =\phi (z)h(x,y)\), where h(xy) is a planform which reflects the shape of the instability cell, cf. Chandrasekhar [48], pp. 43-52, and \(\varDelta ^*h=-a^2h\), for a wavenumber a.

The resulting equations hold on the domain \(z\in (0,1)\), and have form

$$\begin{aligned} \begin{aligned}&\sigma (D^2-a^2)w-\sigma \lambda (D^2-a^2)^2w=(D^2-a^2)^2w+\epsilon (D^2-a^2)^2q_3\\&\quad \qquad -a^2R\theta +{{\mathcal {C}}}a^2\phi ,\\&Pr\sigma \theta =Rw+(D^2-a^2)\theta \,,\\&Pc\sigma \phi ={{\mathcal {C}}}w+(D^2-a^2)\phi \,,\\&(\sigma +\gamma )q_3=w, \end{aligned} \end{aligned}$$
(8)

where \(D=d/dz\).

The boundary conditions are

$$\begin{aligned} w=0,\quad q_3=0, \quad \theta =0,\quad \phi =0,\qquad \mathrm{on}\,\, z=0,1. \end{aligned}$$
(9)

We now consider two stress free surfaces and so we additionally require

$$\begin{aligned} D^2w=0,\qquad z=0,1. \end{aligned}$$
(10)

This allows us to write w, \(q_3\), \(\phi \) and \(\theta \) as a \(\sin \) series in z of form \(\sin (n\pi z)\), \(n=1,2,\ldots \) 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\times 4\) determinant in \(w,q_3,\theta ,\phi \) and \(\varLambda \), where \(\varLambda =\pi ^2+a^2\).

After some calculation one shows the Rayleigh number \(Ra=R^2\) is given by the expression,

$$\begin{aligned} R^2={{\mathcal {C}}}^2\biggl (\frac{\varLambda +Pr\sigma }{\varLambda +Pc\sigma }\biggr ) +\Bigl [\frac{\sigma \varLambda }{a^2}(1+\lambda \varLambda )+\frac{\varLambda ^2}{a^2}\Bigr ](\varLambda +Pr\sigma ) +\epsilon \frac{\varLambda ^2}{a^2}\biggl (\frac{Pr\sigma +\varLambda }{\sigma +\gamma }\biggr )\,.\nonumber \\ \end{aligned}$$
(11)

The stationary convection threshold follows from (11) by taking \(\sigma =0\) and so we find

$$\begin{aligned} R^2={{\mathcal {C}}}^2+\frac{\varLambda ^3}{a^2}\biggl (1+\frac{\epsilon }{\gamma }\biggr )\,. \end{aligned}$$
(12)

Upon minimization in \(a^2\) we find the critical wave number value is \(a_c^2=\pi ^2/2\) and then the critical stationary convection Rayleigh number is

$$\begin{aligned} Ra_{stat}=\frac{27\pi ^4}{4}\biggl (1+\frac{\epsilon }{\gamma }\biggr )+{{\mathcal {C}}}^2. \end{aligned}$$
(13)

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 \(\sigma =i\omega \), \(\omega \in {\mathbb {R}}\), and then the real part of (11) yields

$$\begin{aligned} R^2=\epsilon \frac{\varLambda ^2}{a^2}\Bigl (\frac{\gamma \varLambda +Pr\omega ^2}{\gamma ^2+\omega ^2}\Bigr ) +{{\mathcal {C}}}^2\Bigl (\frac{\varLambda ^2+PrPc\,\omega ^2}{Pc^2\omega ^2+\varLambda ^2}\Bigr ) -(1+\lambda \varLambda )Pr\frac{\varLambda }{a^2}\omega ^2+\frac{\varLambda ^3}{a^2}\,. \end{aligned}$$
(14)

The imaginary part of (11) leads to the equation

$$\begin{aligned} \frac{\epsilon \varLambda }{a^2}\frac{(\gamma Pr-\varLambda )}{(\gamma ^2+\omega ^2)} +\frac{{{\mathcal {C}}}^2(Pr-Pc)}{(Pc^2\omega ^2+\varLambda ^2)} +\frac{\lambda }{a^2}\,(1+\lambda \varLambda +Pr)=0. \end{aligned}$$
(15)

For equation (15) to hold we require, at least,

$$\begin{aligned} \varLambda>\gamma Pr\qquad \mathrm{or}\qquad Pc>Pr\,, \end{aligned}$$
(16)

or both conditions to hold simultaneously. Thus, (16) represent necessary conditions for oscillatory convection.

Equation (15) is rearranged to yield a quadratic equation for \(\omega ^2\), namely,

$$\begin{aligned} \omega ^4+A\omega ^2+B=0, \end{aligned}$$
(17)

where the coefficients A and B are given by

$$\begin{aligned} A=\frac{\gamma ^2Pc^2+\varLambda ^2}{Pc^2}-\frac{{{\mathcal {C}}}^2a^2(Pc-Pr)}{Pc^2\varLambda (1+\lambda \varLambda +Pr)} -\frac{\epsilon (\varLambda -\gamma Pr)}{(1+\lambda \varLambda +Pr)} \end{aligned}$$

and

$$\begin{aligned} B=\frac{\gamma ^2\varLambda ^2}{Pc^2}-\frac{{{\mathcal {C}}}^2a^2\gamma ^2(Pc-Pr)}{Pc^2\varLambda (1+\lambda \varLambda +Pr)} -\frac{\epsilon \varLambda ^2(\varLambda -\gamma Pr)}{Pc^2(1+\lambda \varLambda +Pr)}\,. \end{aligned}$$

The critical values of the oscillatory convection threshold, for \(R^2\), are then found by employing the two solutions, \(\omega _+^2\) and \(\omega ^2_-\), of (17), i.e. \(2\omega ^2=-A\pm \sqrt{A^2-4B}\), in (14), and minimizing \(Ra=R^2\) in \(a^2\) for fixed values of \(\lambda ,\epsilon ,\gamma ,Pr\) and Pc, while simultaneously ensuring \(\omega ^2_+\) and/or \(\omega ^2_-\) are positive. In fact, for all the values we tried \(\omega ^2_-\) is always negative.

The fact that \(\omega ^2_-<0\) means that we have not found the interesting neutral curve behaviour where “instability islands” may occur in the \(Ra,{{\mathcal {C}}},\epsilon \) space (for fixed \(\gamma \)). The stationary convection curve (13) is a plane in this space (for fixed \(\gamma \)) and it means we do not find closed three dimensional structures underneath the stationary convection plane given by (13), where \(\omega ^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,{{\mathcal {C}}},\epsilon ,\gamma ,\lambda ,Pr\) and Pc, and we report numerical findings in Sect. 6.

5 Global Nonlinear Stability

If one writes the perturbation equations (5) in the form of an abstract equation

$$\begin{aligned} Au_t=Lu+N(u) \end{aligned}$$

where ALN are operators mapping into a suitable Hilbert space, where L is the linear operator and N(u) represent the nonlinearities, with

$$\begin{aligned} u\equiv (u,v,w,q_1,q_2,q_3,\theta ,\phi )^T, \end{aligned}$$

then concentrating on the part of the linear operator pertaining to the variable \((w,q_3,\theta ,\phi )^T\) the essential part of the linear operator has form

$$\begin{aligned} L_A= \begin{pmatrix} \varDelta &{} \epsilon \varDelta &{} R &{} -{{\mathcal {C}}}\\ -\epsilon \varDelta &{} \varDelta &{} 0 &{} 0 \\ R &{} 0 &{} \varDelta &{} 0 \\ {{\mathcal {C}}} &{} 0 &{} 0 &{} \varDelta \end{pmatrix} \end{aligned}$$

where we have operated first on equation (5)\(_5\) by \(-\varDelta \) 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 \(\epsilon \). In the absence of the salt field this term yields an increasing stability threshold Ra as \(\epsilon \) 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 \(\epsilon \) term, this salt effect leads to an increase in the Rayleigh number stability threshold as \({{\mathcal {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

$$\begin{aligned} E(t)=\frac{1}{2}\Vert \mathbf{u}\Vert ^2+\frac{\lambda }{2}\Vert \nabla \mathbf{u}\Vert ^2+\frac{\epsilon }{2}\Vert \nabla \mathbf{q}\Vert ^2 +\frac{Pr}{2}\Vert \theta \Vert ^2 +\frac{Pc}{2}\Vert \phi \Vert ^2, \end{aligned}$$
(18)

where \(\Vert \cdot \Vert \) 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

$$\begin{aligned} \frac{dE}{dt}=RI-D, \end{aligned}$$
(19)

where the production term I, and the dissipation D, have form

$$\begin{aligned}&I(t)=2(\theta ,w),\\&D(t)=\Vert \nabla \mathbf{u}\Vert ^2+\epsilon \gamma \Vert \nabla \mathbf{q}\Vert ^2+\Vert \nabla \theta \Vert ^2+\Vert \nabla \phi \Vert ^2 \end{aligned}$$

and \((\cdot ,\cdot )\) denotes the inner product on \(L^2(V)\).

Fig. 1
figure 1

Graph of \(Ra,\epsilon \) with \(\gamma =0.2\), \(Pr=25\), \(Pc=541.75\), \({{\mathcal {C}}}^2=10\). The transition to oscillatory convection is shown for \(\lambda =2\), \(\lambda =1\) and \(\lambda =10^{-3}\). The steady solution is unstable above the appropriate stationary-oscillatory curve. The transition values are approximately \(Ra=909.779,\epsilon =0.0737\) when \(\lambda =10^{-3}\), \(Ra=1107.054,\epsilon =0.1337\) when \(\lambda =1\), and \(Ra=1304.610,\epsilon =0.1938\) when \(\lambda =2\). Kelvin–Voigt fluid of order one

Fig. 2
figure 2

Graph of \(a^2,\epsilon \) with \(\gamma =0.2\), \(Pr=25\), \(Pc=541.75\), \({{\mathcal {C}}}^2=10\). The transition to oscillatory convection is shown for \(\lambda =2\), \(\lambda =1\) and \(\lambda =10^{-3}\). The transition values are approximately \(a^2=5.053,\epsilon =0.0737\) when \(\lambda =10^{-3}\), see AB, \(a^2=5.089,\epsilon =0.1337\) when \(\lambda =1\), see CD, \(a^2=5.125,\epsilon =0.1938\) when \(\lambda =2\), see EF. The stationary convection values of \(a^2\) are \(\pi ^2/2\approx \, 4.9348\). Kelvin–Voigt fluid of order one

Fig. 3
figure 3

Graph of \(Ra,\epsilon \) with \(\gamma =0.05\), \(Pr=25\), \(Pc=1750\), \({{\mathcal {C}}}^2=5\). The transition to oscillatory convection is shown for \(\lambda =1\) and \(\lambda =10^{-4}\). The steady solution is unstable above the appropriate stationary-oscillatory curve. The transition values are approximately \(Ra=693.753,\epsilon =0.002376\) when \(\lambda =10^{-4}\), and \(Ra=729.694,\epsilon =0.00511\) when \(\lambda =1\). Kelvin–Voigt fluid of order one

From equation (19) one may show that a global nonlinear energy stability threshold for Ra is given by \(Ra\le 27\pi ^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 \(\psi =\theta -\delta \theta \) for a suitable constant \(\delta >0\), cf. Joseph [46], Mulone [39], and possibly a combination of terms involving \(u_i,q_i,\nabla u_i,\nabla 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 \(\lambda >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 \(\Vert \nabla \mathbf{u}\Vert \) and \(\Vert \nabla \mathbf{q}\Vert \), 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.

6 Numerical Results

In this section we report numerical results based on the stationary convection curve (13) and the results of minimizing (14) with \(\omega ^2\) given by (17). For the parameters we explore, the values of \(\omega ^2_{+,-}\) we find are always such that \(\omega ^2_-<0\) and we concentrate only on the minimization when \(\omega ^2_+>0\). Figures 1, 2, 3, 4, and 5 and Tables 1, 2, 3, and 4 are based on this strategy.

For fixed \(\gamma >0\), (13) shows the stationary convection surface is a plane in \((Ra,Rs,\epsilon )\) space, increasing in Ra as \(\epsilon \) 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,\epsilon )\) space, for fixed \(\gamma \). When \(\epsilon =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 \(\lambda \). When \({{\mathcal {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 \(\epsilon \), and their position depends on \(\lambda \). However, when \({{\mathcal {C}}}\ne 0\) and \(\epsilon \ne 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=\nu /\kappa ,Pc=\nu /\kappa _c\) and we note the Lewis number \(Le=\kappa /\kappa _c\). In Straughan [47] justification is given for taking Pr in the range 6.99 to 400. Typical values of \(\kappa \) and \(\kappa _c\) given there suggest \(\kappa =1.43\times 10^{-7}\,\,m^2\,\,s^{-1}\), or in the range

$$\begin{aligned} 4.4\times 10^{-8}\le \kappa \le 8\times 10^{-8}\,\,m^2\,\,s^{-1,} \end{aligned}$$

whereas

$$\begin{aligned} 1.286\times 10^{-9}\le \kappa _s\le 2.03\times 10^{-9}\,\,m^2\,\,s^{-1}. \end{aligned}$$

In Straughan [64] \(\gamma \) is taken as 0.2 or 0.5 and \(\lambda \) values of \(10^{-3},10^{-2},0.1\) and 1 are investigated. Given this information we here concentrate on \(Pr=25,Pc=541.75\) or 1750, with \(\lambda \) in the range \([10^{-4},2]\).

Figure 1 shows the stationary convection-oscillatory convection branches when \(\gamma =0.2,Pr=25,Pc=541.75,{{\mathcal {C}}}^2=10,\) with \(\lambda =10^{-3},1\) and 2. We have also computed values for \(\lambda =10^{-2}\) and 0.1, and the trend is in line with that seen here. The transition to oscillatory convection increases in \(\epsilon \) as \(\lambda \) increases but the oscillatory convection branch is not a straight line. Hence the viscoelastic effect of the parameters \(\epsilon \) and \(\gamma \) is influencing the critical Rayleigh number values and the oscillatory convection surface in \((Ra,Rs,\epsilon )\) space is a complex one dependent upon the values of \(Pr,Pc,Rs,\epsilon ,\gamma \) and \(\lambda \). The analogous values of the critcal wavenumber are shown in Fig. 2. We see that before the transition to oscillatory convection the stationary convection value of \(a^2\) is always \(\pi ^2/2\). However, at transition this value jumps discontinuously and thereafter increases with increasing \(\epsilon \) 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 \(\epsilon \) in a nonlinear manner as indicated in Fig. 2.

Figures 3 and 4 show similar transition curves for Ra vs. \(\epsilon \) and \(a^2\) vs. \(\epsilon \) when \(\gamma =0.05,Pr=25,Pc=1750,{{\mathcal {C}}}^2=5\).

Fig. 4
figure 4

Graph of \(a^2,\epsilon \) with \(\gamma =0.05\), \(Pr=25\), \(Pc=1750\), \({{\mathcal {C}}}^2=5\). The transition to oscillatory convection is shown for \(\lambda =1\) and \(\lambda =10^{-4}\). The transition values are approximately \(a^2=4.956,\epsilon =0.002376\) when \(\lambda =10^{-4}\), see AB, and \(a^2=4.960,\epsilon =0.00511\) when \(\lambda =1\), see CD. The stationary convection values of \(a^2\) are \(\pi ^2/2\approx \, 4.9348\). Kelvin–Voigt fluid of order one

Fig. 5
figure 5

Graph of \(Ra,\epsilon \) with \(\gamma =0.2\), \(Pr=25\), \(Pc=541.75\), \({{\mathcal {C}}}^2=10\). The instability and stability regoins are shown for \(\lambda =1\). The transition values are as in Fig. 1. Kelvin–Voigt fluid of order one

For fixed \(\epsilon \) but varying \({{\mathcal {C}}}^2\) the oscillatory convection curve still has a transition but on the oscillatory convection curve the Ra values increase only very slowly. For example, when \(\gamma =0.05,Pr=25,Pc=1750,\epsilon =0.003\) with \(\lambda =10^{-4}\), the Ra value increases from transition at \({{\mathcal {C}}}^2=9.90034,\) where \(Ra=671.357\), to the value \(Ra=672.593\) when \({{\mathcal {C}}}^2=100\). That would indicate that physically the salt field is having a more dominant effect on determining oscillatory convection than the viscoelastic effect of \(\epsilon /\gamma \). 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 \(\gamma =0.2,Pr=25,Pc=541.75,{{\mathcal {C}}}^2=10\), emphasizing the effect of changing \(\lambda \). It is seen that increasing \(\lambda \) leads to an increase in the equivalent \(\epsilon \) value, and an increase in \(a^2\), at the transition point.

Table 2 fixes \(\gamma =0.5,Pr=25,Pc=1750\) and \(\lambda =1\). The transition to oscillatory convection values are reported as \(\epsilon \) is increased from 0.6 to 2. The \({{\mathcal {C}}}^2\) values decrease as \(\epsilon \) 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 \(\epsilon \) would probably help stabilize the layer.

Table 1 Transition from stationary to oscillatory convection values, hence \(Ra_{osc}=Ra_{stat}\). \(\gamma =0.2\), \(Pr=25\), \(Pc=541.75\), \({{\mathcal {C}}}^2=10\). \(\lambda \) varies as shown. Kelvin–Voigt fluid of order one
Table 2 Transition from stationary to oscillatory convection values, hence \(Ra_{osc}=Ra_{stat}\)

Tables 3 and 4 fix \(Pr=25,Pc=1750,\lambda =10^{-4}\) and they employ \(\gamma =0.5\) and \(\gamma =0.05\), respectively. \({{\mathcal {C}}}^2\) is increased from 1 to 8 in both tables and the transition values of \(\epsilon ,Ra\) and \(a^2\) are observed. As \({{\mathcal {C}}}^2\) increases \(\epsilon \) decreases, Ra decreases and \(a^2\) decreases in both cases. The instability values of \(\epsilon \) and Ra are smaller when \(\gamma =0.05\) as opposed to \(\gamma =0.5\). The \(a^2\) values decrease for increasing \({{\mathcal {C}}}^2\) but the smaller value of \(\gamma =0.05\) produces a larger range of \(a^2\).

Table 3 Transition from stationary to oscillatory convection values, hence \(Ra_{osc}=Ra_{stat}\)
Table 4 Transition from stationary to oscillatory convection values, hence \(Ra_{osc}=Ra_{stat}\)

It is clear that the oscillatory convection surface in \((Ra,Rs,\epsilon )\) space is dependent on a complicated interaction between the parameters \(Pr,Pc,\gamma ,{{\mathcal {C}}}^2\) and \(\lambda \). 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, \({{\mathcal {C}}}^2=10,\) but an analogous picture holds for any value of \({{\mathcal {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.

7 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 \(\epsilon \) (with \(\gamma \) 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 \(\lambda \) increases, although the transition values are strongly dependent on the salt Rayleigh number \({{\mathcal {C}}}^2=Rs\).

The increase observed as \(\lambda \) or \(\epsilon \) 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 \(^{\circ }\)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.