Spiral bevel gears are essential components in various mechanical systems, and their reliability is crucial for the proper functioning of the entire system. The reliability analysis of spiral bevel gears is a complex task due to the presence of various uncertainties in the design parameters, boundary conditions, and external environment. In this article, we will explore the structural reliability analysis of spiral bevel gears based on hybrid uncertainties, including random and interval uncertainties.

**Introduction**

The performance of spiral bevel gears can vary or fluctuate due to uncertainties related to the external environment, material properties, geometric dimensions, and loads. Understanding, measuring, and controlling these uncertainties are essential for ensuring the reliability of spiral bevel gears. In practical engineering, traditional reliability modeling and analysis methods based on probability theory may not be applicable due to the presence of imprecise and non-continuous uncertainty variables. Therefore, it is necessary to develop new methods for analyzing the reliability of spiral bevel gears under hybrid uncertainties.

**Classical Reliability Analysis Methods**

In the reliability design and analysis of mechanical structures, the limit state function G(X) is constructed based on the specified functions. X = (x1, x2,…, xn)T represents the independent random variables that affect the function, such as loads, material properties, geometric shapes, sizes, and environmental factors. G(X) determines whether the mechanical structure works normally. If G(X) > 0, it indicates safe operation; otherwise, it indicates failure. The probability of failure Pf can be defined as:

Pf = Pr{G(X) ≤ 0}

where Pr{•} is the probability measure. When G(X) > 0, the reliability probability can be expressed as:

Preliability = 1 – Pr(G(X) ≤ 0) = Pr(G(X) > 0)

This can be calculated by the integral:

Pr(G(X) > 0) = ∫G(X)>0 fX(X) dX

where fX(X) is the joint probability density function of the random variables X.

In the First and second order reliability methods (FORM), the random variables X in the original space are transformed into the variables u in the standard normal space, and the limit state equation G(X) is transformed into:

g(u) = G(X) = G[T(u)]

where T(u) is the transformation from the original space to the standard normal space.

Therefore, the integral can be rewritten as:

Pr(g(u) > 0) = ∫g(u)>0 fu(u) du

where fu(u) is the n-dimensional standard Gaussian density function.

Using the reliability index method, the equation can be further transformed into the following optimization problem:

min β = ||u*||

s.t. g(u) = 0

where || • || represents the norm of the vector, u* is the optimal design point (the design point where the joint probability density function reaches the maximum on the limit state function), and β = ||u*|| is the reliability index (the minimum geometric distance from the origin to the limit state surface in the space).

The first-order Taylor expansion of the limit state function about the mean point u+ can be written as:

g(u) ≈ g(u+) + ∇g(u+)(u – u+)

where ∇g(u+) is the value of ∇g(u) at the point u+.

In the FORM method, the optimal design point is the solution to Equation (6), and it can be expressed as:

u* = βα

where α = ∇g(u*)/||∇g(u*)|| is the unit vector of the design point.

**Hybrid Reliability Analysis Method Based on SORM**

**Polar Coordinate Transformation Method**

In the analysis of mechanical structures with random and interval uncertainties, the limit state equation is g(x, y), where x = (x1, x2,…, xn) is an n-dimensional random variable vector, and y = (y1, y2,…, ym) is an m-dimensional interval variable vector. In the standard normal space, the random variable vector can be represented as u = (u1, u2,…, un), and the normalized independent interval vector can be represented as δ = (δ1, δ2,…, δm) ∈ [-1, 1]. The limit state equation can then be expressed as:

g(x, y) = g[T(u), T(δ)]

where T(u) and T(δ) are the transformations from the original space to the standard normal space.

The transformation from the (u, δ) space to the polar coordinate space can be achieved through two nonlinear features:

v1 = ||(u, δ)||2 = √(∑i=1n ui^2 + ∑j=1m δi^2)

v2 = [(u, δ), a]/(||(u, δ)|| ||a||) = [(u, δ), a]/v1

v1 follows the chi-square distribution:

φ1(v1) = {2^(1 – n/2) v1^(n – 1) / Γ(n/2) e^(-v1^2/2), v1 > 0; 0, v1 ≤ 0}

v2 follows the following distribution:

φ2(v2) = {sin^(n – 2)(arccos v2) + sin^(n – 2)(π – arccos v2) / √(1 – v2^2) ∫0^π sin^(n – 2) α dα, 0 < v2 < 1, n ≥ 2; sin^(n – 2)(arccos(-v2)) + sin^(n – 2)(π + arccos v2) / √(1 – v2^2) ∫0^π sin^(n – 2) α dα, -1 < v2 ≤ 0, n ≥ 2}

These two new variables are independent. To avoid serious errors caused by subjective assumptions in describing the uncertainty characteristics, a new polar coordinate transformation with random and interval uncertainties will be introduced in the next section.

**Random and Interval Uncertainty Polar Coordinate Transformation**

The distribution of the new variables in the polar coordinate space can be derived as follows:

For ∀n, m ∈ N (N is the non-negative integer domain), when n + m ≥ 2, v2 follows the distribution:

φ2(v2) = {sin^(n + m – 2)(arccos v2) + sin^(n + m – 2)(π – arccos v2) / √(1 – v2^2) ∫0^π sin^(n + m – 2) α dα, 0 < v2 < 1; sin^(n + m – 2)(arccos(-v2)) + sin^(n + m – 2)(π + arccos v2) / √(1 – v2^2) ∫0^π sin^(n + m – 2) α dα, -1 < v2 ≤ 0}

For ∀δi ∈ [-1, 1] (i = 1, 2,…, m), it can be deduced that 0 ≤ ∑j=1m δi^2 ≤ m. Since ui (i = 1, 2,…, n) are independent standard Gaussian distribution variables, ∑i=1n ui^2 = v1^2 – ∑j=1m δi^2, and thus ∑i=1n ui^2 follows the χ^2(n) distribution, which can be specifically expressed as:

f(∑i=1n ui^2) = 2^(1 – n/2) (∑i=1n ui^2)^(n – 1) / Γ(n/2) e^(-(∑i=1n ui^2)^2 / 2)

Therefore, for ∀∑j=1m δi^2 ∈ [0, m], v1 follows the distribution:

φ1(v1) = {2^(1 – n/2) (v1^2 – ∑j=1m δi^2)^(n – 1)/2 / Γ(n/2) e^(-(v1^2 – ∑j=1m δi^2)/2), v1 ≥ √∑j=1m δi^2; 0, v1 < √∑j=1m δi^2}

The failure probability of the structural limit state function G(v1, v2) can be expressed as:

Pf = ∫-1^1 ∫-1^1 ∫√∑i=1m δi^2^+∞ φ1(v1, v2) ≤ 0) dv1 dv2

where ∑j=1m δi^2 ∈ [0, m], φ1 follows Equation (19), and φ2 follows Equation (17).

To show the upper and lower bounds of the failure probability or reliability, it only needs to prove that the cumulative distribution function Fv1(v1, ∑j=1m δi^2) of v1 for the variable ∑j=1m δi^2 is non-increasing.

To make the proof more convenient, let Δ = ∑j=1m δi^2, τ = v1^2 – ∑j=1m δi^2, τ ∈ [0, ∞]. Then the partial derivative of Fv1(v1, Δ) with respect to Δ can be expressed as:

∂Fv1(v1, Δ)/Δ = ∂Fv1(τ)/Δ = ∂Fv1(τ)/∂τ ∂τ/Δ = -τ^(n – 1)/2 / (2^(n/2) Γ(n/2)) e^(-τ/2) / √(τ + Δ)

In Equation (21), for ∀∑j=1m δi^2, e^(-τ/2) > 0, 1/√(τ + Δ) > 0, and τ^(n – 1)/2 / (2^(n/2) Γ(n/2)) ≥ 0. Obviously, ∂Fv1(v1, Δ)/Δ ≤ 0. Therefore, the cumulative distribution function of v1 for the variable Δ is non-increasing.

Since v1 and v2 are independent in the polar coordinate space, the failure probability Pf can be expressed as an interval [Pf, Pf], that is:

When ∑j=1m δi^2 = m, Pf takes the lower bound Pf.

Pf = ∫-1^1 ∫√m^+∞ φ1(v1, m) φ2(v2) dv1 dv2

When ∑j=1m δi^2 = 0, Pf takes the upper bound Pf.

Pf = ∫Ωv=[G(v1, v2)≤0]^+∞ φ1(v1, 0) φ2(v2) dv1 dv2

**Improved Hybrid Reliability Analysis Method Based on SORM**

To facilitate the description, the structural reliability analysis problem is defined in an n + m-dimensional ω space using ω = (u, δ), where u is a standard independent Gaussian vector, and δ is a normalized independent interval vector. Therefore, Equation (13) can be simplified to:

g(x, y) = g(ω)

The limit state equation is expanded using the second-order Taylor series:

g(ω) ≈ g(ω*) + ∇g(ω*)(ω – ω*) + λ(ω – ω*)^T(ω – ω*)

where λ is the average curvature. This function can be derived from Equation (7) by using the second-order Taylor expansion formula at the design point ω*.

The function can be further expanded as:

g(ω) ≈ g(ω*) – ∑k=1n+m (∂g/∂ωk|ω* ωk*) + ∑k=1n+m (∂g/∂ωk|ω* ωk) – 2λ(∑k=1n+m (∂g/∂ωk|ω* ωk)) / (||ωk*|| [∑k=1n+m (∂g/∂ωk|ω*)^2]^(1/2)) + λ∑k=1n+m (ωk*)^2 + λ∑k=1n+m (ωk)^2

In the polar coordinates, by substituting Equations (26) and (27) into Equation (28), the limit state function g(ω) can be expressed by v1 and v2:

g(v1, v2) = g(ω*) – ∑k=1n+m (∂g/∂ωk|ω* ωk*) + λv1^2 +…

where v1*^2 = ∑i=1n+m (ωi*)^2, v1* = cos∠(ω+, ω*) = 1.

Let g(v1*, v2*) = 0, then:

λ = (g(ω*) – ∑i=1n+m (∂g/∂ωk|ω* ωk*) + Dv1*)/(2 – 2v1*^2)

Thus, g(v1, v2) can be expressed as:

g(v1, v2) = d + (D – 2λ/v1*) v1v2 + λv1^2

where d = g(ω*) – ∑k=1n+m (∂g/∂ωk|ω* ωk*) + λv1*^2.

When λ > 0, the new safe domain is d + (D – 2λ/v1*) v1v2 + λv1^2 > 0, and the failure domain is d + (D – 2λ/v1*) v1v2 + λv1^2 ≤ 0.

Similar to the derivation process of Equations (22) and (23), the lower and upper bounds of the failure probability interval [Pf, Pf] of the structure can be derived:

Pf = ∫-1^1 ∫√m^+∞ φ1(v1) φ2(v2) dv1 dv2

Pf = ∫-1^1 ∫0^1 ∫0^+∞ φ1(v1) φ2(v2) dv1 dv2

where φ1 follows Equation (19) and φ2 follows Equation (17).

**Case Verification Analysis of Spiral Bevel Gear**

**Basic Parameters of Spiral Bevel Gear**

According to the qualitative analysis of the reliability of the front transmission assembly system of a certain equipment, the failure of the bevel gear has a significant impact, and structural damage caused by contact fatigue is likely to occur. Therefore, the reliability analysis of the spiral bevel gear in the front transmission assembly system is carried out. In the front transmission assembly system, the gear material is 20Cr2Ni4A, and the basic parameters of the spiral bevel gear are shown in Table 1, and the material parameters of the gear are shown in Table 2.

In the case of mass processing, since the various factors such as the environment and operation are independent and there is no dominant tendency, it is considered that their distribution obeys the normal distribution. In the reliability model of contact fatigue of spiral bevel gears, there are many factors affecting the reliability. To more intuitively verify the correctness and effectiveness of the second-order reliability analysis method based on random-interval hybrid uncertainties proposed in this article, only the tooth width b is taken as a random variable in this case, and the rotational speed n and torque T are taken as interval variables. In the gear design process, since the tooth width coefficient is a value within the range of 1.0/3.5 to 1.0/3.0, there is a tolerance range for the tooth width. Assuming that the tolerance of the tooth width is ±Δb (generally b obeys the normal distribution), its statistical parameters are: μb = b¯, σb = Δb / 3. The selection of the simulation operating conditions affects many aspects of the gear design performance. Since the torque and rotational speed are unstable in the simulation operating conditions, the designed simulation operating conditions are n ± Δn and T ± ΔT, and their numerical distributions are shown in Table 3.

**Construction of Limit State Function for Spiral Bevel Gear**

**Contact Stress**

Contact strength is a key consideration in gear design, and the strength calculation method in the standards of various countries for spiral bevel gears generally adopts the method of the Gleason Company in the United States: using the Hertz formula for the contact of cylinders as the basic form, and introducing correction parameters in terms of load, size, stress distribution, etc.

The formula for calculating the contact stress of orthogonal transmission gears is:

σH = ZH · ZE · Zc · Zβ · KA · KV · KHβ · KHα · Ft / (dml · beH · √(u^2 + 1) / u)

where σH is the calculated contact stress, ZH is the node area coefficient, ZE is the elastic coefficient, Zc is the coincidence coefficient for contact strength calculation, Zβ is the spiral angle coefficient for contact strength calculation, KA is the usage coefficient, KV is the dynamic load coefficient, KHβ is the tooth direction load distribution coefficient for contact strength calculation, KHα is the tooth direction load allocation coefficient for contact strength calculation, Ft is the nominal tangential force on the midpoint of the tooth width, dml is the end face of the midpoint of the tooth width of the small gear, beH is the effective tooth width for contact strength calculation, u is the gear ratio of the bevel gear, T1,2 is the acting torque on the small gear, and dm1,2 is the pitch diameter of the midpoint of the tooth width of the small gear.

The calculation methods for several specific values are as follows:

- Zβ is the coincidence coefficient. For the characteristics of the arc tooth in this research object, when the longitudinal coincidence degree εn > 1, its value is only related to the end face coincidence degree Eve. The formulas for the end face and longitudinal coincidence degree and the coincidence coefficient for contact strength calculation are as follows:

εv = gvaR / (πm(R – 0.5b) cosαv)

εva = 0.85bR tanβm / (πm(R – 0.5b))

Zβ = √(1 / εva) - KV is the dynamic load coefficient:

KV = NK + 1

where N is the critical rotation ratio, which is the ratio of the rotation speed n1 of the small gear to the critical rotation speed nE1:- N = 84 * z1 * vm / (10000 * √(u^2 / (u^2 + 1)))

- When N ≤ 0.85, the value of K is:
- K = ((fp1 – ya) * c’) / (KA * Ft / beH) * Cv12 + Cv3
- Here, fpt is the size limit deviation, usually taken according to the large gear; ya is the running-in amount; c’ is the stiffness of a single pair of teeth; Cv12 and Cv3 are the coefficients when N ≤ 0.85.

- The allowable stress is calculated as follows:
- The formula for calculating the allowable contact stress of orthogonal transmission gears is:
- σHp = σHlim / SHm * ZL * Zv * ZR * ZX

- The formula for calculating the allowable contact stress of orthogonal transmission gears is:
- Here, σHp is the contact fatigue strength of the test gear, σHlim is the contact fatigue limit of the test gear, SHm is the minimum safety factor for contact strength calculation, ZL is the lubricant coefficient, Zv is the speed coefficient, ZR is the roughness coefficient, and ZX is the size coefficient for contact strength calculation.

**Construction of Limit State Function**

By assigning values to the relevant parameters: ZH is taken as 2.125 based on the normal tooth shape angle and spiral angle; ZE is taken as 189.8 as both gears are made of steel; ZK is taken as 0.85 when the tooth top and root are properly modified; KA is taken as 1.21 according to the characteristics of the prime mover and the working machine, which is subject to slight impact, uniform and stable operation, and is a speed-increasing transmission; KHβ is taken as 1.5 as both wheels are supported at both ends and used in aircraft and vehicles; KHα is taken as 1.1 as the gear precision level is 6 and it is a hard tooth surface curved gear. Since there is 1 random variable and 2 interval variables in the geometric parameters, these 6 parameters will change accordingly in the calculation of contact stress.

Based on the stress-interference model, the limit state function of the spiral bevel gear can be simplified as:

g(X) = 1602 – 270.258 * Ze * √(2.369 * Kv * KHα * Ft / (dm * beH))

**Calculation Results Analysis**

To verify the calculation efficiency, correctness, and effectiveness of the method proposed in this article, the Monte Carlo sampling analysis method is first introduced in this case for reliability analysis. Most existing MC analysis methods treat interval variables as either a uniform distribution or a truncated normal distribution when dealing with reliability problems involving random and interval mixtures. In this regard, this article treats the interval variables based on both the uniform distribution and the truncated normal distribution. When the sample size is 10^7, the failure probability of this case is calculated respectively. When the interval variables are sampled and processed according to the uniform distribution, the failure probability of this case is 0.1881, and it takes 15.9 seconds. When the interval variables are sampled and processed according to the truncated normal distribution, the failure probability of this case is 0.1989, and it also takes 15.9 seconds.

Since the relevant theories mentioned earlier have proved that the probability density function of v1 is non-increasing with respect to the variable Δ, the upper and lower bounds of the failure probability of this case can be obtained at the endpoints of the interval variables. Therefore, the MC analysis method (with a sample size of 10^7) and the classic SORM method are introduced respectively to obtain the upper and lower bounds of the failure probability of this case at the endpoints of the interval variables. At the same time, the method proposed in this article is applied to obtain the upper and lower bounds of the failure probability of this case without changing the variable information. The results are shown in Table 4.

Method | Range of Failure Probability | ||||
---|---|---|---|---|---|

Lower Bound | Upper Bound | Median | Deviation | Time Taken (s) | |

MC | 0.0031 | 0.5776 | 0.2904 | 0.2873 | 4.2 |

SORM | 0.0051 | 0.5718 | 0.2904 | 0.2834 | 3.2 |

This Article | 0.1874 | 0.3827 | 0.2851 | 0.0977 | 2.3 |

Through comparison, the following can be observed:

- When using the MC method with a sample size of 10^7, the failure probability of this case obtained by sampling the interval variables according to the uniform distribution is 0.1881, and the failure probability obtained by sampling the interval variables according to the truncated normal distribution is 0.1989. Compared with the median value calculated by taking the endpoint values of the interval variables using the same MC method, the relative errors of the uniform distribution sampling and the truncated normal distribution sampling are 35.27% and 31.51% respectively, which are very large. The main reason is that treating the interval variables as different probability distributions leads to the loss of important information, resulting in unreliable probability failure analysis results.
- On the one hand, the failure probability interval of the case obtained using the method in this article is [0.1874, 0.3827], compared with the failure probability interval [0.0031, 0.5776] calculated by taking the endpoint values of the interval variables, the relative error of the median value of the failure probability interval is only 1.83%. Moreover, the failure probability interval obtained by the method in this article completely falls within the failure probability interval calculated by taking the endpoint values of the interval variables and the classic SORM method, verifying the correctness and effectiveness of the method in this article. On the other hand, in terms of the fluctuation range of the failure probability, the fluctuation deviation of the failure probability interval calculated by taking the endpoint values of the interval variables is 0.2873, while the fluctuation deviation of the failure probability interval obtained by the method in this article is only 0.0977, with a reduction of 65.99%. The comparative analysis of the above results shows that the method in this article avoids the problem of interval expansion caused by interval operations to a certain extent and narrows the range of the upper and lower bounds of the failure probability, making the results of this method more accurate and reliable, and thus more instructive for practical engineering design.
- In terms of calculation efficiency, due to the particularity of the MC method, its calculation efficiency is not compared with that of the method in this article. Therefore, the classic SORM method is introduced in this article. The relative errors of the upper and lower bounds of the failure probability calculated by the classic SORM and MC methods are very small, but the fluctuation range of the reliability is large. In terms of time consumption alone, the classic SORM takes 3.2 seconds, while the method in this article takes 2.3 seconds, which improves the efficiency of reliability analysis to a certain extent.

**Conclusion**

This article proposes a second-order reliability analysis method suitable for random-interval hybrid uncertainties. The results show that the concept of the design point in reliability analysis based on probability theory can be applied in this research. Specifically, the unit vector of the design point is the direction in which the output variable implied in the limit state function changes rapidly. To utilize this property and reduce the computational load, this article establishes a two-dimensional mapping with two nonlinear features in the polar coordinate space: ① the distance to the origin; ② the cosine of the angle and the unit vector of the design point. In this case, the hybrid reliability analysis problem is simplified into a classic reliability problem, where the n-dimensional limit state function is approximated by two random variables near the design point, obeying a new probability density function with random and interval variables. This new probability density function is different from the assumption of the uniform density function model of the interval variables, which helps to evaluate lower and higher failure probabilities.

In future research, we can further explore the application of this method in more complex mechanical systems and optimize the algorithm to improve the calculation efficiency and accuracy. Additionally, considering more factors that may affect the reliability of spiral bevel gears, such as the wear and fatigue of the gear surface, can provide a more comprehensive analysis of the reliability of the system.