CC BY
31
Вычислительные технологии
Том 4, № 3, 1999
ON THE STABILITY OF AXISYMMETRIC THERMOCAPILLARY MOTIONS
V. K. Andreev Institute of Computational Modeling SB RAS Krasnoyarsk, Russia e-mail: [email protected]
Исследуются деформации цилиндрических слоев, вызванные термокапиллярными силами, в условиях, близких к невесомости. Соответствующие решения описывают термокапиллярные течения около критических точек и возникающие пространственные стационарные решения. Численно исследуется неустойчивость конвективных стационарных течений в цилиндрических слоях. Моделирование плавающей зоны в виде цилиндрического слоя позволяет исследовать не только конечную фазу, когда зона целиком расплавлена, но и промежуточные стадии. Одним из важных результатов данного исследования является тот факт, что увеличение радиуса твердой фазы слоя влечет расширение области неустойчивости в направлении коротких волн.
The investigation of thermocapillary flows has attracted much attention recently, particularly in the light of results of space experiments. Such surface-driven motions occur in the floating-zone process of crystals growth from the melt. The problems of flow stability of melts, arising here, are of undoubted interest both from practical and theoretical viewpoints.
It has been known since the time of Rayleigh [4] that circular static inviscid jets are unstable with respect to axisymmetric surface waves. A similar instability occurs in viscous jets [1]. At the same time the experimental data are available [2] which prove the existence of stable thermocapillary flows in molten zones.
In [6] it was shown that for a cylindrical float-zone with a unndeformable free surface, an increase of the heat-emission into the environment or an increase of Prandtl number results in stabilization of thermocapillary instability. An analysis of the possibility that the Rayleigh instability is stabilized by the thermocapillary effect was made in [7] for the case of liquid cylinder.
In this paper we considered the case when the fluid flows from the hot wall to the cold one and returns axially. If we neglect influence of the walls this flow can be approximately simulated as the motion in an infinite liquid cylindrical layer. The main problem of investigation then becomes the study of interaction between the capillary, thermocapillary and hydrodynamic instability mechanism and consideration of possibility to suppress Rayleigh instability. The simulation of a float-zone in the form of a cylindrical layer allows us to study not only the finite phase when the whole zone is melted but also the intermediate stages. As found in [5] for the case on non-deformable free surface with small Prandtl number the least stable state exists when the zone is not completely molten.
© V.K. Andreev, 1999.
1. The axisymmetric self-similar flows
Consider a stationary motion in a cylindrical layer. Assume, that on the inner solid body of the prescribed radius r0 a constant linear distribution of the temperature along z — axis is maintained, on the free surface r = ri the heat flux is absent. The length of the inner body and of the cylindrical layer are large enough; due to this we can neglect the boundary effects.
On the solid surface the condition of impenetrability is satisfied while on the free surface the conditions of impenetrability, balance of tangential and viscous stress are satisfied. We also assume that the free surface is undeformable. The reference point in the cylindrical coordinate system is on the solid surface at the point, where the temperature value 90 = const is achieved.
The axisymmetric self-similar solution has been sought in the form
u = -v^(C)/riC, w = nvt'(£)/r1£
9 = 90 + Ar^a(£), A = const,
P = Po - 1 p(v/r1)2[/(£) + An2],
a = a0 + k1(9 — 90)2, k1 = const, (1.1)
where u,w are components of the velocity vector on the axis r and z,9 is the temperature, p is the pressure, n^v/^r1 is a stream function, v,p,k1 are the kinematic viscosity, density and surface tension coefficient, respectively, 90,p0 are the temperature and pressure at the reference point on the solid surface, corresponding to the extreme value of the surface tension, A is unknown constant which must be determined in the course of solution, £ = r/r1,n = z/r1 are dimensionless variables.
The equations of motion and energy lead to the boundary problem for the system of nonlinear ordinary differential equations
/ = t2/£2 + 2f'/£,
t'' = (t' — ft' + t'2)/£ + (tt — t')/£2 — A£,
a'' = [Pr(at' — ta') — a']/£, d<£< 1, (1.2)
t = 0, t'' — t' = Ma2, a' = 0 on £ =1, (1.3)
t = t' = 0, a =1 on £ = d = r0/r1 (1.4)
which has been solved numerically. Here Pr = v/x is Prandtl number and M = k1A2h0/pv2 is modified Marangoni number, x is thermal conductivity. The result of numerical integration of system (1.2)-(1.4) are presented in Fig. 1 for Pr = 1. Parameter A is represented as a function of Marangoni number. We can distinguish four groups of solutions according to the temperature distribution inside the layer.
While solving a dynamic problem two different situation arise. The first one is realized when M > 0, i. e. temperature coefficient of the surface tension is more than zero and, besides, the sign of the Marangoni number doesn’t depend on heating or cooling of the fluid. The second one arose when M < 0. In the first situation (groups 1, 2) the fluid on the free surface flows away from the center n = 0. In the second one (groups 3, 4) the fluid towards to the center. It occurs due to the fact that the fluid on a free surface moves in the direction of the greater value of the surface tension coefficient. With A < 658 the flow (1) always changes the direction on the
Fig.1. The dependence of A on M for Pr = 1.
depth which is equal to 2/3 of the cylindrical layer thickness for different values of Marangoni number, but for the solutions of the group 2 with large M = 7 • 105 the turn point “moves” to the free surface. If A > 658 and M > 0 the second vortex occurs as shown in Fig. 2. With the increase of Marangoni number the secondary vortex domain increases, and with the further increase of M the intensity of the lower vortex begins to prevail the intensity of the upper one.
Fig.2. Streamlines for M = 2.2 • 104, A = 102, Pr = 0.005.
2. Liquid layers with thermocapillarity
Suppose that the constant temperature gradient d9/dz = —A (A > 0 is a specified constant) is imposed along the axis of the free surface. And the dependence of the surface tension coefficient on the temperature is described by the formula a = a0 — k(9 — 90). The length, time, velocity, temperature and pressure are referred to the scales r1, (prf /a0)1/2, (a0/pr1 )1/2, Ar1 and a0/r1.
Then the stationary axisymmetric thermocapillary motion, which occurs due to the change of the surface tension is described as follows:
w
= Re We-1/2
B (£2 — d2) + B2ln (£
pn = 2ReWe,
Re Pr
9 = —n+-— [B1 (£4 —1)/4—(B1 d2 + B2 ln d+B2 )(£2 — 1)+B2(£2 + d2 )ln £ + B1d4ln £ (2.1)
where
1 — d2 + 2ln d
B1
(1 — d2)(3 — d2) + 4ln d’ 1 — d2)2
2=
(1 — d2)(3 — d2) + 4ln d‘
Here Re = rfkA/pv2, We = a0r1/pv2 are Reynolds and Weber numbers respectively. Note that here the Reynolds number Re characterizes the relation of thermocapillary forces to viscosity forces and is modified Marangoni number (the relationship M = PrRe is satisfied).
With d ^ 0 the motion (2.1) reduces to the thermocapillary flow of a liquid jet, the stability of which was investigated in [7].
3. The linearized disturbance equations
We proceed to the investigation of the stability of eq. (2.1) and find the disturbed values of the velocity vector, pressure, temperature and normal component of a free surface in the form
(U, V, W, P,T, R) = (U (£), V (£), W (£), P (£),
T(£), R) exp [ian + im^ — iCr],
where a,m are axial and azimuthal wave numbers, C = Cr + iCi is a complex decrement and t is a dimensional time. In this form waves travel in a direction f = tan-1 (m/a) with respect to the positive n axis. The equations for perturbations take the form
aU + We1/2P? = — iaWe — ^ (£V)?,
aV + y We1/2P
2im
+ -jr
aW + We1/2w? U + iaWe1/2P = 1(£W?)?, (£U )j + imV + ia£W = 0, bT + We1/2 Pr 9? U — We1/2 Pr W = 1(£T?)?,
iWe1/2(aw — C) + m +
2
a
2
b = iWe1/2Pr (aw — C) + m + a2. (3.1)
£ 2
The conditions on the solid boundary £ = d are
U = V = W = T = 0. (3.2)
The conditions on the free surface £ =1 are
V? — V + imU = — imRe We-1/2(T + 9? R), iaU + w?? R + W? = —iaReWe-1/2(T + 9? R),
a
— P + 2We-1/2(U — iaw? R) = —Re We-1/2(T + 9? R) + (1 — a2 — m2)R,
i(aw — C)R = U,
T + Bi T + (9gg + Bi 9^ )R + iaR = 0. (3.3)
Here Bi is Biot number, the length of the disturbance wave is limited by a ^ We/Re.
4. Special case
Note that with We ^ to when v ^ 0 the system (2.1) is formally reduced to an ideal fluid, and the velocity field becomes zero. Therefore, it is interesting to compare the numerical results obtained with the results of the stability investigation of the resting ideal fluid. For this, we
consider a cylindrical layer of an ideal fluid bounded by solid inner and free outer surfaces.
Instead of (3.1)-(3.3) we obtain the equations for the disturbances
im
—iCU + P? = 0, —iCV + —P = 0,
—iCW + iaP = 0, (£U )g + imV + iaW = 0,
U = 0, £ = d,
iCR = U, —P = (1 — a2 — m2)R, £ = 1. (4.1)
We can obtain from eq. (4.1) the expression for the complex decrement
C i f(1 2 2) w (a)K'(ad) — /m(ad)Km(a)i1/2 (42)
C = M(1 — “ — " ’ X I'„M)A„(a) — /„(a)*-,(ad) j ' (4'2)
It is obvious that for m = 0 C is always real and, hence, the rest state is neutrally stable with respect to azimuthal disturbances. But for axisymmetric disturbances (m = 0) the analysis of eq. (4.2) shows that C is complex with 0 < a < 1 and real with a > 1. Note that with d ^ 0 eq. (4.2) changes into the corresponding expression for the complex decrement of a liquid cylinder obtained by Rayleigh [4].
5. Results and discussion
Figure 3 shows Cj plots of a constructed numerically for d = 0 . 1, Re = 200, Pr = 1, We = 104 and Bi = 1. In this case three unstable modes exist. Following [5] we write Wc for the mode with Cj > 0 with a < 1 (curve 1), W + for the mode with Cj > 0 (curve 2) and W- for the
mode with Cr < 0 (curve 3).
As seen in Fig. 3 the mode Wc is unstable at a < 1.026, the mode W + is unstable at 0.623 < a < 1.82, and for the mode W- the instability occurs at 0.945 < a < 1.22. In Fig. 4 the neutral curves Re* of a are represented for the same parameter values. Here curve 2 (mode W +) takes minimum value Re* = 59 with a = 1.05 and curve 3 (mode W-) has a minimum Re* = 116 with a = 1.02. As seen from Fig. 4 with a < 1.1 the disturbance Wc (curve 1) is most critical, and with a > 1.1 the disturbance W +.
In the following we examine the influence of Biot and Prandtl numbers on the disturbance behavior. For this we fix Re = 200 and We = 104, d = 0.1. The decrease of Bi has insignificant
Ci
0.2
0.1
- / > V1
2 , 1 1
3
-0.1
Fig.3. Growth rate Ci, vs. axial wavenumber a for d = 0.1, Re = 200, Pr = 1, Bi = 1, We = 104.
Re*
600
400
200
-
v \ /y
0.5 1 1.5
Fig.4. Neutral curves for d = 0.1, Pr = 1, Bi = 1, We = 104.
influence on modes Wc and W-. The instability domain of the mode W + increases with that, especially in the longwave range and Ci > 0 with a < 1.9. The change of Prandtl number also has an essential influence on the mode W +. With the decrease of Pr the instability domain of the mode Wc increases insignificantly. Modes W + and W- are stabilized with that and there exists some Prandtl number with which these modes are stable for every a. The increase of Pr results in the essential expansion of instability domain for modes W + and W-. Thus, with Pr = 5 for W + the instability occurs with a < 2.1 and for Wc with a < 10. The mode W- does not change practically as compared with the case Pr =1. The conducted analysis showed that Wc is a capillary mode corresponding to Rayleigh instability mechanism caused by the availability of the free surface. The motion in the fluid results in the occurrence of hydrodynamic mode W-, which represents a wave quickly propagating on the surface along the negative z-axis. This instability mechanism is highly dependent on Reynolds number changes. The dependence on Prandtl number is displayed only for Pr < 1, for Pr > 1 the dependence is weak. The third mechanism which can cause instability is the heat mode W + represents a wave quickly propagating along the z-axis in a positive direction. This mode is caused by the presence of the temperature inhomogeneity in the fluid and the thermocapillary effect connected with it. The essential influence on this mechanism effects changes of the Prandtl and Biot numbers. As seen from Fig. 5 the increase of d has a stabilizing influence on hydrodynamic and heat modes.
Thus, with the given parameters values the mode W- is stable for every value of a, and for the mode W + the instability occurs only with 1.12 < a < 1.37. Fig. 5 shows the plots of the a* critical values of a for the mode Wc with Bi =1, curve 1 corresponds to the value Pr = 0.1, curve 2 to Pr =1. As shown in Fig. 5 with the increase of d an expansion of the instability domain in the direction of short waves occurs. For small Pr such expansion is insignificant (e. g. a* = 1.19 with d = 0.95 for Pr = 0.1) while for mean values of Pr the instability domain already increases essentially (e.g. a* = 1.9 with d = 0.95 for Pr = 1).
Consider non-axisymmetric disturbances (m = 1). Figure 6 shows Ci plots of a for all modes with d = 0.1, Re = 200, Pr =1, We = 104 and Bi = 1. We found that for every mode a critical value of a* exists, with a > a* this mode is stable. Thus, for the capillary mode a* = 0.24 the hydrodynamic mode a* = 2.07 and for the heat mode a* = 1.25. With a decrease of the dimensionless layer thickness the stabilization of azimuthal disturbances occurs. As shown in Fig. 7 constructed for a = 0.5 the capillary mode becomes stable for every of a. The stability
a
Fig.6. Growth rate Ci, vs. axial wavenumber a for non-axisymmetric mode (m = 1) for d = 0.1,
Re = 200, Pr = 1, Bi = 1, We = 104.
domain of thermocapillary mode occurs here with 0.47 < a < 0.67 and the hydrodynamic mode for a < 0.54. Thus, thermocapillary flow described by eqs (2.1) is always unstable with respect to long-wave disturbances which is observed in experiments [3].
Ci
0.08 0.04 0
-0.04 -0.08
a
Fig.7. Growth rate Ci, vs. axial wavenumber a for non-axisymmetric mode (m = 1) for d = 0.5,
Re = 200, Pr = 1, Bi = 1, We = 104.
References
[1] Chandrasekhar S. Hydrodynamic and Hydromagnetic Stability. Oxford University Press, Oxford, 1961.
[2] Chun Ch.-H., Wuest W. Experiments on the transition from the steady to the oscillatory marangoni-convection of a floating zone under reduced gravity effect. Acta Astronautica, 6, 1976, 1073-1082.
[3] Feuerbacher B., Hamacher H., Naumann R. J. Materials Sciences in Space. Springer Verlag, New York, 1989.
[4] Rayleigh, Lord. Proc. London Math. Soc., 4, 1879, 10.
[5] RyabitSKII E. A. Instability of thermocapillary motion in cylindrical layer. J. Appl. Mech. Tech. Phys., 4, 1989, 50-52 (in Russian).
[6] Xu J.-J., Davis S. H. Convective thermocapillary instabilities in liquid bridges. Phys. Fluids., 27, No. 5, 1984, 1102-1107.
[7] Xu J.-J., Davis S. H. Instability of capillary jets with thermocapillarity. J. Fluid Mech., 161, 1985, 1-25.
nocmyuuAo, e peda^um 30 umnn 1998 s.
