An autonomous system for chaotic convection in a porous medium using a thermal non-equilibrium model

Chaos, Solitons and Fractals 30 (2006) 672–689 www.elsevier.com/locate/chaos

An autonomous system for chaotic convection in a porous medium using a thermal non-equilibrium model Long-Jye Sheu Department of Mechanical Engineering, Chung Hua University, No. 707 Wufu Road, Sec. 2, Hsinchu 30067, Taiwan Accepted 29 November 2005

Communicated by M.S. El Naschie

Abstract Thermal convection in a fluid-saturated porous medium was analyzed using a thermal non-equilibrium model to take account of the interphase heat transfer between the fluid and the solid. An autonomous system with five differential equations was deduced by applying the truncated Galerkin expansion to the momentum and heat transfer equations. The effects of the porosity-modified conductivity ratio, k, and interphase heat transfer, H, on the routes to chaos for convection in a porous medium were analyzed and discussed. The stability analysis revealed that the interphase heat transfer tends to stabilize the steady convection. The results show that with weak interphase heat transfer, a transition occurs from steady convection to chaos by a period-doubling sequence. However, an abrupt transition to chaos is predicted when interphase heat transfer is moderate and the porosity is small or moderate. As the product of k and H approaches infinity, the present system is reduced to Vadasz’s system which employs a thermal equilibrium model.  2005 Elsevier Ltd. All rights reserved.

1. Introduction The concept of chaos was first introduced by Poincare´ to describe orbits in space mechanics. In 1963, Lorenz [1] found chaos in a simple system of three autonomous ordinary differential equations with only quadratic non-linearity to describe the simplified Rayleigh–Be´nard problem. It may be that small differences in the initial conditions produce very great ones in the final phenomena. Since his work, chaos has been found in mechanical and electrical oscillators, in rotating heated fluids, in chemical reactions, in quantum physics, etc. [2–4]. Chaos and non-linear dynamics have provided new theoretical and conceptual tools that allow us to understand the complex behaviors of systems in essentially every field of science. Various theoretical tools that allow the examination of non-linear phenomena are used in the present study to examine the effect of interphase heat transfer on the thermal convection in a porous medium. The problem of flow and heat transfer in a fluid-saturated porous domain has been the subject of broad interest due to its wide application in different fields such as geothermal energy utilization, oil reservoir modeling, building thermal insulation, and nuclear waste disposal. Comprehensive reviews of applications and the fundamentals of heat transfer in

E-mail address: [email protected] 0960-0779/$ - see front matter  2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.chaos.2005.11.080

L.-J. Sheu / Chaos, Solitons and Fractals 30 (2006) 672–689

673

porous medium are presented by Nield and Bejan [5] and Igham and Pop [6]. For fluid in the porous layer, the problem consists of a fluid layer between two plates, with the lower plate kept at a higher temperature than the upper one. Once the temperature difference exceeds a critical value, the buoyancy causing convection becomes larger than the stabilizing viscous forces, and the fluid begins to move. The Rayleigh number corresponding to the onset of convection is often termed the critical Rayleigh number. Analytical solutions for the linear stability of free convection in the porous layer were presented by Lapwood [7]. Linear stability analysis provides the stability criteria, i.e., the critical Rayleigh number and the critical wave number. However, linear stability analysis cannot provide information on the values of the convection amplitudes nor on the onset of aperiodic or chaotic motions. A comprehensive review paper on non-linear convection in a porous medium was presented by Masuoka et al. [8]. Vadasz and coworkers [9–14] performed a series of investigations on the transition to chaos behaviors in a fluid-saturated porous layer heated from below and subject to gravity. A truncated Galerkin representation of the natural convection problem in a porous layer leads to the familiar Lorenz equations for the evolution of the convection amplitudes with time. In Vadasz’s work, a new dimensionless number that governs the effects of porosity, the Prandtl number, and the Darcy number on flow in the porous medium were defined. The dimensionless number was named the Vadasz number (Va) by Straughan [15]. The works of Vadasz et al. show that for convection with low Prandtl numbers in porous media, the transition from steady convection to chaos is sudden and occurs via a Hopf bifurcation which produces a ‘‘solitary limit cycle’’ and which may be associated with a homoclinic explosion. For moderate Prandtl numbers, the solution shows a transition to chaos via a period-doubling sequence of bifurcations at the Rayleigh number far beyond the critical value associated with the loss of stability of the convection steady solution. The work of Vadasz and coworkers, however, has been investigated under the assumption that the fluid and porous medium are in local thermodynamic equilibrium (LTE) everywhere. However, in many practical applications, the solid and fluid phases are not in thermal equilibrium. Nield and Bejan [5] stated the equations which are generally regarded as modeling unsteady heat transfer in a saturated porous medium where LTE does not apply. Instead of having a single energy equation which describes the common temperature of the saturated medium (the one-field model), the thermal non-equilibrium model utilities two equations to separately model the fluid and solid phases. In the thermal non-equilibrium model, the energy equations are coupled together by terms which account for the heat lost to or gained from the other phase. Rees and coworkers [16–19] in a series of studies investigated thermal non-equilibrium effects on free convective flows in a porous medium. A review by Kuznetsov [20] gives detailed information about the work on thermal non-equilibrium effects. Recently, Banu and Rees [19] studied the effect of thermal non-equilibrium on the onset of convection in a porous layer. It was found that the critical Rayleigh number and wave number are modified by the presence of thermal non-equilibrium. The aim of the present paper was to study the transition to chaos in a porous layer using the non-equilibrium model of microscopic heat transfer between the fluid and solid phases in the porous medium. The truncated Galerkin expansion was applied to the governing equations of thermal convection in a porous medium to deduce an autonomous system with five ordinary differential equations. The system was used to investigate the dynamic behavior of the thermal convection in the porous medium in order to study the effect of interphase heat transfer on the route to chaos. Results show that interphase heat transfer stabilizes steady convection and alters the routes to chaos.

2. Formulations We consider a Boussinesq fluid-saturated porous layer of width l and height d, which is heated from below and cooled from above. The lower surface is held at temperature Tl, while the upper surface is at Tu, and the vertical walls are adiabatic. The Darcy model is employed for the momentum equation. The time derivative term is retained in Darcy’s equation for the conservation of momentum to account for the effect of the Prandtl number in a porous medium. Vadasz and Olem [9,12–14] showed the effect and necessity of retaining this term when investigating the dynamic behavior of convection in a porous medium. Under these assumptions, the respective continuity and momentum equations are [21] r  q ¼ 0; and   K q o q¼ rP  qf g þ f q ; e ot l

ð1Þ ð2Þ

where q ¼ ui þ vj is the fluid velocity vector, g is the gravity, P is the pressure, K is the permeability, l is the fluid viscosity, e is the porosity, and q is the density.

674

L.-J. Sheu / Chaos, Solitons and Fractals 30 (2006) 672–689

We assume 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 [5] oT f þ ðqcÞf q  rT f ¼ ek f r2 T f þ hðT s  T f Þ; ot oT s ¼ ð1  eÞk s r2 T s  hðT s  T f Þ; and ð1  eÞðqcÞs ot q ¼ q0 ½1  bðT f  T u Þ; eðqcÞf

ð3Þ ð4Þ ð5Þ

where T is the temperature. The f and s subscripts 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 (b), the thermal conductivity (k), and the interphase heat transfer coefficient ðhÞ. Eqs. (1)–(5) can be non-dimensionalized by using the following transformations: ðx; y Þ ¼ dðx; yÞ; ðu; vÞ ¼ T s ¼ ðT l  T u Þ/ þ T u ;

ek f ðu; vÞ; ðqcÞf d t ¼

T f ¼ ðT l  T u Þh þ T u ; ð6Þ

ðqcÞf 2 d t kf

The pressure terms can be eliminated by introducing the stream function, u = wy and v = wx. Eqs. (1)–(5) then become   1 o þ 1 ðwxx þ wyy Þ ¼ Rahx ; ð7Þ Va ot ht  wy hx þ wx hy ¼ hxx þ hyy þ hð/  hÞ;

and

a/t ¼ /xx þ /yy þ khðh  /Þ;

ð8Þ ð9Þ

where Va ¼

ePr ; Da

hd 2 h¼ ; ek f

Ra ¼

qf gbðT l  T u ÞKd ; eljf

ðqcÞs k f jf a¼ ¼ ðqcÞf k s js

k¼

ek f ; ð1  eÞk s

ð10Þ

are the Vadasz number (Va), Darcy–Rayleigh number (Ra), a porosity-modified conductivity ratio (k), a scaled interphase heat transfer coefficient (h), and a diffusivity ratio (a), respectively.

3. 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 the temperature distributions of the fluid and porous medium. We represent the stream function, and fluid and porous medium temperature distributions in the respective forms of px w ¼ A11 sin sinðpyÞ; ð11Þ L px sinðpyÞ þ B02 sinð2pyÞ; and ð12Þ h ¼ 1  y þ B11 cos L  px sinðpyÞ þ C 02 sinð2pyÞ. ð13Þ / ¼ 1  y þ C 11 cos L Substituting Eqs. (11)–(13) into Eqs. (7)–(9), multiplying the equations by orthogonal eigenfunctions corresponding to (11)–(13), and then integrating them over the spatial domain yields a set of five differential equations for the time evolution of the amplitudes, in the form of   dA11 Vac Ra B11  A11 ; ¼ 2 ð14Þ p pv ds dB11 1 1 h ¼ B11 þ A11 þ A11 B02 þ 2 ðC 11  B11 Þ; ð15Þ pv v cp ds

L.-J. Sheu / Chaos, Solitons and Fractals 30 (2006) 672–689

dB02 1 h ¼  A11 B11  4cB02 þ 2 ðC 02  B02 Þ; 2v cp ds   dC 11 1 kh ¼ C 11 þ 2 ðB11  C 11 Þ ; and a cp ds   dC 02 1 kh ¼ 4cC 02 þ 2 ðB02  C 02 Þ ; ds a cp

675

ð16Þ ð17Þ ð18Þ

where the time was rescaled and the following notation was introduced: s¼

ðL2 þ 1Þp2 t; L2

v¼

L2 þ 1 ; L

c¼

L2 . L þ1

ð19Þ

2

Although the relationship between the solutions of the governing partial differential and the corresponding truncated ordinary differential systems has not been established, these lower-order spectral models may qualitatively reproduce the convective phenomenon observed in the full system. The results can also be used as the starting values when discussing the fully non-linear problem. It is convenient to introduce the following further notation: R¼

Ra ; p2 v2

r¼

Vac ; p2

H¼

ch ; p2

ð20Þ

and rescale the amplitudes in the form of A11 pRB11 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ; Y ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi; 2v 2cðR  1Þ 2 2cðR  1Þ pRC 11 pRC 02 U ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; W ¼ ðR  1Þ 2 2cðR  1Þ X ¼

Z¼

pRB02 ; ðR  1Þ

ð21Þ

to provide the following set of equations: X_ ¼ rðY  X Þ; Y_ ¼ RX  Y  ðR  1ÞXZ þ H ðU  Y Þ; Z_ ¼ 4cðXY  ZÞ  H ðW þ ZÞ;

ð22Þ

1 U_ ¼ ½U þ kH ðY  U Þ; a 1 W_ ¼ ½4cW  kH ðZ þ W Þ. a

12 λ

10

0 0.01 0.1 1.0 10.0 100.0

Rc1

8 6 4 2 0 0

2

4

H

6

8

10

Fig. 1. The critical Rayleigh number Rc1 as a function of H for specific values of parameter c.

676

L.-J. Sheu / Chaos, Solitons and Fractals 30 (2006) 672–689

System (22) provides a set of non-linear equations with six parameters. Parameter r is a modified Vadasz number, R is a modified Darcy–Rayleigh number, H is the modified interphase heat transfer coefficient between the fluid and porous medium, parameter a is the diffusivity ratio of the fluid and porous medium and k is the porosity-modified conductivity ratio. The value of c 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. This requirement yields a value of c = 0.5 at H = 0 for porous medium convection. We obtained a non-linear autonomous system for the five unknowns: X(t), Y(t), Z(t), U(t), and W(t). System (22) can be reduced to the system provided by Vadasz and Olek [9] at two limiting conditions, H = 0 and kH ! 1. The first condition, H = 0, is provided by a pure mathematical sense. The equations of X(t), Y(t), and Z(t) are decoupled from the equations of Y(t) and W(t) at H = 0. Hence, the first three equations at H = 0 are reduced to the system of Vadasz

Fig. 2. The evolution of the complex eigenvalues with increasing the Rayleigh number, for r = 50, c = 0.5, a = 1.0, k = 1.0 and H = 1.0.

800

λ 0 0.01 0.1 1.0 10.0 100.0

Rc2

600

400

200

0 0

2

4

6

8

10

H Fig. 3. The critical Rayleigh number Rc2 as a function of H for specific values of parameter c.

L.-J. Sheu / Chaos, Solitons and Fractals 30 (2006) 672–689

677

and Olek. With simple transformations, it can be shown that the reduced system is equivalent to the Lorenz system [12]. The second condition, kH ! 1, is provided by both mathematical and physical senses. Mathematically, as kH

Fig. 4. Bifurcation diagrams and Lyapunov spectrum as a function of Rayleigh number: (a) H = 0 and (b) H = 0.2, k = 0.1.

678

L.-J. Sheu / Chaos, Solitons and Fractals 30 (2006) 672–689

approaches infinity, the fourth and fifth equations in system (22) imply the equalities: Y = U and Z = W. Therefore, the last terms in the second and third equations of system (22) are eliminated. Physically, kH approaching infinity means a strong interphase heat transfer between the fluid and solid which achieves equality in the temperatures of the fluid and solid. Hence, the assumption of the thermal equilibrium model by Vadasz and Olem is preserved. A numerical example of thermal equilibrium as kH ! 1 is provided in Section 5.

4. Stability analysis Before proceeding to the numerical solution of system (22), it would be useful to carry out some local stability analyses in an attempt to unravel some of the fundamental differences between thermal equilibrium and non-equilibrium models. System (22) has the following basic properties.

1.3 1.2

Z

1.1 1 0.9 R

0.8

193 194

0.7 -3

-2

-1

0 X (a)

1

2

3

0.4

R

Y(nT)

0.2

193 194

0

-0.2

-0.4 -0.4

-0.2

0 X(nT) (b)

0.2

0.4

Fig. 5. (a) Phase portraits and (b) Poincare section at Z = 1 for R = 193 and 194, H = 0.2, k = 0.1.

L.-J. Sheu / Chaos, Solitons and Fractals 30 (2006) 672–689

679

(1) Dissipation: Since oX_ oY_ oZ_ oU_ oW_ þ þ þ þ oX oY oZ oU oW ¼ ½r þ 2H þ 4c þ 1 þ ð1 þ 2kH þ 4cÞ=a

rV ¼

ð23Þ

with all physically meaningful parameters non-negative, the system is dissipative. (2) Equilibria: The equilibria of the system can be found by solving the following system of equations: rðY  X Þ ¼ 0; RX  Y  ðR  1ÞXZ þ H ðU  Y Þ ¼ 0; 4cðXY  ZÞ  H ðW þ ZÞ ¼ 0; 1 ½U þ kH ðY  U Þ ¼ 0; and a 1 ½4cW  kH ðZ þ W Þ ¼ 0. a

ð24Þ ð25Þ ð26Þ ð27Þ ð28Þ

1.3

-0.0094

0.002

0.2

1.2

-0.0095

0.001

0.1

W

-0.0096 -0.0097

0.9

-0.002

0.8

-0.0098

-0.003

0.7

-0.0099 -3

-0.004

-3

-2

-1

0

X

1

2

3

-2

-1

0

1

2

3

X

-0.1 -0.2 -0.3

-4

-2

0

2

-0.4 -0.4 -0.3 -0.2 -0.1

4

(a)

0

0.1

X(nT)

Y

-0.087

1.3

0

-0.001

U

1

Z

0

Y(nT)

1.1

0.2

0.002 0.001

1.2 -0.088

U

W

Z

-0.089

1

Y(nT)

0

1.1

-0.001

0.9

-0.002

0

-0.2

-0.09 0.8

-0.003 -0.091 -3

-2

-1

0

1

2

3

-0.004 -3

X

-2

-1

0

1

2

3

X

-4

(b)

-0.0094

0.04

1.2

-0.0095

0.03

U

W

Z

0.7 -3

-2

-1

0

1

2

3

-0.0099 -3

0.1

-0.2

-0.02 -2

-1

0

1

2

3

-4

X

-2

0

-0.4 -0.4 -0.3 -0.2 -0.1

4

2

Y

(c) -0.087

X(nT)

0.004

1.2

0.4

0.003 -0.088

0.2

1.1 -0.089

0.9

Y(nT)

1

U

0.002

W

Z

0

0

-0.01

X

1.3

0.1

0.2

0

-0.0098

0.8

0

X(nT)

0.01

-0.0097

0.9

-0.4 -0.4 -0.3 -0.2 -0.1

4

2

0.02

-0.0096

1

0

Y

1.3

1.1

-2

Y(nT)

0.7

0.001

0

0 -0.09

-0.2

0.8

-0.001

0.7 -3

-2

-1

0

X

1

2

3

-0.091 -3

-2

-1

0

X

1

2

-0.002 -4

3

(d)

-2

0

Y

2

4

-0.4 -0.2

0

0.2

0.4

X(nT)

Fig. 6. Phase portraits and Poincare´ sections during period doubling at H = 0.2, k = 0.1: (a) R = 220, period-2; (b) R = 240, period-4; (c) R = 248, period-8 and (d) R = 260, chaos.

680

L.-J. Sheu / Chaos, Solitons and Fractals 30 (2006) 672–689

Fig. 7. Phase portraits and Poincare´ sections at H = 0.2, k = 0.1. (a) R = 261, (b) R = 284, (c) R = 296, (d) R = 325, (e) R = 374 and (f) R = 409.

L.-J. Sheu / Chaos, Solitons and Fractals 30 (2006) 672–689

681

Obviously, S1 = (0, 0, 0, 0, 0) is an equilibrium. Let a¼

ðkH þ 1ÞðR  1Þ  H ; ðkH þ 1ÞðR  1Þ

b¼

kH þ H þ 4c ; kH þ 4c

c¼

kH ; kH þ 1

d ¼

kH kH þ 4c

thus, other non-zero equilibria are given by   pffiffiffiffiffi pffiffiffiffiffi pffiffiffiffiffi d S 2;3 ¼  ab;  ab; a; c ab;  a . b

ð29Þ

ð30Þ

The system has three equilibria. When ab > 0, the three equilibria are all real. For convection in porous media, b is always >0. Thus, when (R  1) > 0 and (kH + 1)(R  1)  H > 0, the three equilibria are all real. The fixed point S1 corresponds to the motionless solution, and S2,3 corresponds to the convection solution. (3) Stability of equilibria: By linearizing system (22), one obtains its Jacobian matrix as follows: 2 3 r r 0 0 0 6 7 H 0 6 R  ðR  1ÞZ ð1 þ H Þ ðR  1ÞX 7 6 7 6 7. 4cY 4cX ð4c þ H Þ 0 H J¼6 ð31Þ 7 6 7 0 kH =a 0 ðkH þ 1Þ=a 0 4 5 0 0 kH =a 0 ð4c þ kH Þ=a Since the matrix is 5 · 5, it is hard to obtain the eigenvalues in a closed form. However, numerical calculations can be performed to discuss the stability at the equilibria. At fixed point S1, the motionless solution loses stability, and the convection solution takes over at the critical value of Rc1. The value of Rc1, which corresponds to the onset of convection, is solved numerically for various values of H and k in Fig. 1 with the values of parameters r = 50 and c = 0.5. The results agree with the linear stability analysis of Buna and Rees [19]. For a wider range of values of H and c, the reader can refer to the work of Bunu and Rees [19]. The stability of fixed point S2,3 is associated with the convection solutions. For r = 50, c = 0.5, k = 1, and H = 1, the evolution of the complex eigenvalues of J is shown in Fig. 2. As the value of

Fig. 8. (a) Bifurcation diagrams and (b) Lyapunov spectrum as a function of Rayleigh number R, H = 2, k = 0.1.

682

L.-J. Sheu / Chaos, Solitons and Fractals 30 (2006) 672–689

R increases, two roots become complex-conjugate at R ffi 1.97. However, they still have negative real parts; therefore the convection fixed points are stable, and it is a spiral node. As the value of R increases further, the real part of the root becomes non-negative at Rc2 = 83.67. In the neighborhood of this point, the convection fixed points lose their stability and another periodic, quasi-periodic, or chaotic solution may become predominant. The influence of H on Rc2 for various values of k is shown in Fig. 3. In the thermal equilibrium model, the value of Rc2 is 58.51. At small values of H, we see from Fig. 3 that Rc2 is close to 58.51 and is independent of k. The physical reason for this is that there is almost no transfer of heat between the phases, and therefore the onset criterion is not affected by the properties of the solid phase. It is obvious that the value of Rc2 increases with H when k is small. A small value of k corresponds to small porosity or small fluid conductivity. For a moderate value of k (k = 10), the value of Rc2 initially increases with H but eventually begins to decrease after reaching a maximum value of Rc2 = 119.7 at H = 2.3. It is obvious that Rc2 is independent of H for the large value of k = 100. In this case, the porosity approaches 1, and interphase heat transfer is diminished.

5. Numerical results By considering the thermal non-equilibrium model of heat transfer between the fluid and solid phase, we deduced system (22) to describe the dynamics of thermal convection in a fluid-saturated porous medium. The objective of this study was to analyze the thermal non-equilibrium effect on the dynamic behavior of thermal convection in a porous medium as the Rayleigh number changes. According to Fig. 3, the important parameters representing this effect are H, the inter-phase heat transfer coefficient, and k, the modified conductivity ratio. In order to simplify the analysis, the other parameters were kept constant. For convection with moderate Prandtl numbers, we considered a value of r = 50 which is consistent with a Vadasz number of Va  987. The value of c used in the computations was chosen to be 0.5, which is consistent with the critical wave number at the marginal stability in porous media convection. The other parameter, a, was chosen to be a = 1 in all computations if not specified. The initial conditions were selected to be in the neighborhood of the positive convection fixed point, i.e., at t = 0: X = Y = Z = 0.9 and U = W = 0. Vadasz

Fig. 9. (a) Bifurcation diagram and (b) Lyapunov spectrum as a function of Rayleigh number R, H = 5.0, k = 0.1.

L.-J. Sheu / Chaos, Solitons and Fractals 30 (2006) 672–689

683

and Olek [13] thoroughly studied the corresponding problem for the parameters chosen in the present study by applying the local thermal equilibrium model. System (22) is solved by applying a fifth- and sixth-order Runge–Kutta–Verner method from the IMSL Library (DIVPRK). All computations were carried out to a value of maximum time smax = 200 with a constant time step Ds = 0.001. In the following discussions, we chose three values of k = 0.1, 10, and 100 to investigate the effects of interphase heat transfer parameter H on the dynamic behaviors of system (22). The three values of k = 0.1, 10, and 100 can respectively be regarded as low, moderate, and large porosity-modified conductivity ratios. 5.1. Low porosity-modified conductivity ratio, k = 0.1 The bifurcation scenarios and Lyapunov spectrum characterized by the parameters k = 0.1 and H = 0 and 0.2 are illustrated in Fig. 4 for the interval 10 6 R 6 600 with a resolution of DR = 1. The Lyapunov exponents were calculated using a method developed by Wolf et al. [22] to identify the chaotic regimes. It can be seen from the Lyapunov spectrum that the first transition to chaotic motion occurred at about R  250 when H = 0.2. However, when H = 0, the first transition to chaos occurred at about R  290. While Vadasz and Olek [13] used peaks and valleys for the bifurcation diagram, Fig. 4(a) and (b) were plotted using the Poincare´ section at Z = 1. Bifurcation diagrams based on this Poincare´ section reveal more physical scenarios based on the peaks and valleys. For example, the bifurcation diagram identifies routes to chaos and the number of periods. In addition, a sudden change in the bifurcation diagram indicates a jump to symmetrically similar phase portraits in shape and location as shown in Fig. 5. It is clear from Fig. 4b that there is a discontinuity when R is 193 and H = 0.2. Fig. 5 shows phase portraits and the Poincare´ section at Z = 1 for R = 193

6

4

X

2

0

-2

-4

-6 0

40

80

120

(a) 1.6

0.8

0 -0.04

1.2

0.4

U

W

Z

-0.08 0.8

0

-0.12 0.4

-0.4

-0.16

0

-0.8

-0.2 -6

-4

-2

0

2

4

6

-6

-4

-2

X

0

X

2

4

6

-6

-4

-2

0

Y

(b) Fig. 10. A transient chaos at R = 122, H = 5, k = 0.1: (a) time history and (b) phase portraits.

2

4

6

684

L.-J. Sheu / Chaos, Solitons and Fractals 30 (2006) 672–689

and 194. Hence, from the bifurcation diagram, the jump between attractors can be realized. By carefully examining Fig. 4, one can find several interesting and important physical phenomena. For small values, e.g., H = 0.2, which indicate weak interphase heat transfer between the fluid and solid, bifurcation from steady convection to period-2 motion was observed in the neighborhood of the critical value of Rc2  60.02. As the Rayleigh number increased, a cascade of a period-doubling route to chaos then ensued. Furthermore, state y shows a band-like shape (a sudden jump in the bifurcation diagram) and eventually becomes wider as R increases. The position of the band-like shape of y implies a jump of symmetrically similar phase portraits in shape and location during the period-doubling route to chaos. Fig. 6 shows the phase portraits and Poincare´ sections at periods 2, 4, and 8 and chaotic motion for various values of the Rayleigh number at H = 0.2 and k = 0.1. One can observe the effect of the non-equilibrium model on the solution in Fig. 6. 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. From Fig. 6, one can observe that the solution departs from this plane. The magnitudes of U and W are 1- or 2-orders of magnitude smaller than those of Y and Z. By comparing the bifurcations of

Fig. 11. Bifurcation diagrams as a function of Rayleigh number R at k = 10: (a) H = 0.2, (b) H = 2.0 and (c) H = 10.0.

L.-J. Sheu / Chaos, Solitons and Fractals 30 (2006) 672–689

685

H = 0 and H = 0.2, one can see that the onset of convection occurs at a larger value of the Rayleigh number while chaos appears at a smaller value of the Rayleigh number. Thus we may say that the non-equilibrium model with weak interphase heat transfer and a low porosity-modified conductivity ratio stabilizes the steady convection while precipitating the onset of chaos. A further examination of Fig. 4(a) reveals as pointed out by Vadasz and Olek [13], possible periodic windows around R  300, 338, 375, 400, 437 and 453 and a transition to a wide periodic window regime above R  535 when H = 0.

(a)

(b)

(c)

(d)

(e) Fig. 12. Phase portraits and Poincare´ sections during inverse period doubling at k = 10, H = 10: (a) R = 78, chaos; (b) R = 81.5, period-16; (c) R = 81.7, period-8; (d) R = 82, period-4 and (e) R = 85, period-2.

686

L.-J. Sheu / Chaos, Solitons and Fractals 30 (2006) 672–689

These windows are also shown in Fig. 4(a) when the largest Lyapunov exponent drops to 0 which indicates a periodic window. With a weak non-equilibrium effect, i.e., H = 0.2, the corresponding periodic windows were also found at around R  261, 284, 296, 325, 374 and 409. The wide periodic regime occurred at around R  470. The interphase heat transfer effects appeared to precipitate all of the periodic windows. The phase portraits and Poincare´ sections of these windows are shown in Fig. 7. The first periodic window was identified as period 6 at R = 261 by the Poincare´ section. The dynamic behavior became chaotic with an increase in the Rayleigh number, and then a periodic window was found at R = 284. The Poincare´ section shows this periodic window to be period 10. The third periodic window appeared at R = 296 and was identified as a period-6 window. The Poincare´ section demonstrated the fourth periodic window to be period 3 at R = 325. By increasing the Rayleigh number further, another period-3 periodic window was found at R = 374. Finally, the periodic window was identified to be period 8 at R = 409. In the broad range of Rayleigh numbers, 261 < R < 409, two of the windows with period 6 were found at R = 261 and 296. It is evident that the shapes of

Fig. 13. Bifurcation diagrams as a function of Rayleigh number at k = 100: (a) H = 0.2, (b) H = 5.0 and (c) H = 100.0.

L.-J. Sheu / Chaos, Solitons and Fractals 30 (2006) 672–689

687

the attractors differ. An attractor which is non-symmetrical in the X–Z plane was found at R = 261, while the period-6 attractor appeared symmetrical in the X–Z plane at R = 296. As R increased to 470, a wide period-4 window was found. The bifurcation diagram and Lyapunov spectrum are shown in Fig. 8(a) and (b) for H = 2.0. When H increased to 2.0, some windows disappeared under the resolution of DR = 1 as can be seen in Fig. 8. The wide window with period 4 was not found. In this case, the route to chaos totally differed from that when H = 0.2. The period-doubling route to chaos was not found. The stability analysis of fixed point S2,3 predicted that the real part of the pair of complex conjugate eigenvalues became positive at R = 81.63. In the neighborhood of this point, at R = 79, a sudden transition from steady convection to chaos was found. As R further increased, two period-8 windows were found at R = 143 and 170. Fig. 9 shows the bifurcation diagram and Lyapunov spectrum when H increased to 5.0 at k = 0.1. A sudden transition from steady convection to chaos was found again at R = 123. As R further increased, no periodic window appeared. Fig. 10 shows the time history of X and phase portraits at R = 122. It can be seen that a transient chaos was found at R = 122 before the transition to chaos. 0.3

1.2

2

-0.9

0.2

1.1

1

-1

0.1

1

0

0.9

-1

0.8

-2

0.7 -3

-2

-1

0 X

1

2

W -1.1 -1.2

-3 -3

3

-2

-1

0 Y

1

2

3

(a)

2

1.2

-0.2 -0.1

1.3

0.4

-0.9

0.3

0

0.1 X(nT)

0.2

0.3

0.2 Y(nT)

W

U

-1.1

0.1 0

-2

-1.2

0.8 0.7

-1.3 0.8

-4 -3

-2

-1

0 X

1

2

-3

3

-2

-1

0 Y

1

2

3

1

0.9

(b)

1.1

1.2

1.2

Z

-4 0 X

1

2

1.3

-2

-1

0 Y

1

2

3

0.8

(c)

0.9

1

1.1

1.2

1.2

0

0.3

0.4

0.4

-0.9

2 1.1

0.2

-1 W

U

-0.2 -0.1

1.3

Z -0.8

4

1

0.1 0.2 X(nT)

0

-1.3 -3

3

0.4

-1.2

0.8

-1

0.3

Y( nT)

-1.1 -2

-2

0.2

0.2

W

U

0

0.9

-3

0.1 X(nT)

0.4

-1

0.7

0

-0.9

2 1.1 1

-0.2 -0.1

1.3

-0.8

4

1.3

-0.1

Y( nT)

Z

1.1

-1 0

0.9

Z

1

-0.8

1.1

Z

0.9

Z

1.2

1

0 -0.1

-1.3 0.8

4

1.3

Y( nT)

-0.8

U

3

Z

1.3

0

-1.1

0.9

0 -2

-1.2

0.8 0.7

-4 -3

-2

-1

0 X

1

2

3

-1.3 -3

-2

-1

0 Y

1

2

3

(d)

0.8

0.9

1

1.1 Z

1.2

1.3

-0.2 -0.2 -0.1 0

0.1 0.2 0.3 X(nT)

0.4

Fig. 14. Phase portraits and Poincare´ sections at k = 100, H = 100: (a) R = 220, period-2; (b) R = 275, period-4; (c) R = 287, period-8 and (d) R = 295, chaos.

688

L.-J. Sheu / Chaos, Solitons and Fractals 30 (2006) 672–689

5.2. Moderate porosity-modified conductivity ratio, k = 10 Fig. 11 shows the bifurcation diagrams as a function of R at k = 10 and H = 0.2, 2.0, and 10.0. When H = 0.2, the first transition from steady convection to periodic motion occurred at R  65. Then, a period-doubling route to chaos was found. At R  264, system (22) became chaotic. The bifurcation scenarios of Y is similar to that in Fig. 4(b) where k = 0.1 and H = 0.2. That is, with weak interphase heat transfer, the modified conductivity ratio does not affect the behavior of the route to chaos. As H further increased to 2.0, as shown in Fig. 11(b), a sudden transition from steady convection to chaos was observed at around R  105. As the Rayleigh number increased, state y shows an band-like shape which indicates a jump of symmetrically similar attractors. The wide period-4 window appeared at R  336. Fig. 11(c) shows the bifurcation of Y at H = 10. It can be seen that a transition from steady convection to period-4 motion occurred at R = 74. Then a short period of chaotic motion followed. As the Rayleigh number continued to increase, an inverse period-doubling bifurcation associated with the transition from chaos towards period-2 motion was observed. Fig. 12 shows details of the phase portraits and Poincare´ sections during the inverse period-doubling bifurcation. It was noted that U and W are one order of magnitude larger than U and W when k = 0.1. From Fig. 12, one finds the orders of magnitudes of U and W are now equivalent to those Y and Z. After the inverse period-doubling bifurcation, a period-doubling route to chaos was found as R further increased. The details of the period-doubling route to chaos were similar to those for the case of k = 0.1 and H = 0.2 and are not addressed here. 5.3. Large porosity-modified conductivity ratio, k = 100 Fig. 13 shows the bifurcation diagrams as a function of Rayleigh number at H = 0.2, 5.0, and 10.0 as k increased to 100. In every case, there was a transition from steady convection to period-2 motion near Rc2. As the Rayleigh number increased, a period-doubling route to chaos was found in all cases. Detailed scenarios of the route to chaos were similar to the case with H = 0. When k = 100, the porosity was approximately unity, and the porous layer approached a pure fluid layer. In such a situation, the effect of the interphase heat transfer, H, on the fluid dynamics was diminished. The case of k = H = 100 was calculated to simulate the situation as kH approaches infinity. Fig. 14 shows the results of the phase portraits and Poincare´ sections at R = 220, 275, 287, and 295 for periods 2, 4, and 8, and chaos, respectively. It is obvious that the solutions lie on the plane of Y = U and Z = W in all cases, which means the temperature distributions of fluid and solid are the same. Hence, the thermal equilibrium between the fluid and porous medium is achieved as kH approaches infinity. Fig. 14 provides numerical evidence, as was mentioned in Section 3, that system (22) is reduced to the system of Vadasz and Olem as kH approaches infinity.

6. Conclusions A five-dimensional autonomous dynamic system was deduced to analyze the thermal convection in a porous medium by applying the thermal non-equilibrium model. The problem was examined by means of stability of equilibria, time history, phase portraits, bifurcation diagrams, Poincare´ sections, and Lyapunov exponents. It was found that interphase heat transfer alters the routes to chaos. Also, the non-equilibrium effect tends to stabilize steady convection. With weak interphase heat transfer and a small-porosity-modified conductivity ratio, a period-doubling route to chaos was predicted. The flow undergoes a heteroclinic bifurcation when H = 5 and k = 0.1. An abrupt transition to chaos is predicted when interphase heat transfer is moderate and the porosity-modified conductivity ratio is small or moderate. As the product of k and H approaches infinity, the present system is reduced to Vadasz’s system which employs a thermal equilibrium model.

Acknowledgements The author would thank Prof. Hsien Keng Chen for valuable suggestions. This research was partially supported by the National Science Council of the Republic of China, under grant NSC94-2622-E-216-014-CC3.

References [1] Lorenz EN. Deterministic non-periodic flows. J Atmos Sci 1963;20:130–41. [2] Hilborn RC. Chaos and nonlinear dynamics. New York: Oxford University Press; 1994.

L.-J. Sheu / Chaos, Solitons and Fractals 30 (2006) 672–689 [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22]

689

Thompson JMT, Stewart HB. Nonlinear dynamics and chaos. 2nd ed. New York: J Wiley; 2002. El Naschie MS. Introduction to chaos, information and diffusion in quantum physics. Chaos, Solitons & Fractals 1996;7:7–10. Nield DA, Bejan A. Convection in porous medium. 2nd ed. New York: Springer-Verlag; 1999. Pop I, Ingham DB. Convection heat transfer: mathematical and computation modeling of viscous fluids and porous media. Oxford, UK: Pergamon; 2001. Lapwood ER. Convection of a fluid in a porous medium. Proc Cambridge Philos Soc 1988;44:508–21. Masuoka T, Rudraiah N, Siddheshwar PG. Nonlinear convection in porous medium: a review. J Porous Media 2003;6:1–32. Vadasz P. Local and global solutions for transitions to chaos and hysteresis in a porous layer heated from below. Transp Porous Media 1999;37:213–45. Vadasz P. The effect of thermal expansion on porous media convection. Part 1. Thermal expansion solution. Transp Porous Media 2001;44:421–43. Vadasz P. The effect of thermal expansion on porous media convection. Part 2. Thermal convection solution. Transp Porous Media 2001;44:445–63. Vadasz P, Olek S. Weak turbulence and chaos for low Prandtl number gravity driven convection in porous media. Transp Porous Media 1999;37:69–91. Vadasz P, Olek S. Route to chaos for moderate Prandtl number convection in a porous layer heated form below. Transp Porous Media 2000;41:211–39. Vadasz JJ, Roy-Aikens JEA, Vadasz P. Sudden or smooth transition in porous media natural convection. Int J Heat Mass Transfer 2005;48:1096–106. Straughan B. A sharp nonlinear stability threshold in rotating porous convection. Proc R Soc Lond A 2001;457:87–93. Rees DAS, Pop I. Free convective stagnation point flow in a porous medium using thermal non-equilibrium model. Int Commun Heat Mass Transfer 1999;26:945–54. Rees DAS, Pop I. Vertical free convective boundary layer flow in a porous medium using a thermal non-equilibrium model. J Porous Media 2000;3:31–44. Rees DAS, Pop I. Vertical free convective boundary layer flow in a porous medium using a thermal non-equilibrium model: elliptic effects. J Appl Math Phys 2002;53:1–12. Banu N, Rees DAS. Onset of Darcy–Benard convection using a thermal non-equilibrium model. Int J Heat Mass Transfer 2002;45:2221–8. Kuznetsov AV. Thermal non-equilibrium forced convection in porous media. In: Ingham DB, Pop I, editors. Transport phenomena in porous media. Oxford, UK: Pergamon; 1998. p. 103–30. Vadasz P. Transitions and chaos for free convection in a rotating porous layer. Int J Heat Mass Transfer 1998;41:1414–35. Wolf A, Swift JB, Swinney HL, Vastano JA. Determining Lyapunov exponents from a time series. Physica D 1985;16:285–317.