Thermal Analysis of High-Speed Spur Gears Under Elastohydrodynamic Lubrication: A Combined Numerical and Experimental Investigation

In the realm of high-performance power transmission systems, the operational reliability of spur gears is paramount. Among the various failure modes, scuffing or adhesive wear stands out as a critical concern for high-speed spur gears. This phenomenon is intrinsically linked to excessive flash temperatures at the gear tooth contact interfaces. Elevated surface temperatures not only precipitate lubrication film collapse leading to scuffing but also induce significant thermal distortions. These distortions adversely affect meshing accuracy, load distribution, transmission efficiency, and ultimately, the service life of the gear pair. Consequently, a precise prediction of the steady-state temperature field within high-speed spur gears is not merely an academic exercise but a fundamental requirement for robust design and operational safety. While experimental measurement offers the most direct evidence, it is often constrained by high costs, limited adaptability, and practical difficulties in instrumenting meshing interfaces. Therefore, developing accurate and efficient theoretical and numerical methodologies for thermal analysis remains a central focus in gear research.

Traditional approaches to gear thermal analysis often rely on empirical or constant friction coefficients for calculating frictional heat generation, which fails to capture the transient and complex tribological conditions at the tooth contact. For high-speed, heavily loaded spur gears, the contact is typically governed by elastohydrodynamic lubrication (EHL), where the lubricant pressure-viscosity effect and surface elastic deformation play dominant roles. Ignoring this refined lubrication state can lead to significant inaccuracies in heat source estimation. Furthermore, many studies utilizing the thermal network method, a powerful tool for system-level thermal modeling, often rely on manual derivation of nodal equations and resistance calculations, which can be tedious, error-prone, and limited in nodal resolution.

This work presents a comprehensive methodology for determining the steady-state temperature field of high-speed spur gears, explicitly accounting for transient EHL conditions. The core of our approach involves a synergistic coupling of tribology and heat transfer analysis. We first establish a line-contact EHL model to accurately compute the friction coefficient at each instantaneous contact point along the path of contact. The frictional heat flux is then derived. Subsequently, we employ an advanced implementation of the thermal network method using a co-simulation framework with Patran (pre/post-processor) and Sinda (solver). This framework automates the generation of a dense nodal network, calculates thermal resistances/capacities, and efficiently solves for the steady-state temperature distribution. The feasibility and accuracy of this combined numerical approach are validated through temperature measurements obtained from a dedicated gear test rig.

Theoretical Foundation: Meshing Parameters and Contact Mechanics

The kinematic and load conditions vary continuously as a pair of spur gear teeth mesh. To analyze the thermal state, we must first quantify these varying parameters along the line of action. The coordinate system is defined along the line of action (N1N2), with the pitch point P as the origin. The coordinate \(s\) defines the position of any contact point.

The boundaries of the contact path are defined by the start of active profile (SAP) and end of active profile (EAP). Their distances from the pitch point are:

$$ PB_1 = \frac{\sqrt{d_{a1}^2 – d_{b1}^2}}{4} – \frac{d’_1 \sin \alpha’}{2} $$
$$ PB_2 = \frac{\sqrt{d_{a2}^2 – d_{b2}^2}}{4} – \frac{d’_2 \sin \alpha’}{2} $$

where \(d_{a}\), \(d_{b}\), and \(d’\) are the addendum circle, base circle, and operating pitch circle diameters, respectively, and \(\alpha’\) is the operating pressure angle. The range for \(s\) is \([x_1, x_2]\), where \(x_1 = -PB_2\) and \(x_2 = PB_1\).

Sliding and Entrainment Velocity: The tangential velocities of the pinion and gear at a contact point are:

$$ U_1(s) = \frac{2\pi n_1}{60\,000} \left( \frac{1}{2}d’_1 \sin\alpha’ + s \right) $$
$$ U_2(s) = \frac{2\pi n_2}{60\,000} \left( \frac{1}{2}d’_2 \sin\alpha’ – s \right) $$

The entrainment velocity \(U\) (average surface speed) and the sliding velocity \(U_s\) are critical for EHL and heat generation calculations:

$$ U(s) = \frac{U_1(s) + U_2(s)}{2} $$
$$ U_s(s) = |U_1(s) – U_2(s)| $$

Contact Geometry and Load: The radii of curvature for the pinion and gear tooth profiles at the contact point are:

$$ R_1(s) = \frac{1}{2}d_1 \sin\alpha + s $$
$$ R_2(s) = \frac{1}{2}d_2 \sin\alpha – s $$

The equivalent radius of curvature \(R\) is given by:

$$ R(s) = \frac{R_1(s) R_2(s)}{R_1(s) + R_2(s)} $$

The normal load distribution along the line of action is not uniform due to load sharing between gear teeth. For a standard non-modified spur gear, the load-sharing factor \(K(s)\) over the zones of single and double tooth contact is:

$$
K(s) =
\begin{cases}
\frac{1}{3}\left(1 + \frac{s – x_1}{x_0}\right) & x_1 \le s \le x_1 + x_0 \\
1 & x_1 + x_0 < s < x_2 – x_0 \\
\frac{1}{3}\left(1 – \frac{s – x_2}{x_0}\right) & x_2 – x_0 \le s \le x_2
\end{cases}
$$

where \(x_0 = \pi m (\epsilon – 1) \cos \alpha\), \(m\) is the module, and \(\epsilon\) is the contact ratio. The normal load \(F_n\) per unit face width on the pinion tooth is:

$$ F_n(s) = \frac{K(s) T_2 z_1}{r_c(s) \cos\alpha\, z_2} $$

Here, \(T_2\) is the gear torque, \(z_1\) and \(z_2\) are tooth numbers, and \(r_c(s)\) is the radial distance on the pinion to the contact point:

$$ r_c(s) = \sqrt{\left( \frac{1}{2}d_1 + s \sin\alpha \right)^2 + (s \cos\alpha)^2 } $$

Elastohydrodynamic Lubrication Model for Friction Coefficient

The meshing of spur gears can be accurately modeled as a transient line-contact EHL problem. At any instantaneous position \(s\), the contact is equivalent to that between two equivalent cylinders, or one cylinder and a plane, with radius \(R(s)\), entrainment speed \(U(s)\), and load \(w(s)=F_n(s)\). The isothermal line-contact EHL governing equations are presented below in dimensional form.

Reynolds Equation: Governs the hydrodynamic pressure generation within the lubricant film.

$$ \frac{d}{dx}\left( \frac{\rho h^3}{\eta} \frac{dp}{dx} \right) = 12U \frac{d(\rho h)}{dx} $$

Boundary conditions: \(p(x_{inlet})=0\) and \(p(x_{exit}) = dp/dx|_{x_{exit}} = 0\).

Film Thickness Equation: Relates the gap between the surfaces to the geometry and elastic deformation.

$$ h(x) = h_0 + \frac{x^2}{2R} + v(x) $$

where \(h_0\) is the central film thickness offset and \(v(x)\) is the elastic deformation.

Elastic Deformation Equation (Boussinesq):

$$ v(x) = -\frac{2}{\pi E’} \int_{x_{inlet}}^{x_{exit}} p(s) \ln(s-x)^2 ds + C $$

where \(E’\) is the effective elastic modulus: \( \frac{1}{E’} = \frac{1}{2}\left( \frac{1-\nu_1^2}{E_1} + \frac{1-\nu_2^2}{E_2} \right) \).

Viscosity-Pressure Equation (Roelands): Models the dramatic increase in lubricant viscosity with pressure.

$$ \eta(p) = \eta_0 \exp\left\{ (\ln\eta_0 + 9.67) \left[ \left(1 + \frac{p}{p_0}\right)^{0.68} – 1 \right] \right\} $$

Density-Pressure Equation (Dowson-Higginson):

$$ \rho(p) = \rho_0 \left( \frac{1 + 0.6p}{1 + 1.7p} \right) $$

Force Balance Equation: Ensures the integrated pressure supports the applied load.

$$ \int_{x_{inlet}}^{x_{exit}} p(x) dx = w $$

Friction Calculation: The shear stress \(\tau\) in a Newtonian lubricant film is:

$$ \tau = \eta \frac{\partial u}{\partial z} = \pm \frac{\partial p}{\partial x}\frac{h}{2} + \eta \frac{U_2 – U_1}{h} $$

The friction forces on the two surfaces per unit width are found by integrating the shear stress at the surfaces (\(z=0\) and \(z=h\)) over the contact domain:

$$ F_{1,2} = \int_{x_{inlet}}^{x_{exit}} \tau|_{z=0,h} \, dx $$

The instantaneous coefficient of friction at contact position \(s\) is then:

$$ f(s) = \frac{|F_1 + F_2|}{w} $$

This system of equations is solved numerically for each contact position \(s\) along the gear mesh cycle using methods like the finite difference method with progressive mesh refinement and the multigrid technique. The solution yields the pressure distribution \(p(x)\), film thickness \(h(x)\), and crucially, the friction coefficient \(f(s)\).

Summary of Key EHL Equations:

Equation Purpose Mathematical Form
Reynolds Pressure Generation $\frac{d}{dx}\left( \frac{\rho h^3}{\eta} \frac{dp}{dx} \right) = 12U \frac{d(\rho h)}{dx}$
Film Thickness Geometry & Deformation $h(x) = h_0 + \frac{x^2}{2R} + v(x)$
Elasticity Surface Deformation $v(x) = -\frac{2}{\pi E’} \int p(s) \ln(s-x)^2 ds + C$
Viscosity-Pressure Lubricant Property $\eta(p) = \eta_0 \exp\left\{ A \left[ (1+p/p_0)^{0.68} -1 \right] \right\}$
Density-Pressure Lubricant Property $\rho(p) = \rho_0 \left( \frac{1 + 0.6p}{1 + 1.7p} \right)$
Force Balance Load Equilibrium $\int p(x) dx = w$

Determination of Frictional Heat Flux

The frictional heat generated at the contact interface of meshing spur gears is the primary internal heat source. Its magnitude depends on the friction coefficient, normal pressure, and sliding velocity. The heat partition between the pinion and gear teeth is unequal and depends on their thermophysical properties and surface speeds.

The instantaneous frictional heat flux density entering the pinion surface at position \(s\) is:

$$ q_1(s) = \varepsilon \kappa(s) f(s) p_h(s) U_s(s) $$

where:
• \(\varepsilon\) is the heat conversion efficiency (typically 0.90–0.95).
• \(\kappa(s)\) is the heat partition coefficient for the pinion.
• \(p_h(s)\) is the maximum Hertzian contact pressure.
• \(f(s)\) is the EHL friction coefficient from the previous section.
• \(U_s(s)\) is the sliding velocity.

The heat partition coefficient \(\kappa\) is derived from the assumption of continuous temperature at the interface and is given by the Blok flash temperature theory:

$$ \kappa(s) = \frac{ \sqrt{\lambda_1 \rho_1 c_1 U_1(s)} }{ \sqrt{\lambda_1 \rho_1 c_1 U_1(s)} + \sqrt{\lambda_2 \rho_2 c_2 U_2(s)} } $$

where \(\lambda\), \(\rho\), and \(c\) are the thermal conductivity, density, and specific heat capacity, respectively.

For steady-state thermal analysis, the transient moving heat source must be averaged over one mesh cycle for a given point on the tooth flank. The time-averaged heat flux density applied to a specific nodal area on the pinion tooth surface is:

$$ q_{1,avg} = \frac{t_b}{T} q_1 = \left( \frac{b}{U_1(s)} \cdot \frac{n_1}{60} \right) q_1(s) $$

Here, \(t_b\) is the time the contact point spends on the tooth flank segment corresponding to the node, \(T\) is the rotational period of the pinion, \(b\) is the face width, and \(n_1\) is the rotational speed (rpm). This averaging transforms the transient local heat source into a steady, distributed heat input for the finite thermal network model.

The Thermal Network Method and Co-Simulation Framework

The Thermal Network Method (TNM), or lumped capacitance method, is an efficient approach for solving complex steady-state and transient heat transfer problems. It discretizes a physical system into a network of nodes, each representing a finite volume with a uniform temperature. These nodes are connected by thermal resistances (conduction, convection, radiation) and possess thermal capacitances.

Governing Principle: For any node \(i\) in the network, the principle of energy conservation (analogous to Kirchhoff’s current law) applies. For steady-state analysis, the sum of all heat flows into the node is zero:

$$ \sum_{j} \frac{T_j – T_i}{R_{i,j}} + Q_i = 0 $$

where:
• \(T_i\), \(T_j\) are the temperatures of nodes \(i\) and \(j\).
• \(R_{i,j}\) is the thermal resistance between nodes \(i\) and \(j\).
• \(Q_i\) is the internal heat source (or sink) at node \(i\) (e.g., frictional heat flux \(q_{1,avg} \times\) area).

Thermal resistances are calculated based on geometry and mode:
• Conduction: \(R_{cond} = L / (\lambda A)\)
• Convection: \(R_{conv} = 1 / (h_c A)\)
• Radiation: Linearized as \(R_{rad} = 1 / (h_r A)\), where \(h_r\) depends on temperature.

Patran-Sinda Co-Simulation Framework: Traditional TNM implementation often struggles with manual node generation, resistance calculation, and equation solving for complex geometries like spur gear teeth. To overcome this, we employ a co-simulation strategy using MSC Patran and Thermal Desktop/Sinda.
1. Geometry and Meshing (Patran): A detailed 3D model of a single spur gear tooth is created and imported into Patran. Instead of a finite element mesh for thermal analysis, a “thermal mesh” of nodes is generated on the geometry surfaces and within the volume. The density of this node distribution can be controlled for accuracy.
2. Property and Boundary Condition Assignment (Patran): Material properties (λ, ρ, c) are assigned. Boundary conditions are applied:
Heat Flux: The calculated \(q_{1,avg}\) is mapped onto the relevant tooth flank nodes.
Convection: Different convective heat transfer coefficients (\(h_c\)) are assigned to tooth flanks (affected by oil jet/spray), tooth tips, side faces (free convection/oil splash), and the root/fillet regions. Correlations for gear convective cooling are used.
Conduction Links: The software automatically calculates conductive resistances between adjacent volume and surface nodes based on material and distance.
Fixed Temperature: The inner bore surface of the gear may be set to a measured or estimated bulk temperature.
3. Model Translation and Solving (Sinda): Patran, configured to use Sinda as the solver, translates the nodal network, all thermal resistances, capacitances (for transient), and boundary conditions into a Sinda input file. Sinda then constructs and solves the large system of linear (steady-state) or differential (transient) equations from the thermal network.
4. Post-Processing (Patran): The temperature results from Sinda are read back into Patran. The temperatures at the discrete nodes are interpolated to produce detailed temperature contour plots on the 3D gear tooth model.

This automated workflow ensures high nodal resolution, eliminates manual calculation errors, and significantly speeds up the analysis process compared to classical TNM.

Summary of Boundary Condition Types:

Surface Type Primary Thermal Process Typical Coefficient / Value
Contact Flank Prescribed Heat Flux, Convection $q_{1,avg}(s)$, $h_{c,flank}$ (Oil Jet)
Tooth Tip Convection $h_{c,tip}$ (Oil Churning/Splash)
Side Faces Convection $h_{c,side}$ (Free/Oil Mist Convection)
Root/Fillet Area Convection $h_{c,root}$ (Lower due to geometry)
Gear Bore Prescribed Temperature / Convection $T_{bulk}$ or $h_{c,bore}$
Internal Volume Nodes Conduction Links via $R_{cond}$

Case Study: Analysis of a High-Speed Spur Gear Test Pinion

We apply the developed methodology to analyze a pinion from a high-speed gear test rig. The primary geometrical and operational parameters are listed below.

Parameter Pinion Gear Unit
Number of Teeth, \(z\) 16 24
Module, \(m\) 4.5 mm
Face Width, \(b\) 20 mm
Pressure Angle, \(\alpha\) 20 deg
Addendum Modification Coeff. +0.85 -0.50
Rotational Speed, \(n\) 1500 1000 rpm
Transmitted Torque, \(T\) 70 Nm
Material (Steel) Elastic Modulus \(E\) = 210 GPa, Density \(\rho\) = 7800 kg/m³, Specific Heat \(c\) = 500 J/(kg·K), Conductivity \(\lambda\) = 40 W/(m·K)
Lubricant ISO VG 32, Dynamic Viscosity \(\eta_0\) = 0.014 Pa·s at 40°C

Step 1: Calculation of Meshing Parameters and Friction Coefficient
Using the formulas in Section 2, key parameters along the path of contact are computed. The results for normal load per unit width \(w(s)=F_n(s)\) and sliding velocity \(U_s(s)\) are plotted. The load shows a characteristic pattern with higher, constant load in the single-pair contact zone and lower, linearly varying load in the double-pair contact zones. The sliding velocity is zero at the pitch point (\(s=0\)) and increases towards the tips and roots.

The EHL model is solved numerically for multiple discrete points \(s\) across the contact path. The calculated friction coefficient \(f(s)\) exhibits a complex profile: it is generally higher near the start and end of contact (low entrainment velocity, high sliding), dips in the mid-regions, and approaches zero at the pitch point where sliding ceases. The corresponding instantaneous frictional heat flux \(q_1(s)\) follows a similar trend but is modulated by the product \(f(s) \cdot U_s(s) \cdot p_h(s)\).

Step 2: Thermal Network Model Setup and Solution
A single-tooth model of the pinion is built. The Patran-Sinda co-simulation framework is employed. The tooth volume and surfaces are populated with a dense distribution of thermal nodes. The time-averaged heat flux \(q_{1,avg}\) is applied to the active flank nodes. Convective boundaries with coefficients representative of forced oil jet cooling on the flanks and free/oil mist convection on other exposed surfaces are applied. The temperature at the gear bore is set as a fixed boundary condition, estimated from test rig measurements. The Sinda solver computes the steady-state temperature for every node in the network.

Step 3: Numerical Results
The steady-state temperature field of the pinion tooth is recovered and visualized. The results indicate a temperature range of approximately 73°C to 77°C on the tooth body. The maximum temperature region is located on the active flank, slightly offset from the pitch line towards the addendum, correlating with the zone of significant heat flux. The temperature distribution is elliptical across the face width, with the center region being hotter than the sides due to better convective cooling at the side faces. A clear temperature gradient exists from the hot flank surface into the cooler tooth core.

Step 4: Experimental Validation
Temperature measurements were conducted on the gear test rig under identical operating conditions. An infrared thermal imaging camera (Fluke Ti32), configured with appropriate emissivity (0.7) and transmission settings for the oil mist environment, was used to capture the steady-state temperature distribution on the accessible side face of the rotating gear. The experimental results showed a maximum tooth flank temperature of approximately 78.3°C.

Comparison and Discussion: The numerical simulation predicted a maximum temperature of 76.7°C. The discrepancy of about 1.6°C is within the typical uncertainty bounds of IR measurements and modeling assumptions (e.g., precise convective coefficients, material properties). More importantly, the spatial pattern of temperature distribution—the location of the hot spot and the trend of temperature variation along the tooth profile—showed excellent qualitative agreement between the simulation and the experiment. This validates the overall methodology, particularly the use of EHL-derived friction for heat source calculation and the efficacy of the Patran-Sinda co-simulation framework for solving the thermal network of spur gears.

Conclusion

This study successfully developed and demonstrated an integrated methodology for analyzing the steady-state temperature field of high-speed spur gears. The key conclusions are:

  1. Importance of EHL-Based Friction: Utilizing a line-contact elastohydrodynamic lubrication model to determine the transient friction coefficient along the gear meshing path provides a physically realistic estimation of the frictional heat source. This approach is superior to using constant empirical values, as it captures the significant variations caused by changing kinematics, load, and lubrication regime.
  2. Advanced Thermal Network Solution: The proposed co-simulation framework using Patran for model preparation and Sinda as the thermal network solver overcomes the major limitations of traditional manual thermal network methods. It enables the automatic generation of a high-resolution nodal network, accurate calculation of thermal resistances, and efficient solution for complex 3D geometries like spur gear teeth.
  3. Validation and Findings: The numerical predictions were in good agreement with experimental infrared thermography measurements. Both methods confirmed that for high-speed spur gears under the given conditions, the maximum temperature occurs on the tooth flank, near the center of the face width and slightly towards the addendum from the pitch line. The temperature distribution exhibits an elliptical pattern across the face width and a gradient from the surface to the core.
  4. Practical Utility: The presented methodology offers a practical and accurate tool for gear designers and engineers to predict operational temperatures during the design phase. This capability is crucial for assessing the risk of thermal failure modes like scuffing, optimizing lubrication and cooling systems, and evaluating the impact of thermal distortions on gear mesh performance, thereby enhancing the reliability of spur gear transmissions.

Future work may focus on extending the analysis to include transient thermal effects during startup or load changes, incorporating non-Newtonian lubricant behavior into the EHL model for even more accurate friction prediction, and applying the method to more complex gear geometries such as helical or hypoid gears.

Scroll to Top