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

G. Gouesbet, J. Maquet, C. Rozé, R. Darrigo

► To cite this version:

G. Gouesbet, J. Maquet, C. Rozé, R. Darrigo. Surface-tension and coupled buoyancy-driven instability in a horizontal liquid layer. Overstability and exchange of stability. Physics of Fluids, 1990, 2 (6), pp.903-911. 10.1063/1.857651. hal-01597366

HAL Id: hal-01597366

https://hal.science/hal-01597366

Submitted on 18 Feb 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Surface-tension- and coupled buoyancy-driven instability in a horizontal liquid layer. Overstability and exchange of stability

G. Gouesbet, J. Maquet, C. Rozé, and R. Darrigo

Laboratoire d'Energétique des Systèmes et Procédés, UA-CNRS 230, INSA de Rouen, BP8, 76131

Mont-Saint-Aignan Cedex, France

(Received 12 April 1988; accepted 7 February 1990)

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 concerning RBM (Rayleigh-Bénard-Marangoni) computations, anchored through a direct filiation in a line of works that goes back to the original theoretical master analysis of Rayleigh.¹ Landmark papers are by Pearson,² who showed that Bénard's cells³,⁴ were not buoyancy driven but actually tension driven, Scriven and Sternling,⁵ who addressed the question of what the roles of flexibility and resistance to deformation 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 authors. The most general work, akin to the present one, is that of Takashima,⁷ who extensively computed overstability and exchange of stability characteristics for surface-tension-driven instability in a horizontal liquid layer with a deformable 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 overstability characteristics would be modified if both Marangoni and Rayleigh numbers were simultaneously included in the analysis. Consequently, the present work can be considered as a significant 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 overstability.

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 crispation numbers, modifications induced by interfacial viscosities and heat transfer at the free surface, and different thermal 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 recovering known formulas and published numerical results; (ii) provide the reader with new expressions completing, in particular, Takashima's work⁷ 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 understanding, 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 other words, microgravity experiments permit the inhibition of exchange of stability but might promote the onset of overstability 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.

## II. THE OVERSTABILITY FORMULATION

### 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 Z<sub>fs</sub> of the free
surface is given by

$$
Z_{\text{fs}}(X,Y,t) = 1 + \delta Z_{\text{fs}}(X,Y,t), \quad (1)
$$

where $\delta Z_{\text{fs}}$ is the perturbation of the location of the free surface. The perturbation $\delta Z_{\text{fs}}$ may be expressed in terms of normal modes:

$$
\delta Z_{\text{fs}} = e^{\beta t} f(X, Y) \delta_s, \quad (2)
$$

in which *f* is the planform function (see Sec. II B) and satisfies a kinematic condition

$$
n_i u_{fs,i} = u_{fs,3} = \frac{\partial \delta Z_{fs}}{\partial t}, \quad (3)
$$

where $u_{fs,i}$ is the dimensionless velocity perturbation at the
free surface and $n_i$ is the normal unit vector to the free sur-
face, with the unit of time $t$ being $d^2/\alpha$ where $\alpha$ is the ther-
mal diffusivity of the liquid.

Assuming that the deformation of the free surface is
infinitesimal, 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,⁸ Scriven and Sternling,⁵ and Scriven,⁹ 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 normal force balance reads

$$
2H = \mathrm{Bo}(\sigma)\delta Z_{\mathrm{fs}} - \frac{\mathrm{Cr}(\sigma)}{\mathrm{Pr}}(p_1 - p_2) + 2 \mathrm{Cr}(\sigma) \frac{\partial U_3}{\partial Z}. \quad (4)
$$

The tangential force balance reads

$$
\Delta_{\text{II}} U_3 - \frac{\partial^2 U_3}{\partial Z^2} = \Delta_{\text{II}} \left( \frac{1}{\text{Cr}(\sigma)} \right) - \text{Vi} \Delta_{\text{II}} \frac{\partial U_3}{\partial Z} . \quad (5)
$$

The dimensionless mean curvature of the surface, *H*, at
a point (*X*, *Y*, *Z*<sub>fs</sub>), is given by

<figure><img src="image_2.png" /><figcaption>FIG. 1. Geometry of the problem.</figcaption></figure>

$$
H = \frac{1}{2} \Delta_{II} Z_{fs}, \quad (6)
$$

where $\Delta_{II}$ is the dimensionless surface Laplacian $\partial^2/\partial X_\alpha^2$, 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($\sigma$), Cr($\sigma$), and Vi are the Prandtl, Bond, crispation, and viscosity numbers, respectively, which are given by

$$
\mathrm{Pr} = \mu / \rho_0 \alpha, \quad (7)
$$

$$
\mathrm{Bo}(\sigma) = \rho_0 g d^2 / \sigma, \quad (8)
$$

$$
\mathrm{Cr}(\sigma) = \mu\alpha/\sigma d, \quad (9)
$$

$$
\mathrm{Vi} = (\kappa + \epsilon) / \mu d, \quad (10)
$$

where μ is the dynamic viscosity of the fluid, ρ₀ is the specific mass (at a reference temperature equal to the nonperturbed steady-state free surface temperature), σ is the surface tension, and (κ + ε) is the sum of the surface dilatational coefficient and of the surface shear viscosity coefficient (interfacial viscosities). Here, P₁ and P₂ are dimensionless pressures in media 1 and 2, respectively (Fig. 1), with ρ₀a²/d² as the pressure unit, and the dimensionless vertical component of the velocity of the liquid U₃ 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.¹⁰,¹¹ The base state is motionless with a dimensionless temperature gradient $\partial T / \partial X_j = (0,0,-k_T)$ in which dimensionless temperatures are produced by using the reference temperature. Then, we follow the standard procedure for linear analysis,⁵⁻⁷ set that space and time are separable in the perturbations ($p$ for pressure, $u_3$ for the $Z$ component of velocity, $\tau$ for temperature), and decompose spatial dependence in terms of normal modes:

$$
[p, u_3, \tau] = e^{\beta t} f(X, Y) [\nu(Z), \zeta(Z), \theta(Z)], \quad (11)
$$

where the planform function $f(X, Y)$ satisfies the Helmholtz equation

$$
\Delta_{\text{II}} f + a^2 f = 0, \quad (12)
$$

with $\beta$ and $a$ being the dimensionless time constant and dimensionless wavenumber of the normal mode, respectively, and obtain the velocity, temperature, and pressure equations given below:

$$
(\mathcal{D}^2 - a^2 r^2)(\mathcal{D}^2 - a^2)\zeta = (1/\text{Pr})\alpha_T g a^2 \theta, \quad (13)
$$

$$
(\mathcal{D}^2 - a^2 q^2)(\mathcal{D}^2 - a^2 r^2)(\mathcal{D}^2 - a^2)\theta = -a^2 \text{Ra } \theta, \quad (14)
$$

$$
(\mathcal{D}^2 - a^2) \nu = \alpha_T g \mathcal{D} \theta. \quad (15)
$$

Using a momentum equation, the pressure equation can
be modified to

$$
a^2 \nu = \mathrm{Pr}(\mathcal{D}^2 - a^2 r^2) \mathcal{D} \xi \quad (16)
$$

and the temperature and velocity Z perturbations are related
by
$$ (\mathcal{D}^2 - a^2 q^2)\theta = -k_T \xi, \quad (17) $$

where the dimensionless acceleration vector $g_i$ has been taken equal to $(0,0,-g)$ and the operator $\mathcal{D}$ is $d/dZ$. The expansion factor $\alpha_T$, the Rayleigh number $Ra$, $r^2$, and $q^2$, are given by

$$ \alpha_T = - \frac{1}{\rho_0} \left. \frac{\partial \rho}{\partial T} \right|_{\rho=\rho_0}, \quad (18) $$

$$ \mathrm{Ra} = (1/\mathrm{Pr})\alpha_T k_T g, \quad (19) $$

$$ r^2 = 1 + \beta/a^2 \mathrm{Pr}, \quad (20) $$

$$ q^2 = 1 + \beta/a^2. \quad (21) $$

### C. The three-constant solution for $\theta$

The $\theta$ solution of Eq. (14) has been researched under the form

$$ \theta = \sum_{k=0}^{\infty} \mathcal{C}_k Z^k, \quad (22) $$

subject to a double kinematic condition at the solid surface ($\xi = \mathcal{D}\xi = 0$, $Z=0$) and to a thermal condition, either $\theta=0$, $Z=0$ (conductor case) or $\mathcal{D}\theta=0$, $Z=0$ (insulating case). Formulas are given in the present paper only for the conductor case (expressions are quite similar for the alternative thermal condition).

We then find that $\theta$ is the sum of an odd and an even series, with three arbitrary constants chosen to be $\mathcal{C}_1, \mathcal{C}_4, \mathcal{C}_5$. The results can be expressed as follows:

$$ \mathcal{C}_{2k} = [-4!(2k)!]e_4^{2k}\mathcal{C}_4, \quad k \ge 0, \quad (23) $$

$$ \mathcal{C}_{2k+1} = [-5!(2k+1)!]e_5^{2k+1}\mathcal{C}_5 - [1!(2k+1)!]e_1^{2k+1}\mathcal{C}_1, \quad k \ge 0, \quad (24) $$

with

$$ e_4^0 = e_4^2 = e_5^1 = e_5^3 = e_5^5 = 0, \quad (25) $$

$$ e_4^4 = e_1^1 = e_5^5 = -1, \quad (26) $$

$$ e_1^3 = -a^2 q^2, \quad (27) $$

and

$$ (e_4^{2k+6}, e_i^{2k+7}) = -a_4(e_4^{2k+4}, e_i^{2k+5}) - a_2(e_4^{2k+2}, e_i^{2k+3}) - a_0(e_4^{2k}, e_i^{2k+1}), \quad i = (1,5), \quad k \ge 0, \quad (28) $$

where

$$ a_4 = -a^2(1 + r^2 + q^2), \quad (29) $$

$$ a_2 = a^4(r^2 + r^2 q^2 + q^4), \quad (30) $$

$$ a_0 = a^2(\mathrm{Ra} - a^4 q^2 r^2). \quad (31) $$

### 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 $\theta$ plus a fourth one for the free surface location. This requires four boundary conditions at the free surface. These conditions can be obtained by combining Refs. 5 and 7. The first condition expresses the kinematic condition (3). The second and third conditions express 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

$$ (\mathcal{D}^2 - a^2 q^2)\theta = -k_T \delta_s \beta, \quad (32) $$

$$ \mathrm{Cr}[\mathcal{D}^2 - a^2(r^2 + 2)](\mathcal{D}^2 - a^2 q^2)\mathcal{D}\theta + a^2 k_T \delta_s (a^2 + \mathrm{Bo}) = 0, \quad (33) $$

$$ (\mathcal{D}^2 + a^2)(\mathcal{D}^2 - a^2 q^2)\theta + a^2 \mathrm{Vi}(\mathcal{D}^2 - a^2 q^2)\mathcal{D}\theta - a^2 \mathrm{Ma}\theta + a^2 k_T \mathrm{Ma}\delta_s = 0, \quad (34) $$

$$ \mathcal{D}\theta + \mathrm{Bi}\theta - \mathrm{Bi}k_T\delta_s = 0, \quad (35) $$

where $\mathrm{Cr}$ and $\mathrm{Bo}$ are taken as constants by writing $\sigma = \sigma_s$ (neglecting higher-order terms), $\mathrm{Bi}$ is the Biot number for heat transfer at the free surface (sometimes called the Nusselt number by some authors), and $\mathrm{Ma}$ is the Marangoni number that may be expressed as

$$ \mathrm{Ma} = -\Gamma k_T, \quad (36) $$

$$ \Gamma = \left. \frac{\partial}{\partial T} \frac{1}{\mathrm{Cr}(\sigma)} \right|_{T=T_s}, \quad (37) $$

where $T_s$ is the dimensionless free surface temperature in the basic state. In conditions (32)–(35), it is now advantageous to replace $k_T\delta_s$ by a single unknown $\delta$.

### E. The characteristic equations and the problem solution for overstability

We introduce the free surface boundary conditions into the three- constant solution for $\theta$. In the case of nontrivial solutions, the characteristic determinant of the resulting homogeneous linear set of equations must be zero, leading to

$$ \mathrm{Ma} = N/D, \quad (38) $$

where

$$ N = \begin{bmatrix} r_j & r_4 & r_5 & \beta/a^2 \\ s_j & s_4 & s_5 & -(\mathrm{Bo} + a^2)/\mathrm{Cr} \\ (t_j + u_j \mathrm{Vi}) & (t_4 + u_4 \mathrm{Vi}) & (t_5 + u_5 \mathrm{Vi}) & 0 \\ (x_j + y_j \mathrm{Bi}) & (x_4 + y_4 \mathrm{Bi}) & (x_5 + y_5 \mathrm{Bi}) & \mathrm{Bi}/a^2 \end{bmatrix}, \quad (39) $$

$$ D = \begin{bmatrix} r_j & r_4 & r_5 & -\beta/a^2 \\ s_j & s_4 & s_5 & (\mathrm{Bo} + a^2)/\mathrm{Cr} \\ v_j & v_4 & v_5 & 1 \\ (x_j + y_j \mathrm{Bi}) & (x_4 + y_4 \mathrm{Bi}) & (x_5 + y_5 \mathrm{Bi}) & -\mathrm{Bi}/a^2 \end{bmatrix}, \quad (40) $$

in which
$$j=1 \text{ for the conductor case,} \tag{41}$$

$j=0$ for the insulating case.

Expressions for $r_0, \dots, y_5$ are not given in this paper. We only mention that they appear similar to series whose convergence is ensured by factorial terms in denominators. For the onset of overstability, $\beta$ is $i\omega$ (marginal oscillatory instability). 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 wavenumber $a$, the critical Marangoni number $Ma_c$ is obtained from the following equations:

$$\mathrm{Im}[\mathrm{Ma}(\beta, a)] = 0, \tag{42}$$

$$Ma_c = Ma(\beta_c, a). \tag{43}$$

Relation (42), where $Ma$ designates the pseudo-Marangoni number, determines a spectrum $\beta_i$ of time constants for which the pseudo-Marangoni numbers are real. The critical time constant $\beta_c$ is the value among the $\beta$-spectrum values for which $|Ma|$ is the smallest. The critical Marangoni number $Ma_c$ is then given by relation (43). Scanning over $a$, we search for the critical wavenumber $a_*$ for which $|Ma_c(a)|$ is a minimum. The final critical time constant $\beta_*$ and Marangoni number $Ma_*$ are $\beta_c(a_*)$ and $Ma_c(a_*)$, respectively.

### III. THE EXCHANGE OF STABILITY FORMULATION

#### A. General formulation

This formulation is derived as a special case of the overstability formulation given in Sec. II. Time constant $\beta$ 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 $N$ and $D$, and expressions for $r_0, \dots, y_5$ must be rewritten by specifying the new relation $r^2 = q^2 = 1$. The relations defining the coefficients $e_i^j$ (25)–(31) also simplify. Finally, (42) becomes an identity because the pseudo-Marangoni number is now a real number. For a given wavenumber $a$, we obtain the critical Marangoni number $Ma_c = Ma(a)$, from which we determine $a_*$ and $Ma_* = Ma_c(a_*)$ by scanning over $a$.

#### B. Pure Marangoni effect ($Ra=0$)

In this case, the formulation simplifies dramatically. We find that coefficients $e_i^j$ can be obtained in plain explicit forms, without any recurrence,

$$e_1^{2k+1} = [\frac{1}{2}k(k-1) - 1]a^{2k}, \quad k \ge 0, \tag{44}$$

$$e_4^{2k} = e_5^{2k+1} = -\frac{1}{2}k(k-1)a^{2(k-2)}, \quad k \ge 0. \tag{45}$$

Injecting these expressions into the expressions for $r_0, \dots, y_5$, we observe that the involved series can be reduced to plain mathematical functions. For instance, $r_0$ becomes $[a^3 \sinh(a)]/2$.

For the conductor case, (38)–(40) then lead to

$$Ma_c = Ma(a) = \frac{8a[(CS - a) + \frac{1}{2}a \mathrm{Vi}(S^2 - a^2)](aC + \mathrm{Bi} S)(\mathrm{Bo} + a^2)}{8 \mathrm{Cr} a^5 C + (\mathrm{Bo} + a^2)(S^3 - a^3 C)}, \tag{46}$$

where $S = \sinh(a)$ and $C = \cosh(a)$. With condition $\mathrm{Vi} = 0$, used by Takashima, (46) becomes identical to relation (3.1) of Takashima⁷ (p. 2745) as it should.

Similarly, for the insulating case (not considered by Takashima), we obtain

$$Ma_c = Ma(a) = \frac{8a[(CS - a) + \frac{1}{2}a \mathrm{Vi}(S^2 - a^2)](aS + \mathrm{Bi} C)(\mathrm{Bo} + a^2)}{8 \mathrm{Cr} a^5 S + (\mathrm{Bo} + a^2)(CS^2 - 2aS + a^2C - a^3S)}, \tag{47}$$

### IV. OVERSTABILITY: NUMERICAL RESULTS

#### A. Generalities

To solve the set (42) and (43) for overstability, we relied 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 algorithm. Outputs are critical quantities $a_*, \omega_*$ (i.e., $\beta_*/i$) and $Ma_*$ by scanning over $a$. Typically, the CPU time needed to obtain a set ($a_*, \omega_*, Ma_*$) is 10 min on a IBM 3090 mainframe and the number of terms in series $r_0, \dots, y_5$ is less than about 150 for an accuracy of about 0.05% on $Ma_*$ although these figures 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 carried

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 Takashima considered the simpler case $\mathrm{Vi} = \mathrm{Bi} = \mathrm{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. However, this is not enough to dismiss the possibility of overstability for a hot rigid wall and coupled Marangoni–Rayleigh effects, because it is impossible to scan all the ranges of thermophysical inputs and program driving parameters numerically. We might simply have been unlucky. Evidence for overstability impossibility could only be obtained from a formal proof. Nevertheless, we strongly believe in it. Because, if we accept Takashima's conclusion in the case of a pure tension mechanism, the occurrence of overstability, when the buoyancy mechanism is also introduced, is even more difficult
<figure><img src="image_3.png" /><figcaption>FIG. 2. Situations for possible investigation (Γ assumed negative). For the sign of Rayleigh number, $\alpha_T$ is assumed positive. For a liquid having a negative coefficient of volume expansion, Ra signs must be reversed.</figcaption></figure>

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 $\omega_*$ and $a_*$ for a value $X_0$ of a parameter to investi- gate (for instance, viscosity number). Critical values for an- other value $X_1 = X_0 + \delta X$, with an increment $\delta X$ not too large, are computed by investigating a range [$a_1, a_2$] con- taining $a_*$ and a range [$\omega_1, \omega_2$] containing $\omega_*$. The use of such an “adiabatic” modification of parameter $X$ from an initial value $X_0$ to a final value $X_f$ to track the critical values is not fully secure because the modification 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 first set of

<figure><img src="image_4.png" /><figcaption>FIG. 3. Overstability. Influence of the Rayleigh number on the critical wavenumber $a_*$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIG. 4. Overstability. Influence of the Rayleigh number on the critical time constant $\omega_*$.</figcaption></figure>

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 $|Ma_*|$, 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 modification of the thermal condition at the rigid wall will be briefly dis- cussed at the end of the next section (for convenience). We now turn to the influence of the Rayleigh number for which the continuation method has been found to fail in some cases.

## B. Influence of the Rayleigh number

Results are shown in Figs. 3–5, for the Fig. 2(b) case, conductor wall, Vi = Bi = 0, Pr = 1, Bo = 0.01, and Cr ranging from $10^{-7}$ to $10^{-1}$.

<figure><img src="image_6.png" /><figcaption>FIG. 5. Overstability. Influence of the Rayleigh number on the critical Mar- angoni number $Ma_*$.</figcaption></figure>
For Ra = 0, they perfectly reproduce results given by Takashima.⁷ 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 $|Ma_c(a)|$ for $a = a_1$ and $a = a_2$. The true critical quantities correspond, for instance, to $a_* = a_1$, for which $|Ma_c(a_1)| < |Ma_c(a_2)|$. When the crispation number increases, $|Ma_c(a_1)|$ increases while $|Ma_c(a_2)|$ decreases. At the point of discontinuity, we have $|Ma_c(a_1)| = |Ma_c(a_2)|$. On the right of discontinuity, the situation is the opposite with $|Ma_c(a_1)| > |Ma_c(a_2)|$. Discontinuity corresponds to a jump of the critical wavenumber from $a_* = a_1$ to $a_* = a_2$ (Fig. 3), accompanied by a jump of the time constants from $\omega_* = \omega_1$ to $\omega_* = \omega_2$ (Fig. 4). Clearly, such a process does not produce any discontinuity in the $Ma_*$ results (Fig. 5).

For $Ra = -10^2$, 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 $\approx -2.5$. For $Ra = -3 \times 10^2$, discontinuity has again drift-ed to the right, but its importance is smaller. Smoothing of this discontinuity is carried out at $Ra = -10^3$, but memory of it is kept through a fairly steep curvature of the profiles for both the wavenumbers and the time constants. This memory is completely lost at $Ra = -10^4$. Nothing special occurs between $Ra = -10^{-4}$ and $Ra = -10^5$, except a decrease in the critical wavenumbers, an increase in the critical time constants, and an increase in the critical $|Ma_*|$. At $Ra = -10^6$, the broken lines in Figs. 3-5 correspond to the critical values we obtain by using the continuation method (Sec. IV A). However, these values are not truly the critical ones because of a structural modification of the stability dia-grams, occurring between $Ra = -10^5$ and $-10^6$. This modification is shown in Fig. 6 for $Cr = 10^{-5}$. We give $Ma_c(a)$ for various values of the Rayleigh number ranging from $-10^5$ to $-10^6$. At $Ra = -10^5$, the curve $|Ma_c(a)|$

<figure><img src="image_7.png" /><figcaption>FIG. 6. Overstability. Modifications of the curve $Ma_c(a)$ vs $Ra$ for $Cr = 10^{-5}$.</figcaption></figure>

presents one minimum corresponding to the points labeled 1 in Figs. 3-5. For $Ra = -2 \times 10^5$, this minimum has drifted to smaller $a$ values and to larger $|Ma_c|$ 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 profile labeled 2. When $Ra$ ranges from $-10^5$ to $-10^6$, the first minimum evolves along the broken line indi-cated in Fig. 6. Somewhere between $Ra = -2 \times 10^5$ and $Ra = -5 \times 10^5$ the critical values jump from the first mini-mum to the second minimum. According to the continu-ation method, the critical values at $Ra = -10^6$ would be found to correspond to the first minimum in Fig. 6 (the curve labeled 4), leading to the points labeled 2 in Figs. 3-5. But the true critical values actually correspond to the second minimum. These values are given for $a_*$ and $Ma_*$ in Figs. 3 and 5 (full lines). The second minimum in Fig. 6 (the curve labeled 4), corresponds to the points labeled 3. The corre-sponding $\omega_*$ values are too high to be indicated in Fig. 4. For instance, for $Cr = 10^{-5}$, $Ra = -10^6$, we have $\omega_* \approx 140$. In Fig. 3, we observe a discontinuity for $log Cr$ between $-5.8$ and $-5.9$. The $a_*$ and $\omega_*$ values to the left of this discontin-uity are also too high to be indicated in the figures. For in-stance, for $log Cr = -6.6$, $Ra = -10^6$, we have $a_* = 2.7$ and $\omega_* \approx 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.¹² They may also be connected with the fact that, according to Davis and Homsy,¹³ 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.⁵,¹⁴ Although this last statement concerns the case of exchange of stability, it indi-cates a conflict 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 Fig. 5 reveals a small region, on the curve $Ra = -10^3$, for $log Cr$ between $\approx (-3)$ and $\approx (-2.7)$, where buoyancy lowers the critical $|Ma_*|$. This region corresponds to the disappearance of discontinuity ob-served in Figs. 3 and 4 when $Ra$ decreases from $0$ to $-10^3$.

Finally, the influence of a modification of the thermal condition on the critical Marangoni numbers is not very sig-nificant for $Pr = 1$, $Bo = 0.01$, $Vi = Bi = Ra = 0$, $Cr rang-ing between $10^{-7}$ and $10^{-1}$. The only noticeable effect for the insulating case is the disappearance of discontinuities shown in Figs. 3 and 4 (detailed results not provided).

## 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 Nield⁶ and Takashima,⁷
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
specifications, 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 specifications since the sys-
tem is stable, or at least is difficult insofar as loss of stability
may arise if the specifications 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 (47)] from
which numerical results can easily be obtained and discussed
[see Ref. 7, p. 2745, for relation (46) with $Vi = 0$. Conse-
quently, results will only be presented for (i) the influence of
the Rayleigh number ($Ra \neq 0$) and (ii) the influence of the
viscosity number, which was not considered in Ref. 7. The
presentation of numerical results is also limited by consider-
ing neither the insulating case nor positive Bond numbers.

For the conductor condition at the rigid wall, Takashima computed the critical Marangoni number $Ma_c(a)$ for $Vi = Bi = Ra = 0$, $Bo = -0.1$, and various values of the crispation number. Under the same specifications, we show our results in Fig. 7 and also in Fig. 8, where they are repeated. The points below each curve represent stable states. The agreement with Takashima's results is perfect (see Fig. 4 in Ref. 7). These results have been obtained from the full computer program but could also be obtained from (46) as done by Takashima. By examining this relation, Takashima observed that all of the neutral stability curves intersect the *a* axis when $a = (-Bo)^{1/2}$ (here $a \approx 0.316$). When $a > (-Bo)^{1/2}$, all of the curves lie in the positive region of $Ma_c$. Conversely, for $0 < a < (-Bo)^{1/2}$, they lie in the negative region of $Ma_c$ as shown in Figs. 7 and 8 ($0 < a < 0.3$). The region below each curve represents stable states; the critical Marangoni number $Ma_c$ 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 specifications. The loss of stability now observed is attributed to the introduction of gravity waves, although the free surface deformation (crispation number) also plays an essential role. Takashima also observed that when $Cr > 0.000085$, the critical $a_*$ is zero. Otherwise, critical $a_*$ becomes nonzero and $Ma_c$ decreases rapidly as $Cr$ decreases (see Figs. 7 and 8). These numerical observations are confirmed 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 < Ma_c$ and Takashima states that cellular convection is also, in principle, possible when $Ma = Ma_c$. Conversely, for

<figure><img src="image_8.png" /><figcaption>FIG. 7. Neutral stability curves for the conductor case, $Bo = -0.1$, $Bi = Ra = 0$, and various values of the crispation number: $Cr = 10^{-3}$ (label 1), $3 \times 10^{-4}$ (label 2), $10^{-4}$ (label 3), $5 \times 10^{-5}$ (label 4), $3 \times 10^{-5}$ (label 5), $2.3 \times 10^{-5}$ (label 6), $2.1 \times 10^{-5}$ (label 7). Influence of the viscosity number: heavy lines: $Vi = 0$; thin lines: $Vi \neq 0$. With $Vi = 10$ (label a), 25 (label b), and 50 (label c). Example: The curve 7c is for $Cr = 2.1 \times 10^{-5}$, $Vi = 50$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 8. Neutral stability curves for the conductor case, $Bo = -0.1$; $Vi = Bi = 0$; and various values of the crispation number: $Cr = 10^{-3}$ (label 1), $3 \times 10^{-4}$ (label 2), $10^{-4}$ (label 3), $5 \times 10^{-5}$ (label 4), $3 \times 10^{-5}$ (label 5), $2.3 \times 10^{-5}$ (label 6), and $2.1 \times 10^{-5}$ (label 7). Influence of the Rayleigh number: heavy lines: $Ra = 0$; thin lines: $Ra \neq 0$ with $Ra = 500$ (label a), 800 (label b), 2500 (label c), and 4000 (label d).</figcaption></figure>
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), interfacial viscosities can provoke a modification of the critical wavenumber $a_*$ from zero to a finite value (label 3c for instance). At small crispation numbers (labels 5–7) the modification of the finite critical wavenumber $a_*$ is real but less and less significant when Cr decreases. The points below each curve represent stable states, and we are led to the conclusion that the interfacial viscosities have (in this case) a destabilizing character, an unexpected result at first sight. The explanation is (for instance) as follows.

If the value of the absolute temperature gradient (or of $|Ma_c|$) 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 gradient. 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 observed 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 modification in the figure, which would, however, remain small. When

the crispation number decreases, the effect of buoyancy becomes more important. For curves labeled 6 and 7, the induced modifications are very important even at $Ra = 500$. This can be physically understood if we interpret the decrease 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 $Ma_c$ changes its sign (becoming positive) in the $a$ range under study. The critical wavenumbers $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 influence of the viscosity number is displayed in Figs. 9 and 11 and the influence of the Rayleigh number in Figs. 10 and 12. Comments are not given because these results merely reflect the data in Figs. 7 and 8.

<figure><img src="image_10.png" /><figcaption>FIG. 10. The critical Marangoni number $Ma_*$ versus the crispation number for the conductor case, $Bo = -0.1$, $Vi = Bi = 0$. Influence of the Rayleigh number.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIG. 9. The critical Marangoni number $Ma_*$ versus the crispation number for the conductor case, $Bo = -0.1$, $Bi = Ra = 0$. Influence of the viscosity number.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIG. 11. The critical wavenumber $a_*$ versus the crispation number for the conductor case, $Bo = -0.1$, $Bi = Ra = 0$. Influence of the viscosity number.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>FIG. 12. The critical wavenumber $a_*$ versus the crispation number for the conductor case, Bo = -0.1, Vi = Bi = 0. Influence of the Rayleigh number.</figcaption></figure>

## VI. CONCLUSION

Overstability and exchange of stability for simultaneous surface-tension and buoyancy-driven instability in a horizontal infinite layer have been theoretically investigated by means of a small disturbance analysis. Formulations and expressions are given in dimensionless forms. Numerical results are obtained relying on computer programming validated by comparisons with the results of Takashima<fnref n="7" /> and Nield.<fnref n="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, except 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 formulation can be simplified 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 wavenumbers $a_*$ versus the crispation number, for a given Bond number, with a discussion of the influence of interfacial viscosities and of the buoyancy instability mechanism. To mention only one point in this conclusion, we observed that the interfacial viscosities may have a destabilizing effect, which, at first sight, is an unexpected result.

We finally state that more details for the formulation and the computing algorithms, listings, and additional results 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'Université de Rouen.

<footnote n="1">¹Lord Rayleigh, Philos. Mag. **32**, 529 (1916).</footnote>

<footnote n="2">²J. R. A. Pearson, J. Fluid Mech. **4**, 489 (1958).</footnote>

<footnote n="3">³H. Bénard, Rev. Gen. Sci. Pures Appl. **11**, 1261 (1900).</footnote>

<footnote n="4">⁴H. Bénard, Ann. Chim. Phys. **23**, 62 (1901).</footnote>

<footnote n="5">⁵L. E. Scriven and C. V. Sternling, J. Fluid. Mech. **19**, 321 (1964).</footnote>

<footnote n="6">⁶D. A. Nield, J. Fluid Mech. **19**, 341 (1964).</footnote>

<footnote n="7">⁷M. Takashima, J. Phys. Soc. Jpn. **50**, 2745, 2751 (1981).</footnote>

<footnote n="8">⁸R. Aris, *Vectors, Tensors and the Basic Equations of Fluid Mechanics* (Prentice-Hall, Englewood Cliffs, NJ, 1982).</footnote>

<footnote n="9">⁹L. E. Scriven, Chem. Eng. Sci. **12**, 98 (1960).</footnote>

<footnote n="10">¹⁰D. D. Joseph, *Stability of Fluid Motions* (Springer, Berlin, 1976), Vols. I and II.</footnote>

<footnote n="11">¹¹*Convective Transport and Instability Phenomena*, edited by J. Zierep, and H. Ortel, Jr. (Braun, Karlsruhe, FRG, 1982).</footnote>

<footnote n="12">¹²A. Pellew and R. V. Southwell, Proc. R. Soc. London Ser. A **176**, 321 (1940).</footnote>

<footnote n="13">¹³S. H. Davis and G. M. Homsy, J. Fluid Mech. **98**, 527 (1980).</footnote>

<footnote n="14">¹⁴W. V. Kayser and J. C. Berg, J. Fluid. Mech. **57**, 739 (1973).</footnote>