Instability thresholds for thermal convection in a Kelvin–Voigt fluid of variable order

We present numerical techniques for calculating instability thresholds in a model for thermal convection in a complex viscoelastic fluid of Kelvin–Voigt type. The theory presented is valid for various orders of an exponential fading memory term, and the strategy for obtaining the neutral curves and instability thresholds is discussed in the general case. Specific numerical results are presented for a fluid of order zero, also known as a Navier–Stokes–Voigt fluid, and fluids of order 1 and 2. For the latter cases it is shown that oscillatory convection may occur, and the nature of the stationary and oscillatory convection branches is investigated in detail, including where the transition from one to the other takes place.


Introduction
The subject of hydrodynamics is one with immense application. The famous Navier-Stokes equations to describe flow of a linearly viscous, incompressible fluid are employed in many branches of applied mathematics and engineering. However, the instantaneous response of the stress to changes in the velocity gradient is not experienced by all types of fluid. There are many fluids which do not react instantaneously and instead remember their history. Such materials are usually put into the category of viscoelastic fluids which incorporate fading memory or other classes of complex fluid. However, the field of complex fluids is intensely studied, as is witnessed by the work of 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], Franchi et al. [8], Franchi et al. [9], Gatti et al. [10], Jordan et al. [11], Jordan and Puri [12], Joseph et al. [13], Payne and Straughan [14], Preziosi and Rionero [15], Yang et al. [16]. By a complex fluid we include not only those displaying viscoelastic behaviour but also such as nanofluid suspensions, cf. Haavisto et al. [17], where 1 3 the stress may need to incorporate dependence on second gradients of the velocity field to describe such physical effects as the flattening of the parabolic profile in Poiseuille flow, cf. Straughan [18]. We stress that this does not preclude complex solution behaviour in a Newtonian fluid where extremely complicated solutions are known for the Navier-Stokes equations and approximations arising therefrom such as the primitive equations, see e.g. Gargano et al. [19], Kukavica et al. [20].
A very interesting class of complex (viscoelastic) materials are those associated with the names of Kelvin and of Voigt, cf. Chirita and Zampoli [21], Layton and Rebholz [22], and interesting models for Kelvin-Voigt fluids have been presented by Oskolkov [23], Oskolkov [24]. The equations of motion for viscoelastic fluids of Maxwell type, of Oldroyd type, and of Kelvin-Voigt type are succinctly summarized in the article of Oskolkov and Shadiev [25]. Generalizations of the Kelvin-Voigt models to incorporate thermal effects are given by Sukacheva and Matveeva [26], Matveeva [27], Sukacheva and Kondyukov [28]. The latter pieces of work mainly concentrate on solvability issues and on establishing properties of solutions such as existence of a suitable solution and uniqueness questions.
In this work we present analysis of thermal convection of a layer of Kelvin-Voigt fluid heated from below employing the models of Sukacheva and Matveeva [26], Sukacheva and Kondyukov [28]. This is the first analysis of these problems we have seen. We present a computational procedure for calculating the neutral curves for instability associated with thermal convection. This is a non-trivial process for Kelvin-Voigt fluids of variable order since it transpires that to calculate the critical Rayleigh number thresholds of instability one has to first calculate the eigenvalues of the variable order theory and then employ these in the minimization process over all wave numbers. It transpires that for the models of orders 1 and 2 both stationary and oscillatory convection may occur depending on what range of parameters one employs. In practical situations this is of paramount importance, since for example in growing crystals for the semi-conductor industry, oscillatory convection may lead to the formation of striations in the crystal which in turn lead to impurities in the final crystal, cf. Jakeman and Hurle [29].
In the following we firstly describe the equations for thermal convection in a Kelvin-Voigt fluid of order 0, 1, 2, … , L, and we then outline the computational procedure in the general case. Following this, the cases of order 0, 1 and 2 are addressed in detail. Practical applications of Kelvin-Voigt fluids are mentioned at the end of Sect. 2.

The Kelvin-Voigt equations for thermal convection
Let v( , t), T( , t), p( , t) be the velocity, temperature and pressure at position and time t of a body of fluid. The Kelvin-Voigt equations of non-isothermal fluid flow of variable order may be written, see Sukacheva and Matveeva [26], Here ̂, , 0 , , g and k are the Kelvin-Voigt coefficient, the kinematic viscosity, a reference density, the coefficient of thermal expansion of the fluid, gravity, and thermal diffusivity. The symbol denotes the Laplacian, k i ≡ = (0, 0, 1) , (observe that k i and k are different), and standard indicial notation together with the Einstein summation convention is employed throughout. For example, . An example of a nonlinear term is The terms W i represent the viscoelastic fading memory effect, and ̂> 0 , are constants, = 1, … , L . In deriving equations (1) a Boussinesq approximation has been employed whereby the density in the buoyancy force is written as = 0 [1 − (T − T 0 )] for a reference temperature T 0 and this gives rise to the gTk i term with the constants appearing as a hydrostatic pressure, cf. Straughan [30], pp. 47-49, Straughan [31], pp. 16-21. As Oskolkov [23], Oskolkov [24] points out, the W i terms are equivalent to a stress relation of form where 3 > 0 and > 0 are constants, d ij is the symmetric part of the velocity gradient, and ij is the Cauchy stress tensor. This representation, Oskolkov [23], Oskolkov [24], also arises from a stress representation of form for suitable coefficients and .
If L = 0 then equations (1) are sometimes called the Navier-Stokes-Voigt equations, cf. Layton and Rebholz [22], or sometimes the Kelvin-Voigt equations of order zero, or the Oskolkov equations, see Oskolkov [24]. For L = 1, 2, … they are known as the Kelvin-Voigt equations of order L.
It is important to note that the Kelvin-Voigt theory is being used in various industrial applications to analyse real materials. For example, Gidde and Pawar [32] use this theory to describe polydimethylsiloxane in a micropump, Jayabal et al. [33] employ the theory in connection with modelling skin in the cosmetics industry, while Jozwiak et al. [34] describe the dynamic behaviour of biopolymer materials. A very important use of Kelvin-Voigt fluids is in viscous fluid dampers which are employed to reduce the effects of vibrations in civil engineering structures, as analysed by Greco and Marano [35] and Xu et al. [36]. Particular cases of this pertain to very high structures such as the tower Taipei 101 in the city Taipei. This is a 1667 ft high structure which

3
is very close to a fault line and has been constructed to withstand typhoons and earthquakes. In order to achieve this Taipei 101 has a 730 ton mass damper fitted and this is connected to 8 viscous fluid dampers which act like shock absorbers when the mass damper moves. The Kelvin-Voigt theory of order 1 has constitutive equation and this contains 2 extra parameters 1 and ̂1 which are not present in the order 0 theory. Kelvin-Voigt theory of order 2 has constitutive equation and this contains a further two parameters 2 and ̂2 and is thus capable of yielding a more accurate fit to experiments than the order 1 or 0 theory, in addition to providing a more refined fit of the fading memory behaviour, which may be required, see e.g. Greco and Marano [35].

Thermal convection in a fluid of variable order
We suppose the fluid occupies the horizontal layer 0 < z < d with gravity acting downward. Then equations (1) are defined on the spatial region The boundary conditions on the planes z = 0 and z = d are given by with constant values T L • K, T U • K, where • K denotes degrees Kelvin, and T L > T U > 0. The steady solution in whose stability we are interested is given by where the temperature gradient, , is given by and the steady pressure p is a quadratic in z which follows from equation (1) To investigate stability of the solution (3) we introduce perturbations (u i , , , q i ) by The equations for the perturbations are then derived and non-dimensionalized with the scalings where T, U, W and T ♯ are the scales for time, velocity, the vectors W i , and temperature, respectively. Define the Rayleigh number, Ra = R 2 , and the Prandtl number Pr, by and rescale the coefficients ̂ , ̂ , by the non-dimensional ones If we drop the *s and now treat x i and t as the non-dimensional variables, the non-dimensional form of the perturbation equations to (1) is Equations (4) are defined on the domain ℝ 2 × {z ∈ (0, 1)} × {t > 0} and the boundary conditions are together with the fact that u i , q i , and satisfy a plane tiling periodicity in the x, y plane.
To analyse linear instability we linearize (4) and seek solutions of the form u i = u i ( )e t , q i = q i ( )e t , = ( )e t , = ( )e t . This yields from (4) 4

Then we utilize this in (4) to derive
We operate on (7) 1 with curl curl and retain the third component of the resulting equation to find (4) , is a planform which reflects the shape of the instability cell, cf. Chandrasekhar [37], pp. 43-52, and * h = −a 2 h , for a wavenumber a. Then (8) may be rearranged as where D = d∕dz and Eq. (9) hold on z ∈ (0, 1).
To progress from this point we must specify further information on boundary conditions. The standard boundary conditions for (9) are If the surfaces are fixed then we additionally require that In this case we must employ a numerical technique like the Chebyshev -tau method, see Dongarra et al. [38], together with the QZ-algorithm of Moler and Stewart [39]. This then requires resolution of an L + 2 order eigenvalue problem for the L + 2 eigenvalues before one seeks the particular eigenvalue leading to a minimum problem for the Rayleigh number. Since we believe this is the first analysis to address the thermal convection problem for a Kelvin-Voigt fluid we restrict attention to two stress free surfaces whence (11) is replaced by In this case we write w and as a sin series in z of form sin(n z) , n = 1, 2, … , cf. Chandrasekhar [37]. Our computations in the cases of Kelvin-Voigt fluids of orders 0, 1 and 2, have always yielded the lowest Rayleigh number instability curves when n = 1 , and so in our computations reported later we employ this value.
For the case of stationary convection where the growth rate is real the instability boundary = 0 yields from (9) the critical Rayleigh number However, we cannot neglect oscillatory convection and then for the instability boundary Re( ) = 0 we may take = i in (9). Converting the terms 1∕( + ) by multiplying top and bottom by the conjugate we find Eq. (9) reduce to in the complex case Insertion of the form sin n z with n = 1 then yields a determinant equation from (14) of form Equation (15) must now be resolved into its real and imaginary parts. Recognise the fact that R 2 is real and then the imaginary part yields a polynomial equation of order L + 2 in 2 . This must be solved numerically for each 2 and then be adaptively employed to calculate the minimum of Ra = R 2 in a 2 from the real part of the equation. In this way a network of eigenvalue problems is solved which lead to the solution of the minimum critical Rayleigh number which yields the instability threshold.
The procedure is described in detail in Sects. 5 and 6 where we find the instability thresholds for the order 1 and order 2 problems.

Stability for thermal convection in a Navier-Stokes-Voigt fluid
When L = 0 and the terms are not present then (4) reduce to those of Navier-Stokes-Voigt theory incorporating thermal effects, cf. Layton and Rebholz [22], Sukacheva and Kondyukov [28]. The nonlinear perturbation equations in this case become on ℝ 2 × {z ∈ (0, 1)} × {t > 0} , together with the boundary conditions and periodicity in the x, y directions.
Let (⋅, ⋅) and ‖ ⋅ ‖ denote the inner product and norm on the Hilbert space L 2 (V) , where V is a period cell for the solution.
Equation (16) may be written in the abstract form cf. Straughan [30], pp. 77-87, where the linear operator L is here symmetric. Then, see general details in Straughan [30], pp. 77-87, one may show that for system (16) the linear instability boundary is equivalent to the nonlinear instability threshold. Thus, we establish the optimum result guaranteeing no subcritical instabilities. One may also show that the instability threshold is the same as that of classical Bénard convection, cf.
For the present problem one finds the decay rate of a solution is different from that observed in classical Bénard convection due to the presence of the Kelvin-Voigt term u i,t and decay is found in the "energy" measure This faster decay rate is also observed by Layton and Rebholz [22], in their Kelvin vortex solutions, see also numerical computations of Matveeva [27].
Remark We may add a Rayleigh friction term to the right hand side of Eq. (16) 1 as in Di Plinio et al. [40]. In this case the perturbation equations (16) are replaced by for a constant > 0 . Again, if one writes (18) in the operator form (17) then the linear operator L is symmetric and one may derive the optimal result that the linear instability and global nonlinear stability boundaries coincide. One may further demonstrate that the Rayleigh number instability threshold is found by minimizing in a 2 the expression where = 2 + a 2 , a being the wavenumber. This threshold is mathematically identical to the one found for thermal convection in a Brinkman porous material. The threshold in the Brinkman case is analysed in depth in Rees [41].

Instability for a Kelvin-Voigt fluid of order one
The Kelvin-Voigt perturbation equations of order 1 are given by (4) with L = 1 . Thus, we must analyse a solution to the system of equations where we have made the identification q 1 i ≡ q i , 1 ≡ , 1 ≡ . The boundary conditions are together with periodicity in x and y.
To develop an energy stability analysis from (19) (22) to arrive at From this equation we may apply the method of energy stability theory, cf. Straughan [30], chapter 4, to find global nonlinear stability provided R does not exceed the classical Bénard threshold. For example, provided Ra < 27 4 ∕4 for two free surfaces, cf. Fig. 1, the region below the line OC. However, as we see from Fig. 1, the linear instability curve increases with beyond 27 4 ∕4 , the lines OA and AB, and there is a region of possible sub-critical instability, in the area OABC. In this respect, the Kelvin-Voigt problem of order 1 behaves very similarly to the thermosolutal convection problem; Fig. 1 is mathematically similar to figure 1, Straughan [46], p. 620.
To determine the instability curves in Fig. 1 we return to Eq. (9) and take L = 1 . For two free surfaces we find that these equations yield the Rayleigh number relation where = 2 + a 2 . For stationary convection = 0 and then one obtains (23) R 2 a 2 = ( + Pr ) ( + 2 ) + 2 + + 2 , and minimizing over a 2 yields a 2 = 2 ∕2 , whence We might expect oscillatory instability for relatively small compared to since in the formal limit → 0 , u i = q i,t and then equation (19) 1 essentially becomes hyperbolic, and such instabilities have been seen in other viscoelastic fluids in similar circumstances, see e.g. Joseph et al. [13].
To determine the oscillatory boundary threshold we follow the method in Chandrasekhar [37], section 29. Since we require the boundary where Re = 0 we put = i , ∈ ℝ , and then from (23) one finds Since R 2 and a 2 are real, (25) yields two equations after separating into real and imaginary parts. We require ≠ 0 and then the imaginary part yields From this relation one obtains necessary conditions for instability of form The real part of (25) yields One now employs 2 from (26) in (27) and minimizes the resulting expression numerically in a 2 for fixed , and Pr. Whether stationary or oscillatory convection occurs depends on the relative sizes of R 2 given by (24) or the minimum of (27). We now present numerical results for the minimization of (27). To present graphical output we fix and allow to vary so that the stationmary convection curve is a straight line as shown as line OA in Fig. 1. For oils the value of the Prandtl number lies typically in the range 25 ≤ Pr ≤ 400 , see Rodenbush et al. [54].
The stationary convection curve and oscillatory convection curves are given for = 0.2, Pr = 25 with taking values = 1, 10 −3 , in Fig. 2, and the transition values as the instability switches from stationary to oscillatory convection are given in Table 1 for = 10 −3 , 10 −2 , 0.1 and 1. The equivalent critical a 2 values are shown in Fig. 3. We see that as the instability switches from stationary to oscillatory convection the a 2 values increase discontinuously and since a ∝ 1∕aspect ratio, (aspect ratio = width/depth of a convection cell), this means that the cells switch from a wider to a relatively narrower cell once the convection is oscillatory. Figures 4 and 5 give the Ra against values and the a 2 against values for = 0.5 , Pr = 25 and the transition values are given in Table 2. It is noticeable that when = 0.5 the a 2 values at = 1 are considerably smaller for oscillatory convection than the = 10 −3 ones. The opposite is true when = 0.2.
We computed instability thresholds also for Pr = 100 . For the parameter ranges we considered only stationary convection was discovered when Pr = 100 . It would appear and ( − 2 ) > Pr + 2 (1 + Pr).

Instability for a Kelvin-Voigt fluid of order two
We now consider the instability of the conduction solution to the thermal convection problem for a Kelvin-Voigt fluid of order 2. From Eq. (4) the governing system of perturbation equations may be written where V i avd W i correspond to W 1 i and W 2 i . The boundary conditions are  where a is the wavenumber and = 2 + a 2 . From (31) we find the stationary convection boundary, = 0 , as To investigate oscillatory convection take = i , ∈ ℝ , and then (31) is written as the real part and the imaginary part. For ≠ 0 the imaginary part yields the quadratic equation , corresponding to the + and -sign in Then, both solutions 2 + and 2 − are employed in the equation for R 2 which arises from the real part of equation (31). This equation has form (33) X 1 4 + X 2 2 + X 3 = 0, Expression (34) is minimized in a 2 for each of 2 + and 2 − and the correct minimum is chosen provided the appropriate 2 is positive. In this manner we find the oscillatory convection instability threshold which is compared to the stationary convection one (32) to assess which type of convection actually occurs.
We present two tables of numerical results taking Pr = 25 in both. Table 3   presented in both tables we found 2 + > 0 and 2 − < 0 , although if and are relatively large then only stationary convection is found.
In Table 3 we see that relatively small values of and compared to and result in oscillatory convection, a fact in keeping with the study in a Kelvin-Voigt fluid of order 1 as seen in Sect. 5. When is kept at fixed value 0.6 stationary convection is witnessed when = 0.5 but once is lowered there is a transition point and thereafter lower values of yield oscillatory convection. For = 0.6 we find an approximate transition value at = 0.474 . However, when we keep fixed at value 0.4 and increase we were unable to find stationary convection, even when = 2.
In Table 4 the value of is larger, but and are smaller. Here we again witness a transition from stationary to oscillatory convection as the pair ( , ) moves from (0.4,0.3) to (0.35,0.25). The critical Rayleigh number values are much higher in Table 4 than those on Table 3. This shows the layer is more stable for the parameters of Table 4.
In Table 3, oscillatory convection is associated to smaller wave numbers and, therefore, to wider convection cells. This trend is also witnessed in Table 4.
In general, we believe there is a three-dimensional picture where if we fix and and plot Ra against and we shall find a transition from stationary to oscillatory convection surface analogous to the Ra against curves in Sect. 5. It is interesting to note that in the order 2 case the oscillatory convection threshold Rayleigh numbers can fall below the transition between the stationary and oscillatory value. Thus, a complete picture in the order 2 case will require substantial computation to create the two dimensional surface as , , and are varied, once real fluid parameter values are identified.

Conclusions
We have analysed the models for viscoelastic flow incorporating thermal effects suggested by Sukacheva and Matveeva [26], Sukacheva and Kondyukov [28], these following from the isothermal models of Oskolkov [23,24]. A procedure for computing the growth rates and then the stability threshold curves is presented for thermal convection in a layer. In general, the computational procedure needs to be adaptive and will benefit by coding over a network.
We have analysed in detail thermal convection in Kelvin-Voigt fluids of order 0, 1 and 2. The order 0 case is special in that we showed the linear instability curve coincides with the global nonlinear stability one. For Kelvin-Voigt fluids of order 1 and 2 this is not the case and there we witness oscillatory convection and more complex neutral curve behaviour. The order 1 case appears to have a transition value and then the oscillatory curve increases with increasing values (for fixed ), as is shown in Figs. 1, 2 and 4. The order 2 case is more complex and yields a two dimensional surface in the Ra, , three dimensional space for fixed , . We stress that the computation of oscillatory curves or surfaces is very important in real life fluid mechanics, since as pointed out in the introduction there are situations such as semi-conductor crystal growth where such instabilities are best avoided.
For higher order Kelvin-Voigt fluids the growth rate will satisfy a polynomial equation of increasing order. The neutral curves must be based on each value of 2 which is found, and the critical Rayleigh number determined taking all 2 into account. There may, therefore, be the opportunity for even greater complex behaviour.