Determination of drag coefficients in automatic ball balancers at low Reynolds numbers

The precise calculation of drag forces in the technical application of automatic balancing of rotating machinery provides important information about efficiency and stability. With increasing geometric complexity of the design, this poses a challenge that can be solved with computer-aided fluid dynamic approaches. In rolling element bearings, the influence of drag induced by the lubricant is predominantly considered in the context of efficiency loss estimations, whereas the movement of the rolling elements is mainly constrained by the contact with the bearing rings and, if present, the cage. Automatic ball balancers, which are installed in rotating machinery to reduce unbalance excitation, are in design very similar to fully lubricated ball bearings missing the cage, the inner ring and the majority of the balls. Inherent to the functional principle, the balancing efficiency and stability are significantly influenced by the choice of lubricant and resulting drag forces. Therefore, the estimation of the drag coefficient based on the geometry and lubrication of automatic ball balancers plays an important role in the engineering process. With a focus on the Stokes flow regime, the drag coefficient for a single sphere in an annular flow domain is determined numerically with finite volume discretization and the SIMPLE steady state solution scheme. Based on a parameter study utilizing the presented solution approach, a simple empirical relation between the design of the automatic ball balancer and resulting drag coefficients is derived. As a result, a drag force formulation based on the balancer geometry and the lubrication fluid properties is presented, which helps to supplement a large number of published kinetic models regarding the analysis of automatic ball balancer stability and transient behavior, giving a better understanding of the influences of design decisions regarding geometry and lubricant.


Introduction
The calculation of drag forces resulting from flows around spheres and cylinders are widely used benchmarks to evaluate computational methods in the area of fluid mechanics due to the simplicity of the geometry, see for example Johnson and Patel (1999), Jones and Clarke (2008), and Bayraktar et al. (2012). When defining a large computational domain compared to the dimensions of the obstacle, the effect of domain boundaries on the solution is attenuated in order to represent the targeted free flow condition. In a further study by Martinez et al. (2009) the single plane contact of a sphere and the resulting influence of the Blasius boundary layer was examined and in the context of process engineering advances in numerical studies on particle flows in pipes are presented by Krishnan and Kaman (2010) and Chen et al. (2019) for example.
In the study at hand, however, the drag of a sphere confined in an annular domain is examined in the context CONTACT Lars Spannan lars.spannan@ovgu.de of automatic balancing of rotating machinery. The concept of automatic balancing describes the arrangement of balancing masses, usually balls, in relation to the axis of rotation of a rotor in order to compensate its unbalance. Depending on the application, the unbalance can vary during operation, for example laundry in washing machines, samples in laboratory centrifuges or abrasion induced imbalance in angle grinders. Since the friction between the balls and the rotor cavity in which they are orbiting around the rotor axis has negative influence on the accuracy of the unbalance compensation, the outer rings of ball bearings are often used in the design to decrease friction as presented by Lindell (1996) for the application in angle grinders. Therefore, the basic design of the automatic balancer can be described as a radial ball bearing missing the inner ring, the cage and the majority of the balls, see Figure 1. In addition to the balls, the cavity is filled with a viscous fluid, which has influence on the dynamics and noise Figure 1. Principle design of automatic ball balancers. When the outer ring of a ball bearing is inserted to function as a raceway with reduced friction, a raceway groove is present which is not shown in this sketch. emission of the automatic balancing unit. As a result, the movement of the balls in the annulus cavity is significantly affected by drag forces, which are to be determined by numerical analysis in the presented study.
Since the first publication by Thearle and Schenectady (1932) on the topic of automatic ball balancers, several research advances have been published on the analysis of the dynamic behavior. For example, Sharp (1975) presents equations of motion, with which the stability of the systems balanced state can be analyzed. More detailed analyses followed by Nesterenko (1984), considering anisotropic supports of the rotor, Van De Wouw et al. (2005), discussing the influence of dry friction on the residual unbalance, Ishida et al. (2012), presenting methods to improve the balancing efficiency and Wright and Peng (2019), examining the influence of varying angular velocities.
Commonly, the viscous forces induced by the fluid on the balancing balls inside the cavity are modeled in linear relation to the orbit velocity of the balls relative to the cavity walls. Even though this enables qualitative studies on the influence of increased viscosity of the used fluid, no methods are proposed to quantitatively determine a suitable linear drag parameter based on the balancer design and the fluid parameters. In the study at hand, a CFD approach is presented delivering the missing modeling aspect in order to describe the viscous forces based on the fluid viscosity and the geometric properties of the annulus and the balancing balls.
The determination of drag on spheres considering contact with domain boundaries is examined in several research publications. Goldman et al. (1967), for example, studied the slow motion of a sphere in contact with a plane wall using asymptotic solutions of the Stokes equations. Jan and Chen (1997) conducted experiments with a sphere rolling down an inclined plane and derived an empirical relation between the drag coefficient c d and the corresponding Reynolds number Re. Miyamura et al. (1981) derived wall correction factors for spheres in contact with cylinders and parallel plates experimentally. Numerical investigations on drag in rolling element bearings were carried out by Marchesse et al. (2014) focusing on the tandem effect of multiple spheres in a row. However, the authors are not aware of any publications presenting quantitative models to examine the viscous drag the scope of automatic balancers.
The simulations presented in the following sections yield the drag coefficient of a single sphere in contact with the geometry of an outer radial ball bearing ring confined in an annular domain. The flow conditions are limited to low Reynolds numbers (Re < 10) and therefore, as will be presented in the results, lead to a steady flow regime in which the product of Reynolds number and drag coefficient equals a constant, reinforcing the assumption of laminar flow. With the transition to unsteady behavior for spheres in contact with a flat surface near Re = 250, see Martinez et al. (2009), sufficient distance to unsteady behavior is assumed for the flow in the curved annulus, considering the restriction to Re < 10.
Section 2 describes the configuration of the simulations carried out with the OpenFOAM 1 software. The resulting drag coefficients for different Reynolds numbers as well as a subsequent parameter study, which analyses the influence of the ball orbit radius, the ball diameter and the cavity dimensions are presented in Section 3. With the chosen parameter range, the design space of commercially available regular radial ball bearings is covered. Section 4 specifies how the identification of the reciprocal relationship between the drag coefficient and Reynolds number in laminar flow regimes helps with the parametrisation of a variety of kinetic models of automatic ball balancers published in the literature. A concluding discussion and outlook is given in Section 5.

Methods
Since the Stokes flow regime at low Reynolds numbers is considered, a symmetric flow regarding the axial geometric symmetry plane is assumed, see Figure 2. The inlet and outlet of the curved domain are defined in a suitable distance to the sphere spanning an orbital angle of 60 • so that the assumption of a periodic flow condition can be implemented while reducing the computational effort. An increase of the orbital angle to 90 • in the calculation of the reference balancer design, Section 3, showed no significant influence on the result, substantiating the proposed approach.
Since the sphere has a single contact point with the raceway of the outer ring, the geometry of the domain is altered by introducing a cylindrical column connecting the sphere and the raceway, see Figure 2. The column diameter is set to 16% relative to the sphere diameter. While the impact on the results is believed to be marginal, this modification improves the meshing and solution process. The marginal influence is attributed to the fact that the height of the cylinder is relatively small, even with a diameter of 16% of the ball diameter, due to the very similar radii of the ball and groove. On the basis of several modeling variants, 16% turned out to be the smallest minimum ratio at which a stable numerical solution could be obtained.
In order to reduce the computational effort, the case of a sphere orbiting inside the annular domain with an orbital angular velocity of , an orbit radius of R orbit and a ball radius of r ball , see Figure 3, is viewed regarding a coordinate system originating in the center of the sphere with the X-axis pointing outwardly in radial direction and the Y-axis pointing in axial direction. Due to The velocity is defined by the product of the ball angular velocity ω about the Y-axis going through the sphere center and the distance to it. ω results from the ideal rolling condition (R orbit + r ball ) = ωr ball .
The OpenFOAM simpleFoam application uses the semi-implicit method for pressure-linked equations (SIMPLE) for a steady-state solution of the continuity equation and the momentum equation for incompressible flow, where u, p k , σ , S u describe the velocity vector, the kinematic pressure, the stress tensor and the momentum source, respectively (Ferziger et al., 2020). Herein, the OpenFOAM incompressible solvers transform the static pressure p s to kinematic pressure p k , i.e. p k = p s /ρ, with ρ being the fluid density. The equations for velocity and kinematic pressure are solved iteratively. As a result, for each mesh element three degrees of freedom for the velocity and one degree of freedom for the pressure are determined. The sought drag coefficient c d is then calculated from the sum of forces in the tangential direction (Z) on the ball surface patches with face area vectors s, the facial reference area A ref and the free stream velocity v ref at the inlet, i.e.
Due to the exploitation of symmetry, the reference area is defined as A ref = πr 2 ball /2. In the conducted calculations, the iterations are stopped when the residuals fall below 2 × 10 −6 . A structured hexaedral mesh is generated with the OpenFOAM blockMesh utility and different mesh refinement levels are realized, see Table 1. The mesh refinement level (Mlvl) scales the number of cell divisions per block in the computational domain. The mesh is defined for half of the computational domain, see Figure 4, and mirrored on the XY-plane. Since the pressure field on the sphere surface is evaluated in order to gain the drag coefficient, a small element size is defined in the adjacent cells. The sphere is surrounded by a radial padding of mesh regions which are discretized gradually with the smallest element height at the sphere surface. The remaining coarser mesh regions are also discretized gradually towards the sphere, resulting in improved results without an increase in degrees of freedom.

Results
The domain geometry for the study is determined by the properties given in Figure 3 and Table 2 representing an automatic ball balancer prototype available to the authors. Following the solution method described in  Section 2 the velocity and pressure fields are obtained and the drag coefficient according to Equation (3) is returned by OpenFOAM. An exemplary normalized result of the flow field in the symmetry plane is shown in Figure 5. The mesh refinement study in Figure 6 shows that the resulting drag coefficient for a Reynolds number Re = 1 does not change significantly for refinement levels above four. The relative deviation with respect to the result of level six is less than 0.4%, which is considered sufficient for the calculation of drag forces presented in Section 4. Therefore, all subsequent simulations are performed with refinement level four. For different Reynolds numbers in the range of Re = 0.01 · · · 100 the drag coefficients are determined and  are shown in Figure 7 along with data from the literature for the free flow condition published by Morrison (2013), the flat surface contact published by Jan and Chen (1997) and experimental results gained for the balancer prototype in the prevalent double logarithmic diagram. In the previous work of Spannan et al. (2017) the utilized experimental method, based on the equilibrium evaluation between drag forces and gravitational forces, is presented. In order to derive the drag coefficient c d for different Reynolds numbers the balancing device prototype is driven on a horizontally aligned axis and the balance of forces is examined. For the considered range of low Reynolds numbers, measurements with silicone oil of kinematic viscosities ν 1 = 50 mm 2 s −1 , ν 2 = 298 mm 2 s −1 and ν 3 = 715 mm 2 s −1 are used as reference for the conducted CFD simulations. For a detailed description of the experimental procedures the reader is kindly advised to the given reference. Figure 7 shows an increased drag coefficient for the considered range of Reynolds numbers compared to the results of the flat surface contact. This result is expected due to the additional confinement of the flow influenced by the domain geometry. Since only the Stokes flow regime is considered, where the product of drag coefficient and Reynolds number result in a constant, the results are transformed to Figure 8 showing the product c d Re for Reynolds numbers up to Re = 100. The simulation results present a constant value of c d Re = 387 and agree well with the experimental data gained for Reynolds numbers Re = 0.1 · · · 0.7 and the oil with a kinematic viscosity ν 3 = 715 mm 2 s −1 . The deviations  of the experimental results acquired with the oil of a kinematic viscosity ν 2 = 298 mm 2 s −1 are believed to be caused by friction induced temperature fluctuations during operation and associated viscosity variations. In addition, the agreement between the numerical and the experimental results for Reynolds numbers around Re = 100 is given so that it can be assumed that the solution approach is also applicable for the transition area after the Stokes flow regime. Since the overall agreement is acceptable, the presented simulation model is considered sufficient to gain useful results for a subsequent parameter study.

Parameter study
Although the presented numerical model can supplement the design process of automatic ball balancers by determining drag coefficients for a given geometry at low Reynolds numbers, the possibility of reducing the calculation time and effort by deriving an empirical relation based on a parameter study is investigated in the following. As described previously, the utilization of radial ball bearing components in automatic ball balancers is advantageous due to low friction. With commonly available radial ball bearings in mind, the variance of the drag coefficient for orbit radii in the range R orbit = 5 mm · · · 300 mm and relative ball sizes r ball /R orbit = 0.05 · · · 0.3 are considered in the parameter study. The annulus height compared to the ball radius varies only slightly between different bearing sizes and is therefore held constant with H = 4r ball . Similarly, also the groove depth ratio W 2 /r ball , see Figure 3, is believed to vary only slightly and is held constant. The distance between the ball center and the inner boundary of the fluid domain may vary for different designs of automatic ball balancer units. With regard to compact designs, low values for the clearance W 1 /r ball to the inner annulus boundary are to be expected. By increasing the distance, it was checked at which point the influence of the inner boundary becomes negligible.
Based on the evaluation of the dimensionless units Reynolds number and drag coefficient, it is apparent that the change in the orbit radius R orbit has no influence on the result, if the remaining proportion quantities in Table 2 remain constant. The check can be carried out quickly by using the converged velocity field of one solution and scale it inversely to the orbit radius magnification. From this initial field solution the convergence is achieved promptly.
Increasing the relative ball radius r ball /R orbit results in an increase of the product of Reynolds number and drag coefficient, see Figure 9. The calculations are conducted with increasing relative ball size while using the same number of cells.The study of increasing distance between the inner boundary of the annular domain and the sphere, Figure 10, shows a decay in drag, which can be approximated by an exponential relation with an asymptotic value for large radial clearances. Clearances W 1 /r ball > 2.5 result in values c d Re < 322, which was identified by Jan and Chen (1997) for the single flat surface contact. The decreased drag for designs with high clearances to the inner wall is attributed to the groove in the outer ring. However, since minimal overall ABB dimensions are aspired, W 1 /r ball is usially less than two. While increasing the domain volume, the number of cell divisions in radial direction are increased accordingly. Furthermore, the mapFields application of OpenFOAM is utilized to copy a previous solution to the coincident part of the domain. Starting with the solution for W 1 /r ball = 5 and consecutively decreasing the domain volume, the collective computational effort is reduced.
Based on the presented simulation results, neglecting interdependence of the two variables, the relation  with a = 248, b = −1.375, c = 245 and d = 813, is proposed in order to approximate the resulting drag in the Stokes regime of a single ball in an automatic ball balancer. The coefficients a, b and c are fitted with trust-region adjusted non-linear least-squares and the coefficients d and e with linear least-squares methods, respectively. Since Equation (4) emerges from the parameter study based on the reference geometry given in Table 2 the best agreement is expected near W 1 /r ball = 1.554 and r ball /R orbit = 0.08.

Application
One of the most researched kinetic models for the application of automatic ball balancers is that of a planar oscillator, Figure 11, with x, y, ϕ rotor , ϕ ball being the two translatoric degrees of freedom, the rotation of the rotor and the orbital position of the ball, respectively. In the Figure 11. Example of a planar oscillator with an automatic ball balancer.
derivation of the equations of motion the viscous damping forces due to the fluid inside the cavity is commonly modeled in a linear manner, e.g. a drag moment proportional to the difference of the rotor angular velocityφ rotor and the ball orbit velocityφ ball , is introduced. See for example Sharp (1975) and Sperling et al. (2004).
In several works the importance of the viscous forces is recognized and simulations are conducted in order to identify values for the coefficient β, which lead to optimal dynamic behavior of the rotor system. See for example Ryzhik et al. (2003) for the analysis of rotor run ups and Bykov and Kovachev (2018) for the identification of stable operating speeds of automatic balancing devices. However, concrete correlations between the introduced constant viscous coefficient, β in Equation (5), and the design of the balancing unit as well as fluid density and viscosity are not discussed. The results of the parameter study described above, yielding a relation between the annulus geometry, the ball size, the fluid properties and the constant c d Re, can be used to remedy this situation as specified in the following.
Considering the familiar equation for drag forces, the discussed drag moment can be written in the form with ρ, ν, A ref , v ref being the fluid density, the fluid kinematic viscosity, the cross-sectional area of the ball and the reference flow velocity, respectively. The Reynolds number is introduced as Under the assumption that v ref = R orbit andφ ref ≈ ϕ rotor −φ ball , meaning that the fluid is rotating with the rotor angular velocity, the coefficient β can be identified as a constant For run up processes of machinery equipped with automatic ball balancing devices the assumption of low Reynolds numbers does not hold generally and velocity differences between the rotor and the fluid occur. But when steady state operations of automatic ball balancers are considered, the stated assumptions are met. As a result, stability analyses, which incorporate Equation (5), see for example Bykov and Kovachev (2018), can be conducted solely based on geometric properties of the balancer, fluid properties and the presented identification of the product of the drag coefficient and Reynolds number, see Equation (4).

Discussion
In the presented study a finite volume approximation solution for the laminar flow around a sphere in an annular domain is presented. The domain geometry is derived from an automatic ball balancer prototype and a parameter study is conducted covering the dimensions of common radial ball bearings. The drag simulation results, in terms of the product of Reynolds number and drag coefficient for a reference geometry, show the anticipated constant numerical value for Reynolds numbers below ten. The increase in drag compared to the case of a sphere in contact with a flat wall discussed in the literature shows the influence of the side walls and the curvature of the domain. It is shown that the knowledge of the product of Reynolds number and drag coefficient is beneficial for the parametrisation of kinetic models for automatic ball balancers. However, due to the dependance of the results on the fluid viscosity the calculation of drag forces in automatic ball balancer models is subject to uncertainties in determining the viscosity of the fluid and possible temperature sensitivity.
In contrast to applications of automatic ball balancers utilizing larger numbers of balancing balls in order to increase the maximum compensable imbalance of the rotor, the presented flow simulations consider only one ball neglecting possible tandem effects. Two balls in tandem arrangement on the orbit will reduce the drag forces on the downstream body (Marchesse et al., 2014) making it worth considering appropriate correction factors to Equation (4).
During run-up of rotors equipped with automatic ball balancers the angular acceleration may be sufficient to leave the Stokes flow regime. As a result, the product of Reynolds number and drag coefficient deviates from a constant value. With regard to computational flow simulations, the assumptions of symmetry and laminar flow will not be met, leading to increased effort to gain results for larger Reynolds numbers.
Concluding from these limitations, future advances in numerical solution of the fluid flow in automatic ball balancers should withdraw the assumption of steady symmetric flow to allow for the consideration of turbulent flow at higher Reynolds numbers with multiple spheres in tandem orientation inside the annular cavity. Consequently, more advanced numerical solution methods must be used.
In addition to the viscous forces, the friction forces and possible geometrical deviations of the raceway are influencing the dynamical behavior of the automatic ball balancing system, see Chao et al. (2005), Green et al. (2006) and Bykov and Kovachev (2019) for example. Therefore, a fully coupled fluid structure interaction of automatic balancers, i.e. in multi body simulation models, should only be pursued if the modeling quality of friction and imprecise geometry is also enhanced. Compared to the increased effort resulting from a fully coupled simulation the presented method based on the widely utilized gas-particle flow approach (Henderson, 1976) represents a good compromise in the view of the authors.