# An Efficient Polynomial Chaos Method for Stiffness Analysis of Air Spring Considering Uncertainties

Authors: Feng Kong; Penghao Si; Shengwen Yin
Journal: Complexity
 (2021)

---

## Abstract

Traditional methods for stiffness analysis of the air spring are based on deterministic assumption that the parameters are fixed. However, uncertainties have widely existed, and the mechanic property of the air spring is very sensitive to these uncertainties. To model the uncertainties in the air spring, the interval/random variables models are introduced. For response analysis of the interval/random variables models of the air spring system, a new unified orthogonal polynomial expansion method, named as sparse quadrature-based interval and random moment arbitrary polynomial chaos method (SQ-IRMAPC), is proposed. In SQ-IRMAPC, the response of the acoustic system related to both interval and random variables is approximated by the moment-based arbitrary orthogonal polynomial expansion. To efficiently calculate the coefficient of the interval and random orthogonal polynomial expansion, the sparse quadrature is introduced. The proposed SQ-IRMAPC was employed to analyze the mechanic performance of an air spring with interval and/or random variables, and its effectiveness has been demonstrated by fully comparing it with the most recently proposed orthogonal polynomial-based interval and random analysis method.


---

## Body

## 1. Introduction

Air springs are widely used for isolating vibrations to enhance comfort for the passengers in railway vehicles. Stiffness is an important characteristic parameter to evaluate the performance of an air spring. Traditional methods to analyze the stiffness of air spring are deterministic methods, in which the input parameters are fixed [1]. However, due to the manufactory error and other factors, the uncertainty related to the material properties is inevitable. It is shown in [1] that the stiffness of the air spring is very sensitive to the input parameters such as cord angle and material of cord. Thus, it is desirable to predict the response of the air spring system with uncertainties.To deal with the uncertainty in the air spring, a suitable uncertain model should be established to describe the uncertain parameters. The probabilistic method can be viewed as the most valuable uncertainty analysis method as it can provide detailed statistical properties of response for engineering design [2]. In the probabilistic method, the uncertainty is modeled as a random variable whose detailed statistical properties are available. Research studies on response analysis of the random uncertain system are rather mature, and different probabilistic methods have been proposed, including Monte Carlo method (MCM) [3], perturbation stochastic method [4], and polynomial chaos method [5, 6]. Among various probabilistic methods, the polynomial chaos method is the most popular uncertainty propagation approach for random analysis due to its good accuracy and efficiency.In the polynomial chaos method, the response of interest is approximated in terms of orthogonal polynomials [7]. By using the optimal random polynomial basis, the polynomial chaos method can converge exponentially for random analysis. In order to construct the optimal random polynomial basis, lots of polynomial chaos methods have been developed [8–10]. The generalized Polynomial chaos (gPC) has become the most widely used polynomial chaos method due to the ease of the construction of polynomial basis [7]. In the gPC, the orthogonal polynomial from the Askey scheme is employed to establish the polynomial chaos expansion. However, the Askey scheme can only provide the optimal polynomial basis for some well known probability distributions. For the random problem with complex probability distribution, a nonlinear variable transformation should be employed in gPC and the accuracy of gPC may be deteriorated. To obtain the optimal random polynomial basis related to arbitrary probability distribution, the Arbitrary Polynomial Chaos (APC) was then developed. In the framework of APC [8], the optimal random polynomial basis related to arbitrary PDFs can be numerically constructed. Originally, the polynomial basis of APC was constructed according to the PDF by using Gram–Schmidt orthogonalization. Subsequently, Oladyshkin et al. found that the APC can be directly based on a matrix of the moment of the random variable [9]. As the Moment-based Arbitrary Polynomial Chaos (MAPC) can avoid the fitting error introduced by the construction of the PDF for random variable, the convergence of the MAPC was significantly better than that of the conventional APC [10].In these above random methods, a huge number of statistical data is needed to obtain the moment or PDF of random variables. However, in practice, the available data may be not sufficient to construct the precise statistical property of the random variable. For the uncertain problem with limited probability information, lots of nonprobabilistic methods and imprecise probabilistic methods have been developed, such as the interval analysis [11, 12], the fuzzy sets [13, 14], the p-box theory [15], and the evidence theory [16–18]. In most of these nonprobabilistic methods, a definition of the imprecise probability of uncertain parameters is required, while in the interval approach, only the bounds of variables should be obtained. In practice, as the bounds for an uncertain parameter is more convenient to be obtained, the interval analysis technique has become a very popular nonprobabilistic method. During last decades, a lot of interval analysis techniques have been developed, such as the first-order perturbation technique [19], the modified first-order perturbation technique [20], the interval factor technique [21], the rational expansion technique [22], the polynomial interval expansion methods [23, 24], and the vertex method [25, 26].As aforementioned, the probabilistic method and interval method are suitable to solve the uncertain problem with pure random variables or pure interval variables. But, in real engineering, the interval and random variables may exist simultaneously if the available information on uncertain parameters is different. Therefore, the engineering application prompts a growing demand for developing the hybrid interval and random analysis method. In particular, lots of hybrid interval and random analysis methods have been proposed based on the perturbation technique and the polynomial chaos theory. The hybrid perturbation method can achieve high efficiency for uncertainty analysis with both the interval and random variable, but it is limited to the uncertain problem with the small uncertainty level [27]. To solve the hybrid uncertain problem with the large uncertainty level, Wu et al. proposed a hybrid polynomial chaos method by integrating the Chebyshev interval expansion and gPC [28]. Based on gPC and dimension-wise analysis technique, Xu et al. developed a hybrid uncertainty analysis method for response analysis of the structure-acoustic system [29]. Based on the parametric Gegenbauer polynomial from the Askey scheme, Yin et al. proposed a unified polynomial chaos method for interval and random analysis of the acoustic system [30]. To improve the accuracy of the polynomial chaos method for the hybrid uncertain structure-acoustic problem with complex probability distributions, the APC has been recently developed for hybrid interval and random analysis [31, 32]. As the APC can obtain the optimal random polynomial basis related to arbitrary PDFs, the APC-based hybrid uncertainty analysis method shows a better accuracy than the traditional gPC-based hybrid uncertainty analysis method.From an overall perspective, the APC has achieved great success in the field of interval and random analysis. In the interval and random arbitrary polynomial chaos method, the Gauss quadrature is used to calculate the expansion coefficient. However, the computational burden of IRAPC by using Gauss quadrature will increase exponentially with the number of variables. The air spring system always involves lots of uncertain parameters, and the computational burden for finite element analysis of the air spring is relatively large. Thus, it is necessary to improve the computational efficiency of IRAPC for uncertainty quantification especially when it is used for response analysis of the air spring.The paper aims to develop a robust orthogonal polynomial expansion approach for interval and/or random analysis of the air spring system with arbitrary statistical data. To avoid the errors introduced by estimating the PDF of random variables, the moment-based polynomial chaos is introduced for hybrid interval and random analysis. Especially, the choice of optimal polynomial basis for interval analysis is justified in the framework of moment-based polynomial chaos theory. To efficiently calculate the expansion coefficient of the moment-based polynomial chaos, the sparse quadrature is introduced. Based on the moment-based interval and random polynomial chaos expansion and the sparse grids quadrature, a new method named as Sparse Quadrature-based Interval and Random Moment Arbitrary Polynomial Chaos (SQ-IRMAPC) method is proposed. The proposed method has been used to analyze the uncertainty of response for the air spring system, and the effectiveness of the proposed method has been also compared with the recently proposed interval and random APC method. Note that the main difference between the proposed method, the recently proposed interval, and random APC method is that different integration methods are used to calculate the expansion coefficient.


## 2. Basic Theory of MAPC

In MAPC, the orthogonal polynomial, which is derived from the moment of the random variable, will be used to approximate the response of interest. The MAPC expansion for the approximation of a function can be expressed as follows:(1)Yx=∑i=0Nyiφix,where N is the retained order of MAPC expansion, yi represents the expansion coefficient to be estimated, φix denotes the polynomial basis of order i, and φix varies for different random variables. In MAPC [9], the polynomial basis can be obtained according the moment of the random variable without knowing the PDF. As a comparison, the polynomial basis in the gPC [7] and the original APC [8] is obtained on the condition that the PDF of the random variable is well defined.The orthogonal polynomial of MAPC can be denoted by the three-term recurrence as follows:(2)φ−1ξ=0,φ0ξ=1,φk+1ξ=ξ−akφkξ−bkφk−1ξ,k=0,1,2,….In the above equation,aj and bj denote the recurrence coefficients to be estimated. In the framework of MAPC, aj and bj can be derived from a Hankel matrix of the random moments. The Hankel matrix of the random moments can be expressed as [9](3)M=μ0μ1…μpμ1μ2μp+1⋮⋱μpμp+1μ2p,where μkk=0,1,…,2p is the kth moment of the random variable and p should be no less than the retained order of MAPC expansion.Under the condition that the Hankel matrix is positive definite, the Cholesky decomposition ofM can be computed, namely, M=RTR. Rutishauser [32] has derived explicit analytic formulas to obtain the recurrence coefficient from the Cholesky matrix entries R. The recurrence coefficients in terms of rij can be expressed as(4)aj=rj,j+1rj,j−rj−1,jrj−1,j−1,bj=rj+1,j+1rj,j,where riji=1,2,…;j=1,2,… is the ith row jth column element of R, r0,0=1, and r0,1=0. The polynomial basis obtained through equations (2)∼(4) can satisfy the following orthogonality relation:(5)∫Ωφixφjxρxdx=δij,where δij denotes the Kronecker delta and ρx denotes the PDF of a random variable. In the framework of MAPC, ρx can be an arbitrary discrete or continuous function.Based on orthogonality of the polynomial basis,yi in equation (1) can be calculated by [32](6)yi=Yx,φixφix,φix=1hi∫ΩYxφixρxdx.The integral in the above equation can be calculated by the Gaussian quadrature as follows [32]:(7)yi=1hi∫ΩYxφixρxdx=1hi∑i=1mYx^iφix^iw^i,where x^i and w^i are the Gaussian nodes and the Gaussian weights, respectively, m is the total number of Gaussian nodes, and x^i and w^i can be obtained from the eigenvalue decomposition of the Jacobi matrix assembled with ai and bi. Particularly, the Jacobi matrix Jn can be expressed as [9](8)Jn=a1b1b1a2b2b2⋱⋱⋱an−1bn−1bn−1an.In particular, ifVTJnV=diagλ1,λ2,…,λn and VTV=I, in which I is the n×n dimension identity matrix, then the desired x^i and w^i can be determined by(9)x^i=λi,w^i=b0vi,12,i=1,2,…,where vi,1 is the first component of the ith column vector of V. The coefficient b0 can be arbitrary and set by convention such that b0=∫wxdx.


## 3. SQ-IRMAPC Method

In this section, the MAPC is introduced for interval and random analysis, and the IRMAPC expansion is established. In addition, the sparse quadrature is used to calculate the expansion coefficient of IRMAPC. By a combination of the IRMAPC and the sparse quadrature, a new hybrid uncertainty analysis method, named as the Sparse Quadrature-based Interval and Random Moment Arbitrary Polynomial Chaos (SQ-IRMAPC) method, is proposed.


### 3.1. IRMAPC Expansion

In this section, both the random and interval variables are considered inY=FxR,xI. The function contains random variable xR and an interval variable xI. Hence, the output of the function will have the characteristics of both random and interval variables, and the proposed IRMAPC method will approximate the response by IRMAPC expansion.The IRMAPC expansion for a function can be expressed as(10)F=FxI,xR=∑0≤i1+⋯+iL1+⋯+iL≤nfi1,…,iLφi1,…,iL1xIφiL1+1,…,iLxR,where(11)φi1,…,iL1xI=φi1x1I×⋯×φiL1xL1I,φiL1+1,…,iLxR=φiL1+1xiL1+1R×⋯×φiLxLR.In traditional IRMAPC, the expansion coefficients are calculated by(12)fi1,…,iL=1hi1×⋯×hiL∑j1=1M1…∑jL=1MLFx^I,x^Rφi1,…,iL1x^IφiL1+1,…,iLx^Rw^i1,…,iL,where x^I=x^j1,…,x^jL1 and x^R=x^jL1,…,x^jL denote the Gaussian nodes of interval variables and random variables, respectively, w^j1,…,jL=∏k=1Lw^jk is the Gaussian weight, Mjj=1,2,…,L denotes the number of Gaussian nodes related to the jth variable, x^jk and w^jk denote the jkth Gaussian node and weight related to the kth variable, respectively, and x^jk and w^jk can be calculated through equations (8) and (9).It can be found from equation (12) that the total number of Gaussian points to determine the coefficient is N=∏k=1LMk. Obviously, the total number of Gaussian points will increase exponentially with the increasing number of uncertain parameters, which may lead to tremendous computational cost. In order to improve the computational efficiency of moment-based polynomial chaos expansion for interval and random analysis, the sparse Gaussian quadrature will be introduced to calculate the expansion coefficient.


### 3.2. Determine the Expansion Coefficient by Using Sparse Quadrature

The sparse quadrature is based on the Smolyak algorithm, which has been widely used in the fields of numerical integration and interpolation and image processing. In this section, the basic principles of sparse quadrature for calculating the expansion coefficient will be deduced.A continuous functionFx defined on x∈−1,1 approximated by polynomial chaos expansion can be denoted as Ql1F, where l denotes the retained order of polynomial chaos expansion. According to the nested hierarchical basis principle of the Smolyak algorithm, the difference format of the approximated function is [33](13)Δk1F=Qk1−Qk−11F,Q01F=0.Furthermore, for ad-dimension problem, the approximated function with the order l of the Smolyak algorithm can be contrasted as(14)QldF=∑k≤l+d−1Δk11⊗⋯⊗Δkd1F,where k denotes the sum of the multidimensional indicators k=∑i=1dki and ⊗ expresses the operation of the tensor product. Through the operation of the tensor product, equation (14) can be expressed as(15)Q1dF=∑l+1≤k≤l+d−1l+d−kd−1l+d−kQk11⊗⋯⊗Qkd1F.Therefore, the integration points in the square grids can be defined as(16)U1d=∪l+1≤k≤l+dUk11⊗⋯⊗Ukd1.The number of the integration points based on the sparse grid method is estimated by(17)N≈2ll!dl.The corresponding coefficient of the weights is(18)wk1…kdi1…id=−1l+d−kd−1l+d−kwk1i1⊗⋯⊗wkdid.


### 3.3. Optimal Choice of the Polynomial Basis for IRMAPC Expansion

According to Witteveen and Bijl [8], the accuracy of MAPC expansion by using different polynomial basis varies greatly. Therefore, when the MAPC expansion is extended for interval analysis and hybrid analysis, a key procedure is to determine the optimal polynomial basis related to each random variable and each interval variable. For the random variable, the polynomial basis derived from the moment of the random variable (through equations (2)∼(4)) can be viewed as the optimal random polynomial basis [11]. For the interval variable, there may also exist an optimal interval polynomial basis. In the following text, the optimal choice of polynomial basis of MAPC for interval analysis will be discussed.The main aim of this section is to choose an optimal interval polynomial basis for interval analysis in the framework of moment-based arbitrary orthogonal polynomial expansion. The interval analysis is to find the maximum and minimum value of the response. Theoretically, the accuracy of moment-based arbitrary polynomial interval expansion can be improved by reducing the error of moment-based arbitrary orthogonal polynomial expansion at the maximum and minimum. However, the maximum or minimum may be an arbitrary point in the range of variation of the interval variable. Therefore, when the moment-based arbitrary orthogonal polynomial expansion is used to approximate the response, we can determine the maximum value of the possible error ofPnx over an interval as follows [33]:(19)ε∞=maxa≤x≤bFx−PnxFx=Fx−PnxFx∞,where Pnx denotes the nth order moment-based arbitrary orthogonal polynomial expansion. Among all the nth order moment-based arbitrary orthogonal polynomial expansions, if Pn∗x satisfied(20)Fx−Pn∗xFx∞=minPnx∈Hnε∞=minPnx∈Hnmaxa≤x≤bFx−PnxFx,then Pn∗x can be viewed as the optimal interval polynomial basis. In equation (14), Hn denotes the universal set of all nth order moment-based arbitrary orthogonal polynomial expansion. According to the theory of Weierstrass [33], Pn∗x is available. However, when n≥1, the computation of Pn∗x is rather difficult. In addition, the polynomial basis of Pn∗x varies for different retained order n. Thus, in general, a polynomial basis whose ε∞ is relatively small is used for numerical approximation instead of Pn∗x. In this paper, a moment-based arbitrary polynomial basis whose ε∞ is relatively small will be used as the interval polynomial basis.According to the theory of orthogonal expansion, the polynomial basis depends on its related weight function or moment. In order to determine interval polynomial basis, the effect of the weight function onε∞ will be investigated through a simple numerical example.Example: suppose thaty=ex2, where the range of variation of x is [−1, 1].Different weight functions are considered, including the weight function of Chebyshev expansion; the weight function whose weight is zero at the midpoint or the bounds of the interval [−1, 1] is shown in Figure1; these two types of weight functions are named as convex shape weight function and concave shape weight function, respectively. For clarity, the expression of each weight function is listed in Table 1. ε∞ of the moment-based arbitrary orthogonal polynomial expansion with different weight functions is shown in Figure 2.Figure 1
Three different types of weight functions.Table 1
Expression of different weight functions.
Type of weight functionExpression of weight functionChebyshev expansion1/π1−x2,−1≤x≤1Constant1/2,−1≤x≤1Concave shapeΓα+β+1/2α+β+1Γα+1Γβ+11−ξα1+ξβ,α=β=0.01,−1≤x≤1Convex shapeΓα+β+1/2α+β+1Γα+1Γβ+11−ξα1+ξβ,α=β=12,−1≤x≤1Figure 2
ε∞ of moment-based arbitrary polynomial interval expansion with different weight functions.It can be found from Figure2 that ε∞ of the moment-based arbitrary orthogonal polynomial expansion is relatively large with the convex or concave shape weight function. The main reason may be that the moment-based arbitrary orthogonal polynomial expansion will introduce large error at the domain where the weight is approaching to zero [32]. As a comparison, ε∞ of the moment-based arbitrary orthogonal polynomial expansion with the Chebyshev weight function or constant weight function is small. In addition, it can be found from Figure 2 that ε∞ of the moment-based arbitrary orthogonal polynomial expansion with the Chebyshev weight function is relatively small than that with the constant weight function. Based on the above analysis, the Chebyshev weight function will be used as the weight function of interval polynomial basis, which can be expressed as(21)wIx=1π1−x2.From the application of MAPC for random analysis, we can find that when the moment is set asμk=∫−11xkρxdx, the polynomial can be orthogonal ρx. It indicates that the polynomial basis will be orthogonal to the weight function of Chebyshev polynomial if the moment is determined by(22)μk=∫−11xkwCx=∫−11xk1π1−x2dx,where wCx=1/π1−x2 is the weight function of the Chebyshev polynomial. Therefore, the moment for each interval variable can be determined according to equation (22).Once the moment is determined, the polynomial basis for interval analysis can then be calculated through equations (2)∼(4).


### 3.4. Calculation of the Bounds of Statistics

To calculate the bounds of statistics, there are two main steps. In the first step, the interval variables are regarded as constant parameters, and the MAPC expansion can be rewritten as the following form:(23)FxI,xR=∑iL1+1=0NL1+1…∑iL=0NL∑i1=0N1…∑iL=0NL1fi1,…,iLφi1,…,iL1xIφiL1+1,…,iLxR,=∑iL1+1=0NL1+1…∑iL=0NLziL1+1,…,iLkφiL1+1,…,iLxR,where(24)ziL1+1,…,iLk=∑i1=0N1…∑iL1=0NL1fi1,…,iLkφi1,…,iL1xI.Due to the orthogonality of the polynomial basis, the expectation and variance of the response can be expressed in terms of the expansion coefficient as follows:(25)μ=E∑iL1+1=0NL1+1…∑iL=0NLziL1+1,…,iLkφiL1+1,…,iLxR=z0,…,0k,σ2=E∑iL1+1=0NL1+1…∑iL=0NLziL1+1,…,iLkφiL1+1,…,iLxR2−μ2,=∑iL1+1=0NL1+1…∑iL=0NLziL1+1,…,iLk2hiL1+1…hiL−z0,…,0k2.Substituting equation (18) into (19) and (20), we can rewrite the expectation and variance of the response as follows:(26)μ=∑i1=0N1…∑iL1=0NL1fi1,…,iL1,0,…,0kφi1,…,iL1xI,(27)σ2=∑iL1+1=0NL1+1…∑iL=0NL∑i1=0N1…∑iL=0NL1fi1,…,iLφi1,…,iL1xI2hiL1+1…hiL−∑i1=0N1…∑iL1=0NL1fi1,…,iL1,0…,0φi1,…,iL1xI2.In the second step, we can calculate the bounds of the expectation and variance through the following process:(28)μ¯,μ¯=minxI∈x¯,x¯μxI,maxxI∈x¯,x¯μxI,σ¯2,σ¯2=minxI∈x¯,x¯σ2xI,maxxI∈x¯,x¯σ2xI.From the above procedure, we can find that the response of interest with hybrid interval and random variables can be approximated by SQ-IRMAPC in a unified form. In this paper, the maximum value of the mean and variance of response is calculated by using the genetic algorithm.


## 3.1. IRMAPC Expansion

In this section, both the random and interval variables are considered inY=FxR,xI. The function contains random variable xR and an interval variable xI. Hence, the output of the function will have the characteristics of both random and interval variables, and the proposed IRMAPC method will approximate the response by IRMAPC expansion.The IRMAPC expansion for a function can be expressed as(10)F=FxI,xR=∑0≤i1+⋯+iL1+⋯+iL≤nfi1,…,iLφi1,…,iL1xIφiL1+1,…,iLxR,where(11)φi1,…,iL1xI=φi1x1I×⋯×φiL1xL1I,φiL1+1,…,iLxR=φiL1+1xiL1+1R×⋯×φiLxLR.In traditional IRMAPC, the expansion coefficients are calculated by(12)fi1,…,iL=1hi1×⋯×hiL∑j1=1M1…∑jL=1MLFx^I,x^Rφi1,…,iL1x^IφiL1+1,…,iLx^Rw^i1,…,iL,where x^I=x^j1,…,x^jL1 and x^R=x^jL1,…,x^jL denote the Gaussian nodes of interval variables and random variables, respectively, w^j1,…,jL=∏k=1Lw^jk is the Gaussian weight, Mjj=1,2,…,L denotes the number of Gaussian nodes related to the jth variable, x^jk and w^jk denote the jkth Gaussian node and weight related to the kth variable, respectively, and x^jk and w^jk can be calculated through equations (8) and (9).It can be found from equation (12) that the total number of Gaussian points to determine the coefficient is N=∏k=1LMk. Obviously, the total number of Gaussian points will increase exponentially with the increasing number of uncertain parameters, which may lead to tremendous computational cost. In order to improve the computational efficiency of moment-based polynomial chaos expansion for interval and random analysis, the sparse Gaussian quadrature will be introduced to calculate the expansion coefficient.


## 3.2. Determine the Expansion Coefficient by Using Sparse Quadrature

The sparse quadrature is based on the Smolyak algorithm, which has been widely used in the fields of numerical integration and interpolation and image processing. In this section, the basic principles of sparse quadrature for calculating the expansion coefficient will be deduced.A continuous functionFx defined on x∈−1,1 approximated by polynomial chaos expansion can be denoted as Ql1F, where l denotes the retained order of polynomial chaos expansion. According to the nested hierarchical basis principle of the Smolyak algorithm, the difference format of the approximated function is [33](13)Δk1F=Qk1−Qk−11F,Q01F=0.Furthermore, for ad-dimension problem, the approximated function with the order l of the Smolyak algorithm can be contrasted as(14)QldF=∑k≤l+d−1Δk11⊗⋯⊗Δkd1F,where k denotes the sum of the multidimensional indicators k=∑i=1dki and ⊗ expresses the operation of the tensor product. Through the operation of the tensor product, equation (14) can be expressed as(15)Q1dF=∑l+1≤k≤l+d−1l+d−kd−1l+d−kQk11⊗⋯⊗Qkd1F.Therefore, the integration points in the square grids can be defined as(16)U1d=∪l+1≤k≤l+dUk11⊗⋯⊗Ukd1.The number of the integration points based on the sparse grid method is estimated by(17)N≈2ll!dl.The corresponding coefficient of the weights is(18)wk1…kdi1…id=−1l+d−kd−1l+d−kwk1i1⊗⋯⊗wkdid.


## 3.3. Optimal Choice of the Polynomial Basis for IRMAPC Expansion

According to Witteveen and Bijl [8], the accuracy of MAPC expansion by using different polynomial basis varies greatly. Therefore, when the MAPC expansion is extended for interval analysis and hybrid analysis, a key procedure is to determine the optimal polynomial basis related to each random variable and each interval variable. For the random variable, the polynomial basis derived from the moment of the random variable (through equations (2)∼(4)) can be viewed as the optimal random polynomial basis [11]. For the interval variable, there may also exist an optimal interval polynomial basis. In the following text, the optimal choice of polynomial basis of MAPC for interval analysis will be discussed.The main aim of this section is to choose an optimal interval polynomial basis for interval analysis in the framework of moment-based arbitrary orthogonal polynomial expansion. The interval analysis is to find the maximum and minimum value of the response. Theoretically, the accuracy of moment-based arbitrary polynomial interval expansion can be improved by reducing the error of moment-based arbitrary orthogonal polynomial expansion at the maximum and minimum. However, the maximum or minimum may be an arbitrary point in the range of variation of the interval variable. Therefore, when the moment-based arbitrary orthogonal polynomial expansion is used to approximate the response, we can determine the maximum value of the possible error ofPnx over an interval as follows [33]:(19)ε∞=maxa≤x≤bFx−PnxFx=Fx−PnxFx∞,where Pnx denotes the nth order moment-based arbitrary orthogonal polynomial expansion. Among all the nth order moment-based arbitrary orthogonal polynomial expansions, if Pn∗x satisfied(20)Fx−Pn∗xFx∞=minPnx∈Hnε∞=minPnx∈Hnmaxa≤x≤bFx−PnxFx,then Pn∗x can be viewed as the optimal interval polynomial basis. In equation (14), Hn denotes the universal set of all nth order moment-based arbitrary orthogonal polynomial expansion. According to the theory of Weierstrass [33], Pn∗x is available. However, when n≥1, the computation of Pn∗x is rather difficult. In addition, the polynomial basis of Pn∗x varies for different retained order n. Thus, in general, a polynomial basis whose ε∞ is relatively small is used for numerical approximation instead of Pn∗x. In this paper, a moment-based arbitrary polynomial basis whose ε∞ is relatively small will be used as the interval polynomial basis.According to the theory of orthogonal expansion, the polynomial basis depends on its related weight function or moment. In order to determine interval polynomial basis, the effect of the weight function onε∞ will be investigated through a simple numerical example.Example: suppose thaty=ex2, where the range of variation of x is [−1, 1].Different weight functions are considered, including the weight function of Chebyshev expansion; the weight function whose weight is zero at the midpoint or the bounds of the interval [−1, 1] is shown in Figure1; these two types of weight functions are named as convex shape weight function and concave shape weight function, respectively. For clarity, the expression of each weight function is listed in Table 1. ε∞ of the moment-based arbitrary orthogonal polynomial expansion with different weight functions is shown in Figure 2.Figure 1
Three different types of weight functions.Table 1
Expression of different weight functions.
Type of weight functionExpression of weight functionChebyshev expansion1/π1−x2,−1≤x≤1Constant1/2,−1≤x≤1Concave shapeΓα+β+1/2α+β+1Γα+1Γβ+11−ξα1+ξβ,α=β=0.01,−1≤x≤1Convex shapeΓα+β+1/2α+β+1Γα+1Γβ+11−ξα1+ξβ,α=β=12,−1≤x≤1Figure 2
ε∞ of moment-based arbitrary polynomial interval expansion with different weight functions.It can be found from Figure2 that ε∞ of the moment-based arbitrary orthogonal polynomial expansion is relatively large with the convex or concave shape weight function. The main reason may be that the moment-based arbitrary orthogonal polynomial expansion will introduce large error at the domain where the weight is approaching to zero [32]. As a comparison, ε∞ of the moment-based arbitrary orthogonal polynomial expansion with the Chebyshev weight function or constant weight function is small. In addition, it can be found from Figure 2 that ε∞ of the moment-based arbitrary orthogonal polynomial expansion with the Chebyshev weight function is relatively small than that with the constant weight function. Based on the above analysis, the Chebyshev weight function will be used as the weight function of interval polynomial basis, which can be expressed as(21)wIx=1π1−x2.From the application of MAPC for random analysis, we can find that when the moment is set asμk=∫−11xkρxdx, the polynomial can be orthogonal ρx. It indicates that the polynomial basis will be orthogonal to the weight function of Chebyshev polynomial if the moment is determined by(22)μk=∫−11xkwCx=∫−11xk1π1−x2dx,where wCx=1/π1−x2 is the weight function of the Chebyshev polynomial. Therefore, the moment for each interval variable can be determined according to equation (22).Once the moment is determined, the polynomial basis for interval analysis can then be calculated through equations (2)∼(4).


## 3.4. Calculation of the Bounds of Statistics

To calculate the bounds of statistics, there are two main steps. In the first step, the interval variables are regarded as constant parameters, and the MAPC expansion can be rewritten as the following form:(23)FxI,xR=∑iL1+1=0NL1+1…∑iL=0NL∑i1=0N1…∑iL=0NL1fi1,…,iLφi1,…,iL1xIφiL1+1,…,iLxR,=∑iL1+1=0NL1+1…∑iL=0NLziL1+1,…,iLkφiL1+1,…,iLxR,where(24)ziL1+1,…,iLk=∑i1=0N1…∑iL1=0NL1fi1,…,iLkφi1,…,iL1xI.Due to the orthogonality of the polynomial basis, the expectation and variance of the response can be expressed in terms of the expansion coefficient as follows:(25)μ=E∑iL1+1=0NL1+1…∑iL=0NLziL1+1,…,iLkφiL1+1,…,iLxR=z0,…,0k,σ2=E∑iL1+1=0NL1+1…∑iL=0NLziL1+1,…,iLkφiL1+1,…,iLxR2−μ2,=∑iL1+1=0NL1+1…∑iL=0NLziL1+1,…,iLk2hiL1+1…hiL−z0,…,0k2.Substituting equation (18) into (19) and (20), we can rewrite the expectation and variance of the response as follows:(26)μ=∑i1=0N1…∑iL1=0NL1fi1,…,iL1,0,…,0kφi1,…,iL1xI,(27)σ2=∑iL1+1=0NL1+1…∑iL=0NL∑i1=0N1…∑iL=0NL1fi1,…,iLφi1,…,iL1xI2hiL1+1…hiL−∑i1=0N1…∑iL1=0NL1fi1,…,iL1,0…,0φi1,…,iL1xI2.In the second step, we can calculate the bounds of the expectation and variance through the following process:(28)μ¯,μ¯=minxI∈x¯,x¯μxI,maxxI∈x¯,x¯μxI,σ¯2,σ¯2=minxI∈x¯,x¯σ2xI,maxxI∈x¯,x¯σ2xI.From the above procedure, we can find that the response of interest with hybrid interval and random variables can be approximated by SQ-IRMAPC in a unified form. In this paper, the maximum value of the mean and variance of response is calculated by using the genetic algorithm.


## 4. SQ-IRMAPC for Stiffness Analysis of Air Spring with Hybrid Interval and Random Variables

### 4.1. Finite Element (FE) Model of an Air Spring

Figure3 shows the FE model of an air spring. The Newton–Raphson method is used for calculation, in which the air spring action is divided into many load-increment steps [34]. At the end of each loaded-increment step, the approximate equation is established. This analysis is completed by using ABAQUS. In ABAQUS, the fluid element named FAX2 is used to simulate the air, while the bellows are simulated with the rebar-reinforced surface element and 4-node bilinear axisymmetric solid element. There are 24183 CAX4R elements and 1512 SFMAX1 elements in the whole FE model. By using the ABAQUS, the stiffness of the air spring can be obtained. More details related to the calculation can be found in [1].Figure 3
FE model of the air spring.


### 4.2. Stiffness of Air Spring with Interval and Random Variables

This paper will apply the proposed SQ-IRMAPC method for hybrid interval and random analysis of the air spring. For uncertain air spring system, the uncertain parameters whose random moments or PDFs are available can be treated as random variables and denoted byxR, while other uncertain parameters can be treated as interval variables and denoted by xI. Accordingly, the stiffness of the air spring can be denoted as FxR,xI. Based on the SQ-IRMAPC expansion, the stiffness can be approximated as(29)F˜=∑i1=0N1…∑iL=0NLfi1,…,iLφi1,…,iL1xIφiL1+1,…,iLxR,k=1,2,…,Ntot,where(30)fi1,…,iL≈1h1×⋯×hL∑j1=1M1…∑jL=1MLFx^I,x^Rφi1,…,iL1x^IφiL1+1,…,iLx^Rw^i1,…,iL.In the above equation,Fx^I,x^R denotes the stiffness of the air spring at the nodes of sparse quadrature. After the polynomial basis and the expansion coefficient are determined, the bounds of expectation and variance of stiffness can be obtained according to equations (26)∼(28).The main procedure of the SQ-IRMAPC method for stiffness analysis of the air spring with hybrid interval and random uncertainties can be summarized asStep 1: determine the moment of the random variableStep 2: determineμk of the interval variable according to equation (22)Step 3: compute the polynomial basis through equations (2)∼(4)Step 4: produce the Gaussian nodes and weights through equations (8)∼(9)Step 5: calculate stiffness of the air spring at the nodes of sparse quadratureStep 6: obtain the expansion coefficient through equation (30)Step 7: compute the bounds of expectation and variance of stiffness through equations (26)∼(28)


## 4.1. Finite Element (FE) Model of an Air Spring

Figure3 shows the FE model of an air spring. The Newton–Raphson method is used for calculation, in which the air spring action is divided into many load-increment steps [34]. At the end of each loaded-increment step, the approximate equation is established. This analysis is completed by using ABAQUS. In ABAQUS, the fluid element named FAX2 is used to simulate the air, while the bellows are simulated with the rebar-reinforced surface element and 4-node bilinear axisymmetric solid element. There are 24183 CAX4R elements and 1512 SFMAX1 elements in the whole FE model. By using the ABAQUS, the stiffness of the air spring can be obtained. More details related to the calculation can be found in [1].Figure 3
FE model of the air spring.


## 4.2. Stiffness of Air Spring with Interval and Random Variables

This paper will apply the proposed SQ-IRMAPC method for hybrid interval and random analysis of the air spring. For uncertain air spring system, the uncertain parameters whose random moments or PDFs are available can be treated as random variables and denoted byxR, while other uncertain parameters can be treated as interval variables and denoted by xI. Accordingly, the stiffness of the air spring can be denoted as FxR,xI. Based on the SQ-IRMAPC expansion, the stiffness can be approximated as(29)F˜=∑i1=0N1…∑iL=0NLfi1,…,iLφi1,…,iL1xIφiL1+1,…,iLxR,k=1,2,…,Ntot,where(30)fi1,…,iL≈1h1×⋯×hL∑j1=1M1…∑jL=1MLFx^I,x^Rφi1,…,iL1x^IφiL1+1,…,iLx^Rw^i1,…,iL.In the above equation,Fx^I,x^R denotes the stiffness of the air spring at the nodes of sparse quadrature. After the polynomial basis and the expansion coefficient are determined, the bounds of expectation and variance of stiffness can be obtained according to equations (26)∼(28).The main procedure of the SQ-IRMAPC method for stiffness analysis of the air spring with hybrid interval and random uncertainties can be summarized asStep 1: determine the moment of the random variableStep 2: determineμk of the interval variable according to equation (22)Step 3: compute the polynomial basis through equations (2)∼(4)Step 4: produce the Gaussian nodes and weights through equations (8)∼(9)Step 5: calculate stiffness of the air spring at the nodes of sparse quadratureStep 6: obtain the expansion coefficient through equation (30)Step 7: compute the bounds of expectation and variance of stiffness through equations (26)∼(28)


## 5. Numerical Example

### 5.1. Stiffness Analysis of the Air Spring with Interval Uncertainties

In this section, the uncertain parameters of the air spring are assumed as the interval variable. To obtain the range of variation of the uncertain parameter, the material properties of five air springs in different use stages have been tested. Figure4 shows the air spring which was used to determine the uncertain range of the uncertain parameters. Table 1 lists the range of variation of different parameters.Figure 4
An air spring.
(a)(b)The SQ-IRMAPC is used to predict the stiffness of the air spring with interval uncertainties. To investigate the computational efficiency and accuracy of SQ-IRMAPC, the traditional IRAPC [31] is also introduced to calculate the response of the air spring. Due to the tremendous computational cost of the Monte Carlo method, the high-order Legendre expansion method will be used as the reference method in this numerical example [31]. The stiffness of the air spring calculated by using different methods is shown in Figure 4. The initial pressure of the air spring is 0.5∼0.68 MPa. In SQ-IRMAPC and IRAPC, the retained order is two.It can be seen from Figure5 that the upper bound of stiffness is much larger than the lower bound of stiffness, which means the stiffness is very sensitive to the uncertain parameters. Thus, it is necessary to consider the uncertainties when to analyze the stiffness of the air spring.Figure 5
Stiffness of the air spring with different initial pressure.By comparing the results calculated by using different methods, it can be found that both IRAPC and SQ-IRMAPC can achieve high accuracy. The relative error of SQ-IRMAPC is slightly higher than that of IRAPC; the main reason is that larger number of polynomial basis is retained in IRAPC. However, increasing the number of polynomial basis will lead to larger computational burden. The computational time of the SQ-IRMAPC and IRAPC is1.31×106s and 1.07×105s. Therefore, compared with the traditional IRAPC, the proposed SQ-IRMAPC can greatly improve the computational efficiency.


### 5.2. Stiffness Analysis of the Air Spring with Hybrid Interval and Random Uncertainties

In this numerical example, Young’s modulus and the cross section of cord is assumed as the random variable, while the other parameters are assumed as the interval variable shown in Table2. In particular, only the raw statistical data of the random variable is available. The raw statistical data of Young’s modulus and cross section of cord is shown in Figures 6 and 7.Table 2
Range of variation of interval variables for the air spring system.
ParametersE (MPa)Cord angleCross section of cordC01 of rubberC10 of rubberThickness of rubberRange of variation[1657, 1996][40.8, 49.1][5.4, 6.5][0.3, 0.9][0.03, 0.09][6, 8] mmFigure 6
The statistical data of Young’s modulus.Figure 7
The statistical data of cross section of cord.Both the proposed SQ-IRMAPC and IRAPCM are used to calculate the stiffness of the air spring. In SQ-IRMAPC and IRAPCM, the retained order is set as 2. Table3 shows the relative error of the bounds of expectation and variance yielded by the traditional IRAPCM and the SQ-IRMAPC. The reference results are obtained by using the high-order Legendre expansion. In addition, the execution time of different methods is also shown in Table 3.Table 3
Relative error and execution time of different methods.
Initial pressure (MPa)IRAPC (traditional)SQ-IRMAPC (proposed)Lower boundUpper boundExecution time (s)Lower boundUpper boundExecution time (s)0.500.7%1.2%5.74×107s1.1%1.9%4.21×106s0.600.6%1.3%1.0%1.4%0.681.7%2.1%1.4%2.7%From Table3, the relative error yielded by the proposed SQ-IRMAPC and the traditional IRAPC is less than 3%. It indicates that the proposed SQ-IRMAPC can achieve high accuracy for stiffness analysis of the air spring with hybrid interval and random uncertainties. However, it can be found from Table 3 that the execution time of the proposed SQ-IRMAPC is much less than that of the traditional IRAPC. Therefore, from the application of the proposed SQ-IRMAPC for hybrid uncertainty analysis of the air spring, we can also conclude that the proposed SQ-IRMAPC can achieve higher efficiency than the traditional IRAPC.


## 5.1. Stiffness Analysis of the Air Spring with Interval Uncertainties

In this section, the uncertain parameters of the air spring are assumed as the interval variable. To obtain the range of variation of the uncertain parameter, the material properties of five air springs in different use stages have been tested. Figure4 shows the air spring which was used to determine the uncertain range of the uncertain parameters. Table 1 lists the range of variation of different parameters.Figure 4
An air spring.
(a)(b)The SQ-IRMAPC is used to predict the stiffness of the air spring with interval uncertainties. To investigate the computational efficiency and accuracy of SQ-IRMAPC, the traditional IRAPC [31] is also introduced to calculate the response of the air spring. Due to the tremendous computational cost of the Monte Carlo method, the high-order Legendre expansion method will be used as the reference method in this numerical example [31]. The stiffness of the air spring calculated by using different methods is shown in Figure 4. The initial pressure of the air spring is 0.5∼0.68 MPa. In SQ-IRMAPC and IRAPC, the retained order is two.It can be seen from Figure5 that the upper bound of stiffness is much larger than the lower bound of stiffness, which means the stiffness is very sensitive to the uncertain parameters. Thus, it is necessary to consider the uncertainties when to analyze the stiffness of the air spring.Figure 5
Stiffness of the air spring with different initial pressure.By comparing the results calculated by using different methods, it can be found that both IRAPC and SQ-IRMAPC can achieve high accuracy. The relative error of SQ-IRMAPC is slightly higher than that of IRAPC; the main reason is that larger number of polynomial basis is retained in IRAPC. However, increasing the number of polynomial basis will lead to larger computational burden. The computational time of the SQ-IRMAPC and IRAPC is1.31×106s and 1.07×105s. Therefore, compared with the traditional IRAPC, the proposed SQ-IRMAPC can greatly improve the computational efficiency.


## 5.2. Stiffness Analysis of the Air Spring with Hybrid Interval and Random Uncertainties

In this numerical example, Young’s modulus and the cross section of cord is assumed as the random variable, while the other parameters are assumed as the interval variable shown in Table2. In particular, only the raw statistical data of the random variable is available. The raw statistical data of Young’s modulus and cross section of cord is shown in Figures 6 and 7.Table 2
Range of variation of interval variables for the air spring system.
ParametersE (MPa)Cord angleCross section of cordC01 of rubberC10 of rubberThickness of rubberRange of variation[1657, 1996][40.8, 49.1][5.4, 6.5][0.3, 0.9][0.03, 0.09][6, 8] mmFigure 6
The statistical data of Young’s modulus.Figure 7
The statistical data of cross section of cord.Both the proposed SQ-IRMAPC and IRAPCM are used to calculate the stiffness of the air spring. In SQ-IRMAPC and IRAPCM, the retained order is set as 2. Table3 shows the relative error of the bounds of expectation and variance yielded by the traditional IRAPCM and the SQ-IRMAPC. The reference results are obtained by using the high-order Legendre expansion. In addition, the execution time of different methods is also shown in Table 3.Table 3
Relative error and execution time of different methods.
Initial pressure (MPa)IRAPC (traditional)SQ-IRMAPC (proposed)Lower boundUpper boundExecution time (s)Lower boundUpper boundExecution time (s)0.500.7%1.2%5.74×107s1.1%1.9%4.21×106s0.600.6%1.3%1.0%1.4%0.681.7%2.1%1.4%2.7%From Table3, the relative error yielded by the proposed SQ-IRMAPC and the traditional IRAPC is less than 3%. It indicates that the proposed SQ-IRMAPC can achieve high accuracy for stiffness analysis of the air spring with hybrid interval and random uncertainties. However, it can be found from Table 3 that the execution time of the proposed SQ-IRMAPC is much less than that of the traditional IRAPC. Therefore, from the application of the proposed SQ-IRMAPC for hybrid uncertainty analysis of the air spring, we can also conclude that the proposed SQ-IRMAPC can achieve higher efficiency than the traditional IRAPC.


## 6. Conclusion

In this paper, the interval and random model is introduced to deal with the uncertainties in the air spring. To efficiently calculate stiffness of the air spring system with interval and random variables, a new interval and random polynomial chaos method named as Sparse Quadrature-based Interval and Random Moment Arbitrary Polynomial Chaos (SQ-IRMAPC) is proposed. In SQ-IRMAPC, the moment-based arbitrary polynomial chaos is used to approximate the response of interest. In particular, the optimal choice of the moment-based polynomial for interval analysis has been investigated. The expansion coefficient is calculated by using the sparse quadrature. The SQ-IRMAPC has been applied for stiffness of the air spring with uncertainties, and main conclusions include the following:(1)
The optimal choice of the moment-based polynomial can improve the accuracy of the moment-based polynomial chaos expansion for interval analysis(2)
The stiffness of the air spring is very sensitive to uncertain parameters; thus, the uncertainties should be considered in the stiffness analysis of the air spring(3)
Compared to the traditional interval and random arbitrary polynomial chaos method which is based on the Gaussian quadrature, the sparse quadrature-based interval and random arbitrary polynomial chaos method can greatly improve the computational efficiency



---

*Source: 5592559-2021-08-26.xml*