# CFD Simulations of Splash Losses of a Gearbox

Authors: Carlo Gorla; Franco Concli; Karsten Stahl; Bernd-Robert Höhn; Michaelis Klaus; Hansjörg Schultheiß; Johann-Paul Stemplinger
Journal: Advances in Tribology
 (2012)

---

## Abstract

Efficiency is becoming a main concern in the design of power transmissions. It is therefore important, especially during the design phase, to have appropriate models to predict the power losses. For this reason, CFD (computational fluid dynamics) simulations were performed in order to understand the influence of geometrical and operating parameters on the losses in power transmissions. The results of the model were validated with experimental results.


---

## Body

## 1. Introduction

Efficiency is becoming more and more a main concern in the design of power transmissions, and appropriate models to predict power losses are fundamental in order to reduce them, starting from the earliest stages of the design phase. Power losses of gearboxes are generally classified according to [1], taking into account the machine elements which are responsible for them and their dependency or nondependency from the load. Gear power losses are strongly related to lubrication, with those load dependent coming from the frictional effects in the lubricant film and those load independent mainly deriving from splashing, churning, and windage effects.Some models, obtained on the basis of experimental tests, can be found in the literature which describe the influence of gear geometric and kinematic parameters on hydraulic losses like for instance those proposed by Mauz [2], who has concentrated on hydraulic losses, or by Dowson [3], who has concentrated on windage losses.Nevertheless, the authors maintain that a deeper understanding of the physical phenomena responsible of gear losses is still needed in order to improve existing models, and CFD simulation can be an effective approach for such investigation.Marchesse et al. [4], on the basis of a state-of-the art on the application of CFD to gear power losses, applied CFD models to study windage losses of gears, or have validated their results by means of experimental tests.A large amount of experimental data on splashing power losses is still available at FZG, resulting from several years of tests [5], with either single discs, single gears and two meshing gears immersed in oil, and this data covers the influence of several parameters like outside diameter, face width, helix angle, temperature, and oil type. In order to improve the understanding of the mechanisms involved in splashing losses, CFD simulations have been run in order to investigate the effect of the same parameters, and the numerical results have been compared with the available test results.In the first phase of the activities, which is presented in this paper, single discs and gears have been considered. The quite good accordance between CFD simulation and experimental tests has confirmed that CFD represents an effective approach to study splashing power losses, and a simulation program with two meshing gears has therefore been started and is under course.


## 2. Composition of Gearbox Power Losses

According to [1], the power losses of gears can be subdivided into load-dependent and load-independent losses:
(1)PV=PVG+PVG0+PVB+PVB0+PVD+PVX.
The total power losses can be further subdivided according to their origin: the subscripts  G,  B,  D, and X are related to gears, to bearings, to contact seals, and to other factors respectively. The subscript 0 indicates load-independent power losses.P
V
G  are the load-dependent power losses of gears and arise primarily from the sliding between the flanks of the gear teeth.P
V
B  are the load-dependent power losses of bearings and are also related to sliding between the rolling elements and the rings.P
V
B
0  are the load-independent power losses of bearings and are, among others, related to viscous effects due to the lubrication.P
V
D  are the power losses of the seals and are related to sliding.P
V
X  are other generic power losses.P
V
G
0  are the load-independent power losses of gears due to the interaction between the lubricant and the rotating/moving elements. These losses are the sum of squeezing, power, and splashing power losses. The lubricant squeezing power losses are related to the fact that the gap at the mesh position is changing its volume during the engagement causing an overpressure that squeezes the oil primarily in the axial direction [6–8]. The power and splashing power losses are related to the viscous and pressure effects of the lubricant on the moving/rotating elements [9, 10]. For dip-lubricated power transmissions, this kind of losses cannot be neglected and are an important part of the total losses. For this reason, the authors have studied this kind of losses for a simple geometry and under different operating conditions.


## 3. Problem Description

In order to understand the influence of the different parameters, the authors have performed some initial simulation with a simplified geometry, consisting of a simple rotating disk in a fully filled case (no free surface is present, both in the experiments and the CFD approach). The presence of one single phase leads to classical mechanical losses generated by moving fluid around a solid. For example, when the fluid is air (or air/oil mixture), these losses are called windage, and the losses approximately evolve with angular velocity power 3. A schematic layout of the analyzed gearbox is shown in Figure1. The disk is mounted on a cantilevered shaft that is supported by two bearings. A second over-hung mounted shaft is placed in the case parallel to the primary shaft. In a second steps a real gear geometry replaced the disk.Figure 1
Geometry of the gearbox with the gear.The simulations were conducted under pressure while the gearbox was developed to operate on the bottom of the sea and this pressure is needed in order to compensate the external water pressure at operating depth.


### 3.1. Oil Properties

In these investigations, different oils were used. The lubricants are mineral based in different viscosity grades. An additional virtual oil FVA3*, with same viscosities as FVA3 and same density as FVA2, was used to evaluate the influence of viscosity on power losses. The detailed oil properties can be seen in Table1.Table 1
Oil properties of the used lubricants.
Oil
Kinematic viscosity at 40°C ν40 (mm²/s)
Kinematic viscosity at 100°C ν100 (mm²/s)
Density at 15°C ρ15 (kg/m³)
FVA2
29.8
5.2
871
FVA3
95.0
10.7
885
FVA3*
95.0
10.7
871


### 3.2. Geometry

The geometry for the numerical model consists only of the oil volume. As shown in Figure2, only a part of the case is modelled as also the bearings were neglected.Figure 2
Geometry for the numerical model.Both disk and gear are not mounted in the shaft’s center; therefore, it is not possible to take advantage of symmetry.The first analyzed component is a disk; its properties are summarized in Table2.Table 2
Geometrical properties of the disk.
Diameterda (mm)
100
Widthb (mm)
40Table3 shows the geometrical parameters for the gear adopted in the different simulations.Table 3
Geometrical properties of the analyzed gears.
Widthb (mm)
Tip diameterda (mm)
Helix angleβ (°)
Pressure angleαn (°)
Number of teethz1()
Normal modulsmn (mm)
Reference case
40
102.5
0
20
23
4
Width influence
20
102.5
0
20
23
4
Tip diameter influence
40
96.5÷ 98
0
20
23
4
Helix angle influence
40
102.5
20
20
23
4


### 3.3. Mesh

For the disk’s case, a simple undeforming mesh was assumed. In order to correctly reproduce the motion of the boundaries, appropriate BC (boundary conditions) were set. Tetrahedral elements have been used. For complex geometries, in fact, the tetrahedral mesh can often be created with far fewer cells than the equivalent mesh consisting of quadrilateral/hexahedral elements. This is because the tetrahedral mesh allows clustering of cells in selected regions of the flow domain. Structured hexahedral meshes will generally force cells to be placed in regions where they are not needed.Tetrahedral cells are not desirable near walls if the boundary layer needs to be resolved because the first grid point must be very close to the wall, while relatively large grid sizes can be used in the directions parallel to the wall. These requirements lead to long thin tetrahedral elements, creating problems in the approximation of diffusive fluxes. For this reason, some grid generation methods generate first a layer of prisms or hexahedra near solid boundaries (inflation), starting with a triangular or quadrilateral discretization of the surface; on top of this layer, a tetrahedral mesh is generated automatically in the remaining part of the domain (Figure3).Figure 3
Detail of the mesh near the wall.For the simulations with the gear, an additional feature was used: in order to model the rotation of the gear without mesh deformation, two separate domains were created and meshed separately. The region around the gear is a cylinder, which was defined as dynamic, it rotates around the gear axis. The other region, corresponding to the rest of the domain, was defined as a static zone. This zone has a cylindrical cavity, in which the dynamic mesh is located (as shown in Figure4).Figure 4
Partition of the mesh: the dynamic zone is marked with light blue, and the static zone is marked with light grey.The two zones have some faces that are intersecting each other but from a numerical point of view they are not connected. In order to link them, an interface zone is defined. The connection ensures that the values of a generic field during the simulations are the same on both sides of the interface.Due to the rotation of one of the two zones, an additional feature is necessary, sliding mesh. The two cell zones will move relative to each other along the grid interface. The sliding mesh condition does not move the fluid directly but moves the mesh only, allowing to set the appropriate motion of the boundaries inside the dynamic zone, in this case the faces of the gear.


### 3.4. Solver Settings

#### 3.4.1. Navier Stokes Equations

The CFD is based on some differential equations. For a generic element, it is possible to write five equations. The first equation is the averaged mass conservation equation for no-stationary incompressible flows and can be written as follows:(2)∂〈ρui〉∂xi=0,
where xi is the Cartesian coordinate and ui  is the velocity component.The second equation is the averaged momentum conservation equation and can be written as follows:(3)∂∂t(ρ〈ui〉)+∂∂xj(ρ〈ui〉〈ui〉)=-∂∂xi〈p〉+∂∂xj[μ(∂〈ui〉∂xj+∂〈uj〉∂xi)-∂∂xjτij],
where  p is the pressure, xi and xj are the Cartesian coordinates, and ρ is the density. The additional term τij that appears in the averaged equation (in comparison to the nonaveraged transport equations) is called unresolved term or Reynolds term. The averaging process produces a set of equations that is not closed. For this reason a turbulence model is needed in order to be able to solve the system of equations.The unresolved or Reynolds term is expressed by using the eddy-viscosity hypothesis as follows:(4)-τij=-ρ〈ui′uj′〉=μt(∂〈ui〉∂xj+∂〈uj〉∂xi)-23ρδijk,
where  k is the turbulent kinetic energy and μt the eddy viscosity.The energy equation was not activated in the given model: the operating temperature was defined as a priori and consequently the properties of the fluid do not change during the calculations.For the pressure-velocity coupling a simple scheme was adopted as suggested for flows in closed domains [11].


#### 3.4.2. Turbulence Models

Turbulent flows are characterized by fluctuating velocity fields. These fluctuations mix transported quantities such as momentum, energy, and species concentration, cause the transported quantities to fluctuate as well. Since these fluctuations can be of small scale and high frequency, they are too much computationally expensive to simulate directly in practical engineering calculations. Instead, the instantaneous exact governing equations can be averaged to remove the small scales, resulting in a modified set of equations that are computationally less expensive to solve. However, the modified equations contain additional unknown variables, and turbulence models are needed to determine these variables in terms of known quantities.The simplest “complete models” of turbulence are two-equation models in which the solution of two separate transport equations allows the turbulent velocity and length scales to be independently determined. This model is a semiempirical model based on model transport equations for the turbulence kinetic energy (k) and its dissipation rate (ε) [12].The RNGk-ε  model assumes that the eddy viscosity is related to the turbulence kinetic energy and dissipation via the following relation:
(5)μt=Cμρk2ε.
The model transport equation for k is derived from the exact equation, while the model transport equation for ε was obtained using physical reasoning and bears little resemblance to its mathematically exact counterpart [11]. In the derivation of the k-ε  model, the assumption is that the flow is fully turbulent, and the effects of molecular viscosity are negligible:
(6)∂∂t(ρk)+∂∂xi(ρkui)=∂∂xj[αkμeff∂k∂xj]+Gk+Gb-ρε-YM+Sk∂∂t(ρε)+∂∂xi(ρεui)=∂∂xj[αεμeff∂k∂xj]+C1εεk(Gk+C3εGb)-C2ερε2k-Rε+Sε.
In these equations, Gk represents the generation of turbulence kinetic energy due to the mean velocity gradients. Gb is the generation of turbulence kinetic energy due to buoyancy effect. YM represents the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate. C1ε,  C2ε, and C3ε are constants. σk and σε are the turbulent Prandtl numbers for k and ε, respectively. The quantities ∝k and ∝ε are the inverse effective Prandtl numbers for k and ε, respectively. Sk and Sε are user-defined source terms. Rε is the additional term that characterizes the RNG-k-ε-model. Rε is derived from a rigorous statistical technique called renormalization group theory and improves the accuracy of the standard-k-ε-model.Some other more complex turbulent models exist, but some workers claim that this model offers improved accuracy in rotating flows [13]. For this reason, this turbulence model was selected to perform the simulations.


#### 3.4.3. Spatial Discretization

By default, the solver stores discrete values of the scalarϕ  at the cell centers. However, face values ϕf  are required for the convection terms and must be interpolated from the cell center values. This is accomplished using an upwind scheme. Upwinding means that the face value ϕf  is derived from quantities in the cell upstream, or “upwind,” relative to the direction of the normal velocity vn.In order to solve the differential equations, the authors have chosen to use second-order upwind schemes. When second-order upwinding is selected, the face valueϕf  is computed on the base of the cell-centered averaged value  ϕ and its gradient Δϕ. For the determination of the gradient, a least squares cell-based evaluation was used. In this method, the solution is assumed to vary linearly between two cell centroids.


#### 3.4.4. Boundary Conditions

All the boundaries were set to a no slip condition. That means that on the wall there is no relative velocity between the walls and the fluid. To describe properly the velocity profile in the normal direction, an enhanced wall treatment was used. This technique calculates the velocity at the elements’ center and then reconstructs the velocity profile starting from these quantities. For this reason, it is important to create an appropriate fine mesh near the walls.In order to evaluate the quality of the results in terms of velocity profile one has to check they+ value that is defined as follows:
(7)y+=u*yν,
and—in case of enhanced wall treatment—should be approximately equal to one. u*  is the friction velocity at the nearest wall, y  is the distance from the center of the first cells to the wall and ν is the local kinematic viscosity of the fluid (Figure 5).Figure 5
Velocity profile near the wall.


#### 3.4.5. Time Steps

The determination of the time step size is based on the estimation of the truncation error associated with the time integration scheme. If the truncation error is smaller than a specified tolerance, the size of the time step is increased; if the truncation error is greater, the time step size is decreased. An estimation of the truncation error can be obtained by using a predictor-corrector type of algorithm according to [14] in association with the time integration scheme. At each time step, a predicted solution can be obtained using a computationally inexpensive explicit method (Adams-Bashford for the second-order unsteady formulation). This predicted solution is used as an initial condition for the time step, and the correction is computed using the nonlinear iterations associated with the implicit formulation. The norm of the difference between the predicted and corrected solutions is used as a measure of the truncation error. By comparing the truncation error with the desired level of accuracy the code, is able to adjust the time step size by increasing it or decreasing it.


## 3.1. Oil Properties

In these investigations, different oils were used. The lubricants are mineral based in different viscosity grades. An additional virtual oil FVA3*, with same viscosities as FVA3 and same density as FVA2, was used to evaluate the influence of viscosity on power losses. The detailed oil properties can be seen in Table1.Table 1
Oil properties of the used lubricants.
Oil
Kinematic viscosity at 40°C ν40 (mm²/s)
Kinematic viscosity at 100°C ν100 (mm²/s)
Density at 15°C ρ15 (kg/m³)
FVA2
29.8
5.2
871
FVA3
95.0
10.7
885
FVA3*
95.0
10.7
871


## 3.2. Geometry

The geometry for the numerical model consists only of the oil volume. As shown in Figure2, only a part of the case is modelled as also the bearings were neglected.Figure 2
Geometry for the numerical model.Both disk and gear are not mounted in the shaft’s center; therefore, it is not possible to take advantage of symmetry.The first analyzed component is a disk; its properties are summarized in Table2.Table 2
Geometrical properties of the disk.
Diameterda (mm)
100
Widthb (mm)
40Table3 shows the geometrical parameters for the gear adopted in the different simulations.Table 3
Geometrical properties of the analyzed gears.
Widthb (mm)
Tip diameterda (mm)
Helix angleβ (°)
Pressure angleαn (°)
Number of teethz1()
Normal modulsmn (mm)
Reference case
40
102.5
0
20
23
4
Width influence
20
102.5
0
20
23
4
Tip diameter influence
40
96.5÷ 98
0
20
23
4
Helix angle influence
40
102.5
20
20
23
4


## 3.3. Mesh

For the disk’s case, a simple undeforming mesh was assumed. In order to correctly reproduce the motion of the boundaries, appropriate BC (boundary conditions) were set. Tetrahedral elements have been used. For complex geometries, in fact, the tetrahedral mesh can often be created with far fewer cells than the equivalent mesh consisting of quadrilateral/hexahedral elements. This is because the tetrahedral mesh allows clustering of cells in selected regions of the flow domain. Structured hexahedral meshes will generally force cells to be placed in regions where they are not needed.Tetrahedral cells are not desirable near walls if the boundary layer needs to be resolved because the first grid point must be very close to the wall, while relatively large grid sizes can be used in the directions parallel to the wall. These requirements lead to long thin tetrahedral elements, creating problems in the approximation of diffusive fluxes. For this reason, some grid generation methods generate first a layer of prisms or hexahedra near solid boundaries (inflation), starting with a triangular or quadrilateral discretization of the surface; on top of this layer, a tetrahedral mesh is generated automatically in the remaining part of the domain (Figure3).Figure 3
Detail of the mesh near the wall.For the simulations with the gear, an additional feature was used: in order to model the rotation of the gear without mesh deformation, two separate domains were created and meshed separately. The region around the gear is a cylinder, which was defined as dynamic, it rotates around the gear axis. The other region, corresponding to the rest of the domain, was defined as a static zone. This zone has a cylindrical cavity, in which the dynamic mesh is located (as shown in Figure4).Figure 4
Partition of the mesh: the dynamic zone is marked with light blue, and the static zone is marked with light grey.The two zones have some faces that are intersecting each other but from a numerical point of view they are not connected. In order to link them, an interface zone is defined. The connection ensures that the values of a generic field during the simulations are the same on both sides of the interface.Due to the rotation of one of the two zones, an additional feature is necessary, sliding mesh. The two cell zones will move relative to each other along the grid interface. The sliding mesh condition does not move the fluid directly but moves the mesh only, allowing to set the appropriate motion of the boundaries inside the dynamic zone, in this case the faces of the gear.


## 3.4. Solver Settings

### 3.4.1. Navier Stokes Equations

The CFD is based on some differential equations. For a generic element, it is possible to write five equations. The first equation is the averaged mass conservation equation for no-stationary incompressible flows and can be written as follows:(2)∂〈ρui〉∂xi=0,
where xi is the Cartesian coordinate and ui  is the velocity component.The second equation is the averaged momentum conservation equation and can be written as follows:(3)∂∂t(ρ〈ui〉)+∂∂xj(ρ〈ui〉〈ui〉)=-∂∂xi〈p〉+∂∂xj[μ(∂〈ui〉∂xj+∂〈uj〉∂xi)-∂∂xjτij],
where  p is the pressure, xi and xj are the Cartesian coordinates, and ρ is the density. The additional term τij that appears in the averaged equation (in comparison to the nonaveraged transport equations) is called unresolved term or Reynolds term. The averaging process produces a set of equations that is not closed. For this reason a turbulence model is needed in order to be able to solve the system of equations.The unresolved or Reynolds term is expressed by using the eddy-viscosity hypothesis as follows:(4)-τij=-ρ〈ui′uj′〉=μt(∂〈ui〉∂xj+∂〈uj〉∂xi)-23ρδijk,
where  k is the turbulent kinetic energy and μt the eddy viscosity.The energy equation was not activated in the given model: the operating temperature was defined as a priori and consequently the properties of the fluid do not change during the calculations.For the pressure-velocity coupling a simple scheme was adopted as suggested for flows in closed domains [11].


### 3.4.2. Turbulence Models

Turbulent flows are characterized by fluctuating velocity fields. These fluctuations mix transported quantities such as momentum, energy, and species concentration, cause the transported quantities to fluctuate as well. Since these fluctuations can be of small scale and high frequency, they are too much computationally expensive to simulate directly in practical engineering calculations. Instead, the instantaneous exact governing equations can be averaged to remove the small scales, resulting in a modified set of equations that are computationally less expensive to solve. However, the modified equations contain additional unknown variables, and turbulence models are needed to determine these variables in terms of known quantities.The simplest “complete models” of turbulence are two-equation models in which the solution of two separate transport equations allows the turbulent velocity and length scales to be independently determined. This model is a semiempirical model based on model transport equations for the turbulence kinetic energy (k) and its dissipation rate (ε) [12].The RNGk-ε  model assumes that the eddy viscosity is related to the turbulence kinetic energy and dissipation via the following relation:
(5)μt=Cμρk2ε.
The model transport equation for k is derived from the exact equation, while the model transport equation for ε was obtained using physical reasoning and bears little resemblance to its mathematically exact counterpart [11]. In the derivation of the k-ε  model, the assumption is that the flow is fully turbulent, and the effects of molecular viscosity are negligible:
(6)∂∂t(ρk)+∂∂xi(ρkui)=∂∂xj[αkμeff∂k∂xj]+Gk+Gb-ρε-YM+Sk∂∂t(ρε)+∂∂xi(ρεui)=∂∂xj[αεμeff∂k∂xj]+C1εεk(Gk+C3εGb)-C2ερε2k-Rε+Sε.
In these equations, Gk represents the generation of turbulence kinetic energy due to the mean velocity gradients. Gb is the generation of turbulence kinetic energy due to buoyancy effect. YM represents the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate. C1ε,  C2ε, and C3ε are constants. σk and σε are the turbulent Prandtl numbers for k and ε, respectively. The quantities ∝k and ∝ε are the inverse effective Prandtl numbers for k and ε, respectively. Sk and Sε are user-defined source terms. Rε is the additional term that characterizes the RNG-k-ε-model. Rε is derived from a rigorous statistical technique called renormalization group theory and improves the accuracy of the standard-k-ε-model.Some other more complex turbulent models exist, but some workers claim that this model offers improved accuracy in rotating flows [13]. For this reason, this turbulence model was selected to perform the simulations.


### 3.4.3. Spatial Discretization

By default, the solver stores discrete values of the scalarϕ  at the cell centers. However, face values ϕf  are required for the convection terms and must be interpolated from the cell center values. This is accomplished using an upwind scheme. Upwinding means that the face value ϕf  is derived from quantities in the cell upstream, or “upwind,” relative to the direction of the normal velocity vn.In order to solve the differential equations, the authors have chosen to use second-order upwind schemes. When second-order upwinding is selected, the face valueϕf  is computed on the base of the cell-centered averaged value  ϕ and its gradient Δϕ. For the determination of the gradient, a least squares cell-based evaluation was used. In this method, the solution is assumed to vary linearly between two cell centroids.


### 3.4.4. Boundary Conditions

All the boundaries were set to a no slip condition. That means that on the wall there is no relative velocity between the walls and the fluid. To describe properly the velocity profile in the normal direction, an enhanced wall treatment was used. This technique calculates the velocity at the elements’ center and then reconstructs the velocity profile starting from these quantities. For this reason, it is important to create an appropriate fine mesh near the walls.In order to evaluate the quality of the results in terms of velocity profile one has to check they+ value that is defined as follows:
(7)y+=u*yν,
and—in case of enhanced wall treatment—should be approximately equal to one. u*  is the friction velocity at the nearest wall, y  is the distance from the center of the first cells to the wall and ν is the local kinematic viscosity of the fluid (Figure 5).Figure 5
Velocity profile near the wall.


### 3.4.5. Time Steps

The determination of the time step size is based on the estimation of the truncation error associated with the time integration scheme. If the truncation error is smaller than a specified tolerance, the size of the time step is increased; if the truncation error is greater, the time step size is decreased. An estimation of the truncation error can be obtained by using a predictor-corrector type of algorithm according to [14] in association with the time integration scheme. At each time step, a predicted solution can be obtained using a computationally inexpensive explicit method (Adams-Bashford for the second-order unsteady formulation). This predicted solution is used as an initial condition for the time step, and the correction is computed using the nonlinear iterations associated with the implicit formulation. The norm of the difference between the predicted and corrected solutions is used as a measure of the truncation error. By comparing the truncation error with the desired level of accuracy the code, is able to adjust the time step size by increasing it or decreasing it.


## 3.4.1. Navier Stokes Equations

The CFD is based on some differential equations. For a generic element, it is possible to write five equations. The first equation is the averaged mass conservation equation for no-stationary incompressible flows and can be written as follows:(2)∂〈ρui〉∂xi=0,
where xi is the Cartesian coordinate and ui  is the velocity component.The second equation is the averaged momentum conservation equation and can be written as follows:(3)∂∂t(ρ〈ui〉)+∂∂xj(ρ〈ui〉〈ui〉)=-∂∂xi〈p〉+∂∂xj[μ(∂〈ui〉∂xj+∂〈uj〉∂xi)-∂∂xjτij],
where  p is the pressure, xi and xj are the Cartesian coordinates, and ρ is the density. The additional term τij that appears in the averaged equation (in comparison to the nonaveraged transport equations) is called unresolved term or Reynolds term. The averaging process produces a set of equations that is not closed. For this reason a turbulence model is needed in order to be able to solve the system of equations.The unresolved or Reynolds term is expressed by using the eddy-viscosity hypothesis as follows:(4)-τij=-ρ〈ui′uj′〉=μt(∂〈ui〉∂xj+∂〈uj〉∂xi)-23ρδijk,
where  k is the turbulent kinetic energy and μt the eddy viscosity.The energy equation was not activated in the given model: the operating temperature was defined as a priori and consequently the properties of the fluid do not change during the calculations.For the pressure-velocity coupling a simple scheme was adopted as suggested for flows in closed domains [11].


## 3.4.2. Turbulence Models

Turbulent flows are characterized by fluctuating velocity fields. These fluctuations mix transported quantities such as momentum, energy, and species concentration, cause the transported quantities to fluctuate as well. Since these fluctuations can be of small scale and high frequency, they are too much computationally expensive to simulate directly in practical engineering calculations. Instead, the instantaneous exact governing equations can be averaged to remove the small scales, resulting in a modified set of equations that are computationally less expensive to solve. However, the modified equations contain additional unknown variables, and turbulence models are needed to determine these variables in terms of known quantities.The simplest “complete models” of turbulence are two-equation models in which the solution of two separate transport equations allows the turbulent velocity and length scales to be independently determined. This model is a semiempirical model based on model transport equations for the turbulence kinetic energy (k) and its dissipation rate (ε) [12].The RNGk-ε  model assumes that the eddy viscosity is related to the turbulence kinetic energy and dissipation via the following relation:
(5)μt=Cμρk2ε.
The model transport equation for k is derived from the exact equation, while the model transport equation for ε was obtained using physical reasoning and bears little resemblance to its mathematically exact counterpart [11]. In the derivation of the k-ε  model, the assumption is that the flow is fully turbulent, and the effects of molecular viscosity are negligible:
(6)∂∂t(ρk)+∂∂xi(ρkui)=∂∂xj[αkμeff∂k∂xj]+Gk+Gb-ρε-YM+Sk∂∂t(ρε)+∂∂xi(ρεui)=∂∂xj[αεμeff∂k∂xj]+C1εεk(Gk+C3εGb)-C2ερε2k-Rε+Sε.
In these equations, Gk represents the generation of turbulence kinetic energy due to the mean velocity gradients. Gb is the generation of turbulence kinetic energy due to buoyancy effect. YM represents the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate. C1ε,  C2ε, and C3ε are constants. σk and σε are the turbulent Prandtl numbers for k and ε, respectively. The quantities ∝k and ∝ε are the inverse effective Prandtl numbers for k and ε, respectively. Sk and Sε are user-defined source terms. Rε is the additional term that characterizes the RNG-k-ε-model. Rε is derived from a rigorous statistical technique called renormalization group theory and improves the accuracy of the standard-k-ε-model.Some other more complex turbulent models exist, but some workers claim that this model offers improved accuracy in rotating flows [13]. For this reason, this turbulence model was selected to perform the simulations.


## 3.4.3. Spatial Discretization

By default, the solver stores discrete values of the scalarϕ  at the cell centers. However, face values ϕf  are required for the convection terms and must be interpolated from the cell center values. This is accomplished using an upwind scheme. Upwinding means that the face value ϕf  is derived from quantities in the cell upstream, or “upwind,” relative to the direction of the normal velocity vn.In order to solve the differential equations, the authors have chosen to use second-order upwind schemes. When second-order upwinding is selected, the face valueϕf  is computed on the base of the cell-centered averaged value  ϕ and its gradient Δϕ. For the determination of the gradient, a least squares cell-based evaluation was used. In this method, the solution is assumed to vary linearly between two cell centroids.


## 3.4.4. Boundary Conditions

All the boundaries were set to a no slip condition. That means that on the wall there is no relative velocity between the walls and the fluid. To describe properly the velocity profile in the normal direction, an enhanced wall treatment was used. This technique calculates the velocity at the elements’ center and then reconstructs the velocity profile starting from these quantities. For this reason, it is important to create an appropriate fine mesh near the walls.In order to evaluate the quality of the results in terms of velocity profile one has to check they+ value that is defined as follows:
(7)y+=u*yν,
and—in case of enhanced wall treatment—should be approximately equal to one. u*  is the friction velocity at the nearest wall, y  is the distance from the center of the first cells to the wall and ν is the local kinematic viscosity of the fluid (Figure 5).Figure 5
Velocity profile near the wall.


## 3.4.5. Time Steps

The determination of the time step size is based on the estimation of the truncation error associated with the time integration scheme. If the truncation error is smaller than a specified tolerance, the size of the time step is increased; if the truncation error is greater, the time step size is decreased. An estimation of the truncation error can be obtained by using a predictor-corrector type of algorithm according to [14] in association with the time integration scheme. At each time step, a predicted solution can be obtained using a computationally inexpensive explicit method (Adams-Bashford for the second-order unsteady formulation). This predicted solution is used as an initial condition for the time step, and the correction is computed using the nonlinear iterations associated with the implicit formulation. The norm of the difference between the predicted and corrected solutions is used as a measure of the truncation error. By comparing the truncation error with the desired level of accuracy the code, is able to adjust the time step size by increasing it or decreasing it.


## 4. Operating Conditions

With CFD analysis, the effect of different parameters was investigated. For the simple case with the disk instead of the gear, the only varied parameter was the velocity in a range between 500 and 8000 rpm (Table4).Table 4
combination of parameters adopted in the different simulations.
Disk or Gear
d
a (mm)
b (mm)
β (°)
n (rpm)
p (bar)
T (°C)
Oil level (%)
Oil
A
D
100
40
/
500
6
40
100
FVA2
B
D
100
40
/
2000
6
40
100
FVA2
C
D
100
40
/
5000
6
40
100
FVA2
D
D
100
40
/
8000
6
40
100
FVA2
1.1
G
102.5
40
0
1500
6
88
*
100
FVA2
1.2
G
102.5
40
0
3000
6
84
*
100
FVA2
1.3
G
102.5
40
0
4500
6
85
*
100
FVA2
1.4
G
102.5
40
0
6000
6
93
*
100
FVA2
1.5
G
102.5
40
0
7500
6
107
*
100
FVA2
2.1
G
102.5
40
0
500
6
90
100
FVA2
2.2
G
102.5
40
0
2000
6
90
100
FVA2
2.3
G
102.5
40
0
5000
6
90
100
FVA2
2.4
G
102.5
40
0
8000
6
90
100
FVA2
3.1
G
102.5
20
0
500
6
90
100
FVA2
3.2
G
102.5
20
0
2000
6
90
100
FVA2
3.3
G
102.5
20
0
5000
6
90
100
FVA2
3.4
G
102.5
20
0
8000
6
90
100
FVA2
4.1
G
98
40
0
500
6
89
+
100
FVA2
4.2
G
98
40
0
2000
6
89
+
100
FVA2
4.3
G
98
40
0
5000
6
89
+
100
FVA2
4.4
G
98
40
0
8000
6
109
+
100
FVA2
5.1
G
98
40
0
500
6
90
100
FVA2
5.2
G
98
40
0
2000
6
90
100
FVA2
5.3
G
98
40
0
5000
6
90
100
FVA2
5.4
G
98
40
0
8000
6
90
100
FVA2
6.1
G
96.5
40
0
500
6
90
100
FVA2
6.2
G
96.5
40
0
2000
6
90
100
FVA2
6.3
G
96.5
40
0
5000
6
90
100
FVA2
6.4
G
96.5
40
0
8000
6
90
100
FVA2
7.1
G
102.5
40
20
500
6
90
100
FVA2
7.2
G
102.5
40
20
2000
6
90
100
FVA2
7.3
G
102.5
40
20
5000
6
90
100
FVA2
7.4
G
102.5
40
20
8000
6
90
100
FVA3
8.1
G
102.5
40
0
500
6
90
100
FVA3
8.2
G
102.5
40
0
2000
6
90
100
FVA3
8.3
G
102.5
40
0
5000
6
90
100
FVA3
8.4
G
102.5
40
0
8000
6
90
100
FVA3
9.1
G
102.5
40
0
500
6
90
100
FVA3
*
9.2
G
102.5
40
0
2000
6
90
100
FVA3
*
9.3
G
102.5
40
0
5000
6
90
100
FVA3
*
9.4
G
102.5
40
0
8000
6
90
100
FVA3
*
+Measured temperature.
*Oil with the viscosity of FVA3 and the density of FVA2.For the case with the gear, the investigated parameters were the velocity, the pressure, and the geometry (face width, tip diameter, and helix angle).Table4 summarizes the combination of parameters adopted in the different simulations.The real gearbox used for the experiments was developed to operate on the bottom of the sea and an internal pressure (6 bar) was needed in order to compensate the external water pressure at operating depth.The experiments have been therefore performed with the operating conditions of 6 bar. The pressurization of the gearbox was made after some experiments with an operating pressure of 0 bar. The choice to operate with 6 bar was taken in order to be sure to have only oil in the gearbox, avoiding the presence of air bubbles.Moreover, the level of the static pressure does not affect the results.


## 5. Experimental Tests

In order to validate the results of the CFD simulations, some experiments were performed. On the driving shaft, a torque meter was mounted. This instrument measures the resistant torque and the rotational velocity so it is possible to evaluate the power losses. Due to the configuration of the system, it is not possible to measure directly the load-independent power losses due to power. For this reason, some additional tests were performed without the gear/disk and with air instead of oil in the gearbox so as to be able to measure the power losses generated by the bearings. Making the difference between the total losses and the losses of the bearings, it is possible to evaluate the losses only due to the churning.


## 6. Results

### 6.1. Results: Disk

Figure6 shows the results in terms of no load loss torque PVZ0,Pl versus circumferential vt  velocity for the case with the disk. The continuous line represents the results of the experiments without the losses of the bearings and the dots the results of the simulations. In this simple case, the results of the simulations appear to underestimate the losses. It is however to say that the experimental results that appear in Figure 6 are obtained like the difference between the total measured losses in the working configuration (including therefore the bearing losses) and the losses measured with only the shaft (no disk) and air instead of oil in the gearbox. In this manner, it is possible to separate the power losses due to the rotation of the disk in the oil from the losses caused by the bearings. It is however to say that in this second measurement, the absence of the lubricant may minimally affect the losses caused by the bearings: the inner race of the bearings, in fact, during the rotation is laterally in contact with air instead of oil, and this leads to a minimal underestimation of the losses of the bearing. This problem becomes less significant in the measurements with the gear where the power losses are significantly higher.Figure 6
Results in terms of resistant torque versus circumferential velocity for the disk.Another reason for the little difference between the numerical and the experimental results can be found in the fact that in the simulations, the temperature (and therefore the oil properties) was fixed to 40°C while in the experiments, it was varying due to the fact that the cooling/heating system was not able to compensate immediately a change in the temperature.


### 6.2. Results: Gear

Figure7 shows the velocity field in the gear symmetry plane for four different tangential velocities. It can be seen that the area in which the velocity is appreciable is expanding with the increase of the tangential velocity.Figure 7
Contour plot of the velocity field for four different tangential velocities (m/s) in the symmetry plane of the gear: (a) 500 rpm; (b) 2000 rpm; (c) 5000 rpm; (d) 8000 rpm.From Figure8, it is possible to have an idea of the fluxes of lubricant inside the gearbox. An ideal particle of fluid escapes the gear radially and then, after a loop, it comes again in contact with the gear from the axial direction. Due to the rotation of the gear, there are no zones in the gearbox where the lubricant is steady but the complete domain is involved in the lubricant flux.Figure 8
Velocitystreamlines;b=40 mm; d=102.5 mm; T=90°C; n=8000 rpm.Figures9 to 13 show the results in terms of splashing resistant torque PVZ0,Pl versus circumferential velocity vt for the cases with the gear. The continuous line represents the results of the experiments without the losses of the bearings, the broken line the results of the calculations according to Mauz [2], and the dotted line the results of the simulations.Figure 9
Results in terms of resistant torque versus circumferential velocity—effect of the temperature.Mauz proposed some equations derived by some experiments. The model proposed by Mauz has the big advantage that it does not need computational resources but, as shown in the next diagrams, it has also some limits. For example, it is not valid for tangential velocities higher than 60 m/s and it does not take into account some parameters like, for example, the helix angle.Figure9 shows the influence of the temperature on the resistant torque. It can be seen that the influence of the temperature on the results cannot be neglected. The experimental tests were performed with a temperature that is not exactly constant, which is due to the cooling/heating system that is not able to compensate immediately a fluctuation of the temperature of the lubrication bath. Some simulations were performed with the same temperature measured in the experiments for each rotational speed: the lubricant properties have been changed in each simulation according to the measured temperature at the different rotational speeds. It can be seen that, the CFD results are in very good agreement with the measured ones. Some other calculations, instead, were performed with a constant temperature of 90°C. The CFD results for the constant temperature are below the measured ones only for low rotational speeds where the effective temperature of the oil bath in the test rig was lower than 90°C. For higher rotational speeds the CFD results calculated with a constant temperature of 90°C are over the measured ones where the temperature of the oil bath in the test rig was higher than 90°C.The diagram shows also the power losses according to Mauz. The results differ significantly from the measured and simulated values.Figure10 shows the similar results as Figure 8 but for another tip diameter. Also in this case, the results of the simulations performed with the same temperature of the experiments overlap the measured results.Figure 10
Results in terms of resistant torque versus circumferential velocity—effect of the temperature.Figure11 shows the effect of the tip diameter on the power losses. Simulations have been performed for a diameter of 102.5 mm, 98 mm, and 96.5 mm. Experiments have been conducted only for the biggest and the smallest diameter (see Figures 9 and 10).Figure 11
Results in terms of resistant torque versus circumferential velocity—influence of the tip diameter.The effect of the tip diameter on the power losses is extremely high. The results according to Mauz did not correctly describe this phenomenon. The results of the simulations, however, are able to well predict the power losses as confirmed by the experiments.An increase of the tip diameter from 96.5 mm to 98 mm causes an increase of the power losses of about 64%.Similarly, an increase of the tip diameter from 96.5 mm to 102.5 mm causes an increase of the power losses of about 142%.Figure12 shows the effect of the face width on the power losses.Figure 12
Results in terms of resistant torque versus circumferential velocity—Influence of the face width.Figure 13
Results in terms of resistant torque versus circumferential velocity—Influence of the helix angle.The effect of the face width on the power losses is extremely high. Also in this case the results of the simulations are better than the results obtained by Mauz even if in this case also for Mauz an increase of the face width of 100% gives a significant increase of the power losses.The doubling of the tip diameter from 20 mm to 40 mm causes an increase of the power losses of about 60%.Figure13 shows the effect of the helix angle on the power losses. It is possible to appreciate how the losses decrease with the helix angle.In particular, the losses decrease to about 45% changing the helix angle from 0° to 20°. The model proposed by Mauz does not take into account the helix angle.Figure14 shows the effect of the oil viscosity on the power losses, while Figure 15 shows the influence of the oil type (density) on the power losses.Figure 14
Results in terms of resistant torque versus circumferential velocity—Influence of the lubricant viscosity.Figure 15
Results in terms of resistant torque versus circumferential velocity—Influence of the lubricant density.Increasing the density from 819.5 Kg/m3 to 824.5 Kg/m3 leads to an increase of the power losses of about 15%.


## 6.1. Results: Disk

Figure6 shows the results in terms of no load loss torque PVZ0,Pl versus circumferential vt  velocity for the case with the disk. The continuous line represents the results of the experiments without the losses of the bearings and the dots the results of the simulations. In this simple case, the results of the simulations appear to underestimate the losses. It is however to say that the experimental results that appear in Figure 6 are obtained like the difference between the total measured losses in the working configuration (including therefore the bearing losses) and the losses measured with only the shaft (no disk) and air instead of oil in the gearbox. In this manner, it is possible to separate the power losses due to the rotation of the disk in the oil from the losses caused by the bearings. It is however to say that in this second measurement, the absence of the lubricant may minimally affect the losses caused by the bearings: the inner race of the bearings, in fact, during the rotation is laterally in contact with air instead of oil, and this leads to a minimal underestimation of the losses of the bearing. This problem becomes less significant in the measurements with the gear where the power losses are significantly higher.Figure 6
Results in terms of resistant torque versus circumferential velocity for the disk.Another reason for the little difference between the numerical and the experimental results can be found in the fact that in the simulations, the temperature (and therefore the oil properties) was fixed to 40°C while in the experiments, it was varying due to the fact that the cooling/heating system was not able to compensate immediately a change in the temperature.


## 6.2. Results: Gear

Figure7 shows the velocity field in the gear symmetry plane for four different tangential velocities. It can be seen that the area in which the velocity is appreciable is expanding with the increase of the tangential velocity.Figure 7
Contour plot of the velocity field for four different tangential velocities (m/s) in the symmetry plane of the gear: (a) 500 rpm; (b) 2000 rpm; (c) 5000 rpm; (d) 8000 rpm.From Figure8, it is possible to have an idea of the fluxes of lubricant inside the gearbox. An ideal particle of fluid escapes the gear radially and then, after a loop, it comes again in contact with the gear from the axial direction. Due to the rotation of the gear, there are no zones in the gearbox where the lubricant is steady but the complete domain is involved in the lubricant flux.Figure 8
Velocitystreamlines;b=40 mm; d=102.5 mm; T=90°C; n=8000 rpm.Figures9 to 13 show the results in terms of splashing resistant torque PVZ0,Pl versus circumferential velocity vt for the cases with the gear. The continuous line represents the results of the experiments without the losses of the bearings, the broken line the results of the calculations according to Mauz [2], and the dotted line the results of the simulations.Figure 9
Results in terms of resistant torque versus circumferential velocity—effect of the temperature.Mauz proposed some equations derived by some experiments. The model proposed by Mauz has the big advantage that it does not need computational resources but, as shown in the next diagrams, it has also some limits. For example, it is not valid for tangential velocities higher than 60 m/s and it does not take into account some parameters like, for example, the helix angle.Figure9 shows the influence of the temperature on the resistant torque. It can be seen that the influence of the temperature on the results cannot be neglected. The experimental tests were performed with a temperature that is not exactly constant, which is due to the cooling/heating system that is not able to compensate immediately a fluctuation of the temperature of the lubrication bath. Some simulations were performed with the same temperature measured in the experiments for each rotational speed: the lubricant properties have been changed in each simulation according to the measured temperature at the different rotational speeds. It can be seen that, the CFD results are in very good agreement with the measured ones. Some other calculations, instead, were performed with a constant temperature of 90°C. The CFD results for the constant temperature are below the measured ones only for low rotational speeds where the effective temperature of the oil bath in the test rig was lower than 90°C. For higher rotational speeds the CFD results calculated with a constant temperature of 90°C are over the measured ones where the temperature of the oil bath in the test rig was higher than 90°C.The diagram shows also the power losses according to Mauz. The results differ significantly from the measured and simulated values.Figure10 shows the similar results as Figure 8 but for another tip diameter. Also in this case, the results of the simulations performed with the same temperature of the experiments overlap the measured results.Figure 10
Results in terms of resistant torque versus circumferential velocity—effect of the temperature.Figure11 shows the effect of the tip diameter on the power losses. Simulations have been performed for a diameter of 102.5 mm, 98 mm, and 96.5 mm. Experiments have been conducted only for the biggest and the smallest diameter (see Figures 9 and 10).Figure 11
Results in terms of resistant torque versus circumferential velocity—influence of the tip diameter.The effect of the tip diameter on the power losses is extremely high. The results according to Mauz did not correctly describe this phenomenon. The results of the simulations, however, are able to well predict the power losses as confirmed by the experiments.An increase of the tip diameter from 96.5 mm to 98 mm causes an increase of the power losses of about 64%.Similarly, an increase of the tip diameter from 96.5 mm to 102.5 mm causes an increase of the power losses of about 142%.Figure12 shows the effect of the face width on the power losses.Figure 12
Results in terms of resistant torque versus circumferential velocity—Influence of the face width.Figure 13
Results in terms of resistant torque versus circumferential velocity—Influence of the helix angle.The effect of the face width on the power losses is extremely high. Also in this case the results of the simulations are better than the results obtained by Mauz even if in this case also for Mauz an increase of the face width of 100% gives a significant increase of the power losses.The doubling of the tip diameter from 20 mm to 40 mm causes an increase of the power losses of about 60%.Figure13 shows the effect of the helix angle on the power losses. It is possible to appreciate how the losses decrease with the helix angle.In particular, the losses decrease to about 45% changing the helix angle from 0° to 20°. The model proposed by Mauz does not take into account the helix angle.Figure14 shows the effect of the oil viscosity on the power losses, while Figure 15 shows the influence of the oil type (density) on the power losses.Figure 14
Results in terms of resistant torque versus circumferential velocity—Influence of the lubricant viscosity.Figure 15
Results in terms of resistant torque versus circumferential velocity—Influence of the lubricant density.Increasing the density from 819.5 Kg/m3 to 824.5 Kg/m3 leads to an increase of the power losses of about 15%.


## 7. Conclusions

The results of the experiments confirm that the CFD represent a valid method to predict power losses. The error in the predictions for the analyzed cases is lower than 5%.The simulations were performed for two different components, disk and gear. The power losses of the dip lubricated disk appear significantly lower than the power losses of dip lubricated gear.Additional simulations were performed in order to understand the influence of different parameters on the power losses. Power losses increase significantly with a larger tip diameter or face depth and decrease with a larger helix angle. Also the temperature has an influence on the results: the losses decrease with a higher temperature and thus reduced lubricant viscosity and density. But while the influence of the density is significant, the viscosity seems to have no significant effects on the power losses.The trend of the results is in line with those presented by Hӧhn et al. [5] for the case with two gears. Moreover, as expected, the losses evolve with angular velocity power 3.Since CFD proves to be a valid tool to predict the load-independent losses, it is planned to investigate other components of the load-independent power losses, like, for example, the oil squeezing losses.



---

*Source: 616923-2012-08-29.xml*