Numerical effectiveness investigation of the automatic ball balancer with a deep chamber

The article describes the simulation results of an unbalanced rotary machine with an automatic ball balancer with a deep chamber (ABB-DC) based on a mathematical model of such a system. The ABB chamber has a cylindrical shape and is filled with balls that can move freely in general motion. The model takes into account mutual collisions between the balls, as well as between the balls and the chamber. The model also includes friction forces and rolling resistance. The comparison of simulation results performed with two different construction types of the chamber, classic single-layer drum (ABB-SC) and ABB-DC, confirmed the better efficiency of the ABB-DC in the optimal range of the number of balls. The machine body movement was limited to planar motion. The correctness of such a limitation has been verified by analyzing waveforms of reaction forces of constraints and their moments. Waveforms were divided into transient and steady states. The results of the analysis indicate the additional constraints do not significantly affect machine body movements. However, this might affect the ball arranging process in the ABB chamber. The paper also presents details of performing numerical simulations based on the developed mathematical model and numerical tests to determine the optimal integration time-step.


Introduction
An automatic ball balancer (ABB) is a device used to automatically balance an unbalanced rotor. It consists of a chamber and a set of balls that move freely inside the chamber mounted on the same axis as the rotor. During the operation of an unbalanced rotary machine, the balls roll around the chamber surface and automatically take an appropriate position to reduce undesirable unbalance using the self-synchronization phenomenon, firstly described by Blekhman [1]. This phenomenon is commonly used in vibratory systems, e.g., to synchronize two [2,3], or even four [4] eccentric rotors (ERs).
The key role in the synchronization process of balls and unbalanced rotor is taken by the tangent components of the inertia forces that result from the chamber vibration. The balls are moving around the rotating raceway until they reach a position where the mass center of the balls and rotor coincides with the axis rotor. When the balls finally take the appropriate position, the rotor becomes balanced, and the chamber vibration successively disappears due to the existing dissipation. Then, only the normal components of the inertia forces remain; thus, the balls do not change positions relative to the chamber. Over the past years, there have been many publications devoted to ABBs-the first being written in 1932 by Thearle [5]. Among these are experimental papers [6][7][8][9][10] and theoretical ones [11][12][13][14]. Due to the complex system of the rotary machine with a ball balancer, the latter is usually based on simplified mathematical models. Many simplifications often include a reduction in the machine body to a material point [15], the absence of resistance in the bearings, and neglecting the collisions between the balls [16][17][18]. Most of the publications concern ABB systems with one layer of S. Pakuła [19 and 20] concerning numerical effectiveness investigation of ABBs with a drum chamber and a ring chamber proved that a higher efficiency was achieved for the first type of chamber where balls were free to move in a radial direction. Therefore, the author assumed that removing the additional constraints imposed on the balls can positively influence the effectiveness of an ABB system. Thus, the author proposed a new type of chamber named a deep chamber (DC). The balls inside such a chamber can freely move and rotate in either direction and this allows them to distribute in many possible ways in many layers. The mathematical model of the system of the machine with the ER and ABB-DC is partially based on the model from paper [21] that concerns the machine with ER and ABB-SC. The equations of motion have been supplemented with equations of motions in the rotor axis direction and with equations of rotational motion relative to two axes perpendicular to the rotor axis. Due to the constraints imposed on the machine body, some of the equations became those of equilibrium. The author for convenience decides to apply a vector notation of differential equations of general motion. The mathematical model used Hertz's model to calculate the contact forces of the mutual collision of the balls as well as the balls with the ABB's chamber. It is described in Sect. 3. The methods for performing the numerical simulations and determining the optimal integration time-step as well as the initial balls distribution algorithm are presented in detail in the Appendices. Simulation results of ABB-DC basing on the developed mathematical model are concluded and compared with the ABB-SC in Sect. 4. The body movement has been restricted to planar motion by the use of the imposed constraints. The reaction forces following from those constraints are the imperfection measure of the model simplifications. Section 4.3 is devoted to investigating a result of imposing additional constraints based on the reaction forces waveforms.

Model
The physical model of a machine with ER and ABB assumes that the machine body is rigid and supported on a set of elastic-damping systems. The ABB chamber is fixed to the rotor axis. The chamber has a cylindrical shape that can contain a few layers of balls along its axis direction. The rotor is driven by an asynchronous motor, which is rigidly mounted in the machine body. The chamber contains a set of freely moving balls. The balls inside collide with each other as well as with the ABB chamber's surface. Figure 1 presents a graphical model of such a system.
Two systems of coordinates have been used to describe the system motion ( Fig. 2): a stationary system, and a system fixed to the movable machine body. The center of the stationary coordinate system (x, y, z) is located in such a manner that the axes' directions coincide with the main axes of inertia of the rotary machine. Besides, it is assumed that the x-axis direction coincides with the rotor axis. The axes of the coordinate system fixed to the body and the ABB is designated by coordinates (ξ, η, ζ ). The axes of this system coincide with the main axes of the body inertia.
It has been assumed that the rotor is rigid and statically unbalanced. The center of mass of ER is located in the vertical symmetry plane of the body machine zy. As a result, the movement of the machine body and ABB is reduced to planar motion (ω ζ ω η ẋ 0). In such a case, the systems of axes (η,ζ ) and (x,y) are always lying on the same plane. To determine the equations of motion, the system is decomposed into subsystems. Then, the differential equations of motion are formulated for the body (2), ABB (7), and each ball (9), respectively.

Machine body
The physical model of the rotary machine body is presented in Fig. 3. This shows the basic geometrical dimensions used in the mathematical model and the vectors of the forces and moments of the reaction forces from the ABB.
The differential equations of body motion are presented as two vector equations (1) and six scalar equations: where: K T ,K ϕ -stiffness matrix (translational and rotational, respectively). B T ,B ϕ -damping matrix (translational and rotational, respectively). P 1 ,P 2 -rolling bearing reaction forces. R x -reaction of forced constraint (only in direction of axis x). r s -vector connecting elasticity center with rotor axis. where: x k , y k , z k -coordinates of the center of body mass relative to the fixed coordinate system (x, y, z).
x 0 , y 0 , z 0 -coordinates of rotor axis relative to the fixed coordinate system (x, y, z). ω x , ω η , ω ζ -rotational velocity of body-machine relative to x, η, and ζ axes, respectively. B-depth of ABB chamber and distance between bearings. P y , P z -components of sum of bearing reactive forces relative to x-and y-axis, respectively: M o -moment of bearing rolling resistance forces expressed by Eq. (4), M-torque of drive according to the Kloss model (5): where: M zn -rated torque of electric motor. p-overload capacity. ω s -synchronous speed. ω u -critical speed. d cz -shaft neck diameter. μ t -reduced coefficient of rolling bearing friction.
The mathematical model from [21] has been supplemented with the action of two rolling bearings-one on each side of the ABB's wall. Due to the arbitrary position of the balls in the ABB's chamber, the values of the rolling bearing reaction (P 1 ,P 2 ) will be different. As a result of such differences, the rotary machine body will tend to revolve relative to the ζ and η central inertia axes. However, those rotations are blocked by the additional forced constraints imposed on the model, which are represented by the action of reaction moment (M R ). The vector M R is parallel to the plane yz, and its value is a measure of the asymmetric distribution of balls in the ABB's chamber relative to the plane yz. The freely moving balls inside the chamber can collide with the front and rear chamber wall, causing a reaction of Bearing 2 in the x-direction. Consequently, an additional constraint is imposed that prevents bearing movement towards the x-axis. The reaction force of this constraint is designated asR x . If the reaction forces and the moment of reaction forces of constraints (M R , R x ) imposed on the body motion are close to zero, this means that the planar model of the machine body and ABB is accurate enough. This is verified in Sect. 4.3.

Automatic ball balancer
A graphical model of the ABB chamber is presented in Fig. 4. The figure also shows the yz plane in detail along with the vector components of the reaction forces P 1 and P 2 .
The ABB's chamber motion equations are presented as two vector equations: As the ABB moves in planar motion, the system of vector equations (6) can easily be presented as six scalar equations: where: x w , y w , z w -coordinates of center of rotor mass relative to fixed coordinate system (x,y,z). F i0 , T i0 -contact forces and friction forces caused by contact of chamber and ith ball. M z -reduced moment of reaction forces of balls M 0i (41) and rolling bearings M P (44) according to The motion equations of the rotary machine and ABB do not include the gravity forces, as it was assumed that the origin of the stationary coordinate system (x,y,z) coincides with the rotor axis in the state of static equilibrium. With these assumptions, the gravity forces are reduced by the pre-tensioning of the springs. The two final equations from Eqs. (7) were equated to zero, because the ABB's rotational motion relative to its main inertia axes η w and ζ w is blocked by the pair of rolling bearing reaction forces M P .

Balls
A ball is a solid body with spherical symmetry; hence, it is possible to write its rotational motion in a coordinate system that is always parallel to the global coordinate system, i.e., (x i , y i , z i ) with the origin in the ball's center mass. The general motion equations of the ith ball are the following system of equations: The system of equations (9) is presented in vector form for convenience: The values designated with subscript "i0" refer to the action of the chamber, and subscript "ij" refers to the action of the jth ball. Section 3 describes methods of determining contact forces F, friction forces T, moments of friction forces M, and moment of rolling resistance M t as a result of the contact of ith ball with the jth ball (Sect. 3.1) and the ABB's chamber (Sect. 3.2).
As the balls are uniform and are central symmetric bodies, a few advantages can be pointed out: • the central moment of inertia has the same value relative to each axis crossing the center of mass, • the scalar value of moment inertia J i in the second of Eqs. (10) can be used instead of a tensor J i , • Euler's angles can be omitted, as the angles of rotation are not a point of interest.

Equations of motion for the system
Assuming that an ABB is fixed to the rigid axis of a rotor and the body rotation angles are small, the constraints concerning the machine body and ABB are formulated according to (11): Additionally, there are three constraints imposed on the body machine that prevent it from moving in the x-direction and from rotating relative to the y and z axes: After substituting the constraints (11) into the systems of equations (2) and (7), we obtain six independent motion equations: The remaining six equations are equilibrium equations that we will use to calculate the values of the reaction forces and moments: After removing the reaction forces components P y and P z from the system of equations (13), we obtain the system of equations in the generalized coordinates of q 0 {y 0 , z 0 , ϕ x , ϕ 0 }. Together with six coordinates q i {x i , y i , z i , ϕ xi , ϕ yi , ϕ zi } for each of n balls, the system respecting the imposed constraints has s 4 + 6n degrees of freedom. The system of the s nonlinear differential equations of second-order was solved by the use of numerical methods described in Appendix 1.

Calculation of contact forces
This section presents the method of calculating the contact forces and friction forces as well as the moments of these forces. These forces arise as a result of balls colliding with each other as well as with the ABB's walls. These interactions are described separately for the pair of balls in Sect. 3.1 and the ball with the chamber walls in Sect. 3.2.

Contact forces between balls
There are a huge number of collisions between the balls inside an ABB's chamber during the normal operation of the system of a rotary machine and ABB. The most intensive of these occur during the starting of the machine. The collisions are elasto-plastic; part of the collision energy is absorbed and converted to heat and plastic deformation. The measure of the energy dissipation is described by the coefficient of restitution R, S. Pakuła which generally depends on the collision intensity. In this paper, a constant value of R 0.54 was assumed. The contact force in the function of the strain was calculated using a modified Hertz model [22] according to Formula (15). The strain depends on the relative position of the balls as shown in Fig. 5.
where δ 2R k − r i j is the strain resulting from contact, with R k the ball radius. The position of the ith ball relative to the jth will be calculated according to the following formula: The relative velocity of the balls' centers is expressed by the following equation: According to Newton's third law of motion, the mutual contact forces between the balls have equal values but opposite directions; hence, F ji − F i j . Consequently, the values of the contact forces between the balls can be calculated only once during the process of the numerical calculations. The values of the friction forces accompanying the balls' collisions are calculated according to the Coulomb model (18). The direction of the friction force vector will depend on the direction of the relative contact point velocity projected on the contact plane. The method of determining the relative contact point velocity is presented below.
Friction force: wherev sτ i j is the versor of the component of relative contact point velocity projected on the contact plane τ (see Fig. 6).
The relative contact point velocity is determined as the difference between the contact point velocity of the ith ball and the contact point velocity of the jth ball: The contact point velocity is calculated according to the principle of the superposition of the ball's center motion and the contact point motion relative to the ball's center: where v τ i , v τ j are the contact point relative velocities resulting from ball rotation around its center. These are calculated according to Finally, after substituting Eq. (20) and (21) into Eq. (19), the contact point velocity will be expressed by Eq. (22) and graphically presented in Fig. 7.
The vector components of the contact point's relative velocity towards the collision direction n (the normal direction of the contact plane) and along the direction parallel to contact plane τ will be calculated according to the following formulas, respectively: Knowing the friction vector T ij from (18), the vector of the friction force moment can be calculated according to Newton's third law of motion also applies to the friction forces; hence, T ij −T ji . In such a case, the vector of the friction force moment for the jth ball is calculated according to the following equation: The next value to be determined is the vector of rolling resistance moment M t ij, which is calculated according to whereω t i j is the versor of rolling of the ith ball on the jth ball. The vector ω t ij determines the rolling direction and is calculated based on the vector of contact point velocity s relative to the ith ball center v τ i (21) and the relative velocity component of the ball's center (17) projected on the plane tangent to the ball at the contact point v τ ij (contact plane) according to the following equation: It is worth mentioning that Formula (27) is not used to calculate the relative velocity of the rolling balls but only to determine the rolling direction. The tangent component of vector v ij is obtained by subtracting the normal component from it: Finally, by substituting the value from Equations (21) through (27), we obtain The relationship between the moments of rolling resistance of the ith and jth balls is as follows:

Contact forces between balls and ABB chamber
This subsection describes the contact forces of the balls on the ABB's chamber as well as the method of calculating them. Description pages are split into two parts (describing the balls' contact with the race and with the flat wall of the chamber, respectively). Similar to the collisions between the balls, a modified Hertz model (36)

General designation
i,ĵ,k-versors of the global stationary coordinate system (x,y,z). | r | r 2 x + r 2 y + r 2 z -vector length r . r r | r | -versor of vector r . r 0 , r i -position vectors of rotor axis and ith ball center relative to the stationary coordinate system (x,y,z). e 0 ·î + e · cos ϕ oĵ + e · sin ϕ 0k -unbalance vector. r i0 r i − r 0 -ball location relative to the geometrical center of the ABB's chamber. v i0 v i − v 0 -ith ball velocity relative to the geometrical center of the ABB's chamber.

Relative contact point tangent velocity with ball rotation
Friction force T i0 All values described above refer to the action of the ABB's chamber on the ith ball. The contact and friction forces of the balls acting on the chamber will satisfy the following equations: The moment of the contact and friction forces acting on the ABB's chamber as a result of contact with the ith ball is calculated according to The manner of the radius r ie calculation depends on the chamber wall which a ball is contacted to as shown in

Reaction forces of rolling bearings
The elimination of the components of rolling bearing reaction forces P y and P z from the system of equations (13) is possible under the assumption that the moment resistance value of rolling bearing M o (4) is known, but it depends on the values of forces P y and P z . However, M o can be approximately estimated based on the reaction forces from the previous integration step in the process of the numerical calculations. The values of force components P y and P z are calculated based on the first two equations from System (13) supplemented with the vector of the ABB's chamber gravity vector m k g: Equations (13) do not include the body machine and ABB weights, because the coordinate system was assumed in the state of the static equilibrium of the machine and an ABB with no balls. However, the ABB weight must be taken into account when calculating the bearing reaction forces. Due to the acceleration present in Eq. (43), such an approach of calculating the bearing reaction forces requires us to assume some values of the rotor axis acceleration at the initial time (ÿ 0 ,z 0 ). When the numerical simulations based on the presented mathematical model only concern the start-up of the machine and the ABB, then it can be assumed that the acceleration will be zero without any significant error. Calculating the forces in successive stages will require knowledge of the accelerations from a previous simulation step. In numerical calculation processes, one should then remember that, in addition to the state vector X, the calculated values of acceleration (ÿ 0 ,z 0 ) should be maintained.
Equations (43) are used to calculate the sum of the components of the reaction forces of both bearings presented in the system. The two final equations from System (7) were used to estimate the share of each bearing in these reaction forces. These equations indicate that the reduced moment M z expressed by Eq. (8) is perpendicular to the plane ζ 0 and η 0 (hence, also to the plane y 0 ,z 0) . To estimate the moment M z , the moment of bearing reaction forces M P (44) must be known. It can be calculated from Eqs. (45). The forces and radius vectors from that equation are presented on the top view of the cross section of the chamber in Fig. 8.
where r O B1 B 2 ·î. Figure 3 shows the vectors of bearing reaction forces P 1 and P 2 , and distance vectors r B1 and r B2 from the mass center of the rotor to the bearings. Considering both Equations (44) and (45), the moment vector of bearing reaction force M P can be expressed as follows: Knowing the components of the vector M P , the components of the vector M z are calculated according to Eq. (8). These components are presented in (47): where P y P 1y + P 2y , P z P 1z + P 2z .
As M z is perpendicular to the plane y 0 z 0 , the two final components of Vector (47) are equal to 0. Along with the first equation from System (7) (whereẍ 0 0), these equations are the three independent equilibrium equations: Using Eqs. (48) and (49), we calculate the values of the components of the reaction forces of both bearings. These are expressed by Eqs. (50) and (51), respectively:

Result of simulations
This section is devoted to the results of the simulation experiments conducted using the mathematical model presented in the previous sections. The details of the model implementation and simulation process are described in Appendix 1. The optimal integration step was selected as a result of the numerical test described in detail in Appendix 2. The initial distribution of balls inside the chamber for all simulations was performed according to the algorithm presented in Appendix 3. In the first step after the simulations were done, a visual verification of the results was performed. To do this, a special visualization application was created in the MATLAB® environment. It allows the user to see the motion of a large number of balls inside the ABB chamber. Next, the simulation results were compared to those performed by using the model from paper [19]. The last subsection is devoted to investigating the effects of the imposed constraints that make the machine body motion planar.

Model verification
The very first verification of the created mathematical model was based on an observation of the balls inside the chamber during the start-up of the machine. The numerical data from the calculation process was used to create an animation of the working machine. Three significant stages of the machine-starting processes can be distinguished. Figure 9 shows the position of the balls before the start as well as during the transient and steady states. The pictures from Fig. 9 show the simulation process in three stages. Before starting the machine, the balls take the initial position inside the chamber (a). After the start-up, the transient process begins, and the balls start moving (b) to new positions to balance the rotor (c). In the steady state, it is noticeable that most of the balls are located on the opposite side of the mass center of the unbalanced rotor. This indicates that the ABB  is working correctly. The balls inside the chamber occupy the correct position using the self-synchronization phenomenon [23,24], where the control signals here are the vibrations of the rotor. Simply, the balls will move inside the chamber until the vibrations stop, whereas vibrations stop due to the dissipative forces when the rotor is balanced statically.
As long as the chamber movements in the rotor axis direction and rotation with respect to two axes perpendicular to the rotor axis are blocked by additional constraints, the x-components of vibrational forces acting on balls don't appear. These forces provide a feedback signal which is necessary to arrange the balls in the correct position along the rotor axis for dynamic balancing. For this reason, imposing additional constraints excludes the possibility of achieving dynamic balance. In this case, however, it does not affect machine vibrations in the xy-plane, but the reaction forces of these constraints. Releasing one of the constraints could make achieving dynamic equilibrium possible but it is not a subject of this paper and demands a mathematical model to be supplemented with at least one additional equation of motion that refers to the released constraint. The next section contains a comparison of the effectiveness of ABB-DC with the classic ABB from paper [19], which concerns a static balance.

Results comparison
After each simulation experiment, the waveforms of the state vector coordinates were saved in matrix form and archived as a binary file. The state vector contains the coordinates of the machine body, each ball and their derivatives (see Table 2 from Appendix 1). From these waveforms, the amplitude value of vertical vibration in steady state was calculated simultaneously and collected in a database. This parameter is a measure of the ABB's effectiveness (the small the amplitude, the better the efficiency). Next, the amplitudes for a number of balls and coefficients of friction were compared with the ones from paper [19], which concerns a model of a planar system of a rotary machine and ABB with a single-layer drum chamber (ABB-SC). The comparison was made for the same physical parameters of both systems and is depicted in Fig. 10.
The charts from Fig. 10 show the amplitude values of vertical vibrations in a steady state-after accelerating the rotor to nominal speed. The charts refer to ABB-SC (Fig. 10a) and ABB-DC (Fig. 5b), respectively. The lowest vibration amplitudes are observed for the ABB-DC; however, certain conditions must be satisfied. The most important factors affecting the work of the ABB are the number of balls and the coefficient of friction. Friction forces push the balls in the tangential direction to the raceway, which allows them to be put into circumferential movement. If there is no sufficient coefficient of friction, the balls get stuck on the bottom of the ABB chamber during the running of the machine. Too small values of the friction forces were not able to Fig. 11 The distribution of three and nine balls inside the ABB-SC and ABB-DC, respectively bring the balls into circumferential motion around the race. This case is shown in Fig. 10a for μ k 0.3 and n 2, where the highest amplitude was observed. A similar phenomenon is observed for ABB-DC in Fig. 10b.
While the number of balls is equal or less than four, the balls inside the ABB-DC behave similarly to one ball inside the classic ABB-SC, as the balls distribute along the x-direction in only one row in such a case. The examples of the ABB chambers filled by three and nine balls are presented in Fig. 11.
The balls placed along the raceway are pushing each other and it helps them to get circumferential motion. Without that, the balls cannot take an appropriate steady position relative to the rotor and thus balance it. The ABB-DC needs at least four balls for μ k 0.7 and eight balls μ k 0.5 to let balls achieve a circumferential motion state. This state was not achieved for any number of balls while μ k 0.3. The ABB-DC only achieves a higher efficiency with a large number of balls except when the coefficient of friction is μ k 0.3.
As the number of balls in ABB-SC increases, the amplitude vibration increases as well (see Fig. 10a). All simulations were made for the same total mass of balls. That means the higher the number of balls, the smaller they are. In steady state, balls tend to be arranged in a single layer along the raceway. As the number of balls increases, the ball filling ratio of the raceway increases. It is shown in Fig. 11 illustrating ABB-SC with three and nine balls. In extreme cases, some of the balls can be situated on the side of the rotor unbalance, thus increasing its eccentricity instead of reducing it. In conclusion, increasing the number of balls decreases the total ability to compensate for the rotor unbalance. The mass optimization inside an ABB-SC chamber is a subject of paper [25]. As the simulation results show in the considered system, the ABB-SC chamber was too small to accommodate a large number of balls with a given total mass. Thanks to the deep chamber in ABB-DC capable to contain the balls in several layers, this problem does not exist.
The highest efficiency was obtained by the ABB-DC for n 5 and μ k 0.7. For such a configuration of parameters, the amplitude of the vertical vibrations is 6 μm. Figure 12 presents the average values of the vibration amplitudes in the steady state for the ABB-DC and ABB-SC for the optimal range of ball numbers ( [8,50] and [3,10]), respectively).
The ABB-DC achieves significantly higher efficiency for the optimal range of the number of balls than the classic ABB-SC does if a friction coefficient μ k ≥ 0.5 is provided. As the results of the simulation show, the average amplitude of the vibrations during the steady state can be reduced by 83% respective to that obtained by the classic ABB if a deep chamber is utilized instead. It is worth noting that the range of the optimal value of the number of balls for the ABB-DC is much broader than with the ABB-SC. A significant disadvantage of the ABB-DC is the need to provide a high coefficient of friction.
A certain imperfection in the presented model is a fixed restitution coefficient R, which in reality is dependent on the intensity of the collisions [26]. Even non-intensive (but numerous) collisions can dissipate a huge amount of energy when a constant value of R is assumed. This may cause the balls to be stuck to the bottom of the eliminator chamber when the coefficient of friction is insufficient. It is possible that, for a big number of balls n > 50, the simulation result will show the eliminator works correctly for the coefficient of friction μ k below 0.5. However, the compiler limitation on the memory allocation did not allow for such simulations to be performed.

Verification
The next stage is focused on the effects of the additional imposed constraints. Imposing additional constraints that restrict body movement to a planar motion can significantly affect the simulation of moving balls arranged inside an ABB's chamber. Therefore, the following reaction waveforms of constraints P x , M y , and M z were used to evaluate the developed mathematical model (Fig. 13). The waveforms were obtained from a rotary machine start-up simulation with five balls and a friction coefficient of μ 0.5.
Nonzero values of reaction forces prove that an unconstrained machine body will move in general motion. To examine the effects of these forces, simulations were conducted of an unconstrained body machine's motions under the imposed reaction forces according to Fig. 13. However, the reactions of these constraints were first analyzed in detail.
Each waveform can be divided into two stages: the transitional process, and the steady state. The first one is the process when the balls are arranged in an ABB's chamber. In the analyzed example, this lasts for about t I 8.6 s. The short impulses of forces are visible during this period as a result of the balls impacting the chamber walls. The second stage is the state when the balls have taken fixed positions relative to an ABB's chamber. This stage is characterized by irregular values of reaction forces within a range of -1500 N to 1500 N at an average value equal to 0.12 N and of moments of reaction force within a range of -350Nm to 350Nm) at an average value equal to 0.19 Nm. At this stage, the P x , M y , and M z values result from the contact forces and friction forces. These in turn depend on the ABB's angular velocity. The model assumed that the balls are in a continuous slippage state. When the balls already take their appropriate positions relative to the ABB's chamber, they perform slight oscillations around the position of balance as a result of the constantly changing friction-force impulses. They always act in the opposite direction of a ball's relative motion. In reality, the friction forces are not the cause of the motion. An analysis of the balls' transition from the slippage stage to the rolling stage without slippage requires many calculations. However, taking into account that the directions of these forces change at nearly each integration step and that their average value is close to zero, it is assumed that they would not affect the body motion to any significant degree. Hence, it has been decided that only those forces occurring in the transitional process will be taken into consideration. This means the values of the reaction force and the moments of reaction forces P x , M y , and M z at the steady state are equal to 0 in further investigations.
A simplified model of the body and ABB was used in the simulations as a system with three non-conjugated degrees of freedom: It was assumed in the simulations that the stiffness and damping coefficients in the horizontal direction are half of those in the vertical direction according to [21]: The remaining parameters are as follows: B 0.2m-body depth; 2L 2m-machine width; M 120kg-reduced system mass; J y 25kgm 2 and J z 20kgm 2 -moments of inertia.
The waveforms of the forces, the moments of the forces, and the system response are presented in Fig. 14.
The waveforms indicate that the maximum displacement of the unconstrained machine body along the xaxis is about 50 μm as a result of being hit by the balls. As a result of the reaction moments of the constraints, maximum rotation angles of 0.015°and 0.0028°were observed along the y and z axes, respectively.

Conclusions
Limiting the body movement to a planar motion in a 3-D model has significantly simplified the equations of motion. Forcing additional constraints that prevent body motion in the x-direction and rotation around the y and z axes has ensured that the main central axes of inertia of the body-ABB chamber system are constantly parallel to the axis of the stationary coordinate system (x,y,z). The consequences of forcing these restraints have been examined by using the reaction forces of these constraints as the input of an unconstrained test model of the machine body. The verification tests indicate that, for typical machines, the displacement along the x-axis is about a few dozen micrometers, and the angular deflection is in the range of a hundredth degree. Such small values are the result of the relatively small mass of the ball as compared to the combined mass of the machine and the ABB chamber. Therefore, it is estimated that limiting the body movement to a planar motion has not caused significant changes in the simulation of a 3-D machine's operation. Please note that there was no conjugation between the actions of the unconstrained body and the balls in the verification test, as a simplified model of a rotary machine and ABB was applied. The displacement waveforms of the unconstrained body might slightly differ from those in Fig. 14. However, the order of magnitude of these displacements will remain the same despite the significant simplification applied in the test model. While the interactions of balls have an inconsiderable influence on a machine's body motion, these small body displacements along the x-axis and slight rotations of the body around the z and y axes could have a significant influence on the balls' arranging process in an ABB's chamber.
The results of the tests from papers [19,20] indicate that the fewer constraints an ABB chamber has, the better the efficiency becomes. A similar conclusion can be made by observing the result of the simulation using the model presented in this paper. The automatic ball balancer with a deep chamber achieves over 80% better efficiency than the classic one in the optimal range of the number of balls; however, it does not work correctly with a low number of balls and low value of the coefficient of friction. Studies have shown that ABB-DC needs to be filled with a sufficient number of balls and the coefficient of friction shouldn't be smaller than 0.3. Although ABB-DC achieves higher performance, there are limitations to its application. Please note that installing ABB with a deep chamber requires sufficient space on the rotor. In many applications, i.e., fans or grinders, this may be a limitation that prevents the usage of this type of chamber. For this reason, in many cases, the ABB-DC installation requires a redesign of the whole rotor device. However, on the other hand, when it is impossible to increase the chamber diameter to provide a sufficient balance mass, a deep chamber with an unchanged diameter can be used instead. Finally, it is worth mentioning that the usage of a deep chamber may also contribute to dynamic balance, but this was not the subject of the study in this publication. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.  where: A-coefficient matrix of size 2s × 2s. X-state vector X X(q 0 , q 1 , ..., q n ,q 0 ,q 1 , ...,q n ) of size 2s × 1. Q-input vector of size 2s × 1.

Fig. 15 Graphical user interface
The indices of the individual coordinates of the machine and ABB as well as their X derivatives are presented in Table 2. The formulas defining the index number of the ith ball coordinate and its derivatives depending on the number of balls n are presented in Tables 3 and 4, respectively. These tables also present examples of index numbers for cases where one, two, or three balls are used in the simulations The first-order differential equations written in matrix form are convenient for the numerical calculations. To calculate the state vector derivative at each time step, the nonlinear equation (54) must be solved. After integration, we obtain a state vector for each step. The LU distribution method and Crout's technique were used to solve Eq. (54), while the integration was performed using the Runge-Kutta algorithm of the fourth-order  16 Testing of optimum integration step based on the experiment of the central collision of balls with a fixed step (see Appendix 2). The fixed-step was used due to the step function sign (42) used in the contact force equation (36). The first-order derivative of such a function is non-continuous; therefore, calculating the time step in each integration step is difficult and time-consuming where the sign function value changes its sign.
The results of the numerical calculations are shown as data tables. Each table is composed of rows that include simulation time t i and the state vector X(t i ) corresponding to that time. The majority of the data from the simulations is saved using a 1-ms time interval in external binary files, which is collected in a specially designed database. The file size mostly depends on the number of balls. The maximal number of balls is limited by the RAM that is allowed to be allocated by the compiler. A Pascal compiler allows us to create a program simulation with a maximum number of balls equal to 50. The file size varies from 60 MB for experiments with 3 balls to 1.2 GB for cases with 50 balls. The time consumption of the simulation calculation using a home-class computer varies from 17 min to 135 h, respectively. The calculation program was designed to work in a single thread.

Post-processing
At this stage, the data results obtained from the numerical calculations are analyzed and visualized. Depending on the needs, the time waveforms are used, i.a., to: • calculate the vibration amplitude in a steady state, • calculate the maximum vibration amplitude in a transitional process, • plot the amplitude graphs in the function of the numbers of balls or friction coefficients, • determine the time when the position of the ball inside the ABB's chamber was established.
The visualization intends to evaluate the model correctness and observe the behavior of the ball inside the ABB's chamber when the system is in operation. To conduct this experiment, a constant coefficient of restitution R k 1 was assumed. Under such conditions, the impact of the balls should be central and elastic. This means that the total kinetic energy of both balls should be the same before and after the collision. If the time step of the integration is too short, differences in the total energy value will be noticed. Consequently, the numerical error δ err related to the integration step will be expressed as the change of total kinetic energy after impact; this is calculated according to the following equation: where v 1 and v 2 -velocity after the collision; v 0 -initial velocity.
The results of the experiments for various relative collision velocities v w 2v 0 are presented in Table 5.
The experiments proved that numerical error δ err also depends on collision velocity v w . All of the observed errors were nonnegative, which means that energy was added to the system during the collisions of the balls and that the internal forces performed work in the system. This is contrary to the principle of energy conservation. During the ABB's operation in the steady state, the balls are in fixed positions relative to the unbalanced rotor and do not collide with each other. The most intensive collisions occur during the start-up. Numerous simulation experiments have proven that the maximum relative velocities of the balls during the start-up collisions do not exceed 4 m/s. Therefore, the optimum integration step was set at h 10 -6 .
All of the generalized coordinates and their derivatives are represented as double-precision floating-point numbers, which have a representation accuracy of 15-16 significant digits in the PASCAL language. Assuming the number representation error is δ L 10 -15 , the maximum representation-related simulation error is as follows (at a typical simulation lasting a few dozen seconds T 10 2 [s] and a chosen integration step of h 10 -6 s): Such an inconsiderable error means that the chosen integration step is right.

Appendix 3: Initial balls distribution algorithm
In the first stage of the simulation experiment, the initial condition for the system of the machine body and drum must be set as well as the initial condition for each of the balls. At the beginning of each simulation experiment, the machine body and ABB chamber stay at rest. This means that those general coordinates and their derivatives corresponding to the body machine and ABB chamber at the initial time are equal to zero. Thus, the first 14 coordinates of the initial space vector are equal to zero (as follows): While the machine body with an ABB chamber stays at rest, all balls should be distributed in the space inside the ABB chamber. Before the start of the actual simulation experiment, the balls need a bit of time to be distributed on the bottom of the ABB chamber. When the balls finally take a steady position, the motor starts. The initial positions of the balls are set according to the algorithm presented in Fig. 17 and depicted Fig.18.
This algorithm involves a multilayer initial distribution of the balls. The maximal number of rows (n x ) and columns (n y ) depends on the radius of ball R k , radius of ABB chamber R 0 , and depth of chamber B according to the following equations: ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ n x 2 · f loor B/2 − R k 2R k + 1 , The following additional functions were used to calculate the initial positions of the balls: • floor-rounds argument to integer number nearest to zero • mod-returns the remainder after division of the first number by the second one • sign-signum function: returns −1, 0, or 1 This distribution algorithm was applied to all the simulation experiments described in this paper.