Pergamon

Adv. Space Rex Vol. 22, No. 8, pp. 1159-l 168, 1998 Q 1998 COSPAR. Published by Elsevier Science Ltd. All rights tr.sctvd Printed in Great Britain 0273-l 177/98 $19.00 + 0.00 PII: SO2734 177(98)00159-8

ON THE BOUSSINESQ APPROXIMATION WITH DEFORMABLE INTERFACES

FOR FLUID SYSTEMS

D.V. LYUBIMOV , T.P. LYUBIMOVA*, J.Iwan D. ALEXANDER** and N.I. LOBOV Perm State University, 15, B&rev Street, 614600 Perm, Russia *Institute of Continuous Media Mechanics UB RAS, 614061 Perm, Russia **Center for Microgravity and Materials Research, University of Alabama in Huntsville, AL 35899, USA.

ABSTRACT The assumptions underlying the Boussinesq approximation place restrictions on its applications to systems with free surfaces and interfaces. In this paper we reconsider the limits in which the Bouesinesq approximation is valid and develop a generalization of the approximation which allows for self-consistent application to such systems. The Boussinesq limit is characterized two parameters, G = gL3/v2, and 6 = /30, where G is a dimensionless measure of gravitational acceleration or system size and c is the product of the fluids coefficient of thermal expansion with a characteristic temperature difference. The Boussinesq limit, in geno;t and E + 0 while the product GE (equal to the familiar Grashof number) eral, corresponds to G remains finite. We consider two problems involving deformable boundaries: the stability of a twolayer fluid system heated from above and below and the influence of buoyancy on long-wavelength Marangoni instability. In the first, two examples, we examine the conditions required to consistently account for the effects of a deformable surface on thermal convection while simultaneously applying the Bousinesq approximation. In particular, the effect of the deformable surface can be included through a term proportional to the product of G6 with the deflection, C, of the interface from planarity, but it is required to treat the density of the two 0). In the second example, that of lo~-~v~~~h fluids as equal in the ~uations of motion (i.e., 6 Marangoni instability, we find that to ~rn~t~~usly consider the effects of thermal buoyancy together with a deformable surface we must either treat the surface as undeformable (G + 00 and 6 -+ 0) or consider finite G and take GE and e to be independent parameters in the equations of motion. This leads to a result which explicitly reveals the role of buoyancy as a destabilizing influence on long-wavelength instability.
01998 COSPAR. Published by Elsevier Science Ltd.

1. INTRODUCTION When studying thermal convection in systems with free surfaces or fluid interfaces, surface deformations must be properly accounted for. In particular, for fluids which can be described using the Navier-StokesBoussinesq (NSB) equations, the presence of a free surface requires that care be taken in the application of the boundary conditions. It is known that for small temperature related density nonuniformities, the deviation of the surface shape from planarity is proportional to e. Within the framework of the Boussinesq approximation, interface deformation is, thus, generally neglected. In this paper the application of the Boussinesq approbation is considered in the context of a ~er~ti~y heated two-layer system with a deformable interface, an encapsulated non-isothermal drop and the interaction of bouyancy with longwavelength Marangoni instability. It is found that the assumptions underlying the approximation place strong restrictions on the way in which such a problem can be approached. The Boussinesq approximation is based on two assumptions (Boussinesq, 1903; Chandrasekhar, 1961). Firstly, that isothermal compressibility can be neglected, and secondly that thermal expansion (manifested

1160

D. V. Lyubimov er al.

through temperature dependence of density) need only be considered in the buoyancy force term in the momentum equations. The first assumption requires that the magnitude of the fluid velocity be small in comp~ison to the speed of sound and that the size of the 3uid system under consideration be less than the adiabatic length. We will suppose that the conditions for the first assumption are satisfied and consider the restrictions imposed by the second condition. The density p = p(T) is taken to be a linear function of temperature, T, and L, L2/u, u/L, r)w/L2, 6 and pe are used to scale for the length, time (t), velocity (q, pressure (p), temperature and density (p) respectively. Here L, v, q, B and po are, respectively, the characterstic length, kinematic viscosity, dynamic viscosity, characteristic temperature difference and reference density. The resulting dimensionless equations representing conservation of momentum, energy and mass are

p(g+(w)3) =-Vp
i3T dt f vs VT = PYAT,

+ Aa -I- ;Vdiu$ + GET?, p = IcT (2)

divpv= rPr-AT,

Note that, it is implicit that we have taken the bulk viscosity to be zero. At the deformable free surface z = {(z, y, t) the force balance conditions yield -p+G<+r,,=CaK, mg=Cl (3)

for the normal, rnn and tangential, rnt, forces. The kinematic boundary condition is

Here j! is the unit vector directed vertically upwards, PP= y/x, Ca = uL/~v, x = ~~~C~ is the thermal diffusivity, K is the thermal conductivity and C, is the heat capacity. To obtain the Boussinesq equations Gem (l)-(2) one should take E equal to zero everywhere except in the last term in (1). Formally, this means that one considers the limit E + 0,G ----t 00, while the product GE and all other qu~tities in (l)-(2) remain finite. For the boundary conditions, it then follows from (4), that the normal velocity vanishes on a flat (C = 0) free surface which is parallel to the zr- y plane. This, together with (3) leads to the usual conditions for an impermeable, tangential stress-free nondeformable boundary. In practice, the magnitudes of G and the capillary number Ca are often comparable. In such cases it is reasonable to take Ca 00 simultaneously with G Y 00, such that the Bond number Bo = G/Ca is finite. The surface shape is then determined by the combined, competitive effects of gravity and surface tension. However, the ~p~cation is unchanged; that is, under the Boussinesq app~~mation, convection should be studied for given free surface shape which is determined by non-convective factors. As emphasized earlier, because of the nature of the limit in which the Boussinesq approximation is valid, care must be taken when using it and simultaneously accounting for the interface deformation. Usually, one accounts for surface d~ormations by expanding the boundary conditions about some flat (or other equilibrium shape) in some small parameter. This parameter is usually proportions to l/G (Davis, 1980; Renardy, 1986; Pavithran, 1991). This means that we must retain terms of (at least) first order with respect to l/G since the surface is flat at the leading order, In this case, consistency requires that we keep terms of order E in the NSB equations because they are of the same order as l/G. To retain terms of order 6 in the boundary conditions while neglecting them in the NSB equations would be inconsistent and can lead to incorrect results. For example, the long wavelength ~stab~ty of a ~ffer~ti~y heated free surface layer (Izakson, 1968) was found to be absent when surface deformation of the layer wag accounted for (Nepomnyaschii, 1983). Order c terms can be readily incorporated into the governing equations (see the pioneering paper of Oberbeck (1879). However, this only leads to minor corrections to solutions obtained using the conventional approach. Some problems do not physically admit such an approximation. For example when the density difference (Ap) across a surface or interface is of the same order of magnitude as therm~y induced density variations, gravity is generally insufficient to keep the surface flat. In this case deformations of the surface are

Boussinesq Approximation

for Fluid Systems

1161

not necessarily small and they must be properly accounted for. In order to be consistent, the Boussinesq approximation can be applied consistently for cases when both Ap and 6 are small by retaining Ap in the force balance boundary condition at the surface and taking p = 1 in equations (l)-(2) in the interior of the fluids. To clarify this point we consider the case where the quantity S = (pr - pz)/pr tends to zero so that the dimensionless density of the two fluids approaches 1. However, the force balance boundary condition contains a term G6 which, like GE in the 6 4 0 limit, can remain finite. In this way the surface deformation can be accounted for within the Boussinesq approximation. In section 2 we develop a generalization of the Boussinesq appro~mation. In sections 3 we illustrate the importance of the r~t~ctions imposed on the application of the Boussinesq approximation for systems with deformable interfaces by employing the equations obtained in this limit to examine the influence of interface deformation on the stability of a two-layer system subject to a vertical temperature gradient. In section 4 we consider a different situation for which the density difference across the interface between two immiscible fluids is large, but interfacial deformations are unrelated to buoyancy-driven convection. This can occur, for example, in the case of propagating capillary waves. In this case the influence of non-isotherms&y cannot be consistently accounted for without a suitable generalization of the Boussinesq approximation, Such a generalization, and its consequences are applied to the problem of coupled Marangoni-Rayleigh instability in the long wavelength limit. 2. GENERALIZED BOUSSINESQ APPROXIMATION

Consider the application of the Boussinesq approximation to case of thermal buoyancy-driven convection in a two-layer system of immiscible fluids with a deformable interface. We assume that the fluids densities depend only on temperature and that the coefficients of the dynamic viscosity, p, and thermal conductivity, K, are constant. According to the general idew of the Boussinesq approach we assume that the density non-uniformities are small besides it is true both for temperature induced non-uniformities of density and for the difference of densities of two fluids. As mentioned in the introduction only in this case the deformability of the interface is consistent with the Boussinesq equations. We represent the state equation in the form
Pj = PjO[l -

PjfTj -

TO)]

(5)

where the subscript j referring to fluid 1 or 2, (1 is the lower fluid) and pje is the density at 2 = TO. Assume that the parameters @#? and 6 = (poi - pez)/(poi + ~02) are small. Applying the conventional procedure of limit transition &3 + 0, S 0, G 00 we come to the following system of equations and boundary conditions for the description of convection in two-layer system of fluids with close densities and deformable interface

2 +(GjV)Cj
3 [q = 0, ITI = 0,

= -kVpj j

+ vjA$j + gpj(Tj - To)?, d~v~j = 0,

+ ($jV)Tj = KjATj,

-Ip]ik + [?I *n+ [po]gCib= aKn [KVT] 0n = o, g + (vV)(; =: $7.

Here [f] denotes the jump of a variable f on the interface, i.e. [f = A - fz; 7 is the Newtonian viscous stress tensor, ii is the unit vector normal to the interface, QIis the surface tension coefficient, and K is the sum of the principal curvatures of the interface. 3. STABILITY OF TWO HORIZONTAL LAYERS OF IMMISCIBLE FLUIDS

Normal mode perturbations

1162

D. V. Lyubimovet al.

Consider a two-layer system of immiscible fluids confined between .two horizontal rigid surfaces z = fh. Each fluid occupies equal part of the total volume. The boundaries are held at fixed different temperatures 01 and 02. For this situation, the equations formulated in section 2 admit a stationary solution that corresponds to zero motion in each fluid (y7i = 0 ) and a flat horizontal interface C = 0. The (conductive) temperature distributions in the layers are
TOI = Or +

AIZ,

To2 = @I+ AO.z, A1 =

m(@z - 01)
&l t 62)

t-42=

&(@a - Or) h(4 +&a)

Here Or is temperature

at the interface in unperturbed

state.

We now consider the stability of the above basic state to small perturbations of the form T = Toj f Tj, F = 6 + i&p = poj + pj,j = 1,2, and consider perturbations of the interface such that it is located at z = C. Substitution of these perturbations into the governing equations (6)-(7) yields the following boundary value problem for the velocity, temperature and pressure perturbations in each fluid: 1 aa. &1= Pr at -Vpj + vj Avfi + RapjTjq, HiVZj = 0, T=O, (11) (12) (13) 04)

&Fj at = XjATj f Aji!jT, at at t= 0:

z = fl : vr=&,
= fi27p
BT2

v=O,

-AI~tTl=-A~~+T2,
K av1, = V2-37
avat

dZ

(15)
@C Z=v.T.

-(w-pa)-GuCi-2 ~$5~~3

=CaA2C,

fW

Here As is the Laplace operator in the variables x,y, and ujt denotes the component of velocity of fluid (j =1,2) tangent to the deformed interface. Equations (ll)-(16) are written in non-dimensional form with h, h2,x*, x,/k x*~*~/~~ and 0 = (@I - f32>/2 as scale factors for length, time, velocity, pressure and temperature, respectively. The scale factors for kinematic viscosities v*, thermal diffusivities x*, thermal conductivities K*, thermal expansion /3*, and for the equilibrium temperature gradients A, are taken to be the mean values for the two fluids. The boundary value problem (ll)-( 16) includes the following non-dimensional parameters: Prandtl number Pr = u*/x*, Rayleigh number Ra = g@,0h3/u,x*, capillary number Ca = ah/u,;y,po, and Gal&i number Go = (ps - ~~)gh3/~*x*. The parameter Ga is determined so that the positive values of Ga correspond to a more dense upper fluid, a potentially unstable stratification. For this choice of units the following relations hold
Y t v2 =

2,

PI t /32 = 2,

XI + x2 = 2,

IEI t

~2 =

2,

AI =

~2,

A2 =

~1,

(171

Thus, the problem is characterized by seven independent nondimensional parameters. For the case of tied physical properties of the fiuids (visc~ity, ~nductivity, etc.) only Ra, Ga, Ca axe independent . We consider normal mode perturbations proportional to exp(& + ikz), where A is the growth rate and k is the wave number. For the perturbation amplitudes u , w , 8 andq of the x-component of velocity, the z-~mponent of velocity, the temperature and the pressure respectively we obtain the spectral problem the eigenvalues of which X are the functions of Ca, &a, Ga and the wavenumber k.
Longwave instability

The eigenvalue problem for normal mode perturbations does not allow analytical solutions for arbitrary values of the governing parameters. Nevertheless, the case of long wavelength perturbations allows a thorough

Boussinesq Approximation for Fluid Systems

1163

analytical investigation. For k = 0 the solution describes neutral monotonous pe~~bations

given by (-1)).

0,

wj

= 0,

dj = yAjC(x-

(18)

At k = 0 all the other perturbations are damped. At small non-zero wavenumbers, the solutions are no longer neutral. Thus, to analyze the stability of eq~b~um with respect to long wavelength perturbations we need to find corrections to the X in terms of k. Thus, we represent the solution of spectral problem in terms of a power series in k for k << 1 we find for X X=k2X2+...,

X2zG

3~~J2[~]~~]~~+ 1 12144 + 5%d2 56[lcfPl 16 - 3[~]2

3( 16 - 4 3f~fr)

Ga

(19)

Thus, & depends linearly on Ra and Ga but the dependencies on the ratios of the viscosities yl/us, thermal conductivities IEI/K~ and thermal expansion coefficients /Jr/p2 are more complicated. Note that, the Prandtl number PT, the capillary parameter CQ, aud x1 and ~2 do not influence the stability of system with respect to long wave perturbations. For equal thermal conductivities, long-wave instabiity is absent (for Gs < 0). This is also true when the viscosities are not equal. Possibly this is due to the fact that the instability mechanism is connected with the difference in temperature gradients in the two fluids. The expression for X (19) determines the boundary of long-wave instability in the plane Gts - &a. This boundary is straight line that passes through the origin and the unstable region is situated to the right of the boundary (at large Ga). Long wave perturbations are damped to the left of the boundary. However, the system may be unstable to cellular perturbations with finite wave length in this region. To demonstrate this mechanism more clearly let us consider the case where the viscosities and the thermal expansion coefficients of two fluids coincide. As one can see from (19), in this case the increment does not depend on the sign of [K]. Let us assume that the thermal conductivity of lower fluid is much higher than that of upper fluid, so that {K]= 2, its maximal possible value. Then (19) reduces to

Under these conditions the temperature in any point of lower fluid coincides with the temperature at lower boundary. Consider the perturbations of the interface sketched in figure 1. At the points where the interface is not horizontal, the upper fluid is locally heated from the side which causes upward motion of fluid. The resulting flow is described qualitatively in the figure. Note that the flow is directed in such a way that the interface perturbations can be amplified. In the lower fluid the buoyancy forces is insuffiecient to cause motion and flow occurs only because of the motion of the upper fluid due to momentum transfer across the interface. The intensity of this flow is much less than that of the upper fluid and it cannot offset the amplification of the interface deformation due to the flow in the upper fluid. A maximum value of vertical component of the velocity is 2.3326 higher in the upper fluid. Its sign is the same aa the interface perturbation in upper fluid, while it has the opposite sign in lower fluid. The case when the liquids differ only by thermal conductivities and the liquid layers are of the same thickness is, in some sense, degenerated. Inspection of (19) reveals that As is proportional to the square of the thermal conductivity jump. As one can see from the same formulae, this degeneracy disappears when at least one more property (thermal expansion coefficient or the viscosity) differs between the two liquids. It is then possible that long wave instability occurs with heating from above. In a recent paper by Smith (1995) it is shown that this degeneracy also disappears when the layers have different thicknesses. Instability with respect to cellular perturbations

To investigate the stability of the system at finite wave numbers requires numerical calculations. To solve the boundary value problem for normal mode perturbations, a set of independent particular solutions was constructed, that satisfied rigid boundary conditions at z = fl. Together with the boundary conditions at

1164

D. V. Lyubimov er al.

the interface this leads to a characteristic equation for X that makes it possible to define stability boundaries of the system for )r = 0. Detailed calculations were carried out using the physical oil. The stabiity of this two-layer system was studied non-deformable interface. The parameters of this system -0,393 and Pr = 176. The boundary problem for normal a differential sweep method. properties of the system formic acid - insulating earlier in (Gershuni, 1982) under assumption of are: [K] = -0.838, [x] = -0.334, [v] = 1.756, p] 5: mode perturbations was solved numerically using

The stability diagram in the Ga - Ra plane for Ca = 0 is presented in figure 2. Stability boundaries with respect to the most dangerous monotonous (solid lines) and oscillatory (dashed lines) perturbations are shown. Curve 1 is the stability boundary with respect to cellular perturbations with finite k. With increasing IGal, curve 1 asymptotes to Ra B 282 . Here the critical wav~umber tends to k M 2.7. Recall that the case fGa/ -+ cc corresponds to au undeformable interface. The above mentioned asymptotic value of Ra is in a good agreement with the results of (Gershuni, 1982). The mechanism of instability along curve 1 appears to be the same as for an undeformable interfaces (Gershuni, 1982). It is interesting to note that the critical Rayleigh number increases with decreasing [Gal ( as interface deformability increases).
400 350 300 250 200 150 100 50 I I I I I I I 0 -40~35~30~25~20~15~100 -50 0 Ga

Ra

For Ga_ < Ga < 2.35 the most dangerous perturbations are travelling waves. Here the critical values of Ra decreases with decreasing IGal ( curve 2 bounds the range of instability from above). Note, that the phase velocity of the travelling waves is several orders less than the phase velocity of gravitational waves on the interface. In the range IGal < 2.35 monotonous perturbations again become the most dangerous. The co~espon~g stability boundary iu Ga - Ra plane is nearly a straight line (curve 3). It passes through the origin of the coordinate system (see inset) and, on the scale of figure 2, coincides with the Ra axis. Its slope is approximately equal to -38.8. For the case of heating from above (Ra < 0), the range of stability is limited from the right side by the &?a axis. This result is trivial: if the upper fluid is less dense (Ga < 0) then the fluid stratification is potenti~y stable. With Ga > 0 Rayleigh- Taylor insta~~ty occurs and the wavelength of the most dangerous perturbations is rather small. The reason for that is evident - at small Gu with heating from above unstable stratification is possible in a small region near the interface. This imposes the spatial scale of growing perturbations. At Ga + 0, the wavelength of most dangerous perturbations tends to zero (k + oo). The capillarity parameter apears in the boundary condition in the term Ga - k2Ca. It is then clear how the neutral curves Ga(k) should change as Ca increases, that is, the critical Galilei number increases more rapidly at larger wave numbers. For any Ru at high enough cap&&ties long-wave pert~bations are the most dangerous. For these, the minimum critical Rayleigh number is Ra w -68.15Ga. (21)

At intermediate values of capillary parameter the most dangerous perturbation is determined by the Rayleigh

Boussinesq Approximation for Fluid Systems

1165

number, At Ra < 80 oscillatory ~stab~ty is absent and at a certain critical value of the capillary parameter the instability to cellular perturbations changes to the long-wave instability. At Ra = 50 the change in the most dangerous perturbation takes place at Ca = 9. At high values of Ra where oscillatory instability can be most d~ge~u8, the long-wave monotono~ instabiity also becomes the most dangerous perturbation for certain values of Ca. At Ra > 282, and Ca = 0 a gap appears on the neutral curve. Here Ca > 0 will result in instability at all values of Ga. For realistic physical properties and laboratory conditions, the capillary parameter is usually large enough that the stability region on the diagram Ga - Ra is limited by two straight lines Ra TV282 and Ra = -68.15&z. Numerical calculations carried out at Ca = lo5 confirm this conclusion. 4. THE INFLUENCE OF BUOYANCY ON A LONGWAVE MARANGONI INSTABILITY

The final problem to be considered i8 the influence of buoy~cy on long-wave Mason insta~~ty. The existence of long-wave Marangoni instability for deformable free surface was first found by Takashima, (1981). The influence of buoyancy on this instability was previously studied by Wilson (1993) for a horizontal fluid layer of infinite extent with a free upper surface and a rigid lower boundary within the framework of formal Boussinesq approximation. Here we investigate this problem and account for non-Roussinesq effects The temperature at the bottom is constant, T = Zr, and the heat flux on the upper boundary is: --_- aT & - w2
- Go),

(22)

where 22is the free surface temperature, adis the thermal conductivity, h is the constant characterizing the heat flux, and T, is the temperature of the su~o~~ng media. We consider the fiuid to be isot~e~~y incompressible, i.e. the density is the function of temperature only. Iu the basic state the fluid is at rest, the upper boundary is flat and situated at a distance d from the lower one. The temperature To is linear function of z and z = 0 correspond8 to the lower boundary. That is, Tz - Tl To = Tr t zd, Ti - T2 = &(Tr - I).

(23)

The linearized equations for monotonous neutral perturbations

(X = 0) to this basic state are (24) (25) (26)

- Vp + Av f ;Vdivv + Ra pJ3gT=: 0, A8 = -w, d&v= -@w, which, together with the boundary conditions

2=

1:

W--O,

%=-Ma 1 PC

t&V

-~+GQj+zg=~~,

-g = sife - C)

define the perturbation amplitude8 at neutral stability, where v, w sre horizontal and vertical components of the velocity (we restrict ourselves to the analysis of 2-B perturbation8~, < is the free surface deviation from planarity. The reference quantities for length, temperature, velocity, pressure, density be giveu by d, B = TI - Ts, x/d, qx/da,and ps, respectively, the reference for the thermal expansion coefficient is P*, its value at the upper boundary temperature. The problem (24)-(29) contains the following dimensionless parameters: Rayleigh number Ra = g&@d3/vx (v = ~po), Bouseinesq parameter t: = j&f?, Marangoni number Ma = aaTgd/?jx, Galilei number Ga = gd3/vx, capillary parameter (Crispation number) C, = qx/crd, Biot number Bi = hd/K. (Yis the surface tension coefficient, CUT = -(l/cr)(dcz/dT), and A = (TI - Ta)/d. There is the relation between the parameters, namely Ra = &a.

1166

D. V. Lyubimov et al.

As noted in the introduction, the Boussinesq approximation is strictly valid in the limit 6 + 0, Ga -+ 00 at fixed Ra. Thus it follows from (29), that for this case one has to accept C = 0, i.e. the free surface must be considered undeformable. Another, often studied, case where the deformations of free surface can be taken into account is in the limit t: --) 0 at finite Ga. In this case Ra -P 0 and thermal convection does not arise. Thus, in this case, the Boussinesq approximation cannot be used for simultaneous consideration of buoyancy effects and the free surface deformation. From a formal point of view this means that the parameters 6 and Ga should be considered as independent. For normal mode perturbations proportions to eikz, where 1% is the wavenumber, and using w, 6 and C to denote the perturbation amplitudes, we obtain

(D2- k2)2w $ ~(0~ - k2)D(@w) - k2Ra p/30 = 0,

(30) (31) (32) (33)

(D2 - k2)(Dw+ q3w) - k2Gac - 2k2Dw=&k4<, where L? = #tk.

For the longwave instability we represent the solution as a series with respect to k2. At zeroth order we obtain we = 0, de = A&, (35) 1tBi At the next order we have the problem D3f = Ra pflt?~, z=o: .z=l: where we introduced the function f
f = Dw, -I-@wl.

(36) (37)

f=O, D2f = Ga <e.

Of=-----1fBi

M% b

(39)

The problem (36)-(38~ defines f uniquely (with accuracy up to the no~~zation factor fo). Equation (39f can be considered as the equation for WI, from which, taking into account the boundary conditions, we obtain the solvability condition

This condition defines the critical value of Marangoni number Mao as a function of the remaining parameters, that is,
Mao = -g

(c Si& t (1+ Bi&),

(41)

where
1

Fj =

eLJpd

- Z, fj(z)dz, f2 = 3, f3 = +,z2

(42)

Boussinesq Approximation for Fluid Systems

1167

and fi is determined by the problem Pfr = ppz, jr(O) = 0, Dfr(l) = D2j-r(1) = 0. (43)

As one can see from (41), Mao is proportional to Ga and is a linear function of Bi. For E = 0 we have F2 = l/2, F3 = -l/3, from which the previously obtained formula (Tahashima, 1981) follows Mao=
zGa(1 3 + Bi).

(44)

To analyze the influence of buoyancy on Marangoni instability we need to specify equations for density and @. Let is consider two different cases. In the first case, we assume the density to-be a linear function of temperature and take p = 1+ c(z - l), Clearly, (45) makes sense &ly if f < 1. Substituti~ (45) into (41) we obtain MTao= Ga 120(1+ Bi) - 3c(15 + 26Bi) + 139Bi , 3-c 60
In particular,

p = (1-t E(.z- 1))-

(45)

(46)

for thermally isolated boundaries (Bi = 0) we have 1 8-3~ Mao = pGU=*

(47)

As one might expect, heating from below (E > 0) leads to a decrease in stability whereas heating from above is stabilizing. In the second case we consider an alternative dependence of density on temperature p = c+-l), p = 1, of the form (48)

Substituting this into (41) we obtain the following formula for the critical Marangoni number

I,Bi
MQo=4 For Bi = 0

(-16 t 6csee + ll& - 14c& + 8c + 5e+) + 2e-$ - 4cec + 4c2 + 4~ * c2 (1 - ef + ce)

(49)

1 2+2r-2ef+c2et Mae== SGU c(l-eeC+ret) For small E (50) yields l-&&c2+... which coincides with (47) to first order in E. As for (46), it follows from (49) that for E > 0 stability decreases.

(50)

(51)

In this case, had we used the Boussinesq approximation in the usual way, that is, not treated Ra and E as independent parameters and, thus, dropped e: from the equations (30)-(34), we would have found the simple, but incorrect, formula:
Mao = %(40(1+

Bi) - 11cBi) = gGa(l+

Bi) - gRaBi

(52)

This formula differs from (49) already in first order of the expansion with respect to buoyancy parameter E.

1168

D. V. Lyubimov et al

For the case of thermally isolated boundaries it follows from (52) that ilfas does not depend on the Rayleigh number. But actually, as one can see it from (51), at Bi = 0 M&-j = ;Ga - &o + O(BQ2). (53)

That is, even for the case of thermaUy isolated boundaries taking into account the buoyancy effect leads to a decrease in stability. 6. CONCLUSIONS We have reconsidered the Boussinesq approximation and its application to systems with free and deformable interfaces and developed a gener~zation of the appro~mation which allows for so-consistent application to these systems. The generalization has been formulated for thermal buoyancy-driven convection in a system of two immiscible fluids separated by a deformable interface at which the densities of the fluids are comparable. Within the framework of this generalization we have examined two different problems. In the first problem we carried out an analysis of the linear stability of two superposed horizontal fluid layers of infinite lateral extent contained between two horizontal rigid boundaries separated by finite distance. Cases corresponding to heating from above and below were considered. At small values of the Bayleigh number, a longwave instab~ty was found for small values of the effective Gal&i number. At moderate vaIues of the GalliIei number is was found that osciUatory instability is the most dangerous, while at high GaILilei numbers the deformability of the interface odes not influence the stability of the system. The second problem concerned the case of a horizontaf layer of fluid with a free surface (at which a temperature dependent surface tension admitted the possibility of Marangoni convection) heated from below. It was found that accounting for the buoyancy effect led to a decrease in stability of the system with respect to the onset of Marangoni convection. It should be emphasized that this result could not be obtained within the framework of the conventional Boussinesq approximation.

ACKNOWLEDGMENTS The study was carried out with the financial support of International REFE~NCES Boussinesq, J., Theo&e analytique de la chafeur, Gauthier-Viiars, Paris Vol2, (1903). Ch~dr~kh~, S., ~yd~y~rn~c and ~yd~magnet~c Sta~~Z~ty, Clarendon Press, Oxford, Chapter 2 (1961). Davis, S. II. and G. M. Homsy, J. F&d Me&., 98, 527 (1980). Gershuni, G. Z. and E. M. Zhukhovitski, DokZ. AN SSSR, 265,302 (1982). Izakson, W. II, and V. I. Yudowich, frv. AN SSSR. Mekhanib zhidk, i gaza, 4,23 (1968). Nepomny~~, A. A., On the long wave convective insta~~ty in a horizon&I Iayer with deformable boundary, in Conuectiue flows, Perm, 25 (1983). Oberbeck, A., Annalen der Physic und Chemie, B7,271 (1879). Pavithran, S. and L. G. Redekopp, Eur. ,I. Mech. B/F&ds, lo,75 (1991). Renardy, Y., P&s. F&ids 29, 356 (1986). Smith, M. K., Phys. FZuids, (1995) . Takaahima, M., J. Phys. Sot. Japan, 60, 2745 (1981). Wilson, SK., J. Eng. Math., 27,161 (1993). Science Foundation (Grant d13 83).