Abstract
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.
1. 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 for convective flow. Later on Yoon et al. [19] studied the onset of oscillatory convection in a horizontal porous layer saturated with viscoelastic fluid by using linear theory. The linear stability of a viscoelastic fluid saturated densely packed horizontal porous layer heated from below and cooled from above was investigated by Malashetty et al. [20], using the Oldroyd-B type fluid. Since linear stability analysis cannot provide the information about the values of convection amplitudes nor about the aperiodic or chaotic motions, nonlinear analysis has to be performed. Kumar and Bhadauria [21–23] studied nonlinear double diffusive convection in a rotating porous layer saturated with viscoelastic fluid. Recently Bhadauria and Kiran [24, 25] studied oscillatory nonlinear thermal instability in viscoelastic fluids under the effect of temperature and gravity modulations.
Lorenz [9] studied nonlinear analysis of Rayleigh-Bénard convection in a fluid layer considering a truncated representation of Fourier series involving only two terms. Later on Khayat [26–28] and Abu-Ramadan et al. [29] followed Lorenz’s work to derive the fourth-order system in Maxwell fluid and Oldroyd’s fluid layers. Akhatov and Chemberisova [30] studied Rayleigh-Bénard convection in Oldroyd’s fluid saturated porous layer with negligible Darcy effect and deduced a three-order autonomous chaotic system. Chaotic Rayleigh-Bénard convection in Oldroyd’s fluid saturated porous medium using a thermal equilibrium model was investigated by Sheu et al. [31]. Recently Bhadauria and Kiran [15] studied double diffusive chaotic and oscillatory magnetoconvection in a viscoelastic fluid under G-jitter and found that modulated gravity field can be used either to delay or enhance the heat and mass transfer in the system.
Further, there can be a large number of practically important situations, where the porous material offers its own source of heat, thus setting up the convective flow in a different way, through the local heat generation within the porous media. This situation may occur through radioactive decay or through, in the present perspective, a relatively weak exothermic reaction which can take place within the porous material. It is well known that internal heating is the main source of energy for celestial bodies, which keeps the celestial objects warm and active. It is due to the internal heating of the earth that there exists a thermal gradient between the interior and the exterior of the earth’s crust, saturated by multicomponents fluids, which helps convective flow, thereby transferring the thermal energy towards the surface of the earth. Therefore, the role of internal heat generation becomes very important in several applications that include geophysics, reactor safety analyses, metal waste form development for spent nuclear fuel, fire and combustion studies, and storage of radioactive materials. However, there are relatively very few studies avaliable in which the effect of internal heating on convective flow has been investigated. Some of these studies are Bhattacharyya and Jena [32], Haajizadeh et al. [33], Rionero and Straughan [34], Rao and Wang [35], Parthiban and Patil [36], Herron [37], Khalili and Huettel [38], Joshi et al. [39], Bhadauria et al. [40], Bhadauria [41], and Bhadauria et al. [42].
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 four-dimensional 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.
2. Mathematical Formulation
An infinitely extended horizontal viscoelastic fluid saturated porous layer, confined between two impermeable boundaries at 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 and , 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 is velocity (), is the dynamic viscosity, is permeability, is the thermal diffusivity, is temperature, is thermal expansion coefficient, and is the density, while and are the reference density and temperature, respectively.
The externally imposed thermal boundary conditions are given by where is the temperature difference across the porous medium.
3. 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 is the Vadasz number, is the Rayleigh number, is the internal Rayleigh number, is the Darcy number, and is the Prandtl number.
The nondimensional basic temperature field , which appears in (14), can be obtained from the expression (10) as
4. 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 . Also In (17) above, time has been rescaled as . 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 .
Some Special Cases(1)For noninternal heat generation case, we have ; then or , , and , where , and we get which is similar to Sheu et al. [31] though with different coefficients.(2)If we let then the above system takes the form which is identical to Khayat’s [26–28] system with the assumption of the parameters: , , , , and .(3)If, in system (21), we consider the case of Newtonian fluid by taking , , we obtain (4)If, in the above system (25), we consider , then we get or , , and , where , which are famous Lorenz [9] equations.(5)In the limit case , that is, Darcy number , the above system (26) reduces to which is identical to the system given by Akhatov and Chembarisova [30].
5. Stability Analysis
5.1. Equilibrium Points
Setting the time derivatives of system (21) to vanish, we obtain the equilibrium points for velocity and temperature fields asSolving 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 is sufficiently large. The other two equilibrium points are the solutions characterize the onset of finite amplitude steady motions, where .
5.2. 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 is governed by the roots of the following characteristic polynomial equation for the eigenvalues, :
Stability depends upon the value of ; for there is an exchange of stability, and for other two steady state solutions origin loses its stability. When , 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 (, , ) is governed by the roots of the following characteristic polynomial equation for the eigenvalues, :
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.
6. 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, and , have been chosen. The time domain is taken as with a constant time step, . 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 , , and with and have been projected onto planes. From Figure 1, we observed that for weak internal heating, that is, at , the motionless solution loses its stability and the convection solution takes over. For , 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 , 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 , 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 the well known routes to chaos, such as through period doubling, quasiperiodicity, or intermittency [43]. The behavior of the system at is chaotic which is confirmed by the broadening of the base in the power spectrum. For the moderate heating, that is, at , the motionless solution loses its stability and the convection solution takes over. For , 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 , it shows the spiral approach before it is attracted towards steady state fixed point. As the value of increased to , the flow becomes homoclinic in nature, and, on further increasing to , a transition to chaos occurs. In the case of strong heating, for , at , the motionless solution loses its stability and the convection solution takes over. For , 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 spiralling approach is even more pronounced, a fact which is consistent with the eigenvalues. From Figure 1 (), for , a homoclinic bifurcation occurs which is more pronounced in this case. Finally at , transition to chaos occurs. Figure 2 () is the plots of plane at , , and with . In Figure 2 (), we find qualitatively similar results as shown in Figure 1 (), except that amplitudes take a little bit larger value than for .
Figure 3 shows the initial supercritical solutions at fixed value of , , and with , projected onto planes. From Figure 3 () we observed that for weak heating, that is, at , the motionless solution loses stability and the convection solution takes over. For 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 () it is clear that there is spiral trajectory for as the flow exhibits an oscillatory decay. In Figure 3 (), spiralling approach of the trajectories towards the steady state fixed point is more pronounced. For , there is a homoclinic pattern of flow, as shown in Figure 3 (). Further increasing the value of , there occurs a transition to chaotic convection as can be seen from Figure 3 (). The transition to chaos in this case is similar to that leading to the Lorenz attractor; one can easily see “two-butterfly” Figure 3 () for . For the case of moderate heating, that is, at , 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 . At trajectories move via a spiral, before they are attracted to their steady state fixed points. At , the flow becomes homoclinic in nature and attracted towards a steady state fixed point. On increasing the value of at , a transition to chaos occurs and “two-butterfly” Figure 3 () are obtained. For the case of strong heating, namely, , the motionless solution loses its stability and the convection solution takes over at . For 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 , spiralling approach is even more pronounced, a fact which is consistent with the eigenvalues. From Figure 3 () for , a homoclinic bifurcation occurs which is more pronounced in this case. Finally at transition to chaos occurs and “two-butterfly” nature occurs. In Figure 4 (), which are for , we found qualitatively similar results to those of Figure 3 () with slightly large amplitudes.
The initial supercritical solutions, calculated at , , and with , projected onto planes have been presented in Figure 5. From Figure 5 (), we observe that for the case of weak heating, that is, at , the motionless solution loses its stability at , and the convection solution takes over. For 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 () it is clear that there are spiral trajectory for as the flow exhibits an oscillatory decay. In Figure 5 (), spiralical approach of the trajectories towards the steady state fixed point is more pronounced. For , there is a homoclinic bifurcation, as shown in Figure 5 (). Further increasing the value of , there occurs a transition to chaotic convection as can be seen from Figure 5 (). 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 . For the case of moderate heating, that is, at , 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 . At trajectories move via a spiral, before it is attracted to its steady state fixed point. At , the flow becomes homoclinic in nature and attracted towards a steady state fixed point. On increasing the value of at , a transition to chaos occurs. For the case of strong heating, namely, at , the motionless solution loses its stability and the convection solution takes over; for 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 , spiralling approach is even more pronounced, a fact which is consistent with the eigenvalues. From Figure 5 () for , a homoclinic bifurcation occurs which is more pronounced in this case. Finally at there occurs transition to chaos. Figure 6 () is the plots of initial supercritical solutions, and it is found that the results are qualitatively similar to those obtained in Figure 5 ().
Figure 7 is the plots of three-dimensional chaotic bifurcation at different values of and . Figures 8(a)–8(f) are the plots of time history of at different values of , at and with initial points for , for , and for . In Figures 9(a)–9(f), we plot the time history of at different values of , at and with initial points for , for , and for . Figures 10(a)–10(f) are the plots of time history of at different values of , at and with initial points for , for , and for .
(a) , , , and ; initial points (7, 7, 19, 0)
(b) , , , and ; initial points (7, 7, 19, 0)
(c) , , , and ; initial points (7.3, 7.3, 20, 0)
(d) , , , and ; initial points (7.3, 7.3, 20, 0)
(e) , , , and ; initial points (7.6, 7.6, 21.9, 0)
(f) , , , and ; initial points (7.6, 7.6, 21.9, 0)
(a) , , , and ; initial points (7, 7, 19, 0)
(b) , , , and ; initial points (7, 7, 19, 0)
(c) , , , and ; initial points (7.3, 7.3, 20, 0)
(d) , , , and ; initial points (7.3, 7.3, 20, 0)
(e) , , , and ; initial points (7.6, 7.6, 21.9, 0)
(f) , , , and ; initial points (7.6, 7.6, 21.9, 0)
(a) , , , and ; initial points (7, 7, 19, 0)
(b) , , , and ; initial points (7, 7, 19, 0)
(c) , , , and ; initial points (7.3, 7.3, 20, 0)
(d) , , , and ; initial points (7.3, 7.3, 20, 0)
(e) , , , and ; initial points (7.6, 7.6, 21.9, 0)
(f) , , , and ; initial points (7.6, 7.6, 21.9, 0)
7. Conclusions
In this paper, chaotic convection in a viscoelastic fluid saturated porous medium with internal heat generation is analyzed. A four-dimensional system of differential equations is deduced by a Fourier series expansion. The effects of internal heat source and viscoelastic parameters are investigated on dynamics of convection. The following conclusions are drawn:(1)An increment in the internal Rayleigh number accelerates the chaotic convection.(2)Effect of increasing retardation time is to increase the amplitude of convection thus accelerating the chaotic convection.(3)Effect of decreasing relaxation time is to accelerate the chaotic convection.(4)As particular case, the present system reduces to the systems of Lorenz [9], Akhatov and Chembarisova [30], Sheu et al. [31], and Khayat [26–28].
Nomenclature
Latin Symbols| Wave number | |
| Critical wave number | |
| Depth of the porous layer | |
| Darcy number, | |
| Prandtl number, | |
| Vadasz number, | |
| Permeability | |
| Pressure | |
| Velocity | |
| Internal heat source | |
| Rayleigh number | |
| Internal Rayleigh number, | |
| Time | |
| Temperature | |
| Temperature difference between the walls | |
| Space coordinate. |
| Relaxation time | |
| Retardation time | |
| Thermal expansion coefficient | |
| Nondimensional relaxation time | |
| Ratio of retardation time to relaxation time, | |
| Thermal diffusivity | |
| Porosity | |
| Density | |
| Dynamic viscosity | |
| Kinematic viscosity, | |
| Stream function. |
| Basic state | |
| Critical | |
| Reference state | |
| Unit normal vector in -direction | |
| Unit normal vector in -direction | |
| Unit normal vector in -direction | |
| , horizontal Laplacian | |
| . |
Conflict of Interests
The author declares that there is no conflict of interests regarding the publication of this paper.