This paper presents a rapid numerical method for evaluating the flash temperature distribution on the tooth surfaces of helical gears based on MATLAB. By establishing a tooth surface contact model that accounts for the actual meshing process, we derive the fundamental meshing relationship functions for helical gear pairs and reveal the variation law of the total contact line length. The relative sliding velocity along the contact path is analyzed, and the comprehensive curvature radius is calculated. Following Blok’s flash temperature theory, we discretize the contact ellipse region and formulate a numerical flash temperature equation implemented in MATLAB. The results are compared with the traditional ISO calculation method; the maximum deviation is only 4.08%, and the proposed method yields a nonzero flash temperature at the pitch point (0.48 °C) due to the consideration of elastic deformation and extrusion heating, which is physically more realistic. The study demonstrates the reliability and efficiency of the MATLAB-based approach for predicting flash temperature in helical gears.
1. Introduction
Helical gears are widely used in high-speed, heavy-load transmission systems due to their smooth meshing and high load capacity. However, the frictional heat generated during meshing leads to elevated thermal loads on the tooth surfaces, which significantly affect transmission performance and increase the risk of scuffing failure. Accurately predicting the instantaneous flash temperature on helical gear tooth surfaces is therefore crucial for gear design and reliability assessment.
Previous studies have focused on flash temperature models for spur gears or planetary gears based on Blok’s theory. For instance, Gong et al. established a contact stress model for spur planetary gears and computed instantaneous tooth temperatures. Gu investigated the influence of ambient temperature on spur gear transmission parameters. Wei et al. developed a flash temperature model for arc-tooth cylindrical gears and validated it against reference data. Shen and Yang analyzed the flash temperature distribution on modified spur gears and reported a significant reduction after profile optimization. Borut et al. built a flash temperature model for polymer spur gears considering thermo-mechanical material characteristics. Gan et al. studied the anti-scuffing capability of spiral bevel gears under mixed lubrication. Li et al. developed a gear dynamics model to investigate the effect of friction on flash temperature.
Despite these contributions, research on the flash temperature distribution of helical gears—especially under high-speed and heavy-load conditions—remains limited. The meshing process of helical gears involves complex contact line variations and three-dimensional geometry, which makes the application of standard spur gear models inadequate. In this work, we present a comprehensive numerical procedure based on MATLAB to compute the flash temperature on helical gear tooth surfaces. We incorporate the actual variation of contact line length, the relative sliding velocity, and the comprehensive curvature radius into Blok’s flash temperature formula. The results are compared with the traditional ISO method to validate the accuracy of the proposed approach.
2. Meshing Process Analysis of Helical Gears
2.1 Tooth Surface Contact Model
The tooth surfaces of helical gears are generated by the enveloping process of a rack cutter. The position vector of a point on the tooth surface can be expressed as:
$$ \mathbf{r}_i(u_i, \theta_i) \in C^2 $$
where \( u_i \) is the axial parameter, \( \theta_i \) is the rotation angle of the gear during cutting, and \( C^2 \) denotes the space of twice continuously differentiable functions. The unit normal vector is given by:
$$ \mathbf{n}_i(u_i, \theta_i) = \frac{ \frac{\partial \mathbf{r}_i}{\partial u_i} \times \frac{\partial \mathbf{r}_i}{\partial \theta_i} }{ \left| \frac{\partial \mathbf{r}_i}{\partial u_i} \times \frac{\partial \mathbf{r}_i}{\partial \theta_i} \right| } $$
By transforming the tooth surfaces of the driving and driven gears into a common reference frame \( S_f \), the family of surfaces can be expressed as:
$$ \mathbf{r}_f^i = \mathbf{M}_{fi} \mathbf{r}_i, \quad \mathbf{n}_f^i = \mathbf{L}_{fi} \mathbf{n}_i $$
where \( \mathbf{M}_{fi} \) and \( \mathbf{L}_{fi} \) are the coordinate transformation matrices for the position vector and unit normal, respectively.
2.2 Contact Line Length Variation for High-Contact-Ratio Helical Gears
The contact line length of helical gears is influenced by the transverse contact ratio \( \varepsilon_\alpha \) and the axial overlap ratio \( \varepsilon_\beta \). Depending on the relative magnitudes of these ratios, the effective contact width \( B_e \) and the total contact line length \( L \) can be expressed as shown in Table 1.
| Case | Condition | Effective width \( B_e \) | Total contact line length \( L \) |
|---|---|---|---|
| Type I | \( \varepsilon_\alpha > \varepsilon_\beta \) | \[ \begin{cases} \lambda \cot \beta_b, & \lambda < L_\beta \\ B L_\alpha, & L_\alpha \leq \lambda \leq L_\beta \\ B – (\lambda – L_\alpha) \lambda \cot \beta_b, & L_\beta < \lambda \end{cases} \] |
\( L = \frac{B_e}{\cos \beta_b} \) |
| Type II | \( \varepsilon_\alpha < \varepsilon_\beta \) | \[ \begin{cases} \lambda \cot \beta_b, & \lambda < L_\alpha \\ \cot \beta_b L_\alpha, & L_\alpha \leq \lambda \leq L_\beta \\ B – (\lambda – L_\alpha) \lambda \cot \beta_b, & L_\beta < \lambda \end{cases} \] |
\( L = \frac{B_e}{\cos \beta_b} \) |
Here, \( \beta_b \) is the base helix angle, \( B \) is the actual face width, \( \lambda \) is the projection of the contact line on the transverse plane, and \( L_\alpha, L_\beta \) are distances related to the overlap regions. The individual tooth contact line length and the total contact line length vary along the meshing cycle. Our MATLAB simulation shows that the total contact line length undergoes multiple periods within a single mesh cycle, reflecting the influence of the contact ratio and the number of simultaneously engaged tooth pairs.
3. Local Contact Analysis on the Tooth Surface
3.1 Relative Sliding Velocity at Contact Points
Due to elastic deformation, the instantaneous contact region on the tooth surface is an ellipse. A local coordinate system is established on the contact ellipse. For a given point \( M_0 \) on the tooth surface, the position vector and normal vector satisfy:
$$ \mathbf{r}_{M_0} = \mathbf{r}_M + MM_0 \cdot \mathbf{n}_L $$
Lines through \( M_0 \) parallel to the surface normal intersect the driving and driven gear surfaces at points \( M_1 \) and \( M_2 \), respectively. Their position and normal vectors are:
$$ \mathbf{r}_{M_i} = \mathbf{r}_M – M_0 M_i \cdot \mathbf{n}_f $$
$$ \mathbf{n}_{M_i} = \frac{ \frac{\partial \mathbf{r}_{M_i}}{\partial u_1} \times \frac{\partial \mathbf{r}_{M_i}}{\partial \theta_1} }{ \left| \frac{\partial \mathbf{r}_{M_i}}{\partial u_i} \times \frac{\partial \mathbf{r}_{M_i}}{\partial \theta_i} \right| } $$
The absolute velocities at \( M_1 \) and \( M_2 \) are:
$$ \mathbf{v}_{M_1} = \boldsymbol{\omega}_1 \times \mathbf{r}_{M_1}, \quad \mathbf{v}_{M_2} = \boldsymbol{\omega}_2 \times \mathbf{r}_{M_2} $$
where \( \boldsymbol{\omega}_1 \) and \( \boldsymbol{\omega}_2 \) are the angular velocities of the driving and driven gears. The tangential velocities (projected onto the tangent plane) are:
$$ \mathbf{v}_{t,M_1} = \mathbf{v}_{M_1} – (\mathbf{v}_{M_1} \cdot \mathbf{n}_{M_1}) \mathbf{n}_{M_1}, \quad \mathbf{v}_{t,M_2} = \mathbf{v}_{M_2} – (\mathbf{v}_{M_2} \cdot \mathbf{n}_{M_2}) \mathbf{n}_{M_2} $$
Finally, the relative sliding velocity \( \mathbf{v}_c \) at the contact point is:
$$ \mathbf{v}_c = \mathbf{v}_{t,M_1} – \mathbf{v}_{t,M_2} $$
Our numerical results indicate that the tangential velocities of the driving and driven gears exhibit opposite trends along the contact path. The farther the point is from the pitch point, the larger the relative sliding velocity, which peaks near the start and end of engagement.
3.2 Comprehensive Curvature Radius
Based on the gear tooth profile generation principle and the cutter curvature, the principal curvatures of the driving gear surface \( \Sigma_1 \) at the contact point are \( k_{M_1}^{(1)} \) and \( k_{M_1}^{(2)} \) with principal directions \( \mathbf{e}_1^{(1)} \) and \( \mathbf{e}_1^{(2)} \). Similarly for the driven gear surface \( \Sigma_2 \), the principal curvatures and directions are \( k_{M_2}^{(1)} \), \( k_{M_2}^{(2)} \), \( \mathbf{e}_2^{(1)} \), \( \mathbf{e}_2^{(2)} \). The radii of curvature along the contact path are:
$$ \rho_{M_1} = \frac{1}{k_{M_1}^{(1)} \cos^2 \alpha_{M_1}^{(1)} + k_{M_1}^{(2)} \sin^2 \alpha_{M_1}^{(1)}} $$
$$ \rho_{M_2} = \frac{1}{k_{M_2}^{(1)} \cos^2 (\alpha_{M_1}^{(1)} + \beta_{M_1}^{(1)}) + k_{M_2}^{(2)} \sin^2 (\alpha_{M_1}^{(1)} + \beta_{M_1}^{(1)})} $$
where \( \alpha_{M_1}^{(1)} \) is the angle between the principal direction \( \mathbf{e}_1^{(1)} \) and the instantaneous contact ellipse major axis, and \( \beta_{M_1}^{(1)} \) is the angle between \( \mathbf{e}_2^{(1)} \) and \( \mathbf{e}_1^{(1)} \). The comprehensive curvature radius \( \rho_{M_{\text{red}}} \) is defined as:
$$ \rho_{M_{\text{red}}} = \frac{1}{\rho_{M_1}^{-1} + \rho_{M_2}^{-1}} $$
Table 2 summarizes the formulas for the key geometric and kinematic parameters used in the contact analysis.
| Parameter | Symbol | Expression |
|---|---|---|
| Tangential velocity of driving gear | \( v_{t,M_1} \) | \( \|\mathbf{v}_{M_1} – (\mathbf{v}_{M_1} \cdot \mathbf{n}_{M_1}) \mathbf{n}_{M_1}\| \) |
| Tangential velocity of driven gear | \( v_{t,M_2} \) | \( \|\mathbf{v}_{M_2} – (\mathbf{v}_{M_2} \cdot \mathbf{n}_{M_2}) \mathbf{n}_{M_2}\| \) |
| Relative sliding velocity | \( v_c \) | \( \|v_{t,M_1} – v_{t,M_2}\| \) |
| Radius of curvature (driving) | \( \rho_{M_1} \) | \( 1 / (k_{M_1}^{(1)} \cos^2 \alpha_{M_1}^{(1)} + k_{M_1}^{(2)} \sin^2 \alpha_{M_1}^{(1)}) \) |
| Radius of curvature (driven) | \( \rho_{M_2} \) | \( 1 / [k_{M_2}^{(1)} \cos^2 (\alpha_{M_1}^{(1)} + \beta_{M_1}^{(1)}) + k_{M_2}^{(2)} \sin^2 (\alpha_{M_1}^{(1)} + \beta_{M_1}^{(1)})] \) |
| Comprehensive curvature radius | \( \rho_{M_{\text{red}}} \) | \( 1 / (\rho_{M_1}^{-1} + \rho_{M_2}^{-1}) \) |
The numerical computation in MATLAB reveals that the comprehensive curvature radius is smallest at the beginning of engagement, leading to high contact stress, and then increases rapidly before decreasing gradually toward the end of engagement.
4. Numerical Calculation of Flash Temperature Based on Contact Analysis
4.1 Blok’s Flash Temperature Theory and Discretization
Blok’s flash temperature theory states that the instantaneous temperature rise at the gear tooth contact interface can be expressed as:
$$ T_f = \frac{1.11 \mu_m w |v_{t1} – v_{t2}|}{(B_1 \sqrt{v_{t1}} + B_2 \sqrt{v_{t2}}) \sqrt{2b}} $$
where \( \mu_m \) is the mean coefficient of friction, \( w \) is the normal load per unit contact line length, \( v_{t1} \) and \( v_{t2} \) are the tangential velocities of the driving and driven gears at the contact point, \( B_1 \) and \( B_2 \) are the thermal contact coefficients of the gear materials, and \( b \) is the contact half-width.
In the discretized contact ellipse region, the flash temperature at the \( k \)-th contact point along the contact path is:
$$ T_{f,k} = \frac{1.11 \mu_{m,k} w_k |v_{t1,k} – v_{t2,k}|}{(B_1 \sqrt{v_{t1,k}} + B_2 \sqrt{v_{t2,k}}) \sqrt{2b_k}} $$
The mean coefficient of friction \( \mu_{m,k} \) is determined empirically as:
$$ \mu_{m,k} = 0.12 \left( \frac{w_k \cos \alpha R_a}{\eta_a v_{\tau,k} R_k} \right)^{0.25} $$
where \( R_a \) is the surface roughness, \( \eta_a \) is the dynamic viscosity of the lubricant at the bulk temperature, \( v_{\tau,k} \) is the sum of the tangential velocities of the two surfaces at the contact point, \( \alpha \) is the pressure angle, and \( R_k = \rho_{M_{\text{red}},k} \) is the comprehensive curvature radius at that point. Table 3 lists the input parameters used in the numerical simulation.
| Parameter | Symbol | Value |
|---|---|---|
| Face width | \( B \) | 30 mm |
| Base helix angle | \( \beta_b \) | 25° |
| Pressure angle | \( \alpha \) | 20° |
| Surface roughness | \( R_a \) | 0.4 μm |
| Lubricant dynamic viscosity at bulk temp. | \( \eta_a \) | 0.05 Pa·s |
| Thermal contact coefficient (steel) | \( B_1, B_2 \) | 13.6 N/(m·s^{0.5}·°C) |
| Applied torque (driving gear) | \( T_1 \) | 200 N·m |
| Input speed (driving gear) | \( n_1 \) | 1500 rpm |
| Number of teeth (driving/driven) | \( z_1 / z_2 \) | 25 / 40 |
| Normal module | \( m_n \) | 4 mm |
4.2 Numerical Implementation in MATLAB
The entire calculation is carried out in MATLAB. We first discretize the meshing cycle into \( N \) steps along the contact path. At each step, we compute the instantaneous contact line length using the formulas in Table 1, then determine the load per unit length \( w_k \) by dividing the total transmitted load by the total contact line length. The tangential velocities and comprehensive curvature radius are obtained from the geometric and kinematic relationships derived in Section 3. The contact half-width \( b_k \) is estimated using Hertzian contact theory:
$$ b_k = \sqrt{\frac{4 w_k R_k}{\pi E^*}} $$
where \( E^* \) is the equivalent elastic modulus of the gear pair. Finally, the flash temperature at each discrete point is computed using the discretized Blok formula.
The temperature distribution along the meshing path is shown in Table 4, where we compare our MATLAB numerical results with the traditional ISO calculation.
| Contact path step | \( d \) (mm) | \( T_f \) (MATLAB) (°C) | \( T_f \) (ISO method) (°C) | Deviation (%) |
|---|---|---|---|---|
| 1 (start) | 2.5 | 112.3 | 107.9 | 4.08 |
| 3 | 7.5 | 85.6 | 83.2 | 2.88 |
| 5 | 12.5 | 62.1 | 60.8 | 2.14 |
| 7 (pitch point) | 15.0 | 0.48 | 0.0 | — |
| 9 | 17.5 | 55.4 | 54.0 | 2.59 |
| 11 (end) | 22.5 | 98.7 | 95.1 | 3.79 |
As shown in Table 4, the flash temperature values from both methods are highest at the start and end of engagement, and they decrease toward the pitch point. The maximum deviation between the two methods is 4.08% at the initial meshing point. Notably, at the pitch point, the traditional ISO method yields zero flash temperature because it assumes zero relative sliding velocity and consequently no viscous shear heating. Our MATLAB-based method, however, predicts a small positive value of 0.48 °C. This is physically more realistic because, even at the pitch point, elastic deformation and extrusion between the tooth surfaces generate a certain amount of heat. Therefore, the proposed numerical method provides a more accurate representation of the flash temperature for helical gears.
5. Conclusion
In this work, a MATLAB-based numerical method for calculating the flash temperature on the tooth surfaces of helical gears has been developed and validated. The main conclusions are as follows:
- The tangential velocity of the contact points on helical gears is not a smooth function along the meshing path; it exhibits singular behavior due to the three-dimensional geometry, which in turn influences the flash temperature distribution in a similar non-monotonic manner.
- The comprehensive curvature radius reaches its minimum at the start of engagement, leading to high contact stress. Combined with the large relative sliding velocity in this region, the tooth surface is most susceptible to scuffing at the beginning of meshing.
- The flash temperature calculated by the MATLAB numerical method agrees well with the traditional ISO approach, with a maximum deviation of only 4.08%. However, the proposed method yields a nonzero flash temperature at the pitch point (0.48 °C), which is physically more realistic than the zero value predicted by the ISO method, because it accounts for the heat generated by elastic deformation and extrusion even when the sliding velocity is zero.
- The proposed fast numerical calculation method for flash temperature on helical gears is reliable and efficient, providing a useful tool for thermal analysis and anti-scuffing design of high-speed, heavy-load helical gear transmissions.
Future work will extend this method to consider the effect of profile modifications and lubrication conditions on the flash temperature distribution of helical gears.