Stochastic dynamic finite element analysis of bridge–vehicle system subjected to random material properties and loadings

Applied Mathematics and Computation 242 (2014) 20–35

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Stochastic dynamic finite element analysis of bridge–vehicle system subjected to random material properties and loadings T.-P. Chang Department of Construction Engineering, National Kaohsiung First University of Science and Technology, Kaohsiung, Taiwan

a r t i c l e

i n f o

Keywords: Bridge–vehicle system Laminated composite beam Random material properties Gaussian random process Finite element method Karhunen–Loéve expansion

a b s t r a c t This paper performs the statistical dynamic analysis on the bridge–vehicle interaction problem with randomness in material properties and moving loads. The bridge is modeled as a laminated composite beam with Gaussian random elastic modulus and mass density of material with random moving forces on top. The mathematical model of the bridge–vehicle system is established based on the finite element model in which the Gaussian random processes are represented by the Karhunen–Loéve expansion. Some statistical response such as the mean value and standard deviation of the deflections of the beam are obtained and checked by Monte Carlo simulation. Ó 2014 Elsevier Inc. All rights reserved.

1. Introduction Recently, the dynamic response of a bridge structure under moving loads has been investigated by many researchers. Fryba [1] obtained the analytical solutions for simply supported and continuous beams with uniform cross-section. Green and Cebon [2] presented the solution on the dynamic response of an Euler–Bernoulli beam in the frequency domain subjected to a ‘‘quarter-car’’ vehicle model by adopting an iterative procedure in conjunction with the experimental verification. Yang and Lin [3] obtained closed-form solution on the dynamic interaction between a moving vehicle and the supporting bridge using the modal superposition technique. Zheng et al. [4] performed the research on a multi-span non-uniform beam subjected to a moving load. Based on the Lagrange equation and modal superposition, the beam bridge model was extended by Zhu and Law [5] to an orthotropic rectangular plate under a series of moving loads. Marchesiello et al. [6] also proposed an analytical approach to the vehicle–bridge dynamic interaction problem with a seven degrees-of-freedom vehicle system moving on a multi-span continuous bridge deck. For more complex bridge–vehicle models, the finite element method has been utilized to perform the dynamic interaction analysis. Henchi et al. [7] presented an efficient algorithm for the dynamic analysis of a bridge discretized into threedimensional finite elements with a stream of vehicles running on top at a prescribed speed. The coupled equations of motion of the bridge–vehicle system are solved directly without an iterative method. Similar work has been conducted by Lee and Yhim [8] and Kim et al. [9] with experimental and field test data respectively. There are also other kinds of finite element methods, such as the ‘‘moving element method’’ [10] and ‘‘moving mass element method’’ [11,12], which are developed to solve the problem of moving forces on frame and plate structures. Generally speaking, the conventional deterministic analysis obtains only an approximation of the actual response due to uncertainties in the structural properties and the external loadings. Therefore, it is quite necessary to perform the stochastic analysis for the bridge–vehicle interaction problem. Recently, some research work has been conducted on the dynamic response of a bridge deck with the road surface roughness modeled as Gaussian random processes. Some researchers E-mail address: [email protected] http://dx.doi.org/10.1016/j.amc.2014.05.038 0096-3003/Ó 2014 Elsevier Inc. All rights reserved.

21

T.-P. Chang / Applied Mathematics and Computation 242 (2014) 20–35

[13–17] considered the randomness in the loadings due to the road surface roughness where the system parameters of both bridge and vehicle were treated as deterministic. Others had the randomness in the mass, stiffness, damping and moving velocity [18–20] of the moving vehicle and perturbation method was adopted to estimate the statistics of the structural dynamic response under random excitations. To model the uncertainties in structural analysis, stochastic finite element method is often used. Fryba et al. [21] computed the statistics of the dynamic response of a beam under a single moving force using the stochastic finite element analysis by using the first order perturbation in which the stiffness and damping were modeled as Gaussian random variables. Due to the fact that the perturbation method tends to lose accuracy when the level of uncertainty increases [22,23], the Karhunen–Loéve (K–L) expansion is used in the present study to represent the Gaussian random processes in the equation of motion of the bridge–vehicle system. The bridge response under Gaussian random vehicular forces which may have non-Gaussian properties will be approximated by Gaussian random processes. This approach is similar to that proposed by Ghanem and Spanos [24] noted as the spectral stochastic finite element method in which the random responses with nonGaussian properties are presented on a polynomial chaos basis. Nevertheless, a large numbers of Karhunen–Loéve components are required to represent the system parameters and excitation [25,26]. Recently, Wu and Law [27] proposed a new method to study the dynamic behaviors of bridge–vehicle system with uncertainties; their mathematical model is established based on finite element model in which the Gaussian random processes are presented by the Karhunen–Loéve expansion. To study the vehicle induced vibrations, a simply supported beam with various constitutive materials is often used. Laminated composite and functionally graded beams due to their appropriate properties are extensively adopted in the engineering structures. Therefore, it is quite essential to know the dynamic characteristics of these structures under the dynamic moving loads. Kadivar and Mohebpour [28] obtained a solution based on a FEM for three deformation theories to investigate the dynamic response of the orthotropic laminated composite beams with shear effect and rotary inertia under the actions of moving loads. Chen [29] derived the dynamic equilibrium equations of composite nonuniform beams made of anisotropic materials considering the effects of transverse shear deformations and structural damping using Hamilton’s principle. Lin and Chen [30] studied the problems of dynamic stability of spinning pre-twisted sandwich beams with a constrained damping layer subjected to periodic axial loads by using the FEM. Kim [31] presented the dynamic stability behavior of the damped laminated beam under the uniformly distributed subtangential forces using the finite element formulation. The effects of various boundary conditions, fiber orientation and external and internal damping have been studied. Simsek and Kocatürk [32] investigated the free and forced vibration behavior of a functionally graded (FG) beam subjected to a concentrated moving harmonic load. They derived the equations of motion by using Lagrange’s equations under the assumptions of Euler–Bernoulli beam theory. Simsek [33] studied the vibration of FG beams under a moving mass considering the different beam theories using the Lagrange’s equations. In addition, he studied the non-linear vibration analysis of a FG Timoshenko beam under action of a harmonic load [34]. Recently, Mohebpour et al. [35] investigated the dynamic response of composite laminated beams subjected to the moving oscillator by using an algorithm based on the finite element method. The first order shear deformation theory (FSDT) is assumed in their beam model. A stochastic finite element model is proposed in the present study for the dynamic response calculation of a bridge structure considering stochastic loading with inherent randomness in a bridge–vehicle system. The algorithm based on the proposed model can handle complex random excitation forces with large uncertainties and relatively small uncertainties in system parameters. The bridge is modeled as a simply supported laminated composite beam with Gaussian random elastic modulus and mass density of material and random moving forces on top. The forces have time-varying mean values and a coefficient of variation at each time instance, and they are considered as Gaussian random processes. The equation of motion of the bridge–vehicle system is presented using the Karhunen–Loéve expansion and the response statistics are obtained by solving the system equation of motion using the Newmark-b method. Some of the numerical results based on the proposed stochastic approach are checked by those obtained from the Monte Carlo simulation. 2. Mathematical formulation 2.1. Beam model A simply-supported laminated composite beam of length L, width b and thickness h is considered, with the coordinate system placed at the mid-plane of the laminate, with the moving loads on the top of the beam, as shown in Fig. 1. z

PN

P2

P1

v x

h

b

L

Fig. 1. Bridge–vehicle system: geometry and coordinate system of the laminated composite beam.

22

T.-P. Chang / Applied Mathematics and Computation 242 (2014) 20–35

All orthotropic layers of the laminate are assumed to be of an equal thickness. The fibers direction is indicated by angle h, which denotes the positive rotation angle of the principle material axes of the layer from the x-axes.u(x, t) and wðx; tÞ are the displacements at any point in the x- and z-directions, respectively. The displacement field for the laminated composite beam based on the first order shear deformation theory (FSDT) can be written as

uðx; z; tÞ ¼ u0 ðx; tÞ þ z/x ðx; tÞ;

v ðx; z; tÞ ¼ z/y ðx; tÞ;

ð1Þ

wðx; z; tÞ ¼ w0 ðx; tÞ; where u0 and w0 are the axial and transverse displacements of a point on the mid-plane in the x- and z-directions, /x and /y are the rotations of the normal to the mid-plane about the y- and x-axes, separately, and t is denoted as time. Based on the assumed displacement field, the material coupling among the extensional, bending and torsional deformation are incorporated in the formulation in order to maximum structural efficiency. The strain–displacement relationships for the theory are given by

 @u0 @/ @u þz x; e0xx ¼ 0 ; @x @x @x @/y @w0 cxy ¼ z ; cxz ¼ þ /x : @x @x

exx ¼

e1xx ¼

 @/x ; @x

ð2Þ

Based on the principle of virtual displacements, the following expression can be written

Z

T

ðdU þ dW  dDÞ dt ¼ 0;

ð3Þ

0

where dU is the virtual strain energy, dW is the virtual work done by applied forces, and dD denotes the virtual kinetic energy. The equation of motion of the beam can be derived from Eq. (3) as follows:

@Nxx @ 2 u0 @2/ ¼ I 0 2 þ I1 2x ; @x @t @t

ð4Þ

@Q x @ 2 w0 þ q þ Pc dðx  x0 Þ ¼ I0 ; @x @t2

ð5Þ

@Mxx @2/ @ 2 u0  Q X ¼ I2 2x þ I1 2 ; @x @t @t

ð6Þ

@ 2 /y @Mxy ¼ I2 ; @x @t 2

ð7Þ

where N xx ; M xx and M xy are the force and moment resultants per unit length, Ii are the mass inertia and Q x is the transverse shear force, q is the resultant distributed force at the surface of the laminate and P c is the concentrated applied force at position x0 : 2.2. Finite element formulation We can approximate the variables (u0 ; w0 ; /x ; /y ) by using the Lagrange interpolation functions as follows

u0 ðx; tÞ ¼ w0 ðx; tÞ ¼ /x ðx; tÞ ¼

n X ui ðtÞaei ðxÞ; i¼1 n X

wi ðtÞaei ðxÞ;

i¼1 n X

ð8Þ

r i ðtÞaei ðxÞ;

i¼1 n X si ðtÞaei ðxÞ; /y ðx; tÞ ¼ i¼1 e i ðxÞ

where a are the Lagrange interpolation functions and ui ; wi ; r i and si denote the nodal values of du0 ; dw0 ; d/x and d/y , respectively. By substituting Eq. (8) into the weak forms of Eqs. (4)–(7), we can obtain

€ e g þ ½K e fW e g  fF e g ¼ f0g; ½M e fW

ð9Þ

T.-P. Chang / Applied Mathematics and Computation 242 (2014) 20–35

23

where fW e g denotes the element nodal displacement vector, ½M e ; ½K e ; ½F e  denotes the element mass matrix, stiffness matrix and load vector respectively. 3. The system equation of motion Based on the element stiffness, mass matrices and element load vector given in Eq. (9) and the assumption of Rayleigh damping, the system equation of motion of the bridge can be written in matrix form as

€ þ Cd W _ þ Kd W ¼ HF; Md W

ð10Þ

_ where Md , Cd and Kd are the deterministic mass, damping and stiffness matrices of the bridge structure, respectively; W, W € are the deterministic nodal displacement, velocity and acceleration vectors of the structure respectively and HF is the and W equivalent nodal load vector of the bridge–vehicle interaction forces. The nodal responses of the bridge model under moving forces can be obtained by directly solving Eq. (10). The displacement of the bridge at position x and time t can then be expressed as:

wðx; tÞ ¼ HðxÞWðtÞ;

ð11Þ

HðxÞTi

HðxÞTi

where HðxÞ ¼ f0    0    0g and is the shape function of the beam structure. HðxÞ is a 1  n vector with zero entries except at the order of the ith beam element on which the position x is located. 4. Karhunen–Loéve expansion The Karhunen–Loéve expansion of a random process v ðx; bÞ is based on its covariance function Cðx1 ; x2 Þ which is bounded, symmetric and positive definite with the following spectral decomposition

Cðx1 ; x2 Þ ¼

1 X kn un ðx1 Þun ðx2 Þ;

ð12Þ

n¼0

where kn and un ðxÞ are the eigenvalue and eigenvector of the covariance kernel, respectively. Then the random process v ðx; bÞ can be written as

v ðx; bÞ ¼ v ðxÞ þ v~ ðx; bÞ ¼ v ðxÞ þ

1 X pffiffiffiffiffi g n ðbÞ kn un ðxÞ;

ð13Þ

n¼1

 ðxÞ denotes the expected value of v ðx; bÞ and g n ðbÞ is a set of uncorrelated Gaussian random variables. b denotes the where v random dimension. An explicit expression for g n ðbÞ can also be expressed as [24]

1 g n ðbÞ ¼ pffiffiffiffiffi kn

Z

v~ ðx; bÞun ðxÞdx:

ð14Þ

Based on the theory of K–L expansion, a m-dimensional stochastic process vector Vðt; bÞ can be denoted as

Vðt; bÞ ¼ fv 1 ðt; bÞ v 2 ðt; bÞ   

vm ðt; bÞgT

ð15Þ

~ i ðt; bÞ can be expressed as After the truncation at the N v th order according to Eq. (13), the random component of the v

v~ i ðt; bÞ ¼

Nv X

ðjÞ

g j ðbÞ xi ;

ð16Þ

j¼1 ðjÞ

where the where xi ðtÞ is of size 1  n representing the jth K–L components of the ith term in Vðt; bÞ, thus Vðt; bÞ becomes

Vðt; bÞ ¼ lv ðtÞ þ

Nv Nv X X g j ðbÞ xðjÞ ðtÞ ¼ g j ðbÞ xðjÞ j¼1

ð17Þ

j¼0 ðjÞ T

 1v 2    v  m gT ; xðjÞ ðtÞ ¼ fx1 x2    xm g . with the mean vector lv ðtÞ, g 0 ðbÞ ¼ 1, xð0Þ ¼ fv ðjÞ ðjÞ

5. System with Gaussian excitation and Gaussian system parameters The mass density qðx; bÞ, Young’s modulus Eðx; bÞ and damping cðx; bÞ of the beam are assumed as Gaussian random pro c and standard deviation rq ; rE ; rc , respectively, and the random components of them can be  ; E; cesses, with mean value q ~ ~c, respectively. The FEM equation of motion of the structure with random material properties and random ~ ; E; denoted as q excitations can be written as follows

€ bÞ þ CWðt; _ bÞ þ þKWðt; bÞ ¼ HFðt; bÞ; MWðt;

ð18Þ

24

T.-P. Chang / Applied Mathematics and Computation 242 (2014) 20–35

_ € bÞ and Wðt; bÞ are the random nodal displacement, velocity and acceleration vectors of the structure, where Wðt; bÞ; Wðt; respectively. M, C, K are the stochastic mass, damping and stiffness matrices of the bridge structure, respectively; ~ K ¼ Kd þ K; ~ K ~ C ¼ Cd þ C; ~ M; ~ C; ~ are the random components of the system mass, damping and stiffness M ¼ Md þ M; matrices, respectively, and they can be generated by assembling the corresponding elemental matrices accordingly. In addition, Rayleigh damping is assumed in the following form:

C ¼ c1 M þ c2 K;

ð19Þ

where c1 and c2 are constants. Based on K–L expansion, the system stiffness matrix K can be expressed as

K ¼ Kd þ

NE X g i ðbÞ Ki

ð20Þ

i¼1

where N E is the number of components in the K–L expansion for the Young’s modulus after truncation, and Ki can be assembled from the element stiffness matrix. Let K0 ¼ Kd , we obtain NE X g i ðbÞ Ki :

K¼

ð21Þ

i¼0

Similarly the system mass matrix can be expressed as

M¼

Nq X

g j ðbÞ Mj :

ð22Þ

j¼0

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

C¼

NC X

g k ðbÞ Ck ;

ð23Þ

k¼0

where N q and N C are the number of the components in the K–L expansion for mass density and damping after truncation, respectively, and N C ¼ N q þ N E . The random excitation force vector Fðt; bÞ can be expressed by the K–L expansion according to Eq. (17) as

Fðt; bÞ ¼

NF X ðnÞ g n ðbÞ f ðtÞ;

ð24Þ

n¼0 ðnÞ

ðnÞ

ðnÞ

ðnÞ T

where f ¼ ff 1 f2    f m g are the K–L components for the random moving forces and N F is the number of K–L components retained after truncation. Since the covariance matrix of the response is not known a priori, the K–L expansion cannot be performed according to Eq. (17) for the nodal displacement vector Wðt; bÞ. However, it is assumed that the randomness in system parameters is not large and the structural response can be approximated to have Gaussian property. Hence the response will take the form introduced in Eq. (17) as

Wðt; bÞ ¼

NW X g n ðbÞ qðnÞ ðtÞ;

ð25Þ

n¼0

where N W is the number of corresponding components qðnÞ ðtÞ which is determined by the number of K–L components for both the excitation forces and system parameters as N W ¼ N q þ N E þ N F Similarly, the nodal velocity vector and nodal acceleration vector can be expressed in similar fashion. By substituting Eqs. (21)–(25) into Eq. (18) and taking the inner product on both sides of the equation with g m ðbÞ and utilizing the orthogonal property, we can obtain

Table 1 pffiffiffiffiffiffiffiffiffiffiffiffiffi  ¼ xðL2 =hÞ q0 =E2 Þ of asymmetric cross-ply (0/90) beams for various boundary conditions. Non-dimensional fundamental frequencies ðx

a

Boundary conditiona

Aydogdu [36]

Mohebpour et al. [35]

Present study

CF SS CS CC FF

2.590 7.218 11.058 15.682 16.210

2.587 7.190 10.972 15.452 16.160

2.588 7.195 10.996 15.556 16.188

CF = clamped–free, SS = simply supported, CS = clamped–simply supported, CC = clamped–clamped, FF = free–free.

25

T.-P. Chang / Applied Mathematics and Computation 242 (2014) 20–35 Table 2 Fundamental natural frequencies (Hz) of symmetric angle-ply (h/h/h/h) composite beams. Solution type

h 0

30

60

90

Clamped–free Jun et al. [37] Mohebpour et al. [35] Present

278.4 278.4 278.4

137.9 137.9 137.9

74.7 74.8 74.8

74.2 74.2 74.2

Simply supported Jun et al. [37] Mohebpour et al. [35] Present

753.2 753.2 753.2

428.6 428.7 428.7

209.2 209.3 209.3

207.4 207.5 207.5

Clamped–simply supported Jun et al. [37] Mohebpour et al. [35] Krishnaswamy et al. [38] Present

1058.5 1058.6 1090.9 1058.8

592.4 592.5 629.2 592.9

323.2 323.2 325.8 323.3

320.6 320.7 323.0 320.8

Clamped–clamped Jun et al. [37] Mohebpour et al. [35] Krishnaswamy et al. [38] Present

1376.4 1376.5 1384.3 1376.9

811.2 811.3 818.3 811.6

462.9 463.0 467.4 463.1

459.1 459.2 463.7 459.3

V*=0.10, Ω *=0.0 1.8

ξi=0 %, PS

1.6

ξi=0 %, MCS ξi=2 %, PS

Mean value of w(L/2,t)/D

1.4

ξi=2 %, MCS

1.2

ξi=5 %, PS ξi=5 %, MCS

1 0.8 0.6 0.4 0.2 0

0

0.1

0.2

0.3

0.4

0.5 T

0.6

0.7

0.8

0.9

1

Fig. 2. Time history of mean value of midspan displacements for V  ¼ 0:10; X ¼ 0:0: PC = Present study, MCS = Monte Carlo simulation.

Nq NC NW X NW X NW X NE X X X € ðnÞ ðtÞ þ fg j ðbÞg n ðbÞg m ðbÞgMj q fg k ðbÞg n ðbÞg m ðbÞgCk q_ ðnÞ ðtÞ þ fg i ðbÞg n ðbÞg m ðbÞgKi qðnÞ ðtÞ n¼0 j¼0

¼ Hf

n¼0 k¼0 ðmÞ

ðtÞ:

n¼0 i¼0

ð26Þ

Rewriting Eq. (26) in matrix form,

^ qðtÞg ^ q ^ ^ € ðtÞg þ ½Cf _ ¼ ½HffðtÞg: ½Mf þ ½KfqðtÞg

ð27Þ

The nodal responses of the bridge are calculated by solving the system Eq. (27) using the Newmark-b method to obtain € ðnÞ ðtÞ for the responses. The mean and the variance of the nodal displacements can then the components qðnÞ ðtÞ, q_ ðnÞ ðtÞ and q be evaluated as

26

T.-P. Chang / Applied Mathematics and Computation 242 (2014) 20–35

V*=0.10, Ω *=0.0 0.18

ξi=0 %, PS

Standard deviation of w(L/2,t)/D

0.16

ξi=0 %, MCS ξi=2 %, PS

0.14

ξi=2 %, MCS

0.12

ξi=5 %, PS ξi=5 %, MCS

0.1 0.08 0.06 0.04 0.02 0

0

0.1

0.2

0.3

0.4

0.5 T

0.6

0.7

0.8

0.9

1

Fig. 3. Time history of standard deviation of midspan displacements for V  ¼ 0:10; X ¼ 0:0: PC = Present study, MCS = Monte Carlo simulation.

V*=0.10, Ω *=0.20 0.6 0.4

Mean value of w(L/2,t)/D

0.2 0 -0.2 -0.4

ξi=0 %, PS

-0.6

ξi=0 %, MCS ξi=2 %, PS ξi=2 %, MCS

-0.8

ξi=5 %, PS

-1 -1.2

ξi=5 %, MCS 0

0.1

0.2

0.3

0.4

0.5 T

0.6

0.7

0.8

0.9

1

Fig. 4. Time history of mean value of midspan displacements for V  ¼ 0:10; X ¼ 0:20: PC = Present study, MCS = Monte Carlo simulation.

MEANfWðtÞg ¼ qð0Þ ðtÞ;

VARfWðtÞg ¼

NW X 2 ðqðnÞ ðtÞÞ :

ð28Þ

n¼1

Thus the mean and variance of displacement at position x and time t can be obtained as

MEANfwðx; tÞg ¼ HðxÞqð0Þ ðtÞ; VARfwðx; tÞg ¼

NW X 2 ðHðxÞqðnÞ ðtÞÞ :

ð29Þ

n¼1

In addition, the mean and variance of random velocity and acceleration at position x and time t can be obtained readily.

27

T.-P. Chang / Applied Mathematics and Computation 242 (2014) 20–35

V*=0.10, Ω *=0.20 0.06

Standard deviation of w(L/2,t)/D

0.04 0.02 0 -0.02 -0.04

ξi=0 %, PS

-0.06

ξi=0 %, MCS ξi=2 %, PS ξi=2 %, MCS

-0.08

ξi=5 %, PS

-0.1 -0.12

ξi=5 %, MCS 0

0.1

0.2

0.3

0.4

0.5 T

0.6

0.7

0.8

0.9

1

Fig. 5. Time history of standard deviation of midspan displacements for V  ¼ 0:10; X ¼ 0:20 PC = Present study, MCS = Monte Carlo simulation.

V*=0.10, Ω *=0.50 1.5

ξi=0 % ξi=2 %

Mean value of w(L/2,t)/D

1

ξi=5 %

0.5

0

-0.5

-1

-1.5

0

0.1

0.2

0.3

0.4

0.5 T

0.6

0.7

0.8

0.9

1

Fig. 6. Mean value of midspan displacements for V  ¼ 0:10; X ¼ 0:50.

6. Numerical results and discussion First of all, in order to check the accuracy and validity of the proposed approach, two deterministic simplified examples are presented and compared to the other available numerical results in the literature. The shear correction factor K S is assumed to be 5/6 in both examples. The following material and geometric properties are adopted in the first example:

E1 ¼ 40 E2 ;

G12 ¼ G13 ¼ 0:6 E2 ;

G23 ¼ 0:5 E2 ;

m12 ¼ 0:25; L=h ¼ 20:

The non-dimensional fundamental frequencies of asymmetric cross-ply (0/90) laminated beams for various boundary conditions are presented in Table 1. Besides, the numerical results of the Ritz method solution by Aydogdu [36] and FEM solution by Mohebpour et al. [35] are cited and compared in the table. It can be found that the numerical results from the present study are in good agreements with those of Refs. [35,36]. In second example, the following material and geometric properties are used:

28

T.-P. Chang / Applied Mathematics and Computation 242 (2014) 20–35

V*=0.10, Ω *=1.00 10

ξi=0 %

8

ξi=2 %

Mean value of w(L/2,t)/D

6

ξi=5 %

4 2 0 -2 -4 -6 -8 -10

0

0.1

0.2

0.3

0.4

0.5 T

0.6

0.7

0.8

0.9

1

Fig. 7. Mean value of midspan displacements for V  ¼ 0:10; X ¼ 1:00.

E1 ¼ 1:448  1011 Pa; E2 ¼ 9:65  109 Pa; G12 ¼ G13 ¼ 4:14  109 Pa;

m12 ¼ m13 ¼ 0:3;

3

q0 ¼ 1389:23 kg=m ; L=h ¼ 15: In Table 2, the present numerical results for the fundamental natural frequencies of the four-layer symmetric angle-ply beams ðh=  h=  h=hÞ for different boundary conditions are compared with those of the FEM solutions by Jun et al. [37] and Mohebpour et al. [35] and the analytical solutions by Krishnaswamy et al. [38]. Once again, the present numerical results match with those of Refs. [35,37,38]. Now we are ready to investigate the statistical dynamic response of the composite beam subjected to moving loads. For simplicity, only one moving concentrated load is considered here, and the mean value of the moving load is assumed as P0 sin Xt; where P 0 is the amplitude of the moving load and X is the driving frequency of the moving load. The effect of the velocity of the moving harmonic load is represented by dimensionless velocity parameter V  as follows:

V  ¼ v P =v cr ;

ð30Þ V*=0.30, Ω *=0.20 1.2

ξi=0 %

1

ξi=2 % ξi=5 %

Mean value of w(L/2,t)/D

0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6

0

0.1

0.2

0.3

0.4

0.5 T

0.6

0.7

0.8

0.9

Fig. 8. Mean value of midspan displacements for V  ¼ 0:30; X ¼ 0:20.

1

29

T.-P. Chang / Applied Mathematics and Computation 242 (2014) 20–35

V*=0.30, Ω *=0.50 2

ξi=0 %

1.5

ξi=2 % ξi=5 %

Mean value of w(L/2,t)/D

1 0.5 0 -0.5 -1 -1.5 -2 -2.5

0

0.2

0.1

0.3

0.4

0.5 T

0.6

0.7

0.8

0.9

1

Fig. 9. Mean value of midspan displacements for V  ¼ 0:30; X ¼ 0:50.

V*=0.30, Ω *=1.00 3

ξi=0 % ξi=2 %

Mean value of w(L/2,t)/D

2

ξi=5 %

1

0

-1

-2

-3

0

0.1

0.2

0.3

0.4

0.5 T

0.6

0.7

0.8

0.9

1

Fig. 10. Mean value of midspan displacements for V  ¼ 0:30; X ¼ 1:00.

where

vP

is the velocity of the moving harmonic load and

v cr ¼ x1 L=p;

v cr

is the critical velocity denoted as

ð31Þ

where x1 is the fundamental natural frequency of the composite beam. The effect of the excitation frequency of the moving harmonic load X is represented by the frequency ratio X , where

X ¼ X=x1 :

ð32Þ

The dimensionless time T is defined as

T ¼ xP =L ¼ v P t=L:

ð33Þ

30

T.-P. Chang / Applied Mathematics and Computation 242 (2014) 20–35

V*=0.50, Ω *=0.20 1.6

ξi=0 %

1.4

ξi=2 %

1.2

ξi=5 %

Mean value of w(L/2,t)/D

1 0.8 0.6 0.4 0.2 0 -0.2 -0.4

0

0.1

0.2

0.3

0.4

0.5 T

0.6

0.7

0.8

0.9

1

Fig. 11. Mean value of midspan displacements for V  ¼ 0:50; X ¼ 0:20.

V*=0.50, Ω *=0.50 1

Mean value of w(L/2,t)/D

0.5

0

-0.5

-1

ξi=0 %

-1.5

ξi=2 % ξi=5 %

-2

0

0.1

0.2

0.3

0.4

0.5 T

0.6

0.7

0.8

0.9

1

Fig. 12. Mean value of midspan displacements for V  ¼ 0:50; X ¼ 0:50.

In the following numerical computations, the laminated composite beam under a harmonic moving load with constant velocity is dealt with using the FSDT and the proposed method. The following properties of the bridge model are adopted in the numerical computations:

L ¼ 40 m; A ¼ 4:8 m2 ;

I ¼ 2:5498 m4 :

The simply supported laminated beam (0/90) is considered having the material properties with mean values as follows:

E1 ¼ 1:448  1011 Pa; E2 ¼ 9:653  109 Pa; G12 ¼ G13 ¼ 4:136  109 Pa; 3

q0 ¼ 1389:23 kg=m :

m12 ¼ m13 ¼ 0:25;

31

T.-P. Chang / Applied Mathematics and Computation 242 (2014) 20–35

V*=0.50, Ω *=1.00 0.5

Mean value of w(L/2,t)/D

0

-0.5

-1

ξi=0 %

-1.5

ξi=2 % ξi=5 % -2

0.1

0

0.2

0.3

0.4

0.5 T

0.6

0.7

0.8

0.9

1

Fig. 13. Mean value of midspan displacements for V  ¼ 0:50; X ¼ 1:00.

V*=1.00, Ω *=0.20 1.5

ξi=0 % ξi=2 %

Mean value of w(L/2,t)/D

ξi=5 % 1

0.5

0

0

0.1

0.2

0.3

0.4

0.5 T

0.6

0.7

0.8

0.9

1

Fig. 14. Mean value of midspan displacements for V  ¼ 1:00; X ¼ 0:20.

The elastic modulus and mass density are assumed as Gaussian random process and they have the spatial correlation represented by an exponential auto-covariance function as:

Cðx1 ; x2 Þ ¼ r2 exp

  jx1  x2 j ;

a

ð34Þ

where r is the standard deviation of the system parameter, a and jx1  x2 j are the correlation length and the positive dislocation of two points in a spatial domain of interest which are set to be 1 m and 5 m, respectively. Both the random elastic modulus and mass density are assumed with the same spatial correlation and the same level of coefficient of variation (COV). In the present study, COV is assumed as 10% for the elastic modulus, mass density and the harmonic moving load. In addition, the damping ratio is assumed to be the same for all the modes.

32

T.-P. Chang / Applied Mathematics and Computation 242 (2014) 20–35

V*=1.00, Ω *=0.50 1.4

ξi=0 % ξi=2 %

1.2

Mean value of w(L/2,t)/D

ξi=5 % 1

0.8

0.6

0.4

0.2

0

0

0.1

0.2

0.3

0.4

0.5 T

0.6

0.7

0.8

0.9

1

Fig. 15. Mean value of midspan displacements for V  ¼ 1:00; X ¼ 0:50.

V*=1.00, Ω *=1.00 0.45

ξi=0 %

0.4

ξi=2 % ξi=5 %

Mean value of w(L/2,t)/D

0.35 0.3 0.25 0.2 0.15 0.1 0.05 0

0

0.1

0.2

0.3

0.4

0.5 T

0.6

0.7

0.8

0.9

1

Fig. 16. Mean value of midspan displacements for V  ¼ 1:00; X ¼ 1:00.

In Fig. 2, the time history of mean value of the non-dimensional dynamic deflections at the midspan are presented for various values of the damping ratio in the case of V  ¼ 0:10; X ¼ 0:0: It is noted that the non-dimensional dynamic deflections are normalized by the static deflection, D, of a beam under a point load P0 at the midspan. As it can be seen from Fig. 2, the non-dimensional dynamic deflections get smaller as the damping ratio increases, in addition, the numerical results based on the present study are checked by the Monte Carlo simulation and they are in excellent agreements. Fig. 3 shows the time history of standard deviation of the non-dimensional dynamic deflections at the midspan, the results from the Monte Carlo simulation are a little higher than those from the present study, especially for zero damping ratio. Nevertheless, they are still in fairly good agreements. Figs. 4 and 5 present the mean value and standard deviation of the non-dimensional dynamic deflections at the midspan individually for the case of V  ¼ 0:10; X ¼ 0:20: Once again, the numerical results from the present study match pretty well with those from Monte Carlo simulation as it can be found from Figs. 4 and 5. As the frequency ratio increases from 0.00 to 0.20, the maximum non-dimensional dynamic deflections at the midspan also increase, which is quite reasonable. For simplicity, the Monte Carlo simulation check won’t be

33

T.-P. Chang / Applied Mathematics and Computation 242 (2014) 20–35

Ω *=0.20 5.5

ξi=5 %

5

ξi=8 %

Mean value of Max.(w(x,t)/D)

4.5

ξi=10 %

4 3.5 3 2.5 2 1.5 1 0.5

0

0.1

0.2

0.3

0.4

0.5 V*

0.6

0.7

0.8

0.9

1

Fig. 17. Variation of mean value of max. dynamic displacements with velocity parameter V⁄ for X ¼ 0:20.

performed from now on. Figs. 6 and 7 present the mean value of the non-dimensional dynamic deflections at the midspan for various values of the damping ratio in the case of V  ¼ 0:10; X ¼ 0:50 and V  ¼ 0:10; X ¼ 1:00: According to the Figs. 2, 4, 6 and 7, the mean values of magnitudes of the non-dimensional dynamic deflections increase as the excitation frequency ratio increases from zero to resonance case. In Figs. 8–16, the time history of mean value of the non-dimensional dynamic deflections at the midspan are presented for various values of the excitation frequency ratio ðX ¼ 0:20; 0:50; 1:00Þ and dimensionless velocity parameter ðV  ¼ 0:30; 0:50; 1:00Þ: As it can be detected from these figures, especially for the relatively low value of the dimensionless velocity ratio, the magnitudes of the non-dimensional dynamic deflections increase as the excitation frequency ratio increases from zero to resonance case. In the case of the moving harmonic load, for a fixed value of the excitation frequency, the number of the oscillations decreases with the increase in the dimensionless velocity parameter. Furthermore, for all values of the excitation frequency ratio considered in the present study, the behaviors of the time histories of the mean values of the dynamic deflections are very similar to each other at the highest velocity parameter, i.e. V  ¼ 1:00: Figs. 17–19 shows the variation of the mean value of the maximum non-dimensional dynamic deflections with the dimensionless velocity parameter for various values of the excitation frequency ratio ð X ¼ 0:20; 0:50; 1:00Þ. Based on Ω *=0.50 5.5

ξi=5 %

5

ξi=8 %

Mean value of Max.(w(x,t)/D)

4.5

ξi=10 %

4 3.5 3 2.5 2 1.5 1 0.5

0

0.1

0.2

0.3

0.4

0.5 V*

0.6

0.7

0.8

0.9

1

Fig. 18. Variation of mean value of max. dynamic displacements with velocity parameter V⁄ for X ¼ 0:50.

34

T.-P. Chang / Applied Mathematics and Computation 242 (2014) 20–35

Ω *=1.00 7

ξi=5 % ξi=8 %

Mean value of Max.(w(x,t)/D)

6

ξi=10 % 5

4

3

2

1

0 0

0.1

0.2

0.3

0.4

0.5 V*

0.6

0.7

0.8

0.9

1

Fig. 19. Variation of mean value of max. dynamic displacements with velocity parameter V⁄ for X ¼ 1:00.

these figures, it is detected that except for the resonance case ð X ¼ 1:00Þ, The mean values of the maximum non-dimensional dynamic deflections get larger as the dimensionless velocity increases, until a certain value of the dimensionless velocity parameter; and after this value, an increase in the dimensionless velocity parameter tends to a decrease in the mean value of the non-dimensional dynamic deflections. It can be expected that the dynamic deflections become very large when the moving load travels at the critical velocity ð V  ¼ 1:00Þ: On the contrary, it can be seen from Figs. 17–19 that the critical velocity defined here does not correspond to the maximum dynamic deflections, and the maximum non-dimensional dynamic deflections occurred at different velocity parameter depending on the excitation frequency ratio. 7. Conclusions In the present study, we perform the statistical dynamic analysis on the bridge–vehicle interaction problem with randomness in material properties and moving loads. The bridge is modeled as a laminated composite beam with Gaussian random elastic modulus and mass density of material with moving forces on top. These forces are time varying with a coefficient of variation at each time instance and they are considered as Gaussian random processes. The first order shear deformation theory is assumed for the laminated composite beam model. The mathematical model of the bridge–vehicle system is established based on the finite element model in which the Gaussian random processes are represented by the Karhunen–Loéve expansion and the equations are solved by the Newmark-b method. Some statistical properties such as the mean value and standard deviation of the structural response are obtained, and the results based on the proposed method match with those from Monte Carlo simulation. Furthermore, the effects on the dynamic deflections due to different vehicle velocity and different driving frequency of the moving loads are investigated. It is concluded that the mean values of magnitudes of the nondimensional dynamic deflections increase as the excitation frequency ratio increases from zero to resonance case, especially for the relatively low value of the dimensionless velocity ratio. Besides, in the case of the moving harmonic load, for a fixed value of the excitation frequency, the number of the oscillations decreases with the increase in the dimensionless velocity parameter. Moreover, the maximum non-dimensional dynamic deflections occurred at different velocity parameter depending on the excitation frequency ratio. It should be noted that the statistical dynamic response based on the proposed method plays an important role in estimating the structural reliability of the bridge–vehicle system. Acknowledgment This research was partially supported by the National Science Council in Taiwan through Grant NSC-100-2221-E-327026. The author is grateful for this support. References [1] L. Fryba, Vibration of Solids and Structures under Moving Loads, third ed., Thomas Telford, London, 1999. [2] M.F. Green, D. Cebon, Dynamic response of highway bridges to heavy vehicle loads: theory and experimental validation, J. Sound Vib. 170 (1) (1994) 51–78. [3] Y.B. Yang, C.W. Lin, Vehicle–bridge interaction dynamic and potential applications, J. Sound Vib. 284 (1–2) (2005) 205–226.

T.-P. Chang / Applied Mathematics and Computation 242 (2014) 20–35

35

[4] D.Y. Zheng, Y.K. Cheung, F.T.K. Au, Y.S. Cheng, Vibration of multi-span non-uniform beams under moving loads by using modified beam vibration functions, J. Sound Vib. 212 (3) (1998) 455–467. [5] X.Q. Zhu, S.S. Law, Dynamic behavior of orthotropic rectangular plates under moving loads, J. Eng. Mech. ASCE 129 (1) (2003) 79–87. [6] S. Marchesiello, A. Fasana, L. Garibaldi, B.A.D. Piombo, Dynamic of multi-span continuous straight bridges subject to multi-degrees of freedom moving vehicle excitation, J. Sound Vib. 224 (3) (1999) 541–561. [7] K. Henchi, M. Fafard, M. Talbot, G. Dhatt, An efficient algorithm for dynamic analysis of bridges under moving vehicles using a coupled modal and physical components approach, J. Sound Vib. 212 (4) (1998) 663–683. [8] S.Y. Lee, S.S. Yhim, Dynamic behavior of long-span box girder bridges subjected to moving loads: numerical analysis and experimental verification, Int. J. Solids Struct. 42 (18–19) (2005) 5021–5035. [9] C.W. Kim, M. Kawatani, K.B. Kim, Three-dimensional dynamic analysis for bridge–vehicle interaction with roadway roughness, Comput. Struct. 83 (19– 20) (2005) 1627–1645. [10] C.G. Koh, J.S.Y. Ong, D.K.H. Chua, J. Feng, Moving element method for train-track dynamic, Int. J. Numer. Methods Eng. 56 (11) (2003) 1549–1567. [11] J.J. Wu, Vibration analysis of a portal under the action of a moving distributed mass using moving mass element, Int. J. Numer. Methods Eng. 62 (14) (2005) 2028–2052. [12] J.J. Wu, Use of moving distributed mass element for the dynamic analysis of a flat plate undergoing a moving distributed load, Int. J. Numer. Methods Eng. 71 (3) (2007) 347–362. [13] T.P. Chang, Y.N. Liu, Dynamic finite element analysis of a non-linear beam subjected to a moving load, Int. J. Solids Struct. 33 (1996) 1673–1688. [14] J.G.S. Da Silva, Dynamic performance of highway bridge decks with irregular pavement surface, Comput. Struct. 82 (2004) 871–881. [15] J.H. Lin, Response of a bridge to moving vehicle load, Can. J. Civil Eng. 33 (1) (2006) 49–57. [16] C.A. Schenk, L.A. Bergman, Response of continuous system with stochastically varying surface roughness to moving load, J. Eng. Mech. ASCE 129 (7) (2003) 759–768. [17] P. Seetapan, S. Chucheepsakul, Dynamic response of a two-span beam subjected to high speed 2DOF spring vehicles, Int. J. Struct. Stab. Dyn. 6 (3) (2006) 413–430. [18] G. Muscolino, S. Benfratello, A. Sidoti, Dynamics analysis of distributed parameter system subjected to a moving oscillator with random mass, velocity and acceleration, Probab. Eng. Mech. 17 (2002) 63–72. [19] T.P. Chang, G.L. Lin, E. Chang, Vibration analysis of a beam with an internal hinge subjected to a random moving oscillator, Int. J. Solids Struct. 43 (2006) 6398–6412. [20] T.P. Chang, M.F. Liu, H.W. O, Vibration analysis of a uniform beam traversed by a moving vehicle with random mass and random velocity, Struct. Eng. Mech. 31 (6) (2009) 737–749. [21] L. Fryba, S. Nakagiri, N. Yoshikawa, Stochastic finite elements for beam on a random foundation with uncertain damping under a moving force, J. Sound Vib. 163 (1) (2003) 31–45. [22] I. Elishakoff, Y.J. Ren, Finite Element Methods for Structures with Large Stochastic Variations, Oxford University Press, 2003. [23] G.I. Schuëller, Developments in stochastic structural mechanics, Arch. Appl. Mech. 75 (10–12) (2006) 755–773. [24] R. Ghanem, P.D. Spanos, Stochastic Finite Elements: A Spectral Approach, Springer-Verlag, New York Inc, 1991. [25] G. Stefanou, The stochastic finite element method: past, present and future, Comput. Methods Appl. M 198 (2009) 1031–1051. [26] C.A. Schenk, G.I. Schuëller, Uncertainty assessment of large finite element systems, Springer-Verlag, Berlin, Heidelberg, 2005. [27] S.Q. Wu, S.S. Law, Dynamic analysis of bridge–vehicle system with uncertainties based on the finite element model, Probab. Eng. Mech. 25 (2010) 425– 432. [28] M.H. Kadivar, S.R. Mohebpour, Finite element dynamic analysis of unsymmetric composite laminated beams with shear effect and rotary inertia under the action of moving loads, Finite Elem. Anal. Des. 29 (1998) 259–273. [29] C.N. Chen, Dynamic equilibrium equations of composite anisotropic beams considering the effects of transverse shear deformations and structural damping, Compos. Struct. 48 (2000) 287–303. [30] C.Y. Lin, L.W. Chen, Dynamic stability of spinning pre-twisted sandwich beams with a constrained damping layer subjected to periodic axial loads, Compos. Struct. 70 (2005) 275–286. [31] N.I. Kim, Dynamic stability behavior of damped laminated beam subjected to uniformly distributed subtangential forces, Compos. Struct. 92 (2010) 2768–2780. [32] M. Simsek, T. Kocatürk, Free and forced vibration of a functionally graded beam subjected to a concentrated moving harmonic load, Compos. Struct. 90 (2009) 465–473. [33] M. Simsek, Vibration analysis of a functionally graded beam under a moving mass by using different beam theories, Compos. Struct. 92 (2010) 904– 917. [34] M. Simsek, Non-linear vibration analysis of a functionally graded Timoshenko beam under action of a moving harmonic load, Compos. Struct. 92 (2010) 2532–2546. [35] S.R. Mohebpour, A.R. Fiouz, A.A. Ahmadzadeh, Dynamic investigation of laminated composite beams with shear and rotary inertia effect subjected to the moving oscillators using FEM, Compos. Struct. 93 (2011) 1118–1126. [36] M. Aydogdu, Vibration analysis of cross-ply laminated beams with general boundary conditions by Ritz method, Int. J. Mech. Sci. 47 (2005) 1730–1755. [37] L. Jun, H. Hongxing, S. Rongying, Dynamic finite element methods for generally laminated composite beams, Int. J. Mech. Sci. 50 (2008) 466–480. [38] S. Krishnaswamy, K. Chandrashekhara, W.Z.B. Wu, Analytic solutions to vibration of generally layered composite beams, J. Sound Vib. 159 (1992) 85– 89.