Numerical Investigation of Combined Buoyancy and Surface Tension Driven Convection in an AxiSymmetric Cylindrical Annulus

Abstract. A numerical study is conducted to understand the effect of su rface tension on buoyancy driven convection in a vertical cylindrical ann ular cavity filled with a low Prandtl number fluid. The inner and outer cylinders are maint ained at different uniform temperatures and the horizontal top and bottom walls are the rmally insulated. The upper free surface is assumed to remain flat and non-deformab le. A finite difference scheme consisting of the Alternating Direction Implicit me thod and the Successive Line Over Relaxation method is used to solve the vorticity stream function formulation of the problem. Detailed numerical results of heat transfer ra te, temperature and velocity fields have been presented for a wide range of physical parame ters of the problem. The flow pattern and temperature distribution in the annular cav ity are presented by means of contour plots of streamlines and isotherms. The rate of he at transfer is estimated by evaluating the average Nusselt number. Further, the pres ent numerical results are compared with the existing results and are found to be in good agreement.


Introduction
It is well known that, when a free liquid surface is present, variations in the liquid surface tension at the free surface due to temperature gradients, can induce motion within the fluid called thermocapillary flows or Marangoni convection.Flows, induced by surface tension gradients, occur, in the processing of materials where small amounts of material are melted, and allowed to re-solidify.In these processes, large temperature differences are created in small amounts of liquid, which may have a free surface.The surface tension gradients, resulting from the large temperature differences in the liquid, can produce fluid velocities as large as one-meter per second.The resulting flow, affects the solute distributions, and crystal structure in the re-solidified material.Since thermocapillary convection is a cause of defects in the crystal growth process, an understanding of the heat transfer, and fluid motion, in the process, is needed, in order to grow defect-free crystals.
The study of fluid motion, and the associated transport processes, generated, either by buoyancy or by surface tension in cavities, have been independently examined by many investigators.Carpenter and Homsy [1] investigated steady thermocapillary flow in a two dimensional square cavity for differentially heated side walls using the finite difference method and concluded that thermocapillary flows are extremely sensitive to Prandtl number.The numerical simulation of thermocapillary convection, in differentially heated rectangular enclosures, have been reported by Strani et al. [2], Zebib et al. [3], Srinivasan and Basu [4], Bergman [5] and Rivas [6].Vrentas et al. [7] have studied surface tension driven convection in a cylindrical geometry using linear stability analysis and perturbation techniques.They showed that, the structure of surface tension driven flow, and heat transfer rate, is significantly modified at super critical Marangoni numbers.An experimental study of oscillatory thermocapillary convection in cylindrical containers was performed by Kamotani et al. [8].From the results, they observed that, beyond a certain temperature difference between the container wall and the heating wire, oscillatory thermocapillary flow patterns appear.Later Kanouff and Greif [9] studied thermocapillary flow in a differentially heated two dimensional square cavity using the control volume technique.They found that, for low Marangoni numbers, the flow was steady, stable, and symmetric about the centre vertical plane, and for Marangoni numbers above 1200, the flow was found to be, given an initial disturbance, oscillatory.
The combined effect of surface tension and buoyancy, on a fluid layer being heated from below, was first studied by Nield [10].He showed, by linear stability analysis, that the destabilizing effects of surface tension gradient and buoyancy, reinforce each other, and the coupling between these two effects is very tight.Vrentas et al. [11] have studied surface tension and buoyancy driven convection in a bounded cylindrical geometry, for a wide range of aspect ratios, to determine the critical Marangoni and Rayleigh numbers.Numerical results are presented, by Chen et al. [12], for steady natural convection in a two dimensional rectangular enclosure with the top wall modelled as an impermeable rigid or free-moving boundary.Their results reveal that, the free surface flow has larger circulation strength, and higher heat transfer rate, relative to the rigid surface case.The influence of Marangoni and buoyancy forces, during solid-liquid phase change in two dimensional melting of a solid slab has been reported by Bergman and Webb [13].
Most of the earlier studies on combined surface tension and buoyancy driven convection are concentrated either on the rectangular cavity or cylindrical geometry.Attention has not been given, to the study of combined buoyancy and surface tension driven convection, in a vertical cylindrical annulus, in spite of its important applications in material processing.The objective of the present numerical study is, to investigate the flow pattern, and heat transfer, of nonlinear convection, driven, by the combined mechanism of buoyancy and surface tension forces, in a cylindrical annulus, filled with a low Prandtl number fluid.

Mathematical formulation
Consider a cylindrical annulus, as shown in Fig. 1, formed by two co-axial cylinders of inner and outer radii r i and r o respectively.The top free surface and bottom wall are thermally insulated, and the side walls are maintained at constant but different temperatures.The cavity is filled with a low Prandtl number fluid, which is assumed to have constant physical properties, and obeys a Boussinesq approximation, according to which, its density is taken as constant, except, in the buoyancy term of the axial momentum equation.Further, we assume that the flow is axi-symmetric, and the surface tension on the upper boundary, is assumed to vary linearly with temperature: where γ = (−1/σ 0 )∂σ/∂θ and the subscript 0 refers to a reference state.The nondimensional equations governing the conservation of mass, momentum and energy for an unsteady laminar flow, after neglecting viscous and ohmic dissipation, in cylindrical coordinates (r, x) are The vorticity ζ is of the form Since, the flow properties depend only on two spatial coordinates, a vorticity-stream function approach, has been chosen.It is well known that the vorticity stream function formulation ensures conservation of volume.The non-dimensional variables are where (u, w) are the dimensional velocity components in (r, x) directions respectively, ζ * and ψ are the dimensional vorticity and stream function, θ is the dimensional temperature, ν is the kinematic viscosity, ρ is the fluid density, α is the thermal diffusivity and β is the coefficient of volumetric expansion.The non-dimensional parameters appearing in the above equations are: να , the Rayleigh number, P r = ν α , the Prandtl number, A = L D , the aspect ratio and λ = ro ri , the radii ratio.The coupled nonlinear equations ( 1)-( 4), namely energy, vorticity, stream function equations and velocity-stream function relation, are solved numerically subject to the following non-dimensional initial and boundary conditions: The boundary conditions on Ψ are obtained from those on the velocity components U and W , such as no through-flow and no-slip conditions.Normally, Ψ on boundary is an arbitrary constant which can be set to zero.The dynamical condition on the upper free surface represents the balance between shear stress and surface tension gradient and is responsible for the establishment of thermocapillary flow in the cavity.The nondimensional parameter is the Marangoni number.With these boundary conditions the system will tend to a steady state having constant properties.
The overall heat transfer rate across the enclosure is expressed by the average Nusselt number at the hot inner cylindrical wall and is defined as where, N u is the local Nusselt number for the hot inner wall which is defined by Nu = − ∂T ∂R R=0 .

Numerical procedure and validation
A numerical technique, based on the two-step ADI (Alternating Direction Implicit) method, has been used to solve the vorticity transport and energy equations.In order to improve stability of the numerical scheme, and to speed up convergence, the non-linear convective terms in the ADI method are approximated by second upwind differences.
The diffusion terms are approximated using central differences.From the known values of U , W , T , ζ and Ψ at time τ = 0 , we determine T and ζ at time τ = △τ using ADI method.The stream function is then computed using newly computed vorticities ζ.The Successive Line Over Relaxation (SLOR) procedure is used to solve the system of equations arising from the stream function equations.This method is preferred in the present study since it converges in fewer iterations than the usual point iteration methods.
The discretized algebraic equations form a tridiagonal matrix which can be solved by the Thomas algorithm.Then the values of U and W at the time τ = △τ are computed using central difference approximation to equation ( 4).For solution of the vorticity equation, vorticity at the boundary, ζ w , is obtained by expanding the stream function in Taylor series [14,Roache].The form used in this study for ζ w is where w denotes the boundary node and ∆η is the spatial interval in the direction normal to the boundary.However, at the top free surface, the vorticity boundary condition is Several grid sizes were used to examine grid dependency.The uniform meshes employed in the present study were: 41 × 41, 61 × 61, and 81 × 81.Because of minor differences between the 61 × 61 and 81 × 81 grids, and for the consumption of computer time, a 61 × 61 grid for A = 1 was used for further calculations.Similar grid dependency studies were carried out for A = 2, and an optimum grid size was obtained.The above computational cycle was then repeated for each of the next time levels and the steady state solution is obtained when the following convergence criterion Here ǫ is a pre-specified constant usually set at 10 −5 .
In the present formulation, the equations are made nondimensional such that the present formulation in cylindrical coordinates can be reduced to rectangular coordinates.Hence equations ( 1) -( 4) of present formulation can be readily converted to those of a rectangular cavity just by putting D = 0 with R and A computed using the width of the rectangle.Hence, to test the applicability and accuracy of the present numerical method, the present numerical results was validated by performing simulations for natural convection in a square cavity (D = 0) and in a cylindrical cavity (D = 0) without surface tension forces.Results are compared with various existing solutions and are shown in Fig. 2. From the figure, it can be observed that the present results for both the cavities are found to be in good agreement with the existing results.

Results and discussion
The numerical study was conducted to determine the effect of surface tension on buoyancy driven convection in a vertical cylindrical annular cavity filled with a low Prandtl number fluid (P r = 0.054).Computations were carried out for a wide range of Rayleigh numbers (10 4 ≤ Ra ≤ 10 6 ) and Marangoni numbers (0 ≤ M a ≤ 10 5 ) for an enclosure of aspect ratios A = 1 and 2. The flow pattern and temperature distribution in the annular cavity are presented by means of contour plots of streamlines and isotherms.The rate of heat transfer was estimated by evaluating the average Nusselt number.Fig. 3 illustrates the streamline and isothermal pattern for Rayleigh number Ra = 10 4 with and without surface tension effect.In the absence of shear at the free surface (M a = 0), the flow exhibits a simple recirculating pattern.The fluid rises along the hot wall, where it is heated, encounters the top adiabatic wall, travels towards the cold wall and recirculates.Also, the centre of rotation is in the middle of the annular cavity.The isotherms show less temperature stratification inside the annular cavity due to the low value of Ra.For low Marangoni number (M a = 10 2 ), the flow pattern and temperature distribution remains the same.This shows the effect of surface tension is unnoticeable for small values of M a.However, when M a = 10 4 , the surface tension effect leads to the development of a pair of counter rotating cells at the top and bottom of the cavity.The top cell is driven by thermocapillary forces, while the bottom cell is due to buoyancy.Also, the maximum stream function value of the main cell is less.This is because thermocapillary energy is utilized to generate a counter rotating cell.The temperature distribution is also strongly affected.There is a pulling of isotherms towards the hot wall.By further increasing the Marangoni number (M a = 10 5 ), surface tension induced shear increases the free surface velocities, and the surface tension driven cell becomes stronger than the buoyancy driven cell.This can be seen from the maximum stream function value.This effect is also reflected in the temperature distributions.The isotherms are more compact near the upper left corner.Also the isotherms are pulled in more in the middle of the cavity.The streamlines and isotherms for moderate Rayleigh number (Ra = 10 5 ) are shown in Fig. 4. A simple recirculating flow pattern and temperature stratification in the vertical direction was observed for M a = 10 4 .When the Marangoni number is increased to M a = 10 5 , two counter rotating cells developed, one being at the top of the annular enclosure due to surface tension.The buoyancy driven cell at the bottom of the annular cavity is bigger in size compared to the surface tension driven cell.
As the Rayleigh number increases further (Ra = 10 6 ), Fig. 5, buoyancy induced convection is dominant and the surface tension effect is noticeable only for large values of M a.For low and moderate values of M a, the isotherms and streamlines resemble that of pure buoyancy driven convection.Also temperature stratification exists in the vertical direction.When M a = 10 5 , there are two counter rotating cells, the smaller cell at top left corner, driven by surface tension, and the another cell, driven by buoyancy, occupying the major portion of the cavity.In this case, pulling of isotherms was found to be more due to the combined buoyancy effect.The influence of aspect ratio is shown in Fig. 6.For low Rayleigh number and M a = 10 3 , a single convective cell driven by buoyancy force exists in the bottom of the cavity.The streamlines indicate that the fluid at the top is attracted towards the top free surface by surface tension pull.The same phenomena can also be observed from the isotherms, as they are parallel near the top of the cavity.When the Marangoni number is further increased to 10 5 , for the same Rayleigh number, thermocapillary effect becomes stronger and two cells are produced in the cavity.The surface tension driven cell at the top of the cavity is stronger than the buoyancy driven cell as is evident from the maximum stream function value.This effect is reflected in isotherms as the isotherms are attracted towards the top free surface.However as the Rayleigh number increases to 10 6 , for M a = 10 3 , a single convective cell occupying the entire cavity was observed.For M a = 10 5 , there is a secondary cell at the top of the annulus, which is comparatively smaller than the buoyancy driven cell.The isotherms show the temperature stratification inside the cavity along the vertical direction.Also, the formation of thermal boundary layers along the side walls can be observed.Fig. 7 shows vertical velocity profiles at the midheight of the enclosure for different values of Ra and M a.At Ra = 10 4 , the vertical velocity increases with Marangoni number.However, there is a reduction in the magnitude of vertical velocity for M a = 10 4 .The same kind of phenomena is also observed as the Rayleigh number is increased to 10 5 .Fig. 8 is the temperature profiles at the midheight of the annular enclosure for two values of Ra and for M a = 0 − 10 4 .It can be observed, from the figures, that the temperature stratification reduces for M a = 10 4 .

Conclusion
The effect of surface tension on buoyancy-driven convection in a vertical cylindrical annular cavity has been numerically studied.Heat transfer and fluid flow results have been presented.Distinct flow regimes in the steady state, have been obtained, for several combinations of the physical parameters.The numerical results indicate that, thermocapillary forces induce multicellular flows even in smaller aspect ratio cavities.For high M a, a bicellular flow is found in the cavity, with a stronger surface tension driven cell, at the top of the cavity.Thermocapillary force has more effect in tall cavities, (A = 2), rather than square cavities (A = 1).Quantitative results were presented in terms of overall Nusselt number.The heat transfer rate increases with Marangoni number.The present numerical results, in the absence of thermocapillary forces, are in good agreement with the existing results for both cylindrical and rectangular cavities.

Fig. 8 .
Fig. 8. Temperature profiles at midheight for A = 1.The variation of average Nusselt number with dimensionless time τ is shown in Fig. 9.During early stages of the flow development, the heat transfer is mainly due to pure conduction.During the transition from conduction to convection heat transfer mode, the Nusselt number curve passes through a minimum.The average Nusselt number then increases slightly and levels off to reach the steady state.The Nusselt number increases with Rayleigh number due to increased convection.