# Vibration Analysis of Driving-Point System with Uncertainties Using Polynomial Chaos Expansion

Authors: Bin Xiao; Yu-long Zhou; Chao Gao; Shuang-Xia Shi; Ze-Xi Sun; Zong-biao Song
Journal: Shock and Vibration
 (2020)

---

## Abstract

A vibration transfer analysis method based on polynomial chaos expansion (PCE) is proposed in this study and is used to analyze the stochastic dynamic compliance of uncertain systems with the Gaussian distribution. The random dynamic compliance is established by utilizing mode superposition on the system as the parameters of system uncertainties are regarded as input variables. Considering the asymptotic probability density function of mode shape, the dynamic compliance is decomposed into the mean of mode shape and the subsystem represented as an orthogonal polynomial expansion. Following this, the vibration transmission analysis approach is proposed for the random vibration. Results of a numerical simulation carried out employing the PCE approach show that broad-band spectrum analysis is more effective than narrow-band spectrum analysis because the former jump of the dynamic compliance amplitude is weakened. This proposed approach is valid and feasible, but since broad-band spectrum analysis loses some important information about the random vibration, both the aforementioned processes need to simultaneously be applied to analyze the random vibration transmission of low-medium frequency systems.


---

## Body

## 1. Introduction

The research of random dynamic system uncertainty originated from the stochastic finite element method (SFEM) applies the Monte Carlo simulation (MCS) method to the structural finite element analysis (FEA) process and repeats FEA to the model through specific samples [1, 2]; consequently, the characteristics of the statistical distribution of the random response are obtained. For determined structural dynamic systems, extensive research regarding vibration analysis has been carried out, such as decoupling analysis on the nonlinear system [3, 4], analysis on vibration characteristics and response of structures, and study on the vibration transmission path.Furthermore, the roles of various boundary conditions [5, 6], for example, classical restraints, elastic supports and their combinations, and geometric properties [7], such as, cylindrical, conical, and spherical, have been taken into account in the research process. It is thus of great interest to understand thoroughly the vibration behaviors of the system. After more accurate and applicable methods for determined vibration analysis are developed, perturbation theory is applied to the transformation equation, including stochastic distribution parameters, stochastic elastic boundary support structures, and stochastic boundary condition problems, and then the stochastic problem is transformed into a series of deterministic differential equations and boundary conditions, and finally, the random vibration analysis is realized. Furthermore, in order to obtain the mean and standard deviations of the results, the uncertainty of parameters and nonparameters should be considered. Fortunately, the stochastic perturbation method is devoted to the theoretical aspects and computational implementation of the generalized stochastic perturbation technique. It is based on any-order Taylor expansions of random variables and enables for determination of up to fourth-order probabilistic moments and characteristics of the physical system response, which provides grounding in the basic elements of statistics and probability and reliability engineering [8]. Therefore, Huang et al. [9] obtained the approximate statistical expressions of stochastic eigenvalues. Wu and Law [10] discussed the statistical prediction of the dynamic response of bridge structures under random moving loads. Xiao et al. [11] analyzed the random vibration characteristics of the Levinson beam on random elastic constraint conditions. Szafran et al. [12] presented a reliability estimation procedure for steel lattice telecommunication towers based on tensioned joint reliability and used the first-order and second-order reliability methods to confirm the elastic-plastic range of joint reliability. Handa and Anderson [13] used the first- and second-order perturbation of the mean of stochastic variables to obtain statistical properties of the response when analyzing the fluctuation of stochastic variables in the system. On this basis, Chen et al. [14] proposed interval variables to describe various uncertain parameters of the system. Adhikari [15] presented variable matrix distributions of the mass, stiffness, and damping of the quantitative linear structure dynamics under uncertainty based on the optimal stochastic matrix theory [16]. Typically, the statistics of the system are obtained via MCS, SFEM, or polynomial chaos expansion (PCE). MCS is often used to obtain reference results. Although simulation techniques can be used for a wide range of structural dynamics problems, it is computationally intractable, particularly for large-scale problems. It is equally evident that while numerical techniques based on Taylor series expansion are efficient, they do not yield sufficient accuracy when there are significant uncertainties.Due to its accuracy and ability to provide powerful representations of vibration transfer capabilities of nonembedded PCE technology, it has become increasingly popular among researchers and has evolved into an attractive approach. Within the framework of orthogonal polynomial expansion, Pascual and Adhikari [17] proposed four kinds of solution forms: Rayleigh quotient, power method, inverse power method, and eigenvalue equation based on the mixed perturbation PCE method. Sarsri et al. [18] obtained the first two statistical moments of the frequency transfer function through chaos expansion and elaborated the coupling problems of the first- and second-order PCE; furthermore, more than three order statistics of the structural response were analyzed by Xiao and Li [19]. Sepahvand et al. [20] used the truncated PCE of arbitrary random basis to express the uncertain elastic modulus, eigenfrequency, and mode of an orthotropic plate in order to study the stochastic free vibration of the plate. Bahmyari et al. [21] combined the meshless Galerkin method with the generalized PCE to analyze the stochastic bending of moderately thick plates with elastic constraint edges. Najlawi et al. [22, 23] proposed a hybrid multiobjective imperialist competitive algorithm (MOICA) and MCS method for the multiobjective robust design optimization of a mechanism and further developed a hybrid MOICA-PCE algorithm to be exploited for robust optimization [24]. In addition, Hadigol and Doostan [25] presented a hybrid sampling method that employs the so-called alphabetic optimality criteria used in the context of the optimal design of experiments in conjunction with coherence-optimal samples, especially for high-order PCE.The main objective of this study was to present an efficient and reliable method for evaluating the vibration transmission and response of a random dynamical system based on PCE of the dynamic compliance when subjected to pure harmonic excitation. For this purpose, mode superposition-based dynamic compliance is introduced for the random dynamical system, which is used in extracting the mean mode when considering random mode elements with the log-normal asymptotic distribution. A new form of PCE is proposed to give the approximate expression of the dynamic compliance, which is combined with embedded Fourier–Hermite polynomial expansion. Thereafter, the proposed approach is established for discussing the vibration transmission and response of the random vibration system.This paper is organized as follows. Section2 briefly reviews the dynamic compliance and orthogonal polynomial expansion theory; thereafter, the PCE approach is presented. In Section 3, numerical simulation is used to demonstrate the PCE approach on both processes of broad-band spectrum analysis and narrow-band spectrum analysis, and analysis of random vibration transmission and response is implemented by comparing with MCS. Finally, in Section 4, the concluding remarks are given.


## 2. Theoretical Method

### 2.1. Dynamic Compliance

Let us consider a multidegree-of-freedom linear structural system with mass, damping, and stiffness matricesM, C, and K, respectively. The equations of motion describing the forced vibration of a linear and damped discrete system are(1)M∂2Y∂t+C∂Y∂t+KY=F,where Y is the nodal displacement vector and F is the external harmonic excitation with angular frequency ω.An undamped or proportional damped system is considered in this study, whose matricesM, C, and K are real symmetric. Hence, the dynamic compliance matrix based on the mode superposition is defined by(2)ρiji:xj:y=∑n=1∞ϕnyϕnxjωZ˜n,where ρij is the vibration transmission from the jth node to the ith node and Z˜n [26] is the mode impedance:(3)Z˜n=jMnω−ωn2ω+Mnηnωn2ω,where Mn, ϕn⋅, and ηn are the mode mass, mode shape, and mode damping ratio with respect to the n-th-order natural frequency ωn, respectively.However, since we are interested in the limit asηn⟶0, consider Ne-order modal truncation; the real dynamic compliance of equation (2) for x=y,(4)ρiii:x=∑n=1Neϕn2xjωZ˜n,may become arbitrarily large ρiii:x⟶+∞ only if the external excitation frequency ω becomes very close to some natural frequency ωn.The weight of eigenfunctions,(5)w=ϕn2x,is defined on 0≤w≤1. The eigenvector of the matrix, ϕnx, at driving point x has to satisfy the orthogonality of eigenfunctions. A convenient way to characterize the statistics of w is via the moments of arbitrary index q>0:(6)Yq=∫01dwΠNwqw,where ΠN⋅ is the probability density function (PDF), and N is the dimension of the square matrices of the system. As a consequence, the weight w is a random variable distributed in a Gaussian orthogonal ensemble (GOE) which obeys log-normal distribution ∏N [27]. The authors in [27] have given the correlation between the log-normal distribution and the multifractal algorithm as follows:(7)ΠNw=1w2πσ2lnNelnw+α0lnN2/2σ2lnN,α≡−lnwlnN.Hence, whenα=α0, the equation holds: the spectrum of fractal dimensions gα0=1; w¯=w satisfies probability distribution ΠNw¯, i.e.,(8)ΠNw¯=2w¯2πσ2lnNe−2lnw¯+α0lnN2/2σ2lnN.If the external excitation frequencyω becomes very close to some natural frequency ωn, then the corresponding term dominates over the others in equation (4), and the dynamic compliance can be approximated by this biggest term; therefore, in such case, the probability distribution of the dynamic compliance is presented as follows:(9)Pρ˜iii:x≃ρ˜iii:x⟶+∞ρωρ˜iii:x2.Otherwise, the dynamic compliance can be approximated by(10)ρ˜iii:x≃Y1∑n=1Neϕn2xjωZ˜nϕn2x≃Y1∑n=1Neln,where Y1 is the mean of w¯ when all the random variables are mean values [28], with Y1=Ew¯⋅ln being the polynomial function.


### 2.2. Orthogonal Polynomial Expansion

#### 2.2.1. Fourier–Hermite Polynomials

Consider a continuous, differentiable, real-valued, multivariate functionyx that is defined on x=x1,x2,…,xMT∈ℝM, where ℝM is an M-dimensional real vector space. A dimensional decomposition of ln=yx, described by(11)yx=y0+∑S=1M∑1≤a1<⋯<aS≤Mya1…aSxa1,…,xaS,can be viewed as a finite hierarchical expansion of an output function in terms of its input variables with increasing dimensions, where y0 is a constant representing the mean response of the transfer function, and ya1…aSxa1,…,xaS is an S-variate component function quantifying the cooperative effects of S input variables xa1,…,xaS.In order to generate a hierarchical and convergent sequence of approximations ofyx, considering input variables following I.I.D Gaussian distribution, a general S-variate approximation of yx is obtained [29]:(12)y˜Sx≃y0+∑S=1M∑1≤a1<⋯<aS≤M∑bS=1m…∑b1=1mCa1…aSb1…bS∏k=1Sψbkxak,where(13)Ca1…asb1…bS=∫Asya1…asxa1,…,xaS∏k=1Sψbkxakfkxkdxk=∫AMyx∏k=1SψbkxakfXxdxis a coefficient associated with the product of b1 through bS basis functions expressed in terms of xa1,…,xaS, where xk follows probability density function (PDF) fkxk, X follows the joint PDF fXx, and ψbxa is the b-th order Hermite polynomial, where b≤m.


#### 2.2.2. Dimension-Reduction Integration for Calculating Coefficients

Following an early idea by Xu and Rahman [30], a lower-variate approximation of the M-variate function yx is considered, which leads to a reduction in the dimensionality of the M-dimensional integral coefficients of the equation. Using Xu and Rahman’s multivariate function theorem [30], it can be shown that a special R-variate approximation of yx, defined by(14)yRx=∑k=0R−1kM−R+k−1k∑1≤k1<⋯<kR−k≤Myxk1,xk2,…,xkR−k,consists of all terms of the Taylor series of yx that have less than or equal to R variables, where yxk1,xk2,…,xkR−k represents the R−kth dimensional component function of yx, with R<M and k=0,…,R.Substituting equation (14) into equation (13), the polynomial expansion coefficients become(15)Ca1…aSb1…bS≃∑k=0R−1kM−R+k−1k∑k1,…,kR−k=1;k1<⋯<kR−kM∫AR−kyxk1,xk2,…,xkR−k∏m=1Sψbmxam∏m=1R−kfkmxkmdxkm,and the mean is(16)y0≃∑k=0R−1kM−R+k−1k∑k1,…,kR−k=1;k1<⋯<kR−kM∫AR−kyxk1,xk2,…,xkR−k∏m=1R−kfkmxkmdxkm.It is obvious that introducing the dimension-reduction integration for calculating coefficients is substantially simpler and more efficient than performing oneM-dimensional integral, particularly when R≪M.


#### 2.2.3. Gauss–Hermite Integration

The Gauss–Hermite multiple quadrature formulas for Hermite polynomials are(17)∫−∞+∞…∫−∞+∞fxa1,…,xaSe−∑k=1Sxk2∏k=1Sdxk=∑aS=1QS…∑a1=1Q1∏k=1SWkfxa1,…,xaS,where Q1,…,QS are the integral points and associated coefficient Wk depends on the probability distribution fxa1,…,xaS of xa1,…,xaS. Considering a special term e−x2, equation (17) can be transformed into the summation form and used to solve the polynomial coefficients and mean values efficiently. Correspondingly, in order to improve its calculation accuracy and efficiency, we can define Q=Q1=⋯QS.


### 2.3. Error Analysis

#### 2.3.1. Mean

For the driving-point vibration analysis, i.e.,x=y, the expectation of the dynamic compliance in equation (4) is obtained:(18)Eρiii:x=E∑n=1Neϕn2xjωZ˜n=∑n=1NeE1jωZ˜nEϕn2x.And the expectation of the dynamic compliance in equation (9) is(19)Eρ˜iii:x≃Y1⋅E∑n=1Neϕn2xjωZ˜nϕn2x≃Y1∑n=1NeE1jωZ˜nEϕn2xϕn2x.Therefore, equation (19) becomes(20)Eρ˜iii:x=∑n,m=1NeE1jωZ˜nEϕnxϕmx.Comparing equation (19) with equation (20) and considering the relationship ∑n=1Neϕn2x≃∑n=1Neϕn2x2, the following relation can be obtained:(21)Eρ˜iii:x=Eρiii:x.


#### 2.3.2. Standard Error

In the statistical analysis procedure of random variablesx1,x2,…,xm, the standard error (SE) is defined as [18](22)Dr.m.sx=1m∑j=1mxj−Exj2.Supposing a quantity is a function of the measured variables, i.e.,y=fx1,x2,…,xm, the error propagation equation is obtained as follows:(23)Dr.m.s2y=∑j=1m∂f∂xjDr.m.sxj2.By using equations (23) and (24), the SE of the dynamic compliance in equation (4) can be obtained:(24)Dr.m.s2ρiii:x=∑n=1Ne2jωZ˜n2ϕn2xDr.m.s2ϕnx+−1jωZ˜n22ϕn4xDr.m.s2Z˜n.And then, the SE of dynamic compliance (10) can be obtained as follows:(25)Dr.m.s2ρ˜iii:x=∑n=1Ne1jωZ˜n2ϕn2x2Y1ϖ+∂Y1∂ϖDr.m.s2ϕnx+−ϖjωZ˜n22Y12Dr.m.s2Z˜n.Using equations (24) and (25), the SE level can be defined as(26)eLd=10lgDr.m.s2ρiii:x,e˜Ld=10lgDr.m.s2ρ˜iii:x.Then, the statistic of the PCE-based dynamic compliance can be yielded as the format of the SE level, i.e.,(27)ΔeLd=e˜Ld−eLd=10lgDr.m.s2ρ˜iii:xDr.m.s2ρiii:x.It follows that the proposed PCE approach is reasonable and feasible, by equation (21), but its error in equation (27) is inevitable, according to equations (24) and (25). So, the curve shape distortion (CSD) is defined based on equation (21) as(28)CSD=∫−∞+∞pMCSx−pPCExdx,where pMCSx and pPCEx are the PDFs of the SE level by employing MCS and PCE.


## 2.1. Dynamic Compliance

Let us consider a multidegree-of-freedom linear structural system with mass, damping, and stiffness matricesM, C, and K, respectively. The equations of motion describing the forced vibration of a linear and damped discrete system are(1)M∂2Y∂t+C∂Y∂t+KY=F,where Y is the nodal displacement vector and F is the external harmonic excitation with angular frequency ω.An undamped or proportional damped system is considered in this study, whose matricesM, C, and K are real symmetric. Hence, the dynamic compliance matrix based on the mode superposition is defined by(2)ρiji:xj:y=∑n=1∞ϕnyϕnxjωZ˜n,where ρij is the vibration transmission from the jth node to the ith node and Z˜n [26] is the mode impedance:(3)Z˜n=jMnω−ωn2ω+Mnηnωn2ω,where Mn, ϕn⋅, and ηn are the mode mass, mode shape, and mode damping ratio with respect to the n-th-order natural frequency ωn, respectively.However, since we are interested in the limit asηn⟶0, consider Ne-order modal truncation; the real dynamic compliance of equation (2) for x=y,(4)ρiii:x=∑n=1Neϕn2xjωZ˜n,may become arbitrarily large ρiii:x⟶+∞ only if the external excitation frequency ω becomes very close to some natural frequency ωn.The weight of eigenfunctions,(5)w=ϕn2x,is defined on 0≤w≤1. The eigenvector of the matrix, ϕnx, at driving point x has to satisfy the orthogonality of eigenfunctions. A convenient way to characterize the statistics of w is via the moments of arbitrary index q>0:(6)Yq=∫01dwΠNwqw,where ΠN⋅ is the probability density function (PDF), and N is the dimension of the square matrices of the system. As a consequence, the weight w is a random variable distributed in a Gaussian orthogonal ensemble (GOE) which obeys log-normal distribution ∏N [27]. The authors in [27] have given the correlation between the log-normal distribution and the multifractal algorithm as follows:(7)ΠNw=1w2πσ2lnNelnw+α0lnN2/2σ2lnN,α≡−lnwlnN.Hence, whenα=α0, the equation holds: the spectrum of fractal dimensions gα0=1; w¯=w satisfies probability distribution ΠNw¯, i.e.,(8)ΠNw¯=2w¯2πσ2lnNe−2lnw¯+α0lnN2/2σ2lnN.If the external excitation frequencyω becomes very close to some natural frequency ωn, then the corresponding term dominates over the others in equation (4), and the dynamic compliance can be approximated by this biggest term; therefore, in such case, the probability distribution of the dynamic compliance is presented as follows:(9)Pρ˜iii:x≃ρ˜iii:x⟶+∞ρωρ˜iii:x2.Otherwise, the dynamic compliance can be approximated by(10)ρ˜iii:x≃Y1∑n=1Neϕn2xjωZ˜nϕn2x≃Y1∑n=1Neln,where Y1 is the mean of w¯ when all the random variables are mean values [28], with Y1=Ew¯⋅ln being the polynomial function.


## 2.2. Orthogonal Polynomial Expansion

### 2.2.1. Fourier–Hermite Polynomials

Consider a continuous, differentiable, real-valued, multivariate functionyx that is defined on x=x1,x2,…,xMT∈ℝM, where ℝM is an M-dimensional real vector space. A dimensional decomposition of ln=yx, described by(11)yx=y0+∑S=1M∑1≤a1<⋯<aS≤Mya1…aSxa1,…,xaS,can be viewed as a finite hierarchical expansion of an output function in terms of its input variables with increasing dimensions, where y0 is a constant representing the mean response of the transfer function, and ya1…aSxa1,…,xaS is an S-variate component function quantifying the cooperative effects of S input variables xa1,…,xaS.In order to generate a hierarchical and convergent sequence of approximations ofyx, considering input variables following I.I.D Gaussian distribution, a general S-variate approximation of yx is obtained [29]:(12)y˜Sx≃y0+∑S=1M∑1≤a1<⋯<aS≤M∑bS=1m…∑b1=1mCa1…aSb1…bS∏k=1Sψbkxak,where(13)Ca1…asb1…bS=∫Asya1…asxa1,…,xaS∏k=1Sψbkxakfkxkdxk=∫AMyx∏k=1SψbkxakfXxdxis a coefficient associated with the product of b1 through bS basis functions expressed in terms of xa1,…,xaS, where xk follows probability density function (PDF) fkxk, X follows the joint PDF fXx, and ψbxa is the b-th order Hermite polynomial, where b≤m.


### 2.2.2. Dimension-Reduction Integration for Calculating Coefficients

Following an early idea by Xu and Rahman [30], a lower-variate approximation of the M-variate function yx is considered, which leads to a reduction in the dimensionality of the M-dimensional integral coefficients of the equation. Using Xu and Rahman’s multivariate function theorem [30], it can be shown that a special R-variate approximation of yx, defined by(14)yRx=∑k=0R−1kM−R+k−1k∑1≤k1<⋯<kR−k≤Myxk1,xk2,…,xkR−k,consists of all terms of the Taylor series of yx that have less than or equal to R variables, where yxk1,xk2,…,xkR−k represents the R−kth dimensional component function of yx, with R<M and k=0,…,R.Substituting equation (14) into equation (13), the polynomial expansion coefficients become(15)Ca1…aSb1…bS≃∑k=0R−1kM−R+k−1k∑k1,…,kR−k=1;k1<⋯<kR−kM∫AR−kyxk1,xk2,…,xkR−k∏m=1Sψbmxam∏m=1R−kfkmxkmdxkm,and the mean is(16)y0≃∑k=0R−1kM−R+k−1k∑k1,…,kR−k=1;k1<⋯<kR−kM∫AR−kyxk1,xk2,…,xkR−k∏m=1R−kfkmxkmdxkm.It is obvious that introducing the dimension-reduction integration for calculating coefficients is substantially simpler and more efficient than performing oneM-dimensional integral, particularly when R≪M.


### 2.2.3. Gauss–Hermite Integration

The Gauss–Hermite multiple quadrature formulas for Hermite polynomials are(17)∫−∞+∞…∫−∞+∞fxa1,…,xaSe−∑k=1Sxk2∏k=1Sdxk=∑aS=1QS…∑a1=1Q1∏k=1SWkfxa1,…,xaS,where Q1,…,QS are the integral points and associated coefficient Wk depends on the probability distribution fxa1,…,xaS of xa1,…,xaS. Considering a special term e−x2, equation (17) can be transformed into the summation form and used to solve the polynomial coefficients and mean values efficiently. Correspondingly, in order to improve its calculation accuracy and efficiency, we can define Q=Q1=⋯QS.


## 2.2.1. Fourier–Hermite Polynomials

Consider a continuous, differentiable, real-valued, multivariate functionyx that is defined on x=x1,x2,…,xMT∈ℝM, where ℝM is an M-dimensional real vector space. A dimensional decomposition of ln=yx, described by(11)yx=y0+∑S=1M∑1≤a1<⋯<aS≤Mya1…aSxa1,…,xaS,can be viewed as a finite hierarchical expansion of an output function in terms of its input variables with increasing dimensions, where y0 is a constant representing the mean response of the transfer function, and ya1…aSxa1,…,xaS is an S-variate component function quantifying the cooperative effects of S input variables xa1,…,xaS.In order to generate a hierarchical and convergent sequence of approximations ofyx, considering input variables following I.I.D Gaussian distribution, a general S-variate approximation of yx is obtained [29]:(12)y˜Sx≃y0+∑S=1M∑1≤a1<⋯<aS≤M∑bS=1m…∑b1=1mCa1…aSb1…bS∏k=1Sψbkxak,where(13)Ca1…asb1…bS=∫Asya1…asxa1,…,xaS∏k=1Sψbkxakfkxkdxk=∫AMyx∏k=1SψbkxakfXxdxis a coefficient associated with the product of b1 through bS basis functions expressed in terms of xa1,…,xaS, where xk follows probability density function (PDF) fkxk, X follows the joint PDF fXx, and ψbxa is the b-th order Hermite polynomial, where b≤m.


## 2.2.2. Dimension-Reduction Integration for Calculating Coefficients

Following an early idea by Xu and Rahman [30], a lower-variate approximation of the M-variate function yx is considered, which leads to a reduction in the dimensionality of the M-dimensional integral coefficients of the equation. Using Xu and Rahman’s multivariate function theorem [30], it can be shown that a special R-variate approximation of yx, defined by(14)yRx=∑k=0R−1kM−R+k−1k∑1≤k1<⋯<kR−k≤Myxk1,xk2,…,xkR−k,consists of all terms of the Taylor series of yx that have less than or equal to R variables, where yxk1,xk2,…,xkR−k represents the R−kth dimensional component function of yx, with R<M and k=0,…,R.Substituting equation (14) into equation (13), the polynomial expansion coefficients become(15)Ca1…aSb1…bS≃∑k=0R−1kM−R+k−1k∑k1,…,kR−k=1;k1<⋯<kR−kM∫AR−kyxk1,xk2,…,xkR−k∏m=1Sψbmxam∏m=1R−kfkmxkmdxkm,and the mean is(16)y0≃∑k=0R−1kM−R+k−1k∑k1,…,kR−k=1;k1<⋯<kR−kM∫AR−kyxk1,xk2,…,xkR−k∏m=1R−kfkmxkmdxkm.It is obvious that introducing the dimension-reduction integration for calculating coefficients is substantially simpler and more efficient than performing oneM-dimensional integral, particularly when R≪M.


## 2.2.3. Gauss–Hermite Integration

The Gauss–Hermite multiple quadrature formulas for Hermite polynomials are(17)∫−∞+∞…∫−∞+∞fxa1,…,xaSe−∑k=1Sxk2∏k=1Sdxk=∑aS=1QS…∑a1=1Q1∏k=1SWkfxa1,…,xaS,where Q1,…,QS are the integral points and associated coefficient Wk depends on the probability distribution fxa1,…,xaS of xa1,…,xaS. Considering a special term e−x2, equation (17) can be transformed into the summation form and used to solve the polynomial coefficients and mean values efficiently. Correspondingly, in order to improve its calculation accuracy and efficiency, we can define Q=Q1=⋯QS.


## 2.3. Error Analysis

### 2.3.1. Mean

For the driving-point vibration analysis, i.e.,x=y, the expectation of the dynamic compliance in equation (4) is obtained:(18)Eρiii:x=E∑n=1Neϕn2xjωZ˜n=∑n=1NeE1jωZ˜nEϕn2x.And the expectation of the dynamic compliance in equation (9) is(19)Eρ˜iii:x≃Y1⋅E∑n=1Neϕn2xjωZ˜nϕn2x≃Y1∑n=1NeE1jωZ˜nEϕn2xϕn2x.Therefore, equation (19) becomes(20)Eρ˜iii:x=∑n,m=1NeE1jωZ˜nEϕnxϕmx.Comparing equation (19) with equation (20) and considering the relationship ∑n=1Neϕn2x≃∑n=1Neϕn2x2, the following relation can be obtained:(21)Eρ˜iii:x=Eρiii:x.


### 2.3.2. Standard Error

In the statistical analysis procedure of random variablesx1,x2,…,xm, the standard error (SE) is defined as [18](22)Dr.m.sx=1m∑j=1mxj−Exj2.Supposing a quantity is a function of the measured variables, i.e.,y=fx1,x2,…,xm, the error propagation equation is obtained as follows:(23)Dr.m.s2y=∑j=1m∂f∂xjDr.m.sxj2.By using equations (23) and (24), the SE of the dynamic compliance in equation (4) can be obtained:(24)Dr.m.s2ρiii:x=∑n=1Ne2jωZ˜n2ϕn2xDr.m.s2ϕnx+−1jωZ˜n22ϕn4xDr.m.s2Z˜n.And then, the SE of dynamic compliance (10) can be obtained as follows:(25)Dr.m.s2ρ˜iii:x=∑n=1Ne1jωZ˜n2ϕn2x2Y1ϖ+∂Y1∂ϖDr.m.s2ϕnx+−ϖjωZ˜n22Y12Dr.m.s2Z˜n.Using equations (24) and (25), the SE level can be defined as(26)eLd=10lgDr.m.s2ρiii:x,e˜Ld=10lgDr.m.s2ρ˜iii:x.Then, the statistic of the PCE-based dynamic compliance can be yielded as the format of the SE level, i.e.,(27)ΔeLd=e˜Ld−eLd=10lgDr.m.s2ρ˜iii:xDr.m.s2ρiii:x.It follows that the proposed PCE approach is reasonable and feasible, by equation (21), but its error in equation (27) is inevitable, according to equations (24) and (25). So, the curve shape distortion (CSD) is defined based on equation (21) as(28)CSD=∫−∞+∞pMCSx−pPCExdx,where pMCSx and pPCEx are the PDFs of the SE level by employing MCS and PCE.


## 2.3.1. Mean

For the driving-point vibration analysis, i.e.,x=y, the expectation of the dynamic compliance in equation (4) is obtained:(18)Eρiii:x=E∑n=1Neϕn2xjωZ˜n=∑n=1NeE1jωZ˜nEϕn2x.And the expectation of the dynamic compliance in equation (9) is(19)Eρ˜iii:x≃Y1⋅E∑n=1Neϕn2xjωZ˜nϕn2x≃Y1∑n=1NeE1jωZ˜nEϕn2xϕn2x.Therefore, equation (19) becomes(20)Eρ˜iii:x=∑n,m=1NeE1jωZ˜nEϕnxϕmx.Comparing equation (19) with equation (20) and considering the relationship ∑n=1Neϕn2x≃∑n=1Neϕn2x2, the following relation can be obtained:(21)Eρ˜iii:x=Eρiii:x.


## 2.3.2. Standard Error

In the statistical analysis procedure of random variablesx1,x2,…,xm, the standard error (SE) is defined as [18](22)Dr.m.sx=1m∑j=1mxj−Exj2.Supposing a quantity is a function of the measured variables, i.e.,y=fx1,x2,…,xm, the error propagation equation is obtained as follows:(23)Dr.m.s2y=∑j=1m∂f∂xjDr.m.sxj2.By using equations (23) and (24), the SE of the dynamic compliance in equation (4) can be obtained:(24)Dr.m.s2ρiii:x=∑n=1Ne2jωZ˜n2ϕn2xDr.m.s2ϕnx+−1jωZ˜n22ϕn4xDr.m.s2Z˜n.And then, the SE of dynamic compliance (10) can be obtained as follows:(25)Dr.m.s2ρ˜iii:x=∑n=1Ne1jωZ˜n2ϕn2x2Y1ϖ+∂Y1∂ϖDr.m.s2ϕnx+−ϖjωZ˜n22Y12Dr.m.s2Z˜n.Using equations (24) and (25), the SE level can be defined as(26)eLd=10lgDr.m.s2ρiii:x,e˜Ld=10lgDr.m.s2ρ˜iii:x.Then, the statistic of the PCE-based dynamic compliance can be yielded as the format of the SE level, i.e.,(27)ΔeLd=e˜Ld−eLd=10lgDr.m.s2ρ˜iii:xDr.m.s2ρiii:x.It follows that the proposed PCE approach is reasonable and feasible, by equation (21), but its error in equation (27) is inevitable, according to equations (24) and (25). So, the curve shape distortion (CSD) is defined based on equation (21) as(28)CSD=∫−∞+∞pMCSx−pPCExdx,where pMCSx and pPCEx are the PDFs of the SE level by employing MCS and PCE.


## 3. Numerical Simulation

### 3.1. FEM-Based Plate Model

Rectangular plates are widely used in our daily life; all kinds of modeling method have been proposed for the vibration characteristics of rectangular plates with system certainty [31]. In this paper, a rectangular plate is fixed at its left boundary and supported by springs at its right boundary. Considering the vibration transmission and response at a driving point, the coordinate of which is (0.8, 0.4), with single-point harmonic excitation, FEM modeling can be given, shown in Figure 1 with parameters listed in Table 1. The stiffness Ki is the i-th random variable with a total of n=19, that is, Ki∼Nμ,σ2,i=1,2,…,n, follows the PDF N⋅,⋅ of a Gaussian distribution with μ=10000/n and σ=2500/n.Figure 1
Schematic of FEM modeling of the rectangular plate.Table 1
Parameters and random boundary conditions of the rectangular plate.
VariablesDistributionMeanVariation coefficientYoung’s modulusEDeterministic2 × 1011 N/m2—Mass densityρDeterministic7860 kg/m3—Poisson’s ratioμDeterministic0.3—LengthLDeterministic1.2 m—WidthBDeterministic0.6 m—ThicknessδDeterministic7.5 × 10−3 m—Spring stiffness∑i=1nKiGaussian10 k N/m0.25Damping coefficientηnUniform0.5%—In order to demonstrate the efficiency of PCE, some benchmark tests are analyzed with random parameters and vertical vibration. The proposed PCE approach with univariateS=1, bivariate S=2, and trivariate S=3 approximations is elaborated and considered to be the main results of this study. For the sake of assessing the accuracy of the PCE method, a comparison MCS approach is used, with its results using 1×105 samples, serving as reference results. The errors between the estimated and reference responses are defined in equation (27), and CSD is defined in equation (28).Table 2
Natural frequencies of the rectangular plate.
Order,n123456789101112131415161718192021fn (Hz)13.4622.5732.2163.1679.34119.93120.68155.94164.76200.17230.54269.92289.84308.85322.80338.42373.41418.32439.28440.96462.04Figure 2
Mean model of the dynamic compliance.Figure 3
Mean of the dynamic compliance.


### 3.2. Narrow-Band Spectrum Analysis

(a)
Employing mode analysis, the natural frequencies within0,500Hz were obtained, as listed in Table 2. Assuming that random variables are mean values, we obtained the mean model of the dynamic compliance, as shown in Figure 2. For comparison purposes, the mean of the dynamic compliance is given in Figure 3. The curves in Figures 2 and 3 are almost identical; these results coincided with equation (21), i.e., the mean model embodies the mean of the vibration transmission and response of the random system, thus implying that the PCE approach is valid and feasible for narrow-band spectrum analysis. ,Figure 4
Error analysis curves of the dynamic compliance spectral level.Figure 5
PDFs of the dynamic compliance amplitude within 18∼22 Hz. (a) 18 Hz. (b) 19 Hz. (c) 20 Hz. (d) 21 Hz. (e) 22 Hz.
(a)(b)(c)(d)(e)Figure 6
PDFs of the dynamic compliance amplitude within 29∼35 Hz. (a) 29 Hz. (b) 30 Hz. (c) 31 Hz. (d) 32 Hz. (e) 33 Hz. (f) 34 Hz. (g) 35 Hz.
(a)(b)(c)(d)(e)(f)(g)Figure4 shows the error analysis curves of the dynamic compliance spectral level. Results reveal that, for both univariate S=1 and bivariate S=2 approximations of the PCE approach, very small differences exist from the results from MCS. That is to say, the error of the PCE approach cannot be eliminated by using higher S-variate approximation, even if S⟶N. After a brief review of the PCE approach, it is not difficult to find that the majority of errors was caused by mode-shape asymptotic estimation, multirandom variable decoupling, and Hermite polynomial truncation in the orthogonal polynomial expansion. Fortunately, the maximum relative error is less than 3 dB, so the statistical analysis results are excellent in the low-frequency band within 100 Hz and are still acceptable at frequencies more than 100 Hz.(b)
In order to further discuss the detailed error of PCE, to consider, for example, just the two frequency bands of 18∼22 Hz and 29∼35 Hz, the PDF curves of the dynamic compliance amplitude are plotted in Figures5 and 6, respectively. The compared results indicate that there exists CSD, but it is very small difference between the estimated and reference results.Although the stiffness of the spring follows the I.I.D. Gaussian distribution, it is evident that the PDF curves of the dynamic compliance do not obey the Gaussian distribution but embody a trend towards the Gaussian distribution. The dynamic compliance curves are in the damping control area, which are very sensitive to damping near the natural frequency. It can be seen from Figure2 and Table 2 that the dominant mode of vertical vibration (resonant frequency) is not evident within the frequency band of 18∼22 Hz, but is evident within the frequency band of 29∼35 Hz. Correspondingly, a comparison of Figures 5 and 6 reveals that the PDFs have CSD. Random fluctuation of the resonant frequency leads to a large jump in the dynamic compliance amplitude when the mode damping is small; this leads to the PDFs having a larger local slope, causing a larger CSD. It is evident that increasing the mode damping is beneficial to eliminate CSD and, subsequently, to reduce its error.


### 3.3. Broad-Band Spectrum Analysis

In view of the plate-structure FEM modeling shown in Figure1, on the basis of narrow-band spectrum analysis, the vibration transmission and response of the 1/3 octave (Oct) spectrum level are further analyzed using the proposed PCE approach.(a)
Figure7 shows the PDFs of the 1/3 Oct spectral vibration response with center frequency fc=20,31.5 Hz. Results indicate that a higher S-variate approximation is better than the univariate one. There is no vertical vibration dominant mode (resonant frequency) within the frequency band of 20 Hz, whose result is significantly better than that of the 31.5 Hz frequency-band result, which does have a vertical vibration dominant mode (resonant frequency). A comparison of Figure 7 with Figures 5 and 6 implies that the PDF of a broad-band spectrum converges towards a Gaussian distribution, and the jump of the dynamic compliance amplitude is weakened due to the 1/3 Oct broad-band energy superposition, which makes the proposed PCE approach more effective.Figure 7
PDFs of the energy-level response of the 1/3 Oct spectrum. (a)fc=20 Hz. (b)fc=31.5 Hz.
(a)(b)Figure 8
Mean of the energy-level response of the 1/3 Oct spectrum.Figure 9
SE of the energy-level response of the 1/3 Oct spectrum.(b)
Figure8 shows the energy-level response mean of the 1/3 Oct spectrum within the frequency range 20,200Hz. The curves in Figure 8 are almost identical, and their results coincide with equation (21), thus implying that the PCE approach is valid and feasible for the broad-band spectrum.Correspondingly, Figure9 shows the SE of the energy-level response of the 1/3 Oct spectrum. Results show that, for the MCS, very small differences exist from the various PCE results; CSD behaves similarly, with similar reason, to the narrow-band spectrum. However, comparing Figure 9 with Figure 4, the relative error of the broad-band spectrum is significantly smaller than that of the narrow-band spectrum. This shows that the PCE approach for broad-band spectrum analysis is more efficient and accurate. Unfortunately, there is a certain degree of small error in the higher-order statistics, but it does not affect the accuracy of the PCE method.


## 3.1. FEM-Based Plate Model

Rectangular plates are widely used in our daily life; all kinds of modeling method have been proposed for the vibration characteristics of rectangular plates with system certainty [31]. In this paper, a rectangular plate is fixed at its left boundary and supported by springs at its right boundary. Considering the vibration transmission and response at a driving point, the coordinate of which is (0.8, 0.4), with single-point harmonic excitation, FEM modeling can be given, shown in Figure 1 with parameters listed in Table 1. The stiffness Ki is the i-th random variable with a total of n=19, that is, Ki∼Nμ,σ2,i=1,2,…,n, follows the PDF N⋅,⋅ of a Gaussian distribution with μ=10000/n and σ=2500/n.Figure 1
Schematic of FEM modeling of the rectangular plate.Table 1
Parameters and random boundary conditions of the rectangular plate.
VariablesDistributionMeanVariation coefficientYoung’s modulusEDeterministic2 × 1011 N/m2—Mass densityρDeterministic7860 kg/m3—Poisson’s ratioμDeterministic0.3—LengthLDeterministic1.2 m—WidthBDeterministic0.6 m—ThicknessδDeterministic7.5 × 10−3 m—Spring stiffness∑i=1nKiGaussian10 k N/m0.25Damping coefficientηnUniform0.5%—In order to demonstrate the efficiency of PCE, some benchmark tests are analyzed with random parameters and vertical vibration. The proposed PCE approach with univariateS=1, bivariate S=2, and trivariate S=3 approximations is elaborated and considered to be the main results of this study. For the sake of assessing the accuracy of the PCE method, a comparison MCS approach is used, with its results using 1×105 samples, serving as reference results. The errors between the estimated and reference responses are defined in equation (27), and CSD is defined in equation (28).Table 2
Natural frequencies of the rectangular plate.
Order,n123456789101112131415161718192021fn (Hz)13.4622.5732.2163.1679.34119.93120.68155.94164.76200.17230.54269.92289.84308.85322.80338.42373.41418.32439.28440.96462.04Figure 2
Mean model of the dynamic compliance.Figure 3
Mean of the dynamic compliance.


## 3.2. Narrow-Band Spectrum Analysis

(a)
Employing mode analysis, the natural frequencies within0,500Hz were obtained, as listed in Table 2. Assuming that random variables are mean values, we obtained the mean model of the dynamic compliance, as shown in Figure 2. For comparison purposes, the mean of the dynamic compliance is given in Figure 3. The curves in Figures 2 and 3 are almost identical; these results coincided with equation (21), i.e., the mean model embodies the mean of the vibration transmission and response of the random system, thus implying that the PCE approach is valid and feasible for narrow-band spectrum analysis. ,Figure 4
Error analysis curves of the dynamic compliance spectral level.Figure 5
PDFs of the dynamic compliance amplitude within 18∼22 Hz. (a) 18 Hz. (b) 19 Hz. (c) 20 Hz. (d) 21 Hz. (e) 22 Hz.
(a)(b)(c)(d)(e)Figure 6
PDFs of the dynamic compliance amplitude within 29∼35 Hz. (a) 29 Hz. (b) 30 Hz. (c) 31 Hz. (d) 32 Hz. (e) 33 Hz. (f) 34 Hz. (g) 35 Hz.
(a)(b)(c)(d)(e)(f)(g)Figure4 shows the error analysis curves of the dynamic compliance spectral level. Results reveal that, for both univariate S=1 and bivariate S=2 approximations of the PCE approach, very small differences exist from the results from MCS. That is to say, the error of the PCE approach cannot be eliminated by using higher S-variate approximation, even if S⟶N. After a brief review of the PCE approach, it is not difficult to find that the majority of errors was caused by mode-shape asymptotic estimation, multirandom variable decoupling, and Hermite polynomial truncation in the orthogonal polynomial expansion. Fortunately, the maximum relative error is less than 3 dB, so the statistical analysis results are excellent in the low-frequency band within 100 Hz and are still acceptable at frequencies more than 100 Hz.(b)
In order to further discuss the detailed error of PCE, to consider, for example, just the two frequency bands of 18∼22 Hz and 29∼35 Hz, the PDF curves of the dynamic compliance amplitude are plotted in Figures5 and 6, respectively. The compared results indicate that there exists CSD, but it is very small difference between the estimated and reference results.Although the stiffness of the spring follows the I.I.D. Gaussian distribution, it is evident that the PDF curves of the dynamic compliance do not obey the Gaussian distribution but embody a trend towards the Gaussian distribution. The dynamic compliance curves are in the damping control area, which are very sensitive to damping near the natural frequency. It can be seen from Figure2 and Table 2 that the dominant mode of vertical vibration (resonant frequency) is not evident within the frequency band of 18∼22 Hz, but is evident within the frequency band of 29∼35 Hz. Correspondingly, a comparison of Figures 5 and 6 reveals that the PDFs have CSD. Random fluctuation of the resonant frequency leads to a large jump in the dynamic compliance amplitude when the mode damping is small; this leads to the PDFs having a larger local slope, causing a larger CSD. It is evident that increasing the mode damping is beneficial to eliminate CSD and, subsequently, to reduce its error.


## 3.3. Broad-Band Spectrum Analysis

In view of the plate-structure FEM modeling shown in Figure1, on the basis of narrow-band spectrum analysis, the vibration transmission and response of the 1/3 octave (Oct) spectrum level are further analyzed using the proposed PCE approach.(a)
Figure7 shows the PDFs of the 1/3 Oct spectral vibration response with center frequency fc=20,31.5 Hz. Results indicate that a higher S-variate approximation is better than the univariate one. There is no vertical vibration dominant mode (resonant frequency) within the frequency band of 20 Hz, whose result is significantly better than that of the 31.5 Hz frequency-band result, which does have a vertical vibration dominant mode (resonant frequency). A comparison of Figure 7 with Figures 5 and 6 implies that the PDF of a broad-band spectrum converges towards a Gaussian distribution, and the jump of the dynamic compliance amplitude is weakened due to the 1/3 Oct broad-band energy superposition, which makes the proposed PCE approach more effective.Figure 7
PDFs of the energy-level response of the 1/3 Oct spectrum. (a)fc=20 Hz. (b)fc=31.5 Hz.
(a)(b)Figure 8
Mean of the energy-level response of the 1/3 Oct spectrum.Figure 9
SE of the energy-level response of the 1/3 Oct spectrum.(b)
Figure8 shows the energy-level response mean of the 1/3 Oct spectrum within the frequency range 20,200Hz. The curves in Figure 8 are almost identical, and their results coincide with equation (21), thus implying that the PCE approach is valid and feasible for the broad-band spectrum.Correspondingly, Figure9 shows the SE of the energy-level response of the 1/3 Oct spectrum. Results show that, for the MCS, very small differences exist from the various PCE results; CSD behaves similarly, with similar reason, to the narrow-band spectrum. However, comparing Figure 9 with Figure 4, the relative error of the broad-band spectrum is significantly smaller than that of the narrow-band spectrum. This shows that the PCE approach for broad-band spectrum analysis is more efficient and accurate. Unfortunately, there is a certain degree of small error in the higher-order statistics, but it does not affect the accuracy of the PCE method.


## 4. Concluding Remarks

A vibration transmission analysis approach is presented based on the PCE of the random dynamic compliance for a system with uncertainties. The numerical illustrations demonstrated that the proposed PCE approach is efficient in comparison with MCS. Some conclusions were drawn as follows:(a)
Compared to the MCS, the results of narrow- and broad-band spectrum analysis under uncertain dynamical systems based on the proposed PCE approach are reasonable and feasible. The jump of the dynamic compliance amplitude caused by the random fluctuation of resonant frequency is weakened because of broad-band energy superposition.(b)
Compared to MCS, the PDF curves obtained by PCE have CSD, which is a relative error of vibration analysis. The error cannot be eliminated because of modal-modeling asymptotic estimation, even if we increase the order of statistics, multirandom variable decoupling, and Hermite polynomial truncation in the orthogonal polynomial expansion. In particular, random fluctuation of the resonant frequency is generated by the system uncertainty, causing a jump in the dynamic compliance amplitude.(c)
The estimated results will converge towards statistical energy analysis results along with an increase in mode overlap numbers. It is noteworthy that some important information about the narrow-band spectrum will be lost when carrying out broad-band spectrum analysis. Therefore, when applying the PCE approach to a random vibration system of low-medium frequency, it is suggested that narrow- and broad-band spectrum analysis be implemented simultaneously.However, the PDF estimation of the mode shape of the random vibration system, the jump in the dynamic compliance amplitude, and nonlinear problems need to be further investigated. Meanwhile, studies should focus on the PCE approach for the cross-point vibration transmission of the dynamical system. In addition, since the system is discretized into FEA models, whose boundary conditions of stochastic elastic constraints are described by finite number of stiffness, PCE can be used for random vibration analysis of the system in this study. However, the problem that PCE is computationally intractable, particularly for numerous random stiffness, even invalid for continuous elastic restraint, offers many challenges; therefore, the methodology of PCE also needs to be discussed in the future.



---

*Source: 8074351-2020-09-28.xml*