Multiaxial Fatigue Analysis of Connecting Bolt at High-Speed Train Axle Box under Structural Subharmonic Resonance

Abstract

Based on the dynamic characteristics of the axle box front cover of high-speed trains in the subharmonic resonance state, the nonlinear single-degree-of-freedom (SDOF) model was proved to be reasonable, and reasons for the ineffectiveness of the common prevention methods for preventing bolt failure were analyzed firstly. Then, dynamic stress of the bolt was simulated by innovatively adopting the linear method based on frequency response analysis. The stress simulation method was verified to be practical under the subharmonic resonance state by analyzing and comparing the experimental and numerical results of the bolted front cover. It was proved that the linear method was accurate enough to simulate the dynamic stress of bolts, which is of great engineering significance. In addition to the transverse resonance stress of bolts caused by drastic vertical vibration of the front cover, the tensile resonance stress at the root of the first engaged thread was too large to be neglected on account of the first-order bending modes of bolts. Next, equivalent stress amplitude of the multiaxial stresses was obtained by means of the octahedral shear stress criterion. Finally, fatigue life of bolts was predicted in terms of S-N curve suitable for bolt fatigue life analysis. It argued that the bolts were prone to multiaxial fatigue failure when the front cover was in subharmonic resonance for more than 26.8 h, and the fatigue life of bolts could be greatly improved when the wheel polygonization was eliminated by shortening the wheel reprofiling interval.

1. Introduction

1.1. Literature Survey

Since bolted joints are widely used in high-speed trains, their reliability is of great significance for the safety operation of vehicles. Loosening and fatigue failure of bolted joints are the two most important problems for vehicles that are subjected to vibration loading [1]. Abundant studies have demonstrated that self-loosening failure of bolts is mainly caused by dynamic transverse loading [2,3,4,5,6], and fatigue failure of bolts is usually caused by dynamic tensile loading [7,8,9,10]. However, transverse vibration can cause fatigue of bolted joints as well. Hashimura [11] proved that the loosening-fatigue life of bolted joints did not depend on the initial clamping force, but it significantly depended on the transverse vibration force. Minguez [12] reported that there was bending as well as shearing of bolted joints because the load transfer path was not symmetric for single lap joints. Increasing bolt tightening torque did not have any effect on the fatigue life for this type of joint because the main cause for fatigue failure was the strong bending moment caused by distortion. Benhaddou [13] and Guo [14] demonstrated that it was easy to generate stress concentration in the bolt hole, and increasing the preload of bolt could alleviate the stress concentration and transfer the cracking position. Wang [9] studied the fatigue behavior of stainless steel bolts under transverse loading and regarded the relative displacements of clamped parts as the deformations of bolts given the rigid clamped parts. Due to the high notch effect, bolts are more prone to fatigue failure than nuts or clamped parts; thus, it is of significance to study the stress state of bolts. However, the above studies mainly discussed the fatigue behavior of clamped parts rather than the transverse capacity of bolts itself. Traditionally, axial load of bolts can be monitored accurately in real time by means of strain gauged bush [15,16], smart washer [17,18], as well as embedding the sensor into the bolt shank [19]. Through the transverse dynamic load test of bolted joints, researchers found that cracks at the root of the first bolt thread could be observed after about 10⁴–10⁶ cycles when the loosening of bolts occurred [11]. But the transverse stress of bolts was less experimentally measured due to the limitation of measuring techniques. Although stress analysis of bolted joints is feasible by FEM [6,9], it is difficult to simulate the complex excitations of the actual operating state, which would result in a large difference between the calculated stress and actual stress. Especially when the structure is in a nonlinear vibration state, it cannot be excited in commercial FEA software such as ABAQUS, so the stress simulation results are not reliable. In short, it is not suitable to simulate the multiaxial stress of bolts under complex operation states by FEM.

Polygonal wear of wheel tread occurs frequently in Chinese railways, which can induce drastic vibration or even resonance and cause fatigue failure of critical structures; it is particularly dangerous for the structures fixed by bolts when the bolt axis is perpendicular to the structural vibration direction. Vibration characteristics and fatigue life under excitation of wheel polygonization were studied a lot for critical structures, such as gearbox housing [20], bearings [21] and lifeguards [22,23], etc.; the connecting bolts are usually subjected to multiaxial dynamic stresses in the real operation condition of high-speed trains. In the past few decades, the issue of damage mechanisms and stress combination of multiaxial loading have been brought to extensive attention, and significant research progress has been made by means of theoretical research and numerical analysis [24,25,26,27]. However, multiaxial fatigue life of connecting bolts was rarely studied because the dynamic stress of bolt threads is difficult to measure. Although metallographic analysis can qualitatively analyze material fatigue failure, the quantitative analysis cannot be realized because the cyclic stress level of the structure is unable to be obtained [28,29]. Consequently, it is of great importance to find a suitable dynamic stress calculation method for bolts.

In recent years, scholars have proposed some novel theoretical models for fatigue life prediction [30,31]. Although theoretical models can yield accurate solutions provided that accurate models can be built, it is difficult for structures under complex loads. There are several dynamic stress simulation methods combining experimental and numerical analysis suitable for bolt dynamic stress simulation. The quasi-static superposition method and the modal superposition method are two common methods in structural dynamic stress simulation. The quasi-static superposition method [32,33] is widely used during the design stage, but this method is not suitable for resonance stress simulation. As for the modal superposition method, the main idea is to obtain the load-stress TF (Transfer Function) of structure, and the external loads are often processed as PSD (Power Spectral Density) in prior research cases [32,34,35,36]; then, the stress PSD can be calculated by multiplying the load PSD by the load-stress TF. Although the modal superposition method is suitable for dynamic stress simulation under resonance conditions, the phase information of the load PSD and TF is ignored. In fact, the effect of the phase between multiaxial excitations on the stress response is so significant that it cannot be ignored; fortunately, the modal superposition method based on FFT of multiaxial excitations and load-stress FRFs takes into account the phase information of multiaxial excitations and was chosen as the bolt dynamic stress simulation method [22].

To simulate the dynamic stress of the bolt in resonance states of structures, the dynamic characteristics of the connected structure must be experimentally analyzed; therefore, processing of the experimental signal is also critical. Because the precision of the fast Fourier transform (FFT) and short-time Fourier transform (STFT) depend on the kind and size of window function, frequency leakage is easy to occur. To avoid leakage, the Hilbert-Huang transform (HHT) was introduced based on an algorithm called empirical mode decomposition (EMD) [37], which is widely used today to recursively decompose a signal into different modes of unknown but separate spectral bands, especially appropriate for nonstationary and nonlinear signals. However, EMD is known for limitations like modal aliasing and edge effects; thus, the algorithm was further improved with the EEMD [38] and CEEMD [39] methods. These methods removed the noise by adding a uniform white noise to the original signal; although mode aliasing can be avoided, amplitude of white noise has great influence on the calculation accuracy. With further improvement, the variational mode decomposition (VMD) method was proposed and proved to be effective in avoiding the shortcomings of the previous methods [40].

1.2. Research Background and Method

As a critical component in high-speed trains, the axle box is mainly composed of the axle box body and axle box front cover; the front cover is connected to the axle box body by six M8 × 40 steel bolts with strength grade of 8.8, as shown in Figure 1. After the high-speed train ran for a period of time on a certain railway line, the connecting bolts would loosen even fracture from time to time, but these failed bolts were theoretically far from reaching the design service life. To solve the problem, the failed bolt was replaced by a new bolt and some prevention methods were used to improve the service life of the connecting bolt, such as adopting a spring washer, replacing the lubricating oil with grease or increasing the preload; these methods are called common prevention methods in this study. However, the new connecting bolts still failed quickly, which obviously greatly threatened the operation safety of the high-speed train; only when the failure mechanism of the connecting bolt is found can the intractable problem be solved.

It was proved that the axle box front cover was in subharmonic resonance of order 1/2 under excitation of 20th-order wheel polygonization, and the resonance caused greater transverse stress of the bolted joint [41]. In this study, an effective method based on numerical and experimental analysis was proposed for multiaxial stress calculation of bolt threads. Firstly, effectiveness of the common prevention methods under structural resonance was assessed based on the nonlinear vibration response of the front cover. Next, feasibility of the method was proved in terms of the axial and transverse stiffness of bolted joints because of the existence of nonlinear characteristics of bolted joints and the first application of this linear stress simulation method on bolts. Then, multiaxial dynamic stresses at the root of the first engaged bolt thread were simulated, the equivalent stress was obtained by the octahedral shear stress criterion and adopted to assess the bolt fatigue life. Finally, a suggestion to improve the fatigue life of bolts was put forward.

2. Assessment for Effectiveness of Common Prevention Method

In this section, firstly, dynamic behavior of the front cover was analyzed according to results of the bench test under 20th-order wheel polygonization. Then, a nonlinear model of the front cover was built by a single degree of freedom (SDOF) modeling strategy, and the feasibility was verified. Finally, effectiveness of common prevention methods was assessed based on the nonlinear vibration response of the front cover.

2.1. Dynamic Characteristics of the Bolted Front Cover under Condition of Wheel Polygonization

Dynamic characteristics of the front cover were experimentally studied on a test bench, as shown in Figure 2. The roller was out-of-round and a 13th-order polygonal wear was formed along the circumference of the tread, the polygonal roller not only simulated the rail but also transformed the roller polygonization into the wheel polygonization by high-speed rotating under the drive of the motor. The weight of the car body was applied by the hydraulic cylinder; thus, the high-speed advance of the polygonal wheel on the rail was simulated successfully at the test bench [41].

In addition, the polygonal order was proportional to the wheel radius, and radiuses of the roller and wheel were 0.3 m and 0.46 m, respectively; thus, the 13th-order polygonal roller meant the wheel had a 20th-order polygonization. Acceleration measuring points were set at the front cover and axle box body, as shown in Figure 1. Then the sweep frequency test was performed with gradual increase of the roller’s rotational speed, and the maximum linear speed was 320 km/h.

2.1.1. Subharmonic Resonance

Figure 3 shows the vertical vibration characteristics of the front cover when the radial deviation of the polygonal roller was 0.075 mm. As shown in Figure 3a, the vibration amplitude was greatly enlarged at region 1 and 2; the corresponding excitation frequency of 20th-order wheel polygonization was also labelled as 1 and 2 in Figure 3b. The relationship between the excitation frequency f and linear speed of wheel v (m/s) is described in Equation (1):

(1)$f=\frac{vn}{2\pi r}$

where n = 20 is the polygonal order of wheel; r = 0.46 m is the wheel radius. The linear speed of the wheel was 288 km/h at 1 and 2, so the excitation frequency of wheel polygon was 550 Hz according to Equation (1). However, the dominant frequency was around 275 Hz at 1 and 2, which was the natural frequency of the bolted front cover and equal to half the excitation frequency; this vibration characteristic indicated that the front cover was in subharmonic resonance of order 1/2, and 288 km/h was a critical speed.

2.1.2. Superharmonic Resonance

Spectral analysis results of the signal of Figure 3c are represented in Figure 4a; obviously, superharmonic resonance of order 2 for the front cover also occurred. Figure 4b represents the spectral analysis results of vertical vibration of the front cover under wheel radial deviation of 0.1 mm and excitation frequency of 280 Hz. According to the main frequencies in Figure 4b, the front cover was not only in primary resonance but also in superharmonic resonance of order 2 and 3.

2.1.3. Jump Phenomenon

Other than subharmonic resonance and superharmonic resonance, the jump phenomenon is also a unique characteristic of nonlinear vibration, in which the response of structural vibration has a sudden change under specific external excitation. The essence is due to stiffness softening or stiffness hardening of a system, which are also called soft spring or hard spring. Figure 5 is a simple illustration of resonant response curves of linear and nonlinear systems; f_e, f₀ and a refer to excitation frequency, natural frequency and response amplitude respectively. The linear system is resonant when f_e ≈ f₀, but the soft spring nonlinearity bends the frequency-response curve to the left, while the hard spring nonlinearity bends the curve to the right. Consequently, the vibration amplitude will jump at certain f_e values (points 1, 2, 3, 4) below and above the natural frequency, respectively.

Figure 3c,d represent the time-domain signals at region 1 and 2 in Figure 3a. Figure 3c shows the jump phenomenon in the acceleration process;. the direction of the arrows represent an increase in speed, the response amplitude jumped from A to B when the wheel speed increased to v₁, which was slightly less than 288 km/h. Speed was constant between B and E, then the wheel speed continued to increase, but the vertical amplitude of the front cover gradually decreased until the maximum test speed of 320 km/h was reached. In contrast, the jump phenomenon in the deceleration process is shown in Figure 3d; while the direction of the arrows represent an decrease in speed, vertical vibration of the front cover gradually increased to maxima at F when the speed was reduced to v₂ (v₂ was slightly smaller than v₁). Speed was constant between F and C, followed by the jump phenomenon from C to D when the speed continued to decrease from v₂. In short, the jump phenomenon of vibration amplitude occurred below the natural frequency, which corresponds to the soft spring curve in Figure 5.

2.2. Nonlinear Modeling of the Bolted Front Cover

The front cover was connected by six bolts evenly distributed along the circumference, which can be simplified to a spring-mass-damper system of SDOF. According to Newton’s second law, dynamics of the front cover can be represented by the second-order ordinary differential equation:

(2)$m\ddot{x}=C(\dot{x},x)\dot{x}+K(\dot{x},x)x+E(t)$

which is a general model for dynamic behavior analysis and nonlinear parameter identification of the bolted joints [42,43], where m = 3 kg is mass of the front cover; the variable $\ddot{x}$ represents acceleration; E(t) = Fcos2πf_et is external force from the periodic excitation of wheel polygonization; f_e is excitation frequency in Hz; $C(\dot{x},x)\dot{x}$ and $K(\dot{x},x)x$ represent the damping force and restoring force, respectively.

According to the analysis results of Section 2.1, vertical vibration of the front cover had the nonlinear characteristics of subharmonic resonances and jump phenomena. It has been proven that subharmonic resonance of order 1/2 and superharmonic resonance of order 2 could be caused by quadratic nonlinear stiffness, while the cubic stiffness nonlinearity took responsibility for occurrence of the superharmonic resonance of order 3 [44], therefore, the restoring force can be expressed as the form of

(3)$K(\dot{x},x)x=-k_{1}x-k_2{x}^2-k_3{x}^3$

The jump phenomenon of a hard spring system can be described by combining the quadratic and cubic nonlinear stiffness as well as linear damping [44]. However, it was found that the jump phenomenon of a soft spring system could occur under the same conditions, which showed that the system with soft spring nonlinearity was more complex than that with hard spring nonlinearity. Therefore, the effect of nonlinear damping was considered in this study, and the damping force had the form of

(4)$C(\dot{x},x)\dot{x}=c_1\dot{x}+c_2{\dot{x}}^2$

The governing equation of vertical vibration of the front cover was obtained by substituting Equations (3) and (4) into Equation (2) and eliminating m:

(5)$\ddot{x}+2\mu \dot{x}+4{\pi }^2{f}_0{}^{2}x+{\alpha }_2{x}^2+{\alpha }_3{x}^3+{\alpha }_4{\dot{x}}^2=k\cos 2\pi f_{e}t$

where f₀ represents the natural frequency of the derived system in Hz, and the derived system is the corresponding linear system when ${\alpha }_2$, ${\alpha }_3$ and ${\alpha }_4$ are zero; k = F/m is the amplitude of excitation from the axle box body; ${\alpha }_2$ and ${\alpha }_3$ refer to the coefficients of quadratic and cubic stiffness, respectively; $\mu$ and ${\alpha }_4$ are coefficients of linear damping and quadratic nonlinear damping, respectively; $\mu$ is proportional to the damping ratio ζ:

(6)$\mu =\frac{c_1}{2m}=\frac{\zeta \omega }{2m}$

ζ usually ranges from 1% to 10% in engineering, and the damping ratio must be relatively small because of the occurrence of nonlinear vibration; let ζ = 1.5%, so $\mu$ = 4.3.

The method of multiple scales used to seek a second approximation [45] for the solution of Equation (5) and yields the frequency-response Equation (7) for the non-trivial solution:

(7)${(\epsilon {\Phi }_1{a}^2-\pi \sigma f_0+\epsilon {\Phi }_2)}^2-{\Phi }_3^2{k}^2+{(2\pi f_0\mu )}^2=0$

(8)${\Phi }_1=\frac{3{\alpha }_3}{8}-\frac{5{\alpha }_2^2}{48{\pi }^2{f}_0{}^2}-\frac{5{\alpha }_2{\alpha }_4}{12}-\frac{{\alpha }_4^2}{6},$

(9)$\begin{array}{c}{\Phi }_2=\frac{k^2}{16{\pi }^4{\left(f_0{}^2-f_e{}^2\right)}^2}\left[3{\alpha }_3+\frac{\left({\alpha }_2+4{\pi }^2{f}_e{}^2{\alpha }_4+4{\pi }^2{f}_0{f}_e{\alpha }_4\right)\left({\alpha }_2-4{\pi }^2{f}_0{f}_e{\alpha }_4\right)}{2{\pi }^2{f}_e\left(f_e+2f_0\right)}\right.\\ -\frac{{\alpha }_2\left({\alpha }_2+4{\pi }^2{f}_e{}^2{\alpha }_4\right)}{2{\pi }^2{f}_0{}^2}\left.-\frac{\left({\alpha }_2+4{\pi }^2{f}_0{f}_e{\alpha }_4\right)\left({\alpha }_2+12{\pi }^2{f}_0{f}_e{\alpha }_4\right)}{8{\pi }^2{f}_0{}^2}\right]-\frac{{\mu }^2}{2}\end{array},$

(10)$\begin{array}{ll}{\Phi }_3=& \frac{1}{16{\pi }^4{\left(f_e{}^2-f_0{}^2\right)}^2}\left\{{\left({\alpha }_2+4\pi f_0{f}_e{\alpha }_4\right)}^2\left(1-\frac{{\epsilon }^2\sigma }{4\pi f_0}\right)+\right.\\ & \left.{(2\pi {\epsilon }^2\mu f_e)}^2{\left[\frac{{\alpha }_2+4\pi f_0{f}_e{\alpha }_4}{2{\pi }^2\left(f_e{}^2-f_0{}^2\right)}-{\alpha }_4\right]}^2\right\}\end{array}$

where a is the amplitude of frequency response; the parameter $\epsilon \ll 1$; σ is the detuning parameter, which quantitatively describes the nearness of f_e to f₀.

In this study, ε = 0.2, f₀ = 275 Hz, and $f_e=2f_0+\frac{\epsilon \sigma }{2\pi }$. As shown in Figure 3b, the jump phenomena occurred at f_e ≈ 545 Hz in the acceleration process and f_e ≈ 540 Hz in the deceleration process; therefore, the difference of σ was about 25 between the two jump points. It is difficult to determine the coefficient of nonlinear term in both theoretical calculation and experimental analysis; a set of values of ${\alpha }_2$, ${\alpha }_3$ and ${\alpha }_4$ is given to make the frequency response (Figure 6) and acceleration response (Figure 7) as close to reality as possible. As shown in Figure 6, the subharmonic resonance amplitude a jumped as the change of f_e. The solid arrows represent the increase process of f_e; a jumped upward from A to B when f_e increased to about 545 Hz, then a gradually decreased as f_e continued to increase. In the decrease process of f_e, which is marked with the dotted arrows, a gradually increased as f_e decreased, and jumped downward from C to D until f_e decreased to about 540 Hz. It can be seen that the two jump phenomena correspond to the jump phenomena of the acceleration and deceleration process in the bench test (Figure 3c,d), respectively.

2.3. Effect of Multiple Parameters on Structural Subharmonic Resonance

In this section, effectiveness of the common prevention method for connecting bolt failure under structural subharmonic resonance was studied by investigating the effect of multiple parameters of Equation (5) on the vertical vibration response of the front cover system.

2.3.1. Parameter Setting

Values of the multiple parameters are shown in Table 1. Curve O is used for comparison, which stands for original values of the parameters that are also the values for Figure 6. In addition, it was found that when the value of ${\alpha }_2$, ${\alpha }_3$ and ${\alpha }_4$ were in the range of (−7.3, −4.3), (−42, 9.9) and (2.76, 4.49), respectively, the jump phenomena would occur.

2.3.2. Comparison between the Test Signal and Numerical Signal

Figure 7a is the locally amplified signal of Figure 3c; there are several frequency components for vertical vibration of the front cover. The test acceleration signal was decomposed by VMD to obtain the main frequencies and reduce noise; the VMD function decomposes a signal into a small number K (K ≥ 2) of narrowband intrinsic mode functions (IMFs) and can effectively avoid modal aliasing and edge effects. To avoid over-decomposition and under-decomposition, the appropriate number K can be predicted by the evaluation index of the mean absolute percentage error (MAPE) [46], which represents the ratio of the residual after VMD decomposition to the original signal. In this paper, the minimum K value with MAPE value less than 1% was selected as the optimal number of mode decomposition. Through verification, as shown in Figure 7b, it was found that K = 5 (MAPE = 0.5657%) was appropriate. As shown in Figure 8, the test signal (Figure 7a) was decomposed to five IMFs; each IMF has a central frequency, and IMF 3 and IMF 4 represent the 550 Hz vibration and 275 Hz vibration, respectively; the former is the excitation frequency of wheel polygonization, while the latter is the subharmonic resonance frequency.

IMF 5 can be regarded as a noise signal due to its irregular vibration and small amplitude. For comparison, the test signal of subharmonic resonance of order 1/2 under excitation of 550 Hz was constructed by adding IMF 3 and IMF 4, as shown in Figure 9. The numerical signal was obtained by substituting the parameters of curve O in Table 1 into Equation (5). Both the experimental and numerical vibrations were steady-state motions and have very similar characteristics.

As shown in Figure 10, the time-frequency spectrums of the two signals were obtained by HHT. The energy distribution showed that both the signals contained frequencies of 275 Hz and 550 Hz, and 275 Hz was the main frequency; in other words, subharmonic resonance of order 1/2 was the most outstanding vibration characteristic, which is also a key point throughout this study. Consequently, it was further proved that the SDOF modeling method is reasonable and feasible to simulate the vertical vibration of the front cover.

2.3.3. Qualitative Analysis for Effect of Multiple Parameters

In this subsection, the effect of multiple parameters on subharmonic resonance response of the front cover was analyzed. As shown in Figure 11a, it was effective to reduce the vibration amplitude by changing the excitation amplitude and the initial condition. The effect of the linear damping coefficient μ on vibration response was closely related to initial conditions; under the nonzero initial condition, the decrease of μ increased the subharmonic resonance response. However, subharmonic resonance was not excited in the zero initial condition; even if μ was further reduced, the vibration amplitude was much smaller than the response in subharmonic resonance. As shown in Figure 11b, subharmonic resonance was more sensitive to ${\alpha }_3$ and ${\alpha }_4$ than to ${\alpha }_2$. In addition, there is no doubt that f₀ and f_e have great influence on the subharmonic resonance, but as long as $f_e\approx 2f_0$, changing the natural frequency of the derived system will only change the cycle period without affecting the occurrence of subharmonic resonance.

2.3.4. Evaluation of the Common Prevention Methods

Essentially, the common prevention methods of adopting spring washers, replacing the lubricating oil with grease, and increasing the preload all increased the friction between interfaces, which is effective for bolt connection under static loads because the friction is more difficult to overcome. However, according to the above analysis results, uncertainty of the effect of damping on a nonlinear system means increasing the preload may not achieve a satisfactory result. Although adequate damping could prevent subharmonic resonance, it is impossible to increase the damping of bolted joints so much by the common prevention methods. In addition, the linear stiffness of the connecting bolt, k₁, can be expressed as

(11)$k_1=\frac{EA}{l}$

where E is Young’s modulus of the bolt material, and A and l are the sectional area and length of bolt, respectively. When the preload increases, the decrease of A and increase of l lead to the decrease of bolt connection stiffness k₁. Meanwhile, the natural frequency of the derived system $\omega$, which is proportional to k₁, also decreases. In this case, the change of f₀ caused by preload is so small due to the small variation of A and l that it is difficult to affect the structural vibration significantly. In other words, the common prevention methods were not effective enough for improving the service life of connection bolts in the front cover because dynamic behavior of the nonlinear system was hard to predict and control. In the following sections, the fatigue life of bolts under structural subharmonic resonance are analyzed and solutions proposed to improve the service life of the connecting bolt.

3. Feasibility Analysis for the Stress Simulation Method

Nonlinear characteristics of structure were ignored automatically at the modal-related analysis stage of dynamic stress simulation. Therefore, feasibility of the calculation method was verified in this section. Firstly, axial load measuring bolts were made and calibrated to measure the axial dynamic stress of bolts; moreover, multiaxial acceleration excitations as well as the operational modes of the front cover were also measured from the vibration test. Then, the FE model of the bolted front cover was established, and the mode-based steady-state dynamics analysis was performed to obtain the FRFs (Frequency Response Functions). Finally, the vertical acceleration of the front cover and axial stress of bolts were simulated based on the multiaxial excitations and FRFs.

3.1. Theoretical Background of the Stress Simulation Method

In this section, the stress simulation method was studied to calculate accurate structure dynamic stress under resonance and multi-load conditions. Based on the structure dynamic theory and previous research results [22,47], the dynamic stress time history calculation method was adopted in this study, in which the phase information between multi loadings was not ignored. Dynamics of multi-degree-of-freedom systems subjected to time-varying excitation can be governed by the following differential equations:

(12)$M\ddot{X}(t)+C\dot{X}(t)+KX(t)=f_{\mathrm{e}}(t)$

where M, C, and K are mass, damping and stiffness matrices, respectively; $X(t)$ and $f_{\mathrm{e}}(t)$ are displacement and excitation vectors, respectively. Assuming that there are totally R modes, the load-stress FRFs can be obtained by Equation (12) [47]:

(13)$H(\Omega )=E·{\displaystyle \sum _{r=1}^R\frac{\left\{{\psi }_r^{\epsilon }\right\}{\left\{{\phi }_r\right\}}^T}{k_r-{\Omega }^2{m}_r+j\Omega c_r}}$

where $m_r$, $c_r$, $k_r$ are the rth modal mass, modal damping, and modal stiffness, respectively; $\Omega$ is the excitation frequency; $\left\{{\phi }_r\right\}$ is the real value mode vector; $\left\{{\psi }_r^{\epsilon }\right\}$ is the rth strain mode corresponding to $\left\{{\phi }_r\right\}$; $j=\sqrt{-1}$; E is elastic module.

Due to the front cover being fixed to the axle box body, so accelerations at the end of the axle box body were regarded as the input load. For the acceleration-stress FRF, the phase information needs to be retained, as shown in Equation (15):

(14)$\begin{array}{ll}H(\Omega )& =E·\left\{{\displaystyle \sum _{r=1}^R\frac{\left\{{\psi }_r^{\epsilon }\right\}{\left\{{\phi }_r\right\}}^T}{{(k_r-{\Omega }^2{m}_r)}^2+{(\Omega c_r)}^2}[(k_r-{\Omega }^2{m}_r)+j(-\Omega c_r)]}\right\}\\ & =E·\left\{{\displaystyle \sum _{r=1}^R\frac{(1-{\lambda }_r^2)}{k_{er}[{(1-{\lambda }_r^2)}^2+{(2{\zeta }_r{\lambda }_r)}^2]}+j·{\displaystyle \sum _{r=1}^R\frac{-2{\zeta }_r{\lambda }_r}{k_{er}[{(1-{\lambda }_r^2)}^2+{(2{\zeta }_r{\lambda }_r)}^2]}}}\right\}\\ & =H^R(\Omega )+j·H^I(\Omega )\end{array}$

where $k_{er}=k_r/\left\{\right\{{\psi }_r^{\epsilon }\left\}{\left\{{\phi }_r\right\}}^T\right\}$, ${\lambda }_r=\Omega /{\omega }_r$, ${\omega }_r{}^2=k_r/m_r$, ${\zeta }_r=c_r/\sqrt{4m_r{k}_r}$.

In addition, assuming the length of a discrete acceleration signal f(n) is N, the FFT of f(n) is

(15)$\begin{array}{ll}F(\Omega )& ={\displaystyle \sum _{n=0}^{N-1}f(n)e^{-j\Omega n}}\\ & ={\displaystyle \sum _{n=0}^{N-1}f(n)\cos n\Omega -j}·{\displaystyle \sum _{n=0}^{N-1}f(n)\sin n\Omega }\\ & =F^R(\Omega )+j·F^I(\Omega )\end{array}$

In the case of multi-load input, that is, when the number of external loadings is more than one, and assuming that the number is k_e, all the loadings should be considered to stimulate the stress response; thus, the frequency domain of concerned point stress under the multi-load input condition can be calculated by Equation (16) below:

(16)$\begin{array}{ll}\sigma (\Omega )& ={\displaystyle \sum _{i=1}^{k_{\mathrm{e}}}[F_i(\Omega )·}H(\Omega )]\\ & ={\displaystyle \sum _{i=1}^{k_{\mathrm{e}}}\{[F_i{}^R(\Omega )H^R(\Omega )-F_i{}^I(\Omega )H^I(\Omega )]+j·[F_i{}^R(\Omega )H^I(\Omega )-F_i{}^I(\Omega )H^R(\Omega )]\}}\end{array}$

Then, the time history stress at the concerned point can be calculated by IFFT:

(17)$\sigma (N)=\frac{1}{N}{\displaystyle \sum _{n=0}^{N-1}[\sigma (\Omega )e^{j\Omega n}]}$

$\sigma (N)$ is the simulated dynamic stress time history, which is used for fatigue life prediction.

3.2. Acquisition of Vibration Signals

3.2.1. Fabrication and Calibration of Axial Load Measuring Bolt

A nonlinear connection such as a bolted joint is usually simplified as a linear system because the accurate nonlinear model is difficult to be established, but the premise is that the accuracy of calculation results can be guaranteed. In this study, damping of the bolted joint was assumed to be linear; therefore, feasibility of the linear stress simulation method was assessed in terms of the axial and transverse stiffness of the bolted joint.

The axial stiffness could be assessed by certifying the bolt axial stress. In order to measure the axial bolt stress in the vibration test, axial stress measuring bolts were made, as shown in Figure 12. Firstly, a hole with diameter of 2.2 mm and an appropriate depth was drilled in the center of the bolt. Then, an LB11 strain gauge was attached to inside of the bolt by a hot-curing epoxy resin adhesive. LB11 is a cylindrical strain gauge to measure strain, force, and vibration in screws or bolts [48].

Drilling reduced the effective cross section and durability of the bolt, which resulted in that the stress-strain relationship could not be depicted with the classical material mechanics formula. Hence, calibration for the drilled bolts was conducted by a strain gauge load cell for tension and compression. Calibration for each bolt was repeated three times, and the results had great consistency. The relationship between the micro strain ε_m (μm/m) and the bolt axial stress S_a (MPa) is shown in Equation (18).

(18)$S_a=\frac{7.9704{\epsilon }_m-2.813}{A_e}$

where A_e is the effective sectional area of the bolt shank, which should minus the area of the drilled hole.

3.2.2. Measurement of Multiaxial Excitations and Bolt Axial Load Variation

As shown in Figure 1, triaxial accelerometers were set at the end of the front cover and the axle box body where near the front cover, two axial load measuring bolts were assembled symmetrically at the front cover and tightened with the torque of 20 Nm. Only vibrations of the axle box body were considered as multiaxial external loads because almost all the external excitations of the front cover came from the axle box body; the multiload can be divided into the longitudinal acceleration, the lateral acceleration, and the vertical acceleration.

The multiaxial excitations of the front cover are shown in Figure 13: the left is the time domain signal at the running speed of 288 km/h, while the right is the FFT result of the left. The largest amplitude of vertical vibration was 40 g, followed by longitudinal vibration, and the smallest was lateral vibration. The dominant frequencies of the three accelerations were all multiples of 27.49 Hz, which is the theoretical rotational frequency of the wheel that is equal to wheel perimeter divided by linear speed of 288 km/h and is also the frequency when n = 1 in Equation (1).

As shown in Figure 14, the dynamic axial stress of the bolt was obtained in accordance with the micro-strain measured by LB11 strain gauge and Equation (18). Obviously, the amplitude of dynamic axial stress also increases in the subharmonic resonance state and changes between ±2 MPa. It is noticeable that the micro-strain measurement started after bolt tightening, so axial stress caused by preloading is not included in Figure 14.

3.2.3. Modal Test of the Front Cover

In this section, the vibration mode of the front cover was obtained by modal test and used to verify the correctness of the subsequent numerical modeling. OMA (Operational Modal Analysis) rather than EMA (Experimental Modal Analysis) was chosen due to the practical advantages [41,49,50,51], which is preferable for in-operation structures because the output-only data is needed and the test structure remains in its normal in-operation condition during the test. More importantly, OMA is quite suitable for structures with nonlinear characteristics.

As shown in Figure 15a, 19 measuring points of acceleration were set on the external surface of the front cover (1–13) and close to the connection bolts (14–19). Then, the frequency sweep test was carried out as in Section 2.1 to obtain the vibration signals in three directions. Finally, modal analysis was carried out in the vibration control system of Simcenter Testlab according to the acceleration signals. As shown in Figure 15b, the solid line model is the undeformed model, while the dotted line model represents the deformed model, and the first-order modal frequency and modal damping were 274.36 Hz and 1.88%, respectively. There are three arrows with different colors at each measuring point; the green, blue, and red arrows represent the vibrations in the x, y, and z directions, respectively. The size of the arrow represents the vibration intensity; obviously, the front cover mainly vibrated in the z direction, that is, perpendicular to the axis of connection bolt, which means that the bolts were mainly subjected to transverse loads.

3.3. Finite Element Analysis

3.3.1. Finite Element Modeling

As shown in Figure 16a, the three-dimensional FE model of the front cover was established in the CAE software of HyperMesh 2019. The numbers 1-6 refer to the six connecting bolts, the external threads of the bolts and internal threads of the axle box body were constructed with the thread pitch of 1.25 mm. The front cover was divided by tetrahedron elements with type of solid 185, while the bolts, washers and cylinders were divided by the mixture of pentahedrons and hexahedrons. Since the axle box body was assumed to be rigid, only a small region around the tapped hole was modeled, and the external nodes were constrained. In addition, contact elements were overlaid on contact regions between the bolt head and washer, the washer and front cover, the front cover and axle box body, as well as between the internal and external threads. The connecting bolts were unlubricated galvanized bolts, which fit into to the friction coefficient class C of Table A5 in VDI 2230 [52], so the friction coefficients range from 0.14 to 0.24. In this study, the friction coefficients were set to 0.18.

Preload of 11,000 N was applied at shank for each bolt in ABAQUS 2021, followed by a general static analysis; convergence of the contact was calculated by the Newton-Raphson method. Stress distribution on the bolt after preloading is shown in Figure 16b; the maximum stress was 408.5 MPa, which is lower than the material yield limit (640 MPa) of grade 8.8 bolt. Moreover, the maximum stress occurred at the root of the first engaged thread, which is reasonable due to the high notch effect at the first load-bearing thread turn of the bolt [52]. The preload calculation results show that FE modeling of the front cover was reasonable.

3.3.2. Modal Analysis

Modal analysis of the bolted front cover was carried out after the preload calculation; the modes were extracted by the block Lanczos method. The first four order modal frequencies are listed in Table 2, and the first two were around 275 Hz, which further proves that subharmonic resonance of order 1/2 did occur when the excitation frequency was 550 Hz.

As shown in Figure 17, the first two modes of the connection bolts were the first-order bending modes; in other words, the front cover vibrated in the z direction and the bolts were subjected to transverse loads. The mode frequencies and mode shapes of the first two modes were in good agreement with the OMA results, so validity of the FE model was further proved.

3.3.3. FRF Calculation

Mode-based steady-state dynamics analysis was conducted to calculate the FRFs, which included the real and imaginary parts. A white noise load was regarded as the base excitation input, the frequency band of white noise was from 200 Hz to 600 Hz, which contains the frequencies of interest in this study. In addition to the accurate modal frequency and modal shape, accurate modal damping was also needed for accurate calculation of FRFs. Therefore, the damping ratio 1.88% of OMA was adopted for the first two calculation modes.

Feasibility of the linear stress simulation method in axial direction of the bolted joint can be verified easily because bolt axial stress can be measured directly, but transverse stress of bolt was difficult to be measured; therefore, the method was verified in the transverse direction of the bolted joint by verifying the vertical vibration of the front cover. The vertical acceleration FRFs of the front cover are shown in Figure 18a–c; the main frequency was about 275 Hz under longitudinal and vertical excitation, which is close to the first-order bending modal frequency of bolts, while the main frequency was about 395 Hz under the lateral excitation. Simultaneously, the axial stress FRFs of the bolt shank were also calculated to simulate the axle stress of bolts, as shown in Figure 18d–f.

3.4. Acceleration and Stress Simulation of the Bolted Front Cover

In this section, dynamic data in the frequency domain are calculated by the frequency domain signals of the multiaxial excitations in Figure 13 and FRFs in Figure 18. Then, dynamic time history can be acquired by Inverse FFT (IFFT) of the frequency domain data. A comparison between the simulation results and test results is shown in Figure 19. As shown in Figure 19a,b, the numerical signals were consistent with the experimental signals. Figure 19c,d show that both the experimental and numerical signals had the dominant frequencies of 274.53 Hz and 548.67 Hz, which corresponded to the first-order bending modal frequency of bolts and the excitation frequency of the 20th-order wheel polygon, respectively.

In summary, the consistency between experimental and numerical results for vertical acceleration and axial stress demonstrates that the linear stress simulation method is able to ensure the calculation accuracy. In other words, this method is feasible for dynamic stress simulation under nonlinear vibration states.

4. Analysis of Multiaxial Fatigue Behavior of Connecting Bolt

After feasibility of the stress simulation method was proved, multiaxial fatigue behavior was analyzed in this section. Transverse and axial resonance stresses at the root of the first engaged thread of connecting bolts were calculated firstly. Then, the multiaxial fatigue analysis criterion was determined. Finally, load change cycles were calculated with the rain-flow counting method, and the bolt fatigue life was predicted.

4.1. Calculation of Dynamic Stress at Bolt Thread

As it is the most dangerous location, resonance stress at the first engaged thread must be evaluated. Figure 20a–c shows the transverse stress FRFs of node 233,773 at the first engaged thread of bolt 2 (Figure 16a). Obviously, the 275 Hz vibration was most sensitive to the vertical excitation. The lateral excitation easily excited the 395 Hz vibration but had little effect on the 275 Hz vibration. In other words, subharmonic resonance of order 1/2 of the front cover was mainly aroused by the vertical vibration.

Figure 20d–f show the axial stress FRFs of node 233,773. It is obvious that the axial stress response at node 233,773 was much larger than the response at the bolt shank (Figure 18d–f). Similarly, 275 Hz vibration was more sensitive to vertical and longitudinal excitations, while 395 Hz vibration was more sensitive to lateral excitation. It means that the axial deformation at node 233,773 under subharmonic resonance was not originated from the axial excitation but from the transverse excitations.

Simulation of the transverse dynamic stress of node 233,773 is shown in Figure 21a. The maximum stress range of the transverse stress was around 20 MPa, which was much larger than the axial stress range at the bolt shank (Figure 19b). As shown in Figure 21b, the tensile stress at the root of the bolt thread could not be ignored due to the bending mode; the maximum range of the axial dynamic stress of node 233,773 was around 45 MPa, which was much larger than the stress range at the bolt shank as well. It is obvious that the axial stress at the bolt thread was more predominant than the transverse stress, which would accelerate fatigue-crack propagation [24,25].

4.2. Analysis of Bolt Fatigue Strength

4.2.1. Criteria of Stresses Combination and Fatigue Life Calculation

There are some stress combination techniques for multiaxial stresses, such as the critical plane approach and the octahedral shear stress criterion; the latter is also often called either the von Mises or the distortion energy criterion. The octahedral stress criterion is work for in-phase or proportional loading, which is suitable for bolt fatigue analysis in this study because the stresses were all caused by the excitation of 20th-order wheel polygonization. There are three main processes to evaluate fatigue life by the octahedral stress criterion: (1) cycle counting for the stress load history, (2) calculation of equivalent stress for multiaxial stresses, and (3) prediction of fatigue life. The fatigue analysis criterion was firstly represented in this subsection.

According to the octahedral stress criterion, the equivalent stress can be expressed as

(19)${\sigma }^{\mathrm{e}}=\frac{1}{\sqrt{2}}\sqrt{{({\sigma }_1-{\sigma }_2)}^2+{({\sigma }_2-{\sigma }_3)}^2+{({\sigma }_3-{\sigma }_1)}^2}$

where σ₁, σ₂, and σ₃ are principal stresses. The amplitudes of the principal stresses, σ_1a, σ_2a, and σ_3a, can then be employed to compute an equivalent stress amplitude ${\sigma }_a{}^{\mathrm{e}}$ by using the similar relationship as Equation (19). In addition, consider the effect of mean stresses; an equivalent mean stress ${\sigma }_m{}^{\mathrm{e}}$ can be calculated from the mean stresses in the three principal directions. Therefore, the cyclic stressing can be specified as

(20)$\left\{\begin{array}{l}{\sigma }_a{}^{\mathrm{e}}=\frac{1}{\sqrt{2}}\sqrt{{({\sigma }_{1a}-{\sigma }_{2a})}^2+{({\sigma }_{2a}-{\sigma }_{3a})}^2+{({\sigma }_{3a}-{\sigma }_{1a})}^2}\\ {\sigma }_m{}^{\mathrm{e}}={\sigma }_{1m}+{\sigma }_{2m}+{\sigma }_{3m}\end{array}\right.$

Equation (20) can be analogized to obtain the equivalent stress amplitude and equivalent mean stress for any convenient coordinate axes:

(21)$\left\{\begin{array}{l}{\sigma }_a{}^{\mathrm{e}}=\frac{1}{\sqrt{2}}\sqrt{{({\sigma }_{xa}-{\sigma }_{ya})}^2+{({\sigma }_{ya}-{\sigma }_{za})}^2+{({\sigma }_{za}-{\sigma }_{xa})}^2+6({\tau }_{xya}^2+{\tau }_{yza}^2+{\tau }_{zxa}^2)}\\ {\sigma }_m{}^{\mathrm{e}}={\sigma }_{xm}+{\sigma }_{ym}+{\sigma }_{zm}\end{array}\right.$

According to the previous analysis, ${\sigma }_y$ and ${\sigma }_z$ at the bolt thread root were obtained in Section 4.1 (Figure 21); ${\sigma }_x$ was neglected for vertical vibration of the front cover, and the shear stresses were also neglected. Thus, equivalent amplitude and mean value of bolt stresses can be expressed as

(22)$\left\{\begin{array}{l}{\sigma }_a{}^{\mathrm{e}}=\frac{1}{\sqrt{2}}\sqrt{{\sigma }_{ya}{}^2+{\sigma }_{za}{}^2+{({\sigma }_{ya}-{\sigma }_{za})}^2}\\ {\sigma }_m{}^{\mathrm{e}}={\sigma }_{ym}+{\sigma }_{zm}\end{array}\right.$

An equivalent stress amplitude of the completely reversed uniaxial stress can be obtained by combining ${\sigma }_a{}^{\mathrm{e}}$ and ${\sigma }_m{}^{\mathrm{e}}$:

(23)${\sigma }_{ar}{}^{\mathrm{e}}=\frac{{\sigma }_a{}^{\mathrm{e}}}{1-\frac{{\sigma }_m{}^{\mathrm{e}}}{{\sigma }_f{}^{\prime }}}$

where ${\sigma }_f{}^{\prime }$ is a constant and approximately equal to the true fracture strength of material. In this study, ${\sigma }_f{}^{\prime }$= 1200 MPa.

The technical committee of the BS standard has determined several S-N curves suitable for fatigue life calculation of bolted joints through fatigue tests under various load conditions and published them in EN1993-1-9:2005 [53]. According to the GL2010 standard [54], for bolts rolled before heat treatment and subjected to tension and bending loads, the S-N curve for category 71 in EN1993-1-9:2005 is applicable. For constant amplitude stress ranges, the S-N curve has the form as follows:

(24)${(\Delta {\sigma }_a)}^m{N}_f={(\Delta {\sigma }_c)}^m\times 2\times {10}^6$

where $\Delta {\sigma }_a$ is cycling stress range, $N_f$ is the number of cycles to fatigue failure, m = 3, and the number $\Delta {\sigma }_c$ = 71 represents the reference value of fatigue strength for the curve of category 71 at 2 million cycles. The S-N curve expressed in stress amplitude (half the stress range) is shown below:

(25)${({\sigma }_{ar}{}^{\mathrm{e}})}^3{N}_f=8.95\times {10}^{10}$

Finally, the fatigue life can be predicted according to Miner’s cumulative damage theory.

4.2.2. Fatigue Life Prediction of Bolt

The number of cycles necessary to break the structure depends on the cycle amplitude and mean value of stress. The rain-flow cycle counting method is an effective method for variable amplitude multiaxial loading which can identify all the stress cycles. According to the rain-flow count, there were totally 2743 different stress cycles and half-cycles for both the transverse stress and axial stress; the number of cycles for each cycle stress was N_i = 1 or 1/2. The stress ranges and mean stresses are shown in Figure 22.

An equivalent stress amplitude ${\sigma }_a{}^{\mathrm{e}}$ and equivalent mean stress ${\sigma }_m{}^{\mathrm{e}}$ were obtained by the 2743 different stress ranges and mean stresses according to Equation (22). Then, the equivalent stress amplitude ${\sigma }_{ar}{}^{\mathrm{e}}$ for fatigue calculation was obtained by Equation (23). According to Equation (25), the failure cycle $N_{fi}$ that corresponded to each N_i was calculated. After getting all the 2743 failure cycles, the accumulation at the root of the first engaged bolt thread caused by transverse and axial cyclic stress can be obtained by Equation (26).

(26)$D={\displaystyle \sum _{i=1}^{2743}\frac{N_i}{N_{fi}}}$

It was calculated that $D=5.18\times {10}^{-6}$, which represents the damage at the thread root caused by the multiaxial excitation signal of 5 s length, so fatigue life of the bolt was $t=5/(5.18\times {10}^{-6}\times 3600)=26.8$ h, which means that fatigue failure of the bolt would occur if only the subharmonic resonance state of the front cover lasted for more than 26.8 h.

5. Method for Improving Fatigue Life of Connecting Bolt

According to the above fatigue analysis results, bolts at the axle box front cover easily failed around the linear speed of 288 km/h after the stable 20th-order polygonal wheel formed. Reprofiling is the most effective way to eliminate the wheel polygonization-induced resonance at present, and 250,000 km is a fixed reprofiling interval for Chinese high-speed trains. Wang et al. [55] reported that vertical vibration of the axle box increased slowly in the early 150,000 km and rapidly in the late 100,000 km, which indicated rapid growth of the polygonization amplitude. Moreover, 280 km/h to 300 km/h is a common running speed range for high-speed trains in China. Therefore, there is a high probability of fatigue failure for connecting bolts of the front cover during a reprofiling interval.

Subharmonic resonance of the front cover must be avoided to prevent the multiaxial vibration fatigue of bolts. There are two methods to improve fatigue life of bolts. One is to optimize the structure of the front cover so that the first-order modal frequency is greater than 290 Hz, which is half the excitation frequency of the 20th-order wheel polygon at the running speed of 300 km/h. Nevertheless, modifying the structure of a widely used component is costly. The other method is to prevent the formation of 20th-order wheel polygonization. In order to explore the fatigue behavior of bolts with no wheel polygonization, a sweep frequency test under the condition of wheel tread defect of wheel-flat but no polygonization was carried out. The reason why the wheel with defects was used instead of the new wheel after reprofiling was that the wheel would not be in a perfect state in most of the running time, so a working condition that was inferior to the new reprofiled wheel but better than the polygonal wheel was used. The same frequency sweep test was carried out as previously. The maximum vibrational amplitude of the multiaxial acceleration loadings from the axle box body was around 20 g at the linear speed range of 280–300 km/h, and the front cover was not in resonance; thus, the stress caused by bending modes of the bolt was ignored under the circumstances, and only the effect of transverse stress on bolt fatigue life was considered. The results of transverse stress simulation at the first engaged thread root and local amplification are shown in Figure 23. Obviously, there was an aperiodic vibration with smaller transverse stress amplitude in comparison with the resonance stress in Figure 21a.

The fatigue life of bolts was calculated by the same method as in Section 4. Accumulation caused by transverse vibration at the root of the first engaged thread was $D_{\tau }=1.95\times {10}^{-9}$, so the fatigue life of the bolt was more than 700 thousand hours, which means the probability of bolt fatigue failure decreases rapidly as long as the front cover is not in resonance, even though the wheel tread has defects. Hence, shortening the reprofiling interval to reduce the polygonal wear could effectively improve the fatigue life of bolts.

6. Conclusions

A linear method based on frequency response analysis was first used to calculate the bolt resonance stress in this work. It was proved that the method was feasible and accurate enough for dynamic stress simulation, which provides a practical way for multiaxial fatigue behavior analysis of bolts under complicated running conditions of high-speed trains. Several conclusions were drawn as follows:

• The SODF nonlinear modeling method is reasonable for the bolted front cover. The common prevention methods for bolt failure were theoretically proved ineffective under structural subharmonic resonance of order 1/2.

• The results of bolt stress measurement show that the dynamic behaviors of bolt under a nonlinear vibration state should be assessed due to the increasing of bolt stress.

• Although the vibration amplitude of the front cover in the direction of the bolt axis was small, the axial resonance stress at the bolt thread’s root could not be neglected due to the first-order bending modes of bolts.

• Feasibility of the linear stress simulation method was proved in terms of the transverse stiffness and axial stiffness of the bolted joint. It indicated that the method was accurate enough to simulate the subharmonic resonance stress of bolts even if the nonlinearity of the bolted joint was ignored.

• Resonance stresses at the root of the first engaged bolt thread were much larger than the resonance stress at the bolt shank, and the axial resonance stress was more predominant at the bolt thread than the transverse resonance stress.

• Since the multiaxial stresses were caused by the homologous excitations, the octahedral shear stress criterion was suitable for equivalent stress calculation. The fatigue life of the bolt was about 26.8 h, which means that the connecting bolts were prone to multiaxial fatigue failure when the front cover was in subharmonic resonance of order 1/2 for a long time.

• The fatigue life of the bolts is greatly improved when the front cover is not in subharmonic resonance. Consequently, the probability of fatigue failure of bolts could be reduced effectively by shortening the reprofiling interval to reduce the wear of the polygon.

The fatigue life analysis method of bolts adopted in this study has great practical significance for safety monitoring in engineering; it is easy to evaluate the reliability of the bolted joint under various vibration conditions. However, this simulation method neglects the effects of nonlinear factors and is not suitable for strongly nonlinear systems; therefore, a more universal stress simulation method should be explored in subsequent studies. In addition, the octahedral shear stress criterion is not suitable for non-proportional loading, and further study should be advanced for a more universal method, such as the critical plane approach.

Author Contributions

Y.F.: Formal analysis, Investigation, Methodology, Validation, Writing—original draft. F.L.: Supervision, Funding acquisition, Writing—review and editing. K.S.: Investigation, Resources, Writing—review and editing. H.D.: Methodology, Resources, Funding acquisition. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Data Availability Statement

The data are not publicly available due to confidentiality agreements of the Institutional.

Conflicts of Interest

The authors declare that they have no known competing financial interest or personal relationships that could have appeared to influence the work reported in this paper.

Nomenclature

EMA	Experimental modal analysis	$H^I(\Omega )$	the imaginary part of H(Ω)
EMD	empirical mode decomposition	HHT	Hilbert-Huang transform
f(n)	the discrete acceleration excitation signal	IMF	intrinsic mode functions
F(Ω)	FFT result of f(n)	k e	number of external loadings
$F^R(\Omega )$	the real part of F(Ω)	MAPE	mean absolute percentage error
$F^I(\Omega )$	the imaginary part of F(Ω)	N	length of f(n)
FE	finite element	OMA	operational modal analysis
FEA	finite element analysis	SDOF	single degree of freedom
FEM	finite element method	STFT	short-time Fourier transform
FRF	frequency response function	S-N	stress-number of cycles
FFT	fast Fourier transform	VMD	variational mode decomposition
H(Ω)	the acceleration-stress FRF	$\sigma (\Omega )$	frequency domain signal of simulation stress
$H^R(\Omega )$	the real part of H(Ω)	$\sigma ({\rm N})$	time domain signal of simulation stress

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Structure of the axle box of high-speed train.

Figure 2. Vibration test bench.

Figure 3. Vertical vibration characteristics of the front cover in bench test. (a) Time-domain signal of acceleration; (b) time-frequency analysis result; (c) jump phenomenon in acceleration process (amplification of the signal in box 1 of Figure 3a); (d) jump phenomenon in deceleration process (amplification of the signal in box 2 of Figure 3a).

Figure 4. Spectral analysis results of vertical vibration signals of the front cover. Excitation frequency of 20th-order polygonization was (a) 550 Hz under wheel radial deviation of 0.075 mm; (b) 280 Hz under wheel radial deviation of 0.1 mm.

Figure 5. Frequency response curve of linear and nonlinear system.

Figure 6. Frequency-response curve of the governing Equation (5).

Figure 7. Vertical acceleration signals of the front cover. (a) Test signal; (b) the optimal number of mode decomposition.

Figure 8. VMD results of the test signal.

Figure 9. Comparison of the experimental and numerical signal.

Figure 10. HHT of (a) test signal and (b) numerical signal.

Figure 11. Effect of multiple parameters on vertical resonance of the front cover. (a) Effect of excitation amplitude and linear damping; (b) effect of nonlinear terms.

Figure 12. Axial stress measuring bolt.

Figure 13. Multiaxial acceleration excitations in time domain and frequency domain. (a) Longitudinal excitation; (b) lateral excitation; (c) vertical excitation.

Figure 14. Time history of the bolt axial stress.

Figure 15. OMA results of the bolted front cover. (a) The analytical model composed of the acceleration measuring points; (b) the first-order mode.

Figure 16. (a) FE model of the front cover; (b) stress distribution of the preloaded bolt.

Figure 17. Modal analysis results of the front cover and bolts. (a) the first-order mode; (b) the second-order mode.

Figure 18. (a–c) represent vertical acceleration FRFs of the front cover under the excitation of longitudinal load, lateral load, and vertical load, respectively; (d–f) represent axial stress FRFs at bolt shank under the excitation of longitudinal load, lateral load, and vertical load, respectively.

Figure 19. Comparison between the experimental and numerical results. (a,b) are time domain signals of the front cover and bolt, respectively; (c,d) are frequency domain signals of the front cover and bolt, respectively.

Figure 20. (a–c) represent transverse stress FRFs at thread root under the excitation of longitudinal load, lateral load, and vertical load, respectively; (d–f) represent axial stress FRFs at thread root under the excitation of longitudinal load, lateral load, and vertical load, respectively.

Figure 21. Dynamic stress simulation results and local amplification. (a) transverse dynamic stress and (b) axial dynamic stress at thread root.

Figure 22. Rain-flow count results for (a) transverse stress time history and (b) axial stress time history of node 233773.

Figure 23. Transverse stress at the thread root and local amplification when there was no wheel polygonization.

Table 1

Parameter table.

Curve Number	Linear Damping Coefficientμ	Excitation Amplitudek (g)	Nonlinear Coefficients			Initial Conditionx0
			α 2	α 3	α 4
O	4.3	20	−4.6	5	3	(0.03 0.1)
μ	3.5	20	−4.6	5	3	(0.03 0.1)
k	4.3	12	−4.6	5	3	(0.03 0.1)
α 2	4.3	20	−7	5	3	(0.03 0.1)
α 3	4.3	20	−4.6	8	3	(0.03 0.1)
α 4	4.3	20	−4.6	5	4.4	(0.03 0.1)
x 0	2	20	−4.6	5	3	(0 0)

Table 2

Modal analysis results.

Mode	1	2	3	4
Frequency (Hz)	274.52	276.51	394.80	614.65

References

1. Hashimura, S. A study of self-loosening and fatigue failure of bolted joint under transverse vibration—Influences of property class of bolt and plastic region tightening. Key Eng. Mater.; 2007; 353–358, pp. 2037-2040. [DOI: https://dx.doi.org/10.4028/www.scientific.net/KEM.353-358.2037]

2. Junker, G. Criteria for self loosening of fasteners under vibration. Aircr. Eng. Aerosp. Technol.; 1972; 44, pp. 14-16. [DOI: https://dx.doi.org/10.1108/eb034949]

3. Yamamoto, A.; Kasei, S. A solution for self-loosening mechanism of threaded fasteners under transverse vibration. Bull. Jpn. Soc. Precis. Eng.; 1984; 18, pp. 261-266.

4. Yang, X.; Nassar, S. Analytical and Experimental Investigation of Self-Loosening of Preloaded Cap Screw Fasteners. J. Vib. Acoust. Trans. Asme; 2011; 133, 031007. [DOI: https://dx.doi.org/10.1115/1.4003197]

5. Baek, K.-H.; Jeong, N.-T.; Hong, H.-R.; Choi, S.-B.; Lee, E.-S.; Kim, H.-M.; Kwon, J.-W.; Song, S.-Y.; Jang, H.-S.; Lee, H.-Y. et al. Loosening mechanism of threaded fastener for complex structures. J. Mech. Sci. Technol.; 2019; 33, pp. 1689-1702. [DOI: https://dx.doi.org/10.1007/s12206-019-0321-2]

6. Gong, H.; Liu, J.; Ding, X. Study on the mechanism of preload decrease of bolted joints subjected to transversal vibration loading. Proc. Inst. Mech. Eng. Part B J. Eng. Manuf.; 2019; 233, pp. 2320-2329. [DOI: https://dx.doi.org/10.1177/0954405419838675]

7. Liu, Y.; Chen, J.; Zhang, X.; Tan, D. Fatigue behaviour of blind bolts under tensile cyclic loads. J. Constr. Steel Res.; 2018; 148, pp. 16-27. [DOI: https://dx.doi.org/10.1016/j.jcsr.2018.05.019]

8. Tizani, W.; Rahman, N.A.; Pitrakkos, T. Fatigue life of an anchored blind-bolt loaded in tension. J. Constr. Steel Res.; 2014; 93, pp. 1-8. [DOI: https://dx.doi.org/10.1016/j.jcsr.2013.10.002]

9. Wang, J.; Uy, B.; Li, D.; Song, Y. Fatigue behaviour of stainless steel bolts in tension and shear under constant-amplitude loading. Int. J. Fatigue; 2020; 133, 105401. [DOI: https://dx.doi.org/10.1016/j.ijfatigue.2019.105401]

10. Wang, Z.-Y.; Wang, Q.-Y.; Xue, H.; Liu, X. Low cycle fatigue response of bolted T-stub connections to HSS columns—Experimental study. J. Constr. Steel Res.; 2016; 119, pp. 216-232. [DOI: https://dx.doi.org/10.1016/j.jcsr.2015.12.009]

11. Hashimura, S.; Socie, D.F. A Study of Loosening and Fatigue Failure of Bolted Joints under Transverse Vibration. Trans. Jpn. Soc. Mech. Eng. Ser. C; 2006; 72, pp. 1297-1304. [DOI: https://dx.doi.org/10.1299/kikaic.72.1297]

12. Mínguez, J.M.; Vogwell, J. Effect of torque tightening on the fatigue strength of bolted joints. Eng. Fail. Anal.; 2006; 13, pp. 1410-1421. [DOI: https://dx.doi.org/10.1016/j.engfailanal.2005.10.012]

13. Benhaddou, T.; Chirol, C.; Daidie, A.; Guillot, J.; Stephan, P.; Tuery, J.-B. Pre-tensioning effect on fatigue life of bolted shear joints. Aerosp. Sci. Technol.; 2014; 36, pp. 36-43. [DOI: https://dx.doi.org/10.1016/j.ast.2014.03.003]

14. Guo, H.; Mao, K.; Liu, Y.; Liang, G. Experimental study on fatigue performance of Q460 and Q690 steel bolted connections. Thin-Walled Structures; 2019; 138, pp. 243-251. [DOI: https://dx.doi.org/10.1016/j.tws.2019.02.011]

15. Chakherlou, T.N.; Abazadeh, B. Investigating clamping force variations in Al2024-T3 interference fitted bolted joints under static and cyclic loading. Mater. Des.; 2012; 37, pp. 128-136. [DOI: https://dx.doi.org/10.1016/j.matdes.2011.12.037]

16. Chakherlou, T.N.; Abazadeh, B.; Vogwell, J. The effect of bolt clamping force on the fracture strength and the stress intensity factor of a plate containing a fastener hole with edge cracks. Eng. Fail. Anal.; 2009; 16, pp. 242-253. [DOI: https://dx.doi.org/10.1016/j.engfailanal.2008.03.002]

17. Yin, H.; Wang, T.; Yang, D.; Liu, S.; Shao, J.; Li, Y. A Smart Washer for Bolt Looseness Monitoring Based on Piezoelectric Active Sensing Method. Appl. Sci.; 2016; 6, 320. [DOI: https://dx.doi.org/10.3390/app6110320]

18. Huo, L.; Chen, D.; Kong, Q.; Li, H.; Song, G. Smart washer-a piezoceramic-based transducer to monitor looseness of bolted connection. Smart Mater. Struct.; 2017; 26, pp. 025-033. [DOI: https://dx.doi.org/10.1088/1361-665X/26/2/025033]

19. Khomenko, A.; Koricho, E.G.; Haq, M.; Cloud, G.L. Bolt tension monitoring with reusable fiber Bragg-grating sensors. J. Strain Anal. Eng. Des.; 2016; 51, pp. 101-108. [DOI: https://dx.doi.org/10.1177/0309324715598265]

20. Hu, W.; Liu, Z.; Liu, D.; Hai, X. Fatigue failure analysis of high speed train gearbox housings. Eng. Fail. Anal.; 2017; 73, pp. 57-71. [DOI: https://dx.doi.org/10.1016/j.engfailanal.2016.12.008]

21. Wang, Z.; Allen, P.; Mei, G.; Wang, R.; Yin, Z.; Zhang, W. Influence of wheel-polygonal wear on the dynamic forces within the axle-box bearing of a high-speed train. Veh. Syst. Dyn.; 2020; 58, pp. 1385-1406. [DOI: https://dx.doi.org/10.1080/00423114.2019.1626013]

22. Li, F.; Wu, P.; Zeng, J.; Liu, C.; Wu, H. Vibration fatigue dynamic stress simulation under multi-load input condition: Application to metro lifeguard. Eng. Fail. Anal.; 2019; 99, pp. 141-152. [DOI: https://dx.doi.org/10.1016/j.engfailanal.2019.02.024]

23. Wei, L.; Sun, Y.; Zeng, J.; Qu, S. Experimental and numerical investigation of fatigue failure for metro bogie cowcatchers due to modal vibration and stress induced by rail corrugation. Eng. Fail. Anal.; 2022; 142, 106810. [DOI: https://dx.doi.org/10.1016/j.engfailanal.2022.106810]

24. Dowling, N.E. Mechanical Behaviour of Materials; Pearson Prentice Hall: Old Bridge, NJ, USA, 2012.

25. Fatemi, A.; Shamsaei, N. Multiaxial fatigue: An overview and some approximation models for life estimation. Int. J. Fatigue; 2011; 33, pp. 948-958. [DOI: https://dx.doi.org/10.1016/j.ijfatigue.2011.01.003]

26. Dantas, R.; Correia, J.; Lesiuk, G.; Rozumek, D.; Zhu, S.-P.; Berto, F. Evaluation of multiaxial high-cycle fatigue criteria under proportional loading for S355 steel. Eng. Fail. Anal.; 2020; 120, 105037. [DOI: https://dx.doi.org/10.1016/j.engfailanal.2020.105037]

27. Macek, W.; Branco, R.; Trembacz, J.; Costa, J.D.; Ferreira, J.A.M.; Capela, C. Effect of multiaxial bending-torsion loading on fracture surface parameters in high-strength steels processed by conventional and additive manufacturing. Eng. Fail. Anal.; 2020; 118, 104784. [DOI: https://dx.doi.org/10.1016/j.engfailanal.2020.104784]

28. Macek, W. Post-failure fracture surface analysis of notched steel specimens after bending-torsion fatigue. Eng. Fail. Anal.; 2019; 105, pp. 1154-1171. [DOI: https://dx.doi.org/10.1016/j.engfailanal.2019.07.056]

29. Goldsmith, N.T.; Wanhill, R.J.H.; Molent, L. Quantitative fractography of fatigue and an illustrative case study. Eng. Fail. Anal.; 2019; 96, pp. 426-435. [DOI: https://dx.doi.org/10.1016/j.engfailanal.2018.10.013]

30. Shaohu, L.; Hao, Z.; Hui, X.; Quanquan, G. A new theoretical model of low cycle fatigue life for coiled tubing under coupling load. Eng. Fail. Anal.; 2021; 124, 105365. [DOI: https://dx.doi.org/10.1016/j.engfailanal.2021.105365]

31. Lee, H.W.; Djukic, M.B.; Basaran, C. Modeling fatigue life and hydrogen embrittlement of bcc steel with unified mechanics theory. Int. J. Hydrog. Energy; 2023; 48, pp. 20773-20803. [DOI: https://dx.doi.org/10.1016/j.ijhydene.2023.02.110]

32. Haiba, M.; Barton, D.C.; Brooks, P.C.; Levesley, M.C. Review of life assessment techniques applied to dynamically loaded automotive components. Comput. Struct.; 2002; 80, pp. 481-494. [DOI: https://dx.doi.org/10.1016/S0045-7949(02)00022-6]

33. Conle, F.A.; Mousseau, C.W. Using vehicle dynamics simulations and finite-element results to generate fatigue life contours for chassis components. Int. J. Fatigue; 1991; 13, pp. 195-205. [DOI: https://dx.doi.org/10.1016/0142-1123(91)90241-P]

34. Han, S.-H.; An, D.-G.; Kwak, S.-J.; Kang, K.-W. Vibration fatigue analysis for multi-point spot-welded joints based on frequency response changes due to fatigue damage accumulation. Int. J. Fatigue; 2013; 48, pp. 170-177. [DOI: https://dx.doi.org/10.1016/j.ijfatigue.2012.10.017]

35. Mršnik, M.; Slavič, J.; Boltežar, M. Frequency-domain methods for a vibration-fatigue-life estimation—Application to real data. Int. J. Fatigue; 2013; 47, pp. 8-17. [DOI: https://dx.doi.org/10.1016/j.ijfatigue.2012.07.005]

36. Palmieri, M.; Česnik, M.; Slavič, J.; Cianetti, F.; Boltežar, M. Non-Gaussianity and non-stationarity in vibration fatigue. Int. J. Fatigue; 2017; 97, pp. 9-19. [DOI: https://dx.doi.org/10.1016/j.ijfatigue.2016.12.017]

37. Huang, N.; Shen, Z.; Long, S.; Wu, M.L.C.; Shih, H.; Zheng, Q.; Yen, N.-C.; Tung, C.-C.; Liu, H. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci.; 1998; 454, pp. 903-995. [DOI: https://dx.doi.org/10.1098/rspa.1998.0193]

38. Wu, Z.; Huang, N. Ensemble Empirical Mode Decomposition: A Noise-Assisted Data Analysis Method. Adv. Adapt. Data Anal.; 2009; 1, pp. 1-41. [DOI: https://dx.doi.org/10.1142/S1793536909000047]

39. Yeh, J.-R.; Shieh, J.-S.; Huang, N. Complementary Ensemble Empirical Mode Decomposition: A Novel Noise Enhanced Data Analysis Method. Adv. Adapt. Data Anal.; 2010; 2, pp. 135-156. [DOI: https://dx.doi.org/10.1142/S1793536910000422]

40. Dragomiretskiy, K.; Zosso, D. Variational Mode Decomposition. IEEE Trans. Signal Process.; 2014; 62, pp. 531-544. [DOI: https://dx.doi.org/10.1109/TSP.2013.2288675]

41. Feng, Y.; Qu, S.; Li, F.; Dai, H.; Shu, K. Nonlinear vibration of the axle box front cover of high-speed train and its effect on connecting bolts. Eng. Fail. Anal.; 2023; 143, 106912. [DOI: https://dx.doi.org/10.1016/j.engfailanal.2022.106912]

42. Alhadid, M.A.; Wright, J.R. Estimation of mass and modal mass in the identification of nonlinear single and multidegree of freedom systems using the force-state mapping approach. Mech. Syst. Signal Process.; 1992; 6, pp. 383-401. [DOI: https://dx.doi.org/10.1016/0888-3270(92)90038-K]

43. Crawley, E.F.; O’Donnell, K.J. Force-state mapping identification of nonlinear joints. AIAA J.; 1987; 25, pp. 1003-1010. [DOI: https://dx.doi.org/10.2514/3.9733]

44. Nayfeh, A.H.; Mook, A.D. Nonlinear Oscillations; Wiley-Interscience: New York, NY, USA, 1995.

45. Nayfeh, A.H. The response of single degree of freedom systems with quadratic and cubic non-linearities to a subharmonic excitation. J. Sound Vib.; 1983; 89, pp. 457-470. [DOI: https://dx.doi.org/10.1016/0022-460X(83)90347-4]

46. Shi, X.; Lei, X.; Huang, Q.; Huang, S.; Ren, K.; Hu, Y. Hourly day-ahead wind power prediction using the hybrid model of variational model decomposition and long short-term memory. Energies; 2018; 11, 3227. [DOI: https://dx.doi.org/10.3390/en11113227]

47. Yam, L.Y.; Leung, T.P.; Li, D.B.; Xue, K.Z. Theoretical and experimental study of modal strain analysis. J. Sound Vib.; 1996; 191, pp. 251-260. [DOI: https://dx.doi.org/10.1006/jsvi.1996.0119]

48. HBM. How to Bond a Strain Gauge into a Bolt. Available online: https://www.hbm.com/en/7901/bond-strain-gauges-to-bolts/ (accessed on 20 March 2023).

49. Zahid, F.B.; Ong, Z.C.; Khoo, S.Y. A review of operational modal analysis techniques for in-service modal identification. J. Braz. Soc. Mech. Sci. Eng.; 2020; 42, 398. [DOI: https://dx.doi.org/10.1007/s40430-020-02470-8]

50. Sadeqi, A.; Amador, S.; Amador, R.; Jacobsen, J.; Ajina, H.; Brincker, R. Structural damage localization using a modified sensitivity analysis. Proc. Dan. Soc. Struct. Sci. Eng.; 2021; 91, pp. 95-115.

51. Parloo, E.; Verboven, P.; Guillaume, P.; Van Overmeire, M. Sensitivity-based operational mode shape normalisation. Mech. Syst. Signal Process.; 2002; 16, pp. 757-767. [DOI: https://dx.doi.org/10.1006/mssp.2002.1498]

52. VDI 2230:2015 Ingenieure, V.D. Systematic Calculation of Highly Stressed Bolted Joints—Part 1: Joints with One Cylindrical Bolt. Ingenieure VD-Gesellschaft Produkt- und Prozessgestaltung: Berlin, Germany, 2015.

53. BS EN 1993-1-9 (2005) Eurocode 3: Design of Steel Structures-Part 1–9: Fatigue. CEO: Brussels, Belgium, 2005.

54. GL 2010 Guideline for the Certification of Wind Turbines. Germanischer Lloyd Wind Energie GmbH: Hamburg, Germany, 2010.

55. Wang, J.; Song, C.; Wu, P.; Dai, H. Wheel reprofiling interval optimization based on dynamic behavior evolution for high speed trains. Wear; 2016; 366–367, pp. 316-324. [DOI: https://dx.doi.org/10.1016/j.wear.2016.06.016]