Regular and chaotic Rayleigh-Bénard convective motions in methanol and water

´ Regular and chaotic Rayleigh-Benard convective motions in methanol and water

Journal Pre-proof

´ Regular and chaotic Rayleigh-Benard convective motions in methanol and water C. Kanchana, Yongqing Su, Yi Zhao PII: DOI: Reference:

S1007-5704(19)30448-4 https://doi.org/10.1016/j.cnsns.2019.105129 CNSNS 105129

To appear in:

Communications in Nonlinear Science and Numerical Simulation

Received date: Revised date: Accepted date:

21 May 2019 12 November 2019 20 November 2019

´ Please cite this article as: C. Kanchana, Yongqing Su, Yi Zhao, Regular and chaotic Rayleigh-Benard convective motions in methanol and water, Communications in Nonlinear Science and Numerical Simulation (2019), doi: https://doi.org/10.1016/j.cnsns.2019.105129

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier B.V.

Highlights • Regular, chaotic and periodic Rayleigh-Benard convective motions in methanol and water are reported. • The comparison of heat transport between water and methanol is presented. • The merits of methanol and the demerits of water as coolants in thermal systems are documented. • From the study it is evident that methanol is stable in regular and periodic regimes and quite vigorously active in the chaotic regime compared to water. • The percentage of enhanced heat transport in water compared to methanol is just around 0.17%. • The possibility of having a transition region between regular convection and chaotic motions is demonstrated.

1

Regular and chaotic Rayleigh-B´enard convective motions in methanol and water Kanchana C.a , Yongqing Sua , Yi Zhaoa,∗ a School

of Science, Harbin Institute of Technology, Shenzhen, Nanshan District, Guangdong Province, Shenzhen-518055, China

Abstract The study of regular and chaotic Rayleigh - B´enard convective motions in methanol and water is made. The stationary mode of convection is shown to be the preferred one at the onset of convection in the case of both the liquids. Using a higher-order truncated Fourier series representation we arrive at the energy-conserving penta-modal Lorenz model and then the tri-modal Lorenz model is obtained as a limiting case of it. To keep the study analytical the Ginzburg-Landau model is derived from the penta-modal Lorenz model. It is shown that the tri- and the penta-modal Lorenz models predict exactly the same results leading to the conclusion that the tri-modal Lorenz model is a good enough truncated model for a weakly nonlinear study of convection. The Rayleigh numbers at which the onset of regular convective and chaotic motions occur are reported for both methanol and water. The behavior of the dynamical system is studied using the spectrum of Lyapunov exponents, the maximum Lyapunov exponent, the bifurcation diagram and the phase-space plots. The Hopf bifurcation Rayleigh number is obtained analytically. It is shown that the thresholds for onset of regular and chaotic motions are smaller in the case of methanol compared to water. Another very important finding of the paper is to show the existence of a developing region for chaos before becoming fully-developed. Keywords: Regular, chaos, periodic, Rayleigh-B´enard convection, Methanol, Water, Lyapunov exponent, Bifurcation diagram.

1. Introduction The convective phenomenon in fluids due to the temperature gradient between two parallel plates with the lower plate hotter than the upper one is known as RayleighB´enard convection(RBC). In the study of flow stability theory, the problem of RBC has prime importance. The RBC is considered as a standard model for studying timeindependent, spatially periodic flows and also for turbulent flow. The research community has been actively working on RBC problems from the last century. It was B´enard ∗ Corresponding

author Email address: [email protected] (Yi Zhao)

Preprint submitted to Communications in Nonlinear Science and Numerical Simulation

[1] who first noticed the hexagonal pattern on the surface of a liquid layer while doing an experiment. Later Rayleigh [2] proposed a theoretical model for it. Rayleigh [2] reported a nondimensional parameter (now referred to as the Rayleigh number) at which the instability occurs and it is defined as: Ra =

ρ0 βg∆T h3 , µχ

where g is the acceleration due to gravity, ∆T is the temperature difference between the upper and lower plates, h is the distance between the plates, β, µ, ρ0 and χ are the thermal expansion coefficient, dynamic viscosity, density and thermal diffusivity of liquid respectively. Thus, the Rayleigh number involves the liquid’s thermophysical properties. Most of the studies on RBC involve Newtonian liquid with Prandtl number, P r = 10, as a working medium ([3]- [14]). In many thermal engineering applications, for instance in heat exchangers such as refrigeration and automobile coolants, the Rayleigh-B´enard system is used as a heat removal mechanism (Carnot’s thermodynamic engine). In these applications the working liquid medium may not be a Newtonian liquid with P r = 10 ([12],[13]). The chaotic motion in a Rayleigh-B´enard convective system for a Newtonian fluid was reported by Lorenz [8]. Lorenz derived a three-dimensional, nonlinear autonomous system of ordinary differential equations which involves three parameters, namely Rayleigh number(Ra), Prandtl number(P r) and geometry ratio(b). He showed that chaotic motion manifests as we increase the Rayleigh number much beyond its critical value, i.e., Ra ≥ 24.74Rac . At this juncture, however, the pioneering works of Saltzman [9] and Lorenz [10] must also be credited. After these pioneering works ([9]-[10]) there have been many reported works on chaotic motion ([11]-[30]). It is a modern day trend to use nanoliquids as a working medium in the study of RBC ([31]-[34]). This is because the presence of nanoparticles/nanotubes in a base liquid significantly increases the stability of the liquid and enhances the heat transport. Siddheshwar et al. [32] proposed a modified two-phase model and reported an enhancement in heat transport due to nanoparticles. Kanchana and Zhao [35] showed that such an enhancement in heat transport is 24.29% for 5% of single-walled carbon nanotubes in water. Recently, Kanchana et al. [36] showed that nanoliquid is more susceptible to chaotic motion than water. Nowadays methanol is used as a heat transfer medium in many industries. This may be because methanol is an energy storage liquid [37]. Hence, in the paper we considered a modern day coolant as well as a traditional one, i.e., methanol and water, and studied onset of regular convective and chaotic motions in these two liquids. The choice of methanol was also dictated by the fact that the study is also of prime importance to the country that we live in. China makes use of methanol extensively and is one of the world’s highest consumers of the liquid [38]. In the paper the critical Rayleigh numbers at which onset of regular and chaotic motions occur are reported. Using the minimal and higher-order truncated Fourier series representations we arrived at the Lorenz and the Ginzburg-Landau models. We studied the dynamical behavior of the Lorenz system using the spectrum of Lyapunov exponents, the bifurcation diagrams and the phase-space plots. In order to quantify heat transport analytically we trans3

formed the Lorenz model into a Ginzburg-Landau model and an analytical solution of which is used to study the heat transport in terms of the Nusselt number. 2. Mathematical formulation Figure (1) is the diagrammatic representation of Rayleigh-B´enard convection in the two-dimensional plane (xz−plane). Figure (1) clearly shows that z = 0 and z = h are lower and upper bounding planes and are maintained at constant temperatures T0 + ∆T (∆T > 0) and T0 respectively. The term ∆T represents the temperature difference between these two planes. The thickness of the layer is h. Newtonian liquid (methanol/water) is confined between the planes. When ∆T reaches a critical value, (∆T )c , convection sets in. The sinusoidal curve represents the intensity pattern in a convective flow.

Figure 1: Physical configuration of the Rayleigh-B´enard convection problem.

In the Rayleigh-B´enard convection problem, the continuity, the momentum and the energy equations govern the system and are given below: ∇ · ~q = 0, ∂~q ρ0 = −∇p − ρ0 [1 − β(T − T0 )]g · kˆ + µ∇2 ~q, ∂t ∂T = χ∇2 T − (~q · ∇)T. ∂t

(1) (2) (3)

Equations (1)-(3) form a consistent system. Small-scale convective motion is assumed and hence |(~q · ∇)~q|  1. The symbols used in the governing equations (1)-(3) are ~q = uˆi + wkˆ for the velocity vector, ρ0 for the density of the liquid at T = T0 , T for the temperature, t for the time, p for the pressure, β, µ and χ for the thermal expansion coefficient, viscosity and thermal diffusivity of the liquid. Equations (1)-(3) are subject to isothermal boundary conditions: ) ~q = 0, T = T0 + 4T, at z = 0 . (4) ~q = 0, T = T0 , at z = h

4

At the quiescent basic state we have q~b = 0, p = pb (z) and T = Tb (z) and hence we get  z Tb = T0 + 4T 1 − . (5) h We now subject the basic state to a perturbation as follows: 0 ~q = q~0 (x, z, t), p = pb (z) + p0 (x, z, t), T = Tb (z) + T (x, z, t).

(6)

Experimentally, perturbation may arise naturally via thermal or mechanical background noise but in the present problem it is due to the externally controllable heating. In a two-dimensional flow, the velocity components can be expressed in terms of the stream function, ψ, as follows: u0 = −

∂ψ ∂ψ 0 ,w = . ∂z ∂x

(7)

We now introduce the following non-dimensional variables (X, Z) =

x z  , , h h

τ=

tχ , h2

Ψ=

ψ , χ

0

θ=

T . ∆T

(8)

Substituting Eqs. (6), (7) and (8) in the governing equations (1)-(3), we get: 1 ∂ ∂θ (∇2 Ψ) = Ra + ∇4 Ψ, P r ∂τ ∂X ∂θ ∂Ψ = ∇2 θ + − J(Ψ, θ), ∂τ ∂X

(9) (10)

ρ0 βg∆T h3 µCp is the Prandtl number, Ra = is the Rayleigh number k µχ and the Jacobian term, J(Ψ, θ), is defined as:

where P r =

J(Ψ, θ) =

∂Ψ ∂θ ∂Ψ ∂θ − . ∂X ∂Z ∂Z ∂X

(11)

The stress-free, isothermal boundary condition for solving Eqs. (9) and (10) are: Ψ = ∇2 Ψ = θ = 0 at

Z = 0, 1.

(12)

The vorticity and the energy equations (9) and (10) characterize the liquid flow. The solution of Ψ and θ provides information about the onset of convection and the heat transport. There are many ways to find the solution of the Eqs. (9) and (10) and the truncated Fourier series representation is one among them. In the next subsection we use the higher modes in the truncated Fourier series representation to derive the pentamodal Lorenz model and to keep the study analytical we transform the penta-modal Lorenz model into the coupled Ginzburg-Landau equation. Apart from the truncated Fourier series model of Lorenz there are other noteworthy truncated models ([39]-[41]). Two recent books that discuss about various aspects of

5

laminar and turbulent RBC are those of Goluskin [42] and Ching [43]. Beyond these main-stream methods there are two non-main-stream studies about RBC by Barna and M´aty´as[44] and Barna et al. [45]. These researchers showed innovation in solving the original hydrodynamic system through analytical means (self-similarity) and gave an explanation for B´enard cells. 2.1. Derivation of the penta-modal Lorenz model and the coupled Ginzburg-Landau equation We perform a weakly nonlinear stability analysis using the following truncated Fourier series representation respectively for the stream function and temperature: √ 2 2δ [A(τ ) sin(πκc X) + A0 (τ ) cos(πκc X)] sin(πZ), (13) Ψ = π 2 κc √ 2 θ = [B(τ ) cos(πκc X) − B 0 (τ ) sin(πκc X)] sin(πZ) πr 1 (14) − C(τ ) sin(2πZ). πr It is important to note here that the penta-model truncated representation[46] in Eqs. (13)-(14) has been chosen keeping in mind that the resulting dynamical system from them have to be energy-conserving. In ensuring that the penta-model is energy-conserving we have followed the guidelines provided by Curry [40] and Roy and Musielak [41]. Further, the bounded nature of the solution of the full system is also captured by the penta-modal truncated model. The two-terms within the square bracket in Eqs. (13) and (14) represent the twocounter-rotating B´enard cells. Substituting Eqs. (13) and (14) in Eqs. (9) and (10) and taking the orthogonality condition with the eigenfunctions associated with the considered modes, we arrive at the penta-modal Lorenz model: dA dτ1 dA0 dτ1 dB dτ1 dB 0 dτ1 dC dτ1

=

P r(B − A),

(15)

=

P r(B 0 − A0 ),

(16)

=

rA − B − AC,

(17)

=

rA0 − B 0 − A0 C,

(18)

=

AB − bC + A0 B 0 .

(19)

6

We now use the following regular perturbation expansion in the Eqs. (15)-(19)            A3 A2 A1 0 A    A03  A02  A01  A0  0               B  = 0 + ε B1  + ε2 B2  + ε3 B3  + · · · ,  0  0  0  0   B3  B2  B1  B  0     C3 C2 C1 0 C    2 4 r = r0 + ε r2 + ε r4 + . . . ,

(20)

and we assume the time variations only at the small time scale, i.e., τ1∗ = ε2 τ1 . Substituting Eqs. (20) in Eqs. (15)-(19) and comparing the like powers of ε on both the sides of the resulting equations, we get the following equations at various orders: First-order system M V1 = 0, (21) Second-order system

Third-order system

 M V2 = R21

where

 M V3 = R31

 −P r  0  M =  r0  0 0

0 −P r 0 r0 0

Pr 0 −1 0 0

0 Pr 0 −1 0

R22

R23

R24

R25

R32

R33

R34

R35

 0 0  0  and 0 −b

T r T r



 Ai  A0   i  Vi =  Bi0  , Bi  Ci

,

(22)

,

(23)

i = 1, 2, 3,

R21 = R22 = 0, R23 = A1 C1 , R24 = A01 C1 , R25 = −A1 B1 − A01 B10 ,  dA1 dA01 dB1  , R32 = , R33 = + A1 C2 + A2 C1 − r2 A1 , R31 =  ∗ ∗ ∗ dτ1 dτ1 dτ1 .  dC1 dB10  0 0 0 0 0  + A C + A C − r A , R = − 2A A − 2A A , R34 = 2 1 35 1 2 1 2 2 1 1 2 dτ1∗ dτ1∗

(24)

(25) (26)

The solution of the first- and the second-order systems subject to appropriate initial conditions, is given by  V1T r = A1

A01

r0 A1

h V2T r = 0 0 0 0

7

r0 A01

T r 0 ,

r0 (A21 +A01 2 ) b

iT r

.

(27)

(28)

Substituting Eq. (26) in the Fredholm solvability condition of the form R31 Aˆ1 + R32 Aˆ01 + R33 Bˆ1 + R34 Bˆ10 + R35 Cˆ1 = 0,

(29)

where (Aˆ1 , Aˆ01 , Bˆ1 , Bˆ10 , Cˆ1 ) = (1, P r, 1, P r, 0), we get

and

h i dA1 (τ1∗ ) Pr r0  3 ∗ ∗ ∗ 02 ∗ r A (τ A (τ = ) − (τ ) + A (τ )A ) 2 1 1 1 1 1 1 1 1 dτ1∗ 1 + r0 P r b

(30)

i h dA01 (τ1∗ ) Pr r0  0 3 ∗ = r2 A01 (τ1∗ ) − A1 (τ1 ) + A1 2 (τ1∗ )A01 (τ1∗ ) . ∗ dτ1 1 + r0 P r b

(31)

Equations (30) and (31) can be combined into a single equation given by:   r0 Pr dA˜1 r2 A˜1 − A˜1 |A˜1 |2 , = ∗ dτ1 1 + r0 P r b

(32)

where A˜1 = A1 + iA01 . In phase-amplitude form, A˜1 can be written as: A˜1 = |A˜1 |eiϕ .

(33)

Substituting Eq. (33) in Eq. (32), we get the amplitude equation:   d|A˜1 | Pr ˜1 | − r0 |A˜1 |3 . = r | A 2 dτ1∗ 1 + r0 P r b

(34)

Reverting to A1 = 1ε A, A01 = 1ε A0 , r0 = 1 and τ1∗ = ε2 τ1 , we get:   ˜ d|A| Pr 1 ˜3 ˜ = (r − 1)|A| − |A| , dτ1 1 + Pr b

(35)

where A˜ = A + iA0 . We note here that the Ginzburg-Landau equation (35) obtained using the pentamodal Lorenz model (15)-(19) is similar to the Ginzburg-Landau equation obtained for the classical Lorenz model (Ref. [47]) by considering the minimal mode truncated representation. This essentially means that the minimal and the higher modes predict the same results and hence we substitute A0 = B 0 = 0 in the penta-modal Lorenz

8

model (15)-(19) in order to arrive at the classical Lorenz model: dA dτ1 dB dτ1 dC dτ1

= P r(B − A),

(36)

= rA − B − AC,

(37)

= AB − bC.

(38)

  dA Pr 1 3 = (r − 1)A − A , dτ1 1 + Pr b

(39)

Equation (35) is now written as:

Solving Eq. (39) subject to the initial condition, A(0) = 2, we get q r 2 (r−1)P 1+P r . A(τ1 ) = r   −2(r−1)P r r (r−1)P r −2(r−1)P τ1 τ1 4P r 1+P r 1+P r + 1+P r e b(1+P r) 1 − e

(40)

In the next section we shall consider the three-mode Lorenz model (36)-(38) to discuss the chaotic motion and its characterization. 2.2. Chaotic regime, maximum Lyapunov exponent and bifurcation diagram It is theoretically known that we have the regular non-chaotic regime for values of r < rH and the chaotic regime for r > rH , where rH =

P r(P r + b + 3) Pr − b − 1

(41)

represents the value of the scaled Rayleigh number called Hopf bifurcation Rayleigh number [48] at which onset of chaos occurs and is obtained by linearizing the Lorenz system (36)-(38) around its critical point, p p ( b(r − 1), b(r − 1), (r − 1)). (42)

Multiplying Eqs. (36) and (37) by A and B, and Eq. (38) by (C − P r − r) and adding the resulting equations, we get: dA dB d dL =A +B + (C − P r − r) (C − P r − r). dτ1 dτ1 dτ1 dτ1

(43)

Integrating the above equation, we get the trapping region in the form L=

1 2 [A + B 2 + (C − P r − r)2 ], 2

9

(44)

where L is the Lyapunov function (a smooth, positive definite function that decreases along the trajectories). Eq. (44)√shows that the trapping region is an ellipsoid with center (0, 0, P r + r) and radius 2: √ A2 + B 2 + (C − P r − r)2 = ( 2)2 . (45) The trajectories traced by post-onset critical points of the Lorenz system (36)-(38) enter and stay within the sphere (45) and it deforms continuously with time hence it is identified as an ellipsoid. It changes its orientation as it evolves and we quantify such an orientation by using the Lyapunov exponents: λi = lim

τ1 →∞

Yi (τ1 ) 1 log2 , τ1 Yi (0)

(46)

where Yi (τ1 ) = [A(τ1 ), B(τ1 ), C(τ1 )], for i = 1, 2, 3, represents the length of the ellipsoid’s principal axis. The conditions λi > 0 and λi < 0 signify the rarefaction and compression of ellipsoid in its principal axis. The dynamical system (36)-(38) is a third-order one and hence we get three Lyapunov exponents, say λ1 , λ2 and λ3 . The maximum one among λ1 , λ2 and λ3 provides information on the chaotic and the periodic flow regimes if it satisfies the following condition: I. If the sum of all the Lyapunov exponents is unique and negative, i.e., if the system is dissipative and II. For λi > 0, i = 1 or 2 or 3, sensitive dependence on the choice of initial condition. We refer the maximum Lyapunov exponent as λmax in the paper. We calculate this maximum Lyapunov exponent (λmax ) and the spectrum of Lyapunov exponents (λ0i s) using the Gram-Schmidt orthonormalization method proposed by Wolf et al. [49] with a time range 0 to 1000 and time step of 0.1 for r ranging from 1 to 200 with step size 0.1. The Wolf algorithm [49] used to obtain λmax is explained in Appendix A. A bifurcation diagram for the present problem was drawn by using the local maxima of C(τ1 ) obtained numerically by solving the Lorenz model (36)-(38) using a Runge-Kutta solver with a time range 0 to 1000 and time step of 0.1 for r ranging from 1 to 200 with step size 0.1. Using the solution of the Ginzburg-Landau equation (40), in the next section we study the heat transport in terms of the Nusselt number at the lower boundary. 3. Estimation of heat transport The thermal Nusselt number, N u(τ1 ), is defined as: Heat transport by (conduction+convection) , Heat transport by conduction  R 2π 
πκc dθb ∂θ 0 dZ + ∂Z dX   , = 2π
R πκ dθb c dX 0 dZ Z=0

N u(τ1 ) =

10

(47)

where θb = (1 − Z). Substituting Eqs. (20), (21) and (22) and taking A0 = B 0 = 0, we get N u(τ1 ) = 1 +

22 C2 (τ1 ), r

(48)

r0 where C2 (τ1 ) = A(τ1 )2 . The amplitude A(τ1 ) is obtained from Eq. (40). b The steady finite amplitude convection of the Lorenz model (36)-(38) yields   1 2 1− . (49) N u(∞) = 1 + b r Using Eqs. (48) and (49) we study heat transport. The results of the paper and their discussion are presented in the next section. 4. Results and discussion The focus of the paper is on the regular convective and chaotic motions in the Rayleigh-B´enard problem by considering two different Newtonian liquids, i.e., methanol and water. The onset of convection in them and the heat transport by them are compared to determine the merits and demerits of using a modern coolant(methanol) and a traditional one(water). The Rayleigh number characterizes the onset of regular convective and chaotic motions and the Nusselt number quantifies the heat transport, and these are explained in the forthcoming subsections. 4.1. Onset of regular convective motion By considering the linear version of the Lorenz model (36)-(38) we arrived at the critical values for Rayleigh and wave numbers as: Rac =

π 4 (1 + κ2c )3 π , πκc = √ . 2 κc 2

(50)

It is well-known that the oscillatory convection is impossible for the present study and thereby the principle of exchange of stabilities is valid (Ref. [14]). Figure (2) represents the marginal stability curve for the stationary Rayleigh number. It is clear from Figure (2) that for stress-free, isothermal boundaries the rolls with the wave number, πκc = √π2 , are rendered unstable under infinitesimal perturbations when Ra takes the value 657.51 which coincides with the value in Eq. (50). This result is true for both methanol and water. The influence of the properties of methanol and water on the onset of convection can be captured using the factor, F , which is defined as: F =

µw χw ρm βm , µm χm ρw βw

11

(51)

Figure 2: Rayleigh number versus wave number.

where the subscripts m and w respectively denote methanol and water. Using F from Eq. (51) and the definition of the Rayleigh number, we may write: (52)

(∆T )cw = F (∆T )cm ,

where the subscripts cw and cm represent the critical values of water and methanol respectively. Using the thermophysical properties of methanol and water (see Table 1), we obtained F = 10.184. This indeed means (53)

(∆T )cw > (∆T )cm .

Equation (53) essentially means that the onset of regular convection begins earlier in methanol than in water. This result is due to the contribution of thermophysical properties of methanol and water. From Table 1 we notice that a) µw > µm , ρw > ρm , kw > km , (Cp )w > (Cp )m , χw > χm and b) βw < βm . Computation reveals that water has 37.08% higher value of dynamic viscosity than methanol. Thus the higher value of viscosity of water than methanol leads to the onset of regular convection earlier in the case of methanol compared to water. µ

ρ

k

(kg/ms)

(kg/m3 )

(K)

Methanol(CH3 OH)

0.00056

784.5

Water(H2 O)

0.00089

997

Newtonian liquid

β × 105

χ × 107 (m2 /s)

(1/K)

(J/kgK)

0.204

119

2540

1.02377

0.613

21

4179

1.47127

Table 1: Thermophysical properties of methanol and water.

12

Cp

4.2. Onset of chaotic motion Figure (3) depicts the AB-, AC- and BC- projections of the ABC - phase-space diagram plotted by solving the Lorenz system (36)-(38) numerically using the initial condition (1,1,1) and taking time range from 0 to 1000 with time step 0.1 for methanol and water.

Figure 3: AB-, AC- and BC-projections of the ABC-phase-space diagram of the Lorenz system (36)-(38) for the methanol and water for initial condition (1,1,1) and for r = 30.

Using in Eq. (41) the thermophysical properties mentioned in the table 1, we obtain rH = 26.7 for methanol and rH = 29.7 for water. Hence to plot Figure (3) the value 13

r = 30 is taken. The Lorenz system (36)-(38) has the following three post-onset critical points:  Λ1 = (0, 0, 0)   p p Λ2 = ( b(r − 1), b(r − 1), (r − 1)) . (54)  p p   Λ = (− b(r − 1), − b(r − 1), (r − 1)) 3

This yields the value (0, 0, 0), (8.794, 8.794, 29), (−8.794, −8.794, 29) at r = 30. Figure (3) shows that the Lorenz attractor has two loops with centers at Λ2 and Λ3 (calculated by taking r=30). The trajectory traces these two loops by approaching the critical point, Λ1 . The trajectory does not cross itself. This essentially confirms that the trajectories confine themselves to the region within the sphere defined by Eq. (45) and also proves uniqueness of the solution of the Lorenz system (36)-(38). The plots in Figures (4) and (5) are respectively those of the three Lyapunov exponents and the maximum one among these three Lyapunov exponents (λmax ) versus the scaled Rayleigh number, r, for methanol and water. It is clear from the maximum Lyapunov exponent plots that onset of chaos occurs at rH = 25.2 for methanol and at rH = 26.5 for water which are different from the values obtained analytically using Eq. (41). The explanation regarding this will be provided a little later in this section by making use of the phase-space plots. Figures (4) and (5) also depict that in the considered range of r, i.e., [0, 200], there exist a number of regimes: i. Regular(non-chaotic) convective regime: 1 < r < rH and λmax < 0 ii. Chaotic regime : r > rH and λmax  0 iii. Mildly-chaotic or nearly-periodic regime: r > rH and 0 < λmax  1 iv. Periodic regime: r > rH and λmax = 0 The attractors in these regimes can be analysed using the sign distribution of the spectrum of Lyapunov exponents (in the present problem there are three Lyapunov exponents, say λ1 , λ2 and λ3 ). From the plots pertaining to the spectrum of Lyapunov exponents, λi , i = 1, 2, 3, versus r in the Figures (4) and (5) it is evident that: ? in the regular convective regime, the attractor is a fixed point since all the three Lyapunov exponents are negative, i.e., have the sign distribution (-,-,-), ? in the chaotic regime the attractor is a strange/chaotic one since one of the Lyapunov exponents (maximum Lyapunov exponent) is positive, i.e., have the sign distribution (+,0,-), ? in the periodic regime the attractor is a limit cycle since the maximum Lyapunov exponent is zero and other two Lyapunov exponents are negative, i.e., have the sign distribution (0,-,-) and

14

? the points mentioned in the Lyapunov exponents’ plots indicate 2-torus since at these points two Lyapunov exponents vanish and one is negative, i.e., have the sign distribution (0,0,-). The values of the scaled Rayleigh number at which the above mentioned regimes exist are documented in the table 2. We need to mention here an observation from the phasespace plots that though we identified from the spectrum Lyapunov exponents that there exist 2-torus points at certain values of r (see table 2) we are unable to depict it using the phase-space plot. Working Medium

Methanol

Water

Interval for FPs

[2, 23.6]

[2, 25.8]

Interval between FPs and rH

[23.7, 25.1]

[25.9, 26.8]

rH

Mildly chaotic or nearly periodic points/intervals

25.2

[64.1, 64.3], [65.9, 66.2], 74.3, 84.7, 84.8, 87.5,94.9,96.3,

26.9

[59.2, 59.4], 60.2, 60.3, 75.7, 76.9, 86.7, [102.7, 103],

Intervals of limit cylcle [89.4,90.6], [96.8,109.6], [127.2,132.9], [133.1, 173.9], [174.1, 200] [81.3,82.6] [87.2, 91.6], [91.8,98], [109.9, 113.7], [113.9, 144.9], [145.1, 200]

2-torus points

101.9, 133, 174

91.7, 113.8, 145

Table 2: The values of the scaled Rayleigh number for FPs(fixed points), interval between FPs and rH (Hopf-bifurcation point), rH , mildly chaotic/nearly periodic points/intervals, intervals of limit cycle, 2-torus point for methanol and water obtained using the Figures (4) and (5).

The features of the Lyapunov exponents can also be studied using the bifurcation diagram as shown in Figure (6). The points mentioned in the bifurcation diagrams are the 2-torus points. They coincide with the 2-torus points obtained using the spectrum of Lyapunov exponents plots in Figures (4) and (5). Actually, excluding values of rH mentioned in Table 2 all other r−values/intervals describing the different regimes match well with the corresponding bifurcation diagrams shown in Figure (6). Thus, both the Lyapunov exponents-plots in Figures (4) and (5) and the bifurcation diagram of Figure (6) provide qualitatively similar features of the Lorenz model. In fact, the bifurcation diagram of Figure (6) also provides certain results that the λi -plots are unable to provide, for example, the exact type of periodicity can be seen in only the bifurcation diagrams and phase diagrams (see the insets in the maximum Lyapunov exponent plots and the same enlarged-phase-space plots in Figures (7) and (8)). We note at this point that the bifurcation diagram also serves the purpose of estimating rH exactly unlike the case of the λmax -plots.

15

Figure 4: Plots of the three Lyapunov exponents, λi , i = 1, 2, 3, and the maximum Lyapunov exponent, λmax , versus r for methanol (Experiment # 1) using the initial condition, (1, 1, 1), with a time range 0 to 1000, time step 0.1 and for r ranging from 1 to 200 with step size 0.1. Inset figures are plots obtained by solving the Lorenz system (36)-(38) numerically.

16

Figure 5: Plots of the three Lyapunov exponents, λi , i = 1, 2, 3, and the maximum Lyapunov exponent, λmax , versus r for water (Experiment # 1) using the initial condition, (1, 1, 1), with a time range 0 to 1000, time step 0.1 and for r ranging from 1 to 200 with step size 0.1. Inset figures are plots obtained by solving the Lorenz system (36)-(38) numerically.

17

Figure 6: Bifurcation diagrams for methanol and water. The amplitude C(τ1 ) (a local maxima of C(τ1 )) as a function of r is plotted by solving the Lorenz system (36)-(38) numerically, with a time range 0 to 1000, time step 0.1 and for r ranging from 1 to 200 with step size 0.1. The shaded regions indicate the periodic regimes.

18

(a) For r=90(Periodicity 5)

(b) For r=91 (Chaotic)

(c) For r=101 (Periodicity 3)

(d) For r=130 (Periodicity 4)

Figure 7: The phase-space plots for methanol obtained by solving the Lorenz system (36)-(38) numerically with a time range 0 to 1000, time step 0.1 and for a particular value of r.

19

(a) For r=82(Periodicity 5)

(b) For r=85 (Chaotic)

(c) For r=89 (Periodicity 3)

(d) For r=116 (Periodicity 4)

Figure 8: The phase-space plots for water obtained by solving the Lorenz system (36)-(38) numerically with a time range 0 to 1000, time step 0.1 and for a particular value of r.

20

Figure 9: Plots of the three Lyapunov exponents, λi , i = 1, 2, 3, and the maximum Lyapunov exponent, λmax , versus r for methanol (Experiment # 2) using the initial condition, Λ2 , which varies with the value of r (see Eq. (54)), with a time range 0 to 1000, time step 0.1 and for r ranging from 1 to 200 with step size 0.1. Inset figures are plots obtained by solving the Lorenz system (36)-(38) numerically.

21

Figure 10: Plots of the three Lyapunov exponents, λi , i = 1, 2, 3, and the maximum Lyapunov exponent, λmax , versus r for water (Experiment # 2) using the initial condition, Λ2 , which varies with the value of r (see Eq. (54)), with a time range 0 to 1000, time step 0.1 and for r ranging from 1 to 200 with step size 0.1. Inset figures are plots obtained by solving the Lorenz system (36)-(38) numerically.

Now coming back to the mismatch between the values of rH obtained analytically by Eq. (41) and by the λmax -plots in Figures (4) and (5), we seek recourse to phase-space plots to provide an explanation. A bit of numerical experimentation was performed to understand this mismatch.

22

Experiment # 1: A(0) = B(0) = C(0) = 1 (In fact, any (A(0), B(0), C(0)) 6= Λ2 ) (the expression of Λ2 is given in Eq. (54) and it’s value varies with r) λmax A BD Result: rH = 25.2 for methanol in place of the analytical value rH = 26.7 = rH A BD rH = 26.9 for water in place of the analytical value rH = 29.7 = rH , where the superscripts λmax , A and BD mean rH values by λmax -plots, analytical value by Eq. (41) and by bifurcation diagrams shown in (6).

Observation : The λmax -plots (see Figure (4) and (5)) match with that of the bifurcaA tion diagrams (see Figure (6)) for r > rH . Experiment # 2: (A(0), B(0), C(0)) = Λ2 , λmax BD A for methanol = rH = 26.7 = rH Result: rH λmax BD A for water = rH rH = 26.9 = rH

Observation : The λmax -plots (see Figure (4) and (5)) do not match with that of the A bifurcation diagrams (see Figure (6)) for r > 120 but match in [rH , 120] and there exist two positive Lyapunov exponents in the interval [120, 200]. From experiment# 2 and Figures (9) and (10) we come to the conclusion that the Lyapunov exponents are sensitive to the initial condition and passing critical points as an initial condition for each value of r is not a valid idea though it predicts the onset of chaos correctly. To understand the implication of these two numerical experimentations, a close look was had at the phase-space plots of Figure (11) in a) the interval [25.2, 26.7] and its neighbourhood for methanol and b) the interval [26.9, 29.7] and its neighbourhood for water. The phase-space plots in Figures (7) and (8), the bifurcation diagrams in Figure (6) and the phase-space plots in a window of the r-interval in Figure (11) explain the implication of the experimentation. It is amply clear that the onset of chaos is not all too sudden (see figure (11)). The beginning of chaos is first signaled by the trajectory starting to veer away from the fixed point to a torus at a value of rH estimated by the λmax -plots when the initial values are other than Λ2 . There is a transition region starting from the points of the trajectory veering away into the developing region to the fully-developed chaos (FDC). The transition regions for methanol and water are [25.2, 26.7] and [26.9, 29.7] respectively. Diagrammatically, one can now understand these results that are shown schematically in Figure (12).

23

Figure 11: Phase-space plots illustrating the shift from regular convection to a transition regime and shift from the transition regime to a fully developed chaotic (FDC) regime.

Figure 12: Illustration of developing and fully-developed chaotic regimes.

By a comparison between the results pertaining to methanol and water using the Lyapunov exponents plots in Figures (4) and (5) and the bifurcation diagram in Figure (6), we noticed that i) Onset of chaos is earlier in methanol compared to that in water. ii) Though onset of chaos occurs earlier in methanol than in water, methanol sustains chaotic motion for a much longer period than that in water before settling down into periodic motions and these periodic motions do not include the ones interspersed in the range where chaotic motion is present. This also means that chaotic motion is more vigorous in methanol than in water. In the next section we discuss the results pertaining to the heat transport between methanol and water using the Nusselt number. 24

4.3. Results on heat transport To keep the study analytical the Lorenz model (36)-(38) obtained from a weakly nonlinear stability analysis is transformed to the Ginzburg-Landau model using the method of multiscales (Ref. [33]). The analytical solution of the Ginzburg-Landau model (Eq. (40)) is used to estimate the heat transport in terms of the Nusselt number. Figure (13) is the plot of Nusselt number versus time, τ1 , for different values of the scaled Rayleigh number for methanol and water. Figure (13) clearly shows that as we increase the value of the scaled Rayleigh number (beyond its critical value), the Nusselt number increases.

Figure 13: Nusselt number, N u, versus τ1 for methanol and water and for different values of r.

Figure (14)a is the plot of the unsteady Nusselt number for methanol and water and for r = 5. The comparison of the heat transport between methanol and water is presented in this figure. From the figure it is clear that there is a slight enhancement in heat transfer in water compared to methanol. The percentage of enhancement of heat

25

transfer in water compared to that in methanol is given by (Ref. [50]): E=

(N uw − N um ) × 100. N uw

(55)

Figure (14)b reveals that when τ1 < 0.47, E increases with τ1 and at τ1 = 0.47, E reaches its maximum value of 0.17 then starts decreasing and at long time (τ1 → ∞) we noticed that N uw = N um . Such a result can also be seen in Figure (15) using the Nusselt number of steady, finite amplitude convection. It is clear from the figure that the Nusselt number increases with increasing r for values of r greater than the critical Rayleigh number. The variation of N u with r is same for both methanol and water.

Figure 14: N u and percentage of N u enhancement, E, for methanol and water for r = 5.

26

Figure 15: N u versus scaled Rayleigh number, r, for methanol and water.

Concluding remark In many heat transfer applications the working medium may not be water. Though water exhibits thermodynamic properties favorable to heat transfer, it can be used in only limited applications. The reason for limited applications of water is summarized in table 3 by considering the merits and the demerits of methanol and water. Hence in Liquid Methanol

Water

Merits Can be produced naturally, Stable, Low freezing point, Good extraction solvent, Used for recycling process

Demerits Reasonably expensive, Highly toxic Harmful, Less flammable

Natural resource, Cheap, Nontoxic, Physiologically harmless, Nonflammable

The freezing point is too high, The vapor pressure rapidly increases with temperature, Poor extraction solvent

Table 3: The merits and demerits of methanol and water.

those unfavorable conditions for water it is important to select an alternate heat transfer medium and this obviously depends upon the application where it is being put to use and what is the maximum operating temperature in the medium. Recent times methanol is being used in many applications, for instance methanol is used as a fuel for stunt and racing cars. The main advantages of methanol among other merits mentioned in table 3 are that it is less inflammable and can be mixed with water. In view of this the study of thermal convection using methanol as a working medium is worthwhile. The present study provides useful information that points towards a preference for methanol over water. The general conclusions from the study are:

27

i) The energy-conserving tri-modal Lorenz model and its higher-modal (the pentamodal Lorenz model) estimate heat transport identically. Hence the tri-modal Lorenz model is a good enough truncated model for a weakly nonlinear study of Rayleigh-B´enard convection. ii) (∆T )cw > (∆T )cm and hence onset of regular motion occurs earlier in methanol than water. iii) Though water has around 3 times higher thermal conductivity than methanol, the enhancement in heat transport by water compared to that by methanol is just 0.17% in any fixed time. At long times the heat transport is same for both the liquids. iv) Onset of chaotic motion sets faster in methanol than in water. v) The onset of chaos is not all too sudden as understood so far. There appears to be a transitional ‘developing-chaos region’ before sighting the ‘fully-developed chaotic region’. vi) Water is definitely much cheaper than methanol but in many applications it has many unfavorable situations whose operating conditions render water unsuitable. vii) It is worthwhile considering the RBC problem in water-methanol mixtures drawing ideas from mixture theory for thermophysical quantities. This is contemplated next for investigation by us with the feeling that such a hybrid liquid (watermethanol mixture) can be a better coolant. Appendix A. Procedure to obtain the maximum Lypunov exponent ([49]) Consider the Lorenz system (36)-(38) as:  f1 = P r(B − A) 

f2 = rA − B − AC f3 = AB − bC

 

.

(A.1)

The trajectories of the above system are defined using the variational equation: f = J
,

(A.2)

where J and
denote the Jacobian and variation of Lorenz system (A.1) and are given by:     −P r P r 0
A1
B1
C1 J = r − C −1 −A and
= 
A2
B2
C2  , (A.3) B A −b
A3
B3
C3

where
A1 is the component of the A variation of f1 . Similarly one can define all the terms in
. Thus, for the third-order Lorenz model we get 3 × 4 equations and we solve these equations using the Runge-Kutta solver with appropriate initial condition and in a time range [τini , τend ] with time step 0.1 for r ranging from 1 to 200 with step size 0.1. After finding the numerical solution of the equations we use the classical Gram-Schmidt reorthonormalization (GSR) procedure to obtain an orthonormal set that 28

serves as the new initial condition for the linearized Lorenz system. We now solve the 0 0 new set of equations using the new initial condition and a new time range [τini , τend ] 0 where τinitial = τend . For the third-order system we will have to do 3 such iterations. Declaration of interests The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgment The authors hereby acknowledge the National Natural Science Foundation of China (NSFC) for financial support under the Project number 61573119 and a Fundamental Research Project of Shenzhen under Project numbers KQJSCX20180328165509766 and JCYJ20170307151312215. The author (YS) is grateful to Harbin Institute of Technology, Shenzhen for providing her financial support. This work is part of the author YS’s Master’s project. The authors thank Professor Pradeep G. Siddheshwar, Bangalore University, India for many helpful discussions. The authors thank the referees for their excellent insights and most educative comments that greatly helped us in refining the paper. The authors shall for long remember this wonderful publishing experience. References [1] H. B´enard, Les tourbillons cellulaires dans une nappe liquide the cellular vortices in a liquid layer, Rev G´en Science Pure Application 11 (1900) 1261–1271. [2] L. Rayleigh, On convection currents in a horizontal layer of fluid, when the higher temperature is on the under side, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 32 (1916) 529–546. [3] W. V. R. Malkus, G. Veronis, Finite amplitude cellular convection, Journal of Fluid Mechanics 38 (1958) 227–260. [4] S. Chandrasekhar, Hydrodynamic and hydromagnetic stability, Clarendon Press, Oxford (1961). [5] J. Platten, J. C. Legros, Convection in liquids, Springer. Berlin (1984). [6] R. B. Bird, C. Armstrong, O. Massager, Dynamics of Polymeric Liquids, volume 1 and 2. Wiley, New York, 2 edition (1987). [7] P. G. Drazin, D. H. Reid, Hydrodynamic Stability, Cambridge University , UK (2004). [8] E. N. Lorenz, The local structure of a chaotic attractor in four dimensions, Physica D: Nonlinear Phenomena 13 (1984) 90–104.

29

[9] E. N. Saltzman, Finite Amplitude Free Convection as an Initial Value Problem-I, Journal of the Atmospheric Sciences 19 (1962) 329–341. [10] E. N. Lorenz, Deterministic nonperiodic flow, Journal of the Atmospheric Sciences 20 (1963) 134–141. [11] M. Giglio, S. Musazzi, U. Perini, Transition to chaotic behavior via a reproducible sequence of period-doubling bifurcations, Physical Review Letters 47(4) (1981) 243–246. [12] P. G. Siddheshwar, D. Radhakrishna, Linear and nonlinear electroconvection under AC electric field, Communications in Nonlinear Science and Numerical Simulation 17 (2012) 2883–2895. [13] D. Laroze, P. G. Siddheshwar, H. Pleiner, Chaotic convection in a ferrofluid, Communications in Nonlinear Science and Numerical Simulation 18 (2013) 2436–2447. [14] P. G. Siddheshwar, P. S. Titus, Nonlinear Rayleigh-B´enard convection with variable heat source, ASME Journal of Heat Transfer 135 (2013) 122502-1–12250512. [15] P. G. Siddheshwar, K. M. Lakshmi, Darcy-B´enard convection of Newtonian liquids and Newtonian nanoliquids in cylindrical enclosures and cylindrical annuli, Physics of Fluids 31 (2019) 084102. [16] J. B. McLaughlin, S. A. Orszag, Transition from periodic to chaotic thermal convection, Journal of Fluid Mechanics 122 (1982) 123–142. [17] W. Decker, W. Pesch, A. Weber, Spiral defect chaos in Rayleigh- B´enard convection, Physical Review Letters 75(3) (1994) 648–651. [18] H. Ishida, S. Kawase, H. Kimoto, The second largest Lyapunov exponent and transition to chaos of natural convection in a rectangular cavity, International Journal of Heat and Mass Transfer 49 (2006) 5035-5048. [19] J. C. Sprott, Simplifications of the Lorenz attractor, Nonlinear Dynamics Psychology and Life Sciences 13(3) (2009) 271–278. [20] R. Barrio, S. Serrano, Bounds for the chaotic region in the Lorenz model, Physica D: Nonlinear Phenomena 238 (2009) 1615–1624. [21] P. Vadasz, Analytical prediction of the transition to chaos in Lorenz equations, Applied Mathematics Letters 23 (2010) 503–507. [22] D. Puigjaner, J. Herrero, C. Sim´o , F. Giralt, From steady solutions to chaotic flows in a Rayleigh-B´enard problem at moderate Rayleigh numbers, Physica D: Nonlinear Phenomena 240 (2011) 920–934. [23] X. Wang, G. Chen, Constructing a chaotic system with any number of equilibria, Nonlinear Dynamics 71 (2013) 429–436. 30

[24] S. Paul, M. K. Varma, Bifurcation analysis of the flow patterns in two dimensional Rayleigh-B´enard convection, International Journal of Bifurcation and Chaos 22(5) (2012) 1230018-1–1230018-15. [25] C. Sparrow, The Lorenz Equations: Bifurcations, Chaos, and Strange Attractors, Springer Verlag, New York (1982). [26] A. J. Lichtenberg, M. A. Lieberman, Regular and Chaotic Dynamics, 2nd ed. Springer-Verlag, New York (1992). [27] S. H. Strogatz, Nonlinear Dynamics and Chaos, Westview Press, Massachusetts (1994). [28] R. C. Hilborn, Chaos and Nonlinear Dynamics, Oxford University Press, Canada (1999). [29] D. D. Joseph, Fluid Dynamics of Viscoelastic Liquids, Springer, New York (1990). [30] R. K. Zeytounian, Convection in Fluids, Springer, Netherlands 90, (2009). [31] P. G. Siddheshwar, N. Meenakshi, Amplitude Equation and Heat Transport for Rayleigh-B´enard Convection in Newtonian Liquids with Nanoparticles, International Journal of Applied and Computational Mathematics 31 (2015) 271-292. [32] P. G. Siddheshwar, C. Kanchana, Y. Kakimoto, A. Nakayama, Steady FiniteAmplitude Rayleigh-B´enard Convection in Nanoliquids Using a Two-Phase Model: Theoretical Answer to the Phenomenon of Enhanced Heat Transfer, ASME Journal of Heat Transfer 139 (2017) 012402-1–012402-8. [33] P. G. Siddheshwar, C. Kanchana, Unicellular unsteady Rayleigh-B´enard convection in Newtonian liquids and Newtonian nanoliquids occupying enclosures: New findings, International Journal of Mechanical Sciences 131-132 (2017) 1061– 1072. [34] P. G. Siddheshwar, C. Kanchana, Effect of trigonometric sine, square and triangular wave-type time-periodic gravity-aligned oscillations on Rayleigh-B´enard convection in Newtonian liquids and Newtonian nanoliquids, Meccanica 54 (2019) 451–469. [35] C. Kanchana, Y. Zhao, Effect of internal heat generation/absorption on RayleighB´enard convection in water well-dispersed with nanoparticles or carbon nanotubes, International Journal of Heat and Mass Transfer 127 (2018) 1031–1047. [36] C. Kanchana, Y. Zhao, P.G. Siddheshwar, Comparative study of individual influences of suspended multiwalled carbon nanotubes and alumina nanoparticles on Rayleigh-B´enard convection, Physics of Fluids 30 (2018) 084101-1–084101-14.

31

[37] B. Roland, A. Hans-Heinrich, Z. Eberhard, Operating resource, i.e., methanol, storage for use in heat transfer device of heat pump that is utilized for thermal use in e.g., stationary technology, has sorbent comprising good heat conducting connection with sheet layers DE102009015102, (2010). [38] Methanol Institute(www.methanol.org), Methanol safe handling manual 4th Edition, (2017). [39] J. H. Curry, A generalized Lorenz system, Communications in Mathematical Physics 60 (1978) 193–204. [40] J. H. Curry, Bounded solutions of finite dimensional approximations to the Boussinesq equations, SIAM Journal on Mathematical Analysis 10 (1979) 71– 77. [41] D. Roy, Z. E. Musielak, Generalized Lorenz models and their routes to chaos. II. Energy-conserving horizontal mode truncations, Chaos, Solitons and Fractals 31 (2007) 747–756. [42] D. Goluskin, Internally heated convection and Rayleigh-B´enard convection, Springer International Publishing, Switzerland (2016). [43] E. S. C. Ching, Statics and scaling in turbulent Rayleigh-B´enard convection, Springer, Singapore (2014). [44] I. F. Barna, L. M´aty´as, Analytic self-similar solutions of the OberbeckBoussinesq equations, Chaos, Solitons and Fractals 78 (2015) 249–255. [45] I. F. Barna, M. A. Poscsai, S. L¨ok¨os, L. M´aty´as, Rayleigh-B´enard convection in the generalized Oberbeck–Boussinesq system, Chaos, Solitons and Fractals 103 (2017) 336–341. [46] P. G. Siddheshwar, G. N. Sekhar, G. Jayalatha, Effect of time-periodic vertical oscillations of the Rayleigh-B´enard system on nonlinear convection in viscoelastic liquids, Journal of Non-Newtonian Fluid Mechanics 165 (2010) 1412–1418. [47] C. Kanchana, P. G. Siddheshwar, Transforming analytically intractable dynamical systems with a control parameter into a tractable Ginzburg-Landau equation: few illustrations, The Nepal Mathematical Sciences Report 35 (2019) 35–45. [48] F. Garcia-Ferrer, E. Roldan, F. Silva, G. de Valcarcel, Didactic application of numerical analysis in nonlinear dynamics:Lorenz model study, Physical Review A 45 (1992) 626–637. [49] A. Wolf, J. B. Swift, H. L. Swinney, J. A. Vastano, Determining Lyapunov exponents from a time series, Physica D: Nonlinear Phenomena 16 (1985) 285–317. [50] H. Ozoe, E. Maruo, Magnetic and graviational natural convection of melted silicon-two dimensional numerical computations for the rate of heat transfer, JSME International Journal 30 (1987) 774–784. 32