A reduced polynomial chaos expansion model for stochastic analysis of a moving load on beam system with non-Gaussian parameters

S. Q. Wu1 , S. S. Law2

1Department of Engineering Mechanics, School of Civil Engineering, Southeast University, Nanjing, P. R. China

1Jiangsu Key Laboratory of Engineering Mechanics, Southeast University, Nanjing, P. R. China

2Department of Civil and Environmental Engineering, Hong Kong Polytechnic University, Hong Kong, P. R. China

1Corresponding author

Journal of Vibroengineering, Vol. 17, Issue 3, 2015, p. 1560-1577.
Received 21 June 2014; received in revised form 22 August 2014; accepted 29 April 2015; published 15 May 2015

Copyright © 2015 JVE International Ltd. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Creative Commons License
Abstract.

Stochastic analysis of a load moving on a beam is conducted in which the Spectral Stochastic Finite Element Method (SSFEM) is adopted to simulation the uncertainties in the system parameters and excitation forces. Since the dimension of Polynomial Chaos Expansion (PCE) of the beam responses will exponentially grow with the number of K-L components of the system and excitation uncertainties, this limits the application of the SSFEM. A reduced PCE model is proposed in this paper to improve the computational efficiency. The non-Gaussian random variables in the Karhunen-Loéve Expansion (KLE) of the non-Gaussian system parameters are assumed “uncoupled”. Numerical simulations show that the computational effort can significantly be reduced while accurate predictions on the response statistics can still be achieved. Studies on the effect of different level of randomness in the system parameters and excitation forces show that the proposed method has good performance even with a high level of uncertainty.

Keywords: moving load, Karhunen-Loéve expansion, reduced polynomial chaos expansion, non-Gaussian, uncertainty.

1. Introduction

The model of “Moving loads on the beam” is widely adopted in simulating engineering problems such as “moving vehicles on bridge deck”, “planes landing on the runway” and “cargos transport by moving overhead crane”, etc. Traditional methods for analyzing the bridge response under moving vehicles can be mainly categorized into two kinds: (1) methods based on modal superposition technique [1, 2]; (2) methods based on finite element method [3, 4].

All the above methods are deterministic approaches in which uncertainties in the bridge structure and the excitation due to moving vehicle are not considered. However, the bridge-vehicle system often exhibits an inherent randomness in the system parameters such as the elastic modulus, mass density, boundary conditions, etc. Furthermore, the interaction forces between the bridge and vehicles are stochastic due to the irregular road surface profile which is usually defined with a power spectral density according to the ISO standard [5]. The samples of the excitation forces are generated from the power spectral density function which will be different from each other. When different samples of the excitation forces are adopted in conventional deterministic analysis, different response samples will be obtained, and this may lead to inaccurate judgment by engineers. In all, the treatment in deterministic analysis only represents an “approximation” of the actual reality in general due to these inherent uncertainties in the structural properties and boundary conditions as well as in the loading process. Stochastic approaches can give a better description on the bridge-vehicle interaction.

The bridge-vehicle interaction problem has been studied considering the uncertainty in the excitation forces with the Monte Carlo Simulation (MCS) method [6], and with stochastic methods in the time domain [7] and in the frequency domain [8]. Others extended the work by introducing Gaussian uncertainties in the system parameters. Fryba et al. [9] performed the stochastic finite element analysis with first-order perturbation to evaluate the variance of the deflection and bending moment of a beam on foundations with uncertain damping and stiffness under a moving force. Dynamic analysis of an Euler-Bernoulli beam excited by a moving oscillator with random mass, velocity and acceleration has been investigated by Muscolino et al. [10] with the improved perturbation approach. Chang et al. [11, 12] studied the dynamic response of an Euler-Bernoulli beam with uncertain parameters such as random mass, stiffness, damping, etc. which is subjected to a moving oscillator/half-car model with random velocity and acceleration. Wu and Law [13] performed dynamic analysis on the bridge-vehicle interaction problem considering Gaussian uncertainties using Karhunen-Loéve expansion. However, the system parameters are always positive in engineering practice, and it is more appropriate to model them with non-Gaussian assumption. Wu and Law [14] later considered the non-Gaussian system parameters in the bridge-vehicle interaction problem where the full dimension of the PCE was adopted in the uncertainty modeling. Liu et al. [15] treated the uncertain system parameters of a bridge-vehicle model as interval variables and perform an interval dynamic analysis on the bridge-vehicle interaction system with these uncertainties.

Though the Spectral Stochastic Finite Element Method (SSFEM) employed by Wu and Law [14] can improve the accuracy of analysis with perturbation technique when the level of uncertainty increases, the SSFEM suffers from the curse of exponentially increased dimension in the Polynomial Chaos Expansion for the solution when the required number of the Karhunen-Loéve (K-L) components for both the system parameters and excitation is large [16]. In the dynamic analysis of a bridge structure, when the excitation force includes high frequency random fluctuations, e.g. vehicle excitation arising from road surface roughness or the ground motion excitation due to an earthquake, a relative large number of K-L components maybe included, e.g. larger than 100 components as shown in the numerical simulation of this paper. It is noted that in the PCE for a random process, the number of K-L components adopted is usually not larger than twenty [17]. This drawback limits the application of the SSFEM in the case with relatively simple excitation forces.

In this paper, a reduced PCE model is proposed to replace the “full” PCE model in SSFEM. A non-Gaussian stochastic process model including high frequency random fluctuations can be represented efficiently with the proposed reduced model. Accurate response statistics of the beam under the moving loads can also be predicted using the proposed stochastic finite element method.

2. System equation of motion with uncertainty

A two-dimensional simply supported Euler-Bernoulli beam with multiple loads moving on top is shown in Fig. 1. The mass per unit length ρ, elastic modulus E and damping c of the beam structure are assumed as non-Gaussian random processes, with mean value ρ-, E-, c- and standard deviation σρ, σE, σc, respectively. The random part of the these parameters, denoted as ρ~, E~, c~, can be obtained by subtracting the mean value from the random variable.

The equation of motion of the structure with random material properties and randomness in the excitations can then be written as:

(1)
ρ x , θ A 2 t 2 u x , t , θ + c x , θ t u x , t , θ + E x , θ I 4 x 4 u x , t , θ
              = j = 1 N F F j t , θ δ x - v j t ,

where A is the cross-sectional area, I is the moment of inertia of the beam and θ denotes the random dimension. NF is the number of the moving forces. Employing the Hermitian cubic interpolation shape functions and with the assumption of Rayleigh damping for the structure, the equation of motion of the beam can be rewritten in matrix form as:

(2)
M θ R ¨ t , θ + C θ R ˙ t , θ + K θ R t , θ = H b F t , θ ,

where M(θ), C(θ), K(θ) are the stochastic mass, damping and stiffness matrices of the beam model with M(θ)=Md+M~(θ), C(θ)=Cd+C~(θ), K(θ)=Kd+K~(θ), respectively. Md, M~(θ), Cd, C~(θ), Kd, K~(θ) are the deterministic and random components of the system mass, damping and stiffness matrices respectively, and they can be obtained by assembling the corresponding elemental matrices as:

(3)
M d e = l H e T ρ - A H e d l ,       M ~ e θ = l H e T ρ ~ θ A H e d l ,
K d e = l B e T E - I B e d l ,       K ~ e θ = l B e T E ~ θ I B e d l ,

where He and Be are respectively the shape function matrix and strain-displacement matrix for each element. l is the length of beam element. F(t,θ) is the stochastic force vector applied on the beam structure with F=F1F2FNFT. R(t,θ), R˙t,θ and R¨t,θ are the stochastic nodal displacement, velocity and acceleration vectors of the structure respectively with R=R1  R2    RnT. Hb is a n×NF location matrix of the moving forces as shown in Eq. (4), where n is the number of degree-of-freedom of the beam structure after considering the boundary condition. We have:

(4)
H b = 0 0 0 H 1 0 0 0 H 2 0 0 0 0 0 H N F 0 0 0 n × N F ,

where the shape function vector Hi can be written in the global coordinate as:

(5)
H i = 1 - 3 x j t - i - 1 l l 2 + 2 x j t - i - 1 l l 3 x j t - i - 1 l x j t - i - 1 l l - 1 2 3 x j t - i - 1 l l 2 - 2 x j t - i - 1 l l 3 x j t - i - 1 l x j t - i - 1 l l 2 - x j t - i - 1 l l T ,

and xj(t) is the location of the jth force on the ith element at time t with (i-1)lxj(t)<il. The two coefficients cM and cK of the Rayleigh damping:

(6)
C = c M M + c K K ,

are assumed constant, and they can be obtained from the mean values of the first two modal frequencies and damping ratios. The stochastic displacement of the beam at position x and time t has the following relationship with the stochastic nodal displacement as:

(7)
w x , t , θ = H x R t , θ ,

where H(x)=0    HxiT  0    0 and HxiT is the shape function of the beam structure. H(x) is a 1×n vector with zero entries except at the ith beam element on which the force is located.

Fig. 1. The beam-moving load system

 The beam-moving load system

3. Karhunen-Loéve expansion

3.1. Theory [18]

Let u(x,θ) be a second-order random process and u-x denote the expected value of u(x,θ). The covariance function C(x1,x2) of u(x,θ) is bounded, symmetric and positive definite and has the following spectral decomposition:

(8)
C x 1 , x 2 = n = 0 λ n ϕ n x 1 ϕ n x 2 ,

where λn and φn(x) are the eigenvalues and the eigenvectors of the covariance function, respectively. They can be proved to be solution of the following integral equation:

(9)
  C x 1 , x 2 ϕ n x 1 d x 1 = λ n ϕ n x 2 .

Due to symmetry and the positive definiteness of the covariance kernel, the eigenfunctions are orthogonal and they form a complete set representing the covariance function. The eigenvectors can then be normalized according to the following criterion:

(10)
  ϕ i x ϕ j x d x = δ i j ,

where δij is the Kronecker delta. The random process u(x,θ) can then be written as:

(11)
u x , θ = u - x + u ~ x , θ = u - x + n = 1 ξ n θ λ n ϕ n x ,

where ξn(θ) is a group of uncorrelated random variables. As a special case, when u(x,θ) is a Gaussian random process, ξn(θ) is a set of standard Gaussian random variables with the following orthogonal properties:

(12)
E ξ n θ = 0 ,       E ξ i θ ξ j θ = δ i j ,

where E() denotes the expectation value.

4. Polynomial Chaos expansion

4.1. Theory [19]

For an element u(θ) from a set in the space spanned by the orthogonal basis, it has been proved that u(θ) can be written in the following form as:

(13)
u θ = i = 0 p j = 0 q b k Ψ i ξ 0 θ , ξ 1 θ , , ξ j θ ,

where Ψi are the Hermite polynomials of ξ1(θ), ξ2(θ),…, ξj(θ), and Ψ0= 1.0, ξ0(θ)= 1.0. bk are the generalized Fourier coefficients with k= (0, 1,…, N) and N=1+p+q!/p!q!. The series in Eq. (13) can be refined along the random dimension θ either by increasing the total number of random variables q or the maximum order of polynomials p included in the Polynomial Chaos (PC) expansion. The first refinement takes into account higher frequency random fluctuations of the underlying stochastic process, while the second refinement captures strong nonlinear dependence of the solution of this stochastic process [20].

The polynomials Ψk are orthogonal satisfying the relationship:

(14)
Ψ i Ψ j = Ψ i 2 δ i j ,

where denotes the inner product and the value of Ψi2 can be calculated analytically [19].

The generalized Fourier coefficients bk in Eq. (13) can be evaluated from:

(15)
b k = u θ Ψ i ξ 1 θ , , ξ j θ Ψ i 2 ξ 1 θ , , ξ j θ .

4.2. Polynomial Chaos coefficients

Numerous approaches have been proposed to obtain the coefficients in the PCE model for a random process. Field et al. [21] proposed two methods with Monte Carlo simulation and numerical integration, respectively, to evaluate the coefficients. Non-intrusive method and Intrusive method were adopted by Reagan et al [22] for the analysis. The method using Latin Hypercube sampling and Fitting Regression model [23] will be introduced as follows due to its high efficiency and great accuracy, and it will be adopted in this paper to evaluate the coefficients in a one-dimensional PCE in the proposed reduced PCE model.

For a one-dimensional problem, i.e. q in Eq. (13) equals to unity, consider the PCE truncated at the pth level:

(16)
u θ = i = 0 p b i Ψ i ξ θ = b 0 + b 1 ξ + b 2 ξ 2 - 1 + b 3 ξ 3 - 3 ξ + + ε
            = b 0 + b 1 f 1 ξ + b 2 f 2 ξ + + b p f p ξ + ε ,

where bi, i= 0, 1, 2,…, p, are the regression coefficients, and ε is the error of the model equation. Functions f1(ξ)=ξ, f2(ξ)=ξ2-1,…, fp(ξ) are the pth one-dimensional Hermite polynomials of ξ(θ).

It is noted that Eq. (16) can be regarded as a linear regression model. To evaluate the coefficients bi, n sample values, e.g. z1~zn, should be drawn from the standard random variable ξ(θ). With these sample values, the matrix form can be formulated as:

(17)
Y = Z B ^ + e ,

where Y=y1y2ynT, B^=b0b1bqT, e=ε1ε2εnT and:

Z = 1 f 1 z 1 f 2 z 1 f p z 1 1 f 1 z 2 f 2 z 2 f p z 2 1 f 1 z n f 2 z n f p z n ,

and y1 to yn are the sample values drawn from u(θ). Vector e includes the error terms. The sample values of fi(zk) (i= 1, 2,…, p) in Eq. (17) can be obtained from the Cumulative Density Function (CDF) of the random variables.

Using the least-squares method, the regression coefficients can be calculated as:

(18)
B ^ = Z T Z - 1 Z T Y .

The fitted model and the residuals are:

(19)
Y ^ = Z B ^ ,       e = Y - Y ^ .

5. Modeling with non-Gaussian uncertainties

5.1. The reduced PCE model

The KLE of a non-Gaussian process results in a set of deterministic coefficients multiplied by the corresponding uncorrelated non-Gaussian random variables {ξi(θ)}. In the PCE of a non-Gaussian process, these uncorrelated non-Gaussian random variables are projected on the Polynomial Chaos basis. The dimension of PCE, denoted as KR,to represent the stochastic response of a beam structure is determined according to the following [20]:

(20)
K R = 1 + k s + p ! k s ! p ! ,

where p is the maximum order of the Polynomial Chaos in PCE; ks is the number of K-L components required in the representation of the response which depends on the number of K-L components adopted in the representations of excitation forces and system parameters. For example, if the randomness in elastic modulus, mass density and excitation forces are independent and the Rayleigh damping is assumed, then ks=kE+kρ+kF. Variables kE, kρ, kF are the number of K-L components in the KLE of the Gaussian distributed elastic modulus, mass density and the excitation forces, respectively.

When a large number of K-L components is required inthe KLE of the excitation force, an extremely large number of Polynomial Chaos will be included in the PCE of the non-Gaussian random responses according to Eq. (20), and this leads to computation problem due to the limited capability of computer. This drawback limits the application of the SSFEM in the case with relatively simple excitation.

For the case with complex excitation forces, a reduced PCE model is proposed to avoid the curse of dimensionity of PCE in the response representation. The uncorrelated random variables in the KLE of the non-Gaussian random processes are assumed to be “uncoupled”. The spaces spanned by these “uncoupled” random variables are assumed to be mutually orthogonal, which results in keeping the terms of Ψi(ξjθ) in Eq. (13) and eliminating the terms of Ψiξ1θ,ξ2θ,,ξjθ in the new reduced model. For example, when p= 2 and q= 2, Eq. (13) has the following form:

(21)
u θ = b 0 + b 1 ξ 1 + b 2 ξ 2 + b 3 ξ 1 2 - 1 + b 4 ξ 1 ξ 2 + b 5 ξ 2 2 - 1
            + + b 6 ξ 1 3 - 3 ξ 1 + b 7 ξ 1 2 ξ 2 - ξ 2 + b 8 ξ 1 ξ 2 2 - ξ 1 + b 9 ξ 2 3 - 3 ξ 2 ,

and in the new reduced model, the terms b4ξ1ξ2, b7ξ12ξ2-ξ2 and b8ξ1ξ22-ξ1 are eliminated while other terms are retained. It is noted that the variance of a non-Gaussian system parameter is derived from the following equation:

(22)
V A R u = k = 1 N b k 2 Ψ i ξ 0 θ , ξ 1 θ , , ξ j θ 2 .

Since the variance of terms retained are larger than those of the eliminated ones when they are of the same order [19], e.g. the variance of ξ13-3ξ1 and ξ12ξ2-ξ2 are equal to 6 and 2, respectively, when the total number of the random variables q increases, the summation of the variance of terms retained will be much larger than that of the eliminated ones. Therefore, the accuracy in the variance representation in Eq. (22) may still be maintained.

The reduced PCE approach on a non-Gaussian random process is implemented as follows: After the decomposition of a non-Gaussian random process u(x,t,θ) along the dimension x and t using the finite element method and KLE respectively, each uncorrelated random variable ξi(θ), which is assumed as “uncoupled” in this reduced PCE model, can be expanded as shown in Eq. (16) with a one-dimensional Polynomial Chaos. The dimension of PCE, Kr, required in a reduced PCE of order p for the stochastic response becomes:

(23)
K r = k s × p .

It is noted from Eqs. (20) and (23) that the dimension of PCE under this “uncoupled” assumption is significantly reduced compared with that in the “full” PCE. This allows the proposed algorithm to include a large number of K-L components in the PCE of the excitation force or the system parameters, with a trade off in the accuracy.

5.2. Modeling of the non-Gaussian system

5.2.1. The system matrices

According to the system equation of motion as described in Section 2, the KLE is applied to the elastic modulus E(x,θ) according to Eq. (11) as:

(24)
E x , θ = E - x + i = 1 k E ξ i θ λ i ϕ i x = i = 0 k E ξ i θ λ i ϕ i x ,

where kE is the number of K-L components. {ξi(θ)} are the normalized non-Gaussian random variables. Then:

(25)
K ~ e θ = i = 1 k E ξ i θ λ i ϕ i x I B e T B e d l = i = 1 k E ξ i θ K i e .

The stochastic system stiffness matrix K(θ) can then be expressed as:

(26)
K θ = K d + i = 1 k E ξ i θ K i ,

where Ki can be assembled from Kie. Let K0=Kd and ξ0(θ)= 1.0, we have:

(27)
K θ = i = 0 k E ξ i θ K i .

According to the reduced PCE for the non-Gaussian random process introduced in Section 5.1, the system matrices can be summarized as follows:

(28)
K θ = i = 0 k E ξ i θ K i = i = 0 k E k = 0 p Ψ i , k θ c 1 , k i K i = j = 0 K E Ψ j θ K j ' .

The stochastic system mass matrix can similarly be expressed as:

(29)
M θ = i = 0 k ρ ξ i θ M i = i = 0 k ρ k = 0 p Ψ i , k θ c 2 , k i M i = j = 0 K ρ Ψ j θ M j ' .

Since the Rayleigh damping matrix is a linear combination of the system mass and stiffness matrices according to Eq. (6), the damping matrix can also be written as:

(30)
C θ = i = 0 k c ξ i θ C i = i = 0 k c k = 0 p Ψ i , k θ c 3 , k i C i = j = 0 K c Ψ j θ C j ' ,

where Kj'=c1,k(i)Ki, Mj'=c2,k(i)Mi, Cj'=c3,k(i)Ci. c1,k(i), c2,k(i) and c3,k(i) are deterministic coefficients corresponding to the ith K-L component and kth order of the Polynomial Chaos for the elastic modulus, mass density and damping, respectively. Variable p is the order of Polynomial Chaos adopted, kE, kρ and kc are the number of the K-L components adopted for the elastic modulus, mass density and damping, respectively. KE, Kρ and Kc are the number of terms of the Polynomial Chaos in the reduced PCE for the elastic modulus, mass density and damping, respectively. According to the assumption of the reduced PCE, the relationship between the number of K-L components and the number of Polynomials for the PCE of the three system parameters are:

(31)
K E = k E × p ,       K ρ = k ρ × p ,       K c = k c × p .

5.2.2. The excitation forces

The stochastic excitation forces can be expanded similarly according to Eq. (28) as:

(32)
F t , θ = i = 0 k F ξ i θ F i t = i = 0 k F k = 0 p Ψ i , k θ c 4 , k i F i t = j = 0 K F Ψ j θ F j ' t ,

where Fj'(t)=c4,k(i)F(i)(t) and c4,k(i) is the deterministic coefficient corresponding to Ψi,k(θ); kF, p and KF are the number of K-L components, the order of Polynomial Chaos and the dimension of the reduced PCE for the excitation force vector, respectively.

5.2.3. The responses

The nodal displacement vector of the system with non-Gaussian property can be represented by the reduced PCE after truncation as:

(33)
R ( t , θ ) = i = 0 k s ξ i ( θ ) y ( i ) ( t ) = i = 0 k s k = 0 p Ψ i , k ( θ ) c 5 , k ( i ) y ( i ) ( t ) = j = 0 K r Ψ j θ y j ' t ,

where Kr is the dimension of the reduced PCE related to the order of the Polynomial Chaos p and the number of K-L components ks for the stochastic displacement vector with ks=kE+kρ+kF and yj'(t)=c5,k(i)y(i)(t).

Taking the first and second derivatives of Eq. (33) with respect to t, the stochastic velocity and acceleration vectors can be obtained respectively as:

(34)
R ˙ t , θ = j = 0 K r Ψ j θ y ˙ j ' t ,
(35)
R ¨ t , θ = j = 0 K r Ψ j θ y ¨ j ' t ,

where y˙j'(t)=c5,k(i)y˙(i)(t) and y¨j'(t)=c5,k(i)y¨(i)(t).

5.2.4. The dynamic equation of equilibrium

Substituting Eqs. (28) to (30) and (32) to (35) into Eq. (2), and taking inner product of both sides of Eq. (2) with respect to Ψk(θ) and employing the orthogonal property of Polynomial Chaos in Eq. (14), we have:

(36)
j = 0 K r i = 0 K ρ Ψ i θ Ψ j θ Ψ k θ M i ' y ¨ j ' t + j = 0 K r i = 0 K c Ψ i θ Ψ j θ Ψ k θ C i ' y ˙ j ' t
          + j = 0 K r i = 0 K E Ψ i θ Ψ j θ Ψ k θ K i ' y j ' t = H b Ψ k 2 F k ' t .

Let:

M k , j = i = 0 K ρ Ψ i θ Ψ j θ Ψ k θ Ψ k 2 M i ' ,       C k , j = i = 0 K c Ψ i θ Ψ j θ Ψ k θ Ψ k 2 C i ' ,
K k , j = i = 0 K E Ψ i θ Ψ j θ Ψ k θ Ψ k 2 K i ' .

Eq. (36) can be written in matrix form as:

(37)
M 0,0 M 1,0 M K r , 0 M 0,1 M 1,1 M K r , 1 M 0 , K r M 1 , K r M K r , K r y ¨ 0 ' t y ¨ 1 ' t y ¨ K r ' t + C 0,0 C 1,0 C K r , 0 C 0,1 C 1,1 C K r , 1 C 0 , K r C 1 , K r C K r , K r y ˙ 0 ' t y ˙ 1 ' t y ˙ K r ' t
              + K 0,0 K 1,0 K K r , 0 K 0,1 K 1,1 K K r , 1 K 0 , K r K 1 , K r K K r , K r y 0 ' t y 1 ' t y K r ' t = H b 0 0 0 H b 0 0 0 H b F 0 ' t F 1 ' t F K r ' t .

The values of inner product of polynomial chaos are constants and they can be obtained analytically [19].

5.3. Response statistics

Eq. (37) can be solved by employing the Newmark-β method to obtain the deterministic coefficients in the reduced PCE for the responses, and the first two response statistics of the stochastic nodal displacement vector can be evaluated as:

(38)
M E A N R t = y 0 ' t ,       V A R R t = j = 1 K r y j ' ( t ) 2 Ψ j 2 ,

where the subscript “R” denotes the random nodal displacement vector of the beam structure. The stochastic displacement of the beam structure at position x and time t can be derived according to Eqs. (7) and (38) as:

(39)
w ( x , t , θ ) = H ( x ) j = 0 K r Ψ j θ y j ' ( t ) = j = 0 K r Ψ j θ H x y j ' t .

The mean and variance of the stochastic displacement at position x and time t can be respectively obtained as:

(40)
M E A N w x , t = H x y 0 ' t ,       V A R w x , t = j = 1 K r H x y j ' ( t ) 2 Ψ j 2 ,

where the subscript “w” denotes the random displacement of the beam. The mean value and variance of the stochastic velocity and acceleration at position x and time t can also be obtained by substituting the first and second derivatives y˙j't and y¨j't obtained from Eq. (37) into Eq. (40). Samples of the stochastic displacement of the beam structure at position x and time t can be generated from Eq. (39) with any available sampling techniques, e.g. the Latin Hypercube sampling [24], to simulate the Gaussian random variables ξi(θ) in the Polynomial Chaos Ψj(θ). The probabilistic density function of the stochastic displacement of the beam structure at position x and time t can then be obtained from Eq. (39) after the coefficients yj't are calculated.

6. Numerical simulations

The following properties of the beam model are adopted in the numerical simulation: length of beam L= 40 m; cross-sectional area A= 4.8 m2; second moment of inertia of cross-section I= 2.5498 m4; damping ratio ζi= 0.02 for all modes; elastic modulus E and the mass density ρ are assumed as random variables with mean value 5×1010 N/m2 and 2.5×103 kg/m3respectively and with a coefficient of variation denoted as COVρ and COVE, respectively. The mean value of the first five natural frequencies of the beam are 3.9, 15.6, 35.1, 62.5 and 97.6 Hz. Rayleigh damping is assumed.

The random time-varying moving forces are generated with the following mean values:

(41)
F 1 = 20000 1 + 0.1 s i n 10 π t + 0.05 s i n 40 π t ,
(42)
F 2 = 20000 1 - 0.1 s i n 10 π t + 0.05 s i n 40 π t ,

and with different Coefficients of Variation (COV) at each time instant. The two loads are moving at a specific speed four meters apart as a group. It should be noted that although only two random forces are considered in the simulation, more forces can be included according to Eq. (1).

The relative error between the statistics of the calculated and reference responses are calculated as:

(43)
R E = R c a l c u l a t e d - R r e f e r e n c e 2 R r e f e r e n c e 2 × 100   % .

Since the system parameters are always positive in engineering practice, both the system parameters and the excitation forces are assumed as lognormal random processes in this study. To describe the uncertainty in the system parameters for the beam structure, covariance kernels for the stochastic system parameters have to be defined. The finite element model of beam is assumed to have eight Euler-Bernoulli beam elements of 5 m long each. In this simulation, the following exponential auto-covariance function is adopted as the covariance function for the stochastic system parameters:

(44)
C x 1 , x 2 = σ 2 e x p - x 1 - x 2 a ,

where σ is the standard deviation of system parameter E or ρ. Expression x1-x2 denotes the positive distance between two points in a spatial domain of interest. a is the correlation length. The positive distance between two points in a spatial domain of interest is set equal to 0.5 m unless it is stated otherwise and a is set as unity for both the elastic modulus E and mass density ρ in the following study. It should be noted that the covariance function for the two system parameters can be totally different with different values of x1-x2 and a. The sampling rate for all the simulations is 200 Hz. Results from the proposed algorithm is compared with those from the Monte Carlo Simulation on 10000 samples for verification.

6.1. Coefficients in the one-dimensional Polynomial Chaos expansion

According to the proposed method, the one-dimensional PCE for each uncoupled non-Gaussian random variable should be performed. A numerical study is conducted to compare the efficiency and accuracy of different existing techniques, such as, the Latin Hypercube sampling and the Fitting Regression model proposed by Choi et al. [23]. Considering the two functions Y1 and Y2 including the normal random variable ξ with a zero mean value and unit standard deviation as:

(45)
Y 1 = 2 + ξ + ξ 2 - 1 + 4 ξ 3 - 3 ξ ,
(46)
Y 2 = e ξ .

The coefficients in the one-dimensional PCE of the two functions are calculated with different existing techniques and with different order of PC up to four. The relative error and the time required in the calculation of coefficients of PCE using MCS method [21] with 10000 samples, the Latin Hypercube Sampling (LHS) method [22] with 300 samples, and Fitting Regression Model (FRM) method [23] using 300 samples respectively are shown in Table 1. The programme runs on a computer with Intel Core 2, Quad CPU 2.40 Hz and 8 GB RAM. Results show that the Fitting Regression Model method has both high efficiency and accuracy comparing with the other two methods especially for random variables which can be approximated as Polynomial functions of Gaussian random variables. Therefore, this method will be adopted in the following study to calculate the coefficients in the one-dimensional PCE in the reduced PCE model.

Table 1. Relative Error (%) in coefficients of one dimensional PCE calculated from three methods

Y 1
Y 2
Methods
MCS
LHS
FRM
MCS
LHS
FRM
b 0
0.07
0.14
0
0.01
0.06
0.03
b 1
0.04
0.98
0
0
0.29
0.29
b 2
0.09
1.84
0
0.05
0.61
0.24
b 3
0.23
2.75
0
0.08
0.80
0.67
b 4
0.02
2.81
0
0.08
0.80
0.25
Time (s)
35.34
0.17
0.14
33.07
0.14
0.11

6.2. Effect of the level of randomness in system parameters

Different levels of randomness in the system parameters, i.e. different COVρ and COVE are investigated. Different orders of Polynomial Chaos are adopted to represent the lognormal random processes of system parameters. To highlight the effect of the level of randomness in system parameters on the accuracy of the non-Gaussian model, the two excitation forces are assumed as deterministic (COVF= 0%) in this Section as defined in Eqs. (41) and (42).

The mean value and variance of the mid-span displacement of the beam structure calculated with the proposed method under different order of PC and the MCS are compared in Figs. 2 and 3, respectively. The detail informations from 0.5 s to 0.7 s are also demonstrated to achieve a better understanding of the results.

The relative errors of the results compared with those from the MCS according to Eq. (43) are listed in Table 2.

Results show that the calculated mean values for all the cases are very accurate and the relative error increases slightly with the increase of level of randomness in system parameters. The order of PC adopted has little effect on the accuracy of calculated mean value in general. For the variance of the mid-span displacement, the use of a low order of PC such as an order of 2, can achieve a high accuracy as noted in Table 2 when the COVE and COVρ are small and below 0.2. The results from large COV values can, however, be improved with the use of 3rd and 4th order PC, and the result from a 4th-order representation is better that that from a 3rd-order representation. Large error exists in the calculated variance when only a 2nd-order PCE is adopted for the cases with COVE and COVρ larger than 0.2.

Table 2. Relative error (%) in mid-span displacement statistics due to different level of PC representation when COVF= 0 %

C O V E = C O V ρ
C O V F =   0 %
Order of PC
5 %
10 %
20 %
30 %
40 %
50 %
Mean value
2
0.01
0.02
0.03
0.06
0.31
0.84
3
0.01
0.02
0.06
0.20
0.56
1.34
4
0.01
0.02
0.06
0.22
0.68
1.85
Variance
2
0.63
1.51
2.58
8.31
14.85
21.66
3
0.06
0.14
1.51
3.27
4.77
5.00
4
0.08
0.11
1.36
2.88
3.01
8.11

Fig. 2. Comparison of the mean mid-span displacement with different order of PC (COVF= 0 %, — MCS, ----- Order 2, –·– Order 3, – – Order 4)

 Comparison of the mean mid-span displacement with different order of PC (COVF= 0 %, — MCS, ----- Order 2, –·– Order 3, – – Order 4)

Fig. 3. Comparison of variance of mid-span displacement with different order of PC (COVF= 0 %, — MCS, ----- Order 2, –·– Order 3, – – Order 4)

 Comparison of variance of mid-span displacement with different order of PC (COVF= 0 %, — MCS, ----- Order 2, –·– Order 3, – – Order 4)

Another observation noted in Table 2 is that for a very large COVE and COVρ, e.g. larger than 0.4, the mean and variance from a 3rd-order representation are better that that from a 4th-order representation. This phenomenon can be explained as follows: the error in the coefficient calculation from the one-dimensional polynomial chaos expansion will be amplified by the K-L components of the random processes, and it will propagate into the reduced PCE model causing larger errors. It may be concluded that when more coefficients in the one-dimensional PCE are included in the calculation with higher order Polynomial Chaos, the computation error will be amplified and accumulated affecting the accuracy of the calculated response statistics.

Via the comparison of response statistics with the MCS method in Table 2, the accuracy of the proposed method is proved. The proposed method also has good efficiency, for example, in the case when COVE=COVρ= 20 %, 10 min is required for 10000 runs of MCS to obtain the response statistics and only 15 s is needed for the proposed method.

6.3. Effect of the distance of interest

It is noted that a refined distance between two points in a spatial domain of interest x1-x2 will increase the number of K-L components in the PCE of a non-Gaussian system parameter. According to Section 5.1, to increase the total number of K-L components in the reduced PCE can improve the accuracy of variance representation of the non-Gaussian parameter. Although the number of the KLE increases in this reduced PCE, yet the total number of Polynomial Chaos adopted under the assumption of “uncoupled” {ξi(θ)} is much smaller than that in a “full” PCE.

Fig. 4. Comparison of variance of mid-span displacement with different distance of interest (Order of PC = 3, COVF= 0 %, — MCS, -*- 0.5 m, ----- 1 m, –·– 2.5 m, – – 5 m)

 Comparison of variance of mid-span displacement with different distance of interest  (Order of PC = 3, COVF= 0 %, — MCS, -*- 0.5 m, ----- 1 m, –·– 2.5 m, – – 5 m)

A numerical study is carried out on the moving load on beam system with COVF= 0 % and the order of PCE equals to three in this sub-section to show the effect of the distance of interest on the accuracy of the predicted variance of the response. A comparison between the predicted variance of response using MCS and the reduced PCE method with different distance of interest are shown in Fig. 4, and the corresponding relative errors are shown in Table 3. The detail informations from 0.5 s to 0.6 s are also demonstrated in Fig. 4 to achieve a better understanding of the results.

Results show that when the level of randomness in the system parameters is not large, e.g. COV 0.2, with a relative large distance of interest x1-x2= 5.0 m, that is 1/8 of the total length of the beam, the relative error in the predicted variances of the response is less than 4 % indicating a good accuracy. When the level of randomness increase, with x1-x2= 5.0 m, larger error exists in the predicted variance of response. By decreasing the distance of interest, the relative error in the predicted variance of the response can be significantly reduced. For example, when the level of randomness in the system parameters has a COV= 0.5, and the distance of interest is decreased from 5.0 m to 1.0 m, the relative error in the predicted variance of response reduced from 57.3 % to 3.68 % according to results in Table 3. The relative error in the predicted variance of response decreases with smaller distance between two points of interest. It is noted from Table 3 that such reduction can be achieved favourably when x1-x2= 1.0 m. Due to the “uncoupled” assumption made in the proposed reduced PCE model, the relative error in the predicted variance of response is not linearly decreasing with the distance of interest.

Table 3. Relative error (%) in variance of mid-span displacement due to different distance of interest when COVF= 0 % and Order of PC = 3

C O V E = C O V ρ
C O V F =   0 %
Distance of interest (m)
5 %
10 %
20 %
30 %
40 %
50 %
Variance
5
0.23
0.64
1.94
3.09
11.9
57.3
2.5
0.32
1.03
3.79
7.38
10.1
10.7
1
0.15
0.33
1.15
1.97
2.04
3.68
0.5
0.06
0.14
1.51
3.27
4.77
5.00

6.4. Effect of the level of randomness in excitation

Investigations on the effect of the level of randomness in excitation is similarly conducted with COVE=COVρ= 20 % in this section. Different levels of randomness of the lognormal distributed excitation forces are adopted with different order in the Polynomial Chaos.

The mean and variance of the mid-span displacement of the beam structure obtained using the proposed method with different order of PCE and from the MCS are compared in Figs. 5 and 6, respectively. The detail informations from 0.5 s to 0.7 s are also demonstrated to achieve a better understanding of the results.The relative errors according to Eq. (43) from both methods are listed in Table 4.

Table 4. Relative error (%) in mid-span displacement statistics due to different level of PC representation when COVE=COVρ= 20 %

C O V F
C O V E = C O V ρ =   20 %
Order of PC
10 %
20 %
30 %
40 %
50 %
Mean value
2
1.81
1.90
1.99
2.11
2.24
3
1.79
1.83
1.84
1.87
1.91
4
1.79
1.83
1.83
1.80
1.88
Variance
2
4.23
11.2
14.0
17.8
23.1
3
3.35
4.44
3.93
4.83
5.21
4
3.59
6.90
3.82
12.5
16.2

Results show that the order of Polynomial Chaos has little effect on the accuracy of the mean mid-span displacement. When the COVF are equal or larger than 0.2, a relatively large error exists in the calculated variance from a 2nd-order representation. When the COVF becomes large, e.g. at 0.4 or 0.5, a 3rd-order representation has better performance than a 4th-order representation. The reason for this phenomenon has been explained in the last paragraph in Section 6.2.

Fig. 5. Comparison of mean mid-span displacement with different level of randomness in excitation and different order of PC used (COVE=COVρ= 20 %, — MCS, ----- Order 2, –·– Order 3, – – Order 4)

 Comparison of mean mid-span displacement with different level of randomness in excitation and different order of PC used (COVE=COVρ= 20 %, — MCS, ----- Order 2, –·– Order 3, – – Order 4)

Fig. 6. Comparison of variance of mid-span displacement with different level of randomness in excitation and different order of PC used (COVE=COVρ= 20 %, — MCS, ----- Order 2, –·– Order 3, – – Order 4)

 Comparison of variance of mid-span displacement with different level of randomness in excitation and different order of PC used (COVE=COVρ= 20 %, — MCS, ----- Order 2, –·– Order 3, – – Order 4)

It may be concluded that when a system with the coefficients of variation of both the system parameters and excitation forces smaller than 0.2, a 2nd-order representation of the Polynomial Chaos should be sufficient to obtain accurate response statistics of the system. In other cases with larger coefficients of variation, a 3rd-order representation of Polynomial Chaos is recommended. When the coefficient of variation of the excitation forces attains a large value of 0.5, the response statistics of the non-Gaussian model with a reduced PCE could still maintain a good accuracy.

7. Conclusions

A reduced PCE model is proposed in this paper to represent the non-Gaussian system parameters in a stochastic moving load on beam system. With the “uncoupled” assumption on the uncorrelated non-Gaussian random variables in the KLE of the non-Gaussian system parameters, the curse of dimensionality in a full PCE model can be alleviated with acceptable predictions on the response statistics of the beam. Numerical simulations with the proposed method and from the MCS show good agreement in the result for cases with different level of uncertainties in the excitation and system parameters and different order of Polynomial Chaos used. The following conclusions can be drawn:

1) The SSFEM with the proposed reduced PCE model can achieve good accuracy in the mean and variance predictions even when the COVρ and COVE are as large as 0.5. When COVE and COVρ are smaller or equal to 0.2, the proposed approach has good accuracy even with large uncertainties in the excitation forces. The relative errors in the first two response statistics of the beam responses increase with an increase in the level of randomness in both the system parameters and excitation forces.

2) When the level of randomness in the system parameters is COV 0.2, with a distance of interest x1-x2 equal to 1/8 total length of the beam, the relative error in the predicted variances of the response are less than 4 %; when the level of randomness increases to 50 %, an error in the predicted variances of the response up to 57.3 % will occur. The relative error in the predicted variances of the response can be significantly reduced to below 5 % when the distance of interest is decreased.

3) For a system with the coefficients of variation of both the system parameters and excitation forces smaller than 0.2, a 2nd-order representation of the Polynomial Chaos should be sufficient to obtain accurate response statistics of the beam structure. In other cases of larger coefficients of variation, a 3rd-order representation of Polynomial Chaos is recommended.

Acknowledgements

The research work was supported by a research grant from the National Natural Science Foundation of China (No. 11402052), the Natural Science Foundation of Jiangsu province (No. BK20140616), the Specialized Research Fund for the Doctoral Program of Higher Education of China (20130092120039) and the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD).

References

  1. Marchesiello S., Fasana A., Garibaldi L., Piombo B. A. D. Dynamic of multi-span continuous straight bridges subject to multi-degrees of freedom moving vehicle excitation. Journal of Sound Vibration, Vol. 224, Issue 3, 1999, p. 541-561. [Search CrossRef]
  2. Zhu X. Q., Law S. S. Dynamic behavior of orthotropic rectangular plates under moving loads. Journal of Engineering Mechanics, Vol. 129, Issue 1, 2003, p. 79-87. [Search CrossRef]
  3. Henchi K., Fafard M., Talbot M., Dhatt G. An efficient algorithm for dynamic analysis of bridges under moving vehicles using a coupled modal and physical components approach. Journal of Sound Vibration, Vol. 212, Issue 4, 1998, p. 663-683. [Search CrossRef]
  4. Lee S. Y., Yhim S. S. Dynamic behavior of long-span box girder bridges subjected to moving loads: Numerical analysis and experimental verification. International Journal of Solids and Structures, Vol. 42, Issues 18-19, 2005, p. 5021-5035. [Search CrossRef]
  5. ISO 8606:1995(E). Mechanical Vibration – Road Surface Profiles – Reporting of Measured Data. 1995. [Search CrossRef]
  6. Sasidhar M. N. V., Talukdar S. Nonstationary response of bridge due to eccentrically moving vehicle at variable velocity. Advances in Structural Engineering, Vol. 6, Issue 4, 2003, p. 309-324. [Search CrossRef]
  7. Seetapan P., Chucheepsakul S. Dynamic responses of a two-span beam subjected to high speed 2DOF sprung vehicles. International Journal of Structural Stability and Dynamics, Vol. 6, Issue 3, 2006, p. 413-430. [Search CrossRef]
  8. Da Silva J. G. S. Dynamical performance of highway bridge decks with irregular pavement surface. Computers and Structures,Vol. 82, Issues 11-12, 2004, p. 871-881. [Search CrossRef]
  9. Fryba L., Nakagiri S., Yoshikawa N. Stochastic finite elements for beam on a random foundation with uncertain damping under a moving force. Journal of Sound Vibration, Vol.163, Issue 1, 2003, p. 31-45. [Search CrossRef]
  10. Muscolino G., Benfratello S., Sidoti A. Dynamics analysis of distributed parameter system subject to a moving oscillator with random mass, velocity and acceleration. Probabilistic Engineering Mechanics, Vol. 17, 2002, p. 63-72. [Search CrossRef]
  11. Chang T. P., Lin G. L., Chang E. Vibration analysis of a beam with an internal hinge subject to a random moving oscillator. International Journal of Solids and Structures,Vol. 43, 2006, p. 6398-6412. [Search CrossRef]
  12. Chang T. P., Liu M. F., O H. W. Vibration analysis of a uniform beam traversed by a moving vehicle with random mass and random velocity. Structural Engineering and Mechanics, Vol. 31, Issue 6, 2009, p. 737-749. [Search CrossRef]
  13. Wu S. Q., Law S. S. Dynamic analysis of bridge-vehicle system with uncertainties based on the finite element model. Probabilistic Engineering Mechanics, Vol. 25, 2010, p. 425-432. [Search CrossRef]
  14. Wu S. Q., Law S. S. Dynamic analysis of bridge with non-Gaussian uncertainties under a moving vehicle. Probabilistic Engineering Mechanics, Vol. 26, 2011, p. 281-293. [Search CrossRef]
  15. Liu N., Gao W., Song C., Zhang N., Pi Y. Interval dynamic response analysis of vehicle-bridge interaction system with uncertainty. Journal of Sound Vibration, Vol. 332, 2013, p. 3218-3231. [Search CrossRef]
  16. Stefanou G. The stochastic finite element method: past, present and future. Computer Methods in Applied Mechanics and Engineering, Vol. 198, 2009, p. 1031-1051. [Search CrossRef]
  17. Eiermann M., Ernst O. G., Ullmann E. Computational aspects of the stochastic finite element method. Computing and Visualization in Science, Vol. 10, Issue 1, 2007, p. 3-15. [Search CrossRef]
  18. Schenk C. A., Schuëller G. I. Uncertainty Assessment of Large Finite Element Systems. Springer-Verlag, Berlin Heidelberg, 2005. [Search CrossRef]
  19. Ghanem R., Spanos P. D. Stochastic Finite Elements: A Spectral Approach. Springer-Verlag, New York, 1991. [Search CrossRef]
  20. Ghanem R. Stochastic finite elements with multiple random non-Gaussian properties. Journal of Engineering Mechanics-ASCE, Vol. 125, Issue 1, 1999, p. 26-40. [Search CrossRef]
  21. Field R. V. Numerical methods to estimate the coefficients of the polynomial chaos expansion. 15th ACSE Engineering Mechanics Conference, Columbia University, New York, 2002. [Search CrossRef]
  22. Reagan M. T., Najm H. N., Ghanem R., Knio O. M. Uncertainty quantification in reacting-flow simulations through non-intrusive spectral projection. Combustion and Flame, Vol. 132, 2005, p. 545-555. [Search CrossRef]
  23. Choi S. K., Grandhi R. V., Canfield R. A., Pettit C. L. Polynomial chaos expansion with latin hypercube sampling for estimating response variability. AIAA Journal, Vol. 42, Issue 6, 2004, p. 1191-1198. [Search CrossRef]
  24. Florian A. An efficient sampling scheme: updated latin hypercube sampling. Probabilistic Engineering Mechanics, Vol. 7, 1992, p. 123-130. [Search CrossRef]