Numerical Calculation of Dynamic Response for Multi-Span Non-Uniform Beam Subjected to Moving Mass with Friction

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 matrices. 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 normal numerical computation. By applying these elements to calculate the non-uniform beam, the computational 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 simulated 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.


Introduction
Continuous beams are general statically indeterminate structures, and have broad applications in civil engineering, 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 beautiful, the holistic structures' stability is well, the spacial span is bigger and on which vehicles can placidly 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 researchers by now.Fryba [1,2] had given an exact solution on dynamic responses of the simple supported beam and continuous beam under moving load.Cai, Cheung and Chan [3] investigated the dynamic responses 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][5][6][7][8][9][10][11][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, Coriolis 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 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 deducted.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 computation.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.

Forced Vibration Differential Equation of Euler-Bernoulli Beam
Forced vibration differential equation of Euler-Bernoulli beam in plane bending with damping takes the form 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, c s is damping coefficient of strain rate, c(x) is damping coefficient of displacement rate.
When a concentrate excited force P(t) was exerted at x i ＝vt of the beam, the force can be expressed, by using the Dirac Delta function δ, as 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 [15].
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 centrifugal force on the beam, that is 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 expressed as follows When a moving vehicle has a brake on the bridge at a moment t n , if the contact friction is big enough, the vehicle will stop at moment t n+1 , the axial force X(ξ,t) can be expressed as ) where a 0 is the initial acceleration of the vehicle, μ and g represent the friction factor and that of gravity, respectively.If the moving vehicle does not brake on the bridge at any moment, the axial force X(ξ,t) will be expressed as The position ξ of the moving vehicle at the bridge can be expressed as follows (see Figure 1) where v 0 is the initial velocity, the position, velocity and acceleration at moment t n are, respectively,

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 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 P v i (ξ, t)=M v i g, and the mass force P a i (ξ, 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 M v i .Due to the location of each moving mass was continuously 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, N j (ξ), while j = 1,6, to get gravity load P v e (ξ, t), which is expressed as The action of the axial action force X i (ξ, t), at current position in each time interval Δt, can also be distributed to the element nodes.Now just the axial action was considered, the shear and moment were ignored.

 
The mass force P a i (ξ, t) can be discretized according to the form of element displacement interposition function y = ΣN j y j = Ny (j = 1,6).Consequently, one can obtain the additive moving mass matrix m a (t), moving damping matrix c a (t) and moving stiffness matrix k a (t) from Equation ( 22) where I is the unit matrix, P a e (ξ, t) is the equivalent nodal force vector inducing by the moving mass force, the shape function vector is a diagonal matrix N=diag(N 1 (ξ) N 2 (ξ) N 3 (ξ) N 4 (ξ) N 5 (ξ) N 6 (ξ)) .A detailed deducing process can be referenced in literature [13].
Instituting the additive moving mass matrix m a (t) and stiffness matrix k a (t) into the entire mass matrix M and stiffness matrix K of the original beam structure, respectively, the new entire matrix The overall damping matrix C of the beam is determined by using the theory of Rayleigh damping, adding the additive damping matrix c a (t) into C, the new entire damping matrix ( ) t C can be gotten by formula ( 25) 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 integration method or with Wilson-θ method.
where the load vector ( ) t F includes the gravity load of the moving vehicle P v (ξ, t) and the axial force X(ξ, t) as well as the adscititious load P(t) induced by other causations.

Stiffness and Mass Matrices of Non-Unform
Beam Elements

Non-Uniform Rectangular Cross Section
For a non-uniform rectangular section beams with both linear and parabolic variable heights in a plane, the stiffness 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 section height, area and the moment of inertia of the beam can be given by expressions (29)-(31), respectively.Copyright © 2010 SciRes.
where h 0 and h 1 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 obtained, respectively, from where B=∂ 2 N / ∂x 2 , while EI(x) and A(x) are the flexural stiffness and the cross section area of the beam, respectively, 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, k ij and m ij, can be obtained as follows

Non-Uniform Box Girder Section
The box section is shown in Figure 3 with an up bottom width of B and a down bottom width of D. The thicknesses of both up and down bottom board are T.The ventral shield thickness is C, and the section height of the box girder is h(x).The centroid distance w from centroidal axis z to z′ axis of self-defined is ( ) So the area and moment of inertia can be denoted by a section height h(x) and a centroid distance w(x) as follows Hypothesis that using a linear and parabolic variable mode such as Figure 2, the beam height h(x) of box girder element can also be denoted by formula 29, 30 and 31.Since w(x) is a composite function of h(x), so the moment of inertia I(x) is also a very complex composite function, moreover have some rational fractions in it.It is difficult to get a fixed form integral result by a manual calculation, so Romberg numerical integral scheme is adopted in this paper.The numerical integral precision is controlled to be 10 -6 .Each coefficient of the stiffness matrix of non-uniform section box girder element can be obtained by means of Formula (37).
When deducing the mass matrix of non-uniform section box girder element, since the area function A(x) is relatively simple, so we adopt a manual calculation fashion and the fixed form integral result is gotten.Observing Formula (35), and comparing the area formula of a rectangular cross section, a superfluous item (B+D-4C) T is found.These coefficients are constants.Form the general Expressions (33) of mass matrix, it is easily known that the mass matrix of non-uniform section box girder element is gotten, as long as substituting the parameter b into 2C and adding an item (such as Formula (38)) in mass matrix of rectangular cross section beam element.0 ( 4 ) For referring and using conveniently, based on before-mentioned three non-uniform section modes, the deduced results for all elements m ij of mass matrix of box girder element are enumerated as follows.

Mass Matrix Coefficients of the Box Girder with
Linear Variable Heights (Figure 2(a))

Mass Matrix Coefficients of the Box Girder with
Parabolic Variable Heights (Fig. 2(c), h 0 <h 1 )

Validation
In order to demonstrate the feasibility of the present stiffness and mass element matrices, a simple example is performed, that is, a non-uniform cross section beam with the length l = 20 m, the width of the rectangular section is b = 0.2 m, the small end and the big end height are h 1 = 0.1 m and h 0 = 0.3 m, respectively.Young's modulus is E = 3.0 × 10 10 N/m 2 , the concentrate load is F = 10 kN, the shape of the beam is shown in  First, taking half length of the beam and computing a cantilever beam with big end is fixed and small end is existed a concentrate force F = 10 kN, the numerical results is shown in Table 1.From the datum in Table 1, it shows that using the element in Figure 2(a) to compute the cantilever beam, one can obtained an analytic solution by taking the overall beam just as one element, whereas using the subsection uniform element one needs to divide the overall beam into 8-10 elements to gain the approximate solutions.Second, taking overall length of the beam and computing a simple supported beam with a concentrate force F = 10 kN at the mid-point of the beam, the computational results are shown in Figures 5 and 6.In above mentioned Figures, it is shown that the numerical results, by using the subsection uniform element, is close gradually to those of by using the paper present element, nevertheless, the convergence rate is decreased gradually with increase of the subsection numbers.According to the present element in Figures 2(b) and (c) to compute the parabolic cantilever beam and simple supported one with non-uniform cross section, the results are shown in Figures 7 and 8, respectively.From the computational results one can know that the accuracy of by using the parabolic elements is not as good as one of by using the tapered ones for the tapered beam, and the results are approximative.Nevertheless, comparing with using the subsection uniform section element, in the condition of ensuring definite computing precision, the needed element numbers for the parabolic beam is also small than that of using the subsection uniform element.By applying the present element to analyze the nonuniform rectangular beams with both linear and parabolic variable heights, the results are being approached the accurate solutions much more.A non-uniform cantilever beam element and a non-uniform simple supported beam are presented to validate the element's reliability, and the calculating results shows that, if using the subsection uniform finite elements one needs to divide more elements to converge the numerical solution to the current exact solution.Therefore, by applying the present element to analysis the non-uniform beam, the analysis can be simplified distinctly and the computational results will approach the accurate solutions.

A Three-Span Continuous Haunched Bridge under a Moving Load
In this example a three-span non-uniform continuous bridge is performed, the height of the beam is changeable in the plane (see Figure 9).A single load value of P v = 100 kN, moving at a speed of v ＝ 17 m/s, is considered.Total length of the bridge is L = 60 m, Young's modulus is E = 3 × 10 10 N/m 2 and the mass density is ρ ＝ 2400 kg/m 3 .The acceleration of the gravity is g = 9.81m/s 2 in all examples.The damping coefficient is ζ 1 = ζ 2 = 0.005 and the associated natural frequencies are ω 1 and ω 2 , which were obtained in the main dynamic program.First, a case of which the haunched bridge with damping subjected to a moving load, is considered.In this case the natural frequencies are not time-variant, because in this example the effects of inertial force, Coriolis force and centrifugal force induced by the moving mass have been ignored, so the mass, stiffness and damping matrices of the entire vibrating system are not time-variant yet.The finite element model of the bridge is composed of 36 uniform beam elements and 24 tapered beam elements.The numerical results for the deflections at each midspan position are in excellent agreement with those available ones from the reference literature [6,8].
Next considering the bridge with damping subjected to a moving mass, since the effect of the inertial force, Coriolis force and centrifugal force induced by the moving mass is existent, the overall matrices(include mass, stiffness and damping) and the associated natural frequencies are all time-variant at each computing time step.The numerical results for the deflections at each mid-span position of the bridge are shown in Figure 10.
From Figure 10 it is shown that the difference of the displacements, for the case with moving load and with moving mass, is no evident.Since the deflections of the bridge are small, the effects of inertial force, Coriolis force and centrifugal force induced by the moving mass are also small to the dynamical responses of the bridge.Nevertheless, the last effect may be significant for other cases.

A Simple Supported Beam under a Moving Mass with Uniform Variable Speeds
Considering a simple supported beam under a moving mass with uniform variable speeds (see Figure 11), the length of the beam is L = 100 m, Young's modulus is E = 2.15 × 10 11 Pa, the mass density of the beam is ρ = 6375 kg/m 3 and the section area is A = 2.4 m 2 , the moment of inertia is I = 0.8 m 4 , the moving mass m v = 61.2× 10 3 kg.We divided the beam into 50 beam elements, and taken the time step as Δt = 0.02 s, the initial velocity of the moving mass is v = 20 m/s, and the acceleration are a = 0, ±3, ±6 and ±9 m/s 2 , respectively.The dynamical displacement results of the mid-point of the beam have been computed by adopting the Newmark method.Under the conditions of being different accelerations (accelerating and decelerating), the dynamical displacement response contrast curves of the midpoint of the beam are given in Figures 12 and 13, respectively.From Figures 12 and 13 it is shown that, under the conditions of being same initial velocity, the dynamical displacements of the midpoint of the beam induced by moving mass with acceleration is bigger than that of with uniform velocity, and the larger the acceleration, the bigger the midpoint displacements of the beam.Whereas the dynamical displacements of the midpoint of the beam induced by moving mass with deceleration is smaller than that of with uniform velocity, and the larger the acceleration, the smaller the midpoint displacements of the beam.It is just opposite to the state with accelerated motion.It is because the moving mass was exerted by a friction induced between the contact interfaces.
When the moving mass motion with acceleration, it is given a friction, which is in line with the movement direction, by the beam.At the same time, the beam is subjected to a reaction force imposed by the moving mass.This force can be regarded as an axial pressure acting on the beam.This pressure will generate an additional bending moment in the beam, so that the deflection of the beam will be increased.On the contrary, when the moving mass deceleration movement, the beam will be subjected to an axial tension, and this tension will generate an additional bending moment within the beam to reduce the deflection.

A Three-Span Continuous Box Girder Bridge with Non-Uniform Section under a Moving Vehicle with Friction
The purpose of this example is to compute a three-span continuous box girder bridge with non-uniform sections under a moving vehicle with friction, the bridge is composed of 7 box girder segments, the height of the boxsection is changeable with both linear and parabolic soffit shape in the plane (see Figure 14).The total length of the bridge is L = 120m, the mass density is ρ ＝ 2400 kg/m 3 and Young's modulus is E = 3 × 10 10 N/m 2 .The acceleration of the gravity is g = 9.81m/s 2 , the moment of inertia I(x) is changeable with the position x.There is a vehicle at left end of the bridge, the weight of the vehicle is P = 20.6 × 10 4 N (while the front-wheel and rear wheel are P 1 = 8.0 × 10 3 N and P 2 = 12.6 × 10 3 N, respectively), the space between the wheels is 2.0 m, at the time of t = 0 the front-wheel is just in the left end of bridge, the vehicle travels at a speed of v ＝15 m/s.When the vehicle moves forward 30 meters from the left end of the bridge at a uniform velocity, then it is broken, so the velocity will be slowed down or stopped by contact friction.In order to impose the braking force, a ramp function is assumed (see Figure 15).This is based on the test results on highway vehicles conduced by the Transport and Road Research Laboratory [18].The braking force increases linearly to a maximum F b max = εP and then stays constant until the vehicle either comes to a stop or crosses the bridge span and is written as where ε is the impact coefficient, P is the vehicle static weight, t n is the moment of begin braking and T is the moment of the braking force increases to the maximum value of F bmax in time interval t n ≤ t ≤ t n+1 .
For this example, 60 linear and parabolic beam elements are divided in overall bridge, the time step is Δt ＝ 0.02 s.We take the friction coefficients as μ = 0.1, 0.2, 0.3 and 0.4, respectively.Taking the impact coefficient ε = 0.2, the numerical results for the midpoint deflections at mid-span position of the bridge are shown in Figure 16.
Figure 16 shows that the dynamical displacements were affected by friction together with the impact coefficient.The larger the friction coefficient, the larger the maximum dynamical displacements.The larger the friction, the shorter the time of braking vehicle.In the condition of being same friction force.The larger the impact coefficient, the bigger the dynamical displacements of midpoint at the mid-span of bridge.It can be observed from

Concluding Remarks
Dynamical solutions of the multi-span non-uniform beam with a moving vehicle can not be obtained by the theoretical means, the mass of the vehicle can not just be regarded as a simple moving load, the effect of inertial force, Coriolis force and centrifugal force is needed to be considered in general.The use of box beam bridge with non-uniform cross section is fairly common in the engineering.Numerical analysis for dynamic response of this kind of beam is beneficial to understand the dynamic characteristics of the bridge, to provide the scientific basis for the safe use of the bridge, and provided with a certain practical significance and application values.Using oblique-shaped and Hparabolic-shaped beam element with non-uniform cross section Hcan improve the accuracy and efficiency in solving.To define the dynamic response effect of the bridge caused by vehicle moving with variable speed as well as the friction and braking impact force on the bridge deck needs to considered the relationships of various factors synthetically, and then to find a design scheme closer to engineering practice combined with the relevant scientific experiment.

Figure 1 .
Figure 1.The position ξ of the moving vehicle at the beam.

Figure 3 .
Figure 3. Geometry size of cross section of the box girder.

Figure 6 .
Figure 6.Relationship between mid-point displacement and element numbers for the tapered simple supported beam.

Figure 7 .
Figure 7. Relationships between end-point displacements and element numbers for the parabolic cantilever beam.

Figure 8 .
Figure 8. Relationship between mid-point displacements and element numbers for the parabolic simple supported beam.

Figure 5 .
Figure 5. Relationship between end-point displacements and element numbers for the tapered cantilever beam.

Figure 9 .
Figure 9.A three-span continuous haunched bridge under a moving load.

Figure 10 .Figure 11 .
Figure 10.Deflections at each mid-point position for each span under a moving load/mass.

Figure 12 .Figure 13 .
Figure 12.Mid-point displacement of the beam with effect of acceleration.

Figure 14 .
Figure 14.A three-span continuous box girder bridge with variable sections.