Mechanism of swirl generation in sink flow

J. Mohammadi*, H. Karimi**, M. H. Hamedi***, A. Dadvand**** *Faculty of Mechanical Engineering, K. N. Toosi University of Technology, Tehran, Iran, E-mail: j_mohammadi@dena.kntu.ac.ir, mohammadijalal@yahoo.com **Faculty of Aerospace Engineering, K. N. Toosi University of Technology, Tehran, Iran, E-mail: karimi@kntu.ac.ir ***Faculty of Mechanical Engineering, K. N. Toosi University of Technology, Tehran, Iran, E-mail: hamedi@kntu.ac.ir ****School of Mechanical Engineering, Urmia University of Technology, Urmia, Iran, E-mail: a.dadvand@mee.uut.ac.ir


Introduction
The bathtub vortex (i.e., sink flow with swirl velocity) is a well-known phenomenon.When water is drained from a tank through a small hole, it experiences a translational motion toward the hole and a rotational movement [1].Due to many effective parameters involved and the sensitivity of the bathtub vortex to external factors, mechanism of swirl generation has not been elucidated sufficiently [2].
The appearance of swirl in the sink flow and the formation of the bathtub vortex can result from many factors.These include small asymmetries in the flow inlet, asymmetric temperature distribution, asymmetric air motion over the water surface, asymmetric initial or boundary conditions, residual fluid motion in the vessel, vessel vibration, and effect of the Coriolis force due to the Earth's rotation [3][4][5].The latter effect has been found to be negli-gibly small provided the diameter of the tank is smaller than six feet [4,5].The basic unanswered questions are "will swirl appear in sink flow if all the external factors are eliminated and the tank size is chosen sufficiently small?" and "will swirl appear if the speed of a swirl-free flow exceeds some threshold value?" The appearance of swirl in the fluid flow without any external factors, namely self-rotation [6], have been observed in many natural systems such as the liquid flow inside Taylor cones [6], natural convection flow in a vertical circular cylinder [7], horizontally oscillating water in a cylindrical container [8], and an electrically driven flow of mercury in a cup [9].
Fundamentally, there have been several experimental and numerical studies to investigate the possibility of self-rotation phenomenon and the related instability in the sink flow.However, there is no consensus among the researchers about either the possibility or impossibility of this phenomenon.The experimental studies showed that the swirl would appear in an initially swirl-free sink flow as the Reynolds number based on the sink flow rate is increased above a critical value [10][11][12][13].In addition, this phenomenon has been investigated numerically for different geometries and conditions as follows.A linear stability analysis of the boundary layer in sink flow was carried out by Fernandez [14].He observed that the instability occurred when the Reynolds number was relatively high and the flow became turbulent.This instability has nothing to do with the formation of a vortex in the sink, a phenomenon that is shown experimentally to occur at much lower Reynolds number.The stability of sink flow was studied numerically based on the axisymmetric and threedimensional (3D) models [11,15,16].It was observed in [15,16] that the flow was stable and swirl-free for all the Reynolds numbers tested.Felice [11] found no swirl in the axisymmetric model.In 3D model, however, both instability and swirl were observed for the Re numbers above a critical value.
In fact, one can see different and conflicting findings from the above mentioned studies, that is, some researchers accept the existence of self-rotation in the sink flow and some refuse.On the other hand, it is deduced from these studies that the circulation generation contradicts the conservation of angular momentum [2].Therefore, this problem needs to be studied in more detail and characterized more accurately.
What we call self-rotation in the current study is either the spontaneous appearance of swirl in an initially swirl-free flow (i.e., swirl generation) or the increase of circulation with respect to its inlet value (i.e., circulation generation).In both cases, the circulation value increases.Note that, increase of the swirl (azimuthal) velocity by decreasing the radius with constant circulation as it occurs in converging flows is not related to self-rotation.This strong growth of local azimuthal velocity near the bathtub drain hole can be explained by Lord Kelvin's circulation theorem.This theorem states the conservation of circulation and it is the hydrodynamic version of a more general principle: the conservation of angular momentum [17].
In the present work, the possibility of circulation generation in the sink flow is studied numerically.We simulate the sink flow or the bathtub vortex in a cylindrical tank with a central drainpipe and a free surface (Fig. 1).The axisymmetric configuration of the problem is simulated using direct numerical simulation (DNS), where the full Navier-Stokes equations are solved without any turbulence model introduced [18].Moreover, we attempt to extract the parameters affecting the circulation and azimuthal velocity and present a more appropriate definition for the Reynolds number that could influence the circulation.The numerical simulations are performed for a wide range of Reynolds numbers and the flow instability in the azimuthal direction as well as the abrupt changes in the circulation related to self-rotation are examined.Finally, it is attempted to discuss the results obtained from the previous experimental studies on the self-rotation phenomena.
The rest of the paper is organized as follows: Sections 2-4 contain the governing equations of motion, boundary condition and the numerical procedure, respectively.In Section 5, numerical results are presented and discussed.Finally, in Section 6, we summarize the main findings and present the conclusions.

Governing equations
We consider the flow inside a cylindrical container as sketched in Fig. 1, a.A flow rate Q of an incompress-ible fluid enters the cylindrical tank horizontally through the circular side.The fluid flows out through a small orifice of radius p R centered at the bottom end wall.The flow is assumed axisymmetric and steady state.We denote the radius of the tank R, the height of the free surface H, and the length of the drainpipe is considered where the asterisk superscripts denote dimensional quantities.The z-axis is chosen as the axis of the cylinder and the bottom end wall lies in the plane z = 0.The corresponding dimensionless variables are: ; ; In order to model the incompressible axisymmetric flow, the stream function-vorticity-circulation formulation [19] is used: 11 , where the stream function ψ, vorticity η, and circulation Γ (related to angular momentum per unit mass) are defined as: ; .
Based on the above formulation, the continuity equation is satisfied and the Navier-Stokes equations for the steady axisymmetric flow can be written as: The Reynolds number (Re) defined in this manner could be an effective parameter influencing the circulation evolution in the solution domain as will be seen later in the results Section.
In the present work, the free surface of the flow inside the cylindrical container is assumed to be flat.This assumption is reasonable when the Froude number is small enough [11].The Froude number of sink flow studies, defined as , is always lower than 0.5 value.For this Froude number, Hocking et al. [20] have shown r z that the maximum depression in the free surface height is less than 0.003 of liquid height.This value is very small compared to the water height.Therefore, in the numerical study the free surface can approximately be assumed flat.For further information on this regard, the reader is referred to [11,21].

Boundary conditions
At the free surface (side A in Fig. 1, a), the axial velocity is set to be zero based on the assumption of constant height of the liquid and there is no shear stress so that the axial differentiation of the circulation is zero: At the entrance (side B in Fig. 1, a), the axial velocity is taken to be zero and we assume a parabolic profile for both the radial velocity and the circulation.These assumptions have been made because they match the boundary conditions at the bottom wall and at the free surface: where in Γ is the average value of Γ at the entrance (r = 1).
can be obtained by considering Eqs. ( 2) and (3).This means that the ratio of the average value of azimuthal velocity to the average value of radial velocity is equal to in Γ .At the solid walls (sides C and D in Fig. 1, a), the velocity vanishes.Thus: At the pipe exit (side E in Fig. 1, a), the velocity profile is assumed to be independent of z, i.e.: .
On the axis of symmetry (side F in Fig. 1, a), we have:

Computational approach
A finite difference approach is employed to discretize the system of Eqs. ( 4)-( 6) subject to the boundary conditions ( 7)- (12).In addition, the successive overrelaxation (SOR) method is used to solve the discretized equations iteratively.The iteration process continues until the maximum difference between two successive iterated values is less than 5 10  in magnitude.The central difference scheme is used to discretize the space derivatives.For higher values of the Reynolds number (Re), central differencing in the convective terms may lead to numerical instability.To solve this problem, we have used the upwind difference scheme in the convective terms while computing for large values of Re number.A non-uniform grid is utilized in the r and z directions and refined near the inlet of the drainpipe where the velocity gradients are large.It has been found that 201 × 101 grids produce a grid independent solution.To substantiate the accuracy of the present numerical approach, the current results have been compared with those found in [21], where identical geometry as ours has been considered.

Results and discussion
In this section, the results are presented and discussed.First, the results corresponding to the stability of the sink flow and the effect of Re number on the azimuthal velocity θ V and circulation Γ are shown in Figs. 2 and 3.
Then, Figs.4-6 show respectively the effects of aspect ratio Λ , drainpipe radius 1 R and average value of circulation at entrance in Γ on θ V and Γ .Finally, the numerical sim- ulation of some related experiments are shown in Fig. 7.
Variations of θ V in terms of Re are shown in frames h-n of Fig. 2. It can be seen that at Re numbers higher than 0.294, the values of θ V in the drainpipe area become higher than the inlet value of θ V .In other words, the concentrated vortex is formed in the drainpipe area (for more information, refer to [21]).As the Re number increases, the azimuthal velocity θ V increases continuously near the drainpipe.For example, at Re = 5.3, the maximum value of θ V approaches 18, which is about 207 times its average value at the inlet.The main reason for this increase in the maximum value of θ V of the concentrated vortex is that the radius of the vortex core decreases as the Re number increases [21].
Fig. 3, b shows the variations of θ V along a specif- ic streamline, which starts from the middle height of the inlet (r = 1, z = 0.5).For the Re numbers lower than approximately 0.622, the values of θ V along this specific streamline are lower than the θ V at the inlet.By increasing the Re number above this value, the values of θ V along the streamline will increase and become higher than its inlet value.Further increase in the Re number, causes the values of θ V on the streamline to approach a limit, which is de- termined from the relation /r V in θ   .The contours of circulation Γ and its variations along the streamline starting from the middle height of the flow at the inlet are shown in Figs. 2, a-g and 3, a, respec-tively.It can be seen from Fig. 2, a-g that the circulation values in the solution domain are always lower than its inlet value.In addition, it can be seen from Fig. 3, a that for all the Re numbers tested, the rate of decrease in Γ with respect to the radius r (i.e., d dr


) reduces as r is reduced.This may be due to the increase in the radial velocity as radius decreases causing the ratio of inertia force to viscous force (i.e., the Re number) to increase.It should be mentioned that the viscous force is applied to the flow from the bottom-end wall.It is worth mentioning that, by increasing the Re number, the decrease in the circulation from the inlet towards the drainpipe will be attenuated so that at the (relatively) high Re numbers the circulation in the whole domain will approach its inlet value.Around r = 0, however, due to the continuity condition, the circulation Γ takes the value of zero.This result is consistent with the results obtained for the sink flow with a rotating body [22].
It can be seen from Fig. 2 that both θ V and Γ take their maximum value at the free surface (due to the stress-free boundary condition in the θ-direction) and zero value at the tank bottom floor due to the no-slip boundary condition.
It is also seen in Figs. 2 and 3 that no sudden changes would occur in the behavior of θ V and Γ as the Re number increases.In other words, there is no critical Re number above which the values of θ V and Γ could suddenly increase.This would imply that the flow is stable in the azimuthal direction.In addition, the value of Γ in the solution domain is less than its inlet value.Therefore, for the Re numbers within the range 0.265 < Re <5.3, the selfrotation phenomenon, i.e., the increase of Γ from its inlet value does not occur.

R/ 
The effects of Λ , 1 R and in Γ on θ V and Γ are represented in Figs.4-6, respectively.In the whole solution domain except near the drainpipe, the variations of Λ and R 1 do not have a noticeable effect on V θ and Γ.It can also be noted that the parameter in Γ has at best negligible ef- fect on


. This shows the fact that in the present study the Re number has been properly defined such that all the parameters affecting the behaviour of V θ and Γ are incorporated into the definition of the Re num- ber.In addition, the value of the maximum velocity of the concentrated vortex (V θ ) max increases with increasing Re, R/H Λ  and in Γ on the one hand and decreasing 1 R on the other hand.This means that      Previously, researchers have studied the selfrotation phenomenon experimentally in various geometries, setups and for various liquids [10][11][12][13].In their studies they have measured the values of the dimensional azimuthal velocity  θ V at various distances from the drain hole.In the present work, it is attempted to simulate numerically the self-rotation phenomenon in geometries identical to (but slightly different from) those used in the experiments to obtain the azimuthal velocity  θ V with respect to the drain flow rate Q.The main discrepancy between the present work and the previous experiments is the fact that in the latter the inlet values of  θ V are not known.Knowing the inlet values of  θ V as is the case in the present numeri- cal study will be useful for understanding the nature of the swirl appearance in the sink flow.Fig. 7 shows four different diagrams of the azimuthal velocity  θ V versus the drain flow rate obtained from the present simulations.Each diagram is associated with one of the four different experimental investigations mentioned above.These diagrams are obtained as follows.In the experiments, the positions at which the measurement of  θ V with respect to the drain hole has been performed are known.On the other hand, the values of θ V at different Re numbers can be obtained from Fig. 3 . The results shown in Fig. 7 obtained in the present numerical simulations are associated with the experimental works found in [10][11][12][13].In all of these simulations, the average value of the inlet circulation in Γ is chosen to be equal to 0.087.It should be noted that the value of in Γ is unknown in the experimental studies.By taking this value for in Γ , the azimuthal velocity deviates by 5° from the radial velocity direction.Thus, as it was pointed out previously, the quantity ratio of the two velocity components becomes 0.087.It should also be mentioned that, the value of in Γ does not affect the behaviors of θ V and Γ in the solution domain (Fig. 6).Fig. 7 shows that the azimuthal velocity is very low when the drain flow rate varies between zero and a certain value.Practically, the swirl appears for the drain flow rates higher than a threshold value, which is called the critical drain flow rate c Q .This may be due to the fact that when the Re number becomes lower than 0.662, the viscous forces will dominate (see Fig. 3) and that the inlet azimuthal velocity will increase as the drain flow rate increases.The phenomenon of swirl appearance for the drain flow rates higher than some critical value is consistent reasonably with the corresponding experimental observations.The main discrepancy between the numerical and experimental results observed in some test cases is due to the differences between the different geometries employed.
The experimental studies [10][11][12][13] claimed that the swirl appears in the sink flow when either the drain flow rate or the Re number reaches some critical value.All of these studies have interpreted this phenomenon based on the flow instability without considering the effect(s) of the external factors.Such interpretations seem to be unconvinced for the following reasons: a) the  θ V values experience very high variations along the flow direction.Since in the experimental studies the variations of  θ V have been measured, the measurement approach can cause erroneous interpretation of the data.Therefore, it seems to be reasonable to consider the circulation *  value rather than the  θ V value; b) the measurements were restricted to some specific points or lines but not to the whole flow domain.Specifically, the  θ V values were not determined at the flow inlet and were assumed to be zero, which seems to be an erroneous assumption; c) in some 2D and 3D numerical studies (see [14,16], for examples), the self-rotation phenomenon in the sink flow was not observed explicitly and thus it was disproved; d) the self-rotation phenomenon in the sink flow is in contradiction with the principle of the angular momentum conservation [2]; e) finally, as it was shown earlier, the sink flow is stable at least in the axisymmetric case and the viscous forces are the main cause for the non-appearance of the swirl for the drain flow rates lower than the critical value.As can be seen, azimuthal velocity (swirl) appears as drain flow rate increases above of a critical value

Conclusions
This study aims to investigate the possibility of circulation generation in the sink flow in the absence of external factors using direct numerical simulation (DNS) based on the axisymmetric Navier-Stokes equations.A new interpretation regarding the origin of swirl appearance in the sink flow is presented.The geometry of the problem consists of a cylinder with the circular hole made on its bottom end.Fluid enters the cylinder horizontally through its lateral wall and is drained from the drainpipe.
The simulations are performed for a range of Re numbers, i.e., 0.265 < Re < 5.3.The results show that the value of circulation in the solution field is always less that its inlet value.By increasing the Re number, the circulation approaches its inlet value   i.e., in  and the azimuthal ve- locity approaches the limit /r Γ in in the whole solution domain, except in the region very close to the symmetry axis.In other words, due to the non-abrupt changes of circulation with respect to the Re number variations, the flow is said to be stable.In addition, the non-increase of circulation with respect to its inlet value indicates that there is no self-rotation in the sink flow.
In the whole solution domain except near the symmetry axis, for a fixed The numerical simulations carried out here in the geometries identical to the previous experiments show that at low drain flow rates, the viscous forces due to the bottom-end wall eliminates the inlet azimuthal velocity.By increasing the drain flow rate above some critical value, the effects of the viscous forces will decrease and swirl appears in the sink flow.This implies that the swirl appearance observed in the previous experiments was not related to the self-rotation and the flow instability.It is worth mention that this new interpretation of the origin of swirl appearance observed in the experiments assists to minimize the present contradictions between the numerical and experimental studies.

Fig. 1 a
Fig. 1 a) schematics of the axisymmetric sink flow with the inlet at the outer edge and the drainpipe at the centre of the bottom-end wall; b) The mesh used in the numerical simulations with R/H = 0, and 20 1 / /R R p 

Fig. 4
Fig. 4 Effect of R / H   Fig. 5 Effect of 1 p R R / R  on circulation  and azimuthal velocity θ V at Re=2.65, 10  Λ , and 087 0. in

Fig. 6
Fig. 6 Effect of in Γ on circulation Γ and azimuthal velocity θ V at Re=2.65, 10  Λ and 1 1 20 R/  , b.Now, from the definition of the Re number, i.e., flow rate corresponding to each value of θ V is obtained.Then, the values of dimen- sional azimuthal velocity  θ V are calculated by using the θ

,
Fig. 7 The dimensional azimuthal velocity  θ V versus the drain flow rate with 087 0. in Γ 

inΓ 1 R and 1 L
both the circulation and the azimuthal velocity are functions of the Re number only and are independent of the geometric parameters  , .The formation of the concentrated vortex at the drainpipe entrance is observed to occur at the Re numbers high-er than 0.294.It is also shown that the strength of the concentrated vortex or the maximum azimuthal velocity