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 KL 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 nonGaussian random variables in the KarhunenLoéve Expansion (KLE) of the nonGaussian 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.
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 bridgevehicle 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 bridgevehicle interaction.
The bridgevehicle 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 firstorder 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 EulerBernoulli 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 EulerBernoulli beam with uncertain parameters such as random mass, stiffness, damping, etc. which is subjected to a moving oscillator/halfcar model with random velocity and acceleration. Wu and Law [13] performed dynamic analysis on the bridgevehicle interaction problem considering Gaussian uncertainties using KarhunenLoéve expansion. However, the system parameters are always positive in engineering practice, and it is more appropriate to model them with nonGaussian assumption. Wu and Law [14] later considered the nonGaussian system parameters in the bridgevehicle 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 bridgevehicle model as interval variables and perform an interval dynamic analysis on the bridgevehicle 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 KarhunenLoéve (KL) 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 KL 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 KL 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 nonGaussian 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 twodimensional simply supported EulerBernoulli beam with multiple loads moving on top is shown in Fig. 1. The mass per unit length $\rho $, elastic modulus $E$ and damping $c$ of the beam structure are assumed as nonGaussian random processes, with mean value $\stackrel{}{\rho}$, $\stackrel{}{E}$, $\stackrel{}{c}$ and standard deviation ${\sigma}_{\rho}$, ${\sigma}_{E}$, ${\sigma}_{c}$, respectively. The random part of the these parameters, denoted as $\stackrel{~}{\rho}$, $\stackrel{~}{E}$, $\stackrel{~}{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:
$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}=\sum _{j=1}^{{N}_{F}}{F}_{j}\left(t,\theta \right)\delta \left(x{v}_{j}t\right),$
where $A$ is the crosssectional area, $I$ is the moment of inertia of the beam and $\theta $ denotes the random dimension. ${N}_{F}$ 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:
where $\mathbf{M}\left(\theta \right)$, $\mathbf{C}\left(\theta \right)$, $\mathbf{K}\left(\theta \right)$ are the stochastic mass, damping and stiffness matrices of the beam model with $\mathbf{M}\left(\theta \right)={\mathbf{M}}_{d}+\stackrel{~}{\mathbf{M}}\left(\theta \right)$, $\mathbf{C}\left(\theta \right)={\mathbf{C}}_{d}+\stackrel{~}{\mathbf{C}}\left(\theta \right)$, $\mathbf{K}\left(\theta \right)={\mathbf{K}}_{d}+\stackrel{~}{\mathbf{K}}\left(\theta \right)$, respectively. ${\mathbf{M}}_{d}$, $\stackrel{~}{\mathbf{M}}\left(\theta \right)$, ${\mathbf{C}}_{d}$, $\stackrel{~}{\mathbf{C}}\left(\theta \right)$, ${\mathbf{K}}_{d}$, $\stackrel{~}{\mathbf{K}}\left(\theta \right)$ 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:
${\mathbf{K}}_{d}^{e}={\int}_{l}{\mathbf{B}}^{eT}\stackrel{}{E}I{\mathbf{B}}^{e}dl,{\stackrel{~}{\mathbf{K}}}^{e}\left(\theta \right)={\int}_{l}{\mathbf{B}}^{eT}\stackrel{~}{E}\left(\theta \right)I{\mathbf{B}}^{e}dl,$
where ${\mathbf{H}}^{e}$ and ${\mathbf{B}}^{e}$ are respectively the shape function matrix and straindisplacement matrix for each element. $l$ is the length of beam element. $\mathbf{F}(t,\theta )$ is the stochastic force vector applied on the beam structure with $\mathbf{F}={\left\{\begin{array}{llll}{F}_{1}& {F}_{2}& \cdots & {F}_{{N}_{F}}\end{array}\right\}}^{T}$. $\mathbf{R}(t,\theta )$, $\dot{\mathbf{R}}\left(t,\theta \right)$ and $\ddot{\mathbf{R}}\left(t,\theta \right)$ are the stochastic nodal displacement, velocity and acceleration vectors of the structure respectively with $\mathbf{R}={\left\{{\mathbf{R}}_{1}\mathrm{}\mathrm{}{\mathbf{R}}_{2}\mathrm{}\mathrm{}\cdots \mathrm{}\mathrm{}{\mathbf{R}}_{n}\right\}}^{T}$. ${\mathbf{H}}_{b}$ is a $n$×${N}_{F}$ location matrix of the moving forces as shown in Eq. (4), where $n$ is the number of degreeoffreedom of the beam structure after considering the boundary condition. We have:
where the shape function vector ${\mathbf{H}}_{i}$ can be written in the global coordinate as:
and ${x}_{j}\left(t\right)$ is the location of the $j$th force on the $i$th element at time $t$ with $(i1)l\le {x}_{j}\left(t\right)<il$. The two coefficients ${c}_{M}$ and ${c}_{K}$ of the Rayleigh damping:
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:
where $\mathbf{H}\left(x\right)=\left\{0\mathrm{}\mathrm{}\cdots \mathrm{}\mathrm{}\mathbf{H}{\left(x\right)}_{i}^{T}\mathrm{}\mathrm{}0\mathrm{}\mathrm{}\cdots \mathrm{}\mathrm{}0\right\}$ and $\mathbf{H}{\left(x\right)}_{i}^{T}$ is the shape function of the beam structure. $\mathbf{H}\left(x\right)$ is a 1×$n$ vector with zero entries except at the $i$th beam element on which the force is located.
Fig. 1The beammoving load system
3. KarhunenLoéve expansion
3.1. Theory [18]
Let $u(x,\theta )$ be a secondorder random process and $\stackrel{}{u}\left(x\right)$ denote the expected value of $u(x,\theta )$. The covariance function $C({x}_{1},{x}_{2})$ of $u(x,\theta )$ is bounded, symmetric and positive definite and has the following spectral decomposition:
where ${\lambda}_{n}$ and ${\phi}_{n}\left(x\right)$ are the eigenvalues and the eigenvectors of the covariance function, respectively. They can be proved to be solution of the following integral equation:
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:
where ${\delta}_{ij}$ is the Kronecker delta. The random process $u(x,\theta )$ can then be written as:
where ${\xi}_{n}\left(\theta \right)$ is a group of uncorrelated random variables. As a special case, when $u(x,\theta )$ is a Gaussian random process, ${\xi}_{n}\left(\theta \right)$ is a set of standard Gaussian random variables with the following orthogonal properties:
where $E(\bullet )$ denotes the expectation value.
4. Polynomial Chaos expansion
4.1. Theory [19]
For an element $u\left(\theta \right)$ from a set in the space spanned by the orthogonal basis, it has been proved that $u\left(\theta \right)$ can be written in the following form as:
where ${\mathrm{\Psi}}_{i}$ are the Hermite polynomials of ${\xi}_{1}\left(\theta \right)$, ${\xi}_{2}\left(\theta \right)$,…, ${\xi}_{j}\left(\theta \right)$, and ${\mathrm{\Psi}}_{0}=$ 1.0, ${\xi}_{0}\left(\theta \right)=$ 1.0. ${b}_{k}$ are the generalized Fourier coefficients with $k=$ (0, 1,…, $N$) and $N=1+\left(p+q\right)!/\left(p!q!\right)$. The series in Eq. (13) can be refined along the random dimension $\theta $ 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 ${\mathrm{\Psi}}_{k}$ are orthogonal satisfying the relationship:
where $\u27e8\bullet \u27e9$ denotes the inner product and the value of $\u27e8{\mathrm{\Psi}}_{i}^{2}\u27e9$ can be calculated analytically [19].
The generalized Fourier coefficients ${b}_{k}$ in Eq. (13) can be evaluated from:
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. Nonintrusive 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 onedimensional PCE in the proposed reduced PCE model.
For a onedimensional problem, i.e. $q$ in Eq. (13) equals to unity, consider the PCE truncated at the $p$th level:
$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}={b}_{0}+{b}_{1}{f}_{1}\left(\xi \right)+{b}_{2}{f}_{2}\left(\xi \right)+\cdots +{b}_{p}{f}_{p}\left(\xi \right)+\epsilon ,$
where ${b}_{i}$, $i=$ 0, 1, 2,…, $p$, are the regression coefficients, and $\epsilon $ is the error of the model equation. Functions ${f}_{1}\left(\xi \right)=\xi $, ${f}_{2}\left(\xi \right)={\xi}^{2}1$,…, ${f}_{p}\left(\xi \right)$ are the $p$th onedimensional Hermite polynomials of $\xi \left(\theta \right)$.
It is noted that Eq. (16) can be regarded as a linear regression model. To evaluate the coefficients ${b}_{i}$, $n$ sample values, e.g. ${z}_{1}$~${z}_{n}$, should be drawn from the standard random variable $\xi \left(\theta \right)$. With these sample values, the matrix form can be formulated as:
where $\mathbf{Y}={\left\{\begin{array}{llll}{y}_{1}& {y}_{2}& \cdots & {y}_{n}\end{array}\right\}}^{T}$, $\widehat{\mathbf{B}}={\left\{\begin{array}{llll}{b}_{0}& {b}_{1}& \cdots & {b}_{q}\end{array}\right\}}^{T}$, $\mathbf{e}={\left\{\begin{array}{llll}{\epsilon}_{1}& {\epsilon}_{2}& \cdots & {\epsilon}_{n}\end{array}\right\}}^{T}$ and:
and ${y}_{1}$ to ${y}_{n}$ are the sample values drawn from $u\left(\theta \right)$. Vector $\mathbf{e}$ includes the error terms. The sample values of ${f}_{i}\left({z}_{k}\right)$ ($i=$ 1, 2,…, $p$) in Eq. (17) can be obtained from the Cumulative Density Function (CDF) of the random variables.
Using the leastsquares method, the regression coefficients can be calculated as:
The fitted model and the residuals are:
5. Modeling with nonGaussian uncertainties
5.1. The reduced PCE model
The KLE of a nonGaussian process results in a set of deterministic coefficients multiplied by the corresponding uncorrelated nonGaussian random variables $\left\{{\xi}_{i}\right(\theta \left)\right\}$. In the PCE of a nonGaussian process, these uncorrelated nonGaussian random variables are projected on the Polynomial Chaos basis. The dimension of PCE, denoted as ${K}_{R}$,_{}to represent the stochastic response of a beam structure is determined according to the following [20]:
where $p$ is the maximum order of the Polynomial Chaos in PCE; ${k}_{s}$ is the number of KL components required in the representation of the response which depends on the number of KL 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 ${k}_{s}={k}_{E}+{k}_{\rho}+{k}_{F}$. Variables ${k}_{E}$, ${k}_{\rho}$, ${k}_{F}$ are the number of KL components in the KLE of the Gaussian distributed elastic modulus, mass density and the excitation forces, respectively.
When a large number of KL components is required inthe KLE of the excitation force, an extremely large number of Polynomial Chaos will be included in the PCE of the nonGaussian 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 nonGaussian 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 ${\mathrm{\Psi}}_{i}\left({\xi}_{j}\left(\theta \right)\right)$ in Eq. (13) and eliminating the terms of ${\mathrm{\Psi}}_{i}\left({\xi}_{1}\left(\theta \right),{\xi}_{2}\left(\theta \right),\dots ,{\xi}_{j}\left(\theta \right)\right)$ in the new reduced model. For example, when $p=$ 2 and $q=$ 2, Eq. (13) has the following form:
${\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}+\cdots +b}_{6}\left({\xi}_{1}^{3}3{\xi}_{1}\right)+{b}_{7}\left({\xi}_{1}^{2}{\xi}_{2}{\xi}_{2}\right)+{b}_{8}\left({\xi}_{1}{\xi}_{2}^{2}{\xi}_{1}\right)+{b}_{9}\left({\xi}_{2}^{3}3{\xi}_{2}\right),$
and in the new reduced model, the terms ${b}_{4}{\xi}_{1}{\xi}_{2}$, ${b}_{7}\left({\xi}_{1}^{2}{\xi}_{2}{\xi}_{2}\right)$ and ${b}_{8}\left({\xi}_{1}{\xi}_{2}^{2}{\xi}_{1}\right)$ are eliminated while other terms are retained. It is noted that the variance of a nonGaussian system parameter is derived from the following equation:
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 ${\xi}_{1}^{3}3{\xi}_{1}$ and ${\xi}_{1}^{2}{\xi}_{2}{\xi}_{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 nonGaussian random process is implemented as follows: After the decomposition of a nonGaussian random process $u(x,t,\theta )$ along the dimension $x$ and $t$ using the finite element method and KLE respectively, each uncorrelated random variable ${\xi}_{i}\left(\theta \right)$, which is assumed as “uncoupled” in this reduced PCE model, can be expanded as shown in Eq. (16) with a onedimensional Polynomial Chaos. The dimension of PCE, ${K}_{r}$, required in a reduced PCE of order $p$ for the stochastic response becomes:
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 KL components in the PCE of the excitation force or the system parameters, with a trade off in the accuracy.
5.2. Modeling of the nonGaussian 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,\theta )$ according to Eq. (11) as:
where ${k}_{E}$ is the number of KL components. $\left\{{\xi}_{i}\right(\theta \left)\right\}$ are the normalized nonGaussian random variables. Then:
The stochastic system stiffness matrix $\mathbf{K}(\theta $) can then be expressed as:
where ${\mathbf{K}}_{i}$ can be assembled from ${\mathbf{K}}_{i}^{e}$. Let ${\mathbf{K}}_{0}={\mathbf{K}}_{d}$ and ${\xi}_{0}\left(\theta \right)=$ 1.0, we have:
According to the reduced PCE for the nonGaussian random process introduced in Section 5.1, the system matrices can be summarized as follows:
The stochastic system mass matrix can similarly be expressed as:
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:
where ${\mathbf{K}}_{j}^{\mathrm{\text{'}}}={c}_{1,k}^{\left(i\right)}{\mathbf{K}}_{i}$, ${\mathbf{M}}_{j}^{\mathrm{\text{'}}}={c}_{2,k}^{\left(i\right)}{\mathbf{M}}_{i}$, ${\mathbf{C}}_{j}^{\mathrm{\text{'}}}={c}_{3,k}^{\left(i\right)}{\mathbf{C}}_{i}$. ${c}_{1,k}^{\left(i\right)}$, ${c}_{2,k}^{\left(i\right)}$ and ${c}_{3,k}^{\left(i\right)}$ are deterministic coefficients corresponding to the $i$th KL component and $k$th order of the Polynomial Chaos for the elastic modulus, mass density and damping, respectively. Variable $p$ is the order of Polynomial Chaos adopted, ${k}_{E}$, ${k}_{\rho}$ and ${k}_{c}$ are the number of the KL components adopted for the elastic modulus, mass density and damping, respectively. ${K}_{E}$, ${K}_{\rho}$ and ${K}_{c}$ 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 KL components and the number of Polynomials for the PCE of the three system parameters are:
5.2.2. The excitation forces
The stochastic excitation forces can be expanded similarly according to Eq. (28) as:
where ${\mathbf{F}}_{j}^{\text{'}}\left(t\right)={c}_{4,k}^{\left(i\right)}{\mathbf{F}}^{\left(i\right)}\left(t\right)$ and ${c}_{4,k}^{\left(i\right)}$ is the deterministic coefficient corresponding to ${\mathrm{\Psi}}_{i,k}\left(\theta \right)$; ${k}_{F}$, $p$ and ${K}_{F}$ are the number of KL 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 nonGaussian property can be represented by the reduced PCE after truncation as:
where ${K}_{r}$ is the dimension of the reduced PCE related to the order of the Polynomial Chaos $p$ and the number of KL components ${k}_{s}$ for the stochastic displacement vector with ${k}_{s}={k}_{E}+{k}_{\rho}+{k}_{F}$ and ${\mathbf{y}}_{j}^{\mathbf{\text{'}}}\left(t\right)={c}_{5,k}^{\left(i\right)}{\mathbf{y}}^{\left(i\right)}\left(t\right)$.
Taking the first and second derivatives of Eq. (33) with respect to $t$, the stochastic velocity and acceleration vectors can be obtained respectively as:
where ${\dot{\mathbf{y}}}_{j}^{\mathbf{\text{'}}}\left(t\right)={c}_{5,k}^{\left(i\right)}{\dot{\mathbf{y}}}^{\left(i\right)}\left(t\right)$ and ${\ddot{\mathbf{y}}}_{j}^{\mathbf{\text{'}}}\left(t\right)={c}_{5,k}^{\left(i\right)}{\ddot{\mathbf{y}}}^{\left(i\right)}\left(t\right)$.
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 ${\mathrm{\Psi}}_{k}\left(\theta \right)$ and employing the orthogonal property of Polynomial Chaos in Eq. (14), we have:
$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}+\sum _{j=0}^{{K}_{r}}\sum _{i=0}^{{K}_{E}}\u27e8{\mathrm{\Psi}}_{i}\left(\theta \right){\mathrm{\Psi}}_{j}\left(\theta \right){\mathrm{\Psi}}_{k}\left(\theta \right)\u27e9{\mathbf{K}}_{i}^{\text{'}}{\mathbf{y}}_{j}^{\text{'}}\left(t\right)={\mathbf{H}}_{b}\u27e8{\mathrm{\Psi}}_{k}^{2}\u27e9{\mathbf{F}}_{k}^{\text{'}}\left(t\right).$
Let:
${\mathbf{K}}^{\left(k,j\right)}=\sum _{i=0}^{{K}_{\mathrm{E}}}\frac{\u27e8{\mathrm{\Psi}}_{i}\left(\theta \right){\mathrm{\Psi}}_{j}\left(\theta \right){\mathrm{\Psi}}_{k}\left(\theta \right)\u27e9}{\u27e8{\mathrm{\Psi}}_{k}^{2}\u27e9}{\mathbf{K}}_{i}^{\text{'}}.$
Eq. (36) can be written in matrix form as:
$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}+\left[\begin{array}{cccc}{\mathbf{K}}^{\left(\mathrm{0,0}\right)}& {\mathbf{K}}^{\left(\mathrm{1,0}\right)}& \cdots & {\mathbf{K}}^{\left({K}_{r},0\right)}\\ {\mathbf{K}}^{\left(\mathrm{0,1}\right)}& {\mathbf{K}}^{\left(\mathrm{1,1}\right)}& \cdots & {\mathbf{K}}^{\left({K}_{r},1\right)}\\ \vdots & \vdots & \ddots & \vdots \\ {\mathbf{K}}^{\left(0,{K}_{r}\right)}& {\mathbf{K}}^{\left(1,{K}_{r}\right)}& \cdots & {\mathbf{K}}^{\left({K}_{r},{K}_{r}\right)}\end{array}\right]\left\{\begin{array}{c}{\mathbf{y}}_{0}^{\mathrm{\text{'}}}\left(t\right)\\ {\mathbf{y}}_{1}^{\mathrm{\text{'}}}\left(t\right)\\ \vdots \\ {\mathbf{y}}_{{K}_{r}}^{\text{'}}\left(t\right)\end{array}\right\}=\left[\begin{array}{cccc}{\mathbf{H}}_{b}& 0& \cdots & 0\\ 0& {\mathbf{H}}_{b}& \cdots & 0\\ \vdots & \vdots & \ddots & \vdots \\ 0& 0& \cdots & {\mathbf{H}}_{b}\end{array}\right]\left\{\begin{array}{c}{\mathbf{F}}_{0}^{\mathrm{\text{'}}}\left(t\right)\\ {\mathbf{F}}_{1}^{\mathrm{\text{'}}}\left(t\right)\\ \vdots \\ {\mathbf{F}}_{{K}_{r}}^{\text{'}}\left(t\right)\end{array}\right\}.$
The values of inner product of polynomial chaos $\u27e8\bullet \u27e9$ are constants and they can be obtained analytically [19].
5.3. Response statistics
Eq. (37) can be solved by employing the Newmark$\beta $ 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:
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:
The mean and variance of the stochastic displacement at position $x$ and time $t$ can be respectively obtained as:
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 ${\dot{\mathbf{y}}}_{j}^{\mathrm{\text{'}}}\left(t\right)$ and ${\ddot{\mathbf{y}}}_{j}^{\mathrm{\text{'}}}\left(t\right)$ 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 ${\xi}_{i}\left(\theta \right)$ in the Polynomial Chaos ${\mathrm{\Psi}}_{j}\left(\theta \right)$. 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 ${\mathbf{y}}_{j}^{\mathrm{\text{'}}}\left(t\right)$ are calculated.
6. Numerical simulations
The following properties of the beam model are adopted in the numerical simulation: length of beam $L=$ 40 m; crosssectional area $A=$ 4.8 m^{2}; second moment of inertia of crosssection $I=$ 2.5498 m^{4}; damping ratio ${\zeta}_{i}=$ 0.02 for all modes; elastic modulus $E$ and the mass density $\rho $ are assumed as random variables with mean value 5×10^{10} N/m^{2} and 2.5×10^{3} kg/m^{3}^{}respectively and with a coefficient of variation denoted as ${\mathrm{C}\mathrm{O}\mathrm{V}}_{\rho}$ and ${\mathrm{C}\mathrm{O}\mathrm{V}}_{E}$, 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 timevarying moving forces are generated with the following mean values:
and with different Coefficients of Variation ($\mathrm{C}\mathrm{O}\mathrm{V}$) 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:
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 EulerBernoulli beam elements of 5 m long each. In this simulation, the following exponential autocovariance function is adopted as the covariance function for the stochastic system parameters:
where $\sigma $ is the standard deviation of system parameter $E$ or $\rho $. Expression $\left{x}_{1}{x}_{2}\right$ 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 $\rho $ 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 $\left{x}_{1}{x}_{2}\right$ 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 onedimensional Polynomial Chaos expansion
According to the proposed method, the onedimensional PCE for each uncoupled nonGaussian 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 ${Y}_{1}$ and ${Y}_{2}$ including the normal random variable $\xi $ with a zero mean value and unit standard deviation as:
The coefficients in the onedimensional 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 onedimensional PCE in the reduced PCE model.
Table 1Relative 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 $\mathrm{C}\mathrm{O}{\mathrm{V}}_{\rho}$ and $\mathrm{C}\mathrm{O}{\mathrm{V}}_{E}$ 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 nonGaussian model, the two excitation forces are assumed as deterministic ($\mathrm{C}\mathrm{O}{\mathrm{V}}_{F}=$ 0%) in this Section as defined in Eqs. (41) and (42).
The mean value and variance of the midspan 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 midspan 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 $\mathrm{C}\mathrm{O}{\mathrm{V}}_{E}$ and $\mathrm{C}\mathrm{O}{\mathrm{V}}_{\rho}$ 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 4thorder representation is better that that from a 3rdorder representation. Large error exists in the calculated variance when only a 2ndorder PCE is adopted for the cases with $\mathrm{C}\mathrm{O}{\mathrm{V}}_{E}$ and $\mathrm{C}\mathrm{O}{\mathrm{V}}_{\rho}$ larger than 0.2.
Table 2Relative error (%) in midspan displacement statistics due to different level of PC representation when COVF= 0 %
$\mathrm{C}\mathrm{O}{\mathrm{V}}_{E}={\mathrm{C}\mathrm{O}\mathrm{V}}_{\rho}$  
${\mathrm{C}\mathrm{O}\mathrm{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. 2Comparison of the mean midspan displacement with different order of PC (COVF= 0 %, — MCS,  Order 2, –·– Order 3, – – Order 4)
Fig. 3Comparison of variance of midspan 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 $\mathrm{C}\mathrm{O}{\mathrm{V}}_{E}$ and $\mathrm{C}\mathrm{O}{\mathrm{V}}_{\rho}$, e.g. larger than 0.4, the mean and variance from a 3rdorder representation are better that that from a 4thorder representation. This phenomenon can be explained as follows: the error in the coefficient calculation from the onedimensional polynomial chaos expansion will be amplified by the KL 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 onedimensional 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 $\mathrm{C}\mathrm{O}{\mathrm{V}}_{E}=\mathrm{C}\mathrm{O}{\mathrm{V}}_{\rho}=\mathrm{}$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 $\left{x}_{1}{x}_{2}\right$ will increase the number of KL components in the PCE of a nonGaussian system parameter. According to Section 5.1, to increase the total number of KL components in the reduced PCE can improve the accuracy of variance representation of the nonGaussian 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” $\left\{{\xi}_{i}\right(\theta \left)\right\}$ is much smaller than that in a “full” PCE.
Fig. 4Comparison of variance of midspan 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 ${\mathrm{C}\mathrm{O}\mathrm{V}}_{F}=$0 % and the order of PCE equals to three in this subsection 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. $\mathrm{C}\mathrm{O}\mathrm{V}\le $ 0.2, with a relative large distance of interest $\left{x}_{1}{x}_{2}\right=$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 $\left{x}_{1}{x}_{2}\right=$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 $\mathrm{C}\mathrm{O}\mathrm{V}=$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 $\left{x}_{1}{x}_{2}\right=$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 3Relative error (%) in variance of midspan displacement due to different distance of interest when COVF= 0 % and Order of PC = 3
${\mathrm{C}\mathrm{O}\mathrm{V}}_{E}={\mathrm{C}\mathrm{O}\mathrm{V}}_{\rho}$  
${\mathrm{C}\mathrm{O}\mathrm{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 ${\mathrm{C}\mathrm{O}\mathrm{V}}_{E}={\mathrm{C}\mathrm{O}\mathrm{V}}_{\rho}=$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 midspan 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 4Relative error (%) in midspan displacement statistics due to different level of PC representation when COVE=COVρ= 20 %
${\mathrm{C}\mathrm{O}\mathrm{V}}_{F}$  
${\mathrm{C}\mathrm{O}\mathrm{V}}_{E}={\mathrm{C}\mathrm{O}\mathrm{V}}_{\rho}=$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 midspan displacement. When the ${\mathrm{C}\mathrm{O}\mathrm{V}}_{F}$ are equal or larger than 0.2, a relatively large error exists in the calculated variance from a 2ndorder representation. When the ${\mathrm{C}\mathrm{O}\mathrm{V}}_{F}$ becomes large, e.g. at 0.4 or 0.5, a 3rdorder representation has better performance than a 4thorder representation. The reason for this phenomenon has been explained in the last paragraph in Section 6.2.
Fig. 5Comparison of mean midspan 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. 6Comparison of variance of midspan 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 2ndorder 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 3rdorder 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 nonGaussian 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 nonGaussian system parameters in a stochastic moving load on beam system. With the “uncoupled” assumption on the uncorrelated nonGaussian random variables in the KLE of the nonGaussian 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 ${\mathrm{C}\mathrm{O}\mathrm{V}}_{\rho}$ and ${\mathrm{C}\mathrm{O}\mathrm{V}}_{E}$ are as large as 0.5. When ${\mathrm{C}\mathrm{O}\mathrm{V}}_{E}$ and ${\mathrm{C}\mathrm{O}\mathrm{V}}_{\rho}$ 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 $\mathrm{C}\mathrm{O}\mathrm{V}\le $ 0.2, with a distance of interest $\left{x}_{1}{x}_{2}\right$ 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 2ndorder 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 3rdorder representation of Polynomial Chaos is recommended.
References

Marchesiello S., Fasana A., Garibaldi L., Piombo B. A. D. Dynamic of multispan continuous straight bridges subject to multidegrees of freedom moving vehicle excitation. Journal of Sound Vibration, Vol. 224, Issue 3, 1999, p. 541561.

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. 7987.

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. 663683.

Lee S. Y., Yhim S. S. Dynamic behavior of longspan box girder bridges subjected to moving loads: Numerical analysis and experimental verification. International Journal of Solids and Structures, Vol. 42, Issues 1819, 2005, p. 50215035.

ISO 8606:1995(E). Mechanical Vibration – Road Surface Profiles – Reporting of Measured Data. 1995.

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. 309324.

Seetapan P., Chucheepsakul S. Dynamic responses of a twospan beam subjected to high speed 2DOF sprung vehicles. International Journal of Structural Stability and Dynamics, Vol. 6, Issue 3, 2006, p. 413430.

Da Silva J. G. S. Dynamical performance of highway bridge decks with irregular pavement surface. Computers and Structures,Vol. 82, Issues 1112, 2004, p. 871881.

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. 3145.

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. 6372.

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. 63986412.

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. 737749.

Wu S. Q., Law S. S. Dynamic analysis of bridgevehicle system with uncertainties based on the finite element model. Probabilistic Engineering Mechanics, Vol. 25, 2010, p. 425432.

Wu S. Q., Law S. S. Dynamic analysis of bridge with nonGaussian uncertainties under a moving vehicle. Probabilistic Engineering Mechanics, Vol. 26, 2011, p. 281293.

Liu N., Gao W., Song C., Zhang N., Pi Y. Interval dynamic response analysis of vehiclebridge interaction system with uncertainty. Journal of Sound Vibration, Vol. 332, 2013, p. 32183231.

Stefanou G. The stochastic finite element method: past, present and future. Computer Methods in Applied Mechanics and Engineering, Vol. 198, 2009, p. 10311051.

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. 315.

Schenk C. A., Schuëller G. I. Uncertainty Assessment of Large Finite Element Systems. SpringerVerlag, Berlin Heidelberg, 2005.

Ghanem R., Spanos P. D. Stochastic Finite Elements: A Spectral Approach. SpringerVerlag, New York, 1991.

Ghanem R. Stochastic finite elements with multiple random nonGaussian properties. Journal of Engineering MechanicsASCE, Vol. 125, Issue 1, 1999, p. 2640.

Field R. V. Numerical methods to estimate the coefficients of the polynomial chaos expansion. 15th ACSE Engineering Mechanics Conference, Columbia University, New York, 2002.

Reagan M. T., Najm H. N., Ghanem R., Knio O. M. Uncertainty quantification in reactingflow simulations through nonintrusive spectral projection. Combustion and Flame, Vol. 132, 2005, p. 545555.

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. 11911198.

Florian A. An efficient sampling scheme: updated latin hypercube sampling. Probabilistic Engineering Mechanics, Vol. 7, 1992, p. 123130.
About this article
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).