 Engineering, 2010, 2, 367-377 doi:10.4236/eng.2010.25048 Published Online May 2010 (http://www.SciRP.org/journal/eng) Copyright © 2010 SciRes. ENG 367Numerical Calculation of Dynamic Response for Multi-Span Non-Uniform Beam Subjected to Moving Mass with Friction Junping Pu, Peng Liu College of Civil Engineering & Architecture, Zhejiang University of Technology, Hangzhou, China E-mail: pjp@zjut.edu.cn, pliu85@126.com Received January 19, 2010; revised February 20, 2010; accepted February 2, 2010 Abstract In order to simulate the coupling vibration of a vehicle or train moves on a multi-span continuous bridge with non-uniform cross sections, a moving mass model is used according to the Finite Element Method, the effect of the inertial force, Coriolis force and centrifugal force are considered by means of the additive ma-trices. For a non-uniform rectangular section beam with both linear and parabolic variable heights in a plane, the stiffness and mass matrices of the beam elements are presented. For a non-uniform box girder, Romberg numerical integral scheme is adopted, each coefficient of the stiffness matrix is obtained by means of a nor-mal numerical computation. By applying these elements to calculate the non-uniform beam, the computa-tional accuracy and efficiency are improved. The finite element method program is worked out and an entire dynamic response process of the beam with non-uniform cross sections subjected to a moving mass is simu-lated numerically, the results are compared to those previously published for some simple examples. For some complex multi-span bridges subjected to some moving vehicles with changeable velocity and friction, the computational results, which can be regarded as a reference for engineering design and scientific research, are also given simultaneously. Keywords: Dynamic Response, Multi-Span Beam, Non-Uniform Section, Friction, Braking Force 1. Introduction Continuous beams are general statically indeterminate structures, and have broad applications in civil engi-neering, mechanism, navigation engineering and so on. Multi-span continuous bridges have been widely used in highway and railway, there is a great deal of merit for the structures, for example, their exterior is beau-tiful, the holistic structures’ stability is well, the spacial span is bigger and on which vehicles can plac-idly pass over. It is of great importance to study the dynamic characteristic of the bridge under moving mass for engineering design and scientific research. Many engineers and scientists have contributed to the solution of the problem with their innovations, and still the subject draws considerable attention from re-searchers by now. Fryba [1,2] had given an exact so-lution on dynamic responses of the simple supported beam and continuous beam under moving load. Cai, Cheung and Chan  investigated the dynamic re-sponses of the infinite continuous beam subjected to a moving force by using the mode superposition method to get an exact solution. However, for a great number of bridge structures in engineering practice, it can not be simply regarded as a perfect state, so their theoretic solutions are not existent in general. Thus, among the numerical methods the finite element method becomes an ideal approximate approach to solve this kind of problems [4-12]. Wu [13,14] performed the dynamic analysis of an inclined beam and a flat plate due to moving loads, and presented the moving mass element by taking account of the effect of inertial force, Corio-lis force and centrifugal force induced by the moving mass. Form mentioned literatures, it is shown that, for the multi-span continuous non-uniform beam, one used a moving load model to obtain the numerical J. P. PU ET AL. 368 results in general, but did not consider the effect of inertial force, Coriolis force and centrifugal force. Some numerical examples dealt with above effect but just analyzed the simple uniform beam. This paper has been performed some complex problems (include multi-span non-uniform beam with moving mass). For a non-uniform rectangular section beam with both linear and parabolic variable heights in a plane, the stiffness and mass matrices of the beam elements have been de-ducted. For a non-uniform box girder which can be widely used in engineering structures, since the integral formula of stiffness matrix is extremely complex, it is difficult to write down the expression of the stiffness coefficients, so Romberg numerical integral scheme is adopted. Each coefficient of the stiffness matrix can be obtained by means of a normal numerical computa-tion. For some complex multi-span bridges subjected to some moving vehicles with changeable velocities and frictions, numerical results of dynamic responses are also obtained, which can be regarded as a reference for engineering design and scientific research. 2. Forced Vibration Differential Equation of Euler-Bernoulli Beam Forced vibration differential equation of Euler-Bernoulli beam in plane bending with damping takes the form ),()()()()(2223222222txftyxctyxmxtyxIxcxyxEIxs (1) where, the mass per unit length of the beam is m(x) = ρA, while ρ is the material density and A is the area of the cross section of the beam, cs is damping coefficient of strain rate, c(x) is damping coefficient of displacement rate. When a concentrate excited force P(t) was exerted at xi＝vt of the beam, the force can be expressed, by using the Dirac Delta function δ, as (,) ()()ifxtPtxvt (2) where v is the moving velocity. When a moving mass M with the constant velocity moved on the beam, the action force on the beam can be considered as a common effect of the gravity and inertia force according to reference literature . )()],([),(vtxtvtygMtxf  (3) Thus the formula (1) will be a set of systems of second order differential equations with time-variant coefficient, it needs to be solved by means of numerical method. Lee’s investigation indicated, in literature [16,17], that when the moving mass M with the constant velocity moved on the Euler-Bernoulli beam, including the grav-ity and inertia force induced by moving mass M, one still needs to consider the effect of Coriolis force and cen-trifugal force on the beam, that is )()],(),(2),([),(2vtxtvtyvtvtyvtvtygMtxf (4) Considering the contact friction between the wheels and the beam, a axial force X(ξ,t) should be exerted on the beam, thus the action force on the beam can be ex-pressed as follows )(),()],(),(2),([),( 2xtXtyvtyvtygMtxf (5) When a moving vehicle has a brake on the bridge at a moment tn, if the contact friction is big enough, the vehi-cle will stop at moment tn+1, the axial force X(ξ,t) can be expressed as 0),( MatX , (6) )0( ntt )(),( 0agMtX , (7) )(1 nn ttt0),(tX . (8) )( 1nttwhere a0 is the initial acceleration of the vehicle, μ and g represent the friction factor and that of gravity, respec-tively. If the moving vehicle does not brake on the bridge at any moment, the axial force X(ξ,t) will be expressed as 0),( MatX (9) )0(tThe position ξ of the moving vehicle at the bridge can be expressed as follows (see Figure 1) 2000 5.0 tatvx  , (10) )( ntt2)(5.0)( nnnnn ttattvx  )(1nn ttt (11) nnn avx2/2 . (12) )( 1nttwhere v0 is the initial velocity, the position, velocity and acceleration at moment tn are, respectively, 2000 5.0 nnn tatvxx  , (13) nn tavv 00  , (14) gaan0. (15) 3. Discrete Model of Vibration Equations under Moving Mass According to the finite element method, the forced vibration differential equation of Euler-Bernoulli beam in plane bending with damping can be written, in matrix form, as follows )()()()( tttt FKyyCyM  (16) where ),(),(),()( tttt avXPPF (17) Copyright © 2010 SciRes. ENG J. P. PU ET AL.369 nttnxP00 ,av nn av,00,0 11 nnavnnnn avtt/1 Figure 1. The position ξ of the moving vehicle at the beam. The action force of the vehicle on the beam can be expressed by the nodal load vector F(t), which is composed of the gravity Pvi(ξ, t)=Mvig, and the mass force Pai (ξ, t) (include inertia force, Coriolis force and centrifugal force) of the vehicle, as well as the axial action force X i(ξ, t), while i=1,n denotes the number of moving mass Mvi. Due to the location of each moving mass was con-tinuously transformed with time, while the moving mass passed along each beam element, the action of the gravity, at current position in each time interval Δt, can be distributed to the element nodes to become the element nodes’ forces by using an interposition function. For Euler-Bernoulli beam element, one can adopt cubic Hermite interposition function of two-nodes, Nj(ξ), while j = 1,6, to get gravity load Pve(ξ, t), which is expressed as 654321),( vvvvvvevFFFFFFt P (18) where )6,,1(),()(),(  jtPNtF ivjjv (19) The action of the axial action force X i(ξ, t), at cur-rent position in each time interval Δt, can also be dis-tributed to the element nodes. Now just the axial action was considered, the shear and moment were ignored. 0000),(41XXteX (20) where (,)() (,)(1,4)jijXtNXtj (21) The mass force Pai(ξ, t) can be discretized according to the form of element displacement interposition func-tion y = ΣNjyj= Ny (j = 1,6). Consequently, one can ob-tain the additive moving mass matrix ma(t), moving damping matrix ca(t) and moving stiffness matrix ka(t) from Equation (22) ykycymINP aaaeat ),( (22) where I is the unit matrix, Pae(ξ, t) is the equivalent nodal force vector inducing by the moving mass force, the shape function vector is a diagonal matrix N=diag(N1(ξ) N2(ξ) N3(ξ) N4(ξ) N5(ξ) N6(ξ)) . A detailed deduc-ing process can be referenced in literature . Instituting the additive moving mass matrix ma(t) and stiffness matrix ka(t) into the entire mass matrix M and stiffness matrix K of the original beam structure, respectively, the new entire matrix)(tMand ()tKare formed and are time-variant. )()( tt amMM  (23) () ()atKKkt (24) The overall damping matrix C of the beam is de-termined by using the theory of Rayleigh damping, adding the additive damping matrix ca(t) into C, the new entire damping matrix ()tCcan be gotten by for-mula (25) )()()()()()( tttttt acKMC  (25) where )/()(2)( 2122122121t (26) )/()(2)(21221122t (27) The coefficients α(t) and β(t) are also time-variant with changing of the natural frequencies ω1 and ω2 of the beam structure at each time steps. Finally, according to Equation (28) one can numerically compute an entire system of vibration equation with Newmark direct inte-gration method or with Wilson-θ method. )()()()()()()( ttttttt FyKyCyM   (28) where the load vector()tFincludes the gravity load of the moving vehicle Pv(ξ, t) and the axial force X(ξ, t) as well as the adscititious load P(t) induced by other cau-sations. 4. Stiffness and Mass Matrices of Non-Unform Beam Elements 4.1. Non-Uniform Rectangular Cross Section For a non-uniform rectangular section beams with both linear and parabolic variable heights in a plane, the stiff-ness and mass matrices are deducted, respectively. So one can analyze the non-uniform beam according to a convenient mode (see Figure 2). According to Figures 2(a), (b) and (c), the cross sec-tion height, area and the moment of inertia of the beam can be given by expressions (29)-(31), respectively.  0130111/hx hhIx bhh12 (29) h0h1 h(x) lh0 h1h(x) lh0h1 h(x) l (a) (b) (c) Figure 2. Variable cross-sectional beam elements. Copyright © 2010 SciRes. ENG J. P. PU ET AL. 370 0130111/hx hhIx bhh12 (30) 01301111111hx hhIx bhh/123l (31) where h0 and h1 are beginning and end height of the beam elements with non-uniform cross sections, respectively. ξ=x/l, while l is the length of the beam element and b is the width of the cross section with rectangular form. The element stiffness and mass matrices can be ob-tained, respectively, from T0dleEI xxKB B (32)   lxxAV 0TTVed)(deNNNNM (33) where B=∂2N / ∂x2, while EI(x) and A(x) are the flexural stiffness and the cross section area of the beam, respec-tively, which are changeable with the height changing of the cross section. Adopting above-mentioned non-uniform beam models and an integral procedure is worked out, the coefficients of the elemental stiffness and mass matrices, kij and mij, can be obtained as follows 4.1.1. Stiffness Matrix Coefficients with Linear Varable Heights (Figure 2(a)) 1144140 132 232225550010113223223350010113223226560010 1132 23330010 113360 0()/2(7337) / 20(522) /20(225)/ 20(11522)/ 60(4kk kEbhhlkkkEbhhhhhhkkEbhhhhhhlkkEbhhhhhhlkEbhhhhh hlkEbhh    223101 1322 36600101112131516 2434 45 464)/60(22511) / 600hhh hlkEbhhhhh hlkkkkkkkk 4.1.2. Mass Matrix Coefficients with Linear Variable Heights (Figure 2(a)) 110 11401220 1223 01250 1226 01333 01235 013/12/12(23/ 5)/ 7(/28/60)9() /140(/ 60/ 70)(/168/ 280)(/ 70/ 60)mblhhmblhhmhhblmh h blmhhblmhhblmh hblmhhbl  3360 1440155()/2803/12(3m hhblmblhhmh 01256 01366 0112 13 15 1624344546/52 )/7(/60/28)(/ 280/168)0hblmh h blmhh blmmmmmmmm  4.1.3. Stiffness Matrix Coefficients with Parabolic Variable Heights (Figure 2(b), h0>h1) 11 44140132 232225550010 113223 223350010 11(2)/3(37423992) / 210(26302128) /210kkkEbh hlkkkEbh hhhh hkkEbhhhhhh l   3l2l 322326560010 1132 2333001011322 3360010 11322 3660010 1112131516 24 34 4546(11121864) / 210(37463324) /420(1514932) / 420(7102796) /4200kkEbhhhhhhkEbhhh hhhlkEbh hhhhhlkEbh hh hhhlkkkkkkkk  4.1.4. Mass Matrix Coefficients with Parabolic Varable Heights (Figure 2(b), h0>h1) 1101140 1220 12230 1250 1226 013330135 01916/1051124/ 210(93077424)/ 45045(21592560)/ 90090(35518032)/ 90090(373/36036928/45045)(173256)/ 45045(103/ 120mblhhmblhhmhhblmhhblmhhblmhhblmhhblmh 213360 1440155 0 1256 013660 112131516 24 34 45 461216/ 715)(391896 )/1801806/21(7194858 )/15015(19/ 200292/ 2145)(31112 )/150150hblmhhblmblhhmhhblmhh blmhhblmmmmmmmm   4.1.5 Stiffness Matrix Coefficients with Parabolic Variable Heights (Figure 2(c), h0h1) 110 1140 14401220 12230 1252 1916354/10521124354/ 210267 4/21293077424167314/4504522159256047194/ 900902 3551mlChh TBDCmlChh TBDCmlChhTBDCmlC hhTBDCmlC hhTBDCmlC012260 13330 12350 1336018032115834/ 9009021865 371255774 /18018021732564294/ 450452 1545403255774/1801802391 89612874hh TBDCmlC hhTBDCmlCh hTBDCmlC hhTBDCmlC hhTBDC  55 02560 13660 11213151624 34 45 46/1801802 719485855774/150152855386447194/ 900902 311121434/150150mlCh h TBDCmlC hhTBDCmlCh hTBDCmmmmmmmm   4.2.3 Mass Matrix Coefficients of the Box Girder with Parabolic Variable Heights (Figure 2(c), h0