Low Reynolds Number Flow Around Tori of Different Slenderness Γ

We present the results of axisymmetric and three-dimensional numerical simulations of the flow around a torus in the low Reynolds number regime [10−2; 3× 103] and aspect-ratio 0 < Γ = 2a/R ≤ 2 (core diameter over toroidal radius). It is shown that as Γ → 0, consistent with intuition and the results from the literature, the flow tends to that of a two-dimensional circular cylinder while for Γ→ 2 the flow is that around a bluff obstacle. It has been observed that in a small region of the Re–Γ phase space, the flow develops an axisymmetric recirculation detached from the torus and with a vorticity distribution that resembles the Hill vortex. In the range 0 < Γ ≤ 2, several different regimes have been observed; the peculiarities of each regime are analyzed and, whenever possible, similarities and differences with other classical flows are discussed.


Introduction
Flows around circular cylinders and spheres have been attracting the interest of researchers for more than a century owing to their relatively easy and controllable conditions, as well as the interesting flow physics [1,2].These problems have often been considered as prototypes of flows around bluff bodies which are relevant to engineering and practical applications.In these cases, the body geometry is described by a single parameter (the diameter d of the cylinder or of the sphere) and the flow is characterized only by the Reynolds number Re d = Ud/ν, with U the free stream velocity and ν the kinematic viscosity of the fluid.
A possible extension of the above flow paradigms consists of considering a geometry that is described by two curvature parameters instead of one.We believe that in this way it is possible to establish the standard for a new benchmark that, even if yet to be well defined and sufficiently simple, shows a richer and more complex flow physics with respect to circular cylinders and spheres.
In the present paper we have considered the flow around a toroidal ring with radius R and core diameter 2a in an attempt to increase the complexity of the problem.In fact, it is possible to change, in addition to the Reynolds number Re = RU/ν, also the slenderness parameter Γ = 2a/R.Apart from the paper by [3] that used a rarefied hypersonic flow, this problem has been considered for an incompressible flow by several studies in the past.Reference [4] considered "fat" tori (Γ = 2) in the Stokes regime, while [5] investigated only slender tori (low Γ) in the very low Reynolds number range.Reference [6] conducted experiments to measure the drag in axisymmetric conditions, and [7][8][9] performed accurate and extensive numerical simulations and stability analyses showing that all the dynamics from the flow around a sphere (Γ → ∞) to that around a circular cylinder (Γ = 0) can be obtained by varying the cylinder slenderness.
In this study we have analyzed the low Reynolds number range [10 −2 ; 3 × 10 3 ], and the aspect ratio has been spanned from Γ = 0.1 up to the limit Γ = 2; the results are summarized in the phase diagram of Figure 1, confirming that indeed the aspect ratio Γ plays a fundamental role in the determination of the flow features.The novel finding of this study is the existence of a new regime in the Re-Γ phase space, indicated by H in Figure 1.In this region the flow develops a spherical vortex in the wake detached from the surface of the torus, whose vorticity is distributed as in the Hill spherical vortex [10].It is worth mentioning that the phase diagram of Figure 1 has been obtained from the results of about 100 numerical simulations, more than 20 of which are fully three-dimensional.The boundaries between the regions have been somehow arbitrarily set by observing the flow behaviour and gathering the simulations into homogeneous categories.Therefore, the lines dividing the different regions should not be interpreted as sharp fronts, but rather as regions where the transition occurs.
The paper is organized as follows: in Section 2, the numerical setup is briefly described, the problem is introduced, and the results are validated.Section 3 contains the results with an analysis of the various flow regimes.The paper is closed by Section 4 with the conclusions.

Problem and Numerical Setup
In this paper, the Navier-Stokes equations for an incompressible viscous flow are numerically integrated; their nondimensional forms read: with u the velocity, p the pressure, D/Dt the material derivative, and f a forcing for the immersed boundary method (specified later).Equation ( 1) have been spatially discretized in a cylindrical coordinate system using staggered central second-order accurate finite-difference approximations.
Details of the numerical method are given in [11], and only the main features are summarized here.
The discretized system is integrated using a fractional-step method where the viscous terms are computed implicitly and the convective terms explicitly.At each time step, the momentum equations are provisionally advanced in time using the pressure at the previous time step, giving an intermediate non-solenoidal velocity field.A scalar quantity Φ is then introduced to project the non-solenoidal field onto a solenoidal one.A hybrid low-storage third-order Runge-Kutta scheme is used to advance the equations in time.
The forcing f of Equation ( 1) is introduced to handle complex geometric configurations without a body conformal mesh; this is the immersed boundary method whose details can be found in [12].In a few words, if Equation ( 1) is rewritten as ∂u/∂t = RHS + f, with RHS containing the advective, pressure, and viscous terms, its time discrete version reads: Here the superscripts indicate the discrete time levels and n + 1/2 indicates the midpoint between times t n and t n+1 .If at the immersed boundary the fluid velocity is imposed as u n+1 = u b , Equation (2) can be explicitly solved for f n+1/2 , yielding which is imposed at all immersed interface nodes while it results f n+1/2 = 0 everywhere else.This procedure allows the imposition of a velocity boundary condition on an arbitrary surface, and is consistent with the centered second-order finite-difference approximation; therefore, the overall accuracy of the scheme remains second-order [12].The problem is sketched in Figure 2, and is summarized as follows: a torus of radius R and core radius a is normally placed along the axis of a cylindrical domain of radius L r and length L x , respectively, in the radial r and axial x directions.The torus is at a distance L i from the inflow (x = 0).The flow boundary conditions are: In the boundary conditions at (x = L x ), all the velocity components are advected outside the domain with a velocity c which is dynamically determined so as to assure the free-divergence of the velocity field at every time step.
The computational domain is discretized by a mesh which is non-uniform in the radial and axial directions to cluster the gridpoints in the torus region and in the wake (see the example in the sketch of Figure 2); for the three-dimensional cases, the mesh is uniform in the azimuthal direction.The axisymmetric cases have been mainly run on two different grids 120 × 201 and 258 × 401 in the radial and axial directions, respectively, obtaining a drag coefficient C D (computed by integrating viscous stresses and pressure over the torus surface and projecting the normalized force in the streamwise direction) that typically showed a difference below 1%.In many cases-especially those with low Γ and high Re-a third grid of 515 × 801 nodes has been used to further assess the grid independence of the results again with differences below 1%, which is within the statistical uncertainty of the data.The standard domain dimensions are L r = 16R and L x = 50R, with the inflow distance L i = 10R; however, since the ideal flow consists of an isolated torus, we have verified that the domain finiteness did not significantly affect the results through the blockage effect.In more detail, the cases at Γ = 1 for Re = 100 and Re = 200 have also been run using a second domain of dimensions L x = 90 and L r = 30, and the number of nodes augmented accordingly: the differences in the drag coefficient now below 0.15%.In a second test it has been verified that the torus was placed far enough from the inflow to have results independent of L i : this was checked by repeating the runs at Γ = 1 and Re = 100 and Re = 200 with L i = 20R and on a grid extended ahead of the torus; in this case, the maximum difference of the drag coefficient was below 2%, which is acceptable for our purposes.
Owing to the large CPU requirements, only a few three-dimensional cases were run, and all of them on the standard domain with a grid 258 × 401 × 65 in the radial, axial, and azimuthal directions.These last runs were mainly aimed at showing the onset of three-dimensional instabilities rather than analyzing the flow details.
A first validation of the results is given by the analytical study of [4], who computed the Stokes flow around a "closed" torus (Γ = 2) for which they found a drag coefficient C D 5.611/Re; the present numerical results yield C D = 5.63 at Re = 1 and C D = 56.04 at Re = 0.1, which are in very good agreement with the prediction.
A further validation is provided in Figure 3 showing the drag coefficient (i.e., the force in the opposite direction with respect to the mean flow normalized by the factor ρU 2 S/2, with S the blockage area) as function of the Reynolds number.The results are compared with the classical paper by [13], and the deviations range from a minimum of 0.53% at Re = 100 up to a maximum of 1.50% at Re = 200.It is worth mentioning that the validation in the range 100 ≤ Re ≤ 1000 is much more challenging than that at low Reynolds, owing to the steeper velocity gradients in the flow and the thinner boundary layers at the sphere surface.Further validation is shown in the next section for different aspect-ratios Γ and Reynolds number ranges.Here it suffices to mention that the results of Figure 3  • present results, filled data from [13].

Results
In order to make the discussion easier, we will briefly illustrate first the flow evolution for a fixed aspect ratio (Γ = 1) as a function of the Reynolds number.For Re ≤ O(1), the dynamics is dominated by viscosity with the flow smooth and steady, without separations; this is the Stokes regime (Figure 4a), where between the upstream and downstream flows there is an approximate symmetry that becomes more exact the more Re → 0. As Re is increased, the flow remains steady and a first recirculation can appear: however, this recirculation is not attached to the torus, but located downstream in the wake.This regime is generated by the competing effects of blockage of the torus and downstream convection of the mean flow.Owing to the vorticity distribution within the bubble, it will be referred to as Hill regime, and it can exist only in a narrow range of aspect ratios (around Γ = 1) and Reynolds numbers (Figure 4b).A further increase of Re produces separations attached to the torus that in some cases can coexist with the Hill recirculation (Figure 4c), but eventually only recirculations attached to the torus surface can survive (Figure 4d).When the Reynolds number is still augmented, the flow becomes unsteady (Figure 4e) with a wake that oscillates in time and sheds compact vortices.The last regime is observed for even higher Re when the flow becomes three-dimensional and also the azimuthal symmetry is lost (Figure 4f).
The above scenario can change considerably depending on the aspect ratio of the torus: in the following sections, we will illustrate the flow dynamics in the various regimes for different values of Γ.

Stokes Regime
The first regime is that of very low Reynolds numbers characterized by dominating viscous effects that determine the flow.In this context, the governing dynamics is a balance between pressure and viscous forces that is expressed by the Stokes equation ∇p = µ∇ 2 u; this does not have inertial terms, and therefore cannot involve the fluid density ρ.On the other hand, the drag force is commonly written as F D = ρU 2 SC D /2, with S a reference surface that-for this problem-is the blockage area S = 4πRa.The only way for F D to get rid of ρ is that C D ∼ 1/Re = µ/(ρUR), and this is a common feature of all very low Re flows.Figure 5 shows that indeed also for the torus the drag coefficient follows this power law; here, in the interest of brevity, we show only the results for Γ = 0.1 and Γ = 1 even if the scaling C D ∼ 1/Re is observed for all values of Γ.For small aspect ratios Γ, the flow tends to that around a two-dimensional cylinder (R/a → ∞) and, accordingly, already for Γ = 0.1 the C D is almost indistinguishable from that for a circular cylinder.It is worth noting that the Reynolds number for the cylinder Re d is defined with the cylinder diameter, while here we use the toroidal radius R; this implies that the two Reynolds numbers are related through Re d = ReΓ, which should be kept in mind when comparing the two flows.
The behaviour is qualitatively similar for the aspect ratio Γ = 1, although the fits have different coefficients than for the two-dimensional cylinder: here we have found in the Stokes regime C D 10/Re.We report for comparison the data for the same flow computed in [9], which show a perfect agreement with the present ones.

Hill Regime
For Γ = O(1), some values of the Reynolds number gave rise to a special steady state which, to the knowledge of the authors, has never been reported before.More in detail, at Γ = 1 the flow at 10 ≤ Re ≤ 70 showed the appearance of a recirculation bubble in the wake of the torus, the bubble consisting of a vortex ring (Figure 6a,b).The shape of this ring is approximately spherical at its onset, but it becomes more elongated in the streamwise direction for increasing Re owing to its interaction with the strain field produced by the torus.In fact, this recirculation shows strong similarities with the Hill spherical vortex [10] that is characterized by a linear azimuthal vorticity distribution ω θ = Ar, with r the radial coordinate and A a constant.In view of this similarity, we have verified whether this relation also holds for the recirculation of the torus, and two representative results are shown in Figure 7.It is clear that at Re = 10, when the vortex is approximately spherical and not strained by flow, the linear relationship is attained very closely (note that in Figure 7 the axes do not have the same range in both panels).On the other hand, at Re = 50 the bubble is significantly elongated in the streamwise direction (Figure 4c) and partially stripped by the external strain field; accordingly, the ω θ (r) relation shows evident deviations from linearity-especially in the upstream part of the vortex that is more subject to the wake of the torus.As a successive step, we analyze how the recirculation is generated: it is expected that the appearance and disappearance of the vortex is somehow related to the relative proportions of fluid passing over and through the torus.As measure of the flux, we have computed the Stokes streamfunction ψ(r, x) = r 0 σu x (σ, x)dσ that, computed for x = L i and r = R − a, yields the flow passing through the torus (inner flux).Some results are reported in Figure 8, showing an increase of the flux with Re and a tendency to saturation at "high" Reynolds number.It is also evident that this flux decreases for increasing Γ, since the torus tends to be "fat" and the surface of the central hole diminishes.Figure 6d clearly shows that when this flux is vanishing the bubble attaches to the torus as for the separation behind a bluff body.On the other hand, for low values of Γ, the inner flux is too intense to allow for an axial steady separation.Therefore, it turns out that this regime can only exist for some combinations of Re and Γ to have an inner flux strong enough to produce a recirculation (through the adverse pressure gradient downstream of the torus) but not too strong to wash it away.The Re-Γ values for which this special steady state was found are contained within the region H of Figure 1.

Steady Regime
Apart from the limited range of Γ where the Hill regime can exist, after the Stokes region of the Re-Γ phase space we find a steady regime in which stationary separation bubbles appear downstream attached to the torus.It has been observed that for Γ = 2, already at Re = 1 a small recirculation is present while for Γ = 0.1 only above Re = 100 does a separation appear.In the first case, the flow is characterized by a single bubble that grows for increasing Re, and it remains steady up to Re = 75 (Figure 9d).As Γ is decreased, two counter-rotating recirculations develop with the inner one (that closer to the axis of symmetry) smaller than the external one (Figure 9c).This difference decreases as Γ is diminished, and for Γ → 0 any difference between the bubbles should disappear since the two-dimensional cylinder case is recovered.Our simulations at Γ = 0.1 show that indeed the flow approaches the cylinder limit even if differences are still present.In particular, Figure 9a reports the length of the separation bubble (measured as the distance from the back of the torus where the streamwise velocity reverts its sign as in Figure 9b) as function of Re compared with the analogous quantity for the plane cylinder.It is observed that the results are very close even if the data for the torus are systematically below those for the cylinder.After having verified that this effect is not related to insufficient numerical resolution or to the finiteness of the computational domain, we have found that this difference is caused by the finite aspect ratio Γ that produces a blockage of the inner flow.In other words, the flow going through the torus is accelerated more than the flow going over, thus pushing the separation towards the rear of the torus surface and consequently decreasing the length of the separation.This argument is in agreement with the recent paper by [14] that computed the flow around a circular cylinder over several domains, and they have observed that indeed the length of the steady separation bubble is reduced for increasing blockage of the flow.These results are qualitatively and quantitatively consistent with those of [7][8][9].

Unsteady Regime
A further increase of the Reynolds number enhances the inertial effects with respect to the viscous ones, and this induces flow transition from steady to unsteady.The recirculations that in the steady regime were attached to the torus are now alternatively shed in the wake and advected by the main flow.A typical case is shown in Figure 10 for the aspect ratio Γ = 0.1 and Re = 1000, and it is similar to the flow around a circular cylinder.
A characteristic quantity of the unsteady vortex street is the Strouhal number St = R f /U, where f is the main frequency of the flow.Here f has been measured from the first peak of the fast Fourier transform of the drag coefficient (Figure 11a), even if any other quantity (for example, velocity or pressure sampled in a point) yielded the same result.Similarly to the Reynolds number, the relation St d = StΓ has been used to compare the present results with those of a circular cylinder, since the Strouhal number for a cylinder St d is computed using the cylinder diameter.Figure 11b shows the normalized oscillation frequencies of the flow as a function of the Reynolds number, and in order to make the comparison easier, all the data are normalized in the same way as for the circular cylinder.For Γ = 0.1, the Strouhal number for the torus appears to be close to the corresponding value for the cylinder, yet systematically below it.Similarly to the separation length of Figure 9a, also in this case the difference is due to the finiteness of the aspect ratio and this is confirmed by the analogous quantities for Γ = 0.25 and Γ = 0.5 showing larger discrepancies and in the same direction as the case at Γ = 0.1.Another fundamental difference between the wake of a cylinder and that of a torus is that the vortices of the latter are in fact vortex rings, and they have a self-induced translation velocity that adds (algebraically) to the convection of the wake.For a vortex ring with thin core and uniform vorticity distribution within the core, the self-induced translation velocity is V = γ/(4πb)[ln(8b/ ) − 1/4], where γ is the core circulation, b the ring radius, and the core radius.Looking at Figure 10, we can see that the vortices shed from the external side rotate clockwise while the others rotate counterclockwise; this implies that the internal vortices must travel downstream faster than the external ones, and accordingly the wake tends to diverge radially as shown by the straight lines connecting the vortex centres.This behaviour is different from that of the cylinder, whose wake spreads symmetrically with respect to the streamwise direction.
In the case at Γ = 0.1, the internal and external vortices have very similar dimensions, and the measured differences of γ were below 2%.However, the situation becomes different as Γ increases because the external vortices increase in size at the expense of the inner vortices that eventually disappear.This completely changes the instability mechanisms, and tends to inhibit the flow unsteadiness.Figure 11b shows that indeed for a given Re the oscillation is slower for increasing Γ, and we have found that for Γ > 1 the axisymmetric unsteadiness is suppressed with a direct transition from steady axisymmetric to three-dimensional flows (Figure 1).

Three-Dimensional Regime
Beyond the unsteady regime the flow loses the axial symmetry and it becomes completely three-dimensional; examples of flow snapshots are given in Figures 4f and 12 for three different aspect ratios.We have seen in the previous section that the flow features are strongly influenced by Γ, and this is also confirmed in the three-dimensional case where the flow instability depends on the topological structure of the separation.In order to clearly distinguish between truly axisymmetric and marginally three-dimensional flows, we have checked all the force and moment components.In the former case, all coefficients except for the drag assumed zero values (the computer round-off error), while for three-dimensional flows-though small and with zero time average-every coefficient was different from zero.Further flow analysis would require the investigation of the azimuthal instabilities of the ring vortices in the wake and their interaction.However, this kind of analysis is quite complex; it would deserve a dedicated study that is beyond the purposes of this paper.Here we only aimed at finding the region of the three-dimensional flow dynamics in the phase diagram of Figure 1 .

Conclusions
This paper has analyzed the flow around a toroidal ring at low Re with the aim of proposing a new benchmark that is more complex than the classical circular cylinder and sphere yet sufficiently simple and well-defined to allow a detailed investigation.The Re-Γ phase diagram reported in Figure 1 shows that several different regimes are possible, and they span from the Stokes dynamics up to the fully three-dimensional flow.For Γ → 0, the flow tends to become that around a circular cylinder, while as the torus aspect ratio increases the differences between internal and external flows become dominant.For Γ = 2, the inner flow is inhibited, only a single separation behind the torus is generated, and the flow resembles that around an obstacle.A particular steady flow referred to as Hill regime has been found for a limited range of Γ and Re in which a recirculation downstream in the wake (but not attached to the torus) is present.
Before concluding this paper we wish to point out that the boundaries of the phase diagram of Figure 1 must be intended only as indicative transitional zones rather than sharp limits.For example, looking at Figure 6a,b it appears that increasing Γ has the effect of moving the recirculation upstream towards the torus until it attaches to its surface; for intermediate Γ, the distinction between a simple massive steady separation from a highly strained Hill vortex is certainly questionable.In addition, the cases delimiting the unsteady regions are considerably more computationally expensive than others, since distinguishing between very slowly decaying disturbances from a marginally sustained perturbation requires an integration over several viscous time scales that-especially for high-Re flows computed on fine grids-becomes infeasible.For example, a typical simulation in the Stokes or steady viscous regime can be run on a single processor desktop computer within 2-4 CPU h, while in the unsteady regime, in order to collect enough statistics, the CPU time increases up to 20-30 h.If a full three-dimensional unsteady case is considered, the CPU time ramps up to 300-500 h and parallel computing is needed.

Figure 1 .
Figure 1.Phase diagram of the regimes for the flow around a torus: • indicate axisymmetric runs, • three-dimensional simulations.V: Stokes regime; S: steady regime; H: Hill regime; U: unsteady regime; 3D: three-dimensional regime.Note that the runs at Re < 0.1 are not reported in the figure to avoid the squeezing of the regions to the right.
radial boundary at r = L r , ∂u ∂t + c ∂u ∂x = 0 for convective outflow at x = L x u = 0 noslip condition at the torus surface.

Figure 2 .
Figure 2. Sketch of the problem.

Figure 3 .
Figure 3. Drag coefficient as function of the Reynolds number for the flow around a sphere (Γ → ∞):• present results, filled data from[13].

Figure 10 .Figure 11 .
Figure 10.Instantaneous snapshot of azimuthal vorticity ( for negative, for positive values) and Stokes streamfunction ( ) for the flow around a torus at Γ = 0.1 and Re = 1000.The dashed lines () connect the centres of the vortices of the same sign.