Uncertainty Quantification of Ride Comfort Based on gPC Framework for a Fully Coupled Human–Vehicle Model

Abstract

We investigated the stochastic response of a person sitting in a driving vehicle to quantify the impact of an uncertain parameter important in controlling defect reduction in terms of ride comfort. Using CarSim software and MATLAB/Simulink, we developed a fully coupled model that simulates a driving vehicle combined with an analytical nonlinear human model. Ride comfort was evaluated as a ride index considering the frequency weights defined in BS 6841. Additionally, to investigate the uncertainty of the ride index, a framework for calculating the ride index was proposed using the generalized polynomial (gPC) method. Further, sensitivity analysis of the ride index was performed for each uncertainty parameter, such as stiffness and damping. The results obtained through the gPC method were in good agreement with those obtained via Monte Carlo simulation (MCS) and were excellent in terms of computation time without a loss of numerical accuracy. Through in-depth investigation, we found that the stochastic distribution of the ride index varies differently for each uncertain parameter in the human model. By comparing linear and nonlinear human models, we also found that the nonlinearity of the human model is an important concern in the stochastic estimation of ride comfort.

1. Introduction

The dynamic properties of the human body have received attention in recent years due to the fact that the response of whole-body vibrations is closely related to the ride comfort of occupants in transportation systems. Meanwhile, a quantitative understanding of how a seated human body responds to whole-body vibrations allows automotive engineers to design more human-oriented vehicles. In addition, whole–body vibration has been recently considered in the design and assessment of road infrastructure (paved roads, highways, bridges, etc.) [1,2,3,4,5]. However, despite many efforts to predict ride comfort, the quantification of ride comfort remains challenging as it relies on the perceptions and emotions that are influenced by many factors, such as driving and environmental conditions.

In the automotive industry, both subjective and objective measurements are commonly used to evaluate the ride comfort of a vehicle. In subjective measures, various kinds of experimental studies have been used to correlate with subjective ratings of occupants exposed to specific driving or seating conditions [6,7]. On the other hand, the quantitative analysis of whole-body vibrations is attracting attention in order to objectively investigate ride comfort because this evaluation can detect small variations in the passenger’s dynamic response compared to subjective evaluations [8]. This vibration method is widely used in various applications, such as conventional vehicles [9], railway vehicles [10], high-speed electric multiplexers [11], and helicopters [12]. Furthermore, the spread of electric or hybrid vehicles has led to a new turning point in the study of vibration analysis in the automotive industry. In general, these vehicles do not mask road-induced vibrations in the low-frequency range owing to the absence of a combustion engine [13]. Therefore, understanding human responses to vibrations will become even more important in the near future.

Because experimental approaches for the evaluation of ride comfort are time-consuming and costly, increasing attention is being paid to the industrial requirements for computational analysis based on dynamic simulations. Brogioli et al. [14] developed a model that combines a seat and a human as well as a five-degrees-of-freedom mathematical model, computed weighted accelerations, and a Seat Effective Amplitude Transmissibility (SEAT) index. Additionally, Mohajer et al. [15,16] introduced a biomechanical human model obtained using a computational multibody system and investigated the ride comfort index with respect to road roughness and the longitudinal velocity of a vehicle. Anandan et al. [17] adopted a human–vehicle model with 14 degrees of freedom to analyze the ride comfort of the proposed active suspension vehicle system. Jun Wu et al. [18,19] developed a train–seat–human model with vertical, lateral, and roll vibrations to study the ride comfort of rail vehicles.

In general, the parameters of a real human body (e.g., weight, length, and stiffness) are highly uncertain compared to typical mechanical structures [20,21]. Furthermore, human parameters result in time-variant systems because of the variations in muscle tension, blood flow, and other environmental circumstances. Consequently, the dynamic response of the human body to whole-body vibrations is uncertain due to the uncertainty of the parameters. Therefore, numerical calculations of dynamic human responses to whole-body vibrations using deterministic parameters are limited in the sense that they do not provide a unique subjective rating of feeling.

Numerous stochastic methods have been introduced to calculate the response of a system in the presence of parameters with uncertainty. Among them, Monte Carlo simulation (MCS) is widely used to estimate system uncertainty [22,23,24,25,26]. In MCS, a set of samples of an input parameter with uncertainty are randomly selected from the corresponding probability distributions that are estimated or assumed to be given, and multiple simulations are taken to compute the corresponding output for each set of samples. However, due to the nature of MSC and random sampling, a huge number of simulations are required, which can be computationally time-consuming. In order to reduce the computational cost, many studies have proposed to use the quasi-Monte-Carlo method with an improved sampling algorithm [27,28]. Despite these efforts, stochastic methods based on MCS can be computationally time-consuming as they typically require a large number of simulations. To deal with modeling uncertainty more efficiently, polynomial chaos expansion has become an important method [29,30], which is based on the homogeneous chaos theory of Wiener [31]. The polynomial chaos method has been applied to various engineering fields, such as vehicle dynamics [32,33,34,35,36].

In this study, we analyzed the ride comfort stochastically to quantify the effect of uncertain human parameters. To calculate the ride comfort of a seated person in a driving vehicle on a rough road, we developed a fully coupled human–vehicle model using a vehicle in CarSim combined with a nonlinear human model in MATLAB/Simulink. The ride comfort was quantified as s ride index according to BS 6841 [37], which considers frequency-weighted accelerations. For stochastic calculations, we adopted an approximate model based on gPC (generalized polynomial chaos) and compared it with MSC. In order to investigate the effect of uncertainty, the stochastic distributions of the ride index were calculated for each uncertain parameter of the human model, such as stiffness and damping. Furthermore, we found significant differences between linear and nonlinear human models in terms of uncertain parameters.

2. Human–Vehicle Model and Road Profile for Stochastic Analysis

This section describes a human–vehicle model, road surface generation, and an approximate model based on gPC. The coupling between the human body and the vehicle model is implemented in MATLAB/Simulink. The road profile is generated using the road roughness classification method based on power spectral density. Finally, an approximate model based on generalized polynomial chaos expansion (gPC) is developed for stochastic analysis [29].

2.1. Model Description and Derivation of Equation of Motion

To calculate the responses of the seated human, the nonlinear five-degree-of-freedom human model, which was used in the previous works [38], was considered in this study. The schematic of the proposed human model is illustrated in Figure 1. The lumped parameter of the human model is composed of a head, a trunk, and a thigh with mass and moment of inertia. Each segment is connected using a torsional spring and a damper, and the trunk and thigh contact the seat through a spring and dampers, which are effective parameters combined with seat cushion and human skin. Regarding the generalized coordinates of the hip joint, x_h and z_h indicate the horizontal and vertical displacement of the hip joint, respectively, and θ₁, θ₂, and θ₃ denote the angular displacement of the thigh, trunk, and head, respectively. Only a normal deflection of the springs and dampers is considered for the calculation of the spring and damping forces, respectively (see Figure 2). In fact, the effective stiffness, k_i, connected to the thigh has different values depending on the contact location according to the experimental results [39]. However, due to the fact that the effect of the uncertainty of the design parameters on human–vehicle ride comfort is the main interest of our study, we adopted an identical value for the simplicity of analysis.

The nonlinear equation of the five-degree-of-freedom human model is derived using Lagrange’s equation. Kinetic energy (T), potential energy (V), and Rayleigh’s dissipation (D) of each segment are derived based on the displacement of the center of gravity, as follows:

$$\begin{array}{c}T=\frac{1}{2}{\displaystyle \sum _{i=1}^3}{m}_i\left({\dot{x}}_i^2+z_i^2\right)+\frac{1}{2}{\displaystyle \sum _{i=1}^3}{J}_i{\dot{\theta }}_i^2\\ V=\frac{1}{2}{\displaystyle \sum _{i=1}^4}{k}_i{\delta }_i^2+\frac{1}{2}{\displaystyle \sum _{i=1}^2}{k}_i^T{\left({\theta }_i^T\right)}^2+{\displaystyle \sum _{i=1}^3}{m}_{i}g{∆}_i\\ D=\frac{1}{2}{\displaystyle \sum _{i=1}^4}{c}_i{\dot{\delta }}_i^2+\frac{1}{2}{\displaystyle \sum _{i=1}^2}{c}_i^T{\left({\dot{\theta }}_i^T\right)}^2\end{array}$$

(1)

where m_i and J_i denote the mass and mass moment of inertia of each segment, respectively. Additionally, k_i and c_i are translational, and $k_i^T$ and $c_i^T$ are torsional stiffness and damping coefficients, respectively. The translational and angular displacements of the springs are represented by ${\delta }_i$ and ${\theta }_i^T$, respectively. In addition, $g$ is the gravitational acceleration and ${∆}_i$ is the vertical displacement of each segment. The positions of the center of gravity at each segment are as follows:

$$\begin{array}{c}{x}_1=x_h+\frac{1}{2}{L}_{th}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_1,\\ x_2=x_h+\frac{1}{2}{L}_{tr}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_2,\\ x_3=x_h+L_{tr}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_2+\frac{1}{2}{L}_h\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_3,\\ z_1=z_h+\frac{1}{2}{L}_{th}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_1,\\ z_2=z_h+\frac{1}{2}{L}_{tr}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_2,\\ z_3=z_h+L_{tr}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_2+\frac{1}{2}{L}_h\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_3\end{array}$$

where x_h and z_h are the horizontal and vertical displacement of the hip joint. The length of the thigh, trunk, and head are denoted as L_th, L_tr, and L_h, respectively. Accordingly, the translational and angular displacements of each spring are given below.

$$\begin{array}{c}{\delta }_1=z_h-\frac{1}{2}{T}_{th}\cos {\theta }_1+L_1\sin {\theta }_1-z_b-{\delta }_1^{*}\\ {\delta }_2=z_h-\frac{1}{2}{T}_{th}\cos {\theta }_1+L_2\sin {\theta }_1-z_b-{\delta }_2^{*}\\ {\delta }_3=x_h\sin {\theta }_2^{*}-z_h\cos {\theta }_2^{*}-L_3\sin \left({\theta }_2-{\theta }_2^{*}\right)-\frac{1}{2}{T}_{tr}\cos \left({\theta }_2-{\theta }_2^{*}\right)-x_3^{*}\sin {\theta }_2^{*}+z_3^{*}\cos {\theta }_2^{*}-{\delta }_3^{*}\\ {\delta }_4=x_h\sin {\theta }_2^{*}-z_h\cos {\theta }_2^{*}-L_4\sin \left({\theta }_2-{\theta }_2^{*}\right)-\frac{1}{2}{T}_{tr}\cos \left({\theta }_2-{\theta }_2^{*}\right)-x_4^{*}\sin {\theta }_2^{*}+z_4^{*}\cos {\theta }_2^{*}-{\delta }_4^{*}\\ {\theta }_1^T={\theta }_2-{\theta }_1-{\theta }_1^{*}\\ {\theta }_2^T={\theta }_3-{\theta }_2-{\theta }_2^{*}\end{array}$$

where the thickness of the thigh and trunk are denoted as T_th and T_tr, respectively. In addition, ${\delta }_i^{*}$ and ${\theta }_i^{*}$ indicate the initial length and angle of the translational and torsional springs, respectively, and ($x_i^{*}$, $z_i^{*}$) indicate the position where the springs are fixed at backrest. The inertial, geometric values, stiffness, and damping coefficients used in this study are listed in Appendix A. Now applying Lagrange’s equation using Equation (1) yields

$$\left[\begin{array}{ccccc}{S}_{11}& 0& S_{13}& S_{14}& S_{15}\\ & S_{22}& S_{23}& S_{24}& S_{25}\\ & & S_{33}& 0& 0\\ & Sym.& & S_{44}& S_{45}\\ & & & & S_{55}\end{array}\right]\left(\begin{array}{c}{\ddot{x}}_h\\ {\ddot{z}}_h\\ {\ddot{\theta }}_1\\ {\ddot{\theta }}_2\\ {\ddot{\theta }}_3\end{array}\right)=\left(\begin{array}{c}{P}_1\\ P_2\\ P_3\\ P_4\\ P_5\end{array}\right)$$

(2)

The elements of each matrix are given in Appendix B.

2.2. Random Road Profile and Human–Vehicle Model

In order to analyze the effect of road surface on human–vehicle vibrations, we generated random road profiles according to ISO 8608, which classifies road roughness using power spectral density (PSD) [40]. An approximate form of the road displacement PSD can be written as follows:

$$\Phi \left(\Omega \right)=\Phi \left({\Omega }_0\right){\left(\frac{\Omega }{{\Omega }_0}\right)}^{-w}$$

(3)

where $\Omega$ is the angular spatial frequency and w is the waviness. The reference spatial frequency of ${\Omega }_0$ is given by 1 rad/m. In addition, the PSD function can be modified as follows [41,42].

$$\Phi \left(\Omega \right)=\frac{2\alpha {\sigma }^2}{{\sigma }^2+{\alpha }^2}$$

(4)

where ${\sigma }^2$ is the variance in the road roughness, and $\alpha$ is the coefficient with respect to the type of road surface. The road roughness variances and coefficients are listed in

Table 1. Standard deviation of the road roughness with respect to the road class at Ω₀ = 1 rad/m [42].

Road Class	σ (10−3 m)	Φ (Ω0) (10−6 m3)	α (rad/m)
A (very good)	2	1	0.127
B (good)	4	4	0.127
C (average)	8	16	0.127
D (poor)	16	64	0.127
E (very poor)	32	256	0.127

. Additionally, the road profile can be approximated using the superposition of sinusoidal functions with respect to longitudinal distance, as follows:

$$z\left(x\right)=\sum _{i=1}^N{A}_i\sin \left({\Omega }_{i}x-{\varphi }_i\right)$$

(5)

where the random phase angle is uniformly distributed in the range of $0\le {\varphi }_i<2\mathsf{\pi }$, and $A_i$ denotes the amplitudes of the road profiles, defined as follows:

$$A_i=\sqrt{\mathsf{\Phi }\left({\Omega }_N\right)\frac{\Delta \Omega }{\pi }}$$

(6)

where

$$\Delta \Omega =\frac{{\Omega }_N-{\Omega }_1}{N-1}$$

In this study, ${\Omega }_N=10\pi$ (rad/m) and ${\Omega }_1=0.02\pi$ (rad/m) are considered to generate a random road profile that covers the frequency range of the weighting functions in BS 6841. We generated the power spectral densities of five road profiles and compared them with eight road roughness references, as shown in Figure 3. Note that the higher the power spectral density, the rougher the road, and the resulting road profile falls into the ISO 8608 classification category. In addition, Figure 4 shows the road profiles generated from the power spectral densities for five different road roughness.

As shown in Figure 5, to simulate the human–vehicle model, the mathematical human model is fully coupled with the vehicle model using CarSim 9.0.2 [43] in MATLAB/Simulink. The vibrations caused by the vehicle movement are transmitted to the human body model as a base excitation caused by the displacement of the vehicle seat floor. Additionally, the reaction forces caused by the human body are transmitted to the seat floor of the vehicle through the springs and dampers. The spring and damping forces are calculated considering the displacements between the thigh and the seat pad of the human model, as follows:

$$F_R=\sum _{i=1}^4{k}_i{\delta }_i+\sum _{i=1}^4{c}_i{\dot{\delta }}_i$$

(7)

where ${\delta }_i$ and ${\dot{\delta }}_i$ indicate the translational displacements and their time rate of changes of each spring, respectively.

Figure 6 illustrates the effects of road roughness and coupling on the reaction forces and the vertical displacement of the seat floor. Figure 6a shows the vertical forces generated between the human model and the vehicle model in the case of full coupling analysis. The average reaction force corresponding to the weight of the human model is about 500 N, regardless of the road class, either A or E. However, the reaction force in the E class, which is classified as a rough road, shows greater fluctuation than when the A class is used, which is classified as a smooth road. Its magnitude is more than 20% greater than average. In addition, Figure 6b shows the vertical displacement of the driver’s seat when the human–vehicle model is connected with full coupling and partial coupling on an ISO E-class road. Full coupling means that the human–vehicle model transmits reaction forces in the mutual direction, and partial coupling means that only the vehicle transmits reaction forces to the human model. The vehicle displacement in the case of full coupling is smaller than the partial coupling case because it ignores the weight of the human model. Although the displacement difference between the full and partial couplings is not significant, coupling analysis is considered in this study because the accuracy of the simulation is crucial to the analysis of the standard deviation when the uncertainty of the human parameters is included.

2.3. Approximate Model Based on Generalized Polynomial Chaos

Suppose the uncertain parameter of the stochastic model is expressed in a set of random variables, $\mathit{\xi }=\left({\xi }_1,{\xi }_2,\dots ,{\xi }_n\right)$, the model response can also be expressed using the same set of variables [44]. Therefore, the output $\tilde{y}\left(t,\mathit{\xi }\right)$ of the stochastic model can be approximated using truncated generalized polynomial chaos expansion (gPC) expressed as the set of $\mathit{\xi }$, as given by [29,45]

$$\tilde{y}\left(t,\mathit{\xi }\right)=\sum _{j=0}^S{a}_j\left(t\right){\varphi }_j\left(\mathit{\xi }\right)$$

(8)

where ${\varphi }_j\left(\mathit{\xi }\right)$ are the generalized Askey–Wiener polynomial chaos basis functions in terms of multi-dimensional random variables $\mathit{\xi }=\left({\xi }_1,{\xi }_2,\dots ,{\xi }_n\right).$ Additionally, $a_j$ are the deterministic coefficients to be estimated. Further, the total number of polynomial terms, $S$, can be obtained as follows:

$$S=\frac{\left(p+n\right)!}{p!n!}-1$$

(9)

here, p indicates the order of the polynomial chaos, and n indicates the number of uncertain parameters. The orthogonal polynomial functions according to the Askey–Wiener polynomial basis are listed in

Table 2. Distributions of random variables and corresponding polynomial basis function.

Distribution	Polynomial Function	Support
Gaussian	Hermite	(−∞, +∞)
Gamma	Laguerre	[0, +∞)
Beta	Jacobi	[a, b]
Uniform	Legendre	[a, b]

.

In this study, we assumed that the uncertain parameter has a Gaussian distribution. For the random variable with a Gaussian distribution, generalized polynomial chaos expansion employs the Hermite polynomials. The general expression of Hermite polynomials of order p is expressed as follows:

$${\varphi }_p\left({\xi }_1,{\xi }_2,\dots ,{\xi }_p\right)={{\left(-1\right)}^{p}e}^{\frac{1}{2}{\mathit{\xi }}^T\mathit{\xi }}·\frac{{\partial }^p}{\partial {\xi }_1·\partial {\xi }_2\cdots \partial {\xi }_p}·e^{-\frac{1}{2}{\mathit{\xi }}^T\mathit{\xi }}$$

(10)

where superscript ^T indicates the transpose of the random variable vector.

The unknown coefficients of polynomial chaos expansion are commonly obtained through the use of the Galerkin method. However, it is not easy to apply the Galerkin method when the stochastic model equation is complicated and difficult to tractable. Therefore, in this study, the unknown coefficient was estimated using the collocation method. The method is accomplished by computing the output of the stochastic model at a set of collocation points. The collocation points of the random variable are the roots of the next higher-order polynomial [36,44]. For example, in a one-dimensional random field, the collocation points of the second-order polynomial chaos are calculated from the root of the third-order Hermite polynomial $\left(-\sqrt{3},0,\sqrt{3}\right)$. In the collocation method, the output of the stochastic model and approximation model based on gPC are the same at the collocation point, as follows:

$$\mathit{y}\left(t,\mathit{\xi }\right)=\tilde{\mathit{y}}\left(t,\mathit{\xi }\right)=\mathbf{T}\left(\mathit{\xi }\right)\mathit{a}\left(t\right)$$

(11)

where $\mathit{y}\left(t,\mathit{\xi }\right)$ represent the output vectors of the stochastic model at the collocation points, ${\mathit{\xi }=(\xi }_1,{\xi }_2,\cdots ,{\xi }_M).$ The subscript of M denotes the number of collocation points. The transform matrix $\mathbf{T}$ and unknown coefficients $\mathit{a}\left(t\right)$ are as follows:

$$\mathbf{T}\left(\mathit{\xi }\right)=\left(\begin{array}{cccc}{\varphi }_0\left({\xi }_1\right)& {\varphi }_1\left({\xi }_1\right)& \cdots & {\varphi }_S\left({\xi }_1\right)\\ {\varphi }_0\left({\xi }_2\right)& {\varphi }_1\left({\xi }_2\right)& \cdots & {\varphi }_S\left({\xi }_2\right)\\ ⋮& ⋮& \ddots & ⋮\\ {\varphi }_0\left({\xi }_M\right)& {\varphi }_1\left({\xi }_M\right)& \cdots & {\varphi }_S\left({\xi }_M\right)\end{array}\right),\mathit{a}\left(t\right)=\left(\begin{array}{c}{a}_1\left(t\right)\\ a_2\left(t\right)\\ ⋮\\ a_M\left(t\right)\end{array}\right)$$

Applying the least square method to Equation (10), the vectors of unknown coefficient can be easily obtained, as follows:

$$\mathit{a}\left(t\right)={\left\{\mathbf{T}{\left(\mathit{\xi }\right)}^T\mathbf{T}\left(\mathit{\xi }\right)\right\}}^{-1}\mathbf{T}{\left(\mathit{\xi }\right)}^T\mathit{y}\left(t,\mathit{\xi }\right)$$

(12)

Using the calculated coefficients of the orthogonal polynomials, the mean and variance can be obtained, as follows:

$$\mu =a_0{\varphi }_0\left(\xi \right)\mathrm{a}\mathrm{n}\mathrm{d}{\sigma }^2=\sum _{j=1}^S{a}_j^2⟨{\varphi }_j^2⟩$$

(13)

where $⟨·,·⟩$ denotes the ensemble average inner product [29].

We used an approximate model based on gPC to analyze the stochastic results of the human–vehicle model instead of Monte Carlo simulation (MCS), which requires a large number of simulations. The proposed approximate model consists of four major steps. First, the selection of the uncertain parameters according to an assumed distribution. Second, the generation of an approximate model through generalized polynomial chaos expansion. Third, the estimation of the coefficients using the collocation method, and finally, a comparison of the stochastic results obtained through gPC with those obtained through MCS. This process of the calculation is illustrated in Figure 7. In essence, gPC can compute the stochastic results of a dynamic system using the polynomials, as opposed to fully integrating the equations of motion through the use of MCS. Therefore, the computational time of gPC is significantly reduced compared to that of MCS. Essentially, in contrast to fully integrating the equations of motion through MCS, gPC can use polynomials to compute the stochastic results of a dynamic system. Therefore, the computational time of gPC is significantly decreased compared to MCS.

3. Stochastic Analysis of Ride Comfort Index

In this section, we analyze the effect of uncertain parameters on ride comfort. We use the generalized polynomial chaos (gPC) method to obtain the stochastic responses of the human model. The gPC results are compared with Monte Carlo simulations (MCS) in terms of accuracy and time efficiency. Finally, a parametric study of stochastic ride comfort is performed.

3.1. gPC Results for Human–Vehicle Model with Uncertain Parameters

A fully coupled human–vehicle simulation was performed by combining CarSim and MATLAB/Simulink in the presence of the parametric uncertainty of the human body model. In this stochastic analysis, it is assumed that there are three uncertain parameters, such as two translational stiffness ($k_1$, $k_3$) and one torsional stiffness ($k_1^T$). It also assumes that the translational springs connected to the same human body are identical to reduce computational costs, that is, $k_1=k_2$ and $k_3=k_4$ (see Figure 1). A single uncertain parameter was selected for the dynamic simulations in turn, and the accelerations of the human model were extracted with respect to the corresponding parameter. The second order of the polynomial chaos was considered to obtain the stochastic responses of the human model. In other words, single-dimensional second-order polynomial chaos expansion was used in the approximate model based on gPC. In addition, those with uncertain stiffness were assumed to be the Gaussian distributed variable. They can be represented as follows:

$$\begin{array}{c}{k}_1\left({\xi }_1\right)={\mu }_{k_1}+{\sigma }_{k_1}{\xi }_1\\ k_3\left({\xi }_2\right)={\mu }_{k_3}+{\sigma }_{k_3}{\xi }_2\\ k_1^T\left({\xi }_3\right)={\mu }_{k_1^T}+{\sigma }_{k_1^T}{\xi }_3\end{array}$$

(14)

where ${\mu }_k$ and ${\sigma }_k$ represent the mean and standard deviation values of the corresponding stiffness, respectively. In this study, the standard deviation of the uncertain stiffness was assumed to be 5% of the corresponding mean value. The conditions of the human–vehicle simulations are as follows. The simulated vehicle model used in CarSim was an E-class sedan. In each segment, the accelerations of the human model were calculated when the vehicle was traveling at a constant speed of 60 km/h in the longitudinal direction. The steering angle of the vehicle was controlled to maintain straight-line driving. The total running time and sampling frequency were set to 50 s and 1000 Hz, respectively. However, in the dynamic simulations, only times between 5 and 50 s are analyzed to eliminate transient effects due to accelerations fluctuating to the target speed.

Since MCS requires a huge number of simulations to achieve satisfactory convergence, we also developed an automatic process for human–vehicle simulations using the COM interface in CarSim. Automation is implemented by generating a set of samples for uncertain parameters and controlling the COM interface with a MATLAB script. A detailed description of the automatic process follows. First, generate the samples corresponding to the probability distribution in the MATLAB script. The number of samples at this stage depends on the number of MCS and gPC runs. Second, launch MATLAB/Simulink in CarSim after executing the COM interface of CarSim in the MATLAB script. Third, provide the generated set of samples as input to the function corresponding to the human model in MATLAB/Simulink. Fourth, once the simulation is performed, save the response of generalized coordinates of the human model for each sample in the workspace of MATLAB. Then, terminate the CarSim and MATLAB/Simulink via MATLAB script. Finally, repeat steps 2 to 4 as many times as the number of simulations of MCS and gPC.

For the three uncertain parameters $k_1$, $k_3$, and $k_1^T$, the variations of the vertical acceleration at the hip joint of the human model were calculated using gPC and MCS, as shown in Figure 8. Vertical accelerations close to zero are displayed in light gray to indicate the distribution of probability densities (see the color bar on the right). Several results are observed. First, it is clear that the accelerations calculated by the approximation model based on gPC show good agreement with those obtained through the use of MCS for all three types of uncertain parameters. To obtain a satisfactory convergent result, gPC requires only three times simulations for three collocation points, while MCS does 500 times simulations. The former is at least 160 times more efficient than the latter. Second, although the average values of the acceleration are the same, the distributions of the probability densities are different depending on the uncertain parameters. Third, when $k_1^T$ is uncertain in Figure 8e,f, the probability density of acceleration appears as a solid line instead of a distribution due to a small standard deviation. This is because that $k_1^T$ is located at the very position of the hip joint that we calculated.

In order to investigate the effect of the variations on vertical accelerations, we calculated the probability densities for the uncertain parameter $k_1$ and $k_3$ in detail. For example, the probability density function of the acceleration at t = 21.1 s (red dashed line in Figure 8) is shown in Figure 9. Here, the shaded blocks and solid lines represent the MCS and gPC results, respectively. In the case of uncertain parameter $k_1$, a wider distribution is observed than that of $k_3$. From a design point of view, it is crucial that the variation of the acceleration shows a different distribution depending on the uncertainty of the parameter. However, it is not practical to examine the variations for all time responses. Therefore, it is necessary to quantify a specific value that takes into account all of the time responses.

Further, we used BS 6841 to quantify subjective ride comfort as a ride index [37]. To calculate the ride index, the acceleration for each axis of the human model is taken from the human–vehicle simulation, as illustrated in Figure 10. Applying Fourier transform, frequency weighting can be imposed using a convolution integral. The inverse Fourier transform of the frequency-weighted acceleration produces the time history of the acceleration, which can be weighted again by multiplying the factor at each axis. Now the root mean square of the time-weighted acceleration on each axis represents the ride index of each axis. Since the developed human model considers two-dimensional motion, in this study, only five accelerations of vertical, horizontal, and rotational directions at the hip joint and vertical and horizontal directions at the back are used in computations. Summing the riding indices of each axis, we can finally quantify the overall ride index of the human model.

Figure 11 shows the overall ride indices with respect to road roughness for three different uncertain parameters with the discomfort scale represented in BS 6841. We computed all ride indices and standard deviations by the approximate model based on gPC for each uncertain parameter. At each road roughness, the solid circle represents the mean of the ride indices, and the error bar represents 3σ, which implies that 99.7% of the data are in this range. The overall ride indices associated with class-A and B road roughness are placed in the ‘not uncomfortable’ zone, whereas those associated with class-E represent the ‘fairly uncomfortable’ zone. In general, class-A indicates a smooth surface, while class-C indicates average road roughness (see Table 1). Therefore, it can be predicted easily that the ride indices at E-class would place in a ‘fairly uncomfortable’ zone. It is also easy to understand that the standard deviation on rough road is greater than on smooth road. However, it is not easy to predict the parameter sensitive to ride index. Now comparing the error bar for uncertain stiffness, it is clear that the ride index is dominantly influenced by $k_1$ and $k_3$, which are the translation stiffness related with thigh and back of the human and seat. However, $k_1^T$, which is the torsional stiffness at hip joint of the human, is not sensitive to the ride index.

3.2. Sensitivity Analysis of Ride Comfort

To improve ride comfort, it is important to understand which segment of the vehicle or human body is closely related or sensitive to ride comfort. In the presence of uncertainties of a parameter, we performed human–vehicle simulation and obtained the probability distributions of ride index as illustrated in Figure 12. The results for translational parameters ($k$ and $c$) are displayed in the left column, and those for torsional parameters ($k^T$ and $c^T$) are in the right column. The simulations were performed when the standard deviation of each uncertain parameter has 5% of its mean value. As it may be expected, the shape of distributions is quite different depending on the uncertain parameters. Surprisingly, the range of the ride index for the translational parameters is at least 20 times wider than those for the torsional parameters. This implies that the translational parameters have much stronger effects on the uncertainty of the ride index than the torsional parameters. It can also be interpreted as the ride index is dominated by the fundamental frequency, which is mainly influenced by the translation parameter of human model [38]. In addition, it can be seen that the distribution for $k_2^T$ and $c_2^T$ have a much shaper shape than $k_1^T$ and $c_1^T$, as shown in Figure 12b,d and

Table 3. Stochastic properties of ride index for the case of Figure 12.

Stiffness	${\mathit{k}}_1$	${\mathit{k}}_3$	${\mathit{k}}_1^{\mathit{T}}$	${\mathit{k}}_2^{\mathit{T}}$
Mean	0.4306	0.4304	0.4312	0.4312
Std	0.85 × 10−3	1.36 × 10−2	1.45 × 10−5	2.27 × 10−6
Damping	${\mathit{c}}_{\mathbf{1}}$	${\mathit{c}}_{\mathbf{3}}$	${\mathit{c}}_{\mathbf{1}}^{\mathit{T}}$	${\mathit{c}}_{\mathbf{2}}^{\mathit{T}}$
Mean	0.4311	0.4311	0.4312	0.4312
Std	0.62 × 10−2	0.52 × 10−2	0.97 × 10−5	2.39 × 10−6

. In other words, the hip joint is more sensitive than the neck part of human body. In a physical sense, it is easy to verify because the acceleration of the hip joint and trunk is more dominant than others in evaluating ride comfort [46]. In conclusion, it was found that the variations in $k_2^T$ and $c_2^T$, which connects the trunk and the head, has little effect on the movement of the hip joint and the trunk.

It is also interesting to investigate the differences between the linear and nonlinear human models concerning human–vehicle simulation in terms of stiffness uncertainties. In this study, linearization was performed using Taylor expansion under a small motion assumption. The mean and standard deviation of the ride index are presented in Figure 13 at the top and bottom for the nonlinear and linear models, respectively. The ride indices were plotted as colormaps to simultaneously look over the effect of the uncertain parameters. For comparison, each parameter varied from 20% to 180% of the reference value and was expressed as a dimensionless value denoted by a superscript *. For example, $k_1^{\mathrm{*}}=1$ means an initial value, and $k_1^{\mathrm{*}}=0.2$ means the 20% of the initial value. The standard deviation of each uncertain parameter was also set to 5% of the mean value, and the human–vehicle simulation was performed on a C-class road. In conclusion, the mean values of the ride index for the linear and nonlinear models show similar trends. The higher the stiffness, the worse the ride comfort. However, the mean of the ride index of the nonlinear model shows a higher value than that of the linear model, and it is probably because the nonlinear model allows for large movements. For example, when $k_1^{\mathrm{*}}$ and $k_3^{\mathrm{*}}$ increases by 40%, the ride index increases by 25.4% from 0.451 to 0.604 in the nonlinear model (Figure 13A), while it increases by 20.5% from 0.421 to 0.529 in the linear model (Figure 13B). On the contrary, for damping, the higher the damping, the better the ride comfort. When $c_1^{\mathrm{*}}$ and $c_3^{\mathrm{*}}$ decreases by 40%, the ride index increases by 21.4% from 0.451 to 0.574 in the nonlinear model (Figure 14A), while it increases by 18.6% from 0.421 to 0.517 in the linear model (Figure 14B).

When the stiffness parameter is uncertain, the standard deviation of the ride index tends to be different depending on the presence of nonlinearity in the human body model. Interestingly, when there is uncertainty only in $k_1$, the maximums of the standard deviation exist at different locations (compare C and D in Figure 13). As mentioned earlier, the higher the standard deviations, the greater the sensitivity, and the results imply that the sensitive ranges of the uncertain parameters are significantly different. Similar results also occur in terms of $k_3$ (compare E and F in Figure 13). In the case of $c_1$ and $c_3$, the effect of standard deviations in the nonlinear and linear models are similar, but their maximums are slightly different (compare C, D and E, F in Figure 14). Based on the above results, it can be concluded that the nonlinear model responds more sensitively to the variations of the seating parameters than the linear model, leading to the need for a nonlinear model for more accurate human–vehicle simulations.

4. Conclusions

In the presence of parameter uncertainty, this study analyzed the ride comfort stochastically through the simulation of a fully coupled human–vehicle model using CarSim software combined with MATLAB/Simulink. Ride comfort was evaluated quantitatively as a ride index using a human–vehicle model according to the BS-6841 standard. In addition, a generalized polynomial chaos extension (gPC) framework was proposed to compute the stochastic response of the human model. Comparing the computational time of gPC with Monte Carlo simulation (MCS), the former proved to be 160 times faster than the latter.

We also carried out a parametric study of stochastic ride comfort in terms of nonlinearities and uncertain parameters in human models. Two main results were observed. First, the translational parameters of the human model are much more important than the torsional parameters when comparing the standard deviations of the ride index. Second, the use of a nonlinear human model shows that the standard deviation of the ride index is significantly different from the linear model. Based on the above observations, it can be concluded that the uncertain quantification method and the nonlinear human model are needed for more accurate and precise ride comfort analysis.

This study also showed the distribution of ride comfort according to the road surface when there is a deviation in the stiffness and damping coefficient for drivers of average body type and weight. In a real driving environment, factors such as posture, road surface, and driving conditions are constantly changing, leading to changes in weight distribution and skin/seat deformation. As a result, even for the same driver or passenger, the response of the human body to vibration becomes unpredictable and uncertain. We expect that the proposed gPC framework can play an important role in automotive engineering as the demand for ride comfort continues to increase. For example, automotive engineers can gain insights from this study to design seats that are robust in terms of ride comfort. In this study, the parameters of k₁, k₃, c₁, and c₃ of the human model represent the springs and dampers that connect the human body to the seat cushion. Therefore, based on the results of the study, it is possible to find the optimal seat design parameters that improve ride comfort with a small standard deviation.