Marangoni attractors

Physica A 270 (1999) 173–181

www.elsevier.com/locate/physa

Marangoni attractors Kausik S. Das, J.K. Bhattacharjee ∗ Department of Theoretical Physics, Indian Association for the Cultivation of Science, Jadavpur, Calcutta - 700032, India

Abstract Aspects of Marangoni convection induced by heating from below are reviewed. A Lorenz like model for Marangoni convection with nite wave number in large aspect ratio situation has been c 1999 Elsevier Science B.V. All rights reserved. constructed. We have got limit cycles. PACS: 47.20.Dr; 47.20. −k Keywords: Convection; Fixed point; Limit cycle

Study of hydrodynamic instabilities [1] acquired a new momentum with the advent of dynamical systems and tremendously sophisticated experimental techniques in the 1970s. One of the oldest instabilities that has been extensively studied is the onset of convection in a thin layer of uid, heated from below. This instability can be buoyancy driven (Rayleigh–Benard) or surface tension driven (Marangoni). In a uid heated from below, the temperature falls as we go upwards, so that the density increases with height making the layer hydrodynamically unstable. However, this does not mean that the layer is dynamically unstable. This is best seen by displacing a small packet of hot, light liquid upwards, where it encounters a cooler, heavier atmosphere and consequently feels the buoyancy upthrust. The packet of displaced

uid need not keep going up, however, because the time it takes the packet to undergo the upward displacement may be larger than the time it takes the packet to loose its heat content. In the case, the instability cannot occur since the displaced packet is no longer warmer than the surrounding. For instability to occur the buoyancy eect must be more eective than the eects of viscosity and thermal conductivity. This is expressed by considering the dimensionless number R (called the Rayleigh number) de ned as (
T )gd3 ; (1) R=  ∗

Corresponding author. E-mail address: [email protected] (J.K. Bhattacharjee)

c 1999 Elsevier Science B.V. All rights reserved. 0378-4371/99/$ - see front matter PII: S 0 3 7 8 - 4 3 7 1 ( 9 9 ) 0 0 1 3 7 - 5

174

K.S. Das, J.K. Bhattacharjee / Physica A 270 (1999) 173–181

where  is the thermal expansion coecient,
T is the temperature dierence across the uid layer, g is the accleration due to gravity, d is the thickness of the uid layer,  is the thermal diusivity and  is the kinematic viscosity. The Rayleigh number expresses the competition between the buoyancy eect (numerator) and the viscous and thermal diusivity eects (denominator). It is seen that R must exceed a critical value Rc (the exact value of Rc depends on the boundary conditions at the top and bottom of the uid layer) for convection to begin. The early experiments where the existance of the convection threshold was demonstrated by Benard were however not on purely buoyancy driven convection. The uid layer was suciently thin to make surface eects important. At the surface of the uid a uctuation in temperature across the surface causes a variation in surface tension which leads to a shear, the eect of which is conveyed to the vertical motion due to the incompressibility constant. A strong surface tension gradient thus induces a strong eect on vertical motion and can cause convection in the uid layer. The analogue of the Rayleigh number is the Marangoni number [2] de ned as M=

ÿ(
T )d 

(2)

@S expresses the strength of the variation of surwhere  is the density and ÿ = − @T face tension with temperature. For a purely surface tension driven convection to occur (Marangoni convection) M must exceed a critical value. From Eqs. (1) and (2) it ÿ 12 ) which demarcates Rayleigh should be clear that there is a critical thickness dc = ( g Benrd and Marangoni eects. For d  dc , it is purely Rayleigh Benard eect, while for d dc , it is purely Marangoni. For d ' dc , the two eects are tied together. We will focus on the two extremes only. The hydrodynamics is governed by Navier–Stokes equation and the heat diusion equation. These read (in the presence of gravity)

∇P @V + (V · ∇)V = − + ∇2 V − gz ; @t 

(3)

@T + (V · ∇)T = ∇2 T (4) @t with the incompressibility constraint appearing as ∇ · V = 0. For the convection free state (conduction only state), V = 0. If the geometry is such that the uid layer exists between the horizontal plates z = 0 and z = d, with the extension in the x–y plane much greater than d (eectively in nite), then the temperature pro le in the conduction state is T = T0 (z) = T1 −

∇T z d

(5)

and the pressure is P = P0 = −0 gz ; where 0 is the constant value of the density.

(6)

K.S. Das, J.K. Bhattacharjee / Physica A 270 (1999) 173–181

175

It is convenient to work with variables de ned about the conduction state. Consequently, we take V =0+u;

(7)

P = P0 + P 0 ;

(8)

T = T0 + T 0 ;

(9)

 =  0 + 0 ;

(10)

where 0 = −0 T 0 and is small since  is small. And in the uctuation variables, we have ∇P0 ∇P 0 ∇P0 @u + (u · ∇)u = − − + 2 0 + ∇2 u − gz @t 0 0 0 =−

∇P 0 + ∇2 u + gT 0 z ; 0

(11)

@T 0 + (u · ∇)T0 + (u · ∇)T 0 = ∇2 T 0 @t or @T 0
T + (u · ∇)T 0 − u3 = ∇2 T 0 @t d

(12)

with ∇·u=

@u2 @u3 @u1 + + =0: @x @y @z

(13)

For a test of the stability of the conduction state, we need to linearize Eqs. (11) and (12) in u and T 0 and ask for the time dependence of convective roll like solutions. These solutions are of the form u = A(t) cos ax f(z) ;

(14)

T 0 = B(t) cos ax g(z) :

(15)

In the above the cos ax indicates the roll-like nature of the solutions, while f(z) and g(z) are functions that would have to satisfy the ordinary dierential equations that follow from Eqs. (11) and (12) with the correct boundary conditions. For the linearized system, the functions A(t) and B(t) must have the structure ept and the original conduction state will be unstable if Re p ¿ 0. The threshold of convection is obtained if Re p = 0. At this point if Im p = 0, we have the onset of stationary convection (exchange of stabilities) but if Im p 6= 0, we have the onset of oscillatory convection via what is known as a Hopf bifurcation. If we are discussing the Rayleigh–Benard convection, then at this point, we can take the uid layer to be con ned between two plates of large thermal conductivity, so that

176

K.S. Das, J.K. Bhattacharjee / Physica A 270 (1999) 173–181

there can be no temperature uctuation at the surfaces. It is easy to check that a couple of curl operations on Eq. (11) produces   1 @ 2 − ∇ w˜ = −R∇21 Â + [∇ × ∇ × (U · ∇)U ]Z ; (16) ∇2 2 @ 2

2

@ @ where ∇1 = @X 2 + @Y 2 , while   @ − ∇2 Â = w˜ − (U · ∇)Â : @

(17)

In the above we have made all the variables, u; T 0 ; x; y; z; and t dimensionless by appropriate scalings and the scaled variables are U ; Â; X; Y; Z; and . The components of U are u; ˜ v; ˜ w˜ and  =  is the Prandtl number. For the linear stability analysis, we simply need to drop the nonlinear terms on the right hand side of Eqs. (16) and (17), insert the solutions indicated in Eqs. (14) and (15) and use the appropriate boundary conditions on f(Z) and g(Z). The conducting plates force g(0) = g(1) = 0, while 2 f(0) = idealized stress free boundary conditions on the velocity give f(0) = f(1) = d dZ 2 d2 f(1) dZ 2

= 0. The critical value of R for the onset of convection is readily found. The diculty arises in describing the nonlinear state above the convection threshold. Analytic treatment of Eqs. (16) and (17) is virtually impossible, and numerical work is dicult. A simpli cation was oered by Lorenz, who said that the qualitative physics could be extracted by working with a few relevant modes. The choice of these modes is to be dictated by the physical constraint. In this case it was important to (i) have modes that exhibit the roll like solution for the velocity as well as the temperature eld and (ii) have modes that exhibit the convective transfer of heat from the bottom plate to the top plate. A minimum model (for idealized boundary conditions) w˜ = A() cos aX sin Z ; Â = A() cos aX sin Z + C() sin 2Z ;  u˜ = − A() sin aX cos Z ; a v˜ = 0 :

(18)

These modes satisfy the equations of continuity and while A() and B() ensure the roll structure for velocity and temperature elds, the mode C() represents the convective heat ow from the bottom to the top plate. What Lorenz [3] did was to introduce the above expression for w; ˜ u; ˜ v˜ and  in Eqs. (16) and(17) and matching the coecients of the dierent Fourier terms, obtained A˙ = (−A + B) ; B˙ = −AC + rA − B ; C˙ = AB − bC ;

(19)

K.S. Das, J.K. Bhattacharjee / Physica A 270 (1999) 173–181

177

where b = 83 for the chosen geometry and boundary conditions and r = RRc ; Rc being the threshold for convection. The above equation constitutes a dynamical system for the Rayleigh–Benard convection. The trivial xed point A = B = C = 0 corresponds to the conduction state since it implies no rolls and no convective heat transfer. The other xed points C = r − 1;

A2 = B2 = b(r − 1)

(20)

correspond to a real convective state for r ¿ 1. Roll amplitudes exist and heat is transferred by convection. The two xed points exchange stability at r = 1, when the trivial xed point loses stability. The instability at r = 1 is stationary. It is interesting to note that the full hydrodynamic equations also yield a stationary instability at r = 1, i.e., R = Rc . The important question is how does the nonlinear state behave as r is increased beyond 1. It is seen from Eq. (20) that the Nusselt number (i.e., the amount of convective heat transfer) represented by the mode C increases as r − 1, while the aplitude 1

of the convection rolls increase as (r − 1) 2 . We now enquire if the state will lose instability as r is increased. It is seen that a Hopf bifurcation occurs at r = r0 = 

+b+3 −b−1

(21)

at which point an oscillatory state (limit cycle) should appear. This oscillatory state is not seen since the bifurcation is backward. The usefulness of the Lorenz model lies in the existence of a strange attractor, which exhibits sensitivity to initial conditions. Hence the benchmark nature of the Lorenz paper. With the ourishing of interest in dynamical system, Lorenz like models were used very eectively to probe the properties of the nonlinear state in somewhat more complicated systems [4 –8]. These studies were quite eective in bringing theory and experiment [11] together. This motivated us to look at a Lorenz like model for Marangoni convection. Unlike the RB convection, Marangoni convection has not been extensively investigated. Here we focus on pure Marangoni convection where we do not yet have an exact answer to fairly simple questions, e.g. the nature of onset of convection – the convection is in general stationary but a principle of exchange of stabilities has not yet been established [9,10]. In the Marangoni limit, we have a uid layer with a free surface and in addition to the velocity uctuation U and the temperature uctuation Â, we need to consider the surface de ection Á(X; Y; ) - the uctuation of the surface about its mean position. The vertical velocity w˜ at the point z = 1, is given by dÁ d w˜ =

dÁ d

at z = 1 :

(22)

In the dynamics of w, ˜ buoyancy no longer aects and hence R = 0 in Eq. (16), which consequently becomes   1 @ − ∇2 w˜ = [∇ × ∇ × (U · ∇)U ]Z (23) ∇2  @

178

K.S. Das, J.K. Bhattacharjee / Physica A 270 (1999) 173–181

while the temperature diusion equation remains the same as before. The roll-like solutions now have the form U = A() cos aX f1 (Z) ; (24)

 = B() cos aX g1 (Z) ; Á = D() cos aX : For the linear stability analysis,   1 @ 2 2 2 2 − (D − a ) f1 (Z) = 0 ; (D − a )  @ 

(25)

 @ − (D2 − a2 ) g1 (Z) = f1 (Z) ; @

(26)

d . where D = dZ It is the boundary conditions which are crucial now. At the lower plate (conducting thermally)

f1 = Df1 = g1 = 0

at Z = 0 :

(27)

At the upper plate, i.e. Z = 1 f1 =

dÁ ; d

Dg1 = 0

(28)

  2    @2 @ @2 Â @2 − w ˜ − M B− @X 2 @Z 2 @X 2 @X 2   2 @2 @ @w˜ @ on Z = 1 ; + 3 − =M · Cr @Z 2 @X 2 @ @Z @Â =0 @Z

on Z = 1 :

(29) (30)

2

ÿd 0 is the Marangoni number (= S10 @S In the above, M = S0 @T , S0 being the mean surface 2

is the Crispation number and B = dS0 g is the Bond number,  being tension), Cr = S 0d the density of the uid. Now comes the crucial issue of the choice of modes. We take two modes for the velocity w, which we write as w˜ = [A()g(Z) + D()f(Z)] cos aX :

(31)

The mode A() corresponds to no surface uctuation and the mode D() contains all the information on the surface uctuations. Thus, we require g(1) = 0, while f(1) = 1. The function g(Z) is chosen such that w˜ = A()g(Z) cos a X sati es ∇4 w˜ = 0 when A() is a constant. This yields g(Z) = sinh aZ − aZ cosh aZ +

aC − S Z sinh aZ ; S

(32)

K.S. Das, J.K. Bhattacharjee / Physica A 270 (1999) 173–181

where C = cosh a and S = sinh a. The function f(Z) needs to satisfy f(0) = and accordingly, we choose f(Z) = Z 2 :

179 d(0) dZ

=0 (33)

Keeping in mind the boundary conditions on  and the lesson of the Lorenz model for Rayleigh convection, we take  = B() sin

3Z Z cos aX + C() sin : 2 2

(34)

To determine the dynamics of the modes, we insert w˜ and  from Eqs. (32) and (35) into Eq. (17) and match the corresponding Fourier coecients. This yields equation for B and C. For the modes A and D, we need to satisfy Eq. (16) in the mean and satisfy the boundary condition given in Eq. (30). There is an unknown variable ‘a’, the wavenumber. This is chosen, as in the Lorenz model, to be the critical wavenumber at the onset of stationary convection. This gives a ' 2. In the dynamics of A; D; B; C all the coecients are known. Approximating the coecients, so that they have a ‘clean’ look, we obtain the dynamical system 3A˙ −

D˙ = 9D ; 2

9A˙ − D˙ = − B˙ =



   3 4 13 − 9 A + 20 + D+ B; M Cr 2M Cr Cr

A 3D 13B + − − AC − 3DC ; 4 4 2

C˙ = −

45 C + AB + 3DB : 2

(35)

(36)

(37) (38)

The trivial xed point (A=B=C =D=0) corresponds to the conduction state. Exchange of stabilities occurs at M = M0 =

169 : 2 + 117Cr

(39)

This is the only possible bifurcation – although we have a cubic equation for the growth rate – a Hopf bifurcation does not occur. The critical Marangoni number for Cr = 0 is M0 = 169 2 = 84:5. The exact answer in the limit is M0 = 81, which goes to show how good the truncation is for the study of the initial bifurcations at least. We can also make the assertion within the model that the rst bifurcation is always an exchange of stabilities. Oscillatory instability does not occur at the onset of convection. Within the complete set of Navier–Stokes’ equation and the heat diusion equation, it has not been possible to establish the principle of exchange of stabilities. However, in all experiments with large aspect ratio, the onset of convection is always stationary. The oscillatory convection occurs only for small aspect ratios.

180

K.S. Das, J.K. Bhattacharjee / Physica A 270 (1999) 173–181

Fig. 1. The simple limit cycle at M = 204.

We now explore the steady convection state. This is given by the non-trivial xed point of Eqs. (36) – (39). A20 =

1 1 169 ( M0 − M ) ; 2 13 M − 9Cr

Bo2 =

169 32

C0 =

169 8





(40)

13 − 9Cr M 1 1 − M0 M



1 1 − M0 M

 ;

(41)

 :

(42)

The structure of the xed point immediately tells us that the model is valid only for M¡

13 1 · : 9 Cr

(43)

Typical values of Cr being 10−3 , this makes the range of validity cover a large range of M . The nontrivial xed point of Eqs. (41) – (43) exists for M ¿ M0 and it is easy to see that the bifurcation is forward. Numerical integration of Eqs. (36) – (39) supports our contention regarding the rst bifurcation. For Cr =10−4 , the trajectory settles down to the trivial xed point for M ¡ 85 and to the nontrivial xed point for M ¿ 85. For Mc ≈ 200, Hopf bifurcation occurs and the resulting limit cycle is shown in Fig. 1. It is in these limit cycles that the role of the surface uctuations become apparent. In the absence of the mode D, there are no stable limit cycles in the system. It is the coupling of the surface uctuations to the heat diusion that gives rise to the stable limit cycle. It is our contention that experiments on Marangoni convection when carried on beyond the initial onset will show onset of oscillations in large aspect ratio systems.

K.S. Das, J.K. Bhattacharjee / Physica A 270 (1999) 173–181

181

Acknowledgements One of the authors (KSD) would like to thank CSIR India for partial nancial assistance and Dr. Anita Gangopadhyay for her cooperation during numerical computation. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11]

S. Chandrasekher, Hydrodynamic and Hydromagnetic Stability, Oxford University Press, London, 1961. D.A. Nield, J. Fluid Mech. 19 (1964) 1941. E.N. Lorenz, J. Atmos. Sci. 20 (1963) 130. G. Ahlers, M. Luecke, Phys. Rev. A 35 (1987) 470. S.J. Linz, M. Luecke, Phys. Rev. A 35 (1987) 3997. E. Knobloch, Phys. Rev. A 34 (1986) 1538. F.H. Busse, K.E. Heikes, Science 208 (1980) 173. J.K. Bhattacharjee, A.J. McKane, J. Phys. A 21L (1988) 555. J.K. Bhattacharjee, Phys. Fluids A 1 (1983) 1938. M. Takashima, J. Phys. Soc. Japan 28 (1970) 810. M.F. Schatz, S.J. VanHook, W.D. McCormick, J.B. Swift, H.L. Swinney, Phys. Rev. Lett. 75 (1995) 1938.