Surface-tension and coupled buoyancy-driven instability in a horizontal liquid layer. Overstability and exchange of stability

Overstability for simultaneous surface-tension- and buoyancy-driven instability in a horizontal infinite liquid layer is theoretically investigated by means of a small disturbance analysis. Formulation and results are given in dimensionless forms. Critical wavenumbers, time constants, and Marangoni numbers are computed. Besides the influence of Prandtl, Bond, and crispation numbers, the modifications induced by interfacial viscosities, heat transfer at the free surface, buoyancy with respect to a pure Marangoni mechanism, and different thermal conditions at the rigid wall, are included in the analysis. The case of exchange of stability is considered as a special case of overstability. This work provides a generalization of Takashima's work [J. Phys. Soc. Jpn. 50, 2745, 2751 ( 1981)] concerning a pure Marangoni mechanism (with less general conditions).


I. INTRODUCTION
The present work follows a long line of studies concern ing RBM (Rayleigh-Benard-Marangoni) computations, anchored through a direct fi liation in a line of works that goes back to the original theoretical master analysis of Ray leigh. 1 Landmark papers are by Pearson, 2 who showed that Benard's cells3. 4 were not buoyancy driven but actually ten sion driven, Scriven and Sternling, 5 who addressed the ques tion of what the roles of flexibility and resistance to deforma tion of the free surface would be, and Nield/ who has the distinction of combining buoyancy and surface tension mechanisms. However, Nield assumed the same restricted free surface conditions as Pearson (nondeformable surface) and did not consider the case of overstability. The pioneering works mentioned above have been expanded by several au thors. The most general work, akin to the present one, is that of Takashima, 7 who extensively computed overstability and exchange of stability characteristics for surface-tension driven instability in a horizontal liquid layer with a deforma ble free surface, including the existence of gravity waves. A minor lack of generality included the facts that only the so called conductor case (solid wall thermal condition) was considered, and that both the Biot number (characterizing the thermal condition at the free surface) and the viscosity number (characterizing the interfacial viscosities) were zero. However, a major lack of generality arose from the fact that bulk buoyancy-driven instabilities were neglected.
More precisely, gravitational acceleration is included in Takashima's analysis but only for gravity wave effects (Bond number), not for bulk effects (Rayleigh number).
Here we are primarily interested in knowing how overstabi lity characteristics would be modifi ed if both Marangoni and Rayleigh numbers were simultaneously included in the anal ysis. Consequently, the present work can be considered as a signifi cant generalization of Takashima's analysis. It may be also considered as a generalization of Nield's work, which also considers simultaneous surface tension and buoyancy agencies, but for the simpler case of exchange of stability, while here primary emphasis is put on the case of overstabi lity.
The analysis is presented under dimensionless form. The mathematical formulation is given in a concise way but all difficulties have been indicated so that the reader wanting to check and use our derivations would be left with only straightforward algebra. Although the primary interest was in the coupling between Marangoni and Rayleigh numbers for overstability, the opportunity has been taken to include other ingredients in the analysis (Prandtl, Bond, and crispa tion numbers, modifi cations induced by interfacial viscos ities and heat transfer at the free surface, and different ther mal conditions at the rigid wall). This results in a large dimension of the parameter space, but only a selected and exemplifying set of results will be displayed.
As a by-product, this paper also considers the case of exchange of stability. The formulation to compute critical quantities for the exchange of stability is obtained as a special case of the overstability formulation. The interest is (i) to supply us with additional checking of our work by recover ing known formulas and published numerical results; (ii) provide the reader with new expressions completing, in par ticular, Takashima's work7 for exchange of stability; and (iii) provide new numerical results of interest such as the fact that interfacial viscosities may be destabilizing in some cases.
Apart from the interest of this work in fundamental un derstanding, we may mention an example of application, namely, the problem of growing large semiconductor single crystals. As exchange of stability is triggered by buoyancy, experiments have been carried out in space laboratories to obtain crystals of higher quality. Unfortunately, we show in this paper that overstability is inhibited by buoyancy. In oth er words, microgravity experiments permit the inhibition of exchange of stability but might promote the onset of oversta bility leading to oscillatory motion. Admittedly, crystal growth typically involves nonmotionless basic states (for ex-ample, in floating zones) but the fact that overstability is promoted by microgravity remains an interesting feature in the crystal growth context.

A. Geometrical setup and free surface equations
The liquid layer of infinite extent is limited by a solid wall and a free surface (Fig. 1 ) . A Cartesian coordinate sys tem (X, Y,Z) is used, where the lengths are made dimension less by using the distance d between the solid wall and the undeformed free surface as a unit of length. The axis Z is directed toward the free surface. The location Zrs of the free surface is given by where 8Zr. is the perturbation of the location of the free surface. The perturbation 8Zrs may be expressed in terms of normal modes: (2) in which/is the planform function (see Sec. II B) and satis fies a kinematic condition a8Zrs n.u, . = u, 3 = --, where Urs.i is the dimensionless velocity perturbation at the free surface and n; is the normal unit vector to the free sur face, with the unit of time t being d 2 I a where a is the ther mal ditfusivity of the liquid.
Assuming that the deformation of the free surface is infi nitesimal, that this free surface is a two-dimensional (2-D) Newtonian fluid with negligible mass, neglecting the vis cous force exerted by the gaseous medium, and following Aris, 8 Scriven and Sternling, 5 and Scriven,9 we obtain a dynamical equation for the free surface. That equation can be split into a normal force balance and a tangential force balance.
The dimensionless mean curvature of the surface, H, at a point (X, Y,Zrs ) , is given by  Phys. Fluids A, Vol. 2, No. 6, June 1990 where a 11 is the dimensionless surface Laplacian a 2 I ax�, with X 1 = X, X 2 = Y, the Einstein rule of summation is used, and a Greek subscript ranging from 1 to 2. Here, Pr, Bo(u), Cr(u), and Vi are the Prandtl, Bond, crispation, and viscosity numbers, respectively, which are given by Pr = J..L I Po a, Bo(u) =p0gd 2 lu, Cr(u) = J..L alud, where f.L is the dynamic viscosity of the fluid, p0 is the specific mass (at a reference temperature equal to the non perturbed steady-state free surface temperature), u is the surface ten sion, and ( K + c) is the sum of the surface dilatational coeffi cient and of the surface shear viscosity coefficient (interfa cial viscosities). Here, P 1 and P 2 are dimensionless pressures in media 1 and 2 , respectively ( Fig. 1 ) , with p0 a 2 I d 2 as the pressure unit, and the dimensionless vertical component of the velocity of the liquid U3 at the free surface is actually a perturbation velocity because the free surface is motionless in the base state. In relation ( 4), the term involving the Bond number is the gravity wave term. A positive g means that the upper side of the layer is a free surface.

B. Bulk equations
We assume that the liquid complies with the Oberbeck Boussinesq approximation.1 0 • 1 1 The base state is motionless with a dimensionless temperature gradient aT laXj = (0,0,-kT) in which dimensionless tempera tures are produced by using the reference temperature. Then, we follow the standard procedure for linear analy sis, [5][6][7] set that space and time are separable in the perturba tions ( p for pressure, u 3 for the Z component of velocity, r for temperature), and decompose spatial dependence in terms of normal modes: where the planform functionf(X, Y) satisfies the Helmholtz equation a uf + aY= o , 0 2 > with /3 and a being the dimensionless time constant and di mensionless wavenumber of the normal mode, respectively, and obtain the velocity, temperature, and pressure equations given below: ( fiJ 2a 2r2 )( fiJ 2a 2 )t = ( l!Pr)aTga20, ( fiJ 2a 2 q 2) ( fiJ 2a 2r 2 ) ( fiJ 2a 2 )0 =a 2 Ra 0, Using a momentum equation, the pressure equation can be modified to (16) and the temperature and velocity Z perturbations are related by (17) where the dimensionless acceleration vector gi has been tak en equal to (0,0,-g) and the operator � is d /dZ. The expansion factor aT, the Rayleigh number Ra, r2, and rf, are We then fi nd that (} is the sum of an odd and an even series, with three arbitrary constants chosen to be � 1 , C(f 4, C(f s . The results can be expressed as follows: r j r4 r s a 2 = a 4 (r2 + r2q2 + q2), a0 = a 2 (Raa 4 q2r2).

D. Free surface boundary conditions
The location of the free surface is not determined for the perturbed state. Consequently, we have four coefficients to determine: the above three for (}plus a fourth one for the free surface location. This requires four boundary conditions at the free surface. These conditions can be obtained by com bining Refs. 5 and 7. The fi rst condition expresses the kine matic condition ( 3) . The second and third conditions ex press the normal force balance ( 4) and the tangential force balance ( 5), respectively. The fourth condition is derived from the Newton heat transfer law. We obtain where Cr and Bo are taken as constants by writing u = u.
(neglecting higher-order terms), Bi is the Biot number for heat transfer at the free surface (sometimes called the Nus selt number by some authors), and Ma is the Marangoni number that may be expressed as where T. is the dimensionless free surface temperature in the basic state. In conditions ( 32 )-( 35) , it is now advantageous to replace kT8s by a single unknown 8.

E. The characteristic equations and the problem solution for overstability
We introduce the free surface boundary conditions into the three-constant solution for (}. In the case of nontrivial solutions, the characteristic determinant of the resulting ho mogeneous linear set of equations must be zero, leading to Vs 1 (x j + y j Bi) (x4 + y4 Bi) ( Xs + Ys Bi) -Bi!a2 in which j = 1 for the conductor case, ( 41) j = 0 for the insulating case.
Expressions for r0,. • • ,y5 are not given in this paper. We only mention that they appear similar to series whose con vergence is ensured by factorial terms in denominators. For the onset of overstability, {Jis i w (marginal oscillatory insta bility). For arbitrary input values, the Marangoni number Ma computed from relation (38) is not a real number. We shall call it a pseudo-Marangoni number. For a given wave number a, the critical Marangoni number Mac is obtained from the following equations: Relation ( 42), where Ma designates the pseudo-Maran goni number, determines a spectrum/3; of time constants for which the pseudo-Marangoni numbers are real. The critical time constant f3c is the value among the /3-spectrum values for which !Mal is the smallest. The critical Marangoni num ber Mac is then given by relation ( 43). Scanning over a, we search for the critical wavenumber a* for which !Mac (a) I is a minimum. The final critical time constant /3 * and Maran goni number Ma* are f3c (a* ) and Mac (a* ), respectively.

A. General formulation
This formulation is derived as a special case of the over stability formulation given in Sec. II. Time constant {3 is now set to zero, leading to r 2 = q 2 = 1 [relations ( 20) and ( 21 ) ) . Relations (39)-( 40) simplify. There is now a zero in the last columns of Nand D, and expressions for r0, • • • ,y5 must be rewritten by specifying the new relation r 2 = q 2 = 1. The relations defi ning the coefficients E!; (25)-( 31) also simplify.
Finally, ( 42) becomes an identity because the pseudo-Mar angoni number is now a real number. For a given wavenum ber a, we obtain the critical Marangoni number Mac = Ma(a), from which we determine a* and Ma* = Mac (a* ) by scanning over a.

B. Pure Marangoni effect (Ra=O)
In this case, the formulation simplifies dramatically. We fi nd that coefficients E!; can be obtained in plain explicit forms, without any recurrence, (47) 8 Cr a5S + (Bo + a 2 ) (CS 2 -2aS + a 2 C-a 3 S)

A. Generalities
To solve the set ( 42) and ( 43) for overstability, we re lied on computer programming. Thermophysical inputs of the program are Pr, Bo, Cr, Vi, Bi, and Ra. Other inputs are driving parameters used to monitor the resolution algo rithm. Outputs are critical quantities a* , w* (i.e.,/3 * /i) and Ma* by scanning over a. Typically, the CPU time needed to obtain a set (a* , w* , Ma* ) is 10 min on a IBM 3090 main frame and the number of terms in series r0, ••• ,y5 is less than about 150 for an accuracy of about 0.05% on Ma* although these fi gures obviously depend on thermophysical values.
As this paper is focused on physics, not on numerics, details of the algorithm are not provided because they would shift the purpose of our discussion. But we mention that a check of the quality of the computer program has been car-ried out by comparing our results with overstability results provided by Takashima/ leading to a perfect agreement. However, this test is clearly not fully complete since Taka shima considered the simpler case Vi = Bi = Ra = 0.
In the case of a hot rigid wall [Figs. 2 (c) and 2 (d) ] , we never observed the possibility of overstability, even when the Rayleigh (buoyancy) mechanism was introduced. How ever, this is not enough to dismiss the possibility of overstabi lity for a hot rigid wall and coupled Marangoni-Rayleigh effects, because it is impossible to scan all the ranges of ther mophysical inputs and program driving parameters numeri cally. We might simply have been unlucky. Evidence for overstability impossibility could only be obtained from a for mal proof. Nevertheless, we strongly believe in it. Because, if we accept Takashima's conclusion in the case of a pure ten sion mechanism, the occurrence of overstability, when the buoyancy mechanism is also introduced, is even more diffi- cult because buoyancy appears to be stabilizing for oversta bility (see Sec. IV B) even when it is destabilizing for the exchange of stability.
When overstability does occur [cold rigid wall, Figs. 2(a) and 2(b) ], we shall limit ourselves to present quantita tive results for the case of the influence of the Rayleigh num ber with a cold floor [ Fig. 2(b)]. Computations are usually carried out by using a continuation method, i.e., knowing the critical w and a for a value X0 of a parameter to investi-* * gate (for instance, viscosity number). Critical values for another value X 1 = X0 + t5X, with an increment t5X not too large, are computed by investigating a range [ a1 ,a 2 ] con taining a* and a range [ w1 ,w 2 ] containing w* . The use of such an "adiabatic" modification of parameter X from an initial value X0 to a fi nal value X1 to track the critical values is not fully secure because the modifi cation of a thermophy sical input can produce internal structural modifications of the involved functions that the continuation method is un able to detect (see the example in Sec. IV B). Consequently, several checks must be carried out by arbitrarily modifying the driving parameters after the determination of a fi rst set of critical quantities. For the influence of the viscosity number (ranging from 0 to 50) and of the Biot number (ranging from 0 to 3 ), the continuation method has been found to apply (quantitative results not given). Qualitatively, the ef fect of increasing Vi is to decrease the wavenumber, i.e., to increase the oscillation wavelengths. This effect can be an ticipated on physical grounds because the interface viscos ities smooth gradients on the surface. Also, an increase of Vi produces an increase of critical IMa* 1. i.e., interfacial vis cosities inhibit overstability as physically expected. Similar ly, we found that an increase of Biot number also inhibits overstability, as again physically expected, since it smooths thermal stresses in the surface. The influence of modifi cation of the thermal condition at the rigid wall will be briefly dis cussed at the end of the next section (for convenience). We now tum to the influence of the Rayleigh number for which the continuation method has been found to fail in some cases. For Ra = 0, they perfectly reproduce results given by Takashima.7 The observed discontinuities in Figs. 3 and 4 originate from the following feature. For a crispation value to the left of discontinuity, we observe two minima in curves I Mac (a) I for a= a1 and a= a2• The true critical quantities correspond, for instance, to a* = a 1 , for which I Mac (a 1 ) I < I Mac ( a2 ) I· When the crispation number in creases, I Mac (a 1 ) I increases while I Mac (a2) I decreases. At the point of discontinuity, we have I Mac (a 1 ) I = I Mac (a2) I· On the right of discontinuity, the situation is the opposite with I Mac (a 1 ) I > I Mac ( a2 ) I· Dis continuity corresponds to a jump of the critical wavenumber from a* = a1 to a* = a2 (Fig. 3 ), accompanied by a jump of the time constants fromw* =w 1 tow* =w 2 (Fig. 4).
Clearly, such a process does not produce any discontinuity in the Ma* results (Fig. 5).
For Ra = -102, discontinuity existing at Ra = 0 has drifted to the right, for both the wavenumbers (Fig. 3) and the time constants (Fig. 4 ). For the critical Marangoni numbers (Fig. 5), the difference with the Ra = 0 case is only marked on the right of discontinuity for log Cr larger than = -2.5. For Ra = -3 X 102, discontinuity has again drift ed to the right, but its importance is smaller. presents one minimum corresponding to the points labeled 1 in Figs. 3-5. For Ra = -2 X lOS, this minimum has drifted to smaller a values and to larger I Mac I values (Fig. 6). A second minimum is detected at very large values of a (a > 1. 5), drifting from the right as indicated by the arrow on the profi le labeled 2. When Ra ranges from -105 to -106, the fi rst minimum evolves along the broken line indicated in Fig. 6. Somewhere between Ra = -2 X 105 and Ra = -5 X 105 the critical values jump from the first mini mum to the second minimum. According to the continu ation method, the critical values at Ra = -106 would be found to correspond to the first minimum in Fig. 6  and w* = 700. Figure 5 shows that the Rayleigh number generally acts as an efficient inhibitor of overstability. In the situation un der study [ Fig. 2(b) ], buoyancy is also stabilizing for the exchange of stability. However, we observed in other com putations that buoyancy is stabilizing even if it is destabiliz ing for exchange of stability. These observations are prob ably connected to the fact that there is no overstability for a pure Rayleigh mechanism. 1 2 They may also be connected with the fact that, according to Davis and Homsy, 1 3 surface deflections for predominantly buoyancy-driven convection are stabilizing. We suggest that the mechanism of inhibition of overstability by buoyancy might be physically understood by remembering that, for tension-driven convection, the free surface is depressed above a hot stream while it is elevated for buoyancy-driven convection.5• 1 4 Although this last statement concerns the case of exchange of stability, it indi cates a confl ict in the direction of deformation of the free surface between the buoyancy case and the surface tension case which might be the cause of the overstability inhibition by buoyancy.
However, examination of

V. EXCHANGE OF STABILITY: NUMERICAL RESULTS
The computer program for exchange of stability is de rived from the one for overstability by internal modifica-tions. Validation of the program is carried out by comparing its results with all results given by Nield6 and Takashima, 7 leading to a perfect agreement.
We shall mention here an interesting observation. For overstability, we numerically observed instability when the rigid wall is cold, in agreement with Takashima's results for a pure Marangoni mechanism. Comparing with exchange of stability results, we notice that instability could set in as overstability only when the layer is stable for exchange of stability via a pure Marangoni mechanism under Nield's specifi cations, i.e., the conductor case at the rigid wall, with Bo = Cr = Vi = 0. In these cases, we may state that loss of stability by exchange of stability via a pure Marangoni mech anism is forbidden under Nield's specifi cations since the sys tem is stable, or at least is difficult insofar as loss of stability may arise if the specifi cations are relaxed. Then, the layer chooses another, easier way to lose stability, namely, over stability.
For a pure Marangoni effect (Ra = 0), the formulation reduces to very simple expressions [ ( 46) and ( 4 7 The region below each curve represents stable states; the critical Marangoni number Mac is always negative when Bo < 0. We note also that the negative half-plane displayed in Figs. 7 and 8 corresponds to a liquid layer suspended from a cold ceiling [ Fig. 2 (a) ] , a case which is stable under Nield's specifi cations. The loss of stability now observed is attribut ed to the introduction of gravity waves, although the free surface deformation (crispation number) also plays an es sential role. Takashima also observed that when Cr > 0.000 085, the critical a* is zero. Otherwise, critical a* becomes nonzero and Mac decreases rapidly as Cr decreases (see Figs. 7 and 8). These numerical observations are con fi rmed by Takashima's examination of relation (46). For the present case of a layer suspended from a cold ceiling (negative Marangoni numbers), stability is observed for Ma < Mac and Takashima states that cellular convection is also, in principle, possible when Ma =Mac. Conversely, for positive Marangoni numbers (layer suspended from a hot ceiling), motionless state and cellular convections cannot be observed and the liquid will, in practice, evolve to pendant drops or fall from the wall.
The influence of the viscosity number is shown in Fig. 7. At high crispation numbers (label 3, for example), interfa cial viscosities can provoke a modifi cation of the critical wavenumber a* from zero to a finite value (label 3c for instance). At small crispation numbers (labels 5-7) the modifi cation of the finite critical wavenumber a* is real but less and less signifi cant when Cr decreases. The points below each curve represent stable states, and we are led to the con clusion that the interfacial viscosities have (in this case) a destabilizing character, an unexpected result at fi rst sight. The explanation is (for instance) as follows.
If the value of the absolute temperature gradient (or of jMac j) is not large enough, the configuration point of the system lies above its neutral curve and the layer is unstable since points below curves represent stable states. Stability is obtained for a high value of the absolute temperature gradi ent. Now, the effect of interfacial viscosities is to smooth out the dynamical effects of the temperature gradients in the surface. Consequently, the absolute temperature gradient to obtain stability must be larger in the presence of interfacial viscosities in order to compensate for their damping effect. Then the domain of the unstable states is expanded as ob served in Fig. 7.
The influence of the Rayleigh number is shown in Fig. 8. A positive Rayleigh number has a destabilizing character as made obvious by Fig. 2(a). At high crispation numbers, we must increase Ra up to about 500 to observe a modifi cation in the fi gure, which would, however, remain small. When the crispation number decreases, the effect of buoyancy be comes more important. For curves labeled 6 and 7, the in duced modifications are very important even at Ra = 500. This can be physically understood if we interpret the de crease of the crispation number as an increase of layer height d and if we remember that we do expect buoyancy to become dominant when d increases.
For the conductor case, Bo = -0.1, Vi-= Bi = Ra = 0, the critical Marangoni numbers Ma* are shown in Fig. 9, and repeated in Fig. 10. Letter P indicates that the Marangoni number Mac changes its sign (becoming positive) in the a range under study. The critical wavenum bers a* are shown in Fig. 11 and are repeated in Fig. 12. These results agree perfectly well with those given by Taka shima. Again, the region below each curve represents stable states. The infl uence of the viscosity number is displayed in Figs. 9 and 11 and the infl uence of the Rayleigh number in Figs. 10 and 12. Comments are not given because these re sults merely reflect the data in Figs. 7 and 8.

VI. CONCLUSION
Overstability and exchange of stability for simultaneous surface-tension and buoyancy-driven instability in a hori zontal infi nite layer have been theoretically investigated by means of a small disturbance analysis. Formulations and ex pressions are given in dimensionless forms. Numerical re sults are obtained relying on computer programming vali dated by comparisons with the results of Takashima7 and Nield.6 For overstability numerical results, emphasis is put on the influence of the Rayleigh number, corresponding to our primary motivation. Buoyancy has a stabilizing effect, ex cept for a very restricted range of input values, even when it is destabilizing for exchange of stability.
For exchange of stability, the formulation is presented as a special case of the overstability formulation. This formu lation can be simplifi ed dramatically in the case of a pure Marangoni effect. For a conductor condition at the rigid wall, and neglecting the interfacial viscosities, our formal result is the same as the one given by Takashima for this case.
Sample cases are then given and discussed. We present the critical Marangoni numbers Ma* and critical wavenum bers a* versus the crispation number, for a given Bond num ber, with a discussion of the influence of interfacial viscos ities and of the buoyancy instability mechanism. To mention only one point in this conclusion, we observed that the inter facial viscosities may have a destabilizing effect, which, at fi rst sight, is an unexpected result.
We fi nally state that more details for the formulation and the computing algorithms, listings, and additional re sults are available on request.

ACKNOWLEDGMENTS
The numerical calculations related in this paper were performed on an IBM 3090-200 at the CIRCE. The authors wish to acknowledge the Service Commun d'Informatique de l'Universite de Rouen.