A comparative study on the fatigue life of the vehicle body spot welds using different numerical techniques: Inertia relief and Modal dynamic analyses

Owing to the expensive and time-consuming nature of durability experiments, finite element based durability analysis is quite prevalent in the automotive industry. Numerical fatigue life analyses are typically divided into two different categories, the quasi-static methods which are faster and the dynamic methods which are more accurate. The aim of this paper is to compare the inertia relief and modal dynamic approaches in terms of formulation, accuracy and computation time. The chosen case study is the fatigue life of the vehicle body which is considered the main load-bearing component in a vehicle. By utilizing multi-body dynamics model and driving the vehicle on different standardized roads and by different velocities, the loadings, which act on the body are calculated and later used for the stress analysis. Then, by using the structural stress method, the fatigue life of the vehicle spot welds is calculated and the results are compared for both approaches. The findings reveal that the modal dynamic method is almost 37 times more time-consuming than the inertia relief approach, but if accuracy is desired, it can be up to 96% more accurate. Also Citation: Ahmadi, A., Farrahi, G.H., Reza Kashyzadeh, K., Azad, Sh., Jahani, K., A comparative study on the fatigue life of the vehicle body spot welds using different numerical techniques: Inertia relief and Modal dynamic analyses, Frattura ed Integrità Strutturale, 52 (2020) 67-81. Received: 31.10.2019 Accepted: 09.01.2020 Published: 01.04.2020 Copyright: © 2020 This is an open access article under the terms of the CC-BY 4.0, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. A. Ahmadi et alii, Frattura ed Integrità Strutturale, 52 (2020) 67-81; DOI: 10.3221/IGF-ESIS.52.06 68 as predicted, at low frequency loading, there is no major difference between the results of both methods.


INTRODUCTION
ue to the competitiveness of the automotive industry, the production of high-quality vehicles is a necessity to sustain the fast-growing auto market. One of the main considerations in vehicle design is the durability of the components and the joints which reflects the quality of the manufactured product. The vehicle body is considered the main load-bearing component in a vehicle and therefore the durability of its components and joints is important. A vehicle Body-In-White (BIW) is referred to a stage in the production line where the body sheets are joined together and the main structure of the vehicle is formed. The BIW does not include the closures i.e., doors, fenders, hood, and trunk lid. BIW sheets are mostly joined by spot welds and a typical BIW contains about 5000 spot welds. Failure of these joints due to fatigue phenomenon changes the vibro-acoustic behavior of the vehicle, leading to unwanted noise and vibrations felt by the passengers [1,2]. Moreover, as it was demonstrated by Xiang et al. who analyzed a spot welded thin-walled component, failure of some spot welds degrades the crashworthiness behavior of the structure [3]. Therefore, knowledge of the BIW fatigue behavior is an important consideration for auto manufacturers. Because of the time consuming and costly nature of experiments (like using the full body fatigue test rigs [4]) and also lack of prototypes in the early stages of the vehicle design, Finite Element Method (FEM) is widely used for stress and fatigue analysis of vehicle components. There are two main methods for the numerical stress analysis of the BIW. The first is the quasi-static method that is mainly implemented using inertia relief approach. In brief, the inertia relief method uses the rigid body accelerations and the d'Alembert's principle to analyze an unconstrained structure like the BIW. An example of this method is the research done by Wen and Du [5] who analyzed the fatigue behavior of a mini car body. A similar research was also done on a truck cab by Chen et al. [6]. There are also other studies on the vehicle body and its parts which used the quasi-static method [7][8][9]. Although this method of analysis is quite easier and faster than other methods [10], it has the big disadvantage of neglecting inertia effects. There are many parts attached to the BIW like the battery, some parts of the powertrain, the seats, the doors, etc. Since these parts are composed of highly concentrated masses, near their mounting locations, they may have very low natural frequencies which lie in the loading frequency spectrum and cause resonance [11]. To take the resonance and inertia effects into account, dynamic methods should be utilized in the analysis. Dynamic analysis of the linear systems with high degrees of freedom like the BIW is mainly done using mode-based dynamic solvers since other methods are not economical in such a case. In brief, the modal dynamic method uses the superposition of the structure's mode shapes to calculate the response of the structure under loading. There are multiple examples in the literature that have applied this method for the analysis of vehicle components. For example, Farrahi and Khalaj used the modal dynamic technique (also known as modal superposition method) to analyze the durability of the rear spindle of the vehicle [12]. Wang [13] and Jordan [14] used this method for the analysis of the vehicle BIW from the fatigue point of view and Lu et al. [15] used the same method on a high-speed train bogie frame. Also, in a more thorough analysis [16], the vehicle body was analyzed using modal dynamic method and the regions having a high probability of fatigue failure were successfully modified to mitigate the stress intensity. In addition, there were other hybrid methods which utilized both static and dynamic analyses simultaneously [17,18]. Based on the results published in a previous research by Anvari and Beigi [19], the use of quasi-static method for the stress analysis of the BIW is only allowed if the loading frequency is less than 10% of the first nonzero natural frequency of the structure. Although the loading applied in their research was fully reversed sinusoidal force at a specific node (not a real loading condition), it gave a good understanding that how much difference inertia effects can make in the results. In order to have a finer conclusion on the effect of incorporating inertia into the durability analysis of BIW, the present research is aimed at comparing inertia relief and modal dynamic analyses formulations and their results for the fatigue life of the vehicle body. To this end, the Multi-Body Dynamics (MBD) model of the full vehicle is driven on different roads that are standardized by the International Organization for Standardization (ISO) at three different speed regimes and the D loadings on the suspension-body interface joints are calculated. Afterward, the loadings are used by the finite element model to calculate the stress time histories using both inertia relief and modal dynamic methods. Finally, the fatigue analysis is performed using the structural stress method (SSM) which uses the structural stress parameter and the master S-N curve of the material to calculate the fatigue life of the spot welds. The numbers of failed spot welds along with the computation time are the outputs of the simulations and are compared for both methods.

EXTRACTING THE SUSPENSION-BODY INTERFACE JOINT REACTIONS USING MBD ANALYSIS
ig. 1 shows the MBD model of the vehicle. It consists of a MacPherson front suspension and a torsion beam rear suspension system. The total mass of the vehicle is 1110 kg and a 75 kg driver is also assumed for the analysis. By driving the model on different Virtual Proving Grounds (VPG), the forces and moments at the joints where the suspension system is attached to the body, are extracted. These joints are represented by 12 points located at the front and rear shock towers bushing, Lower Control Arms (LCA) bushing, anti-roll bar bushings, strut rods bushing, and torsion beam bushings. 36 force time histories and 36 moment time histories are extracted at these 12 points and will be used in the finite element model for the stress analysis. The roads used in this research are based on ISO 8608 standard [20] which categorizes the roads by their roughness Power Spectral Densities (PSD) according to Eq. 1: where n is the spatial frequency and n0 is the reference spatial frequency being set asπcycle/meter. The road classes are defined by the magnitude of Gd(n0) as listed in Tab. 1. The road class A is very smooth and usually does not cause damage to the BIW. In order to gain a better insight into the amount of road roughness attributed to each road class, the vertical displacement of a wheel passing at a 10 km/h velocity for a 60-second long trip is shown in Fig. 2   In this study, three different maneuvers are used in the analysis as follows: The road percentages are rounded numbers taken from Reza Kashyzadeh et al. [22]. The road class F is neglected due to its very small share of the total road (1.6%) and its percentage is allocated to the road class E.

FINITE ELEMENT MODEL OF THE VEHICLE BODY
mesh of 656457 shell elements, of which less than 4% are triangular, is used in the finite element model of the vehicle body. The utilized element types are three-node shell element (S3) and four-node shell element with reduced integration option (S4R). Considering the fact that the vehicle body consists of many parts and the critical locations on the body are not readily recognizable, a mesh sensitivity analysis cannot be carried out on the body. A prevalent approach in the auto industry is to refine the mesh as much as possible around joints and holes on the body since these points are responsible for the stress concentration and are possible failure locations on the body. This approach is adopted in this study. 3498 spot welds are created by using Timoshenko beam elements having circular cross-sections with a diameter equal to that of the actual weld nugget. The reason why Timoshenko beam elements are used is the fact that spot weld nuggets behave much like a short beam and as a result, the shear effects have a significant effect on the analysis results. The detail of the adopted spot weld model is shown in Fig 3. Moreover, instead of assuming some predefined diameters for the spot weld nuggets which is proved to change the fatigue analysis results by orders of magnitude [23], the authors used the following definition for the nugget diameters which is used by the quality control units of some car manufacturers [24].
where tmin is the minimum sheet thickness of each pair in a 2-sheet or 3-sheet connection. According to the utilized standard, spot welds having diameters less than 4 mm will fail the quality control tests. The sprung accessories attached to the vehicle body (front and rear doors, hood, trunk lid, spare tire, front and rear seats, fuel tank, and the driver) are considered as concentrated masses at the respective centroids and are connected to the body via connecting elements [25].

CALCULATION OF STRESS TIME HISTORIES
Inertia relief method olving a Finite Element (FE) problem requires the model to be properly constrained. Some structures like vehicle body and planes do not have specific boundary conditions and in other words, are floating objects. In order to analyze such structures, researchers introduced the inertia relief method which uses a set of self-equilibrating forces and prevents the rigid body motion of the unconstrained structure. By adding the rigid body d'Alembert's forces to the structure, the finite element problem can be stated as: where the term {P}-[M]{ü b } is the so-called self-equilibrating system of forces which puts the structure at rest. The rigid body acceleration of the system is a function of the reference point acceleration vector { r } and the rigid body transformation matrix of the system [T] as: Note that [T] is the response of the system to unit accelerations in all 6 directions (3 translations and 3 rotations). For example, for a node in the 3D space, [T]6×6 is defined as: Where x, y and z are the coordinates of the node and x 0 , y 0 and z 0 are the coordinates of the reference point. For a structure having n number of nodes, the [T] matrix is defined as: ( 2) 6 6 ( ) Using [T], the external and inertia forces on the structure can be transferred to the reference point. This is done as follows: where {P0} is the resultant external force vector on the reference point. If the structure is at rest, the resultant external forces and the inertia forces must be equal as stated by Eq. 10.
Hence, by using Eq. 5, the rigid body acceleration of the system is found to be: Inserting Eq. 12 into Eq. 4 and factoring out {P} in the right side yields: where, where [I] is the unit matrix and [R] is called the inertia relief projection matrix. If [R] is multiplied by any external force vector, it produces a system of self-equilibrating forces. Eventually, any adequate number of boundary conditions that fully constrain the structure can be used to solve the system of Eq. 14. The imposed boundary conditions do not create local stresses since all external forces are neutralized by internal equilibrating forces.

Modal dynamic method
In order to perform a mode-based dynamic analysis, the natural frequencies and mode shapes of the structure are necessary. The modal dynamic method uses these mode shapes and calculates the dynamic solution of the system by superposing them. By extracting the modes and using modal transformation, the nodal displacement vector can be expressed as: where {u} is the displacement vector, [] is the modal matrix and consists of mode shapes and {} is the modal coordinates vector. By ignoring the damping [8,13,26,27], the equation of motion can be written as follows: Inserting Eq. 15 into Eq. 16 and multiplying both sides by [] T , yields: where [n * ] is the modal stress matrix calculated using the nodal displacements in the mode shapes and N is the number of modes utilized in the approximation.
As it was mentioned earlier, the modal dynamic method utilizes the mode shapes of the structure to calculate the dynamic response by superposing the mode shapes. From the mathematical point of view, the response is completely precise if all mode shapes of the structure are extracted and used in the summation. However, since extracting all mode shapes of a structure like BIW with more than a million degrees of freedom is not possible, only a small subset of modes are used in the summation. In order to find the suitable number of modes for reasonable accuracy of the model, some parameters should be defined beforehand. For the -th mode, the generalized mass is defined as follows: In the next step, this parameter is used to calculate the Modal Participation Factor (MPF): where  i is the MPF of the mode  in the i-th direction (i=1, 2, 3, …, 6) and {T i } is the rigid response of the system to a rigid body displacement or infinitesimal rotation in the i-th direction. Actually, {T i } are the columns of [T] defined earlier by Eq. 7 [28]. The MPF is calculated for all extracted modes and in all 6 directions. Modes with higher MPFs are more significant in the global response of the structure. Since the MPF value can be negative or positive, another alwayspositive parameter is defined that enables us to easily compare the significance of the modes. Modal Effective Mass (MEM) is defined as: If MEM is summed for all modes and in all 3 translational directions (i=1, 2, 3) it should be equal to the total mass of the structure. As a rule of thumb, the modes should be extracted until the MEM in all 3 translational directions approaches 90% of the total mass of the structure [29,30]. More modes will enhance the precision but on the other hand, they will significantly increase the time and the hard disk space required for the solution. It is possible to decrease the required number of modes for proper accuracy by using residual modes. These modes are the response of the structure under the application of unit forces and moments exerted on the body. Since the model consists of 72 components of force and moment, it will have 72 residual modes. These modes have high MEMs and will reduce the required number of modes [31]. Hence, in the present study, 423 mode shapes were extracted that lead to 90%, 91% and 90% of MEM to total mass fraction in the X, Y and Z directions, respectively. It should be noted that many of these modes are local modes and do not contribute much to the total MEM.

FATIGUE ANALYSIS OF THE VEHICLE BODY
atigue analysis of the model is performed using two different methods. For areas around spot welds, the mesh insensitive structural stress method is used which was proposed by Dong et al. [32,33] and used in many previous studies [34][35][36][37]. In order to use this method, the equivalent structural stress parameter should be calculated first. To do so, the nodal forces and moments on the weld line around the spot welds are directly calculated by the finite element solver. In the next step, these global forces and moments are transferred into the local coordinate system (x′-y′) defined at each node on the weld periphery lines (l i ) as depicted in Fig. 4(a). This is because the structural stresses are defined as the stress components normal to the spot weld line. The nodal forces in the local coordinate system (i.e., F 1 and F2) are then converted to the linearly distributed force f(x′) as shown in Fig. 4(b). In doing so, it is assumed that the F work done by the nodal forces is equal to the work done by the linearly distributed force. Note that f 1 and f 2 are the values of the distributed force at the nodes. The same procedure is done for the nodal moments. The details of the calculations can be found in Kang et al [38].
where f y′ is the line force in the direction of y′, m x′ is the line moment in the direction of x′ and t is the sheet thickness. In the next step, the equivalent structural stress parameter Ss is calculated using Eq. 24 as follows [38], where r is the bending ratio defined by Eq. 25 and m is an exponent found by experiments and is equal to 3.6 in this case [39]. The function I(r) is calculated according to the fracture mechanics approaches and its diagram for the spot welds can be found in the literature [38].
Now, the calculated structural stresses can be used to assess the fatigue life of the spot welds. To do so, the S s -N curve which is known as the master S-N curve is required. Experimental data on more than 800 steel specimens with different weld types and loadings have shown that utilizing the structural stress to define the S-N curve yields a single line with small scatter of data points around this line [40]. This line can be taken as the master S-N curve for the steel specimens. This is the main idea behind the structural stress method which makes it applicable for all weld types, all materials in the same class (e.g. steels) and all loading types. The master S-N curve is defined by Eq. 26 as follows, where S s is the equivalent structural stress range and coefficients C and h are defined in Tab. 2 for steels. Also, is the variance parameter. As it was explained in the earlier paragraph, experimental data are scattered around a single S-N curve in a narrow band. The variance parameter can be used to take this scatter of data into account. -2 and -3 yield very conservative S-N curves (i.e., lower life for a specific amount of stress) while +2 and +3 are the opposite. The mean curve is used for the analysis in the subsequent sections.  For areas away from spot welds, a multiaxial fatigue criterion is used that consists of the following steps: Step 1: At each material point, the stress tensor is calculated.
Step 2: At each material point, 18 planes which are separated by 10° angles are chosen. Considering the Mohr's circle, all stress states at every material point is now taken into account.
Step 3: Normal stress time histories on the above-mentioned planes are calculated.
Step 4: Using the Rainflow-counting method, the stress cycles are counted on each plane.
Step 5: Fatigue life is calculated on each plane based on the well-known stress-life method.
Step  Table 3: High cycle fatigue and monotonic properties of SAE 1006 steel [41] In order to enhance the precision of the solution, the effect of surface roughness is considered in the fatigue analysis. The arithmetical mean deviation of the assessed surface profile is a parameter that is used to characterize the amount of surface roughness. It is calculated as the mean value of all deviations from a straight line within the evaluation length as shown by Eq. 28, 1 1 n a i i R y n    (28) where n is the number of data points and yi is the amount of deviation from a specific reference line. The average value for the surface roughness of the sheets used in the manufacturing company is found to be R  m. The effect of surface roughness on the fatigue life is characterized by the fatigue stress concentration factor kf. There are many sources in the literature that provide diagrams of k f as a function of R  and material ultimate strength S ut .

Validation of results via a low-frequency loading
s stated earlier, if the loading frequency is less than 10% of the first nonzero natural frequency of the structure, the inertia effects are negligible and therefore, the simulation results for both inertia relief and modal dynamic methods should be identical. Considering the fact that the simulation procedures and the finite element solvers used in these two methods are completely different, if the simulation results for a low-frequency load case are similar, one can make sure that the methods are implemented correctly. In order to calculate the force and moment time histories, first the body is put under the 2g gravitational acceleration as shown in Fig. 5 and all required joint reactions are calculated. Then, a 2 Hz sinusoidal amplitude is applied to the forces and moments acting on the body with a 0.1 sec time delay between the front and rear axles. It should be noted that 2g vertical acceleration is a common load case that is used in some previous studies to assess the strength and durability of the vehicle body and components [42,43].

Durability results for different maneuvers
Previous studies have shown that most of the cracks on a vehicle body are initiated at the areas around the connections like spot welds [44]. Therefore, the authors attempted to present the fatigue results by considering the number of spot welds that fail before 100000 km of service in any of the simulated maneuvers. Note that the failure of a spot weld is characterized by the failure of a single element in any of its two immediately surrounding element layers as depicted in Fig.  3. Fig. 7 shows the location of failed spot welds on the vehicle body for all maneuvers. Also, the number of failed spot welds is depicted in Fig. 8. A Figure 7: Location of the failed spot welds for different maneuvers and two simulation methods.
As it is evident by the simulation results, the relative error for the number of failed spot welds calculated via the quasistatic method is 80, 88, and 96% for maneuvers 1, 2, and 3, respectively. This is a significant error and the quasi-static method of fatigue analysis should be used with great caution. Tab. 4 shows the computation time for both approaches. As it is clear, the time required for the modal dynamic method is almost 37 times higher than that of the inertia relief approach (the workstation used has an 11 core, 3.6 GHz processor with a 64 GB RAM).

CONCLUSION
he aim of this paper is to compare the inertia relief and modal dynamic approaches in terms of formulation, accuracy and computation time. The chosen case study is the fatigue life of the vehicle body structure. To calculate the loadings on the body structure, multi-body dynamics model of the vehicle is created and driven on different standardized roads and by different velocities. The calculated forces and moments are then applied on the BIW and the stress time history on the body is calculated by two different methods which are modal dynamic and inertia relief. These stresses are then used to calculate the fatigue life of the vehicle spot welds by using the structural stress method. The main findings of the present research are as follows:  The inertia relief method is almost 37 times faster than the modal dynamic method.  The modal dynamic method is 96% more accurate than the inertia relief method.  For low frequency loadings, the inertia effects are negligible and the simulations show that there are no major differences between the results of both methods.