Fatigue Life Uncertainty Quantification of Front Suspension Lower Control Arm Design

Abstract

The purpose of this study is to investigate the uncertainty of the design variables of a front suspension lower control arm under fatigue-loading circumstances to estimate a reliable and robust product. This study offers a method for systematic uncertainty quantification (UQ), and the following steps were taken to achieve this: First, a finite element model was built to predict the fatigue life of the control arm under bump-loading conditions. Second, a sensitivity scheme, based on one of the global analyses, was developed to identify the model’s most and least significant design input variables. Third, physics-based and data-driven uncertainty quantification schemes were employed to quantify the model’s input parameter uncertainties via a Monte Carlo simulation. The simulations were conducted using 10,000 samples of material properties and geometrical uncertainty variables, with the coefficients of variation ranging from 1 to 3%. Finally, the confidence interval results show a deviation of about 21.74% from the mean (the baseline). As a result, by applying systematic UQ, a more reliable and robust automobile suspension control arm can be designed during the early stages of design to produce a more efficient and better approximation of fatigue life under uncertain conditions.

1. Introduction

Over the past few decades, structural optimization has become increasingly significant in the automobile sector. The suspension system is one of the car’s most crucial subsystems; it is responsible for protecting the vehicle’s chassis and other components from road shocks, as well as for ensuring stability, comfort, and safety. The suspension system supports the weight of the entire vehicle. One of the independent suspension system designs used in vehicles is the double wishbone suspension system, which employs lower and upper control arms. Forces are transferred from the automobile wheels to the control arm via the ball joint assembly to the wheel as the vehicle crosses over bumps, speed bumps, and other obstacles. Control arms can bend or break when driving over major potholes or can become worn out due to brushing, indicating that the multiaxial loadings associated with many sources of uncertainty are applied to the control arms.

Understanding and managing uncertainty has become increasingly important in the automotive sector to progress the design process, optimize manufacturing procedures, and improve day-to-day technical operations. Areas where uncertainty exists include data collection, processing and interpretation, system modeling, production processes, and operational settings. On the other hand, modern systems require more sophisticated and critical designs that must operate with great dependability, tight safety margins, and high performance. Design approaches for analyzing uncertainty have become more popular as a result and are now being applied in the context of interdisciplinary design [1,2].

The stimulation frequencies of road bumps or speed breakers can cause fatigue damage and dynamic amplification to the suspension’s control arm structure. To ensure service life, it is crucial to take into account the uncertainty in models and parameters. Fatigue damage calculations contain several parameters and combinations of conditions to account for realistic operating scenarios. Although accounting for the variability of some characteristics, such as material qualities and model geometry, is standard practice in the industry, other parameters are typically treated as constants despite data suggesting they may fluctuate dramatically [3]. To understand their significance in terms of the structural response and dependability of the suspension control arm, it is necessary to study these variables. Therefore, the sensitivity of the variables must be addressed before the uncertainty quantification analysis.

The need to identify the most influential factors in the design of support structures, such as suspension control arms, has made sensitivity analysis a popular area of interest in mechanical components. However, determining such values can be challenging due to the complexity of structural design, which includes design requirements, intricate nonlinear models, and numerous input model parameters, all of which significantly increase the uncertainty of the sensitivity results. Researchers from various fields recommend using global sensitivity analysis to account for model uncertainties and nonlinearities in the sensitivity results [4]. Global sensitivity analysis examines the individual and combined effects of the uncertainty of each input parameter on the total variance of the response to determine the most relevant input model parameters. Consequently, a comprehensive examination of the design space that considers interactions between input model parameters can be carried out, providing more reliable sensitivity indices [5]. Global sensitivity analysis has been widely used in the field of structural design due to its advantages. However, previous studies have not provided clear conclusions on which factors to either discard or prioritize (known as model reduction and prioritizing) [6,7] due to the inconsistent selection of input model parameters and the lack of unambiguous convergence metrics. To address this problem, it has been suggested to use multiple global sensitivity analysis methods. For example, rank regression and the Morris methodology, also known as types of elementary effects approaches, have been used with good agreement [4]. However, these methods are not robust enough to identify the most crucial characteristics, and they are limited as researchers will use them before adopting more advanced techniques [5]. Sobol indices, a type of analysis of variance approach, have been successfully used in models with uncorrelated inputs, and are frequently used as a baseline for alternative sensitivity approaches [4,8]. They are able to distinguish between influential and non-influential characteristics, as well as the impacts of their interactions. However, they are computationally expensive and cannot be used for complex models with numerous input parameters as their accuracy highly depends on the number of simulations needed to reach convergence [9].

In order to improve the effectiveness of Sobol’s approach, additional sampling methods and metamodels (surrogate models) have been suggested. Furthermore, the use of metamodels in conjunction with Monte Carlo (MC) simulation has been shown to produce more accurate results than those found in their sampling counterparts [4]. Even though metamodels have been extensively used to quantify uncertainty in mechanical component structures, they have not yet been used in combination with Sobol’s global sensitivity analysis to assist in the design of suspension control arms under fatigue stresses. Therefore, it is necessary to develop sensitivity schemes that enhance the effectiveness and decision making of current sensitivity analyses in the design for fatigue resistance.

Therefore, it is crucial to compute how uncertainties affect model predictions. This study suggested an uncertainty propagation and sensitivity analyses platform in order to measure the variability in fatigue responses of the front-lower suspension control arm. For the sensitivity analysis of the suspension control arm, the Sobol indices approach was suggested, and a Morris one-step-at-a-time (MOAT) parameter screening method [10] was applied for comparison. Surrogate models were compared to replace the computationally expensive finite-element analysis in COMSOL Multiphysics version 6.1 by generating a training sample data using Latin hypercubes (LHCs). A surrogate-based Monte Carlo simulation was then performed to determine the suspension’s lower control arm fatigue life uncertainty. Finally, confidence intervals were computed to quantify the uncertainty in the fatigue life of the lower control arm of the automobile suspension.

The remaining parts of this study are organized as follows: Section 2 explains the details of uncertainty in the fatigue for suspension control arms. Section 3 presents the study methodology, which includes the computational method, fatigue life prediction model, sensitivity analysis scheme, uncertainty propagation model, and surrogate model formulation. The results and discussion are provided in Section 4. Finally, Section 5 provides the study’s conclusion.

2. Uncertainties in Fatigue

Fatigue failure is a common problem in engineering practice, and it can occur randomly. It is particularly challenging to ensure the fatigue resistance of engineering components under specific loading regimes since service conditions and material qualities often fluctuate [11,12,13,14]. Conventional strain-life or stress-life curves that represent only the average response of interest have been shown to be inadequate by [15]. Therefore, probabilistic approaches that take these uncertainties into account are used in practical engineering. In this analysis, probabilistic fatigue modeling is carried out by integrating FE analysis with the Latin hypercube sampling (LHS) technique, which considers the material properties of variability and geometrical uncertainty as capable of producing statistical properties.

The structural engineer is particularly interested in a material’s modulus of elasticity since its value is essential for instability failure scenarios. Even in research that deals with mechanical property variability, Young’s modulus is frequently considered deterministically. However, Galambos and Ravindra [16] addressed the diversity in elastic modulus values and presented important research in this field. For both the tension and compression elastic modulus behavior, they proposed a coefficient of variation (COV) of 0.06 based on the available data at the time. Mansour et al. [17] also researched the variability of Young’s modulus and proposed a COV of 0.031. However, they agreed that the distribution should be normal. Therefore, the COV for this study was set to 0.03 in order to examine fatigue uncertainty within the smaller variability of Young’s modulus under the assumption of a normal distribution; the same COV was also applied to Poisson’s ratio.

Furthermore, the geometrical tolerances between actual and design dimensions are unavoidable due to manufacturing processes or design margins [18,19]. This geometrical tolerance, also known as dimensional tolerance, must be controlled within an authorized range to maintain the interchangeability and assembly of engineering components. As Hu et al. [20] noted, even a small change in dimensional parameters can have a significant impact on how real components respond to stress and strain, thus affecting how well they withstand fatigue.

This study considers geometrical uncertainty when determining the sizes of the suspension control arm. The sampling of dimensions was done in accordance with design margins or manufacturing faults by establishing the tolerance range of each measurement. Following that, the sampled dimensions were used for geometric modeling, and model updates were carried out using finite element (FE) modeling in COMSOL. A few common control arm dimensions were chosen as the geometrical uncertainty inputs from the design and deformation failure analysis perspectives, as shown in Figure 1. In engineering practice, a normal distribution is typically used to model dimensional tolerance [21]. According to DIN EN ISO 286-1 [22], the geometrical tolerances of these dimensions are considered to be ±0.1% changes from the planned size based on the tolerance requirements. Moreover, for this study, the COV for each dimension was considered to be 0.1%, as shown in

Table 1. Uncertain variables of the front lower suspension control arm.

No.	Parameters	Mean Value	Unit	Distribution	COV. %
1	Plate width (w1)	70	mm	Normal	1
2	Cross plate width (w2)	52.5	mm	Normal	1
3	Thickness (ts)	70	mm	Normal	1
4	Supporting plate thickness (tp)	24.5	mm	Normal	1
5	Pin radius (r1)	21.875	mm	Normal	1
6	Bolt 1 radius (r2)	21.875	mm	Normal	1
7	Bolt 2 radius (r3)	17.5	mm	Normal	1
8	Young’s Modulus (E1)	200	GPa	Normal	3
9	Poisson’s ratio (nu1)	0.3	-	Normal	3
10	Fillet_1 (f1)	110.25	mm	Normal	1
11	Fillet_2 (f2)	175	mm	Normal	1
12	Fillet_3 (f3)	21.875	mm	Normal	1
13	Fillet_4 (f4)	110.25	mm	Normal	1
14	Fillet_6 (f6)	17.5	mm	Normal	1
15	Fillet_7 (f7)	8.75	mm	Normal	1
16	Fillet_8 (f8)	11.2	mm	Normal	1
17	Fillet_9 (f9)	8.05	mm	Normal	1
18	Fillet_10 (f10)	8.05	mm	Normal	1
19	Fillet_14 (f14)	11.2	mm	Normal	1
20	Fillet_15 (f15)	11.2	mm	Normal	1
21	Fillet_17 (f17)	5.25	mm	Normal	1
22	Fillet_18 (f18)	4.9	mm	Normal	1

.

3. Methodology

In this study, a computational model and an uncertainty quantification framework were developed to investigate the uncertainty of the design variables of a front suspension lower control arm under fatigue-loading conditions. The principal methodology is depicted in Figure 2. The remainder of this section is allocated as follows: computation modeling, sensitivity analysis, probabilistic modeling, and specification of the analysis. These are the main proposed methods of this study.

3.1. Computational Modeling of the Control Arm

The stresses and deformations caused by mechanical loads can be computed using structural analysis, which helps ascertain how loads affect the physical structure. Since the lower control arm is connected to the wheel and chassis of a car, it is subject to several loads that cause deformation. Failure of the lower control arm is caused by deformation and cyclic stresses, so it is necessary to determine the stresses acting on the lower arm before it is manufactured. Once the stresses for various loads have been computed using structural analysis, fatigue analysis can be performed based on the stress value after the maximum and minimum stresses have been determined.

For this study, the suspension lower control arm model was prepared in COMSOL, and finite element analysis was performed using COMSOL 6.1 software with the structural mechanics module. Figure 3 illustrates the arrangement of the finite elements used for the front suspension lower control arm in this study. In the design process, the appropriate design variables for the parts involves considering the stress distribution (von Mises stress value). In static analysis, the maximum von Mises stress was considered, and fatigue analysis was conducted via the stress-life approach. To achieve the design objective, it is essential to understand the loads assigned to the suspension control arm. For example, as shown in Figure 3a, the forces acting on a vehicle’s tires can be measured in the x, y, and z directions, where x signifies the longitudinal force, y represents the lateral force, and z is the vertical force. The longitudinal force is caused by rolling resistance and traction–compression cycles, which are caused by braking, while the lateral force is produced by camber and toe angles. In this study, the significant force is the vertical force, which is initiated due to the vehicle’s weight and bumping cycles caused by the road bump.

The model is highly dependent on loading and boundary conditions. As previously noted, this study considers a vertical force to be a constant amplitude imposed at the bushing that connects to the tire, while the other two bushings that connect to the body of the car are considered constrained. A tetrahedral mesh was used to approximate the surface contour more accurately, as seen in Figure 3b. To reduce the computational time and to ensure that the FE analysis results were not affected by the mesh size changes, a mesh convergence study was conducted, and the least number of meshes with the highest accuracy was selected.

The mesh in COMSOL Multiphysics always defaults to a physics-controlled mesh with a specified parameter set when you create a model. The predetermined element size parameters are as follows: maximum element size (maxE), minimum element size (minE), maximum element growth rate (maxGR), curvature factor (cf), and the resolution of narrow regions (rnr). The mesh convergency was conducted in a range of extra coarse and extra fine mesh sizes. After the convergency test, the maximum von Mises stress values did not change significantly after the normal mesh size configuration; as such, this study used a fine mesh size (maxE = 8.3 mm, minE = 1.035 mm, maxGR = 1.45, cf = 0.6, and rnr = 0.5) near the necks of the two bushings since these areas are critical for fatigue, and the other part was meshed using the normal mesh size setting (maxE = 10.35 mm, minE = 1.865 mm, maxGR = 1.5, cf = 0.5, and rnr = 0.6) to reduce the computational time.

3.2. Fatigue Life Prediction

Low-cycle fatigue (LCF) and high-cycle fatigue (HCF) are two distinct phases of fatigue characteristics [23]. LCF is typically referred to as “strain-life/strain-based” because strain rather than stress is the most significant descriptive measure of LCF. However, strain is also applicable for damage in HCF, but since HCF occurs in the elastic phase, both stress and strain can be used as the factor. HCF is also regarded as “stress-life/stress-based” since stress is used for practical and historical purposes. The number of cycles determined in this case (suspension control arm) corresponds to an HCF paradigm, which is typically defined as a number of cycles counts needed until reaching a fatigue of more than 10⁴ [24,25].

The type of load pattern applied to a structure under constant load cycles can also have an impact on its performance. The load can be classified as either proportional or non-proportional. In proportional loading, the primary loads and strains remain oriented throughout the load cycle. Another method to distinguish between these two scenarios is to consider the characteristics of the external load. The structural response to a single external load source is represented by a stress tensor, whose components undergo phase changes. If the external load is applied multiple times or, if it is a traveling load, the components of the stress tensor may shift out of phase. Different methods of fatigue evaluation are required for these two distinct types of load cycles [26,27,28,29].

In this research, it was assumed that the load was sinusoidal, which denotes a proportional load. Stress-life models can assess fatigue under proportional loading via the S-N curve model (Wöhler curve model) [30] or Basquin model [31,32]. The S-N curve model (Wöhler curve model) was employed in this study since it is one of the oldest and simplest models for fatigue prediction, and it demonstrates a clear relationship between fatigue life and applied stress. The fatigue properties of the suspension control arm material are summarized in Figure 4. This material has an endurance limit of 110 MPa.

Considering that the suspension arm may have a rough surface due to poor cleanliness or fairly large inclusions, this could affect the value of the endurance limit. For instance, the large inclusions can shorten the lifetime. To ensure safety, the S-N curve data should be modified, and one possible modification is to adjust the amplitude stress as follows [33].

$${\sigma }_a=k·f_{SN}\left(N\right),$$

(1)

where ${\sigma }_a$ is the stress amplitude, which is determined from the difference between the peak stresses and the mean stress $\left(\right({\sigma }_{max}-{\sigma }_{min})/2)$; $N$ is the number of the fatigue lifetime cycle; $k$ is the modification factor; and $f_{SN}$ is the S-N curve. Based on past experience, the modification factor caused by the production process of machining can be set between 0.6 and 0.8 [34].

In fatigue evaluation when using the S-N curve model, the amplitude stress can be modified to account for the mean stress. The mean stress ${\sigma }_m$ is defined as the average of the minimum and maximum stresses during a load cycle. There are three kinds of mean stress correction methods: as per Gerber [35], Goodman [36], and Soderberg [37]. The Gerber and Goodman methods predict failure when the mean stress is at the ultimate tensile stress, whereas the Soderberg method predicts failure when the mean stress is at the yield limit. This study used a yield limit, so we used the Soderberg method. Based on the Soderberg method, the amplitude stress is modified using the mean stress as

$${\sigma }_a\leftarrow {\sigma }_a\left(1-\frac{{\sigma }_m}{{\sigma }_y}\right),$$

(2)

where ${\sigma }_y$ is a yield stress that you need to specify. Figure 5 presents COMSOL’s simulation results for fatigue life, with the mean values of the design variables listed in Table 1. The results indicate that the minimum life cycle occurs around the necks of the two bushings.

3.3. Sensitivity Analysis Model

Sensitivity analysis can identify the parameter or set of parameters that has the greatest impact on the model’s response. As a result, it provides valuable insight into which model input has the most influence on the unpredictability of the model’s response [38]. It is also a useful tool for model developers and users to verify the model’s uncertainty around input parameters, and to provide feedback for model refinement to increase confidence. The results of the sensitivity analysis can help model builders focus on the essential variables that influence the model’s response, especially in the case of a very complex model. As mentioned in the introduction, there are two types of sensitivity analysis: local and global sensitivity analysis. Local sensitivity analysis evaluates variations in the model’s response with respect to variations in a single parameter input [39,40]. In global sensitivity analysis, all parameters are altered simultaneously across the entire parameter space, allowing for the evaluation of both the relative contributions of each individual parameter and the interactions between parameters with respect to the variance of the model’s response.

A variety of global sensitivity methods are available for systems application models, including multiparametric sensitivity analysis, Fourier amplitude sensitivity analysis, the partial rank correlation coefficient, and Sobol’s approach. Zhang et al. [41] attempted to compile the benefits and drawbacks of these approaches, and they concluded that the variance decomposition-based Sobol sensitivity analysis is currently one of the most effective approaches for conducting global sensitivity analysis.

The Sobol method examines the complete input parameter distribution and decomposes the variance of the response into contributions from the input parameters and their interactions. The effects of the parameters are calculated based on the model assessment data. The number of terms in the Sobol indices of the variance of the response with $m$ input parameters grows as $2^m$. To manage this, it is customary to compute only the $m$ first-order effects (first-order Sobol indices) and the $m$ total effects (total Sobol indices). The first-order index shows the contribution of the parameter to the response variance without interaction. The total effect indices show the overall contribution of a parameter to the response variance, including interaction.

Given a scalar response that is defined as $y=M(x_1,x_2,\dots ,x_m)$, the variance-based first-order Sobol index is defined as [42]

$$S_i=\frac{{Var}_{x_i}\left(E_{{\mathit{x}}_{\backslash i}}\right(y\left|x_i\right))}{Var\left(y\right)},$$

(3)

where $x_i$ indicates the $i$th input parameter and ${\mathit{x}}_{\backslash i}$ represents all parameters except $x_i$. The inner expectation value means that the mean of $y$ is calculated over all possible values of ${\mathit{x}}_{\backslash i}$ while $x_i$ is fixed. Then, the outer variance is calculated over all possible $x_i$ values.

The total Sobol index is defined as [42]

$$S_{T,i}=1-\frac{{Var}_{{\mathit{x}}_{\backslash i}}\left(E_{x_i}\left(y|{\mathit{x}}_{\backslash i}\right)\right)}{Var\left(y\right)},$$

(4)

which measures the summation of the first-order effect and the interactions with other parameters of the $i$th input parameter. The second term in the equation can be seen as the first-order effect of all parameters except for the $i$th parameter. Therefore, subtracting the second term from one yields the contribution of all terms related to $x_i$.

In this study, the Sobol indices were computed for the fatigue response of the front lower suspension control arm. With the help of the uncertainty quantification module in COMSOL, there are two different ways through which to compute the Sobol indices. One method is to use the post-processing of a polynomial chaos expansion (PCE) model. The contributions from each input parameter and their interactions can be easily separated given the PCE’s specification. As a result, the Sobol indices can be calculated using the PCE-trained coefficients. In an earlier study from Blatman and Sudret [43], more information on the computation of Sobol indices based on the PCE model was provided. The alternative method for calculating the Sobol indices is the Monte Carlo approach. In this kind of study, a Gaussian process is provided as a surrogate model because performing Monte Carlo simulation directly with model evaluation is computationally expensive. The approach adheres to the recommended method of computing the first and total Sobol indices simultaneously, which was also discussed in Saltelli et al. [5].

3.4. Uncertainty Propagation

Studying the uncertainty propagation of a model is equivalent to approximating the probability density function (PDF) of the response. To construct an accurate PDF of the response, a large number of samples are required. However, this operation is very expensive, so the analysis uses a surrogate model instead. To approximate the PDF, the kernel density estimation (KDE) method is used. This method formulates density in terms of the known kernel functions, which have the following form [44]:

$$\widehat{f}\left(x\right)=\frac{1}{nh}\sum _{i=1}^n{K}_h\left(\frac{x-x_i}{h}\right),$$

(5)

where $K_h$ is a specified symmetric PDF. This study assumed that the kernel function is given as the standard normal distribution, and $h$ is a smoothing parameter, which is termed the bandwidth. Silverman’s rule was also applied for the smoothing parameter $h$. More details related to the KDE can be found in Sheather [44].

To evaluate the uncertainty propagation, a confidence interval was used. A confidence interval gives us an indication of the degree of uncertainty in an estimated probability distribution. In statistics, a probability distribution describes all the possible values that a random variable can take within a range. This range is bound to be within certain minimum and maximum possible values, but precisely where the possible value is likely to be plotted on the probability distribution depends on a number of factors. These factors include the distribution’s mean, standard deviation, and skewness. The most common expression for the mean of a statistical distribution with a discrete random variable is the mathematical average of all the terms. To calculate it, add up the values of all the terms and then divide them by the number of terms. The mean of a statistical distribution with a continuous random variable, also called the expected value, is obtained by integrating the product of the variable with its probability, as defined by the distribution. A greater standard error will result in a wider interval; the wider the interval, the less precise the estimate is. Detailed information on confidence intervals can be found in Dekking et al. [45].

In this study, uncertainty propagation was performed by conducting a number of model evaluations while using Latin hypercube sampling data with the selected parameters after conducting sensitivity analysis. These data were used to create surrogate models by using PCE and the GP model in COMSOL. Next, a Monte Carlo analysis was performed using the surrogate model to make predictions. A KDE was generated with the sampled data to produce an estimation of the probability density for each response. The suggested surrogate models for both sensitivity analysis and uncertainty propagation in this study will be discussed in the next section.

3.5. Surrogate Models

Surrogate models are mathematical models that estimate the computationally expensive simulation models, and they are based on the desired design of experiment samples. In this section, the two most commonly used approaches for surrogate modeling are explained along with their adaptive methods, which are also used for sensitivity analysis and uncertainty propagation analysis in this study.

3.5.1. Polynomial Chaos Expansion

To reduce the computational time, a surrogate model based on PCE was extensively implemented in the statistical analysis of various engineering applications. The main idea behind PCE is to provide a polynomial surrogate for the computational model. In this regard, the polynomial chaos expansion that establishes a relationship between the system response ($Y$) and the independent input parameters of the n-dimensional $\mathit{x}=\{x_1,x_2,...,x_n\}$ is given as follows [46]:

$$Y=\sum c_{\alpha }{\mathsf{\Psi }}_{\alpha }\left(\mathit{x}\right),$$

(6)

where $\alpha$ is a multidimensional index that indicates the degree of the polynomial in each of the input variables, ${\mathsf{\Psi }}_{\alpha }$ denotes the multivariate polynomials, and $c_{\alpha }$ denotes the unknown coefficients to be determined. The multivariate polynomials ${\mathsf{\Psi }}_{\alpha }$ are formed by the tensor products of the univariate orthogonal polynomials, which are as follows:

$${\mathsf{\Psi }}_{\alpha }=\prod _{i=1}^n{\phi }_{\alpha }^{\left(i\right)}\left({\mathit{x}}_i\right),$$

(7)

where ${\phi }_{\alpha }^{\left(i\right)}$ is the univariate orthogonal polynomial in the ith variable of degree α. The probability distribution of the input parameters is taken into consideration when selecting the univariate orthogonal polynomials. Some commonly used classical univariate polynomials with orthonormal distributions can be obtained in [47]. For example, if the input data distribution is uniform or normal, the suggested univariate polynomials are Legendre and Hermite, respectively. If the input parameters do not have a distribution, such as that which is suggested in Xiu and Karniadakis [48], the algorithm defines an isoprobabilistic transformation to a uniform distribution, and then constructs the PCE using the transformed input parameter.

Based on Equation (5), PCE provides an infinite polynomial series, which means it needs to be truncated to a finite sum for computational purposes, resulting in a truncation error. According to the definition of PCE, truncation is defined as the total degree of all polynomials being less than or equal to $p$. Additionally, the polynomial expansions are further truncated using a hyperbolic truncation method that utilizes a parametric q-norm, $q\mathrm{ }[0,1]$. The truncation method truncates all polynomials with total degrees that are less than $p$ when $q=1$. For $q<1$, hyperbolic truncation retains all high-degree terms in each input parameter but discards the equivalent high-order interaction terms. The truncation method truncates all the interaction terms between the input parameters when $q=0$.

A sparse representation of PCE refers to a sparsely selected polynomial basis from all the truncated polynomial expansions. An adaptive algorithm based on the least-angle regression method is used to select the significant sparse coefficients from all the truncated polynomial expansions. During the iterative process, the algorithm only adds new polynomial bases that are mostly correlated with the residuals that are built using the existing polynomial basis. The sparse representation is used for reducing the computational cost and, more importantly, for avoiding overfitting. Here, the leave-one-out cross-validation error estimation is defined with a correction factor that considers overfitting to ensure that the generalization error estimate is not underestimated. The generalization error refers to the error for the sample points with unknown model evaluations. Additional details on the sparse PCE representation and error estimation can be found in [43,48].

3.5.2. Gaussian Process

The Gaussian process, commonly referred to as the Kriging model, is one of the most popular surrogate models that provide an estimation of the responses. A Gaussian process model is a probabilistic model that specifies a distribution over functions and evaluates the variance of the estimation at each sample point in the input parameter space. The Gaussian process’s mean, often referred to as the trend, and covariance function, also known as a kernel or correlation function, are used to define it. While the Gaussian process model is being trained, the mean and covariance hyperparameters are tuned. Consequently, the Gaussian process model may be written as follows (considering $\mathit{x}$ as the input parameters and the responses as vector y):

$$y=m\left(\mathit{x}\right)+f\left(x\right),$$

(8)

where f(x) denotes a Gaussian process with zero mean and covariance $k_{\theta }\left(\mathit{x},{\mathit{x}}_{✶}\right)$, and m(x) is the mean function. Here, $\theta$ denotes the covariance function’s hyperparameters. The joint distribution of the calculated data and the estimated data can then be expressed as [49]

$$\left[\begin{array}{c}f\\ f_{✶}\end{array}\right]=N\left(0,\left[\begin{array}{cc}{k}_{\theta }\left(\mathit{x},\mathit{x}\right)& k_{\theta }\left(\mathit{x},{\mathit{x}}_{✶}\right)\\ k_{\theta }\left({\mathit{x}}_{✶},\mathit{x}\right)& k_{\theta }\left({\mathit{x}}_{✶},{\mathit{x}}_{✶}\right)\end{array}\right]\right),$$

(9)

where $f=y-m\left(\mathit{x}\right)$, $y$ is the calculated response, and ${\mathit{x}}_{✶}$ and $f_{✶}$ are the testing input parameters and estimated response, respectively.

The posterior distribution of the Gaussian process is given as [49]

$$y\left({\mathit{x}}^{\mathbf{\prime }}\right)~N\left(\mu \left({\mathit{x}}_{✶}\right),cov\left(\mathit{x},{\mathit{x}}_{✶}\right)\right),$$

(10)

where $cov(\mathit{x},{\mathit{x}}_{✶})$ and $\mu \left(\mathit{x}\right)$ are the covariance and mean of the model, respectively:

$$\begin{gathered} \mu \left({\mathit{x}}_{✶}\right)=m\left({\mathit{x}}_{✶}\right)+k_{\theta }({\mathit{x}}_{✶},\mathit{x})k_{\theta }{\left(\mathit{x},\mathit{x}\right)}^{-1}f\left(\mathit{x}\right) \\ cov\left({\mathit{x}}_{✶},\mathit{x}\right)=k_{\theta }\left({\mathit{x}}_{✶},{\mathit{x}}_{✶}\right)-k_{\theta }\left({\mathit{x}}_{✶},\mathit{x}\right)k_{\theta }{\left(\mathit{x},\mathit{x}\right)}^{-1}{k}_{\theta }\left(\mathit{x},{\mathit{x}}_{✶}\right) \end{gathered}$$

(11)

Covariance functions are an essential component of Gaussian processes that control the prior and posterior shapes of the model. They encode prior domain knowledge about $f$. These functions, instinctively, allow the model to be generalized by correlating new input data (${\mathit{x}}^{\prime }$) to the existing observation data ($\mathit{x}$). COMSOL provides three commonly used kernel functions: the spectral exponential kernel, and the Marten kernel, where $v=5/2$ and $v=3/2$. The spectral exponential kernel function is given as [49]

$$k_{SE}\left(\mathit{x},{\mathit{x}}^{\prime }\right)=\mathrm{e}\mathrm{x}\mathrm{p}\left(-\sum _{d=1}^m\frac{{(x_d-x_d^{\prime })}^2}{{2l}_d^2}\right),$$

(12)

where each $l_d$ is the length scale parameter for dimension $d$. The Marten kernel also functions when given as [50]

$$k_{M_{\frac{3}{2}}}\left(r\right)=\left(1+\sqrt{3}\mathrm{r}\right)\exp \left(-\sqrt{3}\mathrm{r}\right), \mathrm{a}\mathrm{n}\mathrm{d} k_{M_{\frac{5}{2}}}\left(r\right)=\left(1+\sqrt{5}\mathrm{r}+\frac{5}{3}{r}^2\right)\exp \left(-\sqrt{5}\mathrm{r}\right),$$

(13)

where $r$ is given as function of $\mathit{x}$ and ${\mathit{x}}^{\prime }$:

$$r\left({\mathit{x},\mathit{x}}^{\prime }\right)=\sqrt{\sum _{d=1}^m\frac{{\left(x_d-x_d^{\prime }\right)}^2}{l_d^2}}.$$

(14)

Further details about the covariance function and the Gaussian process can be found in Williams and Rasmussen [49].

The kernel function $k$ also describes the prior on the types of functions that might be represented in the observed data; for instance, it can convey expectations for smoothness or periodicity. The hyperparameters $\theta$ in the parametric kernels influence this prior and, in turn, the posterior prediction. These kernel hyperparameters can be modified in accordance with a dataset’s characteristics to define a prior value over the functions that are suitable for the situation. To maximize the Gaussian process’s marginal likelihood, $p\left(f\right|\mathit{x})$, the kernel hyperparameters are often trained using the gradient-based optimization of [50]

$$\log p\left(f|\mathit{x}\right)=-\frac{1}{2}f{k}_{\theta }\left(\mathit{x},\mathit{x}\right)f-\frac{1}{2}\log \left|k_{\theta }\left(\mathit{x},\mathit{x}\right)+{\sigma }^{2}I\right|-\frac{n}{2}\log 2\pi ,$$

(15)

where ${\sigma }^{2}I$ is a nugget factor, and here $\sigma$ is an additional hyperparameter, which is used while training the Gaussian process model.

In regard to adapting the Gaussian processes, the Gaussian process is created adaptively by incorporating a new input parameter point at the position of the maximum entropy while it is applied to sensitivity analysis or uncertainty propagation. The Gaussian process is trained to accurately replicate the underlying model in the whole input parameter space for sensitivity analysis or uncertainty propagation. The entropy is expressed here as the standard deviation of the Gaussian process model.

Each adaptation phase involves solving a global optimization problem to determine the position corresponding to the maximum entropy or the maximum expected feasibility function in order to explore the region in the input parameter space and to take advantage of the most recently constructed Gaussian process. There are two different categories of global optimization techniques in COMSOL: the Monte Carlo approach and the dividing rectangles method [51]. The Monte Carlo approach was used in this investigation. It should be noted that adaptive Gaussian process models need adequate initial model assessments for the adaptation approach to find good adaptation points that can raise the model’s accuracy. The unsearched region may not be effectively explored given an initial dataset with too few sample points. Adding extra samples is one technique to solve this issue, so this study applied a sequential optimal Latin hypercube sampling, where the ideal LHS samples randomly cover the entire input parameter space.

4. Results and Discussion

Section 3 describes how uncertainty propagation is used to perform computational fatigue damage modeling and the subsequent design assessment. Model sensitivity analysis is first carried out to eliminate the unnecessary factors from statistical calibration as it would be costly to explore uncertainty propagation by using all design variables. The fatigue analysis data are then used in uncertainty propagation to statistically infer the model’s most important parameters. Finally, the calibrated model with measured uncertainty is used to make predictions about the fatigue strength of the material in the specified loaded stress state.

4.1. Sensitivity Analysis of the Model Input Parameters

4.1.1. Polynomial Chaos Expansion

To begin with,

Table 2. Surrogate model comparison for sensitivity analysis.

No.	Parameters
Space polynomial chaos expansion	0.4512
Adaptive polynomial chaos expansion	0.4350
Gaussian process	0.9645
Adaptive Gaussian process	0.4350

describes the validation results of the developed surrogate model for the fatigue life prediction described in Table 1. While we are developing the surrogate model, 400 instances of sample data were generated for the training data via the Latin hypercube sampling method. From the statistics, the comparison between the four types of surrogate models are shown. The root-mean-square error (RMSE) shows the adaptive sparse polynomial chaos expansion, and it shows less errors than the others; thus, for sensitivity analysis, the adaptive sparse polynomial chaos expansion was chosen in this study to proceed with.

4.1.2. Sensitivity Analysis

The results of the influence of the car front suspension lower control arm input parameters, which were based on the Sobol method (also known as the variance-based sensitivity analysis method scheme for fatigue life), are depicted in Figure 6a,b. As can be seen from the figure, there are two different types of Sobol indices: the first-order index and total index. As mentioned in Section 2, the first-order index of a parameter shows the sensitivity by varying this parameter alone. As shown in Figure 6b, the first-order index result demonstrates that the control arm thickness ($ts$) and width ($w$) individually had a significant influence on fatigue life. The total index shows how much a parameter contributes to the overall sensitivity. As shown in the figure, for all parameters, the first and total indices were equal, which indicates there is very little or no interaction between the parameters. The Sobol indices and MOAT plots also show that the thickness ($ts$) of the control arm presents the higher sensitivities, and this is followed by the control arm plate width ($w$); meanwhile, the Young’s modulus ($E$) had the third largest contribution in the third, and fillt-3 ($f3$) had the fourth largest contribution. As a result, small variations in $ts$ and $w$ can significantly increase/decrease the stress concentration of the suspension control arm, which—in turn—increases/decreases the accumulated structural damage (fatigue life) during the suspension control arm operational life. The parameters $tp$ and $f8$ also presented slight contributions to the response, particularly $tp$. The small variance of the other parameters seems to not have any significant effect on the fatigue life of the control arm. In summary, the most influential parameters, in contrast, were found to be $ts, w, E, f3, tp$, and $f8$; thus, to reduce the variance of the response, it is recommended to collect more data and conduct more sophisticated simulations to reduce the uncertainty of those parameters.

4.2. Uncertainty Propagation

The uncertainty quantification used in the fatigue life cycle is described and the model input PDFs are presented; in addition, the confidence interval for the fatigue life cycle response was computed. First,

Table 3. Chosen parameters for the uncertainty propagation.

No.	Parameters	Mean Value	Unit	Distribution	COV. %
1	Plate width ($w1$)	70	mm	Normal	1
2	Thickness ($ts$)	70	mm	Normal	1
3	Supporting plate thickness ($tp$)	24.5	mm	Normal	1
4	Young’s modulus ($E$)	200	GPa	Normal	3
5	Fillet_3 ($f3$)	21.875	mm	Normal	1
6	Fillet_8 ($f8$)	11.2	mm	Normal	1

presents the selected six parameters ($w1, ts, tp, E, f3$, and $f8$) of the twenty-two used in the sensitivity analysis, along with their probability distributions and coefficients of variation.

Then, the surrogate model was developed by generating 300 instances of sample data while using LHS for the chosen parameters.

Table 4. Surrogate model comparison for uncertainty propagation.

No.	Parameters
Space polynomial chaos expansion	0.3251
Adaptive polynomial chaos expansion	0.3165
Gaussian process	0.3314
Adaptive Gaussian process	0.3027

also describes the validation results of the developed surrogate model for the fatigue life prediction. Based on the RMSE result, the adaptive Gaussian process surrogate model showed less errors than the others. A Monte Carlo analysis was then conducted, using an adaptive Gaussian process surrogate model, to evaluate the uncertainty propagation.

Based on the uncertainty quantification study, a kernel density estimation (KDE) plot is shown in Figure 7, along with the corresponding confidence interval information in

Table 5. Quantity of the interest confidence intervals.

Mean	STD	Min.	Max.	Lower 90%	Upper 90%	Lower 95%	Upper 95%	Lower 99%	Upper 99%
22,902.4	4979.2	7917.9	43,473.0	14,711.7	31,093.1	13,143.3	32,661.6	10,056.2	35,748.6

. In the uncertainty propagation study, kernel density estimation is computed. Alternatively, you can think of the KDE as a smoothed version of a histogram, which provides an estimate of the probability density function of the quantity of interest based on the input parameters. Based on the confidence interval table and plot of the quantity of interest, we can see that the mean number of cycles is around 22,902.4 with a standard deviation of 4979.2 or 21.74% deviation from the mean. This is achieved as a result of the small variability of input parameters, i.e., a large deviation from the mean (baseline).

5. Conclusions

In this study, we developed a surrogate-based fatigue uncertainty quantification framework that combines FE analysis with the LHS technique. The objective was to investigate how geometrical uncertainty and material variability affect the fatigue life for a front suspension lower control arm. We compared the fatigue life and examined the failure probability under these uncertainties, along with a sensitivity analysis. For the sensitivity analysis, we used the Sobol method with an adaptive polynomial chaos expansion surrogate model to rank the influential random variables in the fatigue analysis of the front suspension lower control arm. It was found that, out of twenty-two parameters, five geometrical uncertainties and young modulus variabilities had the highest impact on fatigue life and failure probability. We then performed a Monte Carlo analysis using the adaptive Gaussian process surrogate model to examine the uncertainty propagation according to the variability of the selected six parameters.

The uncertainty effect is expressed in terms of fluctuations from the mean value, probability distributions, and percentage deviations from the mean (baseline). We observed a COV of 21.74% as a result of this parameter uncertainty, indicating that uncertainty in parameters has an impact on fatigue life cycle. In this study, the results are based on COV assumptions. However, the outcome of this research clearly illustrates the need to incorporate the uncertainties in car suspension control arm performance analysis for the design and manufacturing process. In the future, researchers can explore more realistic probability distribution functions other than normal distribution functions for uncertainty analysis, and they can use the current results as a benchmark.