#### Abstract

A chaotic system arising from double-diffusive convection in a fluid layer is investigated in this paper based on the theory of dynamical systems. A five-dimensional model of chaotic system is obtained using the Galerkin truncated approximation. The results showed that the transition from steady convection to chaos via a Hopf bifurcation produced a limit cycle which may be associated with a homoclinic explosion at a slightly subcritical value of the Rayleigh number.

#### 1. Introduction

The concept of sensitivity on initial conditions where a small difference on initial conditions may produce large variations in the long-term behaviour of the system is pivotal in chaos theory. This behaviour is also known as the “*butterfly effect*” related to work done by Lorenz [1] where it is already described by Henri Poincare in 1890 in the literature in a particular case of the three-body problem. Chaotic behaviour has been studied intensively in various dynamical systems; see, for example, [2–14].

The investigation of free convection in the Rayleigh-Bénard problem is receiving much attention due to its wide application in different disciplines such as biotechnology for the description of the convection with the microorganisms diffusion, in astrophysics for simulation of the influence of the helium diffusion on convective motions in the stars, in oceanography for the investigation of the salinity influence on the convective motions in the seas, and in engineering and geology. Research in double-diffusive convection begins after the work done by sea-going oceanographers in order to measure the fluctuation of temperature and salinity as a function of depth as stated in the paper of Huppert and Turner [15]. Then, Knobloch et al. [16] and Bhattacharjee [17] studied the transition to chaos in double-diffusive convection with stress-free boundary conditions where oscillatory solution exists. They showed that the instability of fluid becomes oscillatory when thermal Rayleigh number is raised and the truncated model suggests that the appearance of chaos is associated with heteroclinic bifurcations.

Two-dimensional thermosolutal convection between free boundaries was studied numerically by Veronis [18]. From their observation, they found that when the solutal Rayleigh number is large enough, the oscillations underwent a bifurcation to asymmetry as thermal Rayleigh number increased and, for the larger values of solutal Rayleigh number, the transition from chaos to steady motion occurs.

Sibgatullin et al. [19] studied some properties of two-dimensional stochastic regimes of double-diffusive convection in a plane layer. Using the Bubnov-Galerkin method, they obtained that, with the growth of Rayleigh numbers of heat and salinity, the structure of one-dimensional curve becomes more irregular and sophisticated. Transition to chaos in double-diffusive Marangoni convection was studied by Li et al. [20]. It was found that the supercritical solution branch takes a quasiperiodicity and phase locking route to chaos while the subcritical branch follows the Ruelle-Takens-Newhouse scenario. The transitions from regular to chaotic dynamics and analysis of the hyper, hyper-hyper, and spatial-temporal chaos using the Lyapunov exponents of continuous mechanical systems have been studied in [21–24]; they found the Sharkovskii windows of periodicity in the systems investigated.

The objective of the present paper is to study the weak turbulence and chaos in double-diffusive convection involving temperature and concentration as the thermal Rayleigh number increases with rigid, no-slip horizontal boundary condition. Applying the truncated Galerkin approximation to the governing equations yields an autonomous system with five ordinary differential equations which can be used to understand low-dimensional dynamics before moving to studying more complex systems.

#### 2. Mathematical Formulation

Consider a two-dimensional layer of fluid of depth subject to gravity and heated from below as shown in Figure 1. A Cartesian coordinate system is used such that the vertical axis is collinear with gravity, that is, . The two long walls are maintained at temperatures and and solute concentrations and , respectively. A relationship between density, temperature, and solute concentration is assumed linear and can be presented by the following form , where and are volume expansion coefficients due to variations of thermal and solute concentrations. The Boussinesq approximation is applied for the effects of density variations for the gravity term in momentum equation. Therefore, the set of equations governing the conservation of mass, momentum, energy, and concentration for fluid flow is given by the following: We nondimensionalize (1)–(4) using the following transformations: where is the velocity component, is the pressure, and are the temperature and solute concentration variations, is the effective thermal diffusivity, and is fluid’s viscosity.

In this model, all the boundaries are rigid and the solution must follow the impermeability conditions there, that is, on the boundaries, where is a unit vector normal to the boundary. The temperature and solute concentration boundary conditions are at and at .

For convective rolls having axes parallel to the shorter dimension (i.e., ) , by applying the curl operator on (2) to eliminate the pressure and introducing the stream function defined by and , we get where which are, respectively, the Prandtl number, the Rayleigh number, the solutal Rayleigh number, and the Lewis number. The boundary conditions for the stream function are on the horizontal boundaries. Equation (6) forms a nonlinear coupled system which together with the corresponding boundary conditions allows for a basic motionless conduction solution.

#### 3. Diminished Set of Equation

In order to obtain the solution to (6), we represent the stream function, temperature, and solutal distributions in the following form: Substituting (8) into (6), multiplying the equations by the orthogonal eigenfunctions corresponding to (8), and then integrating them over the spatial domain yield a set of five ordinary differential equations for the time evolution of the following amplitudes: In (9), the time, the amplitudes, the Rayleigh number, and the solutal Rayleigh number were rescaled, and the following notations are introduced as follows: Rescaling the equation again in the forms gives the following set of scaled equations which are equivalent to (9): where the dots denote time derivatives .

#### 4. Linear Stability Analysis

In this paper, we investigate the chaotic behaviour with low Prandtl number in double-diffusive convection. We obtained system (12) that provides a set of nonlinear equations with five parameters. The value of has to be consistent with the wave number at the convection threshold, a requirement for the convection cells to fit into the domain and fulfill the boundary conditions. However, the Lorenz equations have been extensively analyzed and solved for parameter values corresponding to convection in pure fluids and, even there, the parameter values most regularly used correspond to and . Therefore, it is of interest to analyze and solve the corresponding equations for parameter values corresponding to the problem under investigation. We employ the MATLAB ODE45 routine for obtaining the numerical solutions.

Before attempting the numerical solution of system (12), it is useful to examine the local stability of equilibrium points. System (12) has the three basic properties which we will discuss in the following: dissipation, fixed points, and stability of fixed points.

##### 4.1. Dissipation

System (12) is dissipative since Therefore, if the set of initial points in the phase space occupies region at , then, after some time, , the endpoints of the corresponding trajectories will fill a volume The expression indicates that the volume decreases monotonically with time.

##### 4.2. Fixed Points

The fixed points for velocity, temperature, and solute concentration can be obtained by setting the derivatives of system (12) to zero: There is one trivial solution, that is, the origin in the phase space which corresponds to the motionless solution. The other nonzero fixed points are given by the following: where The system has five fixed points. When , the five fixed points are all real. Thus, when , is always real and and are always complex; therefore, the three fixed points are all complex and the other two fixed points, , are all real for . The fixed point corresponds to motionless solution and , where corresponds to the convective solution.

##### 4.3. Stability of the Fixed Points

The Jacobian matrix of (12) can be written as follows: Since the matrix is , it is hard to obtain the eigenvalues in a closed form. Hence, the numerical calculation can be performed to discuss the stability at the fixed point. The motionless solution loses stability and the convection solution takes over at the fixed point with the critical value, . Numerical results for the value of , which corresponds to the onset of convection, is obtained for various values of and with the value of parameters and as shown in Figure 2. Increasing the values of and increases the value of .

**(a)**

**(b)**

The stability of the fixed points () is associated with the convective solution. The evolution of the complex eigenvalues of in the case , , and is plotted as shown in Figure 3(a) for the fixed point . These two roots become a complex conjugate at for the case , respectively. At these points exactly, they still have negative real parts; therefore, the convection fixed points are stable, that is, spiral nodes. Of all the cases, both the imaginary and real parts of these two complex conjugate eigenvalues increase as increases and they cross the imaginary axis on the complex plane, so as a result their real part becomes nonnegative at a value of . At these points, the convective fixed points lose their stability and another periodic, quasiperiodic or chaotic solution takes over. Figure 3(b) shows the evolution of the complex eigenvalues of for the case , and for the fixed point . These two roots become a complex conjugate at for the case of , respectively. As mentioned before, the convection fixed points are stable and become spiral nodes. Their real part becomes nonnegative at a value of . Therefore, at these points, the convective fixed points lose their stability and another periodic, quasiperiodic, or chaotic solution takes over.

**(a)**

**(b)**

While for the fixed point , the evolution of the complex eigenvalues of with the same parameter values is always on the positive side of , does not cross the zero axis for , and is not of interest in this study.

#### 5. Results and Discussion

##### 5.1. Bifurcations and Transition to Chaos

In this study, we focused on the dynamic behaviour of thermal convection in double-diffusive fluid layer. The values of and used in all computations are and , respectively, which are consistent with the critical Rayleigh number () and the critical wave number at marginal stability in fluid layer convection. All solutions were obtained using the same initial conditions, which were selected to be in the neighborhood of the positive convection fixed point. The initial conditions are at . All computations were carried out with the value of maximum time, , and a step size using the built-in ODE45 method in MATLAB R2010a.

The bifurcation diagrams illustrated in Figure 4 show the peaks and valleys in the posttransient values of versus . In Figure 4(a), for we have one-point attractors, but the “attracted” value of increases as increases, at least to . Bifurcation occurs at until just beyond , where the system is chaotic. However, the system is not chaotic for all values of , and we will discuss it using phase-portrait diagram. When we fix and increase from to , the range of one-point attractor changes to ; this is shown in Figure 4(b), while in Figure 4(c), one-point attractor dropped to in the case of and . Here we can conclude that increasing the value of the solutal Rayleigh number (with fixed value of Lewis number) will delay the convection process. But increasing the value of Lewis number (with fixed value of solutal Rayleigh number) will enhance the onset of chaos.

**(a)**

**(b)**

**(c)**

Figure 5 shows the projections of the trajectory’s data points on the plane for and . From Figure 5(a), we obtain a solitary limit cycle signifying the loss of stability of the steady convective fixed points. The subcritical value for this transition is . Figure 5(b) shows the projections of the trajectory’s data points for . At this point, the homoclinic explosion occurs and the chaotic regime with the strange attractor takes over. The homoclinic explosion behaviour giving birth to a stable periodic orbit with period-8 at is presented in Figure 5(c). This corresponds to the first wide periodic window within the chaotic regime. In Figure 5(d) we can observe that the data points do align in such a way as to produce an almost clear projection of unconnected points on the projected plane. Increasing the value of further shows the dynamical behaviour’s return to being chaotic again at as shown in Figure 5(e). At and , we have a period-8 and period-4 periodic solutions as shown in Figures 5(f) and 5(g). Figure 5(h) shows that a period-2 periodic solution takes over at and a period-2 periodic type remains when the solutions at higher values of are obtained. We conclude the observation around these regimes of periodic windows within the broadband of chaotic solutions by pointing out a sequence of period-halving as one increases the Rayleigh number.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

Figure 6 shows the projections of the trajectory’s data points on the plane for and . The subcritical values for limit cycle and homoclinic explosion occur at and as shown in Figures 6(a) and 6(b), respectively, while period-8 is observed at and as presented in Figures 6(c) and 6(d). The dynamical behaviour returns to being chaotic again as increases; this happens at as shown in Figure 6(e). Figures 5(f) and 5(g) show a period-8 and period-4 periodic solutions at and . Figure 6(h) shows that a period-2 periodic solution takes over at and remains the same behaviour when the solutions at higher values of are continued.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

##### 5.2. Lyapunov Exponents

The convergence plot of the Lyapunov spectrum for system (12) is shown in Figure 7. The alogrithm as proposed by [24] was employed for this purpose. The values of the Lyapunov exponents for system (12) are tabulated in Table 1. From these results, we can conclude that for eigenvalue the system is always unstable and chaotic with the increasing . For , the system alternates between being stable and dissipative to unstable and chaotic when the value of is increased. For eigenvalues , and , the system is always stable and dissipative with the increasing .

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

#### 6. Conclusion

In this work, chaotic behaviour in double-diffusive convection in a fluid layer has been investigated. A five-dimensional model of chaotic system was obtained using the Galerkin truncated approximation. The transition from steady convection to chaos via a Hopf bifurcation produced a limit cycle which may be associated with a homoclinic explosion at a slightly subcritical value of the thermal Rayleigh number. Both the solutal Rayleigh number and Lewis number affect the stability of the system. Increasing the Rayleigh number shows that the trajectory of the data points pointing out a sequence of period-halving and the behaviour remains the same at a higher Rayleigh number. The different transitions of the system, can be implied by the different values of the Lyapunov exponents. Negative eigenvalues lead to a stable and dissipative system and positive eigenvalues show that the system is always unstable and chaotic, while alternate eigenvalues suggest different transitions of the system (i.e., stable and dissipative to unstable and chaotic) as the value of the Rayleigh number is increased.

#### Acknowledgments

The authors acknowledge the financial support from Universiti Malaysia Terengganu (Grant no. TPM-68006/2012/49), Universiti Kebangsaan Malaysia (Grant no. DIP-2012-12), and University of Malaya (Grant no. RG172-11AFR).