53
AUTHOR’S POST PRINT (Romeo Colour: Green)
Int. J. Num. Meth. Fluids (ISSN: 0271-2091), 36 (1): 53-90 (2001).
DOI: 10.1002/fld.120
Publisher version available at
http://onlinelibrary.wiley.com/doi/10.1002/fld.120/abstract
Three-dimensional numerical simulation of Marangoni
instabilities in liquid bridges: influence of geometrical
aspect ratio
M. Lappa
*
, R. Savino and R. Monti
Università degli Studi di Napoli "Federico II"
Dipartimento di Scienza e Ingegneria dello Spazio "Luigi G. Napolitano", Università degli Studi di Napoli "Federico II",
Napoli, Italy
*
current e-mail address: marcello.lappa@strath.ac.uk
SUMMARY
Oscillatory Marangoni convection in silicone oil liquid bridges with different geometrical aspect ratios is
investigated by three-dimensional and time-dependent numerical simulations, based on control volume
methods in staggered cylindrical non uniform grids. The three-dimensional oscillatory flow regimes are
studied and compared with previous experimental and theoretical results. The results show that the
critical wave number (m), related to the azimuthal spatio-temporal flow structure, is a monotonically
decreasing function of the geometrical aspect ratio of the liquid bridge (defined as ratio of the length to
the diameter). For this function a general correlation formula is found, that is in agreement with the
previous experimental findings. The critical Marangoni number and the oscillation frequency are
decreasing functions of the aspect ratio; however, the critical Marangoni number, based on the axial
length of the bridge, does not change much with the aspect ratio. For each aspect ratio investigated, the
onset of the instability from the axi-symmetric steady state to the three-dimensional oscillatory one is
characterized by the appearance of a standing wave regime that exhibits, after a certain time, a second
transition to a travelling wave regime. The standing wave regime is more stable for lower aspect ratios
since it lasts for a long time. This behaviour is explained on the basis of the propagation velocity of the
disturbances in the liquid phase; for this velocity a general correlation law is found as function of the
aspect ratio and of the Marangoni number.
Key words: liquid bridge, Navier-Stokes calculations, Marangoni flow.
1. INTRODUCTION
Considerable attention has been paid in recent years to surface tension driven (Marangoni)
convection in liquid floating zones and similar configurations under microgravity conditions.
The reason for this interest is that, once buoyancy convection is reduced by several orders of
54
magnitude in space, surface tension-driven flows are the only mechanisms responsible for the
appearance of undesirable instabilities in liquid semiconductors or other liquid melts materials
during container-less floating zone crystal growth methods.
For a better understanding and optimization of these processes, model problems have been
formulated and studied analytically, numerically and experimentally. The most common model
is that of the half zone consisting of two circular, coaxial disks at different temperatures, with a
bridge of liquid suspended between them. Since the end of 1970's (Chun and West
1,2
), on ground
experiments with small size half zones and high Prandtl number liquids (Pr>>1) pointed out that
the Marangoni convection in liquid bridges heated from above or from below shows a transition
from a steady axisymmetric toroidal flow to a three-dimensional oscillatory flow, when a critical
temperature difference between the liquid bridge supports is exceeded. Preisser, Schwabe and
Scharmann
3
, Velten, Schwabe and Scharmann
4
and Frank and Schwabe
5
performed
systematically on ground research over a large range of experimental conditions. They found that
the oscillatory Marangoni flow exhibits different behaviour depending on the Prandtl number of
the liquid employed and on the geometrical aspect ratio of the zone, defined as the ratio of the
length L and of the diameter D of the bridge (A=L/D).
The onset of instability of Marangoni convection was investigated in space using large liquid
bridges with length of the order of several centimeters and in absence of buoyancy effects by
Monti
6
.
Experimental works have been performed recently on ground by Petrov et al.
7
, Muehlner et al.
8
,
Schwabe et al.
9
.
In the last years the development of supercomputers and efficient numerical methods led the
investigators to study the problem through numerical solution of the three-dimensional and time-
dependent Navier Stokes equations.
Rupp, Mueller and Neumann
10
studied by three-dimensional numerical simulations the
instability behaviour for different Prandtl number liquids and for a fixed aspect ratio (A=0.6).
Stability analyses were carried out by several authors
11-18,25
to define, in the non-dimensional
parameter space, sufficient conditions for stability and instability.
The most complete results on the subject are those reported in the linear stability analyses of
Kuhlmann
15
and Kuhlmann and Rath
16
, Chen and Roux
17
and Chen and Hu
18
. The basic
axisymmetric Marangoni flow was obtained by numerical solutions of the Navier-Stokes
equations, together with the appropriate boundary and symmetry conditions, and the eigenvalue
problem for the three-dimensional disturbances was solved over a range of Prandtl numbers and
for aspect ratios close to unity. These results predicted the critical Marangoni numbers and the
form of the most dangerous disturbances, characterized by the appropriate value of the critical
wave number, in the neighbourhood of the neutral stability point (i.e., at the onset).
Kuhlmann and Rath
16
found that the three dimensional supercritical state after the Hopf
bifurcation point, is given by a superposition of two counter-propagating hydrothermal waves,
with axial and azimuthal components. If the two waves have equal amplitude the resulting
disturbance is a "standing wave", with the minimum and maximum disturbances pulsating at
fixed azimuthal positions; while the superposition of waves with different amplitude gives rise to
a "travelling wave", with the minimum and maximum disturbances travelling in azimuthal
direction. Unfortunately, nothing can be predicted by the linear stability analyses about the
amplitude of these disturbances.
55
The experiments performed on ground and in microgravity conditions revealed one or the other
instability mechanism. For instance the travelling wave has been observed by Chun et al.
1
,
Preisser et al.
3
, Frank et al.
5
, Monti
6
,
Muehlner et al.
8
and Schwabe et al.
9
. The standing wave,
instead, has been observed during experimental ground-based activities (Velten et al.
4
, Frank et
al.
5
).
Recent numerical results obtained by Savino, Monti
19
and Monti, Savino, Lappa
20
, have shown
that the standing wave and the travelling wave models correspond to two consecutive transitions
of the Marangoni flow. In particular, for large Prandtl number liquids, the flow exhibits a first
transition from the axi-symmetric steady to the three-dimensional oscillatory state, characterized
by the standing wave instability and, after a certain time, a second transition from the standing
wave to the travelling wave. These results, for a liquid bridge with a fixed aspect ratio (A=1),
have been also validated by experimental results on ground (Monti, Savino, Lappa
21
) where, for
the first time the transition from one regime to the other was clearly observed.
In a previous work
22
a parallel solution numerical method has been introduced and applied to
study three-dimensional Marangoni flow instabilities in liquids with low Prandtl number.
In this paper the method is extended to high Prandtl numbers for which most of the existing
information comes from experiments and linear stability calculations. In particular the attention
is focused on the influence of the aspect ratio on the flow instability, on the critical wave number
(m) and on the spatio-temporal structures (pulsating or rotating).
2. PHYSICAL AND MATHEMATICAL MODEL
2.1 Basic assumptions
Fig.1 shows the geometry of the problem and the boundary conditions. A cylindrical liquid
bridge is suspended between two coaxial disks with constant temperatures (
TT T
O
/2
),
where
T
0
is the ambient temperature and T the overall temperature difference.
T
h
=T
c
+
T
z
r
u
v
v
T
c
Fig.1: sketch of the liquid bridge
56
The liquid is assumed homogeneous and Newtonian, with constant density and transport
coefficients. The bridge is bounded by a cylindrical and undeformable liquid-gas interface with
a surface tension exhibiting a linear decreasing dependence on the temperature:
=
o
-
T
(
T
-T
0
), (1)
where
o
is the surface tension for TT
0
;
T
is the negative rate of change of the surface
tension with temperature
(
T
=-d/d
T
> 0).
The hypothesis of rigid free surface is acceptable if zero-g conditions prevail or, on the ground,
if the volume of the liquid bridge is small and if the geometrical aspect ratio is sufficiently low.
2.2 Nondimensional field equations and boundary conditions
With the above assumptions the flow is governed by the continuity, Navier-Stokes and energy
equations, that in non-dimensional conservative form read :
V0
(2a)
V
t
pVV V Pr
2
(2b)
T
t
VT T
2
(2c)
where V, p and T are the non-dimensional velocity, pressure and temperature, Pr is the Prandtl
number, defined by Pr=/ is the kinematic viscosiy and the thermal diffusivity). The non-
dimensional form results from scaling the lengths by the axial distance between the circular
disks (L) and the velocity by the energy diffusion velocity V
=
/L; the scales for time and
pressure are, respectively, L
/ and
/L
. The temperature, measured with respect to the
initial temperature
T
0
, is scaled by
()
T
:
T= (
TT
0
) /
()T
(3)
The initial conditions are:
t=0: V
(z, r, =0, (z, r, = 0 (4)
i.e. the liquid is motionless and at ambient temperature.
For t > 0, the boundary conditions are the non-slip conditions and the condition of prescribed
temperatures on the circular disks, the kinematic condition of stream surface (zero normal
velocity), the Marangoni conditions (shear stress balance) and the adiabatic condition on the
cylindrical interface:
57
on the cold disk 0
r
RL/ ;
02
V (z=0, r, , t) =
0; (z=0, r, , t) = - 1/2 (5)
on the hot disk
0
r
RL/ ;
02
V (z=1, r, , t) =
0; (z=1, r, , t) = 1/2 (6)
on the cylindrical free surface
01
z
; 02
V
r
(z, r=R/L, , t) =
0
(7a)
V
r
zr R L t Ma
T
z
zr R L t
Z
(, /,,) (, /,,) (7b)
r
V
r
zr R L t V zr R L t Ma
T
zr R L t
(, /,,) (, /,,) (, /,,)
(7c)
T
r
zr R L t(, / , ,)0 (7d)
where the reference Marangoni number Ma is defined as Ma=
T
()
T
L/.
3. NUMERICAL SOLUTION
3.1 Solution method
The equations (2a-c) and the initial and boundary conditions (4-7) were solved numerically in
cylindrical co-ordinates in primitive variables by a control volume method. The domain was
discretized with a non uniform but structured axisymmetric mesh and the flow field variables
defined over a staggered grid.
The axial velocity component V
z
is staggered in axial direction with respect to the point in
which temperature and pressure are computed. In a similar way the radial and azimuthal velocity
components are staggered in radial and azimuthal directions respectively (see Fig. 2).
The finite volume approach relies directly on the application of the integral form of balance
laws. Thus the conservation laws have been written for an arbitrary spatial domain bounded
by a surface . Since the collocation of the variables on the grid is staggered, each variable is
characterized by a different control volume;
