An autonomous system for thermal convection of viscoelastic fluids in a porous layer using a thermal nonequilibrium model

Abstract

Thermal convection of viscoelastic fluids saturating a horizontal porous layer heated from below is analyzed using a thermal nonequilibrium model to take account of the interphase heat transfer between the fluid and the solid. The viscoelastic character of the flow is considered by a modified Darcy’s law. An autonomous system with five differential equations is deduced by applying the truncated Galerkin expansion to the momentum and heat transfer equations. The effects of interphase heat transfer H on the thermal convection of viscoelastic fluids in a porous medium are analyzed and discussed. The results show that the weak interphase heat transfer tends to stabilize the steady convection.

Keywords

thermal convection / porous media / viscoelastic fluid / thermal nonequilibrium model

Introduction

The thermal convection in a saturated porous media is a subject of practical interest for its applications in engineering, such as solar energy storage systems, geothermal reservoirs, passive cooling of nuclear reactors, pollutant transport in underground waters, soil decontamination, storage of chemical or agricultural products, etc. More detailed discussions of its theory and applications are presented by Nield and Bejan [1].

One interesting case of convection in porous media arises when the fluid is viscoelastic. Recently, Kim et al. [2] have studied thermal instability driven by buoyancy forces in an initially quiescent and horizontal porous layer saturated by viscoelastic fluids and found that the overstability is a preferred mode of instability for a certain range of elastic parameters. Yoon et al. [3] have considered the onset of oscillatory convection in a horizontal porous layer saturated with viscoelastic liquid. Bertola and Cafaro [4] have studied theoretically the instability of a viscoelastic fluid saturating a horizontal porous layer heated from below with a dynamical system approach.

However, the work of Bertola and Cafaro has been conducted under the assumption that the fluid and porous medium are in local thermodynamic equilibrium everywhere. However, for many practical applications involving high-speed flows or large temperature differences between the solid and fluid phases, the assumption of local thermal equilibrium is inadequate, and it is important to take account of the thermal nonequilibrium effects. Nield and Bejan [1] have stated that the thermal nonequilibrium model utilizes two equations to separately model the fluid and solid phases instead of using a single energy equation that describes the common temperature of the saturated medium (the one-field model). In the thermal nonequilibrium model, the two energy equations are coupled together by terms that account for the heat lost to or gained from the other phase. As a matter of fact, the thermal nonequilibrium effects on the convection of Newtonian fluids in porous media have been considered already. For example, Banu and Rees [5] have studied the effect of thermal nonequilibrium on the onset of convection in a porous layer and found that the critical Rayleigh number and wave number are modified by the presence of thermal nonequilibrium. Sheu [6] has studied the transition to chaos in a porous layer using thermal nonequilibrium, and the results show that the interphase heat transfer stabilizes steady convection and alters the routes to chaos.

The aim of the present paper is to study the thermal convection of a viscoelastic fluid saturating a horizontal porous layer using a thermal nonequilibrium model to take account of the interphase heat transfer between the fluid and the solid. The viscoelastic character of the flow is considered by a modified Darcy’s law. The truncated Galerkin expansion is applied to the governing equations of the thermal convection in a porous medium to deduce an autonomous system with five ordinary differential equations. The system is used to investigate the dynamic behavior of the thermal convection of viscoelastic fluids in the porous medium in order to study the effect of interphase heat transfer on the convection.

Formulations

Let us consider a fluid-saturated porous layer of width l and height d, which is heated from below and cooled from above. The bottom surface at y=\mathrm{0} is held at constant temperature T_{\mathrm{h}}, the top one at y=d is held at constant temperature T_{\mathrm{c}}, and the vertical walls are adiabatic.

The momentum equation on unsteady flows of a Newtonian fluid in porous media is usually expressed in the form of Darcy’s law, which states that

(1){\displaystyle \frac{{\rho }_{\mathrm{f}}K}{ϵ\mu }}{\displaystyle \frac{\partial \mathbf{u}}{\partial \overline{t}}}+\mathbf{u}=-{\displaystyle \frac{K}{\mu }}(\nabla p-{\rho }_{\mathrm{f}}\mathbf{g}),

where \mathbf{u}=\overline{u}\mathbf{i}+\overline{v}\mathbf{j} is the fluid velocity vector, \mathbf{g} is the gravity, p is the pressure of the fluid, {\rho }_{\mathrm{f}} is the density of the fluid, \mu is the dynamic viscosity of the fluid, K is a quantity called permeability, and ϵ is the porosity.

However, because the pressure drop in viscoelastic flows is nonlinearly related to the filtration velocity, the permeability would also depend on the relaxation time of the fluid. Following a well-established approach [7], the modified Darcy’s law can be obtained by including a relaxation term about the pressure gradient, with a characteristic time \tau depending on viscoelasticity:

(2){\displaystyle \frac{{\rho }_{\mathrm{f}}K}{ϵ\mu }}{\displaystyle \frac{\partial \mathbf{u}}{\partial \overline{t}}}+\mathbf{u}=-{\displaystyle \frac{K}{\mu }}\left(\mathrm{1}+\tau {\displaystyle \frac{\partial }{\partial \overline{t}}}\right)(\nabla p-{\rho }_{\mathrm{f}}\mathbf{g}).

This equation implicitly assumes that the fluid has a constant viscosity and a single relaxation time that is the case, for instance, of an upper convected Maxwell fluid [8].

It is assumed in this paper that the convective fluid and the porous medium are not in local thermodynamic equilibrium, and therefore, a two-temperature model of microscopic heat transfer applies. The governing equations for the fluid and solid temperatures are [1]

(3)ϵ{(\rho c)}_{\mathrm{f}}{\displaystyle \frac{\partial T_{\mathrm{f}}}{\partial \overline{t}}}+{(\rho c)}_{\mathrm{f}}\mathbf{u}·\nabla T_{\mathrm{f}}=ϵk_{\mathrm{f}}{\nabla }^{\mathrm{2}}{T}_{\mathrm{f}}+\overline{h}(T_{\mathrm{s}}-T_{\mathrm{f}}),

(4)(\mathrm{1}-ϵ){(\rho c)}_{\mathrm{s}}{\displaystyle \frac{\partial T_{\mathrm{s}}}{\partial \overline{t}}}=(\mathrm{1}-ϵ)k_{\mathrm{s}}{\nabla }^{\mathrm{2}}{T}_{\mathrm{s}}-\overline{h}(T_{\mathrm{s}}-T_{\mathrm{f}}),

and

(5)\rho ={\rho }_{\mathrm{0}}[\mathrm{1}-\beta (T_{\mathrm{f}}-T_{\mathrm{c}})],

where T is the temperature. The subscripts f and s denote the fluid and solid phases, respectively. The properties of the fluid and porous media include the specific heat (c), the coefficient of cubical expansion (\beta ), the thermal conductivity (k), and the interphase heat transfer coefficient between the fluid and the porous medium (\overline{h}).

Equations (2)-(5) and the continuity equation \nabla \cdot \mathbf{u}=\mathrm{0} constitute the basic equations on convection of a viscoelastic fluid saturating a porous layer. They can be nondimensionalized by using the following transformations:

(6)\{\begin{array}{l}\left(\overline{x},\overline{y}\right)=d(x,y),\hfill \\ \left(\overline{u},\overline{v}\right)={\displaystyle \frac{ϵk_{\mathrm{f}}}{{(\rho c)}_{\mathrm{f}}d}}(u,v),\hfill \\ T_{\mathrm{f}}=(T_{\mathrm{h}}-T_{\mathrm{c}})\theta +T_{\mathrm{c}},\hfill \\ T_{\mathrm{s}}=(T_{\mathrm{h}}-T_{\mathrm{c}})\varphi +T_{\mathrm{c}},\hfill \\ \overline{t}={\displaystyle \frac{{(\rho c)}_{\mathrm{f}}}{k_{\mathrm{f}}}}{d}^{\mathrm{2}}t.\hfill \end{array}

The pressure terms can be eliminated by introducing the stream function, u=-{\psi }_y and v={\psi }_x. In the stream function formulation, the continuity equation is automatically satisfied. Eqs. (2)-(5) then become

(7)\left({\tau }_{\mathrm{1}}{\displaystyle \frac{\partial }{\partial t}}+\mathrm{1}\right)({\psi }_{xx}+{\psi }_{yy})=Ra\left({\tau }_{\mathrm{2}}{\displaystyle \frac{\partial }{\partial t}}+\mathrm{1}\right){\theta }_x,

(8){\theta }_t-{\psi }_y{\theta }_x+{\psi }_x{\theta }_y={\theta }_{xx}+{\theta }_{yy}+h(\varphi -\theta ),

and

(9)\alpha {\varphi }_t={\varphi }_{xx}+{\varphi }_{yy}+\lambda h(\theta -\varphi ),

where

(10)\{\begin{array}{l}Ra={\displaystyle \frac{{\rho }_{\mathrm{f}}g\beta (T_{\mathrm{h}}-T_{\mathrm{c}})Kd}{ϵ\mu k_{\mathrm{f}}}},\hfill \\ \begin{array}{lll}\lambda ={\displaystyle \frac{ϵk_{\mathrm{f}}}{(\mathrm{1}-ϵ)k_{\mathrm{s}}}},& h={\displaystyle \frac{\overline{h}{d}^{\mathrm{2}}}{ϵk_{\mathrm{f}}}},& \alpha ={\displaystyle \frac{{(\rho c)}_{\mathrm{s}}{k}_{\mathrm{f}}}{{(\rho c)}_{\mathrm{f}}{k}_{\mathrm{s}}}}.\end{array}\hfill \end{array}

are Darcy-Rayleigh number (Ra), a porosity-modified conductivity ratio (\lambda ), a scaled interphase heat transfer coefficient (h), and a diffusivity ratio (\alpha ), respectively. In particular, {\tau }_{\mathrm{1}} is the dimensionless retardation time due to the action of the porous matrix, while {\tau }_{\mathrm{2}} is the dimensionless relaxation time depending on viscoelasticity.

Truncated Galerkin expansion

The finite amplitude analysis was carried out using a truncated representation of the Galerkin expansion by considering only one term for stream function and two terms for temperature distributions in the respective forms of

(11)\psi =A_{\mathrm{11}}\sin (ax)\sin (\mathrm{\pi }y),

(12)\theta =\mathrm{1}-y+B_{\mathrm{11}}\cos (ax)\sin (\mathrm{\pi }y)+B_{\mathrm{02}}\sin (\mathrm{2}\mathrm{\pi }y),

and

(13)\phi =\mathrm{1}-y+C_{\mathrm{11}}\cos (ax)\sin (\mathrm{\pi }y)+C_{\mathrm{02}}\sin (\mathrm{2}\mathrm{\pi }y).

Substituting Eqs. (11)-(13) into Eqs. (7)-(9), multiplying the equations by orthogonal eigenfunctions corresponding to (11)-(13), integrating them over the spatial domain, and rescaling the amplitudes yield the following system of nonlinear ordinary differential equations:

(14)\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}X}{\mathrm{d}\tau }}=-{\displaystyle \frac{\mathrm{1}-R{D}_{\mathrm{2}}}{D_{\mathrm{1}}}}Y-(R-\mathrm{1}){\displaystyle \frac{D_{\mathrm{2}}}{D_{\mathrm{1}}}}{\displaystyle \frac{H}{{\gamma }^{\mathrm{2}}}}(U-Y),\hfill \\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}\tau }}=RX-Y-(R-\mathrm{1})XZ+H(U-Y),\hfill \\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}\tau }}=\mathrm{4}\gamma (XY-Z)-H(W+Z),\hfill \\ {\displaystyle \frac{\mathrm{d}U}{\mathrm{d}\tau }}={\displaystyle \frac{\mathrm{1}}{\mathrm{\alpha }}}[-U+\lambda H(Y-U)],\hfill \\ {\displaystyle \frac{\mathrm{d}W}{\mathrm{d}\tau }}={\displaystyle \frac{\mathrm{1}}{\mathrm{\alpha }}}[-\mathrm{4}\gamma W-\lambda H(Z+W)],\hfill \end{array}

where the time was rescaled, and the following notation was introduced:

(15)\{\begin{array}{l}\begin{array}{lll}\tau ={\delta }^{\mathrm{2}}t,& D_{\mathrm{1}}={\delta }^{\mathrm{2}}{\tau }_{\mathrm{1}},& D_{\mathrm{2}}={\delta }^{\mathrm{2}}{\tau }_{\mathrm{2}},\end{array}\hfill \\ \begin{array}{lll}\gamma ={\displaystyle \frac{{\mathrm{\pi }}^{\mathrm{2}}}{{\delta }^{\mathrm{2}}}},& R={\displaystyle \frac{Ra}{{\delta }^{\mathrm{2}}}},& H={\displaystyle \frac{h}{{\delta }^{\mathrm{2}}}},\end{array}\hfill \end{array}

and the amplitudes were rescaled in the form of

(16)\{\begin{array}{l}X={\displaystyle \frac{A_{\mathrm{11}}a\gamma }{\mathrm{2}\mathrm{\pi }\sqrt{\mathrm{2}\gamma (R-\mathrm{1})}}},Y={\displaystyle \frac{\mathrm{\pi }R{B}_{\mathrm{11}}}{\mathrm{2}\sqrt{\mathrm{2}\gamma (R-\mathrm{1})}},}\hfill \\ Z=-{\displaystyle \frac{\mathrm{\pi }R{B}_{\mathrm{02}}}{(R-\mathrm{1})}},U={\displaystyle \frac{\mathrm{\pi }R{C}_{\mathrm{11}}}{\mathrm{2}\sqrt{\mathrm{2}\gamma (R-\mathrm{1})}}},\hfill \\ W={\displaystyle \frac{\mathrm{\pi }R{C}_{\mathrm{02}}}{(R-\mathrm{1})}}.\hfill \end{array}

System (14) provides a set of nonlinear ordinary differential equations with seven parameters. Parameter D_{\mathrm{1}} is a modified dimensionless retardation time, D_{\mathrm{2}} is a modified dimensionless relaxation time, R is the modified Darcy-Rayleigh number, and H is the modified interphase heat transfer coefficient; parameter \alpha is the diffusivity ratio, and \lambda is the porosity-modified conductivity ratio. The value of \gamma has to be consistent with the wave number at the convection threshold, which is required so that the convection cells fit into the domain and fulfill the boundary conditions. System (14) is similar to that given by Eqs. (13) in Sheu [6], the only differences being in the first equation that has two additional terms and different coefficients. System (14) at D_{\mathrm{2}}=\mathrm{0} are reduced to the system of Sheu [6].

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 lies below the critical values for instability to either a marginal state or an overstable infinitesimal perturbation.

The critical (or equilibrium) points correspond to the steady state solutions of the dynamical system (14), from which they are obtained by setting time derivatives at zero. Setting \mathbf{X}=[X,Y,Z,U,W], and {\mathbf{X}}_{\mathrm{1}}=[\mathrm{0},\mathrm{0},\mathrm{0},\mathrm{0},\mathrm{0}] is obviously an equilibrium point. Let

(17)\{\begin{array}{l}a={\displaystyle \frac{(\lambda H+\mathrm{1})(R-\mathrm{1})-[\mathrm{1}+(\mathrm{1}/{\gamma }^{\mathrm{2}}-\mathrm{1})R{D}_{\mathrm{2}}]H}{[(\lambda H+\mathrm{1})-(\mathrm{1}/{\gamma }^{\mathrm{2}}-\mathrm{1})D_{\mathrm{2}}H](R-\mathrm{1})}},\hfill \\ b={\displaystyle \frac{\lambda H+H+\mathrm{4}\gamma }{\lambda H+\mathrm{4}\gamma }},c={\displaystyle \frac{\lambda H}{\lambda H+\mathrm{1}}},\hfill \\ d=-{\displaystyle \frac{\lambda H}{\lambda H+\mathrm{4}\gamma }},e={\displaystyle \frac{(\lambda H+\mathrm{1})-(\mathrm{1}/{\gamma }^{\mathrm{2}}-\mathrm{1})D_{\mathrm{2}}H}{\lambda H+\mathrm{1}}}.\hfill \end{array}

Thus, two nonzero equilibrium points are given by

(18)\{\begin{array}{l}{\mathit{X}}_{\mathrm{2}}=\left[\sqrt{abe},\sqrt{{\displaystyle \frac{ab}{e}}},a,c\sqrt{{\displaystyle \frac{ab}{e}}},-da\right],\hfill \\ {\mathit{X}}_{\mathrm{3}}=\left[-\sqrt{abe},-\sqrt{{\displaystyle \frac{ab}{e}}},a,-c\sqrt{{\displaystyle \frac{ab}{e}}},-da\right].\hfill \end{array}

The solution {\mathit{X}}_{\mathrm{1}}=[\mathrm{0},\mathrm{0},\mathrm{0},\mathrm{0},\mathrm{0}] corresponds to pure conduction, which is known to be a possible solution though it is unstable when R(Ra) is sufficiently large. The remaining solutions {\mathit{X}}_{\mathrm{2}} and {\mathit{X}}_{\mathrm{3}} characterize the onset of finite amplitude steady motions.

Numerical results

By considering the thermal nonequilibrium model of heat transfer between the fluid and solid phase, system (14) was deduced to describe the dynamics of thermal convection of viscoelastic fluids in porous media. The objective of this study was to analyze the thermal nonequilibrium effect on the dynamic behavior of thermal convection of viscoelastic fluids in a porous medium as the Rayleigh number changes. The important parameters representing this effect are H, the interphase heat transfer coefficient, and k, the porosity-modified conductivity ratio. In order to simplify the analysis, the other parameters were kept constant. The value of \gamma used in all computations was 0.5, which is consistent with the critical wave number at marginal stability in porous medium convection. All solutions were obtained using the same initial conditions that were selected to be in the neighborhood of the positive convection fixed point, i.e., at \tau =\mathrm{0}: X=Y=Z=\mathrm{0.9}, and U=W=\mathrm{0}. System (14) was solved by applying Runge-Kutta method. All computations were carried out to a value of a maximum time, {\tau }_{\max }, of 100 with a constant time step \mathrm{\Delta }\tau =\mathrm{0.001}.

In the following discussions, three values of H (0.01, 0.05, 5) were chosen to investigate the effects of interphase heat transfer parameter H on the dynamic behaviors of system (14).

By carefully examining Figs. 1-3, several interesting and important physical phenomena can be found. For small values, e.g., H = 0.05, which indicates weak interphase heat transfer between the fluid and solid. As the Rayleigh number increases, a cascade of a period-2 (R=2), period-4 (R=5) route to multiperiods (R = 40) ensue. Furthermore, it can be predicted that chaos will occur as R increases.

Figures 4-6 show the phase portraits at H=0.01 for various values of the Rayleigh number (R=2, 5 and 40). Steady-state solutions that are different from the initial state can be gained at R=40, as shown in Fig. 6. The waveform of X reaches a stable value as the nondimensional time approaches 40, which means that the numerical solution is independent on the time. The phase trajectories of X-Z, X-W, and Y-U tend to be steady points. It is, however, interesting to observe in Figs. 4(a) and 5(a), that the waveform of X is periodically oscillatory.

The equilibrium model corresponding to temperature equilibrium between the fluid and solid suggests that the solution should lie on the plane Y = U and Z = -W. In Fig. 4, it can be observed that the solution departs from this plane. Z is 1 order of magnitude smaller than -W. In Fig. 5, it can be seen that W and U are four or three orders of magnitude smaller than –Z and Y.

The steady state solution is obtained at R=2 as H further increases to 5, as shown in Fig. 7. Figure 7(a) may be reviewed as the overdamped oscillation with the damping effect being weaker at lower Rayleigh number, and the flows tend to be stable pure conduction state: {\mathit{X}}_{\mathrm{1}}=[\mathrm{0},\mathrm{0},\mathrm{0},\mathrm{0},\mathrm{0}], which means that the numerical solution is independent on initial state \mathbf{X}=[\mathrm{0.9},\mathrm{0.9},\mathrm{0.9},\mathrm{0},\mathrm{0}].

As the Rayleigh number continues to increase to 5, Fig. 8(a) illustrates the oscillations of the time history of X. It is found that a periodic oscillation appears. Moreover, when a period-doubling route to chaos is found at R= 40, system (14) becomes chaotic. The destabilization occurs earlier when R increases. Figures 9(b), (c), and (d) give the chaotic trajectory of the equation. The temporal oscillating flows are quite periodic at lower R number, as shown in Fig. 8(a), but aperiodic at higher R number, as shown in Fig. 9(a).

It is noted in Fig. 8 that U and W are one order of magnitude smaller than Y and Z. In Fig. 9, it can be found that U and W are two or one order of magnitude smaller than Y and Z.

Conclusions

A five-dimensional autonomous dynamic system was deduced to analyze the thermal convection of viscoelastic fluids in a porous medium by applying the thermal nonequilibrium model. The problem was examined by means of stability of equilibria, time history, and phase portraits. It was found that the interphase heat transfer between the fluid and the solid alters the routes to chaos. Also, the very weak interphase heat transfer (H = 0.01) tends to stabilize steady convection. With weak interphase heat transfer (H = 0.05), as the Rayleigh number increases, a cascade of a period-2, period-4 route to multiperiods ensue. Furthermore, it can be predicted that chaos will occur as R increases. As interphase heat transfer H further increases to 5, the destabilization and chaos occurs earlier when R increases.

References

[1] Nield D A, Bejan A. Convection in Porous Media. 2nd ed. New York: Springer-Verlag, 1999

[2] Kim M C, Lee S B, Kim S, Chung B J. Thermal instability of viscoelastic fluids in porous media. International Journal of Heat and Mass Transfer, 2003, 46(26): 5065–5072

[3] Yoon D Y, Kim M C, Choi M C. The onset of oscillatory convection in horizontal porous layer saturated with viscoelastic liquid. Transport in Porous Media, 2004, 55(3): 275–284

[4] Bertola V, Cafaro E. Thermal instability of viscoelastic fluids in horizontal porous layersas initial value problem. International Journal of Heat and Mass Transfer, 2006, 49(21,22): 4003–4012

[5] Banu N, Rees D A S. Onset of Darcy-Benard convection using a thermal non-equilibrium model. International Journal of Heat and Mass Transfer, 2002, 45(11): 2221–2228

[6] Sheu L. An autonomous system for chaotic convection in a porous medium using a thermal non-equilibrium model. Chaos, Solitons and Fractals, 2006, 30(3): 672–689

[7] Alishaev M G, Mirzadjanzade A K. For the calculation of delay phenomenon in filtration theory. Izv Vuzov Neft I Gaz (Russ Chem Rev), 1975, 6(1): 71–77 (in Russian)

[8] Joseph D D. Fluid Dynamics of Viscoelastic Liquids. New York: Springer-Verlag, 1990

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg