Effects of Degrees of Motion Freedom on Free-Fall of A Sphere in Fluid

Free-fall of a sphere in fluid is investigated at a Galileo number of 204 by direct numerical simulations (DNS). We mainly focus on the effects of different degrees-of-freedom (DOFs) of the sphere motion during free-fall. The characteristics of free-fall are compared with those of flow past a fixed sphere. Additional numerical tests are conducted with constraints placed on the translational or rotational DOFs of the sphere motion to analyze different DOFs of sphere motion. The transverse motion contributes significantly to the characteristics of free-fall; it results in the retardation of the vortex shedding, leading to the decrease of the Strouhal number. In addition, the transversal sphere motion exhibits the tendency to promote the sphere rotation. On the contrary, the effects of the sphere rotation and vertical oscillations during free-fall are negligible.


Introduction
The free-fall of particles in fluid is commonly seen in nature and industrial applications. Typical examples are sedimentation of pollutants in the ocean and vertical transportation of ores in deep-sea mining. As the fluid-particle interaction is the basis of the fluid-solid two-phase flow, the studies on freely falling particles are of great importance from the viewpoint of both theory and engineering application.
Generally, the spherical particles in fluid are most widely investigated. The free-fall of a sphere in viscous fluid has been investigated experimentally since the 1900s (Allen, 1900). In the early years, the investigations mainly focused on the drag coefficients of the spheres. Based on a large amount of experimental data, the relations of the drag coefficient and Reynolds numbers were summarized, and typical empirical formulas have been proposed for the applications of engineering (Clift and Gauvin, 1970;Haider and Levenspiel, 1989;Terfous et al., 2013). Subsequently, the Basset force was investigated by Mordant and Pinton (2000) through a series of experiments and the influences of the Reynolds number and sphere density were summarized. Focusing on the path instabilities of freely moving spheres, Jenny et al. (2004) summarized a parametric space (G is the Galileo number and defined as G = , and are the sphere and fluid densit- ies, D is the sphere diameter, g is the gravitational acceleration and is the kinematic viscosity of fluid.) for path regimes and identified five types of path regimes: vertical, steady oblique, oscillating oblique, zigzagging periodic, and chaotic. Similar study was conducted by Horowitz and Williamson (2010). A parametric map for the sphere responses and wake modes was presented. Additionally, the flow patterns have been analyzed to explain for the sphere responses based on numerical simulations (Feng et al., 1994;Rahmani and Wachs, 2014;Uhlmann and Dušek, 2014;Yu et al., 2004;Zhou and Dušek, 2015) and optical measuring techniques (Horowitz and Williamson, 2010;Veldhuis et al., 2005;Veldhuis and Biesheuvel, 2007).
Owing to the strong coupling between the sphere motion and hydrodynamic forces, the characteristics of the free-fall are so complicated that it is challenging to establish a general relation between the sphere motion and hydrodynamic forces. In order to evaluate the effects of different degrees-of-freedom (DOFs) of the motion on the hydrodynamic forces and Strouhal number during free-fall, Namkoong et al. (2008) conducted the constrained numerical simulations on the sedimentation of a two-dimensional (2D) cylinder in fluid by constraining the translational or rotational DOFs of the cylinder motion. The transverse motion was observed to contribute significantly to free-fall and play a crucial role in reducing the Strouhal number. Such a method was adopted to investigate the effects of the Magnus effect owing to the sphere rotation on the lift force by Yu et al. (2004).
To the best of our knowledge, the study on the effects of translational and rotational DOFs of the sphere motion on the free-fall is rare. However, it will help reveal physics of free-fall and is meaningful to the further study on the particle-fluid two-phase flow. The potential application of the work will provide the guidance on dealing with pollutant sedimentation in the ocean and the design of the conveying pipelines for the particle-fluid two-phase flow. In the present study, we investigate the free-fall of a sphere in viscous fluid using the direct numerical simulation (DNS) method. The overset mesh and a moving computational domain are used to simulate the long-distance sphere motion. With the work of Namkoong et al. (2008), the effects of different DOFs of the sphere motion on free-fall are evaluated by constrained numerical simulations. We focus on a simple case in which the sphere moves with the oscillating oblique trajectory at and G=204. The corresponding Reynolds number based on the terminal velocity and the sphere diameter is Re=285. The choice of the parameters is based on the consideration that differences of the sphere responses and the forces can be distinguished easily when some DOFs of the sphere motion are constrained in this case.
The remainder of this paper is organized as follows. The mathematical formula and numerical methods are presented in Section 2. The computational overview and validation studies are provided in Section 3. The results and discussion are described in Section 4. Finally, the concluding remarks are summarized in Section 5.

Mathematical formula
The incompressible flow is governed by the Navier-Stokes (N-S) equations expressed as follows: ∇ · u = 0; (1) ∂u ∂t where u is the velocity of the fluid, p represents the pressure and f g denotes the gravitational force, and t is the time. The N-S equations are discretized using the finite volume method (FVM). The discretization of each term is formulated by integrating the term over a control volume using Gauss's theorem. Volume and surface are then linearized using an appropriate scheme. Here, the second-order upwind scheme is used to discretize the convective term; the implicit unsteady solver is applied to carry out the transient computation and the second-order temporal discretization is utilized to obtain more accurate results compared with the first-order temporal discretization. For each time step, the SIMPLE (Semi-Implicit Method for Pressure Linked Equa-tion) scheme is adopted to solve the governing equations for the velocity component and pressure-correction equation.
The 6-DOF motion solver is employed to simulate the free-fall of the sphere in response to the flow-induced pressure and shear stress as well as the gravity (CD-adapco, 2016). The resultant forces acting on the sphere are calculated first, and used to solve the governing equations of the sphere motion as follows: A dω dt where m represents the mass of the sphere, A is the tensor of the moments of inertia, V and ω are the translational and angular velocities of the sphere, F and M are the resultant forces and moment, respectively.

Computational overview
The free-fall of a sphere with the diameter of D is considered here. The density ratio between the sphere and fluid is and the Galileo number is G=204. Fig. 1 shows the configurations of the simulations in the present study. A cylindrical computational domain with the diameter of 10D is adopted and the effects of the boundaries are verified to be negligible. The local coordinate system (x, y, z) is fixed to the earth with the z-direction opposite to the direction of gravity. A second coordinate system (X, Y, Z) translates with the sphere; its origin is located at the center of the sphere, and the directions of the axes are the same as those of (x, y, z). In addition, a spherical coordinate system is employed with its origin attached to the center of the sphere, as shown in Fig. 1b.
Rather than constructing a long computational domain, the moving computational domain is adopted to obtain a long-term simulation of the free-fall (Asao et al., 2013;Rahmani and Wachs, 2014;Wang et al., 2014). The computa- tional domain is set translating to follow the sphere in the present study. The translating velocity of the computational domain is the same as that of the sphere; the rotation of the computational domain is constrained. Therefore, the position of the sphere related to the computational domain is constant. The bottom and side of the computational domain are set as the inlet boundaries: the flow velocity u x = u y = u z = 0 is prescribed for the quiescent fluid, and the pressure is set as a zero normal gradient. At the top of the domain (the outlet boundary), the velocity is specified as a zero normal gradient boundary condition, and the pressure is set as zero. During the simulation, the translating velocity of the domain will not be superposed on the flow velocity at each volume cell. Thus, the moving computational domain can be used to represent the unbounded and undisturbed fluid. On the surface of the sphere, a no-slip and impermeable boundary condition is represented.
The overset mesh scheme is applied to addressing the rotation of the sphere. As shown in Fig. 2, two sets of meshes are generated: a cylindrical background mesh enclosing the entire domain and a spherical overset mesh containing the moving sphere. The spherical overset mesh is designed with a diameter of 1.8D and the center of the overset mesh coincides with the sphere center. Then, the overset mesh is inserted into the background mesh with the outer boundary defined as the overset condition to separate the background mesh into an active part (the blue mesh in Fig. 2a) and an inactive part (covered by the red overset mesh). The governing equation would not be solved in the inactive part of the background mesh and the flow will be computed in the overset mesh and the active part of the background mesh, similar to a continuous mesh. However, the flow solutions in the overlapping cells (the green cells in Fig. 2a) are obtained from the interpolation of both the overset and background meshes, rather than solving the governing equations. More detailed descriptions of the overset mesh scheme are available in the literature (CD-adapco, 2016;Hadžić, 2006;Koblitz et al., 2017;Romero-Gomez and Richmond, 2016).
A cutaway view of the grid structures is shown in Fig. 3. The grids are initially generated with a base size of 0.1D; refinements are applied in the zones containing the sphere and the vortex to capture the details of the wake structures. In addition, the prism layer cells (see CD-adapco (2016)) are utilized in the vicinity of the sphere surface; this provides a better description of the viscous boundary layer and captures the wake structures accurately at modest computational costs (Koblitz et al., 2017). In the present study, the convective Courant number is no more than 0.5 in most areas of the computational domain.
The non-dimensional forms are used to represent the calculated results. The scales used to non-dimensionalize the results are as follows: L 0 =D for length, V 0 = for velocity, for time, for angular velocity, for force, and for moment. The pressure coefficient is calculated by , where is the pressure of the ambient fluid, and V z denotes the vertical component of the terminal sphere velocity. The drag force coefficient is defined as , where F z denotes the hydrodynamics force in the z-direction.

−
In the present study, the symbols and ' are utilized to represent the time-average value and the amplitude, respectively, of the quantity.

Validation study
Validations of the spatial and temporal resolutions are conducted with four sets of unstructured meshes and three time-step sizes at G=204. Once being released, the sphere accelerates rapidly and falls with an oscillating oblique path trajectory finally. The results in the fully developed states, including the mean vertical and transverse velocities ( and ), mean drag coefficient , oscillating frequency and amplitude of the transverse velocity, and oblique angle between the sphere path and z-direction (see Fig. 5), are summarized in Table 1 and Table 2.
As shown in Table 1, favorable agreement is observed of the results in Cases A2, A3, and A4; marginal difference of the results exists between Cases A3 and A4. With reference to the time-step size, the results in Cases B2(A3) and B3 match adequately, as shown in Table 2. Therefore, based on the validations above, the spatial resolution in Case A3 and the time-step size of are sufficient to simulate the free-fall of a sphere.
To verify the capability of our numerical method, we compare the vertical velocity and drag coefficient with the experimental data of Mordant and Pinton (2000) and numerical results (Rahmani and Wachs, 2014), as shown in Table 3. The calculated and agree well with the results in literature at G=49, 177 and 393; the oblique angle matches well at G=177. Subsequently, an adequate comparison of the sphere motion at G=190 are presented in Table 4, where the marginal differences are observed between the results of Uhlmann and Dušek (2014) and the present simulation. Moreover, the drag coefficients at 150<Re<300 are compared with the results summarized in the literature (Clift and Gauvin, 1970;Terfous et al., 2013;Haider and Levenspiel, 1989), as shown in Fig. 4. A reasonable convergence is observed between the present results and the data in the literature. Based on the validations above, we conclude that the present numerical method is valid to simulate the free-fall of a sphere in fluid.

Results and discussion
4.1 Characteristics of free-fall α The trajectory of the sphere in a viscous fluid at G=204 is shown in Fig. 5. The sphere travels in an oblique direction and exhibits slight transverse oscillations (see Fig. 6a and Fig. 6c). The mean oblique angle from the z-direction is 5.32°. The translational velocities of the sphere presented in Fig. 6a show that the sphere accelerates rapidly once being released and attains the periodic oscillating state finally.     The phase space plot of V x and V y is provided in Fig. 6b, from which we can easily observe that the sphere moves in a single plane. Fig. 6c and Fig. 6d show the transverse (V t ) and angular ( ) velocities of the sphere in this plane. The frequency of the sphere motion is f=0.0723; the Strouhal number defined as St=f/V z corresponds to St=0.0518, which is close to the value of 0.054 predicted by Jenny et al. (2004). For the sake of a clear description, a new coordinate system (X t , X n , Z) is defined with a rotation of the local coordinate system (X, Y, Z) around the Z-direction to position the X axis in the motion plane (i.e., the X t axis). The flow patterns (Fig. 7-Fig. 10), and the hydrodynamic forces (F t and M t , see Fig. 11) are discussed according to the coordinate system (X t , X n , Z).
The instantaneous wake structures at the maximal lift force F tmax in the stable state are presented in Fig. 7. Here, the vortical structures are identified by the iso-surface of Q-criterion proposed by Hunt et al. (1988): where S and Ω denote the strain and rotation tensor, respectively. The Q-criterion has been widely used to represent the three-dimensional vortical structures (Tian et al., 2017a(Tian et al., , 2017b. As shown in Fig. 7, a vortex has just been shed from the previous cycle, and the head of the nascent hairpin is just beginning to form near the sphere and extend downstream. The symmetrical wake structures explain the inplane motion of the sphere. The periodic vortex shedding can be regarded as the "pinch-off" of the wake structure, resulting in the periodic lift force F t and transverse velocity V t of the sphere. The instantaneous distributions of vertical flow velocity u z in the (X t , Z) plane in the cases of free-fall and flow past a fixed sphere are presented in Fig. 8a and Fig. 8b, respectively, at the maximal lift force F tmax . The sway of the wake   flow from the vertical axis during free-fall is larger than that in the fixed case. The wake flow behind the freely moving sphere oscillates around the oblique direction marked by the blue dashed line, whereas the wake behind the fixed sphere oscillates around the vertical axis. It is inferred preliminarily from the wavy shapes of the wake that the frequency of the vortex shedding in the fixed sphere case is higher than that in the freely falling case.
The pressure coefficients C p on the sphere surface (in the (X t , Z) plane) are shown in Fig. 9a; the transverse projection of the resultant pressure coefficient along the vertical axis are revealed in Fig. 9b. Based on the computed data, the sketches of C p and on the freely falling and fixed spheres are plotted in Fig. 10a and Fig. 10b, respect-∆C px θ ively. As changes its direction at ≈-12.5° (see Fig. 9b) during the free-fall, the horizontal projection of the pressure on both sides of the falling sphere cancels out, resulting in a relatively small lift force. As revealed by the sketches in Fig. 10b, the resultant transverse pressure on the fixed sphere along the vertical axis always points to the same direction, leading to a larger lift force.  Fig. 11 shows the hydrodynamic forces on the sphere. The lift force F t on the fixed sphere oscillates around a nonzero mean value . Once the sphere is released, the lift force introduces the transverse sphere migration: contributes to the inclination of the sphere path and oblique angle; whereas the fluctuation (F t -) introduces the oscillation of the transverse velocity. Similarly, the moment M t ∆C px θ Fig. 9. Distribution of (a) C p and (b) on freely falling sphere (represented by FF in the legend) and fixed sphere (represented by Fixed in the legend) in the (X t , Z) plane. The horizontal axis is defined in Fig. 1b. on the fixed sphere introduces the angular velocity shown in Fig. 6d: the non-zero mean value contributes to whereas the fluctuation (M t -) leads to the periodic variation of .
Generally, the frequency of the hydrodynamic forces on a freely falling sphere is lower than that of a fixed sphere; this is related to the lower vortex shedding frequency behind the freely falling sphere. During the free-fall, the flow rolls up on the low-pressure side of the sphere after a previous vortex shedding cycle; a vortex is gradually generated and then shed. The transverse sphere migration acts as a suction flow on the low-pressure side of the sphere, so that the pressure recovers with the separation point retarded, leading to the lower vortex shedding frequency during the free-fall.

Influence of different DOFs of the motion
To have an in-depth understanding of the free-fall, the 6-DOFs of sphere motion are decomposed by constraining the rotational or translational motion of the sphere to evaluate their contributions to the free-fall. Such a method has been adopted to investigate the free-fall of a 2D cylinder (Namkoong et al., 2008).
Numerical tests were conducted by restricting different DOFs of the sphere motion at G=204: (1) the case with the constraints on sphere rotation and transverse migration (Case NN), (2) the case with the constraint on the sphere rotation (Case NR), and (3) the case with the constraints on the transverse migration (Case NT). The results are compared with those of the free-fall (Case FF) and flow past a fixed  First, the marginal oscillations of V z are observed in Case NN, as shown in Fig. 13a. This is owing to the periodic vertex shedding behind the sphere. The stable oscillations slightly reduce the amplitudes of both F t and M t when compared with those in flow past a fixed sphere, as shown in Fig. 14a and Fig. 14c. As the vertical velocity of the sphere oscillates periodically, the relative velocity between the sphere and the surrounding fluid changes slightly, leading to the marginal variation of the pressure and shear stress on the sphere surface, as revealed in Fig. 15b. Nevertheless, In the next step, we focus on the effect of the sphere rotation by comparing the results in Cases NT and NN. As shown in Fig. 13c, the sphere attains a periodic angular velocity in the stable state in Case NT. Owing to the Magnus force which is proportional to , the lift force is reduced slightly (see Fig. 14a). The rotation increases the flow velocity on one side of the sphere and reduces it on the other side. This changes the pressure distribution on the sphere surface; the horizontal projection of the pressure difference is reduced marginally, as shown in Fig. 15a and Fig. 15b. The moment M t impels the sphere to rotate once the constraints on the rotational DOFs of the sphere motion are removed; the shear stress on the sphere surface decreases rapidly. In the stable state, M t varies periodically around zero (Fig. 14c). Similar effect of the sphere rotation on M t is observed by the comparison of Cases FF and NR. When the sphere rotation is constrained (Case NR), M t varies periodically, however, with a non-zero time-averaged value and a larger amplitude. The limitation of the sphere rotation increases the shear stress on the sphere surface, resulting in a larger moment. ω t Finally, we evaluate the effect of the transversal sphere migration, which contributes significantly to the free-fall. When comparing Cases FF and NT, the hydrodynamic forces show entirely different characteristics (see Table 5 and Fig. 14); the Strouhal number in Case FF (St=0.0518) is smaller than that in Case NT (St=0.1284, see Table 5 and Fig. 12). The pressure on the sphere changes significantly, as shown in Fig. 15. Once the constraints on the transverse motion (in Case NT) is removed, the sphere is impelled to move horizontally; the lift force F t decreases sharply. As described in Section 4.1, a similar effect of the suction flow will act on the low-pressure side, which makes the pressure recover with the vortex separation retarded, leading to a lower vortex shedding frequency and a larger amplitude of F t . In addition, in Cases NT and FF, it is observed that the amplitudes of and M t in Case FF are larger than those in  Case NT, which means that the transverse sphere migration promotes the sphere rotation, as shown in Fig. 13c.

Concluding remarks
Based on the engineering concerns of pollutant sedimentation in the ocean and ore transportation in deep-sea mining, the free-fall of a sphere in viscous fluid is simulated by using the DNS method. The Galileo number is 204, and the density ratio between the sphere and fluid is 2. The dynamic responses of the sphere and flow patterns during the free-fall are analyzed. The hydrodynamic forces on the freely falling sphere are compared with the results of flow past a fixed sphere. Moreover, the effects of the translational and rotational DOFs of the sphere motion on the free-fall are evaluated by constraining different DOFs of the sphere motion. Based on the comparisons among the results of different cases, the main conclusions are drawn as follows.
(1) Smaller Strouhal number is observed during the freefall when compared with that of flow past a fixed sphere. Once the fixed sphere is released, the mean values of the lift force and moment contribute to the mean oblique path direction and mean angular velocity, respectively, of the sphere. The fluctuations of the hydrodynamic forces result in the oscillations of the path and rotation of the sphere.
(2) The transverse sphere motion plays a crucial role in the free-fall. Similar to the effect of a suction flow, the transverse sphere motion makes the pressure recover with the vortex separation retarded on the low-pressure side of the sphere, leading to a lower vortex shedding frequency. In addition, the transverse motion promotes the development of the sphere rotation during the free-fall.
(3) The effects of the rotation and vertical oscillation of the sphere are negligible. The Magnus effect owing to the sphere rotation and translation exhibits the tendency to reduce the lift force. However, the effect is marginal.