# A New Integrated Approach for the Prediction of the Load Independent Power Losses of Gears: Development of a Mesh-Handling Algorithm to Reduce the CFD Simulation Time

Authors: F. Concli; A. Della Torre; C. Gorla; G. Montenegro
Journal: Advances in Tribology
 (2016)

---

## Abstract

To improve the efficiency of geared transmissions, prediction models are required. Literature provides only simplified models that often do not take into account the influence of many parameters on the power losses. Recently some works based on CFD simulations have been presented. The drawback of this technique is the time demand needed for the computation. In this work a less time-consuming numerical calculation method based on some specific mesh-handling techniques was extensively applied. With this approach the windage phenomena were simulated and compared with experimental data in terms of power loss. The comparison shows the capability of the numerical approach to capture the phenomena that can be observed experimentally. The powerful capabilities of this approach in terms of both prediction accuracy and computational effort efficiency make it a potential tool for an advanced design of gearboxes as well as a powerful tool for further comprehension of the physics behind the gearbox lubrication.


---

## Body

## 1. Introduction

More and more severe requirements in terms of efficiency have encouraged the gearbox manufacturers to increase the investments for power dissipation reduction. An energy efficient solution, besides a pure power consumption reduction, shows lower operating temperatures and therefore has a higher reliability. The losses are dissipated in form of heat through the housing and shafts. While the exchange area/volume is often limited by external constrains, the internal design of a gearbox is completely in charge of the gearbox manufacturer. Considering that in industrial gearboxes the power losses are more or less equally shared between the bearings, the gear meshing, and the interaction with the oil, it is clear that at least this main contribution (there are other contributions related, e.g., to the seals) should be considered when a new gearbox is developed. Being able to correctly model the power losses enables the capability to predict also the operating temperatures under each operating condition and consequently the so-called thermal limit. While for bearings and for the gear meshing losses accurate models already exist [1–5], for the oil churning/sloshing power losses only basic research was carried out. All the experimentally derived models results are accurate only as long as the operating conditions and geometrical configurations are similar to those used in the original experiments. One of the first works on this topic was carried out by Mauz [6]. Mauz tested several geometries and operating conditions. From the authors experience [7] this model is not capable of taking into account very important parameters like the helix angle. Other studies in this field were carried out by other researchers such as Dawson [8] who concentrated on windage, Seetharaman and Kahraman [9] who concentrated on churning, and several other authors [10–14]. The authors maintain that a deeper understanding of the physical phenomena can only be achieved by means of numerical techniques. Different contributions have already shown the goodness of the CFD approach for such purposes [15, 16]. In the past the authors approached the problem using commercial software [17–22]. The main limitation to the wide diffusion of such approaches is related to the computational effort, the complexity of setting up a model with a general-purpose software, and the costs for the licenses. In this scenario the authors started studying a simplified configuration [23, 24] with the open-source code OpenFOAM [25]. This approach allows both to overcome the license costs (GNU license) but also, thanks to the specificity of the developed algorithm, to reduce the setup and calculation times. In this paper a mesh generation and handling algorithm that enables the simulation of the gear engagement is presented. The developed code was tested and the results show good agreement compared with experimental data [17].


## 2. Power Loss Classification

According to [5] the power losses of gears can be subdivided into load dependent and load independent losses and again according to the mechanical component that is responsible for their generation(1)
P
L
=
P
L
G
+
P
L
G
0
+
P
L
B
+
P
L
B
0
+
P
L
S
+
P
L
X
.The subindexes G, B, S, and X stand for gears, bearings, seals, and other generic components (like clutches or synchronizers). The subindex 0 indicates the load independent losses. The focus of this research is on the load independent gearing power losses P
L
G
0. It results from the interaction of the rotating (or in general moving) mechanical components and the lubricant. This load independent contribution of the gears can further be subdivided into windage, churning, and squeezing/pocketing effects. Windage and churning as in external aerodynamics concern the main interaction of the components with the oil/air or with the oil-air mixture. The pocketing effects instead are related to the volume variation of the gap between the teeth in the mating region and the additional fluxes that take place.


## 3. Objective of the Study

The main objective of this study is to find a reliable method in order to calculate the load independent losses of gears. Beside the accuracy of the results the models should also be manageable in terms of computational times so as to enable systematic studies in reasonable times, a requirement that is particularly stringent for industrial applications.


## 4. Numerical Modeling

### 4.1. Conservation Equations

For a generic volume it is possible to write five conservation equations: the averaged mass conservation equation for no-stationary incompressible flows and the averaged momentum equation (and eventually the energy conservation equation):(2)
∂
ρ
u
i
∂
x
i
=
0
,
∂
∂
t
ρ
u
i
+
∂
∂
x
j
ρ
u
i
u
i
=
-
∂
∂
t
p
+
∂
∂
x
j
μ
∂
u
i
∂
x
j
+
∂
u
j
∂
x
i
-
∂
∂
x
j
τ
i
j
,where x
i
,
j represent the Cartesian coordinate and u
i
,
j is the velocity component, p the pressure, ρ the density, and τ
i
j the Reynolds term.Solving these equations numerically it is possible to describe the fluid behavior in terms of macroscopic properties such as pressure, velocity, and their space and time derivatives. Simulations were run adopting a transient incompressible pressure-velocity coupled solver. In particular, the implementation of the pimpleDyMFoam solver of OpenFOAM (PIMPLE solver, a flexible implementation of a transient solver that allows operation in both PISO [26] and transient SIMPLE mode) was adopted. In this work, the temperature of the lubricant was assumed as uniform in the domain and a priori known; therefore the energy conservation equation was not considered. In a similar way the density and the viscosity values were set on the basis of the temperature.


### 4.2. Boundary Condition

In order to correctly move the gear boundaries according to the kinematic laws, a new boundary condition was developed. The relation on which the boundary condition is based is(3)
ψ
n
=
ψ
n
-
1
+
ω
·
d
tin which ψ
n represents the angular position of the gear at the actual time step, ψ
n
-
1 the angular position at the previous time step, ω the rotational speed of the gear, and d
t the time step. In this manner it was possible to move the gear boundaries according to the kinematic law.


### 4.3. Mesh

OpenFOAM is capable of handling the domain deformations by moving the mesh. The objective of the mesh motion is to accommodate externally prescribed boundary deformations by changing the positions of mesh points (Figure2). During the motion the mesh must remain geometrically valid. The motion u is calculated considering the Laplace smoothing equation (4) and the pseudo-solid equation (5) (linearization of the motion equations for small deformations) (4)
∇
·
γ
∇
u
=
0
,
(5)
x
n
+
1
=
x
n
+
Δ
t
u
,where γ represents the diffusivity, x the position of the nodes of the mesh, and t the actual time.This approach is effective only as far as the deformation does not affect too much the quality of the mesh which should be updated in this case. Different from many commercial codes, OpenFOAM does not allow local mesh regeneration to circumvent the mesh quality degeneration.In order to simulate the complete rotation and meshing of the gears, an approach which employs a multiple number of meshes was applied [27, 28]. In this framework, each mesh is generated for one specific rotational angle and has a certain angular validity. Once a mesh is created, motion is imposed at the boundaries, whereas the inner points of the grid are moved according to the solution of the Laplace equation in order to adapt to the boundary motion. Usually, when the mesh is deformed, its quality in terms of skewness, nonorthogonality and geometrical/topological validity decreases. As a consequence, when the deformation becomes excessive a new mesh must be created and the solution mapped from the old mesh onto the new one resorting to field interpolation techniques. In this case, the interpolation of the computed flow field from one mesh to the next one is performed by means of a second-order, inverse distance weighting method. The conservativeness of the approach was tested in order to verify that the errors induced by interpolation were negligible both in terms of integral qualities (fluid mass) and velocity profiles.In OpenFOAM two main meshing tools are available. The first one, blockMesh, is capable of creating parametric meshes with grading and curved edges; blockMesh is however not suitable to handle complex geometries. SnappyHexMesh, another utility of OpenFOAM, generates 3-dimensional meshes containing hexahedra and split-hexahedra automatically from triangulated surface geometries in stereolithography format. The mesh approximately conforms to the surface by iteratively refining a starting mesh (generated with blockMesh) and morphing the resulting split-hex mesh to the surface. The surface handling is robust with a prespecified final mesh quality but the procedure is also very time consuming and not suitable for the purposes of this study. For this reason, a new and more efficient approach was tested. The geometry is generated with SALOME [29]. A python script allows the parametrization of the geometrical entities that are defined through an analytical formulation. In a second step, the domain is meshed. In order to reduce the computational time a special strategy was implemented. The domain (a box with two mating gears) is cut by a plane that lies on the gear sides. The generated internal faces (Figure 1) are meshed with Netgen [30]. Netgen follows a top down strategy. It starts by computing the corner points. Then the edges are defined and meshed into segments. Next the faces are discretized by an advancing front surface mesh generator. A fast Delaunay algorithm generates most of the elements, but sometimes it fails for the last elements, so a slower back-tracking rule-base algorithm takes over. Basically the Delaunay algorithm subdivides the edges into segments according to the local seed prescribed element size (Figure 3(a)). Then, the surfaces are seeded with internal points (Figure 3(b)). These were connected (triangulation) together. The results of such procedure is often a surface (or volume) mesh that does not match the boundary mesh (Figure 3(c)) so the outer and the intersecting elements are removed (Figure 3(d)).Figure 1
Discretization of the internal planar faces: this mesh will then be extruded in both directions.
(a)
(b)Figure 2
Detail of the first 4 meshes (symmetry plane); before mapping the data, the(
i
-
1
)th mesh is deformed so that its boundary will coincide with the boundaries of the ith mesh.Figure 3
Steps performed by the Delaunay algorithm.
(a)
(b)
(c)
(d)The internal two-dimensional mesh is successively extruded in both directions generating a 2.5D discretization. This procedure reduces significantly the mesh generation times: even if the computational domain is three-dimensional, the meshing algorithms are applied only to two-dimensional surfaces and the resulting meshes copied to create swept elements.The elements belonging to each boundary can be grouped taking advantage of the analytical geometrical entities. The mesh is finally exported and converted into the OpenFOAM format.For this test case the whole geometry and mesh generation takes approximately 30 seconds on a single core 12.8 GFLOPS machine. On the previously tested commercial code [17], meshing the same geometry (with tetrahedrons) takes approximately 6 minutes (same amount of cells—13.2 GFLOPS machine). The algorithm is so effective that also the local mesh regeneration, performed by the commercial code [17] after each iteration to regenerate the collapsed cells, takes more time than the actual procedure takes for the complete remeshing of the domain.


## 4.1. Conservation Equations

For a generic volume it is possible to write five conservation equations: the averaged mass conservation equation for no-stationary incompressible flows and the averaged momentum equation (and eventually the energy conservation equation):(2)
∂
ρ
u
i
∂
x
i
=
0
,
∂
∂
t
ρ
u
i
+
∂
∂
x
j
ρ
u
i
u
i
=
-
∂
∂
t
p
+
∂
∂
x
j
μ
∂
u
i
∂
x
j
+
∂
u
j
∂
x
i
-
∂
∂
x
j
τ
i
j
,where x
i
,
j represent the Cartesian coordinate and u
i
,
j is the velocity component, p the pressure, ρ the density, and τ
i
j the Reynolds term.Solving these equations numerically it is possible to describe the fluid behavior in terms of macroscopic properties such as pressure, velocity, and their space and time derivatives. Simulations were run adopting a transient incompressible pressure-velocity coupled solver. In particular, the implementation of the pimpleDyMFoam solver of OpenFOAM (PIMPLE solver, a flexible implementation of a transient solver that allows operation in both PISO [26] and transient SIMPLE mode) was adopted. In this work, the temperature of the lubricant was assumed as uniform in the domain and a priori known; therefore the energy conservation equation was not considered. In a similar way the density and the viscosity values were set on the basis of the temperature.


## 4.2. Boundary Condition

In order to correctly move the gear boundaries according to the kinematic laws, a new boundary condition was developed. The relation on which the boundary condition is based is(3)
ψ
n
=
ψ
n
-
1
+
ω
·
d
tin which ψ
n represents the angular position of the gear at the actual time step, ψ
n
-
1 the angular position at the previous time step, ω the rotational speed of the gear, and d
t the time step. In this manner it was possible to move the gear boundaries according to the kinematic law.


## 4.3. Mesh

OpenFOAM is capable of handling the domain deformations by moving the mesh. The objective of the mesh motion is to accommodate externally prescribed boundary deformations by changing the positions of mesh points (Figure2). During the motion the mesh must remain geometrically valid. The motion u is calculated considering the Laplace smoothing equation (4) and the pseudo-solid equation (5) (linearization of the motion equations for small deformations) (4)
∇
·
γ
∇
u
=
0
,
(5)
x
n
+
1
=
x
n
+
Δ
t
u
,where γ represents the diffusivity, x the position of the nodes of the mesh, and t the actual time.This approach is effective only as far as the deformation does not affect too much the quality of the mesh which should be updated in this case. Different from many commercial codes, OpenFOAM does not allow local mesh regeneration to circumvent the mesh quality degeneration.In order to simulate the complete rotation and meshing of the gears, an approach which employs a multiple number of meshes was applied [27, 28]. In this framework, each mesh is generated for one specific rotational angle and has a certain angular validity. Once a mesh is created, motion is imposed at the boundaries, whereas the inner points of the grid are moved according to the solution of the Laplace equation in order to adapt to the boundary motion. Usually, when the mesh is deformed, its quality in terms of skewness, nonorthogonality and geometrical/topological validity decreases. As a consequence, when the deformation becomes excessive a new mesh must be created and the solution mapped from the old mesh onto the new one resorting to field interpolation techniques. In this case, the interpolation of the computed flow field from one mesh to the next one is performed by means of a second-order, inverse distance weighting method. The conservativeness of the approach was tested in order to verify that the errors induced by interpolation were negligible both in terms of integral qualities (fluid mass) and velocity profiles.In OpenFOAM two main meshing tools are available. The first one, blockMesh, is capable of creating parametric meshes with grading and curved edges; blockMesh is however not suitable to handle complex geometries. SnappyHexMesh, another utility of OpenFOAM, generates 3-dimensional meshes containing hexahedra and split-hexahedra automatically from triangulated surface geometries in stereolithography format. The mesh approximately conforms to the surface by iteratively refining a starting mesh (generated with blockMesh) and morphing the resulting split-hex mesh to the surface. The surface handling is robust with a prespecified final mesh quality but the procedure is also very time consuming and not suitable for the purposes of this study. For this reason, a new and more efficient approach was tested. The geometry is generated with SALOME [29]. A python script allows the parametrization of the geometrical entities that are defined through an analytical formulation. In a second step, the domain is meshed. In order to reduce the computational time a special strategy was implemented. The domain (a box with two mating gears) is cut by a plane that lies on the gear sides. The generated internal faces (Figure 1) are meshed with Netgen [30]. Netgen follows a top down strategy. It starts by computing the corner points. Then the edges are defined and meshed into segments. Next the faces are discretized by an advancing front surface mesh generator. A fast Delaunay algorithm generates most of the elements, but sometimes it fails for the last elements, so a slower back-tracking rule-base algorithm takes over. Basically the Delaunay algorithm subdivides the edges into segments according to the local seed prescribed element size (Figure 3(a)). Then, the surfaces are seeded with internal points (Figure 3(b)). These were connected (triangulation) together. The results of such procedure is often a surface (or volume) mesh that does not match the boundary mesh (Figure 3(c)) so the outer and the intersecting elements are removed (Figure 3(d)).Figure 1
Discretization of the internal planar faces: this mesh will then be extruded in both directions.
(a)
(b)Figure 2
Detail of the first 4 meshes (symmetry plane); before mapping the data, the(
i
-
1
)th mesh is deformed so that its boundary will coincide with the boundaries of the ith mesh.Figure 3
Steps performed by the Delaunay algorithm.
(a)
(b)
(c)
(d)The internal two-dimensional mesh is successively extruded in both directions generating a 2.5D discretization. This procedure reduces significantly the mesh generation times: even if the computational domain is three-dimensional, the meshing algorithms are applied only to two-dimensional surfaces and the resulting meshes copied to create swept elements.The elements belonging to each boundary can be grouped taking advantage of the analytical geometrical entities. The mesh is finally exported and converted into the OpenFOAM format.For this test case the whole geometry and mesh generation takes approximately 30 seconds on a single core 12.8 GFLOPS machine. On the previously tested commercial code [17], meshing the same geometry (with tetrahedrons) takes approximately 6 minutes (same amount of cells—13.2 GFLOPS machine). The algorithm is so effective that also the local mesh regeneration, performed by the commercial code [17] after each iteration to regenerate the collapsed cells, takes more time than the actual procedure takes for the complete remeshing of the domain.


## 5. Torque Calculation

The total resistant torque around a specified axis is computed by summing the cross products of the pressureF
P and viscous force vectors F
τ for each face and the vector r
A
B, defining the distance of each face center from the specified axis (6)
M
A
=
r
A
B
F
P
+
r
A
B
F
τ
.The pressure and viscosity force vectors are calculated as(7)
F
P
=
p
A
i
,
F
τ
=
μ
∂
u
∂
x
A
i
,where A
i represents the area of the i-face of the mesh.


## 6. Model Validation

The presented model was validated on the basis of some already available experimental results [7, 17]. The tests were performed in 2012 at FZG [31]. Two gears, whose properties are listed in Table 1, are set into rotation with an electric motor whose rotational speed and torque are measured by sensors.Table 1
Gear parameters.
Widthb  [mm]
40
Tip diameterd
a  [mm]
102.5
Helix angleβ  [°]
0
Pressure angleα  [°]
20
Normal modulem
n  [mm]
4
Center distancea  [mm]
91.5
Number of teethz  [
]
23
Chordal thicknesss
j  [mm]
5.35Both gears are mounted on cantilevered shafts each one supported by two rolling bearings. The primary shaft is connected to the motor while the second rotates freely. The power supplied by the motor is therefore completely dissipated by the gears (and shafts), by the bearings, and by the seals. Dedicated tests without the gears were performed in order to quantify the power losses due to bearings, shafts, and seals. In this manner it was possible to separate the pure gear loss contribution.The dimensions of the test rig housing are 290 × 145 × 185 mm. The volume has been completely filled with an FVA2 lubricant and pressurized to 6 bars. The lubricant properties are listed in Table2.Table 2
Fluid properties.
ρ
15  K
g
/
m
3
ν
100  m
m
2
/
s
ν
100  m
m
2
/
s
871
29.8
5.2Measurements were performed up to 4000 rpm (22 m/s of tangential velocity on the pitch diameter). The measured values are the sum of the contributions given by the 4 rolling bearings, the two seals, the shafts interaction with the lubricant, and the gear interaction with the lubricant.This configuration was simulated with the above described approach. Despite the cantilevered shafts, the model was considered symmetric and only half of the domain was modeled.


## 7. Results

Figure4 shows a comparison between the numerical results of present work (OpenFOAM), the experimental ones, and those obtained by the authors with the commercial general-purpose software in a previous research [7, 17]. The results are presented in terms of resistant torque versus rotational speed. The stripes represent a 3rd-order polynomial extrapolation of the numerical results.Figure 4
(a) Gear load independent resistant torque: experimental versus numerical data and (b) convergence monitor.
(a)
(b)The trend of the curve is well reproduced by the simulations even if the absolute values result is underestimated. This can be explained considering the geometrical simplification (the domain was modeled as symmetric despite the cantilevered shafts which are not modeled).Figure5 shows the pressure evolution on the flanks during the engagement. The pressure on the mating flanks increases gradually when the gap between the teeth reduces the volume and squeezes the lubricant generating an axial flow (squeezing/pocketing effect). After the contact point passed through the pitch point, the volume increases again generating a negative pressure and a suction effect that is visible in Figure 6(b) looking to the axial velocity U
z.Figure 5
Relative pressure variation in meshing positions [Pa].
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Figure 6
(a) Path-lines (U [m/s]) and (b, c) detail (U
z [m/s]).
(a)
(b)
(c)These effects that can be evaluated by looking to the path-lines (Figure7) can also be quantified in terms of torque. In a previous study [23] the authors have used the different CFD approaches to calculate the power losses generated by a single rotating gear. The aim of that study was to validate different approaches to reduce the computational effort. The results of that study can be used for a comparison between the losses generated by 2 mating gears and a single rotating gear multiplied by 2.Figure 7
Detail of the contact.Figure8 represents the differences between the torque generated by two mating gears and that generated by one single rotating gear (×2) that is attributable to the squeezing effects. Considering the measurement accuracy and the big simplifications introduced in the numerical model, even if differences in the squeezing torque can be seen with respect to the measured values, the presented approach seems to be able to predict also this phenomena with good accuracy and therefore can be elected as suitable for the prediction of each source of load independent losses of gears.Figure 8
(a) Resistant torque generated by 2 mating gears and by 1 gear (×2): experimental and numerical values and (b) squeezing torque evaluated as the difference between the torque generated by 2 engaging gears and a single rotating gear (×2): experimental and numerical values.
(a)
(b)Regarding computational effort the solution time of each simulation was reduced from approximately 20 h on 16 3.3 GHz CPUs (211 GFLOPS) needed in a general-purpose commercial environment, to 21 h on a single 3.2 GHz CPU (12.8 GFLOPS) with a net reduction of the computational time of 93.5% (considering a perfect scalability). This is mainly due to two aspects. The first one is the reduced time for the mesh-handling and regeneration. The solution of the mass and momentum conservation equations represent just a small part of the total computational time and the presented approach takes less than 1/12th of the time needed by the commercial solver for the mesh generation. The second can be found in the capability of OpenFOAM to reach, with this specific mesh-handling technique/mesh type, convergence also for bigger time steps ensuring a significant reduction of the number of time steps needed to reach the regime condition. The commercial code remeshing algorithm substitutes the original distorted cells with tetrahedrons that can have also very small sizes affecting the maximum possible time step (considering to limit the Co to 1) leading to meshes that are bigger in size with respect to the initial one.


## 8. Conclusions

The results of the experiments confirm that CFD represents a valid method to predict power losses. The error in the predictions with respect to the measured ones is less than 10%. This little difference can be explained with the fact that the geometry was simplified. Simulations with and without the engagement show the capability of the model to predict also the squeezing effects.The main limitation of the wide diffusion of such an approach is traditionally the big computational effort both for modeling and for solving and the license costs. These drawbacks were overcome with the development of specific algorithms capable of handling the geometry generation, the boundary motion, and effective mesh generation and handling. This ensured a reduction of the computational effort of about 93.5%. In other words, with this approach it is possible to perform 15 times more simulations with the same computational effort making such technology suitable for industrial applications.A future work will be the systematic application of this approach for a better understanding of the power loss phenomena and the influence of the different operating conditions.The method can be of course generalized. The first simulations on a planetary configuration [21, 22] have shown good agreement with the experimental evidences and a net computational time reduction of approximately 93% with respect to the simulation performed by the authors with commercial software. For this reason the author maintains that this approach has a big potential and can be successfully extended to each gear geometry and possible gearbox configuration showing the same benefits highlighted with the presented example.



---

*Source: 2957151-2016-04-03.xml*