Numerical simulation of ball bearing flow field using the moving particle semi-implicit method

ABSTRACT The distribution of lubricating oil in the bearing cavity is of great significance to bearing lubrication and cooling. A new idea is provided to further study the ball bearings lubrication, to achieve effective lubrication of bearings. The flow field of oil injection lubrication ball bearings is studied by the moving particle semi-implicit (MPS) method. The accuracy of the numerical calculation method is verified by experiments. The oil distribution of the bearing lubrication flow field and churning torque of the ball and cage at different speeds are analyzed. The research results show the maximum oil content in the bearing cavity is distributed in the range of about 90° to 120°. As the increase of rotating speed, the number of particles in the bearing cavity decreases. The churning torque decreases with the increase of rotating speed and increases with the increase of the oil viscosity. This study provides a new numerical calculation method for the lubrication flow field of ball bearings.


C s
Smagorinsky constant,d The spatial dimension of the solution,d N Nozzle diameter, mm f Body forces (per unit volume), N/m 3 h Effective radius, m l Width, mm m Particle mass, kg n i The rotational speed of the inner ring, r/min n i The particle number density,n 0 Initial particle density,n m Rotation speed of the balls and cage, r/min p Pressure, Pa r The distance between two particles, m r Displacement vector of particle, m r b Ball radius, mm r e Finite space for limiting the interactions, m r in Inner radius, mm r m Pitch radius, mm r out Outer radius, mm S Mean strain rate,t Time, s T Oil temperature,°C u * Predicted velocity, m/s

Introduction
In aircraft engines, motor spindles and drive trains, rolling bearings are widely used (Bossmanns & Tu, 1999;Fondelli et al., 2015;Gorse et al., 2006). Its operating performance is very crucial to the power performance and reliability of the whole transmission system (Gloeckner, 2013). In the automotive industry, power transmission is designed to achieve a higher speed (Nagaoka et al., 2015). The bearing will produce a lot of heat at high-speed. Therefore, effective and reasonable lubrication is very essential for high-speed bearings (Pang et al., 2008). Under the condition of high speed, the bearing temperature increases, which leads to bearing failure (Chang et al., 2017). It is of great significance to evaluate and improve the lubrication status of bearings by studying the distribution characteristics of the internal flow field of bearings under the condition of injection lubrication. The flow field in the bearing cavity is relatively complex, and there is a strong interaction between the fluid and cage and roller (Gao et al., 2018). Some scholars have obtained the internal lubrication characteristics of bearings through experiments. An X-ray computed tomography (CT) imaging system was used to capture the grease distribution in a ball bearing and the internal grease flows (Noda et al., 2016). A specially designed test rig was established, and the flow field data were obtained by particle image velocimetry (PIV) technology (Lee et al., 2005). To study the impact of the oil supply flow rate and rotation speed on the oil distribution, Wu et al. built the bearing temperature and flow visualization test rig (Wu, Hu, Hu, et al., 2016). Aidarinis et al. (2011) established an experimental platform to test the change of the air-flow field in the bearing cavity. The lubricant fluid velocity field in the bearing was measured by the laser Doppler anemometry (LDA) system. It is not easy to test the flow field in the bearing cavity through experiments because the bearing structure is compact and the internal space is small. Computational fluid dynamics (CFD) method for the cooling flow inside the rolling element bearing was proposed. Bristot et al. (2016) used the CFD method to research the gas-liquid two-phase flow and the movement of the liquid droplet. Crouchez-Pillot and Morvan (2014) carried out numerical simulation calculations on the fluid field of bearing chamber in aero-engine, and adopted the volume of fluid (VOF) method and grid adaptive technology to track the flow field. Adeniyi et al. (2017) used the coupled levelset VOF (CLSVOF) method to calculate the lubrication flow field in the bearing cavity. It was found that the variation of the wetted area with shaft speed and the oil flow rate was inconsistent. Gao et al. employed the finite volume method(FVM) to compute the churning loss of roller bearings, and studied the distribution of oil and the change law of volume fraction with rotating speed (Gao et al., 2018;. Then, Gao et al. also researched the distribution of lubricating oil under the bearing ring by CLSVOF multi-phase flow method . Wu and Hu et al. (Hu et al., 2014;Wu et al., 2017) studied the distribution of bearing oil volume fraction and temperature by VOF model and analyzed the influence of oil supply flow, speed and nozzle size on oil distribution and temperature. Yan et al. (2018) defined the rotating coordinates and established the lubrication flow field model of the bearing chamber. The visual flow field test scheme based on PIV technology for bearing chamber flow monitoring was proposed. The fluid distribution and vortex field in the bearing chamber were obtained.
All the above-mentioned researches mainly focused on bearing lubrication flow field with the Euler mesh method. When simulating the characteristics of the bearing lubrication flow field, the results are closely related to the type and quality of mesh generation. For complex geometry and steady and time-varying interacting loads, some results may be affected by the combination of grid modules and the flow information transfer mode. It not only dramatically adds a lot of calculation, but also increases the risk of error expansion and calculation instability.
In 1995, Koshizuka et al. (1995) proposed the MPS method for the first time. The MPS method is to divide the continuous fluid domain into a series of Lagrangian particles. The interaction between particles is realized through the integration of kernel function (1996). The advantage of the MPS method in numerical simulation is that compared with the traditional meshed CFD simulation method, it can save a lot of tedious work of mesh generation when dealing with the flow on the complex interface and large deformation flow on the free surface. Because of the flexibility, robustness and high efficiency of fluid analysis, the simulation of the complex flow field is close to the actual situation. For the past few years, the MPS method has been largely used in numerical calculation of marine engineering and mechanical engineering Duan et al., 2017Duan et al., , 2018Gotoh & Khayyer, 2018;Guo et al., 2020;Sun et al., 2017).
The description of the lubricant flow field distribution of bearing is of great significance to the prediction of bearing temperature. The application of an efficient and high-precision method of calculating bearing lubrication is of inestimable value to bearing optimization design. In this paper, the MPS numerical calculation method is applied to the distribution calculation of the bearing lubrication flow field. The flow field patterns associated with the oil injection lubrication at different speeds are Figure 1. Drawing of oil injection lubrication mechanism (a) Three-dimensional drawing (b) Two-dimensional drawing (Wu, Hu, Yuan, et al., 2016). studied. The aim is to propose a new method for calculating the flow field distribution of bearing lubrication.

System description
The drawing of the two-dimensional and three-dimensio nal structure of ball bearing jet lubrication is shown in Figure 1. Lubricating oil is sprayed from the nozzle to the bearing balls, cage, and inner and outer rings. The oil interacts strongly with the rotating balls, cage and inner ring.
The type of angular contact ball bearing used in this experiment is SKF 7210 BECBP. The structural dimension parameters of the SKF 7210 BECBP are given in Table 1. The lubricating oil is an automatic transmission fluid, which is ATF DEXRON III. The parameters of lubricating oil are shown in Table 2. The spindle rotation causes the bearing inner ring to rotate, and the interaction force between the inner ring and the roller causes the roller and cage to move. The motion of the rolling elements will cause a complex flow of the lubricant oil. In the process of calculation, the rotational speed of bearing balls and cage needs to be determined. The formula for calculating the rotational speed of balls and cage was proposed by Harris and Kotzalas (2006).
where n i is the rotation speed of the inner ring, r b is the rolling element radius, r m is the pitch radius, α is the contact angle.

Governing equations
The numerical calculation method is a meshless method. The fluid will be modeled with particles. The discrete forms of governing equations of CFD include finite difference method(FDM), FVM and finite element method(FEM) (Ghalandari et al., 2019;Salih et al., 2019). Compared with the FVM, the convection term can be neglected. It is easy for solving the equations. Continuity equation: Momentum equation: where t is the time, ∇ is the gradient, f represents the body forces (per unit volume), u is the velocity vector and p is the pressure. Koshizuka and Oka (1996) proposed the kernel function. It defines the weight relation of the interaction between neighboring particles. This function is as follows:

Kernel function
where r e is the finite space for limiting the interactions, r = |r j − r i | is the distance between two particles. When r > r e , there is no interaction between particles, particles only interact with particles within the radius of r e .

Particle number density
Particle density reflects the distribution of particles in the flow field. If the number density of particles is constant, the density of the liquid is constant. That is to ensure the incompressibility of the fluid. The number density of particle is the weighted sum of the number of particles in the radius of action, that is: where r is the displacement vectors of particles. The density of the fluid is as follows: where m is the particle mass.

Gradient operator model
In addition, based on the kernel function, the gradient vector of particle i can be obtained by weighted averaging with the gradient vector of other particles.
where d is the spatial dimension of the solution, φ is the physical parameter of particle, n 0 is the initial particle density.

Laplacian operator model
The diffusion term in the momentum conservation equation is expressed by the Laplacian operator ∇ 2 φ. The formula is as follows: where λ is the correction factor. The calculation formula is:

Numerical scheme for the flow analysis
The momentum equation is solved by a semi-implicit method. In the first step, the effects of viscous force and body force are considered. Then the pressure equation is solved by using the velocity field. Finally, the velocity is corrected by calculating the pressure gradient.
The most important aspect in the calculation method is that the pressure field is estimated based on the Poisson equation. The criterion of the end of calculation is to check whether the end time of calculation has been completed and, if not, go to the next time step to calculate. The process of numerical calculation is shown in Figure 2. The inlet boundary condition is the flow rate, and the outlet boundary condition is atmospheric pressure. As the oil entering the bearing cavity is stirred by the high-speed rotating bearing inner ring, roller and cage, the turbulence model is selected. To improve simulations where turbulence plays an important role, the widely used model is implemented into the method (Smagorinsky, 1963). It is a zero equation turbulence model. The turbulent viscosity coefficient formula is as follows: X is the particle space, C s is the Smagorinsky constant, andS is the mean strain rate. The heat transfer model isn't considered in the calculation process. Considering the viscous term in the momentum equation, the time step needs to meet the Morris condition (Liu & Liu, 2003). The time step is 0.0001 s according to the following formula.
h is effective radius. Generally, the particle radius has a certain influence on the calculation results. In order to ensure that there can be three layers of particles between the cage and the inner or outer ring, the particle radius of 0.35, 0.4 and 0.5 mm are used to verify the independence of particle radius. When the rotating speed of bearing is 600 r/min, the results of churning torque of balls and cage under different particle radius are shown in Table 3. The difference of churning torque between particle radius of 0.4 and 0.35 mm is less than 10%, considering the economy of calculation. In all the following calculations, the particle radius is 0.4 mm. The momentum equation can be discretized into: where Equation (12) can be written as the following equation where predicted velocity u * is And Equation (14) can be formulated as: Equation (15) can be formulated as: Apparently, the left side of the equation is 0. The equation is transformed into The original Poisson Pressure Equation proposed by Koshizuka and Oka (1996) is The following modified Poisson Pressure Equation was proposed by Khayyer and Gotoh (2009).
According to the Poisson Pressure Equation of Equation (8), we can get the following equations:

Free surface boundary condition
Boundary conditions are essential for solving the Poisson pressure equation. Commonly, the ambient pressure is the boundary condition of the free surface. Shibata et al. (2015) proposed the concept of virtual particles. The calculation formula of virtual particles is as follows: Then the Poisson equation is

Inlet and outlet boundary
For the inlet boundary, the new particles with a given velocity are generated over time. The inlet particles are involved in the numerical density calculation, but not in the pressure calculation. Before the inlet particles cross the inlet border, they play as 'inlet particles' with an inlet velocity. When the inlet particles exceed the inlet border, they are transformed into actual fluid particles automatically.
For the outlet boundary, the outlet pressure is given to the virtual particles outside the border. The fluid particles are deleted, when they exceed the outlet border. The particle diagram of the inlet and outlet boundary is shown in Figure 3.

Experimental apparatus
To obtain the oil distribution of bearing lubrication at different speeds, single lens reflex (SLR) camera and other   Figure 4. The motor controller can control the speed of the motor, thus controlling the speed inner ring. The peristaltic pump is composed of gear pump and brushless DC motor, which can regulate the flow rate of the oil supply. The specific parameters of the test apparatus and sensor are shown in Table 4. When taking pictures with SLR camera, in order to obtain a clear picture of the distribution of lubricating oil in the bearing cavity, it is necessary to select an appropriate shutter speed according to the rotating speed of the bearing. The selection basis of shutter speed is shown in Table 5. During the test, first irradiate the parallel light on the bearing end face, then select the appropriate shutter speed, set the appropriate aperture, and finally take photos.

Distribution of lubricating oil
Numerical simulation and visualization experiments are used to study the internal flow field of oil-jet lubrication bearings at speeds of 600, 1200, 2400, 4800, 6000, 9000, 10000 and 12000 r/min. The nozzle radius is 0.5 mm. The flow rate of the oil supply is 500 mL/min. The injection speed of lubricating oil is about 10.6 m/s, which can meet the demand of injection speed. As shown in Figure 5, visualization test results of the flow field at different speeds are presented. Most of the oil is ejected from the gap between the rolling element and the inner and outer raceways, and between the rolling element and the cage when the bearing is stationary. The cage, rolling element and the end face of inner and outer rings can be clearly seen, because the whole bearing cavity is not filled with oil. When the rotation speed is less than 1200 r/min, the bearing cavity is covered with a thick lubricating oil layer. When the rotational speed is greater than 2400 r/min, the oil layer attached to the bearing cavity becomes thinner. With the increase of speed, the contour of the ball and cage becomes more and more evident. Therefore, it can be known that the oil content decreases with the increase of rotating speed. It is because that the oil tends to be thrown centrifugally with a higher rotation speed from the bearing cavity. The storage time of the oil becomes shorter.
When the speed is 4800 r/min, the profile of the cage and balls in the bearing cavity is becoming more and more apparent along the circumferential direction, and there is less oil attached to the surface. There is relatively a smaller amount oil of upstream position of the oil outlet, which is consistent with the calculated results (Hu et al., 2014). In addition, the oil passing through the bearing cavity mainly leaves the bearing from the downstream direction of the nozzle. From the nozzle position along the bearing rotation direction of 120°, the end face oil layer of the bearing outer ring is relatively thick. Following along the direction of rotation, the oil layer on the outer ring face begins to decrease and shrink until the bottom of the outer end face of the bearing.
As shown in Figure 6, the measured and calculated visual flow field results are displayed at speeds of 600 and 1200 r/min. When the speed is 600 r/min, a large amount of oil is observed in the bearing cavity. The clearance between cage and balls, between balls and inner ring, and between balls and outer ring are covered by a thick oil layer, making it difficult to distinguish the boundary between bearing cage and balls. The oil content in the A region of the end face is relatively high. The oil is first thrown into this region by centrifugal force. Mainly due to gravity, part of the oil in the A and C regions is collected in the B region. There is less oil in the C region, because it is far from the jet position. It can be seen in the D region that the oil directly flies away from the bearing, which is closely related to the jet position. When the speed is 1200 r/min, it has a similar effect. The numerical calculation results agree with the experimental results. The reason for the difference between the numerical calculation results and the measured results may be that the effect of air on the fluid is not considered in the numerical calculation.
As shown in Figure 7, the measured and calculated results of the visual flow field are displayed at speeds of 4800 r/min and 12000 r/min. As the rotation speed increases, the centrifugal force acting on the oil entering the bearing cavity becomes stronger and stronger. Therefore, the oil is thrown out, leaving less oil in the bearing cavity attached to the cage and ball surface. Ball surface and the contour of cage and ball become more and more apparent. The difference between the test results and the simulation results may be due to the effect of air on the fluid, cage and ball when the bearing is at high speed during the test. The impact of air on the bearing lubrication flow field is not considered in the simulation computation. In addition, during the experiment, the oil temperature will change to a certain extent. It is assumed that the oil temperature is constant in the numerical calculation, which may also be the reason for the difference. Figure 8 shows the numerical calculation results of the flow field at the outlet side. It is similar to the flow field changes with the other side at different speeds. As the rotation speed increases, the oil layer in the bearing cavity becomes thinner, and the contour of the cage and ball gradually becomes clear. In addition, more broken oil droplets can be seen in the bearing cavity. Therefore, at high-speed, it is difficult to form a complete oil film inside the bearing cavity.
As shown in Figure 9, the enlarged visualization flow field of the test results at different speeds is presented. It shows the evolution process of the bearing cavity internal flow field on both sides of the jet position at different rotating speeds. It can be observed from the enlarged picture that the oil layer in the bearing cavity downstream of the jet position is slightly higher than that in the upstream. In addition, some lubrication oil penetrates the bearing cavity. The remaining oil flows into the bearing cavity and will move with the balls and cage. Due to the small inertia of the oil, it will accelerate rapidly to the same speed as the cage and the inner ring. It is thrown out of the bearing cavity and leaves from the bearing cavity downstream of the jet position. As the rotation speed increases, the region away from the bearing shrinks towards the jet position. When the rotation speed is 12000 r/min, under the irradiation of parallel light source, the surface of the bearing ball has an evidently reflective surface, which indicates that the oil layer thickness is relatively thin.

Particle number distribution and churning torque
To quantitatively analyze the oil content, a cuboid zone is established in the bearing cavity to count the number of particles in each region. When the oil flow rate is 500 mL/min, the rotation speed is 600 and 1200 r/min, and the nozzle diameter is 1 mm, the distribution of the number of particles in the bearing cavity along the rotation direction is shown in Figure 10. The number of particles decreases from the nozzle position to about 30°along the direction of rotation, then increases from about 30°t o 110°, and finally decreases from about 110°to the nozzle position. The number of particles is relatively large between about 90°and 120°, which means that there is more oil in the bearing cavity. It is consistent with the   Figure 6. As the rotation speed increases, the number of particles in the corresponding position along the rotation direction becomes smaller. When the oil flow rate is 500 and 700 mL/min, the injection speed is constant, and the rotation speed is 1200 r/min, the distribution of particle number along the rotation direction is shown in Figure 11. As the increase of flow rate, the number of particle in the corresponding position along the rotation direction increases. The behavior of particle number along the rotation direction is the same.
When the flow rate of oil supply is 500 mL/min, the nozzle radium is 0.5 mm, and the kinetic viscosity is 0.00015 m 2 /s and 0.0001 m 2 /s, the variation of oil churning torque of bearing cage and ball with rotating speed is presented in Figure 12. The churning torque of ball and cage increases linearly with the increase of rotating speed. Peterson et al. (2021) obtained that the drag force increases with the increase of rotation speed. There is the same trend with the results of the study. In addition, at the same speed, the churning torque increases with the increase of kinetic viscosity.

Conclusions and future work
The MPS numerical method is applied to the visualization of ball bearings lubrication flow field. The numerical results are verified by visualization experiments. The oil distribution in bearing cavity under different speeds is discussed. The variation of churning torque between roller and cage at different speeds, flow rates and viscosity are studied. The conclusion is as follows: (1) The oil distribution along the rotation direction obtained by MPS numerical calculation method has the same trend as the experimental results. It is uneven with the distribution of oil in the bearing cavity (2) Along the rotation direction, the number of particles first decreases, then increases and finally decreases. The number of particles is relatively large between about 90°and 120°. (3) As the increase of speed, the number of particles in the bearing cavity decreases. As the increase of flow rate, the number of particles in the bearing cavity increases. (4) The churning torque of the bearing cage and balls changes linearly with the increase of rotation speed. At the same speed, as the decrease of kinematic viscosity, the churning torque of the ball and cage decreases.
Finally, it is worth mentioning that the air has a significant impact on the oil distribution at high speed. Next, the oil-air two-phase flow will be added to the numerical model and we will compare and analyze with Peterson et al. (2021). That will be of significant benefit to the design of ball bearings.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Funding
This work was supported by National Natural Science Foundation of China [grant number 51975045].