Onset of Primary and Secondary Instabilities of Viscoelastic Fluids Saturating a Porous Layer Heated from below by a Constant Flux

We analyze the thermal convection thresholds and linear characteristics of the primary and secondary instabilities for viscoelastic fluids saturating a porous horizontal layer heated from below by a constant flux. The Galerkin method is used to solve the eigenvalue problem by taking into account the elasticity of the fluid, the ratio between the viscosity of the solvent and the total viscosity of the fluid and the lateral confinement of the medium. For the primary instability, we found out that depending on the rheological parameters, two types of convective structures may appear when the basic conductive solution loses its stability: stationary long wavelength instability as for Newtonian fluids and oscillatory convection. The effect of the lateral confinement of the porous medium by adiabatic walls is to stabilize the oblique and longitudinal rolls and therefore selects transverse rolls at the onset of convection. In the range of the rheological parameters where stationary long wave instability develops first, we use a parallel flow approximation to determine analytically the velocity and temperature fields associated with the monocellular convective flow. The linear stability analysis of the monocellular flow is performed, and the critical conditions above which the flow becomes unstable are determined. The combined influence of the viscoelastic parameters and the lateral confinement on the characteristics of the secondary instability is quantified. The major new findings concerning the secondary instabilities may be summarized as follows: (i) For concentrated viscoelastic fluids, computations showed that the most amplified mode of convection corresponds to oscillatory transverse rolls, which appears via a Hopf bifurcation. This pattern selection is independent of both the fluid elasticity and the lateral confinement of the porous medium; (ii) For diluted viscoelastic fluids, the preferred mode of convection is found to be oscillatory transverse rolls for a very laterally-confined medium. Otherwise, stationary or oscillatory longitudinal rolls may develop depending on the fluid elasticity. Results also showed the destabilizing effect of the relaxation fluid elasticity and the stabilizing effect of the viscosity ratio for the onset of both primary and secondary instabilities.


Introduction
The study of viscoelastic fluids have applications in a number of processes that occur in industry, such as the extrusion of polymer fluids, solidification of liquid crystals, suspension solutions and petroleum activities.In contrast to the case of Newtonian fluids, study of thermal convection of viscoelastic fluids in porous media is limited.In rheology, one crucial problem is the formulation of the constitutive equations regarding viscoelastic fluid flows in porous media.Recently, a modified Darcy's law was employed to study the stability of a viscoelastic fluid in a horizontal porous layer using linear and nonlinear stability theory ([1]- [9]).Kim et al. [1] and Yoon et al. [2] performed a linear stability analysis and showed that in viscoelastic fluids such as polymeric liquids, a Hopf bifurcation as well as a stationary bifurcation may occur depending on the magnitude of the viscoelastic parameter.From the nonlinear point of view, Kim et al. [1] carried out a nonlinear stability analysis by assuming a densely packed porous layer and found that both stationary and Hopf bifurcations are supercritical relative to the critical heating rate.The question of whether standing or traveling waves are preferred at onset has been fully addressed by Hirata et al. [4].
The three-dimensional convective and absolute instabilities of a viscoelastic fluid in presence of a horizontal pressure gradient have been analyzed by Hirata and Ouarzazi [5].Alves et al. [6] studied the effect of viscous dissipation of viscoelastic fluids at the onset of convection.In addition to its theoretical interest, Delenda et al [7] have showed that viscoelastic convection in porous media may be useful for industrial applications interested by the separation of species of viscoelastic solutions.
The introduction of a porous packing allows to control the average vertical convective velocity and to generate a homogeneous convection current, improving the separation of species.Fu et al. [8] performed direct numerical simulations on two-dimensional thermal convection of a viscoelastic fluid saturating a porous square cavity.Their numerical experiments revealed the existence of a second transition from oscillatory convection to stationary one followed by a third transition to oscillatory convection for some combinations of rheological parameters while these successive transitions never occur for other combinations of viscoelastic parameters.Taleb et al. [9] used both theoretical and numerical approaches and obtained a global picture and bifurcations diagrams on possible successive bifurcations of convection patterns in a square porous cavity saturated by a viscoelastic fluid.
All the above investigations considered conventional boundary conditions, namely impermeable isothermal horizontal plates and impermeable adiabatic side walls, commonly known as Horton-Rogers-Lapwood convection.However, to the best of our knowledge, no results have been published for thermal convection of viscoelastic fluids when the porous medium is heated from below and cooled from above with a constant heat flux.Therefore, the objective of this work is to fill this gap by investigating the onset of three-dimensional primary and secondary instabilities of a viscoelastic fluids under the assumption that the upper and lower horizontal walls are impermeable and are kept at a constant flux, while the lateral vertical walls are considered impermeable and adiabatic.
For Newtonian fluids, the stability of an infinite porous layer with different boundary conditions was studied by Nield [10] and is well documented in Sect.6.2 of the book by Nield and Bejan [11].For the case of a porous medium heated from the bottom and cooled from the top by a constant heat flux, Nield [10] found that the critical Rayleigh number at the onset of convection is approximately 12 with a vanishing critical wavenumber.Mamou et al. [12] extended the work of Nield [10] by taking into account the effect of the anisotropy of the porous medium.Mojtabi and Rees [13] studied the case where the impermeable boundary walls have a finite thickness.They analyzed the combined influence on the onset of convection of the ratio between the thermal conductivity of the horizontal walls and the thermal conductivity of the porous medium as well as the ratio between the thickness of the horizontal walls and the thickness of the porous layer.
Kimura et al. [14] investigated secondary instabilities for a Newtonian fluid saturating a porous medium heated from below by a constant flux.For Rayleigh number larger that its critical value 12 above which the conduction state looses its stability against long wave instability, these authors used the parallel flow approximation and obtained a nonlinear solution which corresponds to a monocellular flow.Two-dimensional numerical results were also presented to test the validity of the approximated nonlinear solution.In addition, they analyzed its stability against three dimensional disturbances and showed that the monocellular flow is linearly stable to transverse disturbances for This contribution aims at understanding how the viscoelastic character of the fluid influences the properties of convection at the onset of primary and secondary instabilities, when the porous layer is heated from below by a constant flux.Therefore, this work may be viewed as an extension to viscoelastic fluids of the work done by Kimura et al. [14].
The paper is organized as follows.After presenting the governing equations in section 2, the stability of the conductive state is studied in section 3 by considering steady as well as oscillatory three-dimensional perturbations.Section 4 is devoted to the discussion of the combined effects of the viscoelastic parameters and the lateral aspect ratio of the porous medium on the pattern selection at the onset of secondary instabilities.Finally, in section 5, the main conclusions of the present study are presented.

Mathematical formulation
Let us consider an isotropic and homogeneous porous cavity of thickness e, height H, width W (see figure 1).The porous medium is saturated by an Oldroyd-B fluid and we assume that the solid matrix is in local thermal equilibrium with the fluid.The upper and lower horizontal walls are kept at constant flux, while the lateral vertical walls are considered adiabatic.The solid walls of the domain We assume that the Oberbeck-Boussinesq approximation holds.There are several ways to obtain macroscopic laws for polymeric flows in a porous medium: by direct numerical simulations of viscoelastic flows in a specific pore geometry model (a good review of these studies can be found in [15]) or analytical ways.In general, the former is the most commonly used way for the derivation of macroscopic laws.It can be divided in two techniques: the REV method (representative elementary volume method) and the homogenization theory.The starting point for the two techniques is a local description in a pore scale.The pore space is assumed to be saturated by an incompressible viscoelastic fluid.For slow flows, the momentum balance equation can be linearized: where U * is the fluid velocity field, p * is the pressure, τ is the stress tensor and g is the gravity field.
In Newtonian incompressible fluids, the constitutive relation between stress tensor τ and strain The rheological model relating τ and D for viscoelastic fluids, such as a polymeric solution composed of a Newtonian solvent and a polymeric solute of "Newtonian" viscosity µ s and "elastic viscosity" µ p [16] respectively, is given by: with and where λ * 1 represents the relaxation time.Then, by combining 2 -4 we obtain the constitutive equation: where the dynamic viscosity µ and the retardation time λ * 2 are related to µ s and µ p by: An Oldroyd-B fluid may thus be characterized by three parameters: the dynamic viscosity µ, the relaxation λ * 1 and the retardation λ * 2 times.The relation Γ = λ * 2 /λ * 1 may also be used instead of λ * 2 .
In order to derive a macroscopic filtration law based on the Oldroyd constitutive equation, we have to introduce the filtration velocity V * defined by the Dupuit's equation : where φ is the porosity.Substituting Equation 5 into 1 and using the REV method by averaging the resulting equation and taking into account Equation 6 leads to: where K is the permeability.
Under the assumption of low Reynolds number based on the pore dimension, the generalized Darcy's law 7 is also derived by [17] using a homogenization theory.
The fluid density ρ obeys the state law : where ρ 0 is the fluid density at temperature T * 0 which is chosen here as the temperature at the geometric center of the cavity, and β T is the thermal expansion coefficient.Energy and continuity equations can then be written as : The boundary conditions at the impermeable horizontal walls kept at a constant flux q and the impermeable insulated vertical walls are: here, (ρc), µ, ν, k T , α = k T /(ρc) f are respectively the heat capacity per unit volume, the dynamic and kinematic viscosity of the fluid the effective thermal conductivity and the effective thermal diffusivity.Subscript (sf) refers to an effective quantity, while (f) refers to the fluid alone.
We choose ) and qH /k T as reference quantities for lengh, velocity, time, pressure and temperature difference (T * − T * 0 ).With this scaling, the following set of dimensionless equations is obtained: The dimensionless boundary conditions are: The Darcy-Prandtl number P r D is defined as P r D = (φP r)/Da, with Da = K/H 2 and P r = ν/k T .Since in the common porous media the Darcy number is very small, the Darcy-Prandtl number P r D takes quite large values.Therefore, the first term in Equation 13is omitted in what follows.The remaining dimensionless parameters are : the filtration Rayleigh number the horizontal and lateral aspect ratios the relaxation time and the ratio Γ that varies in the interval [0, 1] This model reduces to the Maxwell model in the limit Γ → 0 and to the Newtonian model in the limit Γ → 1.
In the following we will examine the stability of the conductive state (the primary instability) as well as the stability of the monocellular flow (the secondary instability).

Primary stationary and oscillatory instabilities
In the conductive regime, the basic solution is a motionless state V = 0 with a vertical thermal The aim of this section is to perform a temporal stability analysis of the conductive state with respect to both stationary and oscillatory disturbances.

Infinite aspect ratios
To investigate the stability of the basic solution, infinitesimal three-dimensional perturbations are super-imposed onto the basic solution: We first assume very large aspect ratios A(A → ∞) and a(a → ∞).The three-dimensional disturbance quantities are expressed as where k and l are the wave numbers in the x and y directions respectively, and the temporal growth rate of unstable perturbations is given by the imaginary part of the complex frequency ω = ω r +iω i .Therefore, the neutral temporal stability curve is obtained for ω i = 0 which selects dominant modes at the onset of convection.
Substituting Equations ( 20)-( 21) into ( 12)-( 15), linearizing the equations and applying the curl twice to the momentum balance equation, one can obtain where D = d dz and k2 = k 2 + l 2 .The corresponding boundary conditions take the form The system (33) -( 34) is solved by means of the Galerkin method using the following expansions  The number M of modes is chosen so that the quantitative convergence is secured.
As the viscoelastic parameters appear only in front of a time derivative in the momentum equation ( 16), the elasticity of the fluid cannot influence the properties of a stationary instability.
Consequently, the characteristics of the stationary instability are the same as for Newtonian fluids.
For such fluids, linear instability analysis has been considered by Nield [10] and has provided quantitative information on the stability condition when the porous layer is supposed infinite in x and y directions.
We first consider perturbations in the form of stationary convection.Having used the Galerkin expansion (25) -( 26) with M = 5, we obtain results with a very good agreement with those obtained in [10].It is well established that for isothermal horizontal boundaries, competition between the processes of stress relaxation, strain retardation and thermal diffusion may also lead to an oscillatory convective instability as a first bifurcation ([1]- [9]).This feature is also found in the actual study when the viscoelastic fluid saturating the porous layer is heated by a constant flux.
In Figure 2

effect of lateral confinement on pattern selection
This section is devoted to investigate the effect of the lateral confinement of the porous cavity by assuming a very large aspect ratio A(A → ∞) and finite lateral aspect ratio a.The three-dimensional disturbance quantities respecting the boundary conditions 15 are expressed as (u, w, θ, p) = ũ(z), w(z), θ(z), p(z) exp(ikx − iωt)cos(Lπy/a) The governing equations are still the system (33) -(34), except with l now replaced by Lπ/a where the integer L is the number of rolls in the y direction.
We begin the study by considering the stability of the conductive state against stationary rolls with axes parallel to the x direction, called longitudinal rolls (LRs).Steady longitudinal rolls are characterized by k = 0, L 0 and ω r = 0 .The dependence of the critical Rayleigh number at the onset of (LRs) on the lateral aspect ratio a for different number L of rolls is displayed in Figure 4(a).
For comparison we also represent in the same figure the threshold of the steady long wave instability.
The threshold of steady three-dimensional patterns in the form of oblique rolls (i.e.k 0, L 0 and ω r = 0) are bounded by the thresholds of the two limiting cases: the steady long wave instability and steady LRs.
We remark that the mode L = 1 is the most unstable mode for LRs.As it is expected, we note that the critical Rayleigh number increases as a decreases, meaning that the lateral confinement stabilizes the conductive state against longitudinal rolls.We also note that as a → ∞, the limiting value of Ra = 12 is reached monotonically and an infinity of modes may be simultaneously unstable of a, indicating that the true critical Rayleigh number strongly depends on both a and L for fixed rheological parameters.The destabilizing oscillatory LRs mode changes in the intersection points of neural curves from a mono-cellular flow to a two-cellular flow and so on, as the lateral aspect ratio a increases.In addition, the behavior of the critical Rayleigh number is non-monotonic as a increases.We also note that the maximum of critical Rayleigh number decreases as a increases and tends asymptotically to the critical threshold found in the unbounded case (a → ∞).The results are therefore in contrast to the case of stationary LRs where the dominant mode corresponds to L = 1 independently of the lateral confinement.
In Figure 4(b), the critical Rayleigh number at the onset of oscillatory TRs is indicated by the horizontal line.As can be seen from this figure, finite values of a stabilize oscillatory LRs and may select oscillatory TRs as a dominant mode of convection.

Nonlinear solution and formulation of its linear stability
According to above linear stability analysis, we found that a stationary bifurcation occurs giving rise to a convective pattern in the form of a long wave instability in the x direction provided that the elasticity number λ 1 do not exceed a particular value λ f 1 which depends on the viscosity ratio Γ .In that case, the viscoelastic fluid behaves like a Newtonian fluid.Consequently, the nonlinear solution in the regime of steady long wave convection is the same regardless of weather or not the fluid is viscoelastic.
As shown by Bejan [18] for a vertical cavity, and later adopted by Vasseur et al. [19] and Sen et al. [20] for inclined cases, one may assume the existence of a two-dimensional and fully developed counterflow.This may be a good approximation for the mid-region of the horizontally extended space on condition that the unicellular convection is stable.By assuming a shallow cavity A 1 and by using the parallel flow approximation, Kimura et al. [14] found that the analytical solution for the monocellular flow consists of: a horizontal asymmetric velocity with a zero mean along any vertical section, and a vertical as well as a horizontal stratification of the temperature, with and where C is negative or positive according to whether the flow is clockwise or counter-clockwise and both solutions are possible depending on the initial conditions.
From Equation (32) it is seen that no motion may be induced inside the cavity for Ra < 12.For the case of a porous medium heated from the bottom and cooled from the top by a constant heat flux, a critical Rayleigh number of Ra = 12 for the onset of convection was predicted by Nield [10].This result is in agreement with the prediction of Equation (32).
For finite aspect ratio, Kimura et al. [14] performed two dimensional numerical simulations of the full problem.Their numerical results show that the conductive state is stable when the Rayleigh number is smaller than 12. Computations carried out for Ra in excess of 12 were found to agree with analytical solutions (29 -31).
The equations governing the linear stability of the monocellular flow are obtained by the same previous approach used for the stability of the conductive basic solution.By assuming very large aspect ratios A(A → ∞) and a(a → ∞) the following system is obtained where we substitute U 0 and T 0 by their explicit expressions (29) -(31).
The corresponding boundary conditions take the form On the other hand, if we assume a very large aspect ratio A and a finite value of the lateral aspect ratio a, the governing equations are still the system (33) -(34) where k2 is replaced by The resulting linear stability problem is solved by means of the Galerkin method, using the expansion (25) -( 26) at the order M = 30.

Results for Newtonian fluids
To verify the accuracy of our numerical results based on the Galerkin expansion to the order M = 30, we perform a test for the limiting case of a Newtonian fluid and compare the results with those obtained by Kimura et al. [14] .In the first instance, two-dimensional disturbances, corresponding to l = 0, were considered.We found out that for the Newtonian fluid, the base velocity and temperature On the other hand, Kimura et al. [14] considered three-dimensional disturbances with the value of the y-wave number l being gradually increased from zero.For l > 0, the stability analysis indicates that the monocellular flow will be destabilized not by a Hopf bifurcation, but by an exchange of stability for which the x-wave number k vanishes.In that case the threshold of the appearance of steady longitudinal rolls as a secondary instability is found to be R L c2,s ≈ 311.53.Since this critical Rayleigh number is much lower than any of those for the Hopf bifurcations obtained when k 0, Kimura et al. [14] concluded that the monocellular flow will in fact be destabilized by longitudinal, rather than transverse, disturbances.
In the second instance, three-dimensional disturbances, corresponding to k 0 and l 0, were considered in this study.Numerical results performed by assuming infinite aspects ratios A and a indicated that the most unstable mode corresponds to k = 0 and l 0. The corresponding critical Rayleigh number above which this most unstable mode in the form of steady LRs is R L c2,s = 313.107.
Once again, this critical value agrees very well with R L c2,s ≈ 311.53 obtained in [14].
In the third instance, the effect of the confinement of the porous medium in y direction is explored.
We plot in Figure 5(b) the critical Rayleigh number against the aspect ratio for several of the leading modes, from which it is clear that (L = 1) remains the destabilizing mode, ahead of the other modes (L > 1), and that the order of these modes, in the sense that Ra c (L) < Ra c (L + 1), is preserved as a increases.In particular, we also note that as a → ∞, the limiting value of R L c2,s = 313.107 is approached monotonically.We note in this figure that the minimum of neutral stability curves increases when Γ is augmented to reach the critical value for Newtonian fluids in the limit of Γ = 1.Physically, this result means that concentrated polymeric solutions with a small viscosity ratio Γ favor the appearance of oscillatory multicellular flow convection as a secondary instability.On the other hand, for diluted viscoelastic solutions, more heating is needed to trigger the secondary instability.We report in Table 1 the computed results of critical Rayleigh number Ra T c2 , critical frequency ω T c2 and critical wave number k T c2 at the onset of the secondary instability organized as oscillatory TRs for λ 1 = 0.1 and different values of Γ .Table 1 shows a strong stabilizing effect of the viscosity ratio.The values of the critical oscillatory frequency decrease with decreasing Γ .This implies that emerging transversal convection rolls have a larger time-period and move with a larger phase velocity when the polymer concentration is high.
We now present results corresponding to the influence of the fluid elasticity by considering the properties of the emerging oscillatory TRs at different values of λ 1 for a fixed value of Γ .We conclude that the preponderant effect on the properties of the emerging oscillatory TRs is mainly linked to the variations in the viscosity ratio, while the effect of the elasticity remains very weak.Finally, we present in the second part of this section the secondary instability results in the case where disturbances are assumed in the form of longitudinal rolls (LRs).We mention that as for the primary instability, the onset of stationary LRs convection is not affected the two viscoelastic parameters.Consequently, the critical Rayleigh number above which stationary LRs convection develops as a secondary instability is the same as that found for Newtonian fluids, namely R L c2,s = 313.107.However, the computations indicate Hopf bifurcation from steady unicellular flow to oscillatory LRs convection.We emphasize that the Hopf bifurcation to oscillatory LRs is not observed for Newtonian fluids and is due solely to the viscoelastic character of the fluids.The effects of the two viscoelastic parameters on the linear properties of the oscillatory LRs convection are examined in the remainder of this subsection.In order to evaluate the effect of elasticity alone,  In contrast, for the combination (λ 1 = 0.5, Γ = 0.75), the least stable mode of convection changes from steady LRs to oscillatory LRs, meaning that elastic effects become the most important ones in this range.In the same way, the preferred pattern as a secondary instability depends on the viscosity ratio Γ .Table 3 shows that by keeping λ 1 = 0.1 and increasing gradually Γ from Γ = 0.75 (diluted solutions) to Γ = 0.3 (concentrated solutions), the most amplified mode of convection evolves from steady LRs to oscillatory LRs and eventually to oscillatory TRs.
All the results stated in the subsection 4.3 are obtained by assuming infinite aspects ratios in x  and y directions.For the sake of brevity, we exemplify the effect of the lateral aspect ratio a on the pattern selection for two combinations of rheological parameters (Γ = 0.75, λ 1 = 0.3) and (Γ = 0.5, λ 1 = 0.1).We plot in Figures 7(a It is important to note that Figure 7(a) also shows that the curve corresponding to oscillatory longitudinal mode with L = 1 intersects the line representing the critical Rayleigh number of oscillatory TRs Ra T c2 = 355.03at a particular value of the lateral aspect ratio a = a * ≈ 0.4.This means that perturbations promote the appearance of oscillatory TRs if a < a * , oscillatory LRs or a steady monocellular LRs if a * < a < a * * and stationary LRs if a > a * * .In the case of the combination (Γ = 0.5, λ 1 = 0.1), this behavior is not observed since as it can be seen from Figure 7(b), the critical Rayleigh number of oscillatory TRs is much smaller than the critical Rayleigh number for both stationary and oscillatory LRs.For this particular combination, the system selects oscillatory TRs independently of the lateral confinement.

Conclusion
In the present paper, Galerkin method is used to investigate the primary and secondary instabilities of viscoelastic fluids saturating a porous layer heated from below by a constant flux.
The modified Darcy's law based on the Oldroyd-B model was used for modeling the momentum equation.In addition to Darcy-Rayleigh number Ra, two viscoelastic parameters play a key role when characterizing the temporal behavior of the instability, namely, the relaxation time λ 1 which measures the elasticity of the fluid and the ratio Γ between the viscosity of the solvent and the total viscosity of the fluid.In the first part of the paper, three-dimensional disturbances were considered in order to study the stability of the basic motionless solution.For sufficiently elastic fluids, we found that the primary instability is oscillatory.Otherwise, the primary bifurcation gives rise to stationary long wave instability.Results indicated that the lateral confinement of the porous layer by isolated side walls eliminates oblique or longitudinal rolls in favor of two-dimensional transverse rolls.Based on a fully developed parallel flow assumption, a nonlinear analytical solution for the velocity and temperature fields was developed in the range of the rheological parameters where stationary long wave instability develops first.In the second part of the paper, we reported findings on the linear stability analysis of the monocellular flow which is performed with special attention given to the interplay between the viscoelastic parameters and the lateral aspect ratio a of the porous layer.For weakly elastic fluids we determined a second critical value of Rayleigh number above

Figure 1 .
Figure 1.The porous rectangular cavity heated from below by a constant flux.
Fig 2(a)  represents the marginal stability curve in the ( k, Ra) plane and shows that a long wave instability (i.e. the critical wave number kc = 0) may develop if the Darcy-Rayleigh number exceeds the critical value Ra s = 12, 009 in accordance with the critical value Ra s = 12 found in[10].

Γ = 0. 5 Figure 3 .
Fig. 3(a) and in Fig. 3(b) respectively.It is clear from 3(a) that an increase in λ 1 leads to flow destabilization, i.e. to a reduction in the respective critical Rayleigh number.Fig. 3(a) also shows the stabilizing effect of the ratio Γ .Moreover, as it is seen in fig.3(a), for a fixed value of Γ , there exists a particular value of λ 1 = λ f 1 where the critical Rayleigh numbers for the onsets of both oscillatory and stationary convection coincide and therefore a codimension two bifurcation occurs.For λ 1 > λ f 1 , Fig. 3(b) shows that the critical frequency decreases with the decrease of the fluid elasticity or the increase of the viscosity ratio.

Figure 4 .
Figure 4. Critical Rayleigh number against the lateral aspect ratio with different number of rolls (L = 1: red dashed curve , L = 2: green dotted curve, L = 3: black dash dotted curve and L = 4: blue densely dotted curve): (a) steady longitudinal rolls; (b) oscillatory longitudinal rolls for Γ = 0.1 and λ 1 = 0.5.In both figures, the horizontal lines indicate the corresponding critical Rayleigh number for transverse rolls.

PreprintsFigure 5 .
Figure 5. Newtonian fluids: (a) Neutral stability curve at the onset of oscillatory transverse rolls; (b) Critical Rayleigh number at the onset of steady longitudinal rolls against the lateral aspect ratio with different number of rolls (L = 1: red dashed curve , L = 2: green dotted curve, L = 3: black dash dotted curve and L = 4: blue densely dotted curve).The horizontal line corresponds to the threshold of oscillatory transverse rolls.

Figure 5 (
b) also shows that the curve corresponding to steady longitudinal mode with L = 1 intersects the line Ra T c2 = 506.07 at a particular value of the lateral aspect ratio a = a * .This means that perturbations promote the appearance of oscillatory TRs provided that a < a * or stationary LRs otherwise.Preprints (www.preprints.org)| NOT PEER-REVIEWED | Posted: 4 April 2017 doi:10.20944/preprints201704.0013.v1Peer-reviewed version available at Fluids 2017, 2, 42; doi:10.3390/fluids20300424.3.Results for viscoelastic fluids 4.3.1.Hopf bifurcation to transverse rolls In order to study the influence of viscoelastic parameters on the secondary instability, we first computed the bifurcation line from a stationary monocellular convective pattern to oscillatory TRs (l = 0) for either a fixed value of the elasticity number λ 1 with varying values of the viscosity ratio Γ or a fixed value of Γ with varying values of λ 1 .With regard to the question of the influence of the viscosity ratio Γ for a viscoelastic fluid with a relaxation time λ 1 = 0.1 on the onset of a secondary instability in the form of oscillatory TRs, Figure 6(a) illustrates the behavior of neutral stability curves in the (k, Ra) plane for Γ = 0.75, Γ = 0.5 and Γ = 0.3.For a comparison, the Newtonian case (Γ = 1) is also represented on Figure 6(a).

Figure 6 (
b) presents neutral stability curves for Γ = 0.75, a typical viscosity ratio value for Boger fluids and different values of the relaxation time λ 1 = 0.1, λ 1 = 0.35 and λ 1 = 0.5.We note from this figure that the neural stability curves are nearly superposed when λ 1 is increased , meaning that beyond λ 1 = 0.1, the increase in the fluid elasticity has a little influence on the critical Rayleigh number at the onset of oscillatory TRs.

3 and λ 1 =
0.5 cases are investigated for a fixed Γ = 0.75.On the other hand, the effect of viscosity ratio alone is studied by fixing λ 1 = 0.1 and investigating the Γ = 0.75, Γ = 0.6, Γ = 0.5 and Γ = 0.3 cases.The computed results for the six different cases are reported in ) and 7(b) the variation of the critical Rayleigh number for both stationary LRs and oscillatory LRs as a function of the lateral aspect ratio a in the cases (Γ = 0.75, λ 1 = 0.3) and (Γ = 0.5, λ 1 = 0.1) respectively.Computations showed that there is a competition between the two patterns in the sense that depending of the magnitude of lateral confinement, the system may select either stationary LRs or oscillatory LRs.For fixed value of L and by increasing a, the following behavior is observed for the curves representing the critical Rayleigh number for LRs and all values of rheological parameters, (see Figure7(a) and 7(b)): i) the curve associated to the critical Rayleigh number of oscillatory LRs decreases to reach a minimum equal to its value for infinite a.This minimum point moves to the right in the (a, Ra) plane when the number of rolls L is increased; ii) then, the same curve increases to intersect an other branch corresponding to the critical Rayleigh number of steady LRs at a particular value of a; iii) finally, when a exceeds this particular value, the curve associated to the critical Rayleigh number of steady LRs becomes the lower curve, decreases monotonically and tends asymptotically to the critical Rayleigh number R L c2,s = 313.107 of steady LRs found in the case of infinite a.For the particular combination (Γ = 0.75, λ 1 = 0.3), as in the infinite limit of a, the critical Rayleigh number R L c2,s = 313.107 of steady LRs is less than the critical Rayleigh number R L c2,osc = 317.55 of oscillatory LRs, the decreasing curve of the critical Rayleigh number of steady LRs with L = 1 crosses the absolute minimum R L c2,osc = 317.55 of oscillatory LRs at a critical value a = a * * (a * * ≈ 2 in Figure7(a)).Consequently, for all values of a larger than a * * , the dominant mode of convection is a steady monocellular LRs .Otherwise, the system may select oscillatory LRs or a steady monocellular LRs Preprints (www.preprints.org)| NOT PEER-REVIEWED | Posted: 4 April 2017 doi:10.20944/preprints201704.0013.v1Peer-reviewed version available at Fluids 2017, 2, 42; doi:10.3390/fluids2030042depending on a.

Figure 7 .
Figure 7. Critical Rayleigh number for the onset of steady and oscillatory longitudinal rolls as a function of aspect ratio a for different number L of rolls (L = 1: red dashed curve , L = 2: green dotted curve, L = 3: black dash-dotted curve, L = 4: blue densely-dotted curve).(a) Γ = 0.75 and λ 1 = 0.3; (b) Γ = 0.5 and λ 1 = 0.1.The horizontal line corresponds to the threshold of oscillatory transverse rolls

Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 4 April 2017 doi:10.20944/preprints201704.0013.v1 Peer
Rayleigh number as high as 506, at which point a Hopf bifurcation sets in.However, further analysis indicated that an exchange of stability due to longitudinal disturbances will occur much sooner at Rayleigh number equal to 311.53.

Table 2
gathers the results for seven values of λ 1 .It can be observed from Table 2 that critical Rayleigh number Ra T c2 , critical frequency ω T c2 and critical wave number k T c2 at the onset of the secondary instability are weakly dependent on the elasticity number number λ 1 .

Table 1 .
Critical Rayleigh number Ra T c2 , frequency ω T c2 and wave number k T c2 at the onset of moving transverse rolls as a secondary instability for λ 1 = 0.1 and different values of Γ

Table 2 .
Critical Rayleigh number Ra T c2 , frequency ω T c2 and wave number k T c2 at the onset of moving transverse rolls as a secondary instability for Γ = 0.75 and different values of λ 1

Table 3 ,
which indicates the critical Rayleigh number, wave number and oscillatory frequency at the onset of oscillatory LRs secondary instability.As has already been highlighted in the previous sections considering the primary instability and the TRs secondary instability, we recognize the destabilizing effect of the elasticity number λ 1 and the destabilizing effect of the viscosity ratio Γ .Moreover, a comparison between Tables 1, 2 and 3 attests that the frequencies of oscillatory LRs are much smaller An additional remark about Table3is necessary.For comparison purposes, we also indicate in this table, the threshold of both stationary LRs and oscillatory TRs.It is clear that the true critical Rayleigh number depends on the combination of the rheological parameters.The least stable mode of convection is the one with smallest critical Rayleigh number and is identified in Table3with a bold character.For instance, we consider diluted viscoelastic solutions with Γ = 0.75 with different elasticity number λ 1 .For the combination of the rheological parameters (λ 1 = 0.1, Γ = 0.75), the true critical Rayleigh number is R L c2,s indicating that the secondary instability pattern is in the form of steady LRs.In that case, polymeric solutions are almost inelastic and evolve as a Newtonian fluid.

Table 3 .
Critical Rayleigh number Ra L c2 , frequency ω L c2 and wave number k L c2 at the onset of oscillatory longitudinal rolls as secondary instability for different values of Γ and λ 1 s