Natural convection of a viscoelastic Cattaneo–Christov fluid bounded by thick walls with finite thermal conductivity

1 Introduction

Buoyancy convection phenomena occur in nature and in many industrial applications. This is the reason why it has been investigated for many years ago [1]. Particular interest is given to viscoelastic fluids which also have a variety of applications. Natural convection of viscoelastic fluids (polymer melts, polymeric solutions) is an important subject of research because this phenomenon appears in many applications due to the lack of control of temperature gradients when the heating conditions and atmospheric temperature change independently.

Viscoelastic convection has found relevant applications in the replication of DNA molecules in suspension [2–5]. Besides, it is used in the formation of patterned surfaces for electronic devises, by means of the evaporation of convective patterns of polymeric suspensions [6–8]. Geophysics is also a mayor subject where viscoelastic models are used to describe mantle convection [9, 10]. In the process of oil recovery natural convection in porous media is also applicable mainly in the separation of viscoelastic species [11].

The basics of linear viscoelastic convection was investigated by many authors. Vest and Arpaci [12] and Sokolov and Tanner [13] investigated convection in a Maxwell fluid. Takashima [14] considered the Oldroyd model. Kolkka and Ierley [15] first used the fixed heat flux boundary condition. Martínez-Mardones and Pérez-García [16] found a codimension-two point between stationary and oscillatory convection for very good conducting walls. Pérez-Reyes and Dávalos-Orozco [17] found codimension-two points of a Maxwell fluid when the thicknesses and conductivities of the walls are assumed to be finite. The case of an Oldroyd fluid model was explored in [18]. When the free surface is deformable the problem was studied in [19]. A review on this subject can be found in [20].

Nonlinear viscoelastic convection has been investigated in a number of publications. Rosenblat [21] investigated the possibility of supercritical or subcritical instability. Dávalos and Manero [22] considered convection of a second order fluid when the heat conductivity of the walls is very small. In Park and Lee [23, 24] it was assumed that the walls are perfectly conducting. Martínez-Mardones et al. [25] calculate an amplitude equation in rder to determine the preferred pattern of viscoelastic convection. Coupled Ginsburg–Landau equations are found in [26] to describe nonlinear oscillatory convection. Kaloni and Lou [27] used the energy method in a problem where the temperature has an inclined gradient. Li and Khayat [28] used an amplitude equation in three-dimensions to determine the convective pattern in an Oldroyd fluid. The Maxwell model is investigated by Salm and Lücke [29] by means of an amplitude equation and numerical analysis. Pérez-Reyes and Dávalos-Orozco [30] examined the convection of an Oldroyd fluid when the walls are bad conductors and derived a nonlinear partial differential evolution equation to determine the preferred convective pattern. Albaalbaki et al. [31] considered the convection of a shear-thinning viscoelastic fluid. Rotating convection with temperature modulation was investigated by Gopal and Narasimhamurthy [32]. Bistable bifurcation is investigated by Rebhi et al. [33] when double diffusive convection is also taken into account.

The heat transfer with Cattaneo–Christov heat flux in the boundary layer flow of a viscoelastic liquid is investigated in [34, 35] and in nanofluids in [36, 37]. However, the natural convection of a Cattaneo–Christov horizontal fluid layer inside two thick walls of finite thermal conductivity [17, 18, 38, 39], has not been investigated when the liquid satisfies the Maxwell viscoelastic constitutive equation. Therefore, it is the goal of the present research to investigate the natural convection of a liquid layer confined between two thick walls of finite thermal conductivity which presents two relaxations times. One corresponds to that of the Cattaneo–Christov heat flux constitutive equation and the other one corresponds to that of the Maxwell viscoelastic fluid constitutive equation. As will be seen presently, the results are more complex than those obtained before for a viscoelastic fluid [17, 18]. Therefore, the main goal is to find the curves of criticality for different Prandtl numbers. All criticality curves are calculated against log X, where X is the ratio of liquid over wall thermal conductivities. Two magnitudes of the wall over the liquid thicknesses ratio d are used. These two parameters, along with the two relaxation times, give a variety of new and interesting results.

2 Equations of motion and energy

The usual constitutive equation for the heat flux in a continuous medium is the so called Fourier law, where the heat flux depends linearly on the temperature gradient and the heat propagates with an infinity velocity in the medium, which in many cases is in contradiction with experiment. Research has been done to formulate more general constitutive equations in closer agreement with nature, where heat propagates with finite velocity. Some of those formulations change the energy equation of parabolic type into one of hyperbolic type, that is, into a wave equation with finite wave velocity. In this way, Cattaneo [40] proposed a constitutive model which presents a relaxation time in a first order time differential equation for the heat flux, in more agreement with experiments using different materials. This model was generalized with respect to the time derivative in [41] and using different thermodynamic models in [42–44]. This model is also called of Maxwell–Cattaneo [45] because of its similarity with the Maxwell’s model of viscoelasticity. The equation was generalized into the Cattaneo–Christov model with an upper-convected time derivative for vectors, invariant under a change of reference frame by Christov [46]. Reviews on this subject can be found in [47–49]. Natural convection of the so called Maxwell–Cattaneo fluid was investigated for many years ago [50–55]. However, the Cattaneo–Christov model in natural convection was applied first by [56, 57]. A comparison of three kinds of time derivatives in the Maxwell–Cattaneo model was done in [58]. It is important to note that Bissell [57] was the first to investigate in detail, in a wide range of Prandtl numbers, the codimension-two points where stationary and oscillatory instabilities compete to be the first unstable one, when the walls are very good conductors.

In this paper, natural convection of a viscoelastic liquid layer whose heat flux follows the Cattaneo–Christov constitutive model equation is investigated. Further, it is also assumed that the walls are thick and have finite thermal conductivity. A sketch of the system is shown in Figure 1. As can be seen, the liquid layer is confined between two parallel walls perpendicular to gravity. The lower atmosphere below the lower wall has a higher temperature than the upper atmosphere above the upper wall. The thicknesses ratio and thermal conductivities ratio will play an important role on the stability. In Figure 1, the origin is set at the interface between wall 1, of thickness d _W1, and the viscoelastic liquid layer of thickness d _f. The interface between the fluid and wall 2, of thickness d _W2, is located at d _f.

Figure 1:

A fluid confined between two parallel walls at different temperatures and perpendicular to gravity. The origin is set at the interface between wall 1, of thickness d _W1, and the fluid layer of thickness d _f. The interface between the fluid and wall 2, of thickness d _W2, is located at d _f.

Thus, under the Boussinesq approximation, the equations of motion, continuity, energy and the viscoelastic fluid Maxwell constitutive equation are:

where λ _M is the relaxation time of the of the Maxwell viscoelastic fluid [59]. The upper convected derivative for tensors [59] is defined as:

where the superscript T ̲ means transpose of the tensor. It is important to note that when the initial state is hydrostatic, the linear perturbation of Eq. (4) reduces the right hand side to a partial time derivative, that is, ∂τ′/∂t. The energy equation and the heat flux Cattaneo–Christov [46] constitutive equation are:

with the Cattaneo–Christov relaxation time λ _F. Here the upper-convected time derivative for vectors [46] is defined as

Where u ′ ⃗ is the velocity vector, τ′ is the shear stress tensor, e′ is the shear rate tensor, p′ is the pressure, g ⃗ is the acceleration of gravity vector, T′ is the temperature, q ′ ⃗ is the fluid heat flux vector, μ is the dynamic viscosity, ρ is the density and ρ ₀ is a constant reference density. c _p is the specific heat at constant pressure, K _F is the fluid thermal conductivity and ∇′ is the nabla operator. It is assumed that the density satisfies the constitutive relation

where α _t is the coefficient of volumetric expansion. The energy equations of the walls are as follows. For wall 1:

with T _W1 the temperature, ρ _0W1 the density, c _pW1 the specific heat at constant pressure and K _W1 the thermal conductivity of wall 1. For wall 2:

with T _W2 the temperature, ρ _0W2 the density, c _pW2 the specific heat at constant pressure and K _W2 the thermal conductivity of wall 2.

Before onset, the fluid remains in a hydrostatic state. Thus, its main temperature profile satisfies:

In the same way, the main temperature profiles of wall 1 and wall 2 are

Here, T 0 F ′ , T 0W 1 ′ and T 0W 2 ′ are the basic temperatures of the fluid, wall 1 and wall 2, respectively. The boundary conditions to be satisfied are:

where T _0L is the temperature of the lower atmosphere below the lower wall 1 and T _0U is the temperature of the upper atmosphere above the upper wall 2. The fluid is heated from below and it is assumed that T _0L > T _0U.

The solutions of this system of three differential equations, presented in non-dimensional form, are

Here, distances have been made non-dimensional with d _F, time with d F 2 / α , velocity with α/d _F and pressure with ρ 0 α ν / d F 2 , where α = K _F/ρ ₀ c _p is the fluid thermal diffusivity and ν = μ/ρ ₀ is the kinematic viscosity. The temperatures are made non-dimensional by first subtracting T _0U and then dividing by ΔT = (T _0L − T _0U)/(1 + D _L X _L + D _U X _U). The heat flux can be made non-dimensional with K _FΔT/d _F. The parameters D _L and D _U are the wall 1 over fluid thicknesses ratio and wall 2 over fluid thicknesses ratio, respectively. Besides, X _L and X _U are the fluid over wall 1 heat conductivities ratio and fluid over wall 2 heat conductivities ratio, respectively. Thus, T _L = 1 + D _L X _L + D _U X _U at z = −D _L and T _U = 0 at z = 1 + D _U. It is easy to check that T _F satisfies the boundary conditions at z = 0 and z = 1. Notice that these solutions Eq. (15) of the basic temperatures profiles were not presented in the paper by Cerisier et al. [38], but were presented in [39] for a different coordinate system.

In hydrostatic conditions the pressure in the liquid satisfies the non-dimensional equation

where it is assumed that the reference temperature T ₀, in non-dimensional form, is T _F(0) = 1 + D _U X _U. Pr = ν/α, the ratio of the kinematic viscosity over the thermal diffusivity, is the Prandtl number and Ga = g d F 3 / ν 2 , the ratio of the gravity force over viscous forces, is the Galilei number. The Rayleigh number Ra, the ratio of the buoyancy force due to a temperature gradient over viscous forces, is defined as:

Equation (16) will be subtracted from the equations of motion after application of a perturbation to all the variables. In this way, the equations of the linear perturbation in non-dimensional form are

In Eq. (20) use is made of dT _F/dz = −1. Besides, in Eq. (21) use is made of the fact that the main heat flux of the fluid satisfies q ₀ = −dT _F/dz = 1 and dq ₀/dz = −d² T _F/dz ² = 0. Functions without primes are perturbation variables. Here, w is the third component of the velocity perturbation u ⃗ and θ is the perturbation of temperature. The non-dimensional F is the ratio of relaxation time of the Maxwell viscoelastic fluid over the thermal diffusion time across the liquid layer. The non-dimensional ɛ _t is the Cattaneo–Christov relaxation time of the fluid perturbation heat flux q ⃗ over the thermal diffusion time.

The energy equations of the walls satisfy:

respectively, with θ _WL the temperature perturbation of wall 1 and Y _L = ρ _0W1 c _pW1/ρ ₀ c _p. Notice that in fact the product X _L Y _L is the ratio of thermal diffusivities of the fluid over that of wall 1. Besides, θ _WU is the temperature perturbation of wall 2 and Y _U = ρ _0W2 c _pW2/ρ ₀ c _p. Notice again that the product X _U Y _U is the ratio of thermal diffusivities of the fluid over that of the wall 2.

It is possible to write Eq. (20) only in terms of θ applying the operator 1 + ɛ _t ∂/∂t to both sides and using Eq. (21) in the right hand side. That is

where use is made of the continuity equation for the velocity perturbation.

In a similar way as in Eq. (24) it is possible to apply the differential operator 1 + F∂/∂t to both sides of Eq. (17) to replace τ using Eq. (19) and the continuity equation for the perturbation Eq. (18). Then, the third component of the second rotational of this equation is:

where ∇ ⊥ 2 = ∂²/∂x ²+∂²/∂y ² is the horizontal laplacian.

The number of parameters is very large. Thus, a way to understand the convection phenomena in this system is to assume that the two walls have the same thickness and thermal conductivity. Therefore, suppose that X = X _L = X _U, D = D _L = D _U and Y = Y _L = Y _U. It is also assumed that the perturbed variables are in the form of normal modes f(x, y, z, t) = F(z) exp(ik _x x + ik _y y+σt). Here, i is the imaginary number, σ = iω where ω is the frequency of oscillation and F(z) represents the amplitude of any of the variables. k _x and k _y are the x and y components of the wavenumber whose magnitude is k = k x 2 + k y 2 .

In this way, in normal modes these two Eqs. (24) and (25) can be written as:

Notice that D _z = d/dz. When σ = 0, the set of Eqs. (26) and (27) reduces to that of Chandrasekhar [1]. When ɛ _t is zero, the equations reduce to those of Vest and Arpaci [12] for a viscoelastic fluid. In the case F = 0, that set of equations reduces to those of Straugham [56] and Bissel [57] for natural convection with Cattaneo–Christov heat flux. The walls energy equations are:

Here, q = k 2 + σ X Y . The solution of Eq. (28), using Θ_WL = 0 at z = −D, has the form:

where A _WL is an integration constant.

Thus, Θ_WL(z) should satisfy the following boundary conditions:

From these two equations it is possible to eliminate A _WL and obtain the boundary condition for Θ, the amplitude of the temperature of the fluid. That is:

where

is the Biot number, which depends on non-dimensional parameters of the fluid and wall, that is X, Y, D, k and σ. Notice that for very small D, Bi is nearly independent of q and therefore of the wavenumber k and frequency ω, where σ = iω. Notice that this Biot number generalizes the one derived by Cerisier et al. [38] and Pérez-Reyes and Dávalos-Orozco [17, 18].

The boundary condition at the interface z = 1 can be obtained with a similar method. The result for the upper wall is:

The two Eqs. (26) and (27), along with the two boundary conditions Eqs. (31) and (32), and those of W for a solid boundary:

constitute a proper value problem for Ra which will be solved numerically in the next section. Though convective patterns [60] are of great importance in natural convection, this linear proper value problem only allows us to know their magnitude through the wavenumber.

The boundary conditions play an important role in natural convection. The thermal boundary conditions usually are selected as the ideal fixed temperature at the fluid-wall interface [1]. The non-dimensional Biot (Bi) number represents the flow of heat through the fluid-wall interface. In the limit of Bi → ∞ the wall is a very good conductor in comparison with the fluid. This limit corresponds to the ideal fixed temperature at the wall. The other limit of Bi → 0 corresponds to the ideal fixed heat flux boundary condition where the wall is a very bad conductor in comparison with the liquid. This last case also has attracted wide attention [22, 30, 61], [62], [63], [64]. In real applications the Biot number has intermediate values. The thickness and thermal conductivity of the wall play an important role and thus they appear in the definition of the Biot number (see below Eq. (31)). Related theoretical research on natural convection can be found in [38, 39]. Oscillatory viscoelastic convection is investigated in [17, 18] and viscoelastic thermocapillary convection in [65].

3 Numerical analysis and convective stability results

The proper value problem for Ra will be solved numerically applying the Galerkin method as in the book by Gershuni and Zhukhovitskii [66] (see also [67, 68]). In the book [66] only the lowest order is used. Thus, it was necessary to implement and use up to the third order approximation.

The third component of the velocity amplitude W is approximated by the following polynomial function, which satisfies the four boundary conditions Eq. (33):

This approximation is substituted into Eq. (26)

Thus, this second order inhomogeneous ordinary differential equation is solved exactly for Θ using the boundary conditions Eqs. (31) and (32). Let us call this solution Θ _s (z). Then, Θ _s (z) and the approximation of W(z) Eq. (34) are substituted into Eq. (27):

In order to find the solution of A ₁, B ₁ and C ₁, a set of three coupled algebraic homogeneous equations is obtained as follows. First, notice that there is a residue of the approximation in Eq. (36). This residue must be orthogonal to each function z 2 z − 1 2 , z 3 z − 1 3 and z 4 z − 1 4 in Eq. (34) (see page 111 of [68]). Then, the first algebraic homogeneous equation for A ₁, B ₁ and C ₁ is obtained by multiplying Eq. (36) by z 2 z − 1 2 and integrating in the range z = 0 to 1. That is:

The second and third algebraic equations are obtained multiplying Eq. (36) by z 3 z − 1 3 and z 4 z − 1 4 , respectively, and integrating in the range z = 0 to 1. That is:

The determinant of this set of three homogeneous algebraic equations has to be zero in order to have a non trivial solution for A ₁, B ₁ and C ₁. The solutions of this determinant are the proper values of the Rayleigh number Ra. The Rayleigh number appears in the three equations and therefore a cubic equation is expected for Ra. This cubic equation has complex coefficients and consequently the three solutions for Ra are complex. However, the Rayleigh number should be real and the imaginary part of each solution is set equal zero. For a given k, the root of each imaginary part is the frequency of oscillation ω. This root ω is substituted back into the corresponding real part to obtain the marginal Rayleigh number for the given wavenumber k and for all the other parameters fixed. The critical Rayleigh number Ra_c corresponds to the minimum of the marginal curve with the corresponding critical wavenumber k _c and frequency ω _c and all the other parameters fixed. It is important to point out that sometimes the three roots of Ra are needed to trace one curve of criticality in the range of log X. Though the same algorithm may be used for stationary convection where the coefficients of the cubic equation are real, the following results show that it is not possible to find codimension-two points because the curves of criticality Ra_c are located far below the curves of criticality of stationary convection [17].

The numerical method was validated comparing with well known papers found in the open literature, closely related with the problem under investigation. First, comparison was done with papers where use is made of the Cattaneo–Christov model. Very good comparison was obtained (up to three decimals) with Table 1 of the work by Straughan [56] who used the D ² Chebyshev tau numerical method. Next, validation was made with the results of Bissell [57]. Notice that Bissel uses a different notation and that ɛ _t equals his 2 C _T and that ω/Pr equals his γ _T . Bissell [57] used the Chebyshev tau-QZ method to obtain the codimension-two points. Calculations were made for seven of the Prandtl numbers of his Table 1 and very good agreement was found with the two decimals he used. Comparison was also made with respect to classical results of Maxwell viscoelastic convection. Kolkka and Ierley [15] made calculations for fixed temperature and fixed heat flux at the walls. They used the spectral Chebyshev tau method and a QZ algorithm for the numerical analysis. A very good agreement was found up to three decimals. Validation was also made with the work by Martínez-Mardones and Pérez-García [16] who converted the boundary value problem into an initial value problem which was solved with the fourth-order Runge–Kutta method. The marginal curve in their Figure 4 corresponding to their Λ = 0 (Maxwell viscoelastic fluid) was reproduced and the agreement was very good.

Therefore, the discussion of the results will be restricted to oscillatory convection of a viscoelastic Maxwell fluid with Cattaneo–Christov heat flux. First the case of Pr = 10 is discussed. The curves of criticality for Ra_c against log X are shown in Figure 2. In Figure 2(a) it is assumed that the viscoelastic relaxation time is F = 0.1. The left sub-figure corresponds to the relative wall thickness D = 0.1 and the right one to D = 100. In both sub-figures the important ratio of densities and heat capacities is Y = 0.1 (solid), 1 (dashed), 10 (dotted). As can be seen, the increase of the heat conductivities ratio X stabilizes convection in the case of a purely viscoelastic fluid (ɛ _t = 0) by increasing the critical Rayleigh number. However, when the liquid is viscoelastic with Cattaneo–Christov heat flux (ɛ _t > 0) the contrary occurs. This effect is very small when ɛ _t = 0.07. Moreover, an increase of ɛ _t increases all the curves of criticality. An important characteristic of the curves of criticality is that, when ɛ _t > 0, a bump is found in the middle range of log X. The height of the bump increases with ɛ _t . Clearly this is a relevant stabilizing effect in a range around X = 1 before log X destabilizes. As can be seen, the increase of Y destabilizes to the left of log X = 0 but stabilizes above this values. Thus, this parameter Y has an important influence in the middle range of log X. In the right sub-figure for D = 100, it is seen that now the curve for Y = 10 (dashed) is more stable than that of Y = 1 in the same middle range of log X. Notice that the bump now has a slightly smaller height than that in D = 0.1. The criticality curves follow the same destabilizing (ɛ _t = 0) and stabilizing (ɛ _t > 0) effect with respect to log X.

Figure 2:

Ra_c versus log X. Pr = 10. (a) F = 0.1, left: D = 0.1, right: D = 100. Y = 0.1 (solid), 1 (dashed), 10 (dotted). (b) F = 100. For D = 0.1 and Y = 0.1, 1, 10 (dashed). For D = 100: Y = 0.1 (solid), 1 (dotted), 10 (dash-dotted). ɛ _t = 0 (purely viscoelastic fluid [17]). ɛ _t > 0 viscoelastic Cattaneo–Christov fluid. Ra_c destabilizes with log X. However, for small F, Ra_c stabilizes in the middle range of log X.

When F is very large the viscoelastic fluid has more elastic properties than viscous and it is very easy to destabilize. This is shown in Figure 2(b) where F = 100. In this case all the curves have a very small Ra_c and show a destabilizing effect with respect to log X, the same as for the purely viscoelastic fluid (ɛ _t = 0). The criticality curves have a small difference with respect to ɛ _t . However it is clear that the system is more unstable for ɛ _t = 0.07, in contrast to the case of F = 0.1 presented in Figure 2(a). The purely viscoelastic fluid is the more stable one. With respect to the parameter D, it is shown that the dashed curves for D = 0.1 (and the corresponding three magnitudes of Y) are the more stable in the middle range of log X. Below follow the three curves of D = 100.

The behavior of the critical wavenumber corresponding to the curves of criticality of Figure 2(a) for Pr = 10 are presented in Figure 3(a). The left sub-figure for F = 0.1 and D = 0.1 shows that the critical wavenumber increases with log X when ɛ _t = 0 and ɛ _t = 0.07. However, it decreases when ɛ _t = 0.03 and 0.04298. Here it is found a bump in the middle range of log X for ɛ _t > 0 but not when ɛ _t = 0. Notice that the height of the bump increases with ɛ _t . The dotted curve for Y = 10 has always the smallest k _c in the middle range of log X, but that is not the case for larger log X. A similar situation is found in the right sub-figure for D = 100. However, in this case the bump in k _c has a smaller height than when D = 0.1.

Figure 3:

k _c versus log X. Pr = 10. (a) F = 0.1, left: D = 0.1, right: D = 100. Y = 0.1 (solid), 1 (dashed), 10 (dotted). (b) F = 100. For D = 0.1 and Y = 0.1, 1, 10 (dashed). For D = 100: Y = 0.1 (solid), 1 (dotted), 10 (dash-dotted). ɛ _t = 0 (purely viscoelastic fluid [17]). ɛ _t > 0 viscoelastic Cattaneo–Christov fluid. For small F, k _c decreases with log X up to a certain ɛ _t , then k _c increases with log X. A strong increase occurs in the middle range of log X. When F is large, k _c decreases with log X.

In Figure 3(b) the curves of k _c follow a very similar behavior as the curves of Ra_c in Figure 2(b). That means that the largest k _c corresponds to ɛ _t = 0 and D = 0.1 for the three Y’s. The smallest one corresponds to ɛ _t = 0.07, D = 100 and Y = 0.1 (solid).

The curves of the critical frequency of oscillation are given in Figure 4. The behavior of these curves for F = 0.1 in Figure 4(a) is in contrast to those of k _c. Here, in both the left and right sub-figures, ω _c increases with respect to log X for all ɛ _t ’s. In both sub-figures the lowest frequency in the middle range of log X is obtained when Y = 10 (dotted). However, it is interesting that to the left of this range the highest frequency is attained with Y = 1 (dashed). In the case for F = 100 of Figure 4(b), the critical frequency of oscillation decreases with respect to log X. However, it is found another important difference with respect to the results of Figure 3(b). Here, the curves of ω _c decrease with the increase of ɛ _t in the region from negative log X upto the middle range, but they increase with respect to ɛ _t in the region from the middle range to positive log X.

Figure 4:

ω _c versus log X. Pr = 10. (a) F = 0.1, left: D = 0.1, right: D = 100. Y = 0.1 (solid), 1 (dashed), 10 (dotted). (b) F = 100. For D = 0.1 and Y = 0.1, 1, 10 (dashed). For D = 100: Y = 0.1 (solid), 1 (dotted), 10 (dash-dotted). ɛ _t = 0 (purely viscoelastic fluid [17]). ɛ _t > 0 viscoelastic Cattaneo–Christov fluid. For small F, ω _c increases with log X, mainly in the middle range. However, for large F, ω _c decreases with log X.

The results of Ra_c for Pr = 50 are given in Figure 5. In Figure 5(a) the curves of F = 0.1 show an interesting behavior. First, clearly the curve of ɛ _t = 0 of the purely viscoelastic fluid remains far below the curves for ɛ _t > 0. That means that the relaxation time of the Cattaneo–Christov heat flux has a very strong stabilizing effect. Notice the cut made in the middle of the ordinates axis. Second, the curves destabilize with respect to log X for ɛ _t = 0.03 and 0.0412. However, when ɛ _t = 0.07 the curves stabilize when log X increases. This behavior occurs for both the left (D = 0.1) and right (D = 100) sub-figures. Third, the height of the stabilizing bump in the middle range of log X still increases with the increase of ɛ _t , despite the previous results. Fourth, the curves of Y = 10 (dotted) are the more unstable ones in the middle range of log X but for larger X, they are the more stable in a short range. In Figure 5(b) for F = 100, Ra_c decreases with log X which has a destabilizing effect. Moreover, the increase of ɛ _t also has a destabilizing influence on the system.

Figure 5:

Ra_c versus log X. Pr = 50. (a) F = 0.1, left: D = 0.1, right: D = 100. Y = 0.1 (solid), 1 (dashed), 10 (dotted). (b) F = 100. For D = 0.1 and Y = 0.1, 1, 10 (dashed). For D = 100: Y = 0.1 (solid), 1 (dotted), 10 (dash-dotted). ɛ _t = 0 (purely viscoelastic fluid [17]). ɛ _t > 0 viscoelastic Cattaneo–Christov fluid. For small F, Ra_c decreases with log X up to a certain ɛ _t , then Ra_c increases with log X. A strong increase occurs in the middle range of log X. When F is large, Ra_c decreases with log X.

The critical wavenumbers for Pr = 50 are given in Figure 6. Corresponding to the behavior of Ra_c, here too the behavior of k _c in Figure 6(a) for F = 0.1 is in contrast with that in Figure 3(a). In this case, the magnitude of the curve of ɛ _t = 0 is located far above those of ɛ _t > 0. Notice the cut made in the middle of the ordinates axis. The critical wavenumber decreases with respect log X for ɛ _t = 0.3 and 0.0412 but increases for ɛ _t = 0.07. This process makes the curve of ɛ _t = 0.0412 the lowest one to the right of the middle range of log X. This behavior is similar for both D = 0.1 and 100. As can be seen, the parameter Y = 10 (dotted) not always gives the lowest magnitude of k _c in the middle range of log X. An interesting phenomenon is found for ɛ _t = 0.07 for both D = 0.1 and 100. In the middle range of log X instead of a bump, it is found a fast increase in the magnitude of k _c. This does not happen for F = 100 in Figure 6(b). However, now a depression is found in the middle range of log X. Therefore, k _c decreases with log X, but that decrement is even more important in that particular range. As can be seen, the curves also decrease with the increase of ɛ _t .

Figure 6:

k _c versus log X. Pr = 50. (a) F = 0.1, left: D = 0.1, right: D = 100. Y = 0.1 (solid), 1 (dashed), 10 (dotted). (b) F = 100. For D = 0.1 and Y = 0.1, 1, 10 (dashed). For D = 100: Y = 0.1 (solid), 1 (dotted), 10 (dash-dotted). ɛ _t = 0 (purely viscoelastic fluid). ɛ _t > 0 viscoelastic Cattaneo–Christov fluid. For small F, k _c decreases with log X up to a certain ɛ _t , then k _c increases with log X. A strong increase occurs in the middle range of log X. When F is large, k _c decreases with log X.

In Figure 7 can be found the results of the critical frequency of oscillation. The case of F = 0.1 is found in Figure 7(a). In both sub-figures, for D = 0.1 and 100, it is found that ω _c increases with respect to log X. However, it is shown that ɛ _t produces a very important decrease in the critical frequency of oscillation. Notice the cut in the ordinates axis. The lower magnitude of ω _c is found for Y = 10 (dotted) in the middle range of log X. It is interesting that ω _c also shows a sharp increase when ɛ _t = 0.07, corresponding to that found for in both sub-figures of Figure 6(a). In the case of F = 100, Figure 7(b) shows that ω _c decreases with log X but all the curves present a smooth depression in the middle range of log X. The increase of ɛ _t decreases the magnitude of ω _c but smoothly in comparison with the decrease found in Figure 7(a). In the middle range of log X the higher ω _c is found for D = 100 and Y = 10, for all ɛ _t . The lowest one corresponds to D = 100 and Y = 0.1. The curve of D = 0.1 and Y = 0.1 to 10 goes in between.

Figure 7:

ω _c versus log X. Pr = 50. (a) F = 0.1, left: D = 0.1, right: D = 100. Y = 0.1 (solid), 1 (dashed), 10 (dotted). (b) F = 100. For D = 0.1 and Y = 0.1, 1, 10 (dashed). For D = 100: Y = 0.1 (solid), 1 (dotted), 10 (dash-dotted). ɛ _t = 0 (purely viscoelastic fluid [17]). ɛ _t > 0 viscoelastic Cattaneo–Christov fluid. For small F, ω _c increases with log X, mainly in the middle range. However, for large F, ω _c decreases with log X.

The awkward behavior of Ra_c with respect to ɛ _t in Figure 5(a) for Pr = 50 and F = 0.1 in comparison to that in Figure 2(a) for Pr = 10, was the motivation to check what happens for intermediate values of Pr. It was found that the behavior of Ra_c with respect to ɛ _t started to change at Pr = 25. This laborious search was an incentive to look for a more simple way to understand what might be going on with the curves of Ra_c at different Pr’s. It was found that the curves of Ra_c against ɛ _t shown in Figure 8 were an interesting way to appreciate that behavior. Two magnitudes of log X were selected as representative of both regions, to the left of the middle range log X = −3.5 (dashed) and to the right log X = 3.5 (solid). Only two Prandtl numbers, Pr = 10 and Pr = 50, were used. Notice that, to the left hand side of the minima, Ra_c decreases with respect to ɛ _t , after the minima Ra_c increases and after the maxima Ra_c decreases again. Observe that the dashed curve is below the solid one for both Pr = 10 and 50 in a wide range of ɛ _t . However, after a certain magnitude, to the right, this order is reversed. This means that the increase of log X stabilizes in the middle range of ɛ _t and destabilizes in the range located in the far right side of the curves. Note that the increase of log X also destabilizes in the region located to the left of the two minima. In the present paper ɛ _t was limited to the range 0 ≤ ɛ _t ≤ 0.07. However, Figure 8 shows that interesting behavior still could be found outside that range.

Figure 8:

Ra_c versus ɛ _t . Pr = 10 and 50. Graphical explanation of why the curves of criticality not always increase with respect to the Cattaneo–Christov relaxation time parameter ɛ _t and with respect to log X. Two representative magnitudes of log X have been selected: log X = −3.5 (dashed) and log X = 3.5 (solid). Notice how the pair of curves corresponding to each Pr intersect to each other twice changing the relative magnitude of Ra_c for each log X.

4 Conclusions

The natural convection of a viscoelastic Maxwell liquid with Cattaneo–Christov heat flux confined between two thick walls with finite thermal conductivity has been investigated. It is found the impossibility of having codimension-two points because the curves of criticality of oscillatory convection appear very far below those of stationary convection where the lowest critical Rayleigh number is 720, corresponding to the limit of very large X. Therefore, the research is restricted to oscillatory convection alone. The Prandtl numbers investigated were Pr = 5, 10, 25 and 50. In this paper only plots of the results corresponding to the magnitudes Pr = 10 and 50 are presented. It was found that with Pr = 5 the Cattaneo–Christov and Maxwell relaxation times were not able to stabilize the curves in such way to have a codimension-two point. The results are helpful to understand that the main influence of the Maxwell relaxation time is to destabilize the flow. However, the Cattaneo–Christov relaxation time always stabilizes the curves in the middle range of log X with a bump. The ratio of thermal diffusivities XY is relevant too. When Y increases the height of the bump in the middle range of log X decreases. The behavior of the curves of criticality of F = 0.1 and for Pr = 5 and 10 is as follows. The curves decrease with log X, then they show a bump and finally they decrease again to reach a limit when X is very large. Besides, these curves increase in all the range of log X with ɛ _t . However, the Cattaneo–Christov relaxation time also has an important influence on the convective stability of the viscoelastic liquid. It is found that for Pr = 50 the curves of criticality do not increase with ɛ _t as in Pr = 5 and 10. In this sense, the case Pr = 25 played an important role. It is found that at this Prandtl number the behavior of the criticality curves of Pr = 5 and 10 with respect to ɛ _t start to change into the behavior of the curves of Pr = 50 with respect to this parameter. That is, the dependence on ɛ _t started to reverse. This interesting behavior led to the idea of plotting the critical Rayleigh number against ɛ _t . It is found that, in fact, the critical Rayleigh number has a minimum against the Cattaneo–Christov relaxation time. In other words, Ra_c does not depend monotonically on ɛ _t when the liquid is viscoelastic with F = 0.1. It is also shown that these curves intersect to each other for different Prandtl numbers. However, for a large Maxwell viscoelastic relaxation time F = 100, it is reveled that Ra_c decreases very slightly, but monotonically, with respect to ɛ _t . Finally, it is clear that the presence of these two particular relaxation times modifies considerably the curves of pure viscoelastic oscillatory convection (corresponding to ɛ _t = 0).