Fatigue Analysis of Sub-frame Based on Load Spectrum Extrapolation

: The extrapolation of the whole life load spectrum from the strengthened road load spectrum in the test site is the key factor for the accuracy of the fatigue durability evaluation of automobile components. In this paper, aiming at the limitation of parameter method extrapolation single distribution estimation, the non parametric extrapolation estimation is introduced to describe the comprehensive characteristics of the load spectrum distribution, and the two extrapolation methods are compared parametrically, and the optimal fitting extrapolation spectrum is selected to carry out the fatigue damage analysis of the sub-frame. Based on the real vehicle test in the proving ground strengthened road, the stress spectrum of the dangerous point of the sub-frame is extracted by co-simulation of multi-body dynamics and finite element method. The parametric method and non-parametric method are used to calculate the goodness of fit of stress spectrum, the optimal mean value and amplitude value probability density of stress spectrum are selected according to the goodness of fit, and the stress spectrum cycle was extrapolated to 10 6 times, and then the stress spectrum was divided into eight levels. According to Miner criterion, Haibach and NASA S-N curve correction models considering small load damage were respectively used to estimate fatigue damage of front sub-frame. The results show that the fitting effect of non-parametric extrapolation is better than that of parametric extrapolation, and the goodness of fit of equivalent amplitude is 98.6%, which improves the description accuracy of stress spectrum distribution characteristics and provides more accurate extrapolation spectrum for fatigue durability evaluation of components. Based on the extrapolated spectrum and using Haibach model and NASA model to study the fatigue damage of sub-frame, it is shown that NASA model is more suitable for reliability design of components, while Haibach model is suitable for lightweight design of components.


Introduction
The analysis and verification of the durability and reliability of automobile components is a key link in product development. The measured load spectrum of real cars on the strengthening road is the data basis [1]. Due to the limitation of factors such as test cost and development cycle, only a small amount of payload data can be collected, and on this basis, the full life load spectrum can be obtained by extrapolation technology [2].
The extrapolation technique is based on mathematical statistics to expand the load range and scientific expand of the number of cycles. Initially, normal distribution or Weibull distribution is commonly used by foreign experts to fit and extrapolate load spectrum amplitudes. This method only considers the effect of amplitude on fatigue damage and ignores the effect of the mean value. The extrapolated load spectrum obtained is not perfect. Normal distribution and Rayleigh distribution were first attempted by Kowalewski to fit the load spectrum average value, amplitude value and extrapolation. This method takes into account the influence of load mean value, realizes the conversion from one dimension to two dimensions, and improves the rationality of extrapolation load spectrum. In China, Gao and Yan [3] assumes that the amplitude obeys the log-normal distribution and the mean obeys the normal distribution . After gradual improvement, he is pointed out that the fatigue load amplitude of the mechanical structure obeys the Weibull distribution and the mean amplitude obeys the normal distribution. Based on this, Wan et al. [4] extrapolated the loader working device load and compiled the test load spectrum; Liu, Lv, and Li et al. [5] extrapolated the loader drive train load and established the experimental load spectrum; Zou, Li and Peng et al. [6] extrapolated the axle load of the vehicle in the test field conditions and compiled bench test load spectrum of transmission according to the power transmission system test specification; Gao, Xu and Fang [7] extrapolated the vehicle body torque load and a notation strategy for simplifying accelerated body bench fatigue tests is proposed, which achieved good results. Li et al. [8] extrapolated the compressed automobile transmission load by using the statistical extremum theory and compiled the bench test acceleration spectrum. In the above application research, the load distribution characteristics are fitted with a single distribution estimation method, which is more suitable for extrapolation with single peak load characteristics, but for some load spectrum of multi-peak or complex distribution, the load characteristics are difficult to reflect, and the mean value, amplitude value and corresponding frequency of load spectrum extrapolation are directly affected. However, the non-parametric extrapolation does not need to assume that the load obeys a certain distribution in advance, thus getting rid of the dependence on the load distribution itself, fully reflecting the distribution characteristics of the load spectrum to be extrapolated, so as to improve the accuracy of the extrapolated load spectrum [9][10][11][12] .
Based on the above discussion, in order to explore the difference and applicability of different extrapolation methods for the front sub-frame under the enhanced road conditions in the test field, in this paper, real vehicle test wheel center six-component force of test field strengthened road is used as basic data, and the whole vehicle multi-body model and the sub-frame finite element model as the data carriers, the above models were adjusted by the indoor bench test and the free mode test respectively； Based on virtual iterative method and inertial release method, multi-body dynamics and finite element simulation are carried out to solve the stress spectrum of sub-frame dangerous points; The parametric method and non-parametric method are used to fit and extrapolate the distribution characteristics of stress spectrum, and the extrapolation methods are comprehensively compared and analyzed from the perspective of goodness of fit and probability map distribution, so as to explore the technical route to improve the accuracy of extrapolation of enhanced road load spectrum. In addition, considering the large number of small load cycles in the extrapolation spectrum, the influence of S-N curve correction model, Haibach model and NASA model on the fatigue damage of the sub-frame is studied, which provides reference for the lightweight and reliability design of automobile chassis components.

Sub-frame Model Establishment
The sub-frame, an important load-bearing part of the chassis, has the function of vibration reduction and support. It is connected with the body, lower arm and other components through bushing, which is the support device of various major chassis systems. The sub-frame of the test vehicle is plate stamping and welding structure, and the steel plate material is SAPH440, which is composed of front beam assembly (plate thickness of 3 mm), left and right longitudinal beam assembly (plate thickness of 2 mm), and rear cross beam assembly (plate thickness of 2.5 mm). The sub-frame solid model was obtained by 3D scanning reverse modeling, and Hypermesh was used for mesh generation, quality inspection and assign material properties. The welding seam is simulated by welding element (acm), and the posture of the finite element model of the sub-frame is adjusted by hard point coordinates, as shown in Figure 1. The sub-frame finite element model can be verified by free mode test, as shown in Figure 2 [12]. Firstly, the modal analysis of the finite element model of the sub-frame is carried out by using Hypermesh software, and the natural frequencies and modal shapes of the first six non rigid modes are obtained. The corresponding modal force hammer is selected according to the natural frequency obtained by simulation, and the appropriate suspension point, excitation point and response point are selected according to the modal analysis cloud image, the suspension point should select the position with small response and convenient suspension, the excitation point and test point should avoid the part with small response and select the part with larger response as far as possible.
The single point excitation method is used to test, the excitation point is knocked by modal force hammer and the excitation force is collected, the test point is arranged with unidirectional acceleration sensor to collect the response vibration signal, DHDAS software is used to calculate the frequency response function of the sub-frame and get the modal parameters. According to the natural frequency range of simulation is less than 400 Hz, so rubber hammer (frequency range: 0 ~500 Hz) is selected. Determine the appropriate excitation force by trial tapping, according to Nyquist-Shannon sampling theorem, the sampling frequency must be more than 2 times of the maximum natural frequency of the sampled samples [13], otherwise it will lead to sampling distortion. In general practical application, the sampling frequency is guaranteed to be 5~10 times of the highest signal frequency, so the sampling frequency of the test is set as 5 kHz. By comparing the statistical data of the first six modes (as shown in Table 1), the relative error of the natural frequency of simulation and test modes is less than 10%, and the finite element model of the sub-frame meets the accuracy requirements [14].

Establishment and verification of Vehicle
Multi-body Dynamic model A SUV was selected for the test vehicle. The weight of the whole vehicle is 1531 kg, the wheelbase is 2680 mm, the track is 3574 mm, front suspension is macpherson, rear suspension is double wishbone. Creaform hand held 3D scanner was used to obtain the whole chassis and component point clouds; Based on CATIA software, the chassis hard point parameters and parts quality characteristic parameters are obtained; the full load and no-load real vehicle centroid coordinates are calculated through the static load test and the longitudinal lifting method uniaxial lifting test. Based on the above parameters, the vehicle multi-body dynamics rigid-flexible coupling model is built, as shown in Figure 3. The communicator is used to connect the subsystems and components. During the movement, the stabilizer bar has large deformation and needs flexible processing. The car body is simplified to spherical geometry, and the actual mass characteristic parameters are given respectively.
In order to ensure the accuracy of the vehicle multi-body model, the suspension K&C characteristics and real vehicle dynamic characteristics are verified respectively. The suspension K&C characteristics of the vehicle multi-body model are verified by using the wheel center vertical displacement, wheel center vertical force, toe angle, me-domain curve of toe angle changing with wheel center vertical displacement. By adjusting the damper damping, spring stiffness, limit block stiffness and other parameters, the calibration results are shown in Figure 4. Compared with the time domain curves of test and simulation, the variation trend of the two is consistent, which shows that the property parameter of elastic element and the connection relation of parts in the multi-body model are correct.       Table 2. For vehicle signal acquisition, wheel center six component force analyzer, acceleration sensor and displacement sensor are used to collect test signals such as wheel center six component force, axle head vertical acceleration, suspension upper fulcrum vertical acceleration, sub-frame acceleration, suspension spring displacement, etc, the sensor layout is shown in Figure 5 Based on ADAMS multi-body dynamics software, the test conditions in Table 2 are simulated, and the signals except the six component force of wheel center are taken as simulation monitoring signals. According to the comparison results of simulation data and monitoring signals, as well as the influence trend of parameters of the motion equation (Eq. (1)) of the single-wheeled vehicle model, the main parameters of the multi-body model, such as damper damping and suspension spring stiffness, are adjusted Where: Fz is the wheel center vertical force (N); mw is the un-spring mass (kg); K is the spring stiffness (N/m); C is the shock absorber damping (N·m/s); mb is the spring upper mass (kg); W is the suspension spring displacement (mm); aw is the shaft head acceleration (mm/s 2 ); ab is the suspension fulcrum acceleration (mm/s 2 ). When the time history curves of simulation signal and monitoring signal approach to the same, and the relative error of RMS value between them is within 30% [15], the initial adjustment of the multi-body dynamic model is completed. The bench test marking results of the right front wheel monitoring signal of the vehicle under the condition of random white noise in the bench test are shown in Figure 6. The relative error statistics of RMS values in the time-history curve are shown in Table 3 Fig. 5 Indoor bench vibration test

Calculation of Dangerous Point Load Spectrum of Control Arm
In this paper, the durability road test in the test field is taken as the data source, the road test was conducted by the test driver in accordance with the Specifications for the Testing of Automobile Product Stability and Reliability in Highway Transportation Proving Ground of the Ministry of Communications, the road condition information is shown in Table 4. The vehicle counterweight scheme is based on driver's seat driver 80 kg, co-driver's seat occupant 75 kg, and two 70 kg occupants in the rear seat, the vehicle weight is 1831 kg, and the front and rear axle mass distribution coefficient is 1.1. The test vehicle is equipped with sensors and signal acquisition system, such as wheel center six component force meter, wheel center acceleration, suspension spring displacement, shock absorber upper pivot acceleration, and the acquisition frequency is set 1000 Hz. Combined with the six component force of wheel center, acceleration of wheel center, displacement of suspension spring and acceleration of upper fulcrum of shock absorber and other test data measured by the above real vehicle test, the dynamic simulation of vehicle multi-body model is carried out based on ADAMS multi-body dynamics software. Since the measured excitation frequency of 0 ⁓ 50 Hz is close to the first natural frequency of the sub-frame 58.9 Hz, it is necessary to consider the resonance of the components caused by the excitation frequency and aggravate the fatigue damage. In this paper, based on the displacement reverse loading method, the multi-body dynamic simulation is carried out by using the virtual iteration method, and the modal displacement of the front sub-frame is obtained, and the virtual iteration formula is Where: n is the number of iterations; un is the driving signal for n iterations; un-1 is the driving signal for n-1 iterations; F -1 is the inverse of the transfer function of the multi-body dynamics model; yd is the benchmark monitoring signal: yn-1 is the output response signal for iteration n-1 times. The problem of simulation distortion when multi-body dynamic load simulation wheel center applies vertical force is displacement signal. Forming five signal combination of force and a vertical displacement was used as the driving signal un-1; the measured six-point force and acceleration of the wheel center, the spring displacement, and the acceleration of the sub-frame are used as the benchmark monitoring signal yd. Iterative operation makes yn-1 and yd approximate. If the pseudo damage ratio between the two reaches 0.5 ~ 2, indicating that the iteration is convergent. Then, un-1 will be used as the final driving signal after iteration, the modal displacement of the sub-frame can be obtained, as shown in Figure 7, which is the modal displacement of the sub-frame of 7~9 order. After the modal displacement is obtained, the fatigue simulation of the sub-frame can be carried out based on the modal displacement recovery method. The modal displacement recovery method takes the undamped modes of the parts as the spatial base, calculates the response of the system by superimposing the contributions of the modes of each order, takes into account the influence of the excitation frequency on the parts, and can better solve the vibration fatigue problem of the vehicle parts. The modal stress of each node of the part can be obtained by applying the mode superposition method. The solution formula is as follows: Where: Eσ is the modal stress matrix; Φ is the modal displacement vector; σ is the modal stress of each node. Based on the nCode software platform, the fatigue simulation analysis of the sub-frame is carried out under the interface of Designlife. The fatigue damage dangerous point is the area where fatigue failure occurs the earliest and the accumulated damage is the largest. The fatigue life of the component itself is often evaluated by the fatigue failure time at this point in engineering. The damage danger points of the front sub-frame finite element model are identified through the Hot Spot Detetion hot spot detection function, and the stress load of the dangerous points is extracted by using the unique Virtual Gauge function of nCode software, as shown in Figures 8 and 9.

Parameter Estimation of Load Spectrum Distribution
Using statistical principles to process the load time history is an important part of load spectrum compilation. The common counting methods are divided into two categories: single-parameter counting and dual-parameter counting. The single-parameter counting method considers only one variable in the load cycle. It cannot describe the characteristics of the load cycle, but only an approximate description of random loads. The two-parameter counting method can record two parameters in the load cycle, and can record all the information of the load cycle, it is a better counting method, and its representative is the rain flow count method, this method can obtain all the information of the mean and amplitude, the counting principle is similar to the cyclic stress-strain of metal parts due to the actual working load, and has a solid mechanical basis [17] rainflow count statistics on the signed Von-Mises stress time history to get the mean, range and corresponding frequency, according to the relationship between the range and the amplitude: amplitude (σa) = range (Sa) / 2, the average and amplitude can be obtained three-dimensional histogram, as shown in Figure 10. Stress spectrum compression is achieved by deleting loads less than 10% of maximum amplitude [18], and obtaining the compression signed Von-Mises stress spectrum, the mean and amplitude frequency information of the compression stress spectrum can be calculated, and a two-dimensional histogram of the mean and amplitude can be obtained and carry out independence test. According to Fisher's theorem [19], if the mean and amplitude are independent of each other, they follow a  2 distribution with degrees of freedom (r-1)× (s-1), that is Where: k is the total number of times of stress spectrum cyclic load; i, j are the levels of the average and amplitude divided by the application of the rain flow counting method, and they are divided into 64 levels; ki is the frequency of the average value at the i-th level; kj is the frequency of the amplitude at the j-th stage; kij is the frequency of the cyclic load corresponding to the average of the i-th stage and the j-th level.
From the nature of the  2 distribution, when the degrees of freedom (r-1)×(s-1)≥45, the  2 distribution approximately follows the normal distribution of N(m,2 m), m=(r-1) ×(s-1). Therefore, the upper a quantile of the  2 distribution can be obtained according to formula (5).
Where: m is the degree of freedom of the chi-square distribution, Ua can be obtained from the standard normal distribution, and the degree of freedom of the  2 distribution is 3969. According to equation (5), the distribution quantile 2 0.05 (3969)  =4115.6. By formula (4), the  2 value is less than 4115.6, and it can be judged that the average and amplitude are independent of each other at the 95% confidence level. After determining the independent conditions of mean and amplitude, two-dimensional bar graphs of mean and amplitude are fitted with parameters respectively. The mean-amplitude rain flow matrix is obtained from the rain flow counts. The rain flow matrix is transformed into mean statistical histogram and amplitude statistical histogram, and then load average and amplitude statistical histogram are transformed into corresponding probability density histogram. According to the histogram probability density solution formula (6), the probability density function corresponding to the load spectrum mean and amplitude can be solved. Then, the joint probability density function is established. Then, according to the parameter extrapolation method, the normal distribution and Weibull distribution are respectively used to fit the probability density histogram of the load spectrum mean and amplitude, and the parameter estimation is carried out Where, f represents the probability density corresponding to the load; h represents the group spacing of load division, with mean group spacing of 7.5 and amplitude group spacing of 3.8. The random variables x and y are used to represent the load average and amplitude respectively. The expressions of the probability density function of the mean and amplitude are shown in Formulas (7) and (8). Since the mean and amplitude are independent, the joint probability density function of the load spectrum can be expressed in Formula (9) Fatigue Analysis of Sub-frame Based on Load Spectrum Extrapolation (9) where, f(x) is the mean probability density function; μ is the normal distribution mean; σ is the normal distribution standard deviation; f(y) is the amplitude probability density function; β and α are the shape parameters and scale parameters of the Weibull distribution respectively; f(x,y) is the combined probability density function of the load.
Based on the self-programming of damped least squares method, the mean and amplitude are fitted by normal distribution and Weibull distribution respectively. The fitting effect is shown in Figure 11. The goodness of mean fitting is 96.6% and that of amplitude fitting is 97.5%. The results of parameter estimation are obtained, mean value μ=-13.82, standard deviation σ=30.27, shape parameter β=0.96 and scale parameter α=26. 22.
By observing Figure 11, it is found that the fitting effect of mean and amplitude is well. In order to further evaluate and analyze the fitting effect, normal distribution fitting test and Weibull distribution fitting test are carried out for mean and amplitude respectively. The results are shown in Figure 12, load amplitude can well obeys the Weibull distribution except for lower amplitude, but the load mean only obeys the normal distribution in the [-80,54] interval. Therefore, there is deviation in fitting load spectrum by single distribution parameter method, and the load data should not be extrapolated by parameters.

Non-parametric load spectrum extrapolation
Non-parametric extrapolation is also called non-parametric estimation extrapolation, of which the most typical one is nuclear density estimation method extrapolation. The previous study shows that there is a deviation between the mean value of load spectrum extrapolated by parameter method and the mean value of actual load. Therefore, this paper continues to use nuclear density estimation method to extrapolate loads. This method has high data fitting accuracy and gets rid of the dependence on parent samples. According to the different dimensions of the kernel density function, it can be divided into one-dimensional kernel density estimation method and two-dimensional kernel density estimation method. In this paper, the one-dimensional kernel density estimation function is mainly used to carry out non-parametric extrapolation of the load spectrum.

Dimensional reduction of load spectrum
The mean-amplitude rain flow matrix obtained by rain flow counting includes two parameters, mean and amplitude. In order to avoid discussing the correlation between load mean and amplitude, the two-dimensional mean-amplitude load spectrum can be converted into one dimensional equivalent amplitude-frequency load spectrum by using Goodman average stress correction formula (10) 1 UTS m eq a   S S S S (10) Where: Sa is the stress amplitude; Sm is the stress mean value; Seq is equivalent stress amplitude; SUTS is the strength limit. Then, the one-dimensional equivalent amplitude -frequency load spectrum is compressed and 10% of the maximum amplitude is taken as the deletion basis of invalid amplitude. The equivalent stress amplitude-frequency histogram is shown in Figure 13.

One-dimensional nuclear density estimation
The kernel density estimation method is fitted according to the sample data, which can retain the distribution characteristics of the data itself, describe the dispersion data well, and approach the probability density of different distributions. Assuming that X is a random variable, [X1,X2,X3...XN] is the sample point of X, then the probability density function at any point Xi is: Where: h is the bandwidth; N is the number of Xi in the interval [x-h, x+h]. To define the core function K(x), K(x) needs to satisfy the following three conditions: Then the probability density function of Xi is estimated by nuclear density: (13) When estimating nuclear density, both the type of core function and the bandwidth have an effect on the results. The type of core function controls the utilization of data points considered when estimating the probability density f(x) of point x. The types of kernels control the utilization of data points considered when estimating the probability density f(x) of point X. There are many types of kernels, common types of kernels include triangular kernels, Gaussian kernels and Epanechnikov kernels, etc. At present, there is no recognized optimal kernels function [20]. When the bandwidth is constant, the difference of probability density functions generated by different kernels is very small. However, bandwidth has a great influence on the estimation results. When bandwidth is too small, the convergence of probability density function is poor and error peak is generated; When the bandwidth is too large, it will make it difficult for the probability density function to fit the multi-peak distribution features. The obtained probability density function is usually too smooth and important features of the load distribution are lost [21]. Therefore, it is necessary to use the optimal bandwidth kernel function to fit the load spectrum histogram. The methods to solve the optimal bandwidth are generally divided into the cross validation method and the insertion method. In this paper, based on the insertion method, formula (14) is used to solve the optimal bandwidth of the Gaussian kernel function. 5 1 06 σn h (14) Where: h is the bandwidth, σ is the standard deviation of the equivalent amplitude data, and n is the total number of load cycles. According to the equivalent amplitude load data, σ=12.5, h=5413, the optimum bandwidth h=2.3 of Gauss core function is obtained by substituting the above formula. Based on the optimal bandwidth, the equivalent amplitude load is fitted by using Gauss core function, and the fitting goodness R 2 is used as the test index of fitting effect. The difference between the mathematical model and the measured data is quantitatively described from the aspect of fitting accuracy. As R 2 approaches 1, the fitting accuracy gradually increases. The fitting test results of Gaussian kernel function are shown in Figure. 14, when the load value is more than 20 MPa, the Gauss core function can better fit the equivalent amplitude load, and the goodness of fitting is 98.6%.

Equivalent Amplitude Load Extrapolation
Due to the limited test data, it is necessary to use the extreme load derived from the statistical characteristics of the load spectrum to reasonably reflect the load distribution characteristics of the test field. According to the probability criterion proposed by Conover, the probability P of occurrence of extreme load in 10 6 load cycles is 10 -6 , and the equivalent amplitude load extreme value Xmax is 175MPa through Eq (15). The original load cycle times are extrapolated to 10 6 times, and the corresponding load cycle times are obtained by formula (16). The corresponding load cycle frequency cumulative curve is plotted by MATLAB programming as shown in Figure 15.
x dx x f P (15)  dx Where: Ni is the number of cycles per stage load.

Load spectrum classification
Due to the large amplitude density of the extrapolated spectrum obtained, the load spectrum can be graded under the premise that the damage amount of the load spectrum remains unchanged. The continuous load cycle frequency accumulation curve can be transformed into a frequency accumulation curve with a stepped type. Conver found that the load spectrum can be divided into 8 stages to accurately reflect its fatigue effect, which has been applied in engineering practice [22]. Since the damage caused by large amplitude is much higher than that caused by low amplitude, the amplitude is divided according to unequal distances, and the amplitude ratio coefficient G is 1, 0.95, 0.85, 0.725, 0.575, 0.425, 0.275 and 0.125 respectively, that is, xi=g xmax for each amplitude. The load spectrum of 8 stages is shown in Table 5.

Extrapolation Spectrum Fatigue Damage Analysis
The fatigue of sub-frame belongs to high cycle fatigue, so the nominal stress method is used to estimate its fatigue life. In order to estimate the fatigue life of the sub-frame, the corresponding S-N curve must be obtained. Generally, the S-N curve of component is obtained by modifying the S-N curve of material. In this paper, the part of fatigue life N>10 3 is modified, that is, the high-cycle fatigue zone is corrected, and the correction effect changes linearly in the high-cycle fatigue zone. The fatigue limit of component decreases by KD times compared with the material fatigue limit, and its expression is as follows: Where: CD is the dimension factor; Cs is the surface processing factor; Kf is the fatigue notch factor, which is obtained from the following formula (18) Where: q is the notch sensitivity coefficient; Kt is the theoretical stress concentration factor. By consulting the references and taking CD=0. 7 (19) Where: σ is nominal stress (MPa), N is fatigue life corresponding to nominal stress, and A and B represent parameters related to material properties.
Assuming that the fatigue limit of the front accessory frame is σ-1,-1, the corresponding life is 10 7 , σ-1,-1=0.5σb/KD, B is the tensile strength of SAPH430, σb=430MPa, the corresponding coordinate point of the fatigue limit of the front accessory frame on the S-N curve is obtained (79.6, 10 7 ). When the fatigue life is 10 3 , the corresponding load σ1000=0.85CRσb, and CR=1 represents the level coefficient of reliability. From this, the corresponding coordinate point of σ1000 on the S-N curve of the component is (365.5, When carrying out fatigue analysis on components, it is first necessary to select an appropriate fatigue cumulative damage theory to calculate the total load damage. The existing fatigue cumulative damage theory can be divided into four categories: linear cumulative damage theory, bilinear cumulative damage theory, nonlinear cumulative damage theory and other cumulative damage theory (empirical or semi-empirical formula induced by experimental analysis) [23] Miner law, as the most representative theory in the linear cumulative damage theory, is simple in form and easy to calculate [24], and has been widely applied in modern engineering practice. Its expression is as follows Where: D is the cumulative total damage; i is the load level; ni is the actual number of cycles of the stress in stage i; Ni is the fatigue life corresponding to the stress in stage i. Considering that the extrapolation spectrum contains a large number of loads lower than the fatigue limit of components, while Miner law neglects the damage caused by loads below the fatigue limit to components, it is necessary to modify the traditional S-N curve model and obtain two S-N curve correction models considering small load cycle: Haibach model and NASA model, as shown in Figure 16.
Based on Miner damage accumulation theory and the corresponding S-N curve correction model, the fatigue damage corresponding to each stage load in the eight-stage load spectrum is obtained, as shown in Table 6.
The extrapolated eighth-order spectrum is defined as a load block , and the cyclic loading times of the spectrum block are used as the evaluation measurement of the fatigue life of the front sub-frame. According to Table 6 above, the damage amount of a single spectrum block is calculated, and the corresponding total damage amount and the fatigue life of the front sub-frame are obtained by combining the standard durability test, as shown in Table 7. It can be found that the damage amount of the standard durability test of the front sub-frame is less than the critical damage specified by Miner theory and meets the durability requirements.
In combination with Tables 6 and 7, it can be seen that in the two S-N curve correction models considering small load cycles, load damage amount of below the fatigue limit calculated from NASA model accounted for 17.1% of the spectrum block damage, while Haibach model only accounted for 1.54%. This shows that the damage weight given by NASA model to small load is higher, the damage estimated during fatigue analysis is larger, and the fatigue life obtained is more conservative than that of Haibach model. It is beneficial to reduce the risk of fatigue failure and is suitable for component reliability design. Haibach model gives a moderate weight to minor damage, and estimates minor damage, which is beneficial for component lightweight design.

Conclusion
(1) The boundary loads of sub-frame which are difficult to collect by sensors are obtained by virtual iteration and multi-body dynamics simulation method, and the position of dangerous damage points of sub-frame is accurately calculated by modal displacement recovery method, taking into account the influence of excitation frequency on components. (2) The real-vehicle measured strain is synthesized into Signed Von-Mises stress time history, which is more in line with the true stress state of the component under complex loading conditions. The rain-flow counting method is used to convert the stress-time history into the mean-amplitude rain-flow matrix which can be used for data analysis. The two-dimensional load spectrum is transformed into one-dimensional equivalent load spectrum by Goodman mean stress correction method, thus avoiding the discussion of the correlation between load average and amplitude. The applicability of parametric and non-parametric extrapolation to durable road load spectrum is compared and studied. Based on parametric extrapolation, the mean data do not fully follow the normal distribution, the fitting effect of load spectrum extrapolated by Gaussian kernel function non-parametric method with equivalent amplitude load is better, and the goodness of fit is 98.6%, which can better describe the distribution characteristics of load spectrum. (3) Extrapolation spectrum are divided into eight levels by amplitude ratio coefficient method, and S-N curves of sub-frame components are obtained by correcting material S-N curves. Considering the effect of a large number of small loads lower than fatigue limit on 6fatigue damage in the eighth stage spectrum, based on Miner law, sub-frame fatigue analysis is carried out by using Haibach model and NASA model respectively.