Copyright c 2007 Tech Science Press FDMP, vol.3, no.4, pp.317-327, 2007 Instabilities and Pattern Formation in Thermocapillary Liquid Pools U. Schoisswohl 1 and H. C. Kuhlmann 2 Abstract: The flow in thermocapillary liquid pools heated or cooled from above can exhibit various flow patterns depending on the thermal conditions and the geometrical constraints. This pattern formation and the respective physical mechanisms are studied numerically by means of a linear-stability analysis. We focus on the transition from the steady axisymmetric to a threedimensional flow. Keyword: linear stability, Marangoni effect, thermocapillary flow, pattern formation. 1 Introduction Fusion welding is of increasing importance in industrial applications. There is a large body of work on the overall flow field in weld pools (see e.g. DebRoy and David, 1995) in which the effects of material composition, surface-active agents, vaporization, and others on the flow has been considered. The influence of the beam diameter and the beam power has been considered by Kamotani and Ostrach (1994) and Limmaneevichitr and Kou (2000). The shape of the axisymmetric weld pool due to the movement of the work piece as well as the influence of plasma effects have been studied by Do-Quang (2004). Yet, relatively little is known about the pattern formation in weld pools during the welding process and its dependence on the heating mode, in particular, under conditions of weightlessness. In this paper we investigate the melt flow for zero-gravity conditions for the case in which the 1 e-mail:ulrich.schoisswohl@tuwien.ac.at, Corresponding author, Institute of Fluid Mechanics and Heat Transfer, Vienna University of Technology, Resselgasse 3, A-1040 Vienna, Austria 2 e-mail: h.kuhlmann@tuwien.ac.at, Institute of Fluid Mechanics and Heat Transfer, Vienna University of Technology, Resselgasse 3, A-1040 Vienna, Austria workpiece is at rest. This problem is relevant to space applications of this technology and it is a step towards those terrestrial applications in which thermocapillary forces dominate buoyancy. More specifically, we search for sufficient conditions for an axisymmetric basic steady flow to become unstable to a non-axisymmetric one and we analyze the underlying physical mechanisms driving the instability. The model to be described below has previously been studied by Wagner, Friedrich, and Narayanan (1994) and others, and it has similarities with the model of Kamotani, Ostrach, and Masud (2000). A scaling analysis of the boundary layer structure in the corresponding infinite system has been carried out by Pumir and Blumenfeld (1996). However, an accurate prediction and a systematic study of the stability boundaries of the basic axisymmetric flow has not yet been made. 2 Statement of the Problem We consider a cylindrical volume of liquid with a height d and a radius R bounded by solid walls from the side and from below and by a free surface from above (cf. fig.1). The aspect ratio is defined as Γ = R/d. In this study we shall consider only Γ = 1. The analysis is carried out for the limit of large mean surface tension. Hence, static and dynamic deformations of the free surface are absent. The bottom and the walls of the cylinder represent the liquid solid interface and are assumed to be at the constant melting-point temperature. The temperature distribution on the free surface is given by an imposed axisymmetric heat flux Q(r) with a maximum value Q(0) =Q max at the center of the free surface. To minimize the governing parameters we assume a parabolic profile for Q(r).
318 Copyright c 2007 Tech Science Press FDMP, vol.3, no.4, pp.317-327, 2007 f ree sur f ace solid wall d d/2 R z ϕ r solid bottom For the boundary conditions we assume no-slip andnopenetrationontherigidwallsandonthe bottom of the cylindrical domain (U = 0), and a constant wall temperature corresponding to the melting point, i.e. T = T melt = T 0 = 0. The thermal boundary condition on the free surface is determined by the heat flux z T = (1 r) 2. Note that we use the temperature scale ΔT = Q max /k,where k is the thermal conductivity of the liquid. The scale derives from the maximum heat flux which is prescribed at the free surface. Neglecting stresses in the ambient gas the stress balance on the flat non-deformable liquid gas interface is given by (see, e.g. Kuhlmann, 1999) Figure 1: Geometry and coordinate system. S n +Re (I nn) T = 0. (4) The non-uniform heat flux causes a non-uniform distribution of the temperature on the free surface which gives rise to a non-uniform surface tension σ(t ). Owing to the thermocapillary effect a surface stress σ drives a flow from the hotter near-axis region into the colder peripheral regions. This thermocapillary convection is also called Marangoni flow (see Scriven and Sternling, 1960; Mills, Keene, Brooks, and Shirali, 1998). The temperature dependence of the surface tension σ can be approximated by σ(t )=σ(t 0 ) γ(t T 0 )+O [ (T T 0 ) 2], (1) Here S is the viscous stress tensor, I the unity matrix, and n = e z the unity normal vector of the free surface. Proper boundary conditions on the axis complete the problem. The two dimensionless parameters are the thermocapillary Reynolds and Prandtl numbers Re = γδtd ρν 2 and Pr = ν κ, (5) where γ is the linear Taylor expansion coefficient γ = σ T. (2) Employing cylindrical coordinates (r,ϕ,z) the Navier-Stokes-equations are formulated in dimensionless form using the scalings d, d/u, γδt/ρν,andδt for length, time, velocity, and temperature to obtain ( t +U )U = P + 2 U, (3a) Pr( t +U )T = 2 T, (3b) U = 0. (3c) With U =(U,V,W) T, P, andt denoting the dimensionless velocity, pressure and temperature fields. where κ is the thermal diffusivity and ν the kinematic viscosity. 3 Numerical Solution Techniques 3.1 Basic state The basic state of the above problem features a steady axisymmetric toroidal vortex flow ( t 0, ϕ 0,V 0). The solution is calculated, using a formulation in primitive variables and a
Instabilities and Pattern Formation in Thermocapillary Liquid Pools 319 cylindrical coordinate system, i.e., U U r U W r +W U z = P r + 1 r r ( r U r ) + 2 U z 2 (6a) W +W z = P z + 1 ( r W ) + 2 W r r r z 2 U T T +W r z = 1 [ 1 Pr r 1 r (ru) + W r z = 0. r ( r T ) + 2 ] T r z 2 (6b) (6c) (6d) This system of partial differential equations is replaced by a system of difference equations by means of a discretization with finite volumes (Albensoeder, Kuhlmann, and Rath, 2001) on a staggered grid (Ferziger and Perić, 2002) with a resolution of 70 70 grid cells (see also Sec.3.4). The resulting system of difference equations is solved by means of a Newton-Raphson-Method (Deuflhard, 2004) and an efficient linear-systems solver from mathematical subroutine libraries (BLAS, GSL, LAPACK). 3.2 Linear stability analysis Once the 2D-axisymmetric basic state is calculated its stability is investigated by means of a linear-stability analysis. To that end we decompose the 3D non-axisymmetric flow state into the basic state and a perturbation (U,P,T) T =(U 0,P 0,T 0 ) T +(ũ, p, T ) T. (7) Substitution into the above equations and linearization with respect to the perturbation quantities yields t ũ +ũ U 0 +U 0 ũ = p + 2 ũ, (8a) Pr ( ) t T +U 0 T +ũ T 0 = 2 T, (8b) ũ = 0. (8c) A general solution of this system can be written as ũ û p (r,ϕ,z,t)= ˆp (r,z)e imϕ e λt +c.c., (9) T ˆT where m is an azimuthal wave number. Owing to the periodicity in ϕ-direction m is an integer. Note that the boundary conditions on the axis r = 0for the perturbations depend on the the wave number m (see, e.g. Xu and Davis, 1984). The complex eigenvalues λ = σ +iω depend on the wave number m. The real part σ = R(λ ) represents the temporal growth rate while the imaginary part ω = I(λ ) is an angular frequency. Employing this normal-mode representation of the solution and utilizing the same discretization as for the basic state we obtain the generalized eigenvalue problem A x = λ B x, (10) with eigenvector x and eigenvalue λ. A and B are the matrix representations of the linear difference equations and the boundary conditions. The generalized eigenvalue problem is solved using inverse iteration. In order to find a neutral stability boundary one has to find an eigenvalue λ with a real part σ = 0. To find this neutral boundary the roots of σ are searched for by means of the secant method. During this zero search the Reynolds number Re has to be varied, and the basic state calculation as well as the linear-stability analysis have to be carried out repeatedly. The search is performed until an eigenvalue λ with a real part σ = σ neut = 0 respectively a neutral Reynolds number Re neut is found. The critical Reynolds number Re c is defined as the minimum envelope of the neutral Reynolds numbers Re neut. 3.3 Energy analysis For a deeper insight into the physical mechanisms driving the transition process we analyze the transfer rates of kinetic and thermal energy between the basic state (U 0,P 0,T 0 ) T and the neutral mode (ũ, p, T ) T. The rate of change of kinetic energy Ė kin is governed by the Reynolds Orr equation E kin t = 1 2 t V ũ 2 dv = D +M r +M ϕ +I v. (11) It can be derived by scalar multiplication of the momentum equation (3a) by ũ andbyintegration
320 Copyright c 2007 Tech Science Press FDMP, vol.3, no.4, pp.317-327, 2007 over the volume V. The terms on the right-handside are the rate of viscous dissipation D = ( ũ) 2 dv (12) V and total energy production I v = ũ [(ũ )U 0 ]dv. (13) V Energy is produced by advection of basic state momentum U 0 by the perturbation flow ũ, thus adding to the perturbation flow itself. The quantities M r = ũ z ũds, M ϕ = ṽ z ṽds (14) S represent the work done by Marangoni forces acting on the free surface in the radial and azimuthal directions. In a similar way we define a thermal energy. Multiplying the temperature equation (3b) by T and integrating over the volume V we find E th t Here = 1 2 t V S T 2 dv = D T +H +I T. (15) D T = 1 ( T ) 2 dv (16) Pr V is the rate of heat diffusion, H = 1 1 Pr S 2 ( z T 2) ds (17) is a measure for the supply of thermal energy E th through the free surface, and I T = T (ũ )T 0 dv (18) V is the rate of production of thermal energy. It is generated by convection of basic-state temperature T 0 by the perturbation flow field ũ thus adding to the perturbation temperature field. For a further analysis we shall also need the local rates of change of energy, i.e. the densities of the rates of change of energy. These are the integrands of the above integrals and will be denoted, henceforth, by lower-case letters, e.g. i T = T (ũ )T 0. Table 1: Critical Reynolds numbers Re c as function of the grid resolution N r N z. The critical wave numbers for Pr = 0.02 and Pr = 4.0 are m = 3andm = 2, respectively. Re c N r N z Pr = 0.02 Pr = 4 30 30 39,556 119,731 40 40 36,647 117,413 50 50 35,528 114,571 60 60 34,992 112,453 70 70 34,709 110,363 80 80 34,547 109,969 90 90 34,488 109,283 100 100 34,386 108,818 3.4 Grid and grid convergence The calculations have been performed on a staggered grid with a resolution of 70 70 cells, i.e. N r = N z = 70. In order to resolve the developing boundary layers the grid is compressed towards the top and side walls, the stretching factor being 0.98 in both directions. The grid convergence is demonstrated in Table 1 for two representative Prandtl numbers. 4 Results 4.1 Basic state for low Prandtl numbers The temperature distribution on the free surface together with the thermocapillary effect drives a radial fluid motion from the center to the rim of the free surface. Owing to continuity a return flow is created below the free surface flow, resulting in a toroidal vortex. For a representative low Prandtl number of Pr = 0.02 this basic flow is shown in fig. 2a for a Reynolds number Re = Re c = 34,709 at which this flow becomes unstable. The corresponding basic temperature field, displayed in fig. 2b, is almost conducting at criticality, because the relevant critical Peclet number, i.e. the critical Marangoni number Ma c = Re c Pr = 694 is small for this low value of the Prandtl number. Here and in all cross sections shown in the following the direction of rotation of the basic vortex is clockwise
Instabilities and Pattern Formation in Thermocapillary Liquid Pools 321 0.5 0.3 0.1 z 0.1 0.5 0.3 0.1 z 0.1 0.3 0.5 0.0 0.2 0.4 0.6 0.8 1.0 (a) r 0.3 0.5 0.0 0.2 0.4 0.6 0.8 1.0 (a) r Pro 0.5 0.3 0.1 z 0.1 0.3 0.5 0.0 0.2 0.4 0.6 0.8 1.0 r (b) Figure 2: Basic state for Pr = 0.02 at critical conditions Re = Re c = 34,709. (a) Stream function ψ 0 with isolines at Δψ 0 = 2.046. (b) Temperature field with isotherms at ΔT 0 = 0.0263. For the relatively high Prandtl number of Pr = 4 and at the critical Reynolds number Re c = 110,362 the structure of the toroidal vortex (fig. 3a) differs little from the one at low Prandtl number Pr = 0.02 (fig. 2a), even though the critical Reynolds number Re c is more than three times as high. However, the temperature field given in fig. 3b differs significantly from the low-prandtlnumber case. This is due to the high critical Marangoni number Ma c = Re c Pr = 441,448 indicating a strong convective effect on the temperaof 0.3 0.1 z 0.1 0.3 0.5 0.0 0.2 0.4 0.6 0.8 1.0 r (b) Figure 3: Basic state for Pr = 4 at critical conditions Re = Re c = 110,362. (a) Stream function roo0.5 ψ 0 with isolines at Δψ 0 = 0.588. (b) Temperature field with isotherms at ΔT 0 = 0.0013. (ψ 0 < 0). 4.2 Basic state for high Prandtl numbers ture field. In fact, the isotherms significantly deviate from the conducting state and thermal boundary layers are about to be established. 4.3 Stability boundary Neutral Reynolds numbers Re neut have been calculated for Prandtl numbers in the range 10 10 Pr 8 and for wave numbers m = 1,2, and 3. The dependence of the stability boundary on Pr and m is shown in fig. 4 in form of neutral curves Re neut (Pr,m). The relatively high values of the critical Reynolds numbers are due to the temperature scale used. If we would use the actual temperature difference as the temperature scale crit-
322 Copyright c 2007 Tech Science Press FDMP, vol.3, no.4, pp.317-327, 2007 ical Reynolds numbers of the order of O(10 3 ) would result. However, the actual temperature difference is not known a priori, because it is part of the solution and can thus be determined only a posteriori. Two ranges can clearly be distinguished: in the low-prandtl-number range (Pr 1) the 2D axisymmetric steady flow is unstable to a 3D nonaxisymmetric steady flow. The critical wave numbers for Γ = 1 are either m c = 2orm c = 3 depending on the Prandtl number. In the high-prandtlnumber range (Pr 1) the basic flow for Γ = 1 is unstable to a 3D non-axisymmetric oscillatory flow with a critical wave number m c = 2 throughout. neutral Reynolds numbers 1e+06 1e+05 m1 m2 m3 3d stationary state 2d axisymmetric state 0.01 0.1 1 Prandtl number 3d oscillatory state Figure 4: Neutral curves for wave numbers m = 1,2, and 3. The type of the neutral mode (steady or time-dependent) is indicated. 4.4 Low-Prandtl-number instability mechanism Figure 5 shows the critical velocity field ũ and the critical perturbation temperature T on the free surface at z = 0.5forPr= 0.02. We find three strong perturbation temperature maxima and three corresponding minima. In addition, three weak maxima and minima arise. The free surface perturbation flow between adjacent extrema is either from hot to cold or from cold to hot. Most significant, of course, are the strong perturbation temperature extrema and the flow connecting them. For the strong spots, the perturbation flow is always directed from cold to hot. Such a motion, however, Figure 5: Perturbation flow (arrows, interpolated from the numerical data) and temperature field (lines) on the free surface at z = 0.5forPr= 0.02. Negative values are indicated by gray lines. The parameters are m c = 3andRe c = 34,709. cannot be created by the normal Marangoni effect, since γ > 0. In fact, the Marangoni stresses produced by the strong surface perturbation temperature spots are counteracting the surface perturbation flow of the critical mode. Hence, the Marangoni effect is acting stabilizing in this situation. For that reason the perturbation velocity field must be driven by a different mechanism. For small Prandtl numbers at which heat diffusion dominates heat convection such a mechanism should be inertial. For the inertial instability of the axisymmetric toroidal thermocapillary vortex flow in low- Prandtl-number liquid bridges Nienhüser and Kuhlmann (2002) have shown that vortex straining as well as centrifugal effects may contribute to an inertial destabilization of the basic flow (for the lid-driven cavity, see Albensoeder, Kuhlmann, and Rath, 2001). In the following we shall argue that centrifugal mechanisms are dominating for the present low-prandtl-number instability for Γ = 1. Bayly (1988) derived a generalized Rayleigh cri-
Instabilities and Pattern Formation in Thermocapillary Liquid Pools 323 terion which states that the flow of an inviscid fluid is subject to a centrifugal-type instability if it consists of closed convex streamlines, and if the magnitude of circulation decreases outwards. This criterion has been reformulated by Sipp and Jacquin (2000) as follows. A two-dimensional inviscid flow is centrifugally unstable if and along the solid sidewall until the basic flow separates from the sidewall to form the radial return flow still in the upper half of the cylindrical domain to the central region where the return flow approaches the free surface again. Φ(r) := U 0 ω 0 < 0 (19) R all along a closed convex streamline. Here ω 0 is the vorticity of the basic flow and R is the local radius of curvature of a streamline which can be calculated as (see Sipp and Jacquin, 2000) R = U 0 3 ( ψ 0 ) (U 0 U 0 ). (20) Figure 7: Vertical cut along the axis of the cylinder showing the total local kinetic-energyproduction rate i v (integrand of I v ) restricted to positive values and the regions for which Φ(r) < 0 holds (gray-shading). The cut is shown at an azimuthal angle for which the maximum local kinetic energy production rate takes its absolute maximum. The parameters are Pr = 0.02, m c = 3, and Re c = 34,709. Note that i v has twice the azimuthal period of the critical mode. Figure 6: Vertical cut along the axis of the cylinder showing the basic-state stream function ψ 0 (solid lines) and the regions for which Φ(r) < 0 holds (gray-shading). The data are evaluated at the critical conditions for Pr = 0.02 (m c = 3, Re c = 34,709). Even though the criterion provided by Sipp and Jacquin (2000) is valid for inviscid flows only, we have evaluated (19) for the present viscous basic flow. The result is shown in fig. 6 for the Prandtl number Pr = 0.02, the critical Reynolds number Re c = 34,709, and the critical wave number m c = 3. The criterion (19) holds true in the gray-shaded areas of fig. 6. Most notably, the regions which would favor a centrifugal instability in an inviscid flow are aligned with the outer streamlines of the toroidal vortex. The region extends from the cold corner from which the accelerated free-surface flow is deflected into the bulk A comparison of the regions favoring a centrifugal instability with the total local production rate of kinetic energy i v is provided in (fig. 7). In fact, the regions in which a significant amount of kinetic energy is produced lie well within the regions in which (19) is satisfied. I.e., most of the kinetic energy transfer from the basic flow to the perturbation mode takes place in a region that would be subject to a centrifugal-type instability if the flow were inviscid. Hence we conclude, that the low-prandtl-number flow in a cylindrical weld pool of aspect ratio one and a parabolic heat flux is unstable to a centrifugal-type instability. This behavior is very similar to the centrifugal instabilities in lid-driven cavities (Albensoeder, Kuhlmann, and Rath, 2001) and the Taylor Görtler instability of the boundary layer flow along convex walls Drazin and Reid (1981).
324 Copyright c 2007 Tech Science Press FDMP, vol.3, no.4, pp.317-327, 2007 4.5 High-Prandtl-number instability mechanism As a representative case for high Prandtl numbers we consider Pr= 4. In this case the critical wave number is m c = 2 and the corresponding critical Reynolds number is Re c = 110,362. The perturbation temperature on the free surface at z = 0.5 (fig. 8) exhibits four temperature extrema, two maxima and two minima. Since the perturbation flow field is essentially directed from the hot to the cold perturbation temperature spots, thermocapillary forces contribute to the driving of the perturbation flow. In fact, they may be the only effective driving force of the perturbation flow if no other inertial mechanisms are operating. The question then arises: how are the surface temperature extrema created? Since the rate of diffusion of perturbation temperature (16) is much smaller for high Prandtl numbers than for low ones, the surface spots could possibly be created by the vertical component of the perturbation flow which must arise due to continuity (similar as in the classical Marangoni problem, see Pearson, 1958). Such a mechanism, however, cannothold, because the vertical temperature gradient has the wrong sign: the free surface is hotter than the fluid below it (cf. fig. 3b). Figure 9: (color online). Perturbation flow (arrows), perturbation temperature field (color) and local thermal energy production rate i T (lines) at a vertical cut taken at an azimuthal angle for which the temperature perturbation takes its absolute maximum, for parameters Pr = 4, Re = Re c = 110,362 and m = m c = 2. Figure 8: Perturbation flow (arrows, interpolated from the numerical data) and perturbation temperature field (lines) on the free surface for z = 0.5 and Pr = 4. Negative values are indicated by gray lines. The pattern rigidly propagates in clockwise direction. The parameters are m c = 2and Re c = 110,362. The only other possibility that remains is heat conduction from even stronger temperature perturbation extrema in the bulk. In fact, the maximum of the perturbation temperature arises in the bulk, as can be seen from fig. 9 showing a vertical cut along the axis for an azimuthal angle for which the perturbation temperature takes its absolute maximum. The azimuthal angle is indicated by the solid line in fig.11. Again the question arises how these high and low perturbation temperature spots in the bulk are created. The strong perturbation-temperature spots in the bulk are generated by thermal energy production. The extrema of the local thermal energy production rate i T (the integrand of I T ) are located in close vicinity of those of the perturbation temperature. This is shown in fig. 10, which demonstrates how the thermal energy is extracted from the basic state by convection of basic-state temperature due to the perturbation velocity field. The vertical plane through the axis of the cylinder in which the maximum energy production arises is indicated as a dotted line in fig. 11. The
Instabilities and Pattern Formation in Thermocapillary Liquid Pools 325 Figure 10: (color online). Perturbation flow (arrows), local thermal energy production rate i T (lines) and basic state temperature field (color) shown at a vertical cut taken at an azimuthal angle for which the maximum local thermal energy production rate takes an absolute maximum, for parameters Pr = 4, Re = Re c = 110,362 and m = m c = 2. The color scale is a range from blue to red, corresponding to a range cold to hot. Figure 11: (color online). Perturbation flow (arrows), perturbation temperature field (color) and local thermal energy production rate i T (lines) at the midplane for z = 0.0, for parameters Pr = 4, Re = Re c = 110,362 and m = m c = 2. Red (blue) colors indicate perturbation temperature maxima (minima). The straight solid line gives the location of the cut in fig.9, and the straight dotted line the location of the cut given in fig.10. same mechanism applies to the orthogonal vertical plane parallel to the axis in which the flow direction and the temperature perturbations are inverted. From fig. 11, which shows the fields in the midplane z = 0 as viewed from above, we find that the production extrema arise slightly ahead in clockwise direction of the temperature extrema. This is an indicator for the clockwise rotation of the pattern and consistent with the negative phase velocity which, for m > 0 and together with (9) is determined by the positive critical angular frequency ω c = 54.54 for the case presented. Of course, the critical modes arise as pairs with ω = ±ω c. 5 Discussion and Conclusion The linear stability of the axisymmetric thermocapillary flow in a model weld pool has been investigated for Prandtl numbers in the range 10 10 Pr 8. We found two different types of instabilities. For Pr 1 the instability is stationary and the critical flow is three-dimensional and steady. For Pr 2 the instability is timedependent leading to three-dimensional waves. We did not investigate the intermediate region, because both branches of the critical curve increase to very high Reynolds and Marangoni numbers such that the numerical calculations become increasingly demanding. Most likely, however, both critical curves will intersect at some point in the interval Pr [1;2]. For low Prandtl numbers, the temperature field is nearly conducting and it practically decouples from the stability problem. The basic temperature field simply serves to drive the toroidal ring vortex. The latter becomes unstable due to centrifugal effects arising in the region of curved streamlines which originate from the cold corner. The mechanism has been demonstrated for the specific case of Pr = 0.02. The increase of the critical curve with the Prandtl number can be attributed to a thermocapillary effect: The critical mode creates surface temperature perturbations which result in Marangoni stresses counteracting the flow that produces the surface-temperature perturbations.
326 Copyright c 2007 Tech Science Press FDMP, vol.3, no.4, pp.317-327, 2007 For Pr 2 the instability is oscillatory. The perturbation flow extracts thermal energy from the basic temperature field leading to strong temperature extrema in the bulk. By conduction the strong bulk extrema create weak temperature perturbations on the free surface which drive a thermocapillary perturbation flow. The corresponding return flow reinforces the bulk flow which creates the internal extrema. The basic flow structures and also the instability mechanism found for low and high Prandtl numbers are very similar to those in thermocapillary liquid bridges (half-zone model) (Wanschura, 1996; Nienhüser and Kuhlmann, 2002). We conclude that the present low- and high-prandtlnumber instabilities are due to centrifugal and straining effects and hydrothermal waves, respectively, just as in the half-zone model. Acknowledgement: Computing time provided by the Zentraler Informatik Dienst (ZID) of the Vienna University of Technology is gratefully acknowledged. References Albensoeder, S.; Kuhlmann, H. C.; Rath, H. J. (2001): Three-dimensional centrifugal-flow instabilities in the lid-driven cavity problem. Phys. Fluids, vol. 13, pp. 121 135. Bayly, B. J. (1988): Three-dimensional centrifugal-type instabilities in inviscid twodimensional flows. Phys. Fluids, vol. 31, pp. 56 64. DebRoy, T.; David, S. A. (1995): Physical processes in fusion welding. Rev. Mod. Phys., vol. 67, pp. 85 112. Deuflhard, P. (2004): Newton Methods for Nonlinear Problems. Affine Invariance and Adaptive Algorithms, volume 35 of Series Computational Mathematics. Springer, Heidelberg, Berlin. Do-Quang, M. (2004): Melt convection in welding and crystal growth. PhD thesis, KTH Stockholm, 2004. Drazin, P. G.; Reid, W. H. (1981): Hydrodynamic Stability. Cambridge University Press, Cambridge. Ferziger, J. H.; Perić, M. (2002): Computational Methods for Fluid Dynamics. Springer, Berlin. Kamotani, Y.; Ostrach, S. (1994): Analysis of velocity data taken in surface tension driven convection experiment in microgravity. Phys. Fluids, vol. 6, pp. 3601 3609. Kamotani, Y.; Ostrach, S.; Masud, J. (2000): Microgravity experiments and analysis of oscillatory thermocapillary flows in cylindrical containers. J. Fluid Mech., vol. 410, pp. 211 233. Kuhlmann, H. C. (1999): Thermocapillary Convection in Models of Crystal Growth, volume 152 of Springer Tracts in Modern Physics. Springer, Berlin, Heidelberg. Leypoldt, J. (1999): Dreidimensionale numerische Simulation thermokapillarer Strömungen in zylindrischen Flüssigkeitsbrücken. PhD thesis, Universität Bremen, 1999. Limmaneevichitr, C.; Kou, S. (2000): Visualization of Marangoni convection in simulated weld pools. Welding J. (Res. Suppl.), vol. 79, pp. 126 135. Mills, K. C.; Keene, B. J.; Brooks, R. F.; Shirali, A. (1998): Marangoni effects in welding. Phil. Trans. Roy. Soc. London A, vol. 356, pp. 911 925. Nienhüser, C.; Kuhlmann, H. C. (2002): Stability of thermocapillary flows in non-cylindrical liquid bridges. J. Fluid Mech., vol. 458, pp. 35 73. Pearson, J. R. A. (1958): induced by surface tension. 4, pp. 489 500. On convection cells J. Fluid Mech., vol. Pumir, A.; Blumenfeld, L. (1996): Heat transport in a liquid layer locally heated on its free surface. Phys. Rev. E, vol. 54, pp. R4528 R4531.
Instabilities and Pattern Formation in Thermocapillary Liquid Pools 327 Scriven, L. E.; Sternling, C. V. (1960): The Marangoni effects. Nature, vol. 187, pp. 186 188. Shtern, V.; Hussain, F. (1993): Azimuthal instability of divergent flows. J. Fluid Mech., vol. 256, pp. 535 560. Sipp, D.; Jacquin, L. (2000): Threedimensional centrifugal-type instabilities of twodimensional flows in rotating systems. Phys. Fluids, vol. 12, pp. 1740 1748. Smith, M. K. (1986): Instability mechanisms in dynamic thermocapillary liquid layers. Phys. Fluids, vol. 29, pp. 3182 3186. Wagner, C.; Friedrich, R.; Narayanan, R. (1994): Comments on the numerical investigation of Rayleigh and Marangoni convection in a vertical circular cylinder. Phys. Fluids, vol.6, pp. 1425 1433. Wanschura, M. (1996): Lineare Instabilitäten kapillarer und natürlicher Konvektion in zylindrischen Flüssigkeitsbrücken. PhD thesis, ZARM Universität Bremen, 1996. Xu, J.-J.; Davis, S. H. (1984): Convective thermocapillary instabilities in liquid bridges. Phys. Fluids, vol. 27, pp. 1102.