Simulation of torque caused by the
lubrication fluid in a ball bearing
David Westsson
March 11, 2015
Abstract
The purpose of the study has been to get a first estimation of how well adapted Computational Fluid Dynamics (CFD) based on the Fi-nite Element Method (FEM) is for studying the large scale flow of lubrication fluid within a ball bearing and the torque arised thereby. Focus has been on calculating the velocities and forces of the lubri-cation fluid in all three directions in order to estimate the unwanted torque.
Sammanfattning
M˚a let med studien har varit ha f˚a en f¨o rsta uppfattning om hur v¨a l anpassad finita elementmetoden ¨a r f¨o r att studera det storskaliga fl¨o det hos sm¨o rjv¨a tskan i ett kullager samt det moment som detta ger upphov till. Tyngdpunkten i arbetet har varit att ber¨a kna hastigheterna och krafterna orsakade av sm¨o rjv¨a tskan i alla tre dimensioner och d¨a rigenom uppskatta det o¨o nskade momentet. Under hela studiens g˚a ng har en fullst¨a ndigt parallell kod
-utvecklad p˚a Ume˚a universitet - anv¨a nts och resultaten visar att den ¨
Contents
1 Introduction 4
1.1 Background . . . 4
1.2 Finite Element Method . . . 5
1.3 Ball Bearings . . . 6 2 Simulation set-up 6 2.1 Included . . . 6 2.2 Excluded . . . 7 2.3 Solver Method . . . 7 2.4 Geometry . . . 7 2.5 Mesh . . . 11 2.6 Boundary Conditions . . . 11 2.7 Model Parameters . . . 12 3 Results 13 3.1 Runtime and Computational Performance . . . 13
3.2 Forces and Torque Inside the Ball Bearing . . . 13
3.2.1 Torque . . . 13
3.2.2 Forces in the Different Directions . . . 15
3.2.3 The velocity field . . . 18 4 Conclusions and Discussion 19
5 Future and On Going Work 20
1
Introduction
1.1 Background
This is the final report of a Bachelor’s Thesis done at the Department of Mathematics and Mathematical Statistics Ume˚a University. The purpose of the study has been to get a first estimation of how well adapted
Computa-5
tional Fluid Dynamics (CFD) based on the Finite Element Method (FEM) is for studying the large scale flow of lubrication fluid within a ball bearing and the torque arised thereby. Focus has been upon finding a good test model and then calculating the velocities and forces of the lubrication fluid in all three directions in order to estimate the unwanted torque.
10
My part of the project has mainly been to design a ball bearing model, and based thereupon supply suitable meshes. The design work has been done in SolidWorks (see Figure 1 and 2 below) and the meshes in ANSYS CFX. All the calculations have been performed in the Akka-cluster of the HPC2N with help of a fully parallel code, developed at the Department of Mathematics
15
and Mathematical Statistics.
Figure 2: Showing the numbers of the balls in relation to the coordinate system.
1.2 Finite Element Method
The main idea of all numerical methods is to consider a continous problem, to descretize it into a finite problem and solve it with the help of a computer. The Finite Element Method (FEM) is a numerical method for solving partial differential equations (PDE’s). Particularly for the viscous incompressible Navier-Stokes equation the problem to be solved can be formulated as: Find the velocity (u) and the pressure (p) such that:
˙ u − ν∆u + (u · ∇)u + ∇p = f ∈ Ω × I (1.1) ∇ · u = 0 ∈ Ω × I (1.2) u = g ∈ ΓD× I (1.3) νn · ∇u − p · n = 0 ∈ ΓN × I (1.4) u(·, 0) = 0 ∈ Ω (1.5)
on the domain Ω, which boundary δΩ is divided into a Dirichlet part (ΓD)
with a ”non-slip” condition and a Neumann part (ΓN) with a ”do-nothing”
condition. g is a given function, n is the unitary vector outward normal to
the boundary, ν is the kinematic viscosity and f a known force. I = (0, T∗) is the applied time interval.
1.3 Ball Bearings
The main purpose of a ball bearing is to decrease the rotational friction and to support radial and axial loads. Though depending upon a particular
5
application the design of the ball bearing can be adapted in enumerous different ways, in order to optimize the behaviour and life length of the bearing. Another important aspect for deminishing the friction forces is the lubrication fluid used. Depending upon the rotational speed of the bearing in combination with the viscosity of said fluid, the frictional forces will vary
10
and in order to minimize them it is important to choose a lubrication fluid with e.g. a vicosity suitable for the rotational speed of interest.
Yet another important feature is the number of balls within the bearing and their arrangement. Like in the self-aligned classical ball bearing (SABB) by Sven Wingquist, the balls can also be lined up in more than one plane and
15
thus allow a certain misalignment between the outer and the inner races. Further, the balls are normally held in place by a so called cage, which presences of course also effects the ease of lubrication and thus the suitable viscosity of the lubrication fluid.
2
Simulation set-up
20
The project is a first trial in preparation for a more detailled study of the frictional forces in ball bearings, thus the geometry and conditions defined below have been kept simple enough for not requiring a too complicated set-up and adaptation of the code, but yet thourough enough for testing the methodology and evaluating the quality of the results.
25
2.1 Included
The trial was set to include: 1. a simple six balls geometry 2. ridgid balls
3. perfect rotation of the balls
30
4. a low Reynold’s number
2.2 Excluded
No particular regards have been taken to:
1. true geometrical models of existing ball bearings 2. non-perfect rotation of the balls
3. Reynold’s numbers with highly turbulent flow
5
4. radial force from work load 5. linear deformation of the balls
6. a two phase flow - with two lubrication oils or lubrication oil mixed with air
2.3 Solver Method
10
There are different possible approches for solving the Navier-Stokes equa-tion, however all of them have advantages and disadvantages emerging from the complexity of a time dependent, non-linear equation.
The method applied in the solver is based on Chorin’s Projection Method which is based on the introduction of the addition and subtraction of an intermediate velocity, ˆu, in each time step, n of Equation 1.1:
un+1+ ˆu − ˆu − un
dt − ν∆u
n+ (un· ∇)un+ ∇˜p = fn (2.1)
decouples the problem into one Crank-Nicolson part for the velocities and one pressure correction part. This splitting introduces extra calculation in form of intermediate velocities, but since the code is completely parallel it is a problem easily overcome by using some extra computational power.
2.4 Geometry
15
The geometrical model was first made in SolidWorks and then exported as an IGS-file to the ANSYS CFX mesh generator, but any common computer aided design program could have been used. The model was constructed in a right-handed coordinate system with the z-axis along the rotational axis of the inner axis.
20
Considering the geometry of the ball bearing (see Figure 3 below), two lay-outs were tested for avoiding mesh problems at the area between the ball and the inner and outer ring respectively. As seen in Figure 4 the balls all have a small collar at the contact zone to make the transition smoother and avoid poorly shaped mesh elements. Figure 5 shows the second layout,
25
which has a small gap between the balls and the respective rings.
flow close to the collars and considering the purpose of the project, it was decided that free hanging balls should be used for further simulations. Another simplification made in order to allow for simple boundary condi-tions on the side walls were the elongated sides of the bearing. In Figure 6 the longitudinal dimensions are shown.
5
Figure 6: A wireframe side view of the ball bearing showing the dimensions in millimeters.
2.5 Mesh
After importation of the IGS-file into the ANSYS CFX module, a fluid body, consisting of the hollow space between the inner wall of the outer shell, the balls and the inner axis, was created and meshed. Throughout the whole project unstructured tetrahedral elements were used, due to there ability to
5
adapt to different types of geometries. All hanging nodes were deleted and the mesh was saved in the TetGen-format in order to be compatible with the solver.
Although most of the mesh sizes created were around 100k-800k nodes, test runs were made with as many as 2.2M nodes.
10
2.6 Boundary Conditions
Since the mesh generator and solver are not automatically connected, a static mesh had to be used. Thus to impose a rotation of the inner axis on the surface nodes of the balls, the point of gravity of the ball bearing was set to coicide with origo and the surface nodes of the balls were given
15
corresponding velocity.
with the z-axis along the rotational axis (ω) of the inner axis of the ball bearing.
2.7 Model Parameters
During the simulations the only physical parameters varied were the (kine-matic) viscosity (ν) and the number of rotations per minute (RP M ). Since
5
the rotation of both the outer and inner ring had to be considered, RP M will equal the sum of their number of rotations per minute.
The values used for the viscosity and the rotation per minute were: 1. ν = 4.0 · 10−3 until 4.0 · 10−3m2s−1 and
2. RP M = 500 to 5000 RP M .
10
Other parameters of interest - that either have been kept constant or have been calculated from set parameters - are:
1. The dimensionless Reynolds number (Re) defined as: Re = L · U
ν , (2.2)
where L is the characteristic length scale of the bearing - i.e. the distance from the center of the bearing to the centre of the balls; U is
15
the characteristic velocity of the bearing - i.e. the velocity of the inner ring at the contact surface with the fluid; ν is the kinematic vicosity of the fluid.
2. Number of nodes (nno) of the mesh.
3. Total simulated time (T∗) - set to T∗= 10 s
20
4. Mesh size (h) defined as the shortest side of the smallest element in the mesh.
5. Constrained Delaunay tetrahedralisation (CDT ) is a constant that is selected to assure the quality of the mesh through fulfilling the Delaunay criteria. Throughout the simulations the vale of CDT = 0.05
25
has been used.
6. Time step dt is calculated for each iteration through: dt = CDT · h
2
8. Finally the ramp parameter (ramp) - for gradually accelerating the system during the first 5th of the simulation time - was defined as:
ramp = 5 · iter niter
, (2.5)
where iter is the number of the current iteration.
The ramping allows a smoother, more realistic acceleration and since it is only applied for a relatively short period in the beginning of the simulation, any eventual effects due to the change from this linear acceleration to a
5
constant speed of the system will be small or even neglectible compared to a system where the full speed is applied directly. Moreover, during this phase it is also possible to save computational power without losing accuracy, by only updating the velocity vectors every 5th iteration.
3
Results
10
The results for the two separate purposes of this study - evaluating the performance of the method and code in the aspects of runtime and accuracy - are presented in the this section below.
3.1 Runtime and Computational Performance
In order to get an understanding for the speed of which the simulation is
15
performed with different numbers of CPUs for the different mesh sizes, the runtimes for the runs where measured (see table 1).
3.2 Forces and Torque Inside the Ball Bearing
For all the different mesh sizes, the torque on the inner axis and the forces on the balls were calculated using a drag function in the code. Below, a
20
part of these results is plotted in order to more comprehensively see the effect the mesh sizes have on the result and later on use that to balance the total runtime and the accuracy. Where the different runs have had similar characteristiques, the results from the 797k run have for the sake of consistency and comparision been taken as example.
25
3.2.1 Torque
Among the different z-torques of the different mesh sizes there is a minor difference and also a certain periodic fluctuation over time. All graphs show instability during the ramp-up phase.
Figure 8: The torque [Nm] in the z-direction on the inner axis due to the lubrication fluid for the different mesh sizes.
The means and the standard deviation (taken after the ramp-up phase) of the different torques shown in Figure 8 are:
5
1. mean32k= 1.07e-6 Nm and std32k = 3.36e-8 Nm
2. mean81k= 1.16e-6 Nm and std81k = 2.70e-8 Nm
3. mean462k= 1.18e-6 Nm and std462k = 3.46e-8 Nm
4. mean797k= 1.23e-6 Nm and std797k = 3.93e-8 Nm
3.2.2 Forces in the Different Directions
10
forces acting on the balls in the respective directions in order to evaluate to consistency of the results. In Figure 9, Figure 10 and Figure 11 this is given for the 797k run.
Figure 10: The forces [N] in the y-direction on the balls in the bearing due to the lubrication fluid for the 797k mesh.
3.2.3 The velocity field
In the search for unexpected behavior it is also informative to depict the velocity field and to some extent visualize the velocity field and the major modes of the solution of the posed problem.
Figure 13: The velocity field for a 232k run with ν = 8e−3 of the lubrication fluid between two of the balls in the bearing.
4
Conclusions and Discussion
Studying the result presented in Figure 8, there are a couple of conclusions that can be directly drawn; the torque of the 32k run differs significantly from the other mesh sizes and the higher number of nodes, the higher value for the top peak of the oscillation. It also seems as the 797k results in two
5
resonance peaks with one of them slowly increasing while the amplitud of the other three rather seem to flaten out. The 462k mesh appears to be the most periodic one. The torque in x- and y-direction respectively seen in Figure 7 is most likely an effect of a not perfectly symmetric mesh, but the effect of it - if needed - is expected to decrease by even a finer meshes.
10
Shifting to figures 9 and 10, one notices a particularity while comparing the forces in the x- and y-direction respectively. For the x-direction the forces seem to pair up with (see Figure ?? in the appendix) balls 1 and 6, balls 2 and 5, and balls 3 and 4 having comparable values. However in the y-direction this symmetry seems broken, since there the pairs with
compa-15
As shown in Figure 11 the force in the z-direction also appears to have one anomali consisting of the force on ball 1 which is not zero as bluntly ex-pected. The reason therefore is assumed also to be due to a non-symmetric mesh, possibly due to an asymmetry introduced while creating the bearing in SolidWorks, even if whenever possible mirroring functions in SolidWorks
5
were used.
Figure 13 shows an appearently smooth velocity field between two of the balls and intuitively opposite balls should have opposite forces acting on them, but since the solution has different eigenmodes (with different wavelengths and directions), what appears to be asymmetries in one or two dimensions
10
in fact is symmetric on the whole in three dimensions. Part of the result of these eigenmodes are pictured in Figure 12, wherefrom it also (especially in slice number 2 and 3 from the top) is clear that the z-component of the velocity field is nonzero with a certain periodicity.
Probably the amplitud of this z-component would be smaller when using
15
a true geometrical model of a ball bearing with a cage, since there is less free space to build up on. The cage might also even out the force between the balls making the x- and y-forces more symmetrical.
Combining the run times with the qualitative results of the toruqes and the forces, one realizes that the 81k run is almost 5 times faster than the 462k
20
run with a fourth of the number of CPUs. Yet the result is comparable to the much finer mesh, what definitely makes it suitable as a test mesh during a forthcoming study, where short runtimes are preferred.
5
Future and On Going Work
25
There is a lot of loose ends in this project, but a natural continuation would be to run the code with SKF’s reference geometry (see Figure 14 below) and settings. After confirmed results yet another step would be to adapt the code and boundary conditions in order to be able to run the simulation with a true geometrical model and conditions. With these preparations everything
30
Figure 14: Drawing of SKF’s benchmark problem with a possibility to add a two phase flow through not filling the whole space with lubrication fluid.
6
References
1. Concepts And Applications Of Finite Element Analysis - Robert Davis Cook and Michael E. Plesha (2001)
2. Finite Element Mesh Generation - B.H.V. Topping and J. Muylle and R. Putanowicz (2002)
5
3. Applied Finite Element Analysis - Larry J. Segerlind (1985)
4. An Introduction to Computational Fluid Dynamics - H. K. Versteeg (2007)
5. Applied Computational Fluid Dynamics Techniques: An Introduction Based on Finite Element Methods - Rainald Lohner (2008)
10
6. An Introduction to Parallel Computational Fluid Dynamics - Francesco Papetti and Sauro Succi
7. Simulation of Lubricant Flow in Ball Bearings - Mats G. Larson and Per Vesterlund (2010)
8. Multiphase Flow Dynamics 1: Fundamentals - Nikolay I. Kolev (2004)
9. Beta Mathematics Handbook: Concepts, Theorems, Methods, Algo-rithms, Formulas, Graphs, Tables - Lennart Rade and Bertil Wester-gren (1997)
10. Schaum’s Mathematical Handbook of Formulas and Tables - Murray R. Spiegel (1998)
