Thermal instability in a three-dimensional rigid container with prescribed heat flux at lower boundary

International Journal of Engineering Science 39 (2001) 1315±1325

www.elsevier.com/locate/ijengsci

Thermal instability in a three-dimensional rigid container with prescribed heat ¯ux at lower boundary Purna Chandra Biswal a, A. Ramachandra Rao b,* a

Department of Computer Science, Birla Institute of Technology and Science, Pilani 333031, India b Department of Mathematics, Indian Institute of Science, Bangalore 560012, India Received 7 July 2000; received in revised form 7 August 2000; accepted 8 September 2000

Abstract The Benard±Marangoni convection is studied in a three-dimensional container with thermally insulated lateral walls and prescribed heat ¯ux at lower boundary. The upper surface of the incompressible, viscous ¯uid is assumed to be ¯at with temperature dependent surface tension. A Galerkin±Tau method with odd and even trial functions satisfying all the essential boundary conditions except the natural boundary conditions at the free surface has been used to solve the problem. The critical Marangoni and Rayleigh numbers are determined for the onset of steady convection as a function of aspect ratios x0 and y0 for the cases of Benard±Marangoni, pure Marangoni and pure Benard convections. It is observed that critical parameters are decreasing with an increase in aspect ratios. The ¯ow structures corresponding to the values of the critical parameters are presented in all the cases. It is observed that the critical parameters are higher for case with heat ¯ux prescribed than those corresponding to the case with prescribed temperature. The critical Marangoni number for pure Marangoni convection is higher than critical Rayleigh number corresponding to pure Benard convection for a given aspect ratio whereas the reverse was observed for twodimensional in®nite layer. Ó 2001 Elsevier Science Ltd. All rights reserved. Keywords: Benard±Marangoni convection; Thermal instability

1. Introduction When a ¯uid layer is heated from below, it is well known that convection sets in after a critical temperature di€erence is reached between the bottom and the top of the layer. This type of convection is called thermal convection. The appearance of convection results either from

*

Corresponding author. Tel.: +91-80-309-2265; fax: +91-80-360-0146. E-mail addresses: [email protected] (P.C. Biswal), [email protected] (A. Ramachandra Rao).

0020-7225/01/$ - see front matter Ó 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 0 2 0 - 7 2 2 5 ( 0 0 ) 0 0 0 9 6 - 3

1316

P.C. Biswal, A. Ramachandra Rao / International Journal of Engineering Science 39 (2001) 1315±1325

buoyancy e€ects, Rayleigh±Benard convection or from thermocapillary e€ects, Marangoni convection. The trend setting papers by Rayleigh [1] and Pearson [2] on thermogravitational and thermocapillary instabilities have created great interest for many researchers to work in this area. In most of the works the ¯uid layer is assumed to be of in®nite horizontal extent. Benard convection in three-dimensional boxes with perfectly conductive side walls has been analysed by Davis [3], Catton [4,5] and Kirchartz and Oertel [6]. They have utilized the Galerkin method to determine the critical Rayleigh number for the onset of Benard convection. In selecting the trial function for the Galerkin method, Davis [3] and Catton [4,5] have established that any three-dimensional motion can be generated by superposing two two-dimensional motions that individually satisfy equation of continuity, if sucient number of terms are included in the expansions. Thermocapillary convection in ®nite boxes has been ®rst studied by Rosenblat et al. [7,8] and Dauby et al. [9]. These authors have presented a linear and non-linear study of thermocapillary convection in circular and rectangular containers. However, these works are based on the slippery boundary conditions at the lateral walls. This boundary condition for velocity is not realistic but enables the linear problem to be solved by using the method of separation of variables, at least for adiabatically insulated side walls. Winters et al. [10], van de Vooren and Dijkstra [11] and Dijkstra [12] have examined the pure Marangoni convection with a no-slip condition at the lateral boundaries. The velocity actually vanishes on the lateral walls and the method of separation of variables does not work. Their approaches are based on ®nite-element or ®nite-volume methods but are restricted to two-dimensional containers. For a critical review of the thermocapillary convection in ®nite size systems, one may refer Kuhlmann [13]. Recently, Dauby and Lebon [14] have investigated the thermocapillary convection in threedimensional rectangular container with realistic no-slip boundary conditions at the rigid lateral walls by prescribing temperature at the lower boundary. They assumed that the upper surface of the layer is non-deformable and that the lateral walls are adiabatically insulated. A Galerkin±Tau method is used to determine the critical Marangoni number and the convective pattern at the threshold is obtained as a function of the aspect ratios of the container. In many situation it is interesting to study the thermocapillary convection when heat ¯ux is prescribed at the lower boundary. Wilson [15] has studied the e€ect of a uniform magnetic ®eld on the onset of Marangoni convection in a layer of conducting ¯uid with prescribed heat ¯ux at the lower boundary. In this paper steady Marangoni convection is studied with or without buoyancy e€ects in a three-dimensional rectangular ®nite container, with thermally insulated rigid lateral walls, by prescribing a heat ¯ux at the lower boundary. The upper surface of the ¯uid layer is assumed to be ¯at with temperature dependent surface tension. The critical parameters for the onset of convection are obtained as a function of aspect ratios of the container by using a Galerkin±Tau method. The in¯uence of buoyancy and Biot number on the critical parameter is analysed. The ¯ow structure and the isotherms at the threshold are presented in di€erent cases. 2. Mathematical formulation The system consists of an incompressible viscous Newtonian ¯uid ®lling a three-dimensional container as shown in Fig. 1. The thickness of the container is d and the length and the width of

P.C. Biswal, A. Ramachandra Rao / International Journal of Engineering Science 39 (2001) 1315±1325

1317

Fig. 1. Geometry of the three-dimensional system.

the container are x0 d and y0 d, respectively (x0 and y0 are the aspect ratios). A Cartesian coordinate system is chosen with x-, y-axes in the plane of the lower rigid surface and the z-axis vertically upwards. The ¯uid is bounded below by z ˆ 0 and above by the free surface at z ˆ d. In undisturbed state, the lower rigid surface and the free surface are maintained at constant temperature h0 and hd , respectively by applying constant heat ¯ux b. The free surface is in contact with a passive gas at constant pressure. In the initial state, the temperature varies linearly in the z direction and the surface tension at the free surface is uniform. It is assumed that the surface tension is large enough to make the free surface ¯at when motion occurs. The surface tension r, supposed to be a linear function of temperature h, is given by r…h† ˆ r…h0 †

rh …h

h0 †;

…1†

where rh ˆ oh r is a constant which is positive for most ¯uids. Temperature condition at the free surface is given by Newton's law of cooling foz h ˆ jg …h

h1 †;

…2†

where jg is the heat transfer coecient between the free surface, f the thermal conductivity of the ¯uid and h1 is the constant temperature of the passive gas lying above the free surface. The lateral walls of the container are thermally insulated. ~ is the velocity, ~ ˆ …u; v; w† ˆ 0 and h ˆ h0 ‡ b z, where U The basic state solutions are U b ˆ …hd h0 †=d. Introducing the non-dimensional variables (d, j=d, b and lj=d 2 as the appropriate scales for the unit of length, velocity, temperature gradient and pressure, where j, l being the thermal di€usivity, the dynamic viscosity respectively) into the governing equations and then linearizing the steady disturbed state equations over the basic state solution using Boussineq approximation, q ˆ q0 …1 a…h h0 ††, where q is the density of the ¯uid, and a is the volume expansion, we get:

1318

P.C. Biswal, A. Ramachandra Rao / International Journal of Engineering Science 39 (2001) 1315±1325

~ ˆ 0; r
U

…3†

~ ˆ rP ‡ Rahez ; r2 U

…4†

r2 h ‡ w ˆ 0;

…5†

where Ra ˆ agbd 4 q0 =…lk† Rayleigh number, ez ˆ …0; 0; 1† unit vector in z direction. The corresponding non-dimensional boundary conditions are: ~ ˆ 0; U

ox h ˆ 0

at x ˆ 0; x0 ;

…6†

~ ˆ 0; U

oy h ˆ 0

at y ˆ 0; y0 ;

…7†

h‡wˆ0

at z ˆ 1;

oz u ‡ Ma ox h ˆ 0;

…8† oz v ‡ Ma oy h ˆ 0;

oz h ‡ Bi h ˆ 0

at z ˆ 1;

…9†

where ox , oy and oz denote derivatives with respect to x, y and z, Ma ˆ rh b d 2 =…lj† Marangoni number. Dauby and Lebon [14] have presented a detailed derivation of Eq. (9). On the rigid bottom boundary for a prescribed heat ¯ux, we have ~ ˆ 0; U

oz h ˆ 0

at z ˆ 0;

…10†

where Bi ˆ jg d=f Biot number. 3. Method of solution The solution of the problem given in (3)±(5), satisfying the boundary conditions (6)±(10), is obtained by Galerkin±Tau method which is a re®nement of the well-known Galerkin method [16,17]. For computational convenience, we introduce the following coordinate transformation ~x ˆ 2x0 1 x 1, y~ ˆ 2y0 1 y 1 and ~z ˆ 2z 1, which allows us to replace the real ¯uid volume by the cube ‰ 1; 1Š  ‰ 1; 1Š  ‰ 1; 1Š. As a consequence, the nabla operator becomes r ˆ …2x0 1 o~x ; 2y0 1 oy~; 2o~z †, while the boundary conditions are expressed at ~x ˆ 1; y~ ˆ 1; ~z ˆ 1. In the following analysis, we have droped tilde ` ' for convenience. The unknown quantities velocity and temperature are expanded in a series of trial functions which form a complete set satisfying the boundary conditions (6)±(8) and (10) the so called `essential' boundary conditions. The unknown ~ and h are written in the form U 

~ U h



      Ny X Nx X Nz  X ~X ~Y 0 U U jkl jkl ; ˆ Ajkl ‡ Bjkl ‡ Cjkl hjkl 0 0 jˆ1 kˆ1 lˆ1

…11†

~ X ˆ …uX ; 0; wX †, U ~ Y ˆ …0; vY ; wY † are the trial functions, Ajkl , where Nx ; Ny ; Nz are integers; U jkl jkl jkl jkl jkl jkl Bjkl and Cjkl are unknown constants. The total number of trial functions or degrees of freedom is

P.C. Biswal, A. Ramachandra Rao / International Journal of Engineering Science 39 (2001) 1315±1325

1319

thus given by 3  Nx  Ny  Nz . Eq. (11) is introduced into (4), (5) and (9) which are not satis®ed a priori by the trial functions. The unknown constants in Eq. (11) are determined by minimizing the error with usual residual method. ~ Y ˆ …0; vY ; wY † correspond to the `®nite rolls' so~ X ˆ …uX ; 0; wX † and U The expressions U jkl jkl jkl jkl jkl jkl lutions of Davis [3], i.e., the modes of convection for which one of the two horizontal components of the velocity vanish but whose non-zero components are actually functions of the three spatial coordinates. The superscripts X and Y de®ne the `x-roles' and `y-roles' of Davis which are parallel to y- and x-axes, respectively. Note that, owing to the symmetry between the `length' and `width' of the box, the solution of the eigenvalue problem may be separated into four classes which are characterized by the parity of the velocity and temperature with respect to the coordinates x and y. The four classes of solutions are denoted by EE, EO, OE and OO where E stands for even and O for odd. The ®rst and the second letters of each of these symbols refer to the parity of x- and ydependences, respectively of the temperature ®eld. The velocity and temperature distributions are chosen as 8 X < ujkl ˆ x0 fj …x† gk …y† oz hl …z†; …12† x-rolls vXjkl ˆ 0; : X wjkl ˆ ox fj …x† gk …y† hl …z†; 8 Y < ujkl ˆ 0; y-rolls vYjkl ˆ y0 gj …x† fk …y† oz hl …z†; : Y wjkl ˆ gj …x† oy fk …y† hl …z†;

…13†

temperature hjkl ˆ mj …x† mk …y† nl …z†:

…14†

The choice of de®nite form of the expressions for the functions f …
†, g…
†, h…
†, m…
† and n…
† are explained below. It is easily veri®ed that, with the choice of (12) and (13), the incompressibility condition (3) is automatically satis®ed. Since the box and the boundary conditions are symmetric with respect to the planes x ˆ 0 and y ˆ 0, it is sucient to consider horizontal trial functions f …
†, g…
† and m…
† which are either symmetric or anti-symmetric with respect to these planes, i.e., even or odd horizontal trial functions with respect to x or y. The functions w, h, ox u and oy v have the same parity in x and y as seen from (3) and (5) and the classi®cation of the trial functions may be done in terms of the parity of the mj …
† only. Thus, we choose: fj …n† ˆ …n2

1†2 Tj 1 …n†;

hl …z† ˆ …z ‡ 1†2 …z 8Rn < 0 …12 mj …n† ˆ 1 :Rn 2 …1 0

gk …n† ˆ …n2

1† Tk 1 …n†;

…15†

1† Tl 1 …z†;

nl …z† ˆ …z ‡ 1†2 Tl 1 …z†;

…16†

1† T2r …1† d1

for j ˆ 2r ‡ 1; for j ˆ 2; for j ˆ 2r …j > 2†;

…17†

1† T2r 3 …1† d1

1320

P.C. Biswal, A. Ramachandra Rao / International Journal of Engineering Science 39 (2001) 1315±1325

where n stands either for x or y and T …
† are the Chebyshev polynomials given in [18]. An explicit decomposition of the unknown pressure P is not given because that does not appear in the ®nal equations. It is easy to verify that all the essential boundary conditions are satis®ed automatically by the trial functions chosen above in (15)±(17). Field equations and natural boundary conditions are projected on the same trial functions, i.e., the governing equations and natural boundary conditions are multiplied by the trial functions before being integrated over the ¯uid volume. The purpose of this procedure is to replace the set of di€erential equations by a set of algebraic equations to determine the unknown constants in (11). By introducing the velocity and the temperature expressions given in (11) into (4), (5) and (9) and projecting on the same trial functions, and then integrating over the ¯uid volume, we get, **

X

~X
U pqr

**

j;k;l

X

uXpqNz

j;k;l

** ~Y U pqr

** vYpqNz

** hpqr

~ X ‡ Bjkl U ~Y r2 Ajkl U jkl jkl

2oz Ajkl uXjkl

! 

‡ Ra

X

!!++ ˆ 0;

Cjkl hjkl ez

…18†

j;k;l

1

‡ 2 x0 MaCjkl ox hjkl



!++ ˆ 0;

…19†

zˆ1

X






j;k;l

X j;k;l

r

2



~X Ajkl U jkl

2oz Bjkl vYjkl

‡

~Y Bjkl U jkl

! 

‡ Ra

X

!!++ Cjkl hjkl ez

ˆ 0;

…20†

j;k;l

1

‡ 2 y0 MaCjkl oy hjkl

!++ 

ˆ 0;

…21†

zˆ1

!++  X 2 X Y ˆ 0; r Cjkl hjkl ‡ Ajkl wjkl ‡ Bjkl wjkl

…22†

j;k;l

** hpqNz

X j;k;l

2oz Cjkl hjkl ‡ Bi Cjkl hjkl

!++


ˆ 0;

…23†

zˆ1

where p ˆ 1; . . . ; Nx , q ˆ 1; . . . ; Ny and r ˆ 1; . . . ; Nz 1; 


 denotes the integral over the ¯uid volume. The set of Eqs. (18)±(23) are 3  Nx  Ny  Nz in number with that many unknowns form an algebraic eigenvalue problem. Formally, this system may be written as

P.C. Biswal, A. Ramachandra Rao / International Journal of Engineering Science 39 (2001) 1315±1325

22

3

332

2

1321

3

Ajkl D1 Ra D2 44 D3 5 ‡ 4 Ma D4 554 Bjkl 5 ˆ 0; Cjkl D5 O

…24†

where D1 , D2 are matrices coming from Eqs. (18) and (20), D3 , D4 are matrices coming from Eqs. (19) and (21), and D5 , O (zero matrix) are matrices coming from (22) and (23). 4. Discussion of results System (24) is solved to determine Ra or Ma using NAG F02BJF. The solution EO is same as OE due to the symmetry of x and y and therefore it is enough to consider EE, OO and OE. The critical Mac or Rac characterizing the onset of steady convection is obtained as the minimum of the three Ma or Ra corresponding to three classes of solution EE, OO and OE. First the critical parameters Mac and Rac are evaluated with Nx  Ny  Nz ˆ 5  5  5 and 4  4  4 and found the results di€er only by 0.045%. The subsequent numerical calculations are performed only with Nx  Ny  Nz ˆ 5  5  5 which reduces considerably the number of integrations involved in calculating the various elements of the coecient matrices in Eq. (24). We ®rst consider the Benard±Marangoni convection for a special case in which Ma and Ra are equal, which implies both the buoyancy and the surface tension forces are of equal magnitude. This can be achieved for a given ¯uid by a proper choice of d as RaMa 1 ˆ q0 agrh 1 d 2 . The general case of unequal non-zero Ma and Ra is possible but it is not considered here. However, we present the results for pure Marangoni convection …Ra ˆ 0† and pure Benard convection …Ma ˆ 0†. The critical values Mac with Bi ˆ 0 for di€erent aspect ratios are presented in Table 1. Zero Biot number at the free surface implies oz h ˆ 0, i.e., the free surface is insulated. From Table 1, we observe that Mac decreases with an increase in x0 or y0 (the aspect ratios). The onset of convection occurs earlier for the case of a rectangular box, i.e., when the length and width are more than the height. This happens due to the absorption of more heat at the lower boundary. In the case of pure Marangoni convection, i.e., Ra ˆ 0 with Bi ˆ 0, the critical Mac are calculated for di€erent aspect ratios and presented in Table 2. Comparing these values with the values given by Dauby and Lebon [14] for a prescribed temperature on the lower boundary, we Table 1 Critical Mac for Benard±Marangoni convection for di€erent x0 and y0 when Bi ˆ 0 x0 =y0

1

2

3

4

5

6

7

8

1 2 3 4 5 6 7 8

203.445 166.408 149.044 135.013 129.553 124.382 116.866 113.739

166.408 104.350 94.5163 90.8192 89.2513 88.4967 88.0900 87.8498

149.044 94.5163 84.0449 82.6410 82.2020 81.6619 80.8694 80.3832

135.013 90.8192 82.6410 79.7486 79.1426 78.5484 77.5490 77.1821

129.553 89.2513 82.2020 79.1426 77.6055 77.3628 76.7888 76.5786

124.382 88.4967 81.6619 78.5484 77.3628 77.0938 76.5854 76.4655

116.866 88.0900 80.8694 77.5490 76.7888 76.5854 76.5016 76.4445

113.739 87.8498 80.3832 77.1821 76.5786 76.4655 76.4445 76.0727

1322

P.C. Biswal, A. Ramachandra Rao / International Journal of Engineering Science 39 (2001) 1315±1325

Table 2 Critical Mac as function of x0 and y0 when Bi ˆ 0 and Ra ˆ 0 x0 =y0

1

2

3

4

5

6

7

8

1 2 3 4 5 6 7 8

220.059 181.432 143.236 142.975 142.611 142.208 141.068 140.681

181.432 115.356 104.550 100.494 98.7757 97.9490 97.5036 97.2406

143.236 104.550 92.9903 91.9233 90.9827 90.4021 90.3447 90.0875

142.975 100.494 91.9233 88.3331 88.7702 86.9993 85.9109 85.6109

142.611 98.7757 90.9827 88.7702 85.9720 85.7064 85.6789 85.5095

142.208 97.9490 90.4021 86.9993 85.7064 85.4151 84.8558 84.5445

141.068 97.5036 90.3447 85.9109 85.6789 84.8558 84.6544 84.3010

140.681 97.2406 90.0875 85.6109 85.5095 84.5445 84.3010 84.2906

observe that the critical Mac is increased for all aspect ratios. Thus, the onset of convection is delayed for the case with prescribed heat ¯ux at the lower boundary. In this case we observe that the critical values are increased for all aspect ratios compared with the previous case in which Ma ˆ Ra due to the absence of buoyancy force. Hence the convection with non-zero Ma and Ra sets in earlier due to the presence of both the destabilizing buoyancy and thermocapilary forces and a similar observation was made by Nield [19]. Table 3 shows the critical Mac for pure Marangoni convection with conducting free surface (Bi ˆ 1) for di€erent aspect ratios. The critical values are found to be more than those with insulated free surface (Bi ˆ 0) for all aspect ratios. The critical values of Rayleigh number Ra for pure Benard convection with Bi ˆ 0 is computed for di€erent aspect ratios and they are presented in Table 4. The critical Rayleigh numbers are quite higher than the critical Marangoni numbers for all aspect ratios (see Tables 2±4). It appears that the onset of convection is delayed for pure Benard convection compared to pure Marangoni convection whereas the reverse was observed for two-dimensional in®nite layer [15]. A regular pattern of critical numbers decreasing with an increase in aspect ratio is observed because the system absorbs more heat at the lower boundary due to bigger length and width compared to height, in all the cases. From the critical Mac calculated earlier, the corresponding eigenvector satisfying Eq. (24) is determined. Using this, the velocity components u; v; w and the temperature h are calculated from Table 3 Critical Mac as function of x0 and y0 when Bi ˆ 1 and Ra ˆ 0 x0 =y0

1

2

3

4

5

6

7

8

1 2 3 4 5 6 7 8

287.053 236.326 224.239 212.306 210.570 209.666 208.926 205.964

236.326 174.208 151.252 149.787 148.534 147.850 146.597 145.750

224.239 151.252 131.868 130.685 130.506 129.259 129.160 128.990

212.306 149.787 130.685 128.760 128.039 126.518 125.387 125.156

210.570 148.534 130.506 128.039 125.302 124.719 124.187 124.134

209.666 147.850 129.259 126.518 124.719 123.261 123.230 122.751

208.926 146.597 129.160 125.387 124.187 123.230 122.795 122.303

205.964 145.750 128.990 125.156 124.134 122.751 122.303 122.024

P.C. Biswal, A. Ramachandra Rao / International Journal of Engineering Science 39 (2001) 1315±1325

1323

Table 4 Critical Rac as function of x0 and y0 when Bi ˆ 0 and Ma ˆ 0 x0 =y0

1

2

3

4

5

6

7

8

1 2 3 4 5 6 7 8

2089.80 1606.10 1443.44 1310.13 1232.01 1196.26 1132.33 1106.26

1606.10 1010.86 902.371 878.608 864.995 858.277 855.717 853.948

1443.44 902.371 793.336 779.191 778.474 769.741 767.663 753.245

1310.13 878.608 779.191 772.975 765.322 742.392 738.572 736.596

1232.01 864.995 778.474 765.322 738.607 735.534 733.480 732.398

1196.26 858.277 769.741 742.392 735.534 731.570 729.989 729.306

1132.33 855.717 767.663 738.572 733.408 729.989 729.337 726.736

1106.26 853.948 753.245 736.596 732.398 729.306 726.736 722.540

Eq. (11). Figs. 2 and 3 show the isovalues of velocity and isotherms at z ˆ 0:5 when Ma ˆ Ra and Bi ˆ 0 for di€erent aspect ratios. The isovalues of velocity and isotherms at z ˆ 0:5 (for the original volume) with Ra ˆ 0, Bi ˆ 0 for di€erent aspect ratios show the same behaviour as in

Fig. 2. Isovalues of velocity and isotherms at z ˆ 0:5 with x0 ˆ 4; y0 ˆ 3 when Ma ˆ Ra, Bi ˆ 0.

1324

P.C. Biswal, A. Ramachandra Rao / International Journal of Engineering Science 39 (2001) 1315±1325

Fig. 3. Isovalues of velocity and isotherms at z ˆ 0:5 with x0 ˆ 8; y0 ˆ 8 when Ma ˆ Ra, Bi ˆ 0.

Figs. 2 and 3. It is observed that, the Marangoni convection dominates the Benard convection. For a rectangular cavity (x0 6ˆ y0 ) or square cavity (x0 ˆ y0 ) the isovalues of velocity u; v are not symmetric with respect to x; y but they are symmetric about the mid planes. But isovalues of w and isotherms are symmetric with respect to x; y for square cavity (Fig. 3). Further, we observe that the number of cells increase with an increase in the aspect ratios. Acknowledgements P.C. Biswal thanks CSIR, India for ®nancial assistance and IISc, Bangalore for providing local hospitality under IISc-TIFR short term visitors program. Authors thank the referee for useful suggestions. References [1] L. Rayleigh, On the convection currents in a horizontal layer of ¯uid when the higher temperature is on the under side, Philos. Mag. 32 (1916) 529±546.

P.C. Biswal, A. Ramachandra Rao / International Journal of Engineering Science 39 (2001) 1315±1325

1325

[2] J.R.A. Pearson, On convection cells induced by surface tension, J. Fluid Mech. 4 (1958) 489±500. [3] S.H. Davis, Convection in a box: linear theory, J. Fluid Mech. 30 (1967) 465±478. [4] I. Catton, Convection in a rectangular region: the onset of motion, Trans. ASME J. Heat Transfer 92 (1970) 186± 188. [5] I. Catton, E€ect of wall conduction on the stability of a ¯uid in a rectangular region heated from below, Trans. ASME J. Heat Transfer 94 (1972) 446±452. [6] K.R. Kirchartz, H. Oertel Jr., Three-dimensional thermal cellular convection in rectangular boxes, J. Fluid Mech. 192 (1988) 249±286. [7] S. Rosenblat, S.H. Davis, G.M. Homsy, Non-linear Marangoni convection in bounded layers. I. Circular cylindrical containers, J. Fluid Mech. 120 (1982) 91±122. [8] S. Rosenblat, G.M. Homsy, S.H. Davis, Non-linear Marangoni convection in bounded layers. II. Rectangular cylindrical containers, J. Fluid Mech. 120 (1982) 123±138. [9] P.C. Dauby, G. Lebon, P. Colinet, Hexagonal Marangoni convection in a rectangular box with slippery walls, Quart. J. Mech. Appl. Math. 46 (1993) 683±705. [10] K.H. Winters, Th. Plesser, K.A. Cli€e, The onset of convection in a ®nite container due to surface tension and buoyancy, Physica D 29 (1988) 387±401. [11] A.I. van de Vooren, H.A. Dijkstra, A ®nite-element stability analysis of the Marangoni problem in a twodimensional container with rigid sidewalls, Comput. Fluids 17 (1989) 467±485. [12] H.A. Dijkstra, On the structure of cellular solutions in Rayleigh±Benard±Marangoni ¯ows in small-aspect-ratio containers, J. Fluid Mech. 243 (1992) 73±102. [13] H.C. Kuhlmann, Thermocapillary ¯ows in ®nite size systems, Math. Comput. Modelling 20 (1994) 145±173. [14] P.C. Dauby, G. Lebon, Benard±Marangoni instability in rigid rectangular containers, J. Fluid Mech. 329 (1996) 25±64. [15] S.K. Wilson, The e€ect of a uniform magnetic ®eld on the onset of steady Marangoni convection in a layer of conducting ¯uid with a prescribed heat ¯ux at its lower boundary, Phys. Fluids 6 (1994) 3591±3600. [16] C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, in: Spectral Methods in Fluid Dynamics, Springer, New York, 1988. [17] B.A. Finlayson, in: The Method of Weighted Residuals and Variational Principles, Academic Press, New York, 1972. [18] M.A. Snyder, in: Chebyshev Methods in Numerical Approximation, Prentice-Hall, Englewood Cli€s, NJ, 1966. [19] D.A. Nield, Surface tension and buoyancy e€ects in cellular convection, J. Fluid Mech. 19 (1964) 341±352.