Mechanism of swirl generation in sink flow/Sukurio susidarymo nutekamajame vamzdyje mechanizmas.

Mohammadi, J. ; Karimi, H. ; Hamedi, M.H. 等

1. 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-5]. The latter effect has been found to be negligibly 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-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 three-dimensional (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.

[FIGURE 1 OMITTED]

2. Governing equations

We consider the flow inside a cylindrical container as sketched in Fig. 1, a. A flow rate Q of an incompressible fluid enters the cylindrical tank horizontally through the circular side. The fluid flows out through a small orifice of radius [R.sub.p] 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 [L.sub.p]. The cylindrical polar coordinates ([r.sup.*], [theta], [z.sup.*]), with the velocity field ([V.sub.r.sup.*], [V.sub.[theta].sup.*], [V.sub.z.sup.*]) 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:

r = [r.sup.*]/R; z = [z.sup.*]/H; (1)

[V.sub.r] = [V.sub.r.sup.*] 2[pi] RH/Q; [V.sub.[theta]] = [V.sub.[theta].sup.*] 2[pi] RH/Q; [V.sub.z] = [V.sub.z.sup.*] 2[pi] RH/Q. (2)

In order to model the incompressible axisymmetric flow, the stream function-vorticity-circulation formulation [19] is used:

[V.sub.r] = [LAMBDA] 1/r [partial derivative][psi]/[partial derivative]z; [V.sub.z] = 1/r [partial derivative][psi]/[partial derivative]r; [V.sub.[theta]] = [GAMMA]/r; (3)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (4)

where the stream function [psi], vorticity [eta], and circulation [GAMMA] (related to angular momentum per unit mass) are defined as:

[psi] = [[psi].sup.*] 2[pi]H/RQ; [GAMMA] = [[GAMMA].sup.*] 2[pi]H/Q; [eta] = [[eta].sup.*] 2[pi][R.sup.2]H/Q.

Based on the above formulation, the continuity equation is satisfied and the Navier-Stokes equations for the steady axisymmetric flow can be written as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]; (5)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (6)

where Re = QH/2[pi] [R.sup.2]v and [LAMBDA] = R/H.

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 Fr = [square root of [Q.sup.2]/g[H.sup.5]], is always lower than 0.5 value. For this Froude number, Hocking et al. [20] have shown 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].

3. 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:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (7)

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:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (8)

where [[bar.[GAMMA]].sub.in] is the average value of [GAMMA] at the entrance (r = 1). The relation [([[bar.V].sub.[theta].sup.*]).sub.in] / [([[bar.V].sub.r.sup.*]).sub.in] = [[bar.[GAMMA]].sub.in] and [([[bar.V].sub.[theta]]).sub.in] = [[bar.[GAMMA]].sub.in] 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 [[bar.[GAMMA]].sub.in].

At the solid walls (sides C and D in Fig. 1, a), the velocity vanishes. Thus:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (9)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (10)

At the pipe exit (side E in Fig. 1, a), the velocity profile is assumed to be independent of z, i.e.:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (11)

On the axis of symmetry (side F in Fig. 1, a), we have:

[psi] = - 1/[LAMBDA]; [eta] = 0; [GAMMA] = 0; r = 0; - [L.sub.1] < z < 1. (12)

4. 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 [10.sup.-5] 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 x 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.

5. 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.sub.[theta]] and circulation [GAMMA] are shown in Figs. 2 and 3. Then, Figs. 4-6 show respectively the effects of aspect ratio [LAMBDA], drainpipe radius [R.sub.1] and average value of circulation at entrance [[bar.[GAMMA]].sub.in] on [V.sub.[theta]] and [GAMMA]. Finally, the numerical simulation of some related experiments are shown in Fig. 7.

Variations of [V.sub.[theta]] 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.sub.[theta]] in the drainpipe area become higher than the inlet value of [V.sub.[theta]]. 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.sub.[theta]] increases continuously near the drainpipe. For example, at Re = 5.3, the maximum value of [V.sub.[theta]] 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.sub.[theta]] 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.sub.[theta]] along a specific 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.sub.[theta]] along this specific streamline are lower than the [V.sub.[theta]] at the inlet. By increasing the Re number above this value, the values of [V.sub.[theta]] along the streamline will increase and become higher than its inlet value. Further increase in the Re number, causes the values of [V.sub.[theta]] on the streamline to approach a limit, which is determined from the relation [V.sub.[theta]] = [[GAMMA].sub.in]/r.

The contours of circulation [GAMMA] 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, respectively. 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 [GAMMA] with respect to the radius r (i.e., d[GAMMA]/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 [GAMMA] takes the value of zero. This result is consistent with the results obtained for the sink flow with a rotating body [22].

[FIGURE 2 OMITTED]

It can be seen from Fig. 2 that both [V.sub.[theta]] and [GAMMA] take their maximum value at the free surface (due to the stress-free boundary condition in the [theta]-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.sub.[theta]] and [GAMMA] as the Re number increases. In other words, there is no critical Re number above which the values of [V.sub.[theta]] and [GAMMA] could suddenly increase. This would imply that the flow is stable in the azimuthal direction. In addition, the value of [GAMMA] in the solution domain is less than its inlet value. Therefore, for the Re numbers within the range 0.265 < Re <5.3, the self-rotation phenomenon, i.e., the increase of r from its inlet value does not occur.

[FIGURE 3 OMITTED]

The effects of [LAMBDA], [R.sub.1] and [[bar.[GAMMA]].sub.in] on [V.sub.[theta]] and [GAMMA] are represented in Figs. 4-6, respectively. In the whole solution domain except near the drainpipe, the variations of [LAMBDA] and R1 do not have a noticeable effect on [V.sub.[theta]] and [GAMMA]. It can also be noted that the parameter [[bar.[GAMMA]].sub.in] has at best negligible effect on [V.sub.[theta]] / [([[bar.V].sub.[theta]]).sub.in] and [GAMMA]/[[bar.[GAMMA]].sub.in]. 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.sub.[theta]] and [GAMMA] are incorporated into the definition of the Re number. In addition, the value of the maximum velocity of the concentrated vortex [([V.sub.[theta]]).sub.max] increases with increasing Re, [LAMBDA] = R/H and [[bar.[GAMMA]].sub.in] on the one hand and decreasing [R.sub.1] on the other hand. This means that [([V.sub.[theta]]).sub.max] / [[bar.[GAMMA]].sub.in] is a function of [Re.sub.d] = Q/ (2[pi] v[R.sub.p]) = Re [LAMBDA] / [R.sub.1]. From the fact that [[bar.[GAMMA]].sub.in] = [([[bar.V].sub.[theta]]).sub.in], it be said that [([V.sub.[theta]]).sub.max] / [([[bar.V].sub.[theta]]).sub.in] is a function of either one of the following new defined Reynolds number, i.e., [Re.sub.d] = Re [LAMBDA]/[R.sub.1] or [Re.sub.d] = Q / (2[pi] v[R.sub.p]). In other words, the strength of the concentrated vortex, which forms close to the drainpipe, is a function of the drain flow rate, the drain-hole size, the liquid kinematic viscosity and the average inlet azimuthal velocity [([[bar.V].sub.[theta]]).sub.in].

[FIGURE 4 OMITTED]

[FIGURE 5 OMITTED]

[FIGURE 6 OMITTED]

Previously, researchers have studied the self-rotation phenomenon experimentally in various geometries, setups and for various liquids [10-13]. In their studies they have measured the values of the dimensional azimuthal velocity [V.sub.[theta].sup.*] 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.sub.[theta].sup.*] 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.sub.[theta].sup.*] are not known. Knowing the inlet values of [V.sub.[theta].sup.*] as is the case in the present numerical 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.sub.[theta].sup.*] 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.sub.[theta].sup.*] with respect to the drain hole has been performed are known. On the other hand, the values of [V.sub.[theta]] at different Re numbers can be obtained from Fig. 3, b. Now, from the definition of the Re number, i.e., Re = QH / (2[pi][R.sup.2]v), the drain flow rate corresponding to each value of [V.sub.[theta]] is obtained. Then, the values of dimensional azimuthal velocity [V.sub.[theta].sup.*] are calculated by using the [V.sub.[theta]] definition [V.sub.[theta]] = [V.sub.[theta].sup.*] 2[pi] RH/Q.

The results shown in Fig. 7 obtained in the present numerical simulations are associated with the experimental works found in [10-13]. In all of these simulations, the average value of the inlet circulation [[bar.[GAMMA]].sub.in] is chosen to be equal to 0.087. It should be noted that the value of [[bar.[GAMMA]].sub.in] is unknown in the experimental studies. By taking this value for [[bar.[GAMMA]].sub.in], the azimuthal velocity deviates by 5 [degrees] 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 [[bar.[GAMMA]].sub.in] does not affect the behaviors of [V.sub.[theta]] and [GAMMA] 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 [Q.sub.c]. 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-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.sub.[theta].sup.*] values experience very high variations along the flow direction. Since in the experimental studies the variations of [V.sub.[theta].sup.*] have been measured, the measurement approach can cause erroneous interpretation of the data. Therefore, it seems to be reasonable to consider the circulation [[GAMMA].sup.*] value rather than the [V.sub.[theta].sup.*] value;

b) the measurements were restricted to some specific points or lines but not to the whole flow domain. Specifically, the [V.sub.[theta].sup.*] 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.

[FIGURE 7 OMITTED]

6. 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., [[bar.[GAMMA]].sub.in]) and the azimuthal velocity approaches the limit [[bar.[GAMMA]].sub.in]/r 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 [[bar.[GAMMA]].sub.in] both the circulation and the azimuthal velocity are functions of the Re number only and are independent of the geometric parameters [LAMBDA], [R.sub.1] and [L.sub.1]. The formation of the concentrated vortex at the drainpipe entrance is observed to occur at the Re numbers higher than 0.294. It is also shown that the strength of the concentrated vortex or the maximum azimuthal velocity [([V.sub.[theta]]).sub.max] for a specific [[bar.[GAMMA]].sub.in] is a function of the Reynolds number defined as [Re.sub.d] = Q / (2[pi] v[R.sub.p]).

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.

Nomenclature

H--height of the free surface, m; [L.sub.p]--drainpipe length, m;

[L.sub.1] - dimensionless drainpipe length, [L.sub.1] = [L.sub.p]/H; Q - volumetric drain flow rate, [m.sup.3][S.sup.- 1]; [R.sub.p] - drainpipe radius, m;

[R.sub.1] - dimensionless drainpipe radius, [R.sub.1] = [R.sub.p]/R; R - radius of the cylinder, m; [r.sup.*] -radial coordinate, m; r - dimensionless radial coordinate, r = [r.sub.*]/R; Re - Reynolds number, Re = QH/(2[pi] [R.sup.2] v); [V.sub.r.sup.*] - radial velocity, m[s.sup.-1];

[V.sub.[theta].sup.*] - azimuthal velocity, m[s.sup.-1]; [V.sub.z.sup.*] - axial velocity, m[s.sup.- 1];

[V.sub.r.sup.*] - dimensionless radial velocity, [V.sub.r] = [V.sub.r.sup.*] 2[pi] RH/Q;

[V.sub.[theta]] - dimensionless azimuthal velocity, [V.sub.[theta]] = [V.sub.[theta].sup.*] 2[pi] RH/Q;

[V.sub.z] - dimensionless axial velocity, [V.sub.z] = [V.sub.z.sup.*] 2[pi] RH/Q;

[([[bar.V].sub.r]).sub.in] - average value of [V.sub.r] at entrance; [([[bar.V].sub.[theta]]).sub.in] - average value of [V.sub.[theta]] at entrance; [([V.sub.[theta]]).sub.max] - maximum value of [V.sub.[theta]] of concentrated vortex; [z.sup.*] - axial coordinate, m; z - dimensionless axial coordinate, z = [z.sup.*]/H.

Greek symbols--

[upsilon] - kinematic viscosity, [m.sup.2][s.sup.-1]; [LAMBDA] - aspect ratio of cylinder, [LAMBDA] = R/H; [[bar.[GAMMA]].sub.in] - average value of [GAMMA] at entrance; [[GAMMA].sup.*] - circulation, i.e., angular momentum per unit mass, [m.sup.2] [s.sup.-1]; [GAMMA] - dimensionless circulation, [GAMMA] = [[GAMMA].sup.*] 2[pi] H/Q .

cross ref http://dx.doi.org/ 10.5755/j01.mech.19.6.5998

Received October 19, 2012

Accepted November 11, 2013

References

[1.] Klimenko, A.Y. 2001. Moderately strong vorticity in a bathtub-type flow, Theoretical and Computational Fluid Dynamics 14(4): 243-257. http://dx.doi.org/10.1007/s001620050139.

[2.] Vladimir, S.; Fazle, H. 1999. Collapse, symmetry breaking and hysteresis in swirling flows, Annual Review of Fluid Mechanics 31(1): 537-566. http://dx.doi.org/10.1146/annurev.fluid.31.1.537.

[3.] Lugt, H.J. 1983. Vortex Flow in Nature and Technology, John Wiley & Sons, New York, 297 p.

[4.] Shapiro, A.H. 1962. Bath-tub vortex, Nature, 196(4859): 1080-1081. http://dx.doi.org/10.1038/1961080b0.

[5.] Trefethen, L.M.; Bilger, R.W.; Fink, P.T.; Luxton, R.E.; Tanner, R.I. 1965. The bath-tub vortex in the southern hemisphere, Nature 207(5001): 1084-1085. http://dx.doi.org/10.1038/2071084a0.

[6.] Herrada, M.A.; Barrero, A. 2002. Self-rotation in electrocapillary flows, Physical Review E 66(3): 036311-036320. http://dx.doi.org/10.1103/PhysRevE.66.036311.

[7.] Torrance, K.E. 1979. Natural convection in thermally stratified enclosures with localized heating from below Journal of Fluid Mechanics 95(3): 477-495. http://dx.doi.org/10.1017/S0022112079001567.

[8.] Funakoshi, M.; Inoue, S. 1988. Surface waves due to resonant horizontal oscillation, Journal of Fluid Mechanics 192: 219-247. http://dx.doi.org/10.1017/S0022112088001843.

[9.] Bojarevics, V.; Freibergs, Y.A.; Shilova, E.I.; Shcherbinin, E.V. 1989. Electrically Induced Vortical Flows: Kluver Academic, 380 p.

[10.] Tanaka, D.; Mizushima, J.; Kida, S. 2004. The origin of the bathtub vortex, Kyoto University Research Information Repository 1406: 166-177 (in Japanese).

[11.] De Felice, V.F. 2008. The free surface vortex due to instability, PhD Thesis, Degli Studi di Salerno University (in French).

[12.] Kawakubo, T.; Tsuchiya, Y.; sugaya, M.; Matsumura, K. 1978. Formation of a Vortex around a Sink, a kind of phase transition in a nonequilibrium open system, Physics Letters A 68(1): 65-66. http://dx.doi.org/10.1016/0375-9601(78)90759-4.

[13.] Fernandez-Feria, R.; Sanmiguel-Rojas, E. 2000. On the appearance of swirl in a confined sink flow, Physics of Fluids, 12(11): 3082-3085. http://dx.doi.org/10.1063/1.1313566.

[14.] Fernandez-Feria, R. 2002. Stability analysis of boundary layer flow due to the presence of a small hole on a surface, Physical Review E 65(3): 036307. http://dx.doi.org/10.1103/PhysRevE.65.036307.

[15.] Sanmiguel-Rojas, E. 2002. On The Phenomenon Of Self-Rotation, Ph.D. thesis, Malaga University (in French).

[16.] Sanmiguel-Rojas, E.; Fernandez-Feria, R. 2006. Nonlinear instabilities in a vertical pipe flow discharging from a cylindrical container, Physics of Fluids 18(2): 024101. http://dx.doi.org/10.1063/L2168445.

[17.] Tyvand, P.A.; Haugen, K.B. 2005. An impulsive bathtub vortex, Physics of Fluids 17(6): 062105. http://dx.doi.org/10.1063/1.1938216.

[18.] Drazin, P.G. 2002. Introduction to hydrodynamic stability: Cambridge University Press, 258 p. http://dx.doi.org/10.1017/CBO9780511809064.

[19.] Hoffmann, K.A.; Chiang, S.T. 2000. Computational Fluid Dynamics, Vol. I: Wichita.

[20.] Hocking, G.C.; Vanden-broeck, J.M.; Forbes, L. K. 2002. A note on withdrawal from a fluid of finite depth through a point sink, ANZIAM 44: 181-191. http://dx.doi.org/10.1017/S1446181100013882.

[21.] Bohling, L.; Andersen, A.; Fabre, D. 2010. Structure of a steady drain-hole vortex in a viscous fluid, Journal of Fluid Mechanics 656: 177-188. http://dx.doi.org/10.1017/S0022112010001473.

[22.] Yukimoto, S.; Niino, H.; Noguchi, T.; Kimura, R. 2010. Structure of a bathtub vortex: importance of the bottom boundary layer, Theoretical and Computational Fluid Dynamics 24: 323-327. http://dx.doi.org/10.1007/s00162-009-0128-3.

J. Mohammadi, Faculty of Mechanical Engineering, K. N. Toosi University of Technology, Tehran, Iran, E-mail: j_mohammadi@dena.kntu.ac.ir, mohammadijalal@yahoo.com

H. Karimi, Faculty of Aerospace Engineering, K. N. Toosi University of Technology, Tehran, Iran, E- mail: karimi@kntu.ac.ir

M. H. Hamedi, Faculty of Mechanical Engineering, K. N. Toosi University of Technology, Tehran, Iran, E-mail: hamedi@kntu.ac.ir

A. Dadvand, School of Mechanical Engineering, Urmia University of Technology, Urmia, Iran, E-mail: a.dadvand@mee.uut.ac.ir