1 Introduction

Relative to the research activity on thermal convection, the work on thermosolutal convection is quite limited. It is now well known that the two mutually opposing diffusing mechanisms give rise to a wide range of interesting phenomena. In such systems with opposing influence on instability, the existence of a net negative density distribution does not ensure stability. Salt fingering (stationary) and diffusive instabilities (overstable) are of importance in these systems. A linear stability analysis of the two-component system was first performed by Nield and Bejan ^[1].

Most of the literature on double diffusive convection, however, did not consider coupled molecular diffusion, namely, the cross-diffusion effects of Soret and Dufour. These coupled fluxes are due to irreversible thermodynamic effects. If one of the two properties is temperature, the diffusion of solute due to an applied temperature gradient is called the Soret effect or the thermo-diffusion effect, and the diffusion coupled with it is the Dufour effect or diffusion-thermo effect. In liquids, the Dufour effect is negligibly small.

In the presence of chemically inert micron-sized particles, the suspension can be modelled by Stokes' couple stress liquid whose governing equations can be obtained from Eringen's equations for a micropolar continuum. Siddheshwar and Pranesh ^[2] were the first to study thermoconvective instabilities in such liquids with the single-buoyancy effect. Malashetty et al. ^[3] extended this study to a system with a second diffusing component and the Soret effect.

Rudraiah and Siddheshwar ^[4] considered both Soret and Dufour effects on two-component convection of Newtonian liquids saturating a sparsely-packed porous medium. A very comprehensive account of the cross-diffusion coefficient on finite-amplitude instability and heat and mass transport was presented in the paper. There have been other works that consider thermoconvective instability in the couple-stress liquid with the single-buoyancy effect ^[5-12].

Thermo-diffusion (Soret effect) and diffusion-thermo (Dufour effect) effects on convective instabilities in nanofluids have been theoretically investigated by Kim et al. ^[13]. They noticed that both the Soret and Dufour effects make nanofluids unstable, and the heat transfer enhancement by the Soret effect in binary nanofluids is more significant than that in normal nanofluids. Siddheshwar et al. ^[14] made a detailed study of the effect of Soret coefficient on heat transport in twenty nanoliquids. Wang and Tan ^[15] investigated stability analysis of double-diffusive convection in a viscoelastic fluid with the Soret effect occupying a porous medium using a modified-Maxwell-Darcy model. They showed that the relaxation time also enhances the instability of the system along with the Soret parameter. Altawallbeh et al. ^[16] analytically studied double-diffusive convection in an anisotropic porous layer heated and salted from below with an internal heat source and Soret effect using both linear and nonlinear stability analyses. They found that increasing the mechanical anisotropy parameter, the Soret parameter, and the internal Rayleigh number enhances heat and mass transfer. The instability of a horizontal layer of a binary nanofluid in a vertical magnetic field was investigated by Gupta et al. ^[17] using the normal mode analysis and weighted residual method. Complex expressions for the Rayleigh number were simplified by valid approximations for an analytical study, and numerical investigations were made for the alumina-water nanofluid. They showed that the critical wave number increases when the Chandrasekhar number increases, and it is independent of the solute and nanoparticles. Agarwal and Rana ^[18] studied the onset of periodic and aperiodic convection in a binary nanofluid saturated rotating porous layer. They obtained the Rayleigh numbers for stationary and oscillatory convections in terms of various non-dimensional parameters. The effect of the physical parameters on the aperiodic convection was analyzed graphically. The thermal instability using the linear stability analysis in a rotating anisotropic porous medium saturated by a nanofluid was analytically studied by Agarwal et al. ^[19]. The expression of the Rayleigh number for both stationary and oscillatory convections in the case of the bottom-heavy arrangement was obtained. A reversed trend between the stationary and oscillatory modes for the bottom-heavy and top-heavy arrangements was explained in the paper. Ren and Chan ^[20] studied double-diffusive convection in a vertical cavity with horizontal temperature and concentration gradients using the lattice Boltzmann method. They showed that an increase in the buoyancy ratio from 0.01 to 2 causes a decrease in the average Nusselt and Sherwood numbers. However, a further increase in the buoyancy ratio from 2 causes an increase. In addition, the double-diffusive flow was observed to be unsteady at a small Prandtl number (Pr=0.1) and large Lewis numbers (Le > 6). Hu et al. ^[21] studied transient growth due to non-normality for the Poiseuille-Rayleigh-Bénard flows of binary fluids with the Soret effect. They considered two cases (ⅰ) negative separation factors and (ⅱ) positive separation factors. They showed that for negative separation factors, the transient growth is strong, whereas for positive separation factors, it is weak. They used the least-stable mode for computation. Other noteworthy works considering Soret and/or Dufour effects in their problems were those of Ibrahim et al. ^[22], Nawaz et al. ^[23], and Al-Odat and Al-Ghamdi ^[24].

The works on convection in couple-stress liquids cited above are, however, silent about the following subjects:

(ⅰ) Takens-Bogdanov bifurcation,

(ⅱ) codimension-two bifurcation,

(ⅲ) inverted bifurcation,

(ⅳ) Eckhaus and zigzag instabilities,

(ⅴ) heat and mass transport.

These aforementioned unconsidered aspects of the two-component thermoconvective instability problem in couple-stress liquids with the Soret effect are studied in the paper using the Lorenz and Newell-Whitehead-Segel equations ^[25].

2 Basic equations

Consider an infinite extent horizontal couple stress liquid layer of thickness d. The upper and lower boundaries are held at a constant temperature T₀ and T₀ + △T (△T > 0), respectively. The solutal concentrations are maintained at S₀ and S₀+△S, respectively. The bounding surfaces of the layer are further assumed to be stress-free, isothermal, and isohaline. The Boussinesq approximation is assumed to be valid. The thermal and solutal gradients lead to simultaneous heat and mass transfer. In mixtures, temperature and concentration gradients induce mass flow, called the Soret effect. With the above effect, heat and mass flows can be expressed as follows:

(1)

(2)

where κ and D_m are the thermal conductivity and the mass diffusivity of species, respectively. The quantity D₁ is the Soret coefficient that arises due to cross diffusion. Here, T^' is the first diffusing component, and S^' is the second diffusing component. Equation (2) signifies the cross-diffusion phenomenon wherein there can be a flux of S^' due to ▽T^'. Similarly, there can be a flux of T^' due to ▽S^', but in the present case, this is assumed to be of negligible magnitude.

Using the above expression, the governing equations ^[2-3] are written as

(3)

(4)

(5)

(6)

(7)

The conduction state is characterized by

(8)

(9)

We use the Cartesian system of coordinates whose dimensional coordinates x^', y^', and z^' are scaled by d. The time t^' is scaled by . The velocity vector V (u^', v^', w^'), the temperature T^', the concentration S^', and the pressure P^' are non-dimensionalized by the scales κ/d, △T, △S, and ρ₀κ²/d². The dimensionless equations for the perturbed quantities of primary and secondary thermosolutal convective instabilities in two-component couple-stress liquids with the Soret effect are

(10)

(11)

(12)

(13)

where

and the Soret coefficient arises due to cross diffusion.

We apply curl two times on the momentum equation (11) and take the z-component of the resulting equation,

(14)

where

From Eqs. (12) - (14), we get

(15)

where

(16)

(17)

In the case of small scale convective motions, we may neglect the inertial term (see Siddheshwar and Pranesh ^[2]), and this is assumed in the paper from here on.

We study the primary instabilities in the next section using the linear stability analysis.

3 Linear stability analysis and study of primary instabilities

We consider stress-free, isothermal, isohaline boundary conditions. In terms of w, the conditions translate into the following ones:

(18)

Substituting the normal mode solution

(19)

into the linearized version of Eq. (15), viz., , we obtain the dispersion relation,

(20)

where

Performing the classical analysis on Eq. (20), we get the expression of the Rayleigh numbers R_Tsc and R_Toc for the stationary and oscillatory modes of convection and the expression for the frequency ω in the form of

(21)

where

(22)

(23)

(24)

and the subscripts sc and oc denote stationary and oscillatory convections, respectively.

The minimum value of R_Tsc is obtained at a minimum (critical) wave number q_scc, which is given below,

(25)

Since the equation governing the critical value of wave number for oscillatory convection does not come in an elegant form as that of q_sc above, the same is not recorded here. In Sections 4 and 5, a local nonlinear stability analysis is presented using the Lorenz equation and the Newell-Whitehead-Segel equation.

Using these equations, the inverted bifurcation and Eckhaus and zigzag instabilities are studied.

4 Derivation of fifth-order Lorenz equation at threshold of stationary convection and study of inverted (subcritical) bifurcation

Eliminating the pressure P from Eq. (11), assuming two-dimensional motions, and using the following minimal representation of Fourier series for u, w, θ, and S:

(26)

(27)

(28)

(29)

in the resulting equations from Eqs. (10), (11), (12), and (13), we get the fifth-order Lorenz model in the forms of

(30)

(31)

(32)

(33)

(34)

We now use the following scaling, keeping in mind the form of the classical Lorenz model which can be obtained as a particular case obtained in the paper,

(35)

Using the scaling (35), the Lorenz system of Eqs. (30) - (34) reduces to

(36)

(37)

(38)

(39)

(40)

where

and the overdot denotes t₁-derivative.

It is a well known fact in the context of the classical Lorenz model that its trajectories remain within the confines of a sphere. The nonlinear terms are responsible for keeping the trajectories confined. Following Siddheshwar and Titus ^[26], the trapping region of the trajectories of the solution of the Lorenz model in Eqs. (36) - (40) can be obtained in the form of

(41)

From the above equation, we note that the trapping region is a four-ellipsoid in a five-dimensional phase-space.

The Lorenz model of Eqs. (36) - (40) is intractable, but in the steady state, it does possess a solution. The steady-state solution of the Lorenz model of Eqs. (36) - (40) is

(42)

(43)

(44)

(45)

(46)

where

(47)

(48)

(49)

Equating the discriminant in Eq. (42) to zero, the finite-amplitude Rayleigh number R_T^f is obtained as a quadratic equation as follows:

(50)

Solving Eq. (50), we get

(51)

where

Equation (51) will be used to explore the possibility of inverted (subcritical) bifurcation, and the same is discussed in the results and discussion section. In the next section, we derive the Newell-Whitehead-Segel equation with the intention of using it to study two secondary instabilities, namely, Eckhaus and zigzag instabilities.

5 Derivation of nonlinear two-dimensional Newell-Whitehead-Segel equation at threshold of stationary convection and study of secondary instabilities

In this section, we follow Newell and Whitehead to obtain the Newell-Whitehead-Segel equation governing the amplitude of the imposed disturbance. We assume the manifestation of cylindrical rolls with the axis parallel to the y-axis so that the y-derivatives need not be considered in Eq. (15). The z-dependence comes through the trigonometric functions which ensures that the stress-free boundary condition is satisfied. We use the expansion parameter ϵ as

(52)

to obtain the system of equations of several orders in powers of ϵ. However, the variables x, y, z, and t are scaled as follows:

(53)

The inherent symmetry breaking of instability is reflected through the difference of scaling in the two horizontal directions. In view of Eq. (52), the differential operators can be expressed as

(54)

Due to the fact that we are making a local nonlinear stability analysis, R_Tsc is assumed to take values close to the threshold value, R_Tscc, i.e., ϵ ≪ 1. The solution to Eqs. (10) - (13) is now assumed as follows:

(55)

(56)

(57)

(58)

(59)

The first approximation is given by the eigenvector of the linearized problem,

(60)

(61)

(62)

(63)

(64)

where c.c. stands for complex conjugate. Using Eqs. (53) and (54), the linear operator in Eq. (16) can be written as

(65)

where

(66)

(67)

(68)

Similarly, the nonlinear term in Eq. (17) is taken to be

(69)

where

(70)

(71)

(72)

Substituting Eqs. (54), (57), (65), and (69) into Eq. (15), we get

(73)

(74)

(75)

Equation (73) gives us the expression for the critical Rayleigh number for the onset of stationary convection, viz., Eq. (21).

Using Eqs. (62), (67), and (70) in Eq. (74) and solving them, we get w₁=0. Using the continuity equation (10), we get u₁=0. The equations for θ₁ and S₁ are given by

(76)

(77)

Substituting Eqs. (54), (57), (58), (59), and (65) into Eqs. (76) and (77), taking u₁=w₁=0, and solving the resulting equations, we get

(78)

(79)

Substituting Eqs. (62), (67), (68), (71), (78), and (79) in Eq. (75) with w₁=0 and using the Fredholm alternative condition on the resulting equation, we get the Newell-Whitehead-Segel equation in the form of

(80)

where

(81)

(82)

(83)

(84)

(85)

Since we consider cylindrical rolls with the axis parallel to the y-axis, we neglect the y-derivatives and the time-derivatives. Equation (80) takes the form of

(86)

Since λ₁>0, the solution to Eq. (80) turns out to be

(87)

(88)

In the succeeding section, Eq. (80) will be used to study the secondary instabilities of Eckhaus and zigzag.

5.1 Secondary instabilities

In terms of fast variables x, y, z, and t, Eq. (80) may be written as

(89)

We now consider the phase-winding solution to Eq. (89) in the form of

(90)

Substituting Eq. (90) into Eq. (89), we get

(91)

The steady-state uniform solution to Eq. (91) is

(92)

To study the secondary instabilities, we impose an infinitesimal perturbation on the uniform steady state solution in the form of

(93)

Substituting Eq. (93) into Eq. (91) and equating real and imaginary parts on either side of the resulting equation, we get

(94)

(95)

We analyze Eqs. (94) and (95) by using the normal modes in the form of

(96)

Substituting Eq. (96) into Eqs. (94) and (95), we get

(97)

(98)

where

For a non-trivial solution to Eqs. (97) and (98), we require

(99)

We are not interested in a negative root of Eq. (99), since the mode would then be stable. The positive root of Eq. (99) is

(100)

For this positive root, the rolls can become unstable.

5.2 Eckhaus instability

Substituting q_y=0 in Eq. (99), we get

(101)

The possibilities regarding the roots of Eq. (101) are as follows:

(ⅰ) both roots are negative, or

(ⅱ) one root is positive, and the other is negative.

When both roots of Eq. (101) are negative (Case (ⅰ) above), the product of roots is positive. Both the roots being negative refers to a stable situation.

In Case (ⅱ), the product of roots is negative, and this refers to an unstable situation due to the positive root. By the Descartes rule of sign, the following condition ensures the existence of a positive root:

(102)

Since q_x²>0, we require

(103)

and this condition defines the domain of Eckhaus instability in the plane.

5.3 Zigzag instability

Let us now take q_x=0 in Eq. (99) to get

(104)

where

Following the procedure of Subsection 5.2, we can obtain the domain of the zigzag instability in the form of

(105)

Having thus far studied the possibility of different types of primary and secondary instabilities, we now quantify the heat and mass transport by considering the stationary instability (primary instability).

6 Heat and mass transports at lower boundary

The heat transport is quantified through the thermal Nusselt number Nu which is defined as

(106)

where

Substituting Eqs. (8), (52), (58), (63), and (78) into Eq. (106) and completing integration, we get

(107)

Similarly, mass transport is quantified through the Sherwood number Sh which is defined as

(108)

Substituting Eqs. (8), (9), (52), (58), (59), (63), (64), (78), and (79) into Eq. (108), we get

(109)

In computing Nu and Sh, we make use of Eq. (87) in Eqs. (107) and (109). The next section documents the results and discussion extracted from the computation.

7 Results and discussion

Primary and secondary instabilities are discussed in the paper concerning double diffusive convection in a couple stress liquid, in the presence of the cross diffusion effect of Soret. We first discuss the variation of the critical Rayleigh number for the stationary convection R_Tscc and the oscillatory convection R_Tocc with the couple stress parameter C_S for different values of the solutal Rayleigh number R_S, the Lewis number Le, the Soret parameter D_S, and the Prandtl number Pr, and those are shown in Tables 1 and 2. In the above parameters, C_S represents the suspended particles, R_S represents the Rayleigh number corresponding to the second diffusing component, and D_S represents the coupled molecular diffusion or cross-diffusion.

From the tables, one may make the inference that in the presence of suspended particles (C_S≠0), the system is stabilized in the cases of both stationary and oscillatory convection.

The effect of R_S on the onset of stationary and oscillatory convection is also shown in Tables 1 and 2. From the tables, we may conclude that when the magnitude of R_S increases, the critical Rayleigh numbers of stationary and oscillatory modes increase, thereby indicating that the solutal concentration (second diffusing component) has a stabilizing effect on the onset of convection.

The tables show the effect of the Lewis number Le on the onset of convection for fixed values of the other parameters. With an increase in Le, one can see that the critical Rayleigh number increases in the cases of both stationary and oscillatory convection.

For positive values of D_S, one can observe from the tables that R_Tscc increases with the increase in D_S, while for negative D_S, R_Tscc decreases with the increase in |D_S|. A similar observation on the effect of D_S on R_Tocc can be made. In the case of D_S taking negative values, oscillatory convection is shown to be possible. In the case of positive values of D_S, however, the PES is valid provided D_S≥0.3. Overstability is seen only for the case when the sign of D_S is such as to make a stabilizing contribution to the density gradient.

The tables illustrate the fact that the Prandtl number has a destabilizing effect on oscillatory convection and has no effect on stationary convection.

We now discuss about the possibility of pitchfork, Hopf, Takens-Bogdanov, and codimension-two bifurcations (primary instabilities) using the results of linear stability analysis. The schematic of pitchfork and Hopf bifurcations are shown in Fig. 1.

Based on the sign of D_S, we may consider the following two cases:

(ⅰ) D_S>0 and

(ⅱ) D_S < 0.

The pitchfork and Hopf bifurcations are characterized by stationary and oscillatory Rayleigh numbers given by Eq. (21). The finding on the effect of D_S on R_Tsc and R_Toc is in tune with the results of Malashetty et al. ^[6] and Hu et al. ^[21].

From Eq. (84), it is clear that the supercritical pitchfork, subcritical and tricritical bifurcations are, respectively, given by bifurcation λ₃>0, λ₃ < 0, and λ₃=0 (see Fig. 2). The Takens-Bogdanov bifurcation point is an equilibrium point (here, it is a wave number), at which the neutral curve of oscillatory convection intersects the neutral curve of stationary convection, i.e.,

(110)

where q_oc=q_sc. If q_oc ≠ q_sc, it is called the codimension-two bifurcation point.

Analysing the expression of ω² in Eq. (21), we note that when I₃=0, the frequency on the neutral curve of oscillatory convection approaches zero, i.e.,

(111)

Equation (111) is a double zero corresponding to the Takens-Bogdanov bifurcation point at R_S(q_oc), where

(112)

The schematic of the Rayleigh number versus the wave number plot corresponding to Takens-Bogdanov and codimension-two bifurcations is shown in Fig. 2.

We now proceed to discuss the results of a nonlinear stability analysis using the Lorenz and Ginzburg-Landau models. The Lorenz model is used to discuss about the possibility of inverted (subcritical) bifurcation. Our findings from the weakly nonlinear stability analysis concur with those of Platten and Chavepeyer ^[28] and Malashetty et al. ^[7] that subcritical instability is possible for D_S < 0.

Subsection 4.1 discusses the condition under which the secondary instabilities exist. The Ginzburg-Landau equation is used in studying secondary instabilities such as Eckhaus and zigzag.

It is found from our computations that the regions of Eckhaus and zigzag instabilities increase by increasing |D_S| for fixed values of the other parameters (see Fig. 3). However, D_S does not significantly alter the region of both Eckhaus and zigzag instabilities.

We evaluate the Nusselt number Nu and the Sherwood number Sh to study heat and mass transport. Figure 4 shows the plots of Nu versus the scaled thermal Rayleigh number r₁ for different values of the Soret parameter D_S and keeping the other parameters fixed. From the figure, it is clear that for D_S>0, heat transport is diminished, while for D_S < 0, heat transport is enhanced compared with the value of Nu of D_S=0.

Figure 5 shows that as we increase Le, Nu increases for D_S < 0, while the opposite result is true for D_S>0. Thus, Le enhances heat transport in the case of D_S < 0 and diminishes it in the case of D_S>0. A similar effect can be observed in Fig. 6 which is the plot of Nu versus r₁ for different values of C_S.

Figure 7 throws up an interesting result in the case of D_S=0.3, and the similar result applies for other positive values of D_S. In this case, the second diffusing component fails to alter heat transport of a single component system, and this is essentially due to the cross-diffusion effect signified by the Soret coefficient nullifying the solutal concentration effect on heat transfer. However, in the case of D_S < 0, heat transport is significantly altered by the second diffusing component. We also find from the figure that the Nusselt number increases with the increase in R_S.

We next discuss the result on the Sherwood number. From Fig. 8, it is obvious that the solutal transport decreases with an increase in the value of D_S for D_S>0, a result similar to the result on Nu in Fig. 4. Likewise, the results of Figs. 4 and 9 are similar for the case of D_S < 0. The effect of Le on Sh for the cases D_S>0 and D_S < 0 are qualitatively similar to those of the Nusselt number depicted in Fig. 5 but the only difference is that when Le significantly influences Nu, there is comparatively insignificant influence of Le on Sh. From Eqs. (107) and (109), we may write

(113)

This aspect can be explicitly seen by comparing the corresponding plots of Nu and Sh in Figs. 4-11.

An important related paper to the current one is the excellent paper by Hu et al. ^[21] who not only iterated that Poiseuille-Bénard flows are subcritically unstable for negative separation factors but most importantly established that, for very small negative or large positive separation factors and large Rayleigh numbers, the maximum transient growth rate gets enlarged. In the absence of a Poiseuille flow in their Bénard configuration, the results qualitatively match with that of Hu et al. ^{[21, 27]}. To compare the current paper with the paper of Hu et al. ^[21], we note the following one-to-one correspondence between the papers:

(ⅰ) The separation ratio ψ of Hu et al. ^[21] is of our paper and their Lewis number Le is reciprocal of our Le.

8 Conclusions

(ⅰ) The second diffusing component stabilizes the system due to the fact that the solutal concentration increases the density. The Soret effect characterized by D_S makes a stabilizing or destabilizing contribution to the density gradient depending on the sign of D_S. Overstability is seen only for the case when the sign of D_S is such as to make a stabilizing contribution to the density gradient.

(ⅱ) The suspended particles increase the viscosity of the fluid and thereby stabilize the system. The presence of suspended particles is characterized by C_S.

(ⅲ) In addition to the pitchfork and Hopf bifurcations, Takens-Bogdanov and codimension-two bifurcations are also possible in the system.

(ⅳ) The inverted (subcritical) bifurcation is possible only for D_S < 0.

(ⅴ) D_S does not alter the regions of the Eckhaus and zigzag instabilities.

(ⅵ) for D_S>0, and for D_S < 0.

(ⅶ) Variations of R_Tsc, R_Toc, Nu, and Sh with Le, C_S, and R_S for positive and negative values of D_S are tabulated in Table 3.

Acknowledgements The second author (C. KANCHANA) is grateful to the University Grants Commission (UGC), New Delhi, India for supporting her research work with a Rajiv Gandhi National Fellowship. The authors are grateful to the two reviewers for their educative comments that helped us refine the paper to the present form.