Gear transmission, as one of the most widely used mechanical transmission methods, is extensively applied in various mechanical fields such as aerospace, aviation, and marine engineering. Its performance plays a crucial role in the reliability of the entire transmission system. At the same time, the reliability of gear meshing is closely related to the safety of personnel. Therefore, conducting a reasonably comprehensive analysis of gear meshing contact stress, thereby guiding the design and development of components, is of great significance.
Many domestic and international experts and scholars have conducted extensive research on the contact stress during the meshing process of straight spur gears. However, most of this research is based on deterministic analysis methods, which ignore the influence of uncertainties in related variables on contact stress. The results of such studies lack universality and cannot effectively guide the design and development of gear components. Therefore, to more objectively and accurately describe the distribution of contact stress during the gear meshing process and improve the reliability of gear transmission, the randomness of factors influencing meshing contact stress should be considered. This necessitates transforming the deterministic analysis of gear meshing contact stress into reliability analysis.
Probabilistic analysis methods have been widely applied in many fields, but they have not yet been utilized for the reliability analysis of meshing gear contact stress. Conducting a probabilistic analysis of meshing gear contact stress allows us to not only obtain the probabilistic distribution characteristics of the contact stress based on the distribution features of random variables but also determine the variation of random variables based on these distribution characteristics. This facilitates the design optimization of meshing gears and helps improve the reliability of the gear transmission system. Based on finite element analysis, this paper employs a hybrid response surface method that combines the Monte Carlo Method (MCM) and the Response Surface Method (RSM). Considering the influence of random variables on contact stress within the nonlinear time domain during the meshing phase of straight spur gears, a nonlinear reliability sensitivity analysis of straight spur gear contact stress is conducted.
Hybrid Response Surface Method
The hybrid response surface method integrates the Monte Carlo method with the response surface method. It rationally selects test points and iteration strategies, conducts a series of deterministic experiments via the Monte Carlo method to obtain a series of random variable and output response sample points, and fits a polynomial response surface function. This function converges to the true implicit limit state function within the failure probability domain. To describe the mathematical model of the hybrid response surface method, let \( Y \) and \( \mathbf{X} = [X_1, X_2, \ldots, X_i, \ldots, X_r]^T \) (where \( X_i \sim N(\mu_i, \sigma_i^2) \)) represent the structural output response and random variables, respectively. Using the Monte Carlo random sampling method (specifically, Latin Hypercube Sampling), we obtain s sets of sample values for the random variables \( \mathbf{X} \), namely \( [\mathbf{x}_1, \mathbf{x}_2, \ldots, \mathbf{x}_i, \ldots, \mathbf{x}_s] \). Subsequently, through a series of deterministic experiments, we obtain a set of output response sample values \( [y_1, y_2, \ldots, y_i, \ldots, y_s] \). These are fitted into a polynomial function as shown in Equation (1), which then replaces the complex finite element model for probabilistic analysis.
$$ Y = a_0 + \sum_{i=1}^{r} b_i X_i + \sum_{i=1}^{r} \sum_{j=1}^{r} c_{ij} X_i X_j \tag{1} $$
In Equation (1), \( a_0 \), \( b_i \), and \( c_{ij} \) (where \( i = 1, \ldots, r; j = 1, \ldots, r \)) are the coefficients for the constant, linear, and quadratic terms, respectively. Equation (1) is referred to as the polynomial response surface function.
This paper employs the Box-Behnken matrix sampling method. For each random variable, three levels are selected. The sample points are located at the center and edge centers according to a specific rule. The level points \( q^n \) for \( X_i \) satisfy the following condition:
$$ \int_{-\infty}^{q^n} f(q) dq = p_n \quad (n = 1, 2, 3) \tag{2} $$
Here, \( f(q) \) is the probability density function of the random variable, and \( p_n \) represents the level points, with values \( p_1 = 0.01 \), \( p_2 = 0.5 \), and \( p_3 = 0.99 \). The variable \( q \) is the normally distributed variable, satisfying:
$$ q^n = \mu + \sigma \psi^{-1}(p_n) \tag{3} $$
In Equation (3), \( \mu \) is the mean, \( \sigma \) is the standard deviation, \( \psi(\cdot) \) is the standard normal distribution function, and \( \psi^{-1}(p_n) \) can be obtained from standard normal distribution tables.
Regression analysis is performed on the s sets of sample values for the random variables and output response to obtain:
$$ s = \sum_{i=1}^{s} \left[ y_i – \left( a_0 + \sum_{i=1}^{r} b_i X_{ii} + \sum_{i=1}^{r} \sum_{j=1}^{r} c_{ij} X_{ii} X_{ji} \right) \right]^2 \tag{4} $$
Solving the following partial derivative equations:
$$ \frac{\partial s}{\partial a_0} = 0; \quad \frac{\partial s}{\partial b_i} = 0; \quad \frac{\partial s}{\partial c_{ij}} = 0 \quad (i = 1, \ldots, r; j = i, \ldots, r) $$
We can estimate the undetermined coefficients \( a_0 \), \( b_i \), and \( c_{ij} \) in Equation (1), thereby determining the functional relationship between the system output response \( Y \) and the random variables \( \mathbf{X} \).
The hybrid response surface method combines the accuracy of the Monte Carlo method with the efficiency of the response surface method, effectively improving the precision and efficiency of reliability analysis for complex structures.
Nonlinear Analysis of Straight Spur Gear
A straight spur gear from a certain aero-engine reducer is selected for analysis. The model of the driving and driven gears is appropriately simplified according to the needs of the analysis. To improve computational efficiency, a surface model is used instead of a solid model to reflect the contact stress distribution during the meshing of two gear teeth. The finite element model is established using ANSYS software. The PLANE182 surface element is selected for meshing the gear model, resulting in 13,407 nodes and 11,344 elements.
The time domain [0, 1] (in seconds) is selected as the analysis range, covering the entire meshing process from the start to the completion of engagement for a pair of teeth on the driving and driven gears. The meshing process within this time domain is treated as a nonlinear problem. The material for both the driving and driven gears is 20CrMnTi carburized steel, with a yield strength of 850 MPa (according to GB/T3077-1988). The gear dimensional parameters and material properties are detailed in Table 1. By inserting a contact pair between the two gear pairs for dynamic contact analysis, the contact stress contour of the meshing gears is obtained. During the meshing process of the two gear pairs, the maximum contact stress is 813.859 MPa.
Table 1. Random Variables and Their Statistical Characteristics
| Random Variable | Mean \( \mu \) | Standard Deviation \( \delta \) | Distribution Type |
|---|---|---|---|
| Modulus \( m_n \) (mm) | 8 | 0.16 | Normal Distribution |
| Number of Teeth of Driving Gear \( Z_1 \) | 21 | 0.42 | Normal Distribution |
| Number of Teeth of Driven Gear \( Z_2 \) | 36 | 0.72 | Normal Distribution |
| Helix Angle \( \beta \) (°) | 13.6 | 0.272 | Normal Distribution |
| Pressure Angle \( \alpha \) (°) | 20 | 0.40 | Normal Distribution |
| Elastic Modulus \( m \) (MPa) | 2.12e5 | 4.24e3 | Normal Distribution |
| Poisson’s Ratio \( p \) | 0.298 | 0.00596 | Normal Distribution |
| Density \( d \) (kg/m³) | 7860 | 157.2 | Normal Distribution |
Reliability Sensitivity Analysis of Straight Spur Gear Contact Stress
Reliability Sensitivity Mathematical Model
Assuming the maximum ultimate stress of the gear pair is \( \sigma_{max} \), the limit state function for contact stress is:
$$ h(\mathbf{X}) = \sigma_{max} – Y = \sigma_{max} – \left( a_0 + \sum_{i=1}^{r} b_i X_i + \sum_{i=1}^{r} \sum_{j=1}^{r} c_{ij} X_i X_j \right) \tag{5} $$
From the above equation, \( h(\mathbf{X}) \leq 0 \) indicates a failure mode, while \( h(\mathbf{X}) > 0 \) indicates a safe mode. If the random variables are independent of each other, the mean vector and variance matrix are \( \boldsymbol{\mu} = [\mu_1, \mu_2, \ldots, \mu_r] \) and \( \mathbf{D} = [D_1, D_2, \ldots, D_r] \), respectively. The expectation and variance for the random variables are:
$$
\begin{aligned}
E[X_i^2] &= D_i + \mu_i^2 \\
E[X_i X_j] &= \mu_i \mu_j \\
D[X_i^2] &= 2D_i^2 + 4\mu_i^2 D_i \\
D[X_i X_j] &= D_i D_j + D_i \mu_j^2 + D_j \mu_i^2
\end{aligned} \tag{6}
$$
The mean and variance of \( h(\mathbf{X}) \) can be expressed as functions of the means and variances of the random variables:
$$
\begin{aligned}
\mu_h &= E[h(X_1, X_2, \ldots, X_r; \mu_1, \mu_2, \ldots, \mu_r, D_1, D_2, \ldots, D_r)] \\
D_h &= D[h(X_1, X_2, \ldots, X_r; \mu_1, \mu_2, \ldots, \mu_r, D_1, D_2, \ldots, D_r)]
\end{aligned}
$$
If \( h(\mathbf{X}) \) follows a normal distribution, the reliability index \( \beta \) and reliability \( R \) of the limit state equation are:
$$ \beta = \frac{\mu_h}{\sqrt{D_h}}, \quad R = \psi(\beta) \tag{7} $$
Here, \( \psi(\cdot) \) can be evaluated using the Monte Carlo method. The reliability sensitivity with respect to the mean and variance matrices is:
$$
\begin{aligned}
\frac{\partial R}{\partial \boldsymbol{\mu}^T} &= \frac{\partial R}{\partial \beta} \left( \frac{\partial \beta}{\partial \mu_h} \frac{\partial \mu_h}{\partial \boldsymbol{\mu}^T} + \frac{\partial \beta}{\partial D_h} \frac{\partial D_h}{\partial \boldsymbol{\mu}^T} \right) \\
\frac{\partial R}{\partial \mathbf{D}^T} &= \frac{\partial R}{\partial \beta} \left( \frac{\partial \beta}{\partial \mu_h} \frac{\partial \mu_h}{\partial \mathbf{D}^T} + \frac{\partial \beta}{\partial D_h} \frac{\partial D_h}{\partial \mathbf{D}^T} \right)
\end{aligned} \tag{8}
$$
Where:
$$
\frac{\partial R}{\partial \beta} = \psi'(\beta); \quad \frac{\partial \beta}{\partial \mu_h} = \frac{1}{\sqrt{D_h}}; \quad \frac{\partial \beta}{\partial D_h} = -\frac{1}{2} \mu_h D_h^{-3/2}
$$
The gradient vectors of \( \mu_h \) and \( D_h \) are defined as:
$$
\begin{aligned}
\frac{\partial \mu_h}{\partial \boldsymbol{\mu}^T} &= \left[ \frac{\partial \mu_h}{\partial \mu_1}, \frac{\partial \mu_h}{\partial \mu_2}, \ldots, \frac{\partial \mu_h}{\partial \mu_r} \right]^T \\
\frac{\partial \mu_h}{\partial \mathbf{D}^T} &= \left[ \frac{\partial \mu_h}{\partial D_1}, \frac{\partial \mu_h}{\partial D_2}, \ldots, \frac{\partial \mu_h}{\partial D_r} \right]^T \\
\frac{\partial D_h}{\partial \boldsymbol{\mu}^T} &= \left[ \frac{\partial D_h}{\partial \mu_1}, \frac{\partial D_h}{\partial \mu_2}, \ldots, \frac{\partial D_h}{\partial \mu_r} \right]^T \\
\frac{\partial D_h}{\partial \mathbf{D}^T} &= \left[ \frac{\partial D_h}{\partial D_1}, \frac{\partial D_h}{\partial D_2}, \ldots, \frac{\partial D_h}{\partial D_r} \right]^T
\end{aligned} \tag{9}
$$
Reliability Analysis
Based on the nonlinear analysis of the driving and driven gears, a probabilistic analysis of the contact stress during the meshing process of the gear pair is performed. Considering the key factors influencing the contact stress distribution, the gear’s critical dimensional parameters (modulus, number of teeth, helix angle, pressure angle) and material properties (elastic modulus, Poisson’s ratio, density) are set as mutually independent random variables. Their statistical characteristics are listed in Table 1.
The distribution characteristics of each random variable are expressed in the form shown in Equation (10):
$$
\begin{aligned}
m_n’ &= 6.37332 \cdot m_n – 50.9865 \\
Z_1′ &= 2.36683 \cdot Z_1 – 49.7128 \\
Z_2′ &= 1.35352 \cdot Z_2 – 48.7195 \\
\beta’ &= 3.55913 \cdot \beta – 48.3987 \\
\alpha’ &= 2.47221 \cdot \alpha – 49.4346 \\
m’ &= 2.37379 \times 10^{-4} \cdot m – 50.3243 \\
p’ &= 1.66354 \times 10^{-2} \cdot p – 49.5723 \\
d’ &= 6.43615 \times 10^{-4} \cdot d – 50.5859
\end{aligned} \tag{10}
$$
Using the Monte Carlo method and the Latin Hypercube Sampling method, 80 sets of random variable and output response sample points are obtained. These sample points are used to fit the response surface functions. To investigate the effect of the form of the response surface function (linear function, quadratic function without cross terms, quadratic function with cross terms) on the distribution characteristics of the output response, the linear response surface function, the quadratic response surface function without cross terms, and the quadratic response surface function with cross terms are fitted, as shown in Equations (11), (12), and (13), respectively:
$$ Y_1 = 813.863 – 5.05715 \times 10^{-3} \cdot m_n’ + 16.1654 \cdot Z_1′ – 0.759638 \cdot p’ \tag{11} $$
$$ Y_2 = 813.859 – 3.94320 \times 10^{-3} \cdot m_n’ + 16.1647 \cdot Z_1′ – 0.759602 \cdot p’ + 4.08190 \times 10^{-3} \cdot p’^2 \tag{12} $$
$$
\begin{aligned}
Y_3 = & 813.860 + 16.1724 \cdot m’ – 0.756951 \cdot p’ + 2.48338 \times 10^{-5} \cdot p’^2 – 1.85084 \times 10^{-5} \cdot m’ \cdot p’ \\
& + 1.88130 \times 10^{-5} \cdot Z_1′ \cdot p’ – 2.01456 \times 10^{-5} \cdot d’ \cdot p’ – 3.07290 \times 10^{-5} \cdot m’ \cdot p’^2 – 1.50093 \times 10^{-5} \cdot \beta’ \cdot p’
\end{aligned} \tag{13}
$$
Using these response surface functions, further simulations are performed on the meshing gears, yielding 10,000 simulation output samples. Based on the response surface functions from Equations (11), (12), and (13), the corresponding output sample histories and output response probability distributions are obtained.
With \( \sigma_{max} = 850 \) MPa, the reliability corresponding to the different forms of the response surface functions is obtained using Equations (5), (11), (12), and (13), as shown in Table 2. At the same time, the Monte Carlo method is used to sample the output response 10,000 times as a reference for evaluating the accuracy of different response surface functions in fitting the relationship between random variables and the output response. The results are also presented in Table 2.
Table 2. Contact Stress Distribution Characteristics and Reliability for 10,000 Samples from Different Response Surface Functions
| Response Surface Function | Mean (MPa) | Standard Deviation (MPa) | Reliability | Time (s) |
|---|---|---|---|---|
| Monte Carlo Method | 813.81 | 16.292 | 98.6341% | 126,000 |
| Linear Function | 813.86 | 16.303 | 98.6112% | 1,036 |
| Quadratic Function without Cross Terms | 813.86 | 16.302 | 98.6116% | 1,035 |
| Quadratic Function with Cross Terms | 813.86 | 16.310 | 98.6045% | 1,038 |
From Table 2, it is observed that the quadratic response surface function without cross terms exhibits high accuracy. The contact stress characteristics and reliability obtained from this method are very close to those from the Monte Carlo method. This indicates that during the meshing process of the driving and driven gears, the reliability of the contact stress for the straight spur gear is 98.6116%, which meets the design requirements.
Sensitivity Analysis
The results of the reliability analysis demonstrate that the quadratic response surface function without cross terms can effectively represent the relationship between the random variables and the output response. Using Equations (5) and (12), the influence of each random variable on the contact stress is determined, as presented in Table 3.
Table 3. Sensitivity and Probability of Random Variables
| Variable | Sensitivity (×10⁻³) | Probability (%) |
|---|---|---|
| Elastic Modulus \( m \) (MPa) | 998.82 | 88.15 |
| Poisson’s Ratio \( p \) | -62.419 | 5.51 |
| Modulus \( m_n \) (mm) | -23.571 | 2.08 |
| Helix Angle \( \beta \) (°) | 16.485 | 1.45 |
| Density \( d \) (kg/m³) | -12.453 | 1.10 |
| Number of Teeth of Driven Gear \( Z_2 \) | 6.9679 | 0.61 |
| Number of Teeth of Driving Gear \( Z_1 \) | -6.8919 | 0.61 |
| Pressure Angle \( \alpha \) (°) | -5.4633 | 0.48 |
From Table 3, it is clear that the elastic modulus \( m \) is the primary factor influencing the contact stress of the straight spur gear pair, with an influence probability of 88.15%. The impact of other random variables on the contact stress distribution is relatively minor. This conclusion is consistent with actual experimental results, providing a theoretical foundation for the optimization and redesign of gear structures.
Conclusion
Through the analysis of the nonlinear contact stress during the meshing process of a straight spur gear pair, the distribution characteristics of the contact stress were initially investigated. The maximum contact stress was selected as the output response for the probabilistic analysis of the gear pair’s contact stress. The analysis of the contact stress probability and the sensitivity of random variables for the straight spur gear pair reveals the following: For the maximum ultimate stress \( \sigma_{max} = 850 \) MPa, the probability of safety is 98.6116%, which essentially meets the design requirements. Furthermore, the key factors influencing the contact stress distribution were identified, providing a theoretical basis for the optimization of gear structures. By comparing different forms of response surface functions with the Monte Carlo method, the accuracy and efficiency of the hybrid response surface method, which combines the Monte Carlo method with a quadratic function without cross terms, were evaluated. This demonstrates the effectiveness of the hybrid response surface method and provides an effective approach for the nonlinear probabilistic analysis of the output response of complex structures, particularly in the context of straight spur gear analysis.