Random Vibration Fatigue Life Analysis of Aeroengine Bolts Using the Frequency-Domain Method

The power spectral density function S_x(ω) and its autocorrelation function of any sample function x(t) in a random process X(t) together form a Fourier transform pair: 1Sx(ω)=∫−∞∞Rx(τ)e−jωτdτRx(τ)=12π∫−∞∞Sx(ω)ejωτdω\matrix{ {{S_x}(\omega )} \hfill & { = \int_{ - \infty }^\infty {{R_x}} (\tau ){e^{ - j\omega \tau }}d\tau } \hfill \cr {{R_x}(\tau )} \hfill & { = {1 \over {2\pi }}\int_{ - \infty }^\infty {{S_x}} (\omega ){e^{j\omega \tau }}d\omega } \hfill \cr }

Among them, the autocorrelation function R_x(τ) is defined as follows: 2Rx(τ)=E[ x(t)·x(t+τ) ]{R_x}(\tau ) = E\left[ {x(t)\cdotx(t + \tau )} \right]

The correlation function describes the time-varying characteristics of a stationary stochastic process, also referred to as its time-domain characteristics. The power spectral density (PSD) function, on the other hand, describes the frequency-dependent characteristics of a stationary stochastic process, also referred to as its frequency-domain characteristics.

The PSD function can be expressed as either a bilateral (two-sided) or unilateral (one-sided) function. Formula (1) defines the bilateral PSD, S_x(ω). In practical engineering applications, because frequency is considered only in the non-negative range, and due to the even-function property of the bilateral PSD, the amplitude of the negative frequency range is symmetrically mapped to the positive frequency range. This leads to the unilateral power spectral density function, G(ω), expressed as: 3G(ω)={ 2Sx(ω),ω>0Sx(ω),ω=0G(\omega ) = \left\{ {\matrix{ {2{S_x}(\omega ),\omega > 0} \cr {{S_x}(\omega ),\omega = 0} \cr } } \right.

The frequency unit commonly used in engineering practice is Hz, rather than rad/s as in mathematics. By letting f = ω / 2π, the form of the power spectral density function form commonly used in engineering can be obtained.

The n-th moment of inertia of the power spectral density function is defined as follows: 4mn=∫0∞fnG(f)df,n=1,2,…{m_n} = \int_0^\infty {{f^n}} G(f)df,n = 1,2, \ldots

When n = 0, then 5m0=∫0∞G(f)df{m_0} = \int_0^\infty G (f)df

Formula (5) shows that the 0-order moment of inertia m₀ is the area enclosed below the power spectral density function curve. Thus, the Root Mean Square (RMS) value of the random process can be obtained.

The moment of inertia also characterizes the important time-domain characteristics of stochastic processes x(t). The relationship between some time-domain parameters and their derivatives and the moment of inertia in the frequency domain is as follows: 6m0=σx2=E[ x2(t) ]m2=σx˙2=E[ x˙2(t) ]m4=σx¨2=E[ x¨2(t) ]\matrix{ {{m_0}} \hfill & { = \sigma _x^2 = E\left[ {{x^2}(t)} \right]} \hfill \cr {{m_2}} \hfill & { = \sigma _{\dot x}^2 = E\left[ {{{\dot x}^2}(t)} \right]} \hfill \cr {{m_4}} \hfill & { = \sigma _{\ddot x}^2 = E\left[ {{{\ddot x}^2}(t)} \right]} \hfill \cr }

In the frequency domain, the above moments of inertia are commonly used to approximate the zero mean forward crossing frequency E[0] and the peak frequency E[P]. 7E[0]=m2m0E[0] = \sqrt {{{{m_2}} \over {{m_0}}}} 8E[P]=m4m2E[P] = \sqrt {{{{m_4}} \over {{m_2}}}}

The bandwidth factor α_n is also an important parameter for describing the power spectral density function, defined as: 9αn=mnm0m2n{\alpha _n} = {{{m_n}} \over {\sqrt {{m_0}{m_{2n}}} }}

As shown in formula (9), α_n is a dimensionless number, 0 ≤ α_n ≤ 1, where n can take a non-integer value.

Specifically, when n = 2 in formula (9), the spectral irregularity factor γ and spectral width coefficient ε that describe the bandwidth distribution of the power spectral density function with frequency can be obtained as follows: 10γ=αn=m2m0m4=E[0]E[P]∈(0,1)\gamma = {\alpha _n} = {{{m_2}} \over {\sqrt {{m_0}{m_4}} }} = {{E[0]} \over {E[P]}} \in (0,1) 11ε=1−γ2∈(0,1)\gamma = {\alpha _n} = {{{m_2}} \over {\sqrt {{m_0}{m_4}} }} = {{E[0]} \over {E[P]}} \in (0,1)

When ε → 0 (or γ → 1), it is represented as an ideal narrowband process (similar to a harmonic process). When ε → 1 (or γ → 0), it is represented as an ideal broadband process (similar to a white noise process).

Assuming the critical stress level is S, the expected critical crossing frequency is expressed as E[S]. The peak distribution function of a stationary Gaussian stochastic process at time t can be written as: 12F[S]=1−E[S]E[P]F[S] = 1 - {{E[S]} \over {E[P]}}

The probability density function of the peak can be obtained by taking the differential equation of formula (12), as follows: 13p(S)=−1E[P]∂∂SE[S]p(S) = - {1 \over {E[P]}}{\partial \over {\partial S}}E[S]

Formula (13) can be used to obtain a more general representation of the peak probability density function, which is: 14p(S)=12πσ(1−γ2)12exp(−S2{ 2σ2(1−γ2) }−1)+S2σ2γ[ 1+erf{ Sσ(2γ−2−2)−12 } ]exp(−S22σ2)p(S) = {1 \over {\sqrt {2\pi } \sigma }}{\left( {1 - {\gamma ^2}} \right)^{{1 \over 2}}}\exp \left( { - {S^2}{{\left\{ {2{\sigma ^2}\left( {1 - {\gamma ^2}} \right)} \right\}}^{ - 1}}} \right) + {S \over {2{\sigma ^2}}}\gamma \left[ {1 + erf\left\{ {{S \over \sigma }{{\left( {2{\gamma ^{ - 2}} - 2} \right)}^{ - {1 \over 2}}}} \right\}} \right]\exp \left( { - {{{S^2}} \over {2{\sigma ^2}}}} \right)

Where erf(x) is the error function: 15erf(x)=2π∫0xexp(−t2)dterf(x) = {2 \over {\sqrt \pi }}\int_0^x {\exp } \left( { - {t^2}} \right)dt

When γ → 1, the peak probability density function of formula (14) can be expressed as: 16p(S)=Sσ2exp(−S22σ2)p(S) = {S \over {{\sigma ^2}}}\exp \left( { - {{{S^2}} \over {2{\sigma ^2}}}} \right)

Where 0 ≤ S < ∞. This peak probability density function exhibits a Rayleigh distribution.

When γ → 0, the peak probability density function of formula (14) can be expressed as: 17p(S)=12πσexp(−S22σ2)p(S) = {1 \over {\sqrt {2\pi } \sigma }}\exp \left( { - {{{S^2}} \over {2{\sigma ^2}}}} \right) where –∞ < S < ∞. The peak probability density function exhibits a Gaussian distribution.

According to the theory of structural fatigue damage, it is the amplitude of stress cycles – rather than the peak value – that causes structural damage. Therefore, in fatigue life analysis, the first step is to obtain the corresponding stress amplitude probability density function based on the above stress peak probability density function.

When the probability density function of the amplitude of the stress cycle is obtained, it can be concluded that n (S_i) in the Miner linear damage accumulation formula is: 18n(Si)=E[P]·T·p(Si)·ΔSin\left( {{S_i}} \right) = E[P]\cdotT\cdotp\left( {{S_i}} \right)\cdot\Delta {S_i}

Where n (S_i) represents the actual number of stress cycles during time T when the stress amplitude is S_i, and p (S_i) represents the probability density function when the stress amplitude is S_i.

By substituting formula (18) into the S-N curve equation, the estimated value of damage can be obtained as follows: 19D=E[P]TC−1∑i=1+∞Simp(Si)·ΔSiD = E[P]T{C^{ - 1}}\sum\limits_{i = 1}^{ + \infty } {S_i^m} p\left( {{S_i}} \right)\cdot\Delta {S_i}

When ΔS_i → 0, formula (19) can be written in the following integral form: 20D=E[P]TC−1∫0+∞Smp(S)dSD = E[P]T{C^{ - 1}}\int_0^{ + \infty } {{S^m}} p(S)dS

Under random loading, the sequential effect of load cycles is generally negligible. It is commonly assumed that fatigue failure occurs when the accumulated damage reaches a critical value, D_CR = 1. Therefore, the fatigue life of the structure can be derived as: 21T=CE[P]∫0+∞Smp(S)dST = {C \over {E[P]\int_0^{ + \infty } {{S^m}} p(S)dS}}

As shown in formula (21), the key to estimating the fatigue life of a structure using frequency domain analysis is to convert the power spectral density function of the stress response spectrum into the amplitude probability density function of the stress cycle.

In engineering, it is generally believed that when the spectral width parameter ε < 0.3, the stochastic process can be regarded as a narrowband process, and its amplitude probability density function is generally approximated as a Rayleigh distribution; When the spectral width parameter ε > 0.3, the random process is considered a broadband process, and its amplitude probability density function is relatively complex, with multiple descriptive models.

A system with an irregular factor close to 1 has a load time history that looks like an amplitude-modulated sine wave. In fact, there are no troughs above the average value of the signal, so for any given level S, the expected number of peaks above S is equal to the expected value of the positive intersection with S level. Therefore, for narrowband processes, the estimation of the fraction of cycles with peaks greater than x(t) = S is given by the following equation: 22pp(S)=pS(τ)p(τ)e−S2/2Rx(0){p_p}(S) = {{{p_S}(\tau )} \over {p(\tau )}}{e^{ - {S^2}/2{R_x}(0)}}

The probability [1 – p_p(S)] that the peak value is less than S is defined as the distribution function of the peak value, and the derivative of this quantity yields the peak probability density function: 23pp(S)=SR(0)e−S2/2Rx(0){p_p}(S) = {S \over {R(0)}}{e^{ - {S^2}/2{R_x}(0)}}

Therefore, for the special case of narrowband processes, the probability density function of the peak amplitude obtained will be a Rayleigh distribution function.

Since the peak frequency of the narrowband process is similar to the zero mean forward crossing frequency (i.e. the center frequency of the narrowband process), it can be assumed that E[P] = E[0]. Substituting formula (23) into formulas (19) and (20) yields the fatigue damage and life as follows: 24DNB=C−1σ−2E[0]·∫0+∞Sm+1e−S2/2σ2dS{D_{NB}} = {C^{ - 1}}{\sigma ^{ - 2}}E[0]\cdot\int_0^{ + \infty } {{S^{m + 1}}} {e^{ - {S^2}/2{\sigma ^2}}}dS 25TNB=Cσ2E[0]·∫0+∞Sm+1e−S2/2σ2dS{T_{NB}} = {{C{\sigma ^2}} \over {E[0]\cdot\int_0^{ + \infty } {{S^{m + 1}}} {e^{ - {S^2}/2{\sigma ^2}}}dS}}

When the expected damage is equal to 1, fatigue failure occurs.

Dirlik conducted simulations using the Monte Carlo method to study the power spectral density functions of 70 different spectral types. In Dirlik (1985), he proposed a quasi-empirical model that no longer uses a Rayleigh distribution to represent the probability density function of stress amplitude in broadband processes. It is believed that the probability density function of stress amplitude in broadband processes is a combination of an exponential distribution and two Rayleigh distributions. Then, based on the minimum regularization error between the model and the rainfall amplitude, relevant fitting parameters are extracted. The Dirlik equation is: 26p(S)=1σ(D1Qe−ZQ+D2ZR2e−Z22R2+D3Ze−Z22)p(S) = {1 \over \sigma }\left( {{{{D_1}} \over Q}{e^{ - {Z \over Q}}} + {{{D_2}Z} \over {{R^2}}}{e^{ - {{{Z^2}} \over {2{R^2}}}}} + {D_3}Z{e^{ - {{{Z^2}} \over 2}}}} \right)

Where Z = S/σ is the regularized stress amplitude, and the other parameters in the equation are: D1=2(xm−γ2)1+γ2D2=1−γ−D1+D121−RD3=1−D1−D2Q=1.25(γ−D3−D2R)D1=1.25D1R=γ−xm−D121−γ−D1+D12xm=m1m0m2m4\matrix{ {{D_1}} \hfill & { = {{2\left( {{x_m} - {\gamma ^2}} \right)} \over {1 + {\gamma ^2}}}} \hfill \cr {{D_2}} \hfill & { = {{1 - \gamma - {D_1} + D_1^2} \over {1 - R}}} \hfill \cr {{D_3}} \hfill & { = 1 - {D_1} - {D_2}} \hfill \cr Q \hfill & { = {{1.25\left( {\gamma - {D_3} - {D_2}R} \right)} \over {{D_1}}} = 1.25{D_1}} \hfill \cr R \hfill & { = {{\gamma - {x_m} - D_1^2} \over {1 - \gamma - {D_1} + D_1^2}}} \hfill \cr {{x_m}} \hfill & { = {{{m_1}} \over {{m_0}}}\sqrt {{{{m_2}} \over {{m_4}}}} } \hfill \cr }

As can be seen from the above equation, the Dirlik equation involves four moments of inertia of the power spectral density function. Although the Dirlik equation did not take into account the influence of the mean values of each stress cycle, practice has shown that it usually provides relatively accurate life estimation results.

By substituting formula (26) into formulas (19) and (20), the fatigue damage and life can be obtained, and then expressed in the form of a Gamma function through mathematical transformation, that is: 27DDK=E[P]C·σm·[ D1QmΓ(1+m)+(2)mΓ(1+m2)·(D2|R|m+D3) ]{D_{DK}} = {{E[P]} \over C}\cdot{\sigma ^m}\cdot\left[ {{D_1}{Q^m}\Gamma (1 + m) + {{(\sqrt 2 )}^m}\Gamma \left( {1 + {m \over 2}} \right)\cdot\left( {{D_2}|R{|^m} + {D_3}} \right)} \right] 28TDK=CE[P]·σm·[ D1QmΓ(1+m)+(2)mΓ(1+m2)·(D2|R|m+D3) ]{T_{DK}} = {C \over {E[P]\cdot{\sigma ^m}\cdot\left[ {{D_1}{Q^m}\Gamma (1 + m) + {{(\sqrt 2 )}^m}\Gamma \left( {1 + {m \over 2}} \right)\cdot\left( {{D_2}|R{|^m} + {D_3}} \right)} \right]}} where C is the fatigue strength of the engine bolt material, m is the fatigue index of the engine bolt material, and Γ is the Gamma function.

The environmental spectrum of random engine vibration consists of a combination of broadband and narrowband components, as shown in Figure 5. The acceleration power spectral density (PSD) function is defined as follows: L1 is taken as 1.0 g²/Hz, F1 as 188 Hz, F2 = 2F1, F3 = 3F1, F4 = 4F1, Peak bandwidth = ± 5%. Note: when inputting this curve, the acceleration PSD dimension should be converted from g²/Hz to (mm/s²) 2/Hz, that is, multiplying the values in the figure by 98,002 = 9.604 × 107. Apply acceleration excitation to the X, Y, and Z axes during analysis and calculation.

Figure 5.

Vibration excitation spectrum

Following the procedure outlined in the diagram, a dynamic analysis of the structure was carried out, consisting sequentially of modal analysis, frequency response analysis, and random response analysis. The root mean square (RMS) stress results for the connecting bolts are shown in Figure 6.

Figure 6.

RMS stress cloud maps of connecting bolts: (a) X-direction, (b) Y-direction, (c) Z-direction.

The stress power spectral density (PSD) curves at the critical points of the structure under different excitation conditions are presented in the following figure.

Figure 7.

Stress PSD curve of dangerous points in the X direction.

Figure 8.

Stress PSD curve of dangerous points in the Y direction.

Figure 9.

Stress PSD curve of dangerous points in the Z direction.

The S–N curve of TC4 material is shown in Fig. 10, where the fatigue strength coefficient is 1973 Mpa, the first fatigue strength exponent is -0.104, and the fatigue limit is 468.95 Mpa.

Figure 10.

S–N curve of TC4 material.

Under random vibration analysis, the RMS stress in the Y-direction was found to be the highest, resulting in the shortest fatigue life. Under these most critical conditions, the estimated fatigue life is 1.598 × 10⁴ seconds, or approximately 4.44 hours, as shown in the figure. This service life is relatively low and does not meet design requirements. Therefore, a bolt model with higher fatigue strength should be selected during the design stage.

Figure 11.

Fatigue life results of connecting bolts.