Chaotic Convection in a Viscoelastic Fluid Saturated Porous Medium with a Heat Source

Chaotic convection in a viscoelastic fluid saturated porous layer, heated from below, is studied by using Oldroyd’s type constituting relation and in the presence of an internal heat source. A modified Darcy law is used in the momentum equation, and a heat source term has been considered in energy equation. An autonomous system of fourth-order differential equations has been deduced by using a truncated Fourier series. Effect of internal heat generation on chaotic convection has been investigated. The asymptotic behavior can be stationary, periodic, or chaotic, depending upon the flow parameters. Construction of four-scroll, or “two-butterfly,” and chaotic attractor has been examined.


Introduction
Transport phenomenon in a fluid saturated porous medium is of great practical importance in many areas such as geothermal energy utilization, oil reservoirs, solar energy storage systems, passive cooling of nuclear reactors, pollutant transport in ground water, and storage of chemical and agricultural products, to mention a few.The problem of convection in a fluid saturated porous medium has been studied during the past few decades due to its applications in thermal and engineering sciences.An interesting problem was studied by Horton and Rogers [1] and independently by Lapwood [2], who addressed the Rayleigh-Bénard convection in porous media.Katto and Masuoka [3] employed Darcy's law to express the fluid characteristics in porous layer and experimentally showed the effect of Darcy's number on the onset conditions of buoyancy-driven convection.Important reviews of most of the findings on convection in porous medium are given by Ingham and Pop [4], Vafai [5], and Nield and Bejan [6].
The concept of chaos was first introduced by Poincaré [7,8], who investigated orbits in celestial mechanics and realized that the dynamical system generated by the three-body problem is quite sensitive to the initial conditions exhibiting chaotic behavior.Since the introduction of the chaotic attractors by Lorenz [9] to study atmospheric convection, many chaotic systems have been introduced, such as Rössler [10], Chen and Ueta [11] systems.Because of their potential applications in engineering, the study of chaotic systems has attracted the interest of many researchers.Recently Vadász et al. [12] have investigated the effect of vertical vibrations on chaotic convection in porous medium employing Darcy model.Their results show that periodic solutions and chaotic solutions alternate as the value of the scaled Rayleigh number varies, when forced vibrations are present.Very recently Kiran and Bhadauria [13] have studied chaotic convection in a Newtonian fluid saturated porous medium under temperature modulation at the boundaries.They found that the effect of temperature modulation is to enhance the behavior of chaotic motion.Very recently Bhadauria et al. [14] and Bhadauria and Kiran [15] have studied the chaotic convection using different models.
Although viscoelastic fluid in porous media has been considered many years before by Marshall and Metzner [16] and James and McLaren [17], it is only recently that attention has been given to convection in viscoelastic fluid saturated porous media.Kim et al. [18] studied thermal instability of viscoelastic fluid in porous media by making linear and nonlinear stability analyses and obtained the stability criteria The motive for the present work is to study the effect of internal heat source on dynamics of convection in a viscoelastic fluid saturated porous medium.In the momentum equation, the viscoelastic model of the Oldroyd type is considered by a modified Darcy law.To deduce fourdimensional system, the Fourier series expansion has been applied to the governing equations of the thermal convection in a viscoelastic fluid saturated porous medium.The effects of internal heat source, relaxation, and retardation parameters on the dynamics of the system have been examined in detail.Further, the present system has been reduced to some famous systems provided in the literature.

Mathematical Formulation
An infinitely extended horizontal viscoelastic fluid saturated porous layer, confined between two impermeable boundaries at  = 0 and  = , heated from below and cooled from above, has been considered.A Cartesian frame of reference is chosen in such a way that the origin lies on the lower plane and the -axis as vertical upward.Adverse temperature gradient is applied across the porous layer and the lower and upper planes are kept at temperatures  0 + Δ and  0 , respectively.Oberbeck-Boussinesq approximation is applied to account the effect of density variations.Under these postulates the governing equations for thermal convection in a viscoelastic fluid saturated porous medium are given by where q is velocity (, V, ),  is the dynamic viscosity,  is permeability,   is the thermal diffusivity,  is temperature,   is thermal expansion coefficient, and  is the density, while  0 and  0 are the reference density and temperature, respectively.The externally imposed thermal boundary conditions are given by  =  0 + Δ, at  = 0,  =  0 , at  = , (5) where Δ is the temperature difference across the porous medium.

Basic State
The basic state is assumed to be quiescent, and the quantities in this state are given by Substituting ( 6) in ( 1)-( 4), we get the following relations, which helps us to define basic state pressure and temperature: The solution of ( 8), subject to the boundary conditions (5), is given by The finite amplitude perturbations on the basic state are superposed in the following form: We introduce (11) and the basic state temperature field given by ( 10) in ( 1)-( 4).The resulting equations are then nondimensionalized using the following transformations: Introducing the stream function  such that  * = −/ * and  * = / * in the resulting momentum equation (1) and then eliminating the pressure  by operating curl on it, we get the following momentum and energy equations in nondimensionalized form (dropping the asterisks) as where Va = Pr/Da is the Vadasz number, Ra =   (Δ)/]  is the Rayleigh number,   =  2 /  is the internal Rayleigh number, Da = / 2 is the Darcy number, and Pr = ]/  is the Prandtl number.
The nondimensional basic temperature field   , which appears in (14), can be obtained from the expression (10) as

Mathematical Solution
To obtain the solution of the nonlinear coupled system of partial differential equations ( 13)-( 14), we represent the stream function and temperature in the following Fourier series expressions [12,13,15]: Substituting expressions ( 16) in ( 13)-( 14), using the orthogonality condition with the eigenfunctions associated with expressions (16), and integrating over the domain, we get a set of ordinary differential equations for the time evolution of the amplitudes, in the form where Γ is nondimensional relaxation time, Λ is ratio of retardation time to relaxation time, and  2 = ( 2 +  2 ).Also In ( 17) above, time  has been rescaled as  = ( 2 +  2 ).The above low-order spectral model may qualitatively reproduce convective phenomenon observed in the full system.Further, the solution of the system can be used as initial values in studying the fully nonlinear convection problem.Now, for our convenience, we use the following notations: Then after rescaling the amplitudes in the form we get the following set of equations: where "⋅" denotes the derivative with respect to the scaled time .

Equilibrium Points.
Setting the time derivatives of system (21) to vanish, we obtain the equilibrium points for velocity and temperature fields as Solving the above algebraic equations, we get the trivial solution which corresponds to pure heat conduction solution.This is known to be a possible solution though it is unstable when (Ra) is sufficiently large.The other two equilibrium points are the solutions ( 2,3 ,  2,3 ,  2,3 ,  2,3 ) characterize the onset of finite amplitude steady motions, where  > /.

Stability of the Equilibrium Points.
The Jacobian matrix of system (21) may be written as The stability of the fixed point corresponding to pure conduction solution ( 1 =  1 =  1 =  1 = 0) is governed by the roots of the following characteristic polynomial equation for the eigenvalues,   ( = 1-4): Stability depends upon the value of Γ; for Γ < ( − )/[(Λ − 1)] there is an exchange of stability, and for other two steady state solutions origin loses its stability.When Γ > ( − )/[(Λ − 1)], there is a pair of pure imaginary roots of (32).The oscillatory or overstable solutions arise at a critical value of Rayleigh number given by The stability of the fixed point corresponding to convection solution ( =  = ±√( − ),  = ( − ),  = 0) is governed by the roots of the following characteristic polynomial equation for the eigenvalues,   ( = 1-4): The steady state solutions are useful because they predict that a finite amplitude solution to the system is possible for subcritical values of the Rayleigh number and that the minimum values of Ra for which a steady solution is possible lie below the critical values for instability to either a marginal state or an overstable infinitesimal perturbation.

Results and Discussions
In the previous section, the set of results were obtained for supercritical values of  with the effect of internal heat generation.All the calculations were taken over using Mathematica's inbuilt forth-order Runge-Kutta method.Solutions were obtained using the same initial conditions which were selected to be in the neighbourhood of the positive convective fixed point.The common initial conditions, (0) = (0) = (0) = 0.9 and (0) = 0.1, have been chosen.The time domain is taken as 0 ≤  ≤ 100 with a constant time step, Δ = 0.001.In the paper, we demonstrate the effects of internal heat source and viscoelastic parameters on the system, in the form of space projections of trajectories onto the -, -, and - planes, as the value of  increases.
In Figure 1, the initial supercritical solutions are presented at fixed value of   = 0.1[1( 1 - 4 )],   = 2[1( 5 - 8 )], and   = 5[1( 9 - 12 )] with Λ = 0.7 and have been projected onto - planes.From Figure 1, we observed that for weak internal heating, that is,   = 0.1 at   = 0.987361, the motionless solution loses its stability and the convection solution takes over.For  = 2, trajectories move towards steady convection stability point on a straight line for a Rayleigh number slightly above the loss of stability of the motionless solution.From Figure 1, it is evident that there is phase (spiral) trajectory for  = 8, as the flow exhibits an oscillatory decay.In Figure 1, spiralling approach of the trajectories towards the steady state fixed point is quite pronounced.This behavior is also inferred from the roots of (32).The flow undergoes a homoclinic bifurcation, as shown in Figure 1 for  = 9.294, similar to that predicted by the Lorenz [9] equations.This is known as a global bifurcation which cannot be detected through local stability analysis around the fixed point.On further increasing the value of , the flow becomes completely chaotic as can be depicted from the spatiotemporal structure of flow.The transition to chaos in this case is similar to that leading to the Lorenz attractor.It does not follow any of (a 1 ) the well known routes to chaos, such as through period doubling, quasiperiodicity, or intermittency [43].The behavior of the system at  = 28 is chaotic which is confirmed by the broadening of the base in the power spectrum.For the moderate heating, that is,   = 2 at   = 0.756963, the motionless solution loses its stability and the convection solution takes over.For  = 2, trajectories move towards steady convection stability point on a straight line for a Rayleigh number slightly above the loss of stability of the motionless solution.At  = 8, it shows the spiral approach before it is attracted towards steady state fixed point.As the value of  increased to  = 9.294, the flow becomes homoclinic in nature, and, on further increasing to  = 28, a transition to chaos occurs.In the case of strong heating, for   = 5, at   = 0.430905, the motionless solution loses its stability and the convection solution takes over.For  = 2, trajectories move towards steady convection stability point on a straight line for a Rayleigh number slightly above the loss of stability of the motionless solution.In this case for  = 8 spiralling approach is even more pronounced, a fact which is consistent with the eigenvalues.From Figure 1 ( 11 ), for  = 9.294, a homoclinic bifurcation occurs which is more pronounced in this case.Finally at  = 28, transition to chaos occurs.Figure 2 ( 1 - 12 ) is the plots of - plane at   = 0.1[2( 1 - 4 )],   = 2[2( 5 - 8 )], and   = 5[2( 9 - 12 )] with Λ = 0.9.In Figure 2 ( 1 - 12 ), we find qualitatively similar results as shown in Figure 1 ( 1 - 12 ), except that amplitudes take a little bit larger value than for Λ = 0.7.
Figure 3 shows the initial supercritical solutions at fixed value of   = 0.1[1( 1 - 4 )],   = 2[1( 5 - 8 )], and   = 5[1( 9 - 12 )] with Λ = 0.7, projected onto - planes.From Figure 3 (c 1 -c 12 ) we observed that for weak heating, that is,   = 0.1 at   = 0.987361, the motionless solution loses stability and the convection solution takes over.For  = 2 trajectories move towards steady convection stability point through a spiral for a Rayleigh number slightly above the loss of stability of the motionless solution.From Figure 3 ( 2 ,  6 ,  10 ) it is clear that there is spiral trajectory for  = 8 as the flow exhibits an oscillatory decay.In Figure 3 (c 2 , c 6 , c 10 ), spiralling approach of the trajectories towards the steady state fixed point is more pronounced.For  = 9.19, there is a homoclinic pattern of flow, as shown in Figure 3 (c 3 , c 7 , c 11 ).Further increasing the value of , there occurs a transition to chaotic convection as can be seen from Figure 3 (c 4 , c 8 , c 12 ).The transition to chaos in this case is similar to that leading to the Lorenz attractor; one can easily see "two-butterfly" Figure 3 ( 4 ) for  = 28.For the case of moderate heating, that is,   = 2 at   = 0.756963, the motionless solution loses stability and the convection solution takes over.Trajectories move towards steady convection stability point through a spiral for a Rayleigh number slightly above the loss of stability of the motionless solution for  = 2.At  = 8 trajectories move via a spiral, before they are attracted to their steady state fixed points.At  = 9.19, the flow becomes homoclinic in nature and attracted towards a steady state fixed point.On increasing the value of  at  = 28, a transition to chaos occurs and "two-butterfly" Figure 3 ( 4 ) are obtained.For the case of strong heating, namely,   = 5, the motionless solution loses its stability and the convection solution takes over at   = 0.430905.For  = 2 trajectories move towards steady convection stability point through a spiral for a Rayleigh number slightly above the loss of stability of the motionless solution.In this case for  = 8, spiralling approach is even more pronounced, a fact which is consistent with the eigenvalues.From Figure 3 ( 11 ) for  = 9.19, a homoclinic bifurcation occurs which is more pronounced in this case.Finally at  = 28 transition to chaos occurs and "twobutterfly" nature occurs.In Figure 4 ( 1 - 12 ), which are for Λ = 0.9, we found qualitatively similar results to those of Figure 3 ( 1 - 12 ) with slightly large amplitudes.
The initial supercritical solutions, calculated at   = 0.1, 2, and 5 with Λ = 0.7, projected onto - planes have been presented in Figure 5. From Figure 5 ( 1 - 4 ), we observe that for the case of weak heating, that is, at   = 0.1, the motionless solution loses its stability at   = 0.987361, and the convection solution takes over.For  = 2 trajectories move towards steady convection stability point through a spiral for a Rayleigh number slightly above the loss of stability of the motionless solution.From Figure 5 ( 6 ) it is clear that there are spiral trajectory for  = 8 as the flow exhibits an oscillatory decay.In Figure 5 ( 5 ), spiralical approach of the trajectories towards the steady state fixed point is more pronounced.For  = 9.45, there is a homoclinic bifurcation, as shown in Figure 5 ( 7 ).Further increasing the value of , there occurs a transition to chaotic convection as can be seen from Figure 5 ( 8 ).The transition to chaos in this case is similar to that leading to the Lorenz attractor, and one can easily see the transition to chaos for  = 28.For the case of moderate heating, that is,   = 2 at   = 0.756963, the motionless solution loses its stability and the convection solution takes over; the trajectories move towards steady convection stability point through a spiral for a Rayleigh number slightly above the loss of stability of the motionless solution at  = 2.At  = 8 trajectories move via a spiral, before it is attracted to its steady state fixed point.At  = 9.45, the flow becomes homoclinic in nature and attracted towards a steady state fixed point.On increasing the value of  at  = 28, a transition to chaos occurs.For the case of strong heating, namely,   = 5 at   = 0.430905, the motionless solution loses its stability and the convection solution takes over; for  = 2 trajectories move towards steady convection stability point through a spiral for a Rayleigh number slightly above the loss of stability of the motionless solution.In this case for  = 8, spiralling approach is even more pronounced, a fact which is consistent with the eigenvalues.From Figure 5 ( 11 ) for  = 9.45, a homoclinic bifurcation occurs which is more pronounced in this case.Finally at  = 28 there occurs transition to chaos.Figure 6 ( 1 - 12 ) is the plots of initial supercritical solutions, and it is found that the results are qualitatively similar to those obtained in Figure 5 ( 1 - 12 ).

Figure 1 :
Figure 1: Projections and evolution of trajectories over the planes -.

Figure 7 :
Figure 7: Projections and evolution of trajectories.