# Coupled Dynamics of Vehicle-Bridge Interaction System Using High Efficiency Method

Authors: Lu Sun; Xingzhuang Zhao
Journal: Advances in Civil Engineering
 (2021)

---

## Abstract

Vehicle-bridge interaction is the core for a variety of applications, including vehicle vibration, bridge vibration, bridge structural health monitoring, weight-in-motion, bridge condition inspection, and load rating. These applications give rise to a great interest in pursuing a high-efficiency method that can tackle intensive computation in the context of vehicle-bridge interaction. This paper studies the accuracy and efficiency of discretizing the beam in space as lumped masses using the flexibility method and as finite elements using the stiffness method. Computational complexity analysis is carried out along with a numerical case study to compare the accuracy and efficiency of both methods against the analytical solutions. It is found that both methods result in a similar level of accuracy, but the flexibility method overperforms the stiffness method in terms of computational efficiency. This high efficiency algorithm and corresponding discretization schema are applied to study the dynamics of vehicle-bridge interaction. A system of coupled equations is solved directly for a simply supported single-span bridge and a four-degree-of-freedom vehicle modeling. Pavement roughness significantly influences dynamic load coefficient, suggesting preventative maintenance or timely maintenance of pavement surface on a bridge, to reduce pavement roughness, is of significant importance for bridge’s longevity and life-cycle cost benefit. For class A and B level pavement roughness, the dynamic load coefficient is simulated within 2.0, compatible with specifications of AASHTO standard, Australian standard, and Switzerland standard. However, the Chinese code underestimates the dynamic load coefficient for a bridge with a fundamental frequency of around 4 Hz. The proposed method is applicable to different types of bridges as well as train-bridge interaction.


---

## Body

## 1. Introduction

Vehicle-bridge interaction is the key to obtaining the interaction force between the tire and pavement. It has many applications in the weigh-in-motion system [1–32], bridge dynamics and vibration control [33–50], bridge structural health monitoring, bridge ambient testing, bridge damage identification [51–64], and high-speed railway [29, 65–68].Due to the complexity of models of vehicle-bridge interaction, an analytical solution of vehicle-bridge interaction is normally beyond achievement for most cases. Therefore, numerical methods such as the finite difference method and finite element method (FEM) become indispensable for any realistic applications [1, 2, 69–76].One of the essential components in modeling vehicle-bridge interaction is the bridge subsystem. In FEM, a bridge is typically discretized as finite elements and solved by the stiffness method [2, 4, 41, 73, 76–84]. Rieker et al. studied the influence of discretization of beams on the dynamic response of elastic beams under moving load [1]. And, they concluded that an acceptable level of accuracy (i.e., relative error within 1%) could be achieved when the number of elements in discretization for dynamic analysis should be 16 and be at least two to eight times greater than that used in static analysis. Law et al. [77] investigated the influence of the number of finite elements and argued that at least 8 elements were required to discretize the beam bridge.The abovementioned applications require intensive computation of time-varying vehicle-bridge interaction, demanding high-efficiency methods. What is missing in the literature on vehicle-bridge interaction is the investigation of an efficient and accurate method for solving the vehicle-bridge interaction problem.The innovative contributions of this study are three-fold. First of all, this study studies the accuracy and efficiency of discretizing the beam in space as lumped masses using the flexibility method and as finite elements using the stiffness method. Using computational complexity analysis, it is shown that both methods result in a similar level of accuracy, but the flexibility method is two to eight times faster than the stiffness method.Secondly, based on the new finding of the advantage embedded in discretizing the beam in space as lumped masses using the flexibility method, this study develops a high efficiency method for solving the coupled dynamics of vehicle-bridge interaction using a 4-DOF vehicle model, a simply supported single-span bridge model and ISO pavement roughness model.Lastly, this study reveals that pavement roughness significantly influences dynamic load coefficient (DLC) and Chinese code [85] underestimates the DLC for a bridge with fundamental frequency around 4 Hz, suggesting the importance of preventative and timely maintenance of the pavement surface on a bridge to reduce pavement unevenness for bridge’s longevity and life cycle cost benefit.


## 2. Equation of Beam Motion

Without forfeiting the validity of the analysis, a moving vehicle is simplified as a dynamic loadq(x, t) in this study [86]. This simplification allows analytical solutions to be developed for a simple form of the dynamic load, therefore serving as a benchmark for comparison against the numerical solution.In this section, the equation of motion of a beam is derived. Consider an Euler–Bernoulli beam, with unit width and lengthL subjected to an arbitrarily distributed force qx,t, as shown in Figure 1. The parameters of the beam are E which denotes Young’s modulus, Iz which denotes the cross-section moment of inertia, ρ which denotes the density, and A which denotes the area of cross section. It is assumed that the beam bridge is uniform in cross section.Figure 1
Euler–Bernoulli beam.Taking out an infinitesimal segment from the beam, with which the force applied on this segment is shown in Figure2. Here, M denotes the moment applied on left cross section, V denotes the shear force applied on left cross section, M+∂M/∂xδx denotes the moment applied on right cross section, V+∂V/∂xδx denotes the shear force applied on right cross section, qx,t denotes the distributed force applied, m denotes the mass of the infinitesimal segment volume, m=ρAδx,a denotes the acceleration of the infinitesimal segment volume, a=∂2w/∂t2, where w is the deflection of the beam and t is time variable, fd is the unit damping force caused by the ambiance and material viscosity, which is proportional to the velocity, fd=c∂w/∂t, where d is the damping coefficient of the unit length.Figure 2
Infinitesimal segment analysis.Similar to [87], it is plausible to assume that the distributed force be uniform and the inertial force is applied at the center of the segment because of the infinitesimal length of the segment. The equation of beam motion can be established by Newton’s second law. Different from [87], an additional term of damping force proportional to the speed of beam motion is added to the equation.Apply Newton’s second law on the segment in the vertical direction:(1a)ρAδx∂2w∂t2=V+∂V∂xδx−V+qx,tδx−c∂w∂tδx,or equivalently(1b)ρA∂2w∂t2=∂V∂x+qx,t−c∂w∂t.Take moment with respect to point O. The moment balance equation can be obtained as follows:(2)M+Vδx−M−qx,tδx22−∂M∂xδx−maδx2=0.Substitutingm=ρAδx and a=∂2w/∂t2 into (2) leads to(3)Vδx−qx,tδx22−∂M∂xδx−ρAδx22=0.Divide both sides byδx and take the limit of δx⟶0, obtaining(4)V=∂M∂x.Substitution ofV into (1b) leads to(5)ρA∂2w∂t2=∂2M∂x2+qx,t−c∂w∂t.According to Euler–Bernoulli beam theory,(6)EIz∂2w∂x2=−M.The equation of motion for the beam bridge is obtained as follows:(7)ρA∂2w∂t2+c∂w∂t+∂2∂x2EIz∂2w∂x2=qx,t,where ρ,A,d, and EIz could be functions of both position x and time t, that is, ρ=ρx,t,A=Ax,t,c=cx,t, and∂2/∂x2EIz∂2w/∂x2=rx,t∂4w/∂x4,rx,t≠0. This generalization leads to(8)px,t∂2w∂t2+cx,t∂w∂t+rx,t∂4w∂x4=qx,t.The coefficient of the third term∂4w/∂x4 in the continuous equation of motion (8) could be set to 1 by dividing rx,t into both sides of (8), resulting in(9)p˜x,t∂2w∂t2+b˜x,t∂w∂t+∂4w∂x4=q˜x,t.For a uniform beam with a constantEIz, (7) can be further simplified as(10)ρA∂2w∂t2+c∂w∂t+EIz∂4w∂x4=qx,t,where p˜x,t=px,t/rx,t, b˜x,t=bx,t/rx,t, and q˜x,t=qx,t/rx,t.When neglecting the damping, i.e.,c=0, (8) degrades to Euler–Lagrange equation [88]. For a beam bridge and the initial time is set at t = 0, equations (7), (8), and (10) are defined for t≥0,0≤x≤L, and L is the length of the bridge. The coefficients of ∂4w/∂x4, ∂w/∂t, and ∂2w/∂t2 present stiffness, damping, and mass, respectively.


## 3. Discretization of Beam as Lumped Masses

A beam bridge can be discretized into lumped masses [89]. Each lumped mass is assigned to the node with an interval of Δx, as shown in Figure 3.Figure 3
Discretization schema of the beam as lumped masses.In the lumped mass discretization schema, only the vertical displacement at each of the lumped masses is considered. The deflection of the lumped mass can be obtained for both determinate and indeterminate structures [90]. For instance, the deflection of a simply supported beam due to static force is governed by the following formula given in the AISC steel construction manual [91]:(11)w=−Fbx6EIzLL2−x2−b2,w=−Fb6EIzLLbx−a3+L2−b2x−x3,where F is the concentrated force acting on the beam, x,E,Iz,andL are defined the same way as before, a is the distance from the origin, and b=L−a.The deflection of every lumped mass is obtained by applying a unit force to the mass to achieve the corresponding deflection. The stiffness matrix K in equation (13) can be computed as the reverse of the flexibility matrix Fd [84]. The mass matrix is obtained by constructing a diagonal matrix of size n×n, with the element of value ρxi,tAxi,tΔx. The force vector Ft applied on the lumped mass could be directly represented by qxi,tΔx.


## 4. Discretization of Beam as Finite Elements

The beam bridge could also be discretized asn finite elements, with Δx being the length of the elements, as shown in Figure 4. The two ends of each element are defined as the nodes of the element. From the figure, it is observed that the distances between the nodes are Δx.Figure 4
Discretization schema of the beam as finite elements.MultiplyingΔx on both sides of equation (1a), discretizing the system into(12)ρxi,tAxi,tΔx∂2wi∂t2+cxi,tΔx∂wi∂t+Exi,tIzxi,t∂4wi∂x4Δx=qxi,tΔxi=1,…,n.Equation (12) can be written in a matrix form as(13)Mtw¨+Ctw˙+Ktw=Ft.Here,ρxi,tAxi,tΔx∂2wi/∂t2 is equivalent to Mtw¨, cxi,tΔx∂wi/∂t is equivalent to Ctw˙, qxi,tΔx is equivalent to Ft, EtIzt∂4wi/∂x4Δx is equivalent to Ktw, Mtis the mass matrix of the discretized system, Ct is the damping coefficient matrix, Kt is the stiffness matrix, Ft=qxi,tΔx is the excitation force that is applied on each lumped mass, and w is the deflection vector of the beam.In Figure5, the nodes of the element are located on two clamped ends. An external exciting force applied to the element is transferred onto the nodes of the elements. For a concentrated load, as shown in Figure 5, the equivalent node force vector in the local coordinate is −Pb/Δx,−Pab2/Δx,−Pa/Δx,Pab2/ΔxT. For a distributed load, as shown in Figure 6, the equivalent node force vector in the local coordinates is −qΔx/2,−qΔx2/12,−qΔx/2,qΔx2/12T, Here, P and q are respective the concentrated force and the distributed force applied locally on the element.Figure 5
The concentrated force.Figure 6
The distributed force.The local stiffness matrix and mass matrix can be written as(14)ki=2EiIiΔx363Δx−63Δx3Δxi2Δx2−3ΔxΔx2−6−3Δx6−3Δx3Δxi2Δx2−3Δx2Δx2,(15)mi=ρiAiΔx42015622Δx54−13Δx22Δx4Δx213Δx−3Δx25413Δx156−22Δx−13Δx−3Δx2−22Δx4Δx2.Besides, the damping matrix is constructed in a similar way as the mass matrix in (17):(16)ci=ciAiΔx42015622Δx54−13Δx22Δx4Δx213Δx−3Δx25413Δx156−22Δx−13Δx−3Δx2−22Δx4Δx2.The exciting force can be written accordingly as follows:(17)Fi=−qΔx2,−qΔx212,−qΔx2,qΔx212Tfor a distributed  load,(18)Fi=−PbΔx,−Pab2Δx,−PaΔx,Pab2ΔxTfor a concentrated load.The global stiffness matrix, the global mass matrix, and the global force vector could be obtained by assembling all the local stiffness matrices (15), the local mass matrices (16), and the local force vectors [92].


## 5. Computational Complexity

The following equation can be constructed for the dynamic response of the beam bridge:(19)Mtw¨+Ctw˙+Ktw=Ft.For the discretization schema of the beam bridge as lumped masses, the stiffness matrix is obtained by inverting the flexibility matrix. For the discretization schema of the beam bridge as finite elements, the stiffness matrix is obtained by assembling the local stiffness matrices of elements. The coefficients using the flexibility method and the stiffness method are given in Table1. Here, dij is the deflection at lumped mass j when the unit force is applied to the beam at lumped mass i.Table 1
Two different discretization schemas of beam.
Discretization of the beam as lumped masses (the flexibility method)Discretization of the beam as finite elements (the stiffness method)Ktn×n=Fd−1t=d11d12…d1nd21d22…d2n⋮⋮⋱⋮d11d11⋯d11−1Kt2n×2n=ki,i=1,2,…,nki=2EiIi/Δx363Δx−63Δx3Δxi2Δx2−3ΔxΔx2−6−3Δx6−3Δx3Δxi2Δx2−3Δx2Δx2Mtn×n=ρAΔx/EI10…001…0⋮⋮⋱⋮00⋯1Mt2n×2n=mi,i=1,2,...,nmi=ρiAiΔx/42015622Δx54−13Δx22Δx4Δx213Δx−3Δx25413Δx156−22Δx−13Δx−3Δx2−22Δx4Δx2Ctn×n=cAΔx/EI10…001…0⋮⋮⋱⋮00⋯1Ct2n×2n=ci,i=1,2,…,nci=ciAiΔx/42015622Δx54−13Δx22Δx4Δx213Δx−3Δx25413Δx156−22Δx−13Δx−3Δx2−22Δx4Δx2Fdt=qxi,tΔx/EI11⋮1n×1Ft2n×1=Fi,i=1,2,…,nFi=−qΔx/2,−qΔx2/12,−qΔx/2,qΔx2/12TorFi=−Pb/Δx,−Pab2/Δx,−Pa/Δx,Pab2/ΔxTEquation (17) could be solved using classical numerical methods, such as the central finite difference method, Houbolt method, Wilson theta method, and Newmark’s Beta method [93]. In this study, the Newmark method [94] is implemented. The procedures for solving (20) are given as follows.Step 1: compute coefficient matrix,Mt,Ct,andKt using the lumped masses and the finite elements, respectivelyStep 2: compute the excitation force vectorFtStep 3: solve equation (19) using Newmark’s methodThe size of the coefficient matrixM,C,andK for the discretization of the beam as lumped masses is n×n, whereas the size of the coefficient matrix M,C,andK for the discretization of the beam as finite elements is 2n×2n. Here, n is the number of lamped masses or finite elements after discretization. The matrices involved in the computation of the lumped masses are four times the size of the matrices involved in the computation of the finite elements.Neglecting the initialization step described in Table2, the computational operations (addition/subtraction and multiplication/division) for the algorithm are approximately 2n3/3 and 16n3/3, respectively. Therefore, the difference between the flexibility method and the stiffness method is due to the size of their coefficient matrices, which mainly account for the number of operations of the flexibility method and stiffness method. Theoretically, the efficiency of discretizing the beam as lumped masses in the flexibility method is roughly eight times more efficient than discretizing the beam as finite elements using the stiffness method.Table 2
Computational schema of Newmark’s method using (9).
AlgorithmDiscretization of the beam as lumped masses (the flexibility method)Discretization of the beam as finite elements (the stiffness method)1.0 initialization 1.1 u¨0=p0−cu˙0−ku0/m 2.2 set Δt 3.3 calculate a1,a2,a3 4.4 calculate k^O (1)O (1)2.0 loop for i = 1 to Tt5.5 calculatep^i+1OnO2n6.6 calculateui+12n3/3+On216n3/3+O4n27.7 calculateu^i+1OnO2n8.8 calculateu¨i+1OnO2nEnd


## 6. Performance Comparison

A numerical case study is provided in this section to illustrate and compare the computational accuracy and efficiency of two different discretization schemas.The case study involves a forced vibration of the beam bridge without the damping term. The equation of motion is shown in (19). The beam bridge is excited by an external force which varies in time and space; the coefficients of the term ∂2w/∂t2 and ∂4w/∂x4 are constant. The parameters of the bridge are adopted from Law et al. [77]; the total length of the beam bridge: L=30m, EI=2.5×1010N⋅m2, and ρA=5.0×103kg/m:(20)ρA∂2w∂t2+EI∂4w∂x4=1102ρA+EIt2πL2sinπLx,(21)w0,t=0,wL,t=0,w″0,t=0,w″L,t=0.The boundary conditions and the initial conditions are given in equations (21) and (22):(22)wx,0=0,w˙x,0=0.The theoretical solution for this beam bridge is provided in (23), which can be verified by directly substituting (23) into (19):(23)w=110t2sinπLx.Dynamic displacement for the first 1 second is computed at three different time steps: 0.01 s, 0.005 s, and 0.001 s and plotted in Figures4 and 5. Due to the symmetry of the geometry and the boundary condition, only two locations on the beam bridge are computed and compared: 0.3L and 0.5L. The relative error εx is defined as(24)εx=abscalculated value–theoretical valuetheoretical value×100%,where the subscript x of εstands for the method,f stands for the flexibility method, and s stands for the stiffness method.Figures7 and 8 show the relative error of both methods at 0.3L and 0.5L. The horizontal axis of both figures is computation time duration from the beginning into the simulation. Newmark’s method is a stable computational method when the time step is small enough, say, smaller than 10% of the overall time duration, [94]. Therefore, the cumulated error would not increase without bound as the computation time duration increases.Figure 7
The relative error at 0.3L.Figure 8
The relative error at 0.5L.According to Figures7 and 8, the relative errors of both methods are approximately at the same level for the same number of elements and the same time step, though the accuracy of the flexibility method is slightly better than that of the stiffness method. The variations of flexibility methods for a different number of elements are very small when the time step is the same.Increasing the number of elements or nodes not only slightly improves the accuracy of the flexibility method but also improves considerably the accuracy of the stiffness method.The time step plays a significant role in improving accuracy. Decreasing the time step improves significantly the accuracy for both methods at the cost of increased computation time. When the time step is small and the element number is over 50, the increase of element number does not improve the accuracy considerably for both methods. When the time step is small enough, even a small number of elements (e.g., 20 elements or nodes for a beam bridge of the length of 30 m) could achieve satisfactory accuracy.Figure9 plots the computation time of two discretization schemas versus the number of elements or nodes. It is evident that the computation time for the flexible method increases slightly with the increase in the number of elements, while the computation time for the stiffness method increases significantly with the increase in the number of elements.Figure 9
Efficiency comparison between the flexibility method and the stiffness method.Figure10 shows the time ratio, defined as the ratio between the time used by the flexibility method over the stiffness method, versus the number of elements or nodes. As the number of elements increases from 10 to 100, the time ratio increases from roughly 1.8, 2.2, and 3.2 for time steps 0.01 s, 0.005 s, and 0.001 s, respectively, to 4.9, 5.6, and 5.7. This is in general consistent with and slightly lower than the theoretical prediction on the computational complexity of the two methods, which suggests 8.0 computational efficiency differences between the two methods. The reason for such discrepancy between theoretical prediction and realization could be caused by coefficients embedded in computational complexity analysis in Table 2 and the actual computer used to execute the computation.Figure 10
The time ratio between the stiffness method over the flexibility method.


## 7. Vehicle Model

Figure11 shows a vehicle traversing a single span simply supported beam bridge. A four-degree-of-freedom (4-DOF) half-car model is utilized for analysis in this paper. Taking advantage of Lagrange formulation, the equation of motion of the vehicle can be derived as follows.Figure 11
A 4-DOF model of the vehicle-bridge interaction system.The kinetic energy equation of the vehicle model is(25)K=12mvy˙v2+12mfy˙12+12mry˙r2+12Ivθ˙2.The potential energy equation of the vehicle model is(26)V=12ktfyf−ybf2+12ktryr−ybr2+12ksfyv−yf+a1θ2+12ksryv−yr−a1θ2.The energy dissipation equation is(27)D=12ctfy˙f−y˙bf2+12ctry˙r−y˙br2+12csfy˙v−y˙f+a1θ˙2+12csry˙v−y˙f−a1θ˙2.Substituting into Lagrange formation.(28)ddt∂K∂q˙i−∂K∂qi+∂d∂q˙i+∂V∂qi=fii=1,2,3,4,equation of motion of the vehicle model is obtained, presenting in the matrix form,(29)Mvy¨+Cvy˙+Kvy=Fv,where(30)y=yvθyfyr,(31)Mv=mv0000Iv0000mf0000mr,(32)Cv=csf+csr−csfa1+csra2−csf−csr−csfa1+csra2csfa12+csra22csfa1−csra2−csfcsfa1csf+ctf0−csr−csra20csr+ctr,(33)Kv=ksf+ksr−ksfa1+ksra2−ksf−ksr−ksfa1+ksra2ksfa12+ksra22ksfa1−ksra2−ksfksfa1ksf+ktf0−ksr−ksra20ksr+ktr,(34)Fv=00ktfybf+ctfy˙bfktrybr+ctry˙br,(35)ybft=wxft,t+rxft,t,(36)ybrt=wxrt,t+rxrt,t,Here,ybft and ybrt are, respectively, the absolute vertical displacements of the front wheel and the rear wheel evaluated at the contact points along the route of the vehicle, which are determined based on pavement roughness and deflection of the bridge. In this study, it is assumed that the vehicle and bridge always keep full contact during the entire interaction process, which is reflected through equations (35) and (36). In these two equations, w(x, t) and r(x, t) stand for bridge deflection and pavement roughness, respectively. The values of ybft and ybrt are dependent on the location on the bridge x during the motion of the vehicle. For a moving vehicle at a constant speed, x=vt; therefore, ybft and ybrt are also functions of time.


## 8. Bridge Model

The bridge considered herein is a uniform cross-section beam. For varied cross-section beams, the same method is applied. The beam is discretized intomn lumped mass, as shown in Figure 3. In this study, only the vertical displacement of each lumped mass is considered. In other words, the total number of degree-of-freedom of the bridge is mn, excluding boundary conditions.The mass matrixMb for the bridge model can be conveniently obtained as a diagonal matrix. For the stiffness matrix Kb, the flexibility method is implemented. Taking advantage of knowledge from mechanics, the vertical deflection at each of the aforementioned lumped masses is calculated, resulting in the flexibility matrix. Taking the inverse of the flexibility matrix, the stiffness matrix Kb is obtained.In establishing the damping matrix, the classical Rayleigh damping for structure dynamics is adopted. The damping matrix isCb, which is a combination of both mass matrix Mb and stiffness matrix Kb:(37)Cb=αMb+βKb,where α=2ξω1ω2/ω1+ω2, β=2ξ/ω1+ω2, ξ is the damping ratio and ξ=0.02 is used in this study, and ω1 and ω2 are the first and second fundamental frequency, respectively.The equation of motion of the beam can be derived by Newton’s second law and reads(38)Mbw¨+Cbw˙+Kbw=Fbt,where Mb is the mass matrix of the discretized system, which is a mn×mn diagonal matrix, Cb is damping coefficient matrix, Kb is the stiffness matrix ybft and ybrt, Fbt is excitation force, concentrated on each lumped mass, which is induced by the motion of the vehicle, and w is the deflection vector at each of the lumped mass of the beam, the size of which is mn×1.The excitation forceFbt applied on the bridge model at each of the vehicle wheels is consisted of static and dynamic interaction load, which can be written as(39)Fbt=00ktfyf−ybft+ctfy˙f−y˙bft−mf+a2mvagktryr−ybrt+ctry˙r−y˙brt−mr+a1mvag.


## 9. Pavement Roughness Model

Pavement roughness is an important parameter affecting the dynamic response of the vehicle-bridge interaction system [95–97]. When the wheels are separated from the road surface [98], the dynamic response of the vehicle-bridge interaction system is dramatically different. While many studies ignore pavement roughness in analyzing vehicle-bridge interaction systems due to its complexity, we do take into account pavement roughness in this study.Dodds and Robson [99] proposed a method to generate pavement roughness by inverse Fast Fourier Transformation (FFT) based on displacement power spectral density (PSD) function. Pavement roughness in the spatial domain can be determined by(40)rx=∑i=1m4Gdni/Δmcos2πnix+θi.Here,rx is the pavement surface profile function, x=vt is the location indicator, which is time dependent, θi is the random phase angle uniformly distributed from 0 to 2π, ni is the wavenumber (cycle/m), ni=i/Δm, φni is the PSD function for the pavement surface profile (m3/cycle/m), Δ is the interval of successive points of the pavement surface profile, m is the total number of sample points of the surface profile, and Gdni is the PSD function.In this paper, the PSD function recommended by ISO-8608 [100] is employed:(41)Gdni=Gdn0nin0−wnl<ni<nu,where n0 is the reference spatial frequency (= 0.1 cycles/m), Gdn0 is the roughness coefficient (m3/cycle), nuandnl are upper and lower cut-off frequencies, w is the exponent of the fitted PSD, and w=2 is adopted according to ISO-8608 [100].Roads are classified into eight categories from A to H (very good to very poor) in ISO-8608 [100] according to pavement roughness, which is represented by Gdn0. Of the eight categories, the first four categories are used for describing pavement road surface while the last four are used for unpaved road surface [101]. The first three categories of pavement roughness shown in Table 3 are examined in this study.Table 3
ISO-8608 road roughness classification [100].
Road classGdn0×10−6mLower limitGeometric meanUpper limitA—1632B3264128C128256512


## 10. Equation of Motion of Vehicle-Bridge Interaction System

For a coupled vehicle-bridge interaction system, the vehicle and the bridge models are coupled by the relative displacement. To make sure that the vehicle is coupled with the corresponding lumped mass on the bridge, there exists an implicit bound,dx=vdt, in the discretization of the beam if velocity is constant though not mandatory. In order to obtain the dynamic response of the vehicle-bridge interaction system, the vehicle model and bridge model should be efficiently coupled together through contact points between the wheel and bridge pavement surface.The vehicle model needs to be transformed using more convenient variables before coupling with the bridge model. The dynamic motion of the vehicle is described by new variables, as shown in Figure12. The relation between the variables is as follows:(42)θ=u3−u4a,yv=a1u4+a2u3a,yf=u1,yr=u2,where uii=1,2,3,4 are the absolute displacement of the bodies and a1anda2 are shown in Figure 1, and a=a1+a2.Figure 12
A 4-DOF model of the vehicle-bridge interaction system with different representations.The conversion is conducted through the following transformation matrix:(43)yvθyfyr=a2aa1a00a−1a0000100001u3u4u1u2.After transformation, the equation of motion is written as(44)Mvtu¨+Cvtu˙+Kvtu=Fv,where(45)u=u3u4u1u2,Mvt=a2mvaa1mva00Iva−Iva0000mf0000mr,Cvt=csfcsr−csf−csrcsfa1−csra2−csfa1csra2−csf0csf+ctf00−csr0csr+ctr,Kvt=ksfksr−ksf−ksrksfa1−ksra2−ksfa1ksra2−ksf0ksf+ktf00−ksr0ksr+ktr.By introducing relative displacement variableszii=1,2,3,4 between the bodies and the corresponding bridge contact points, the equation of relative motion is acquired. The relative displacement variables are defined as(46)z3z4z1z2=u3u4u1u2−ybftybrtybftybrt,where zii=1,2,3,4 is the relative displacement with respect to the absolute displacement of the two contact points ybft and ybrt; uii=1,2,3,4 are the absolute displacement of the bodies.Substituting (46) into (44), the relative equation of motion is written as(47)Mvzz¨+Cvzz˙+Kvzz=Fvzt,where z=z3z4z1z2T, Mvz=Mvt, Cvz=Cvt, and Kvz=Kvt,(48)Fvzt=−a2mvay¨bft+a1mvay¨brtIvay¨bft−Ivay¨brtmfy¨bftmry¨brt.Applying (35) and (36), (49) is equivalent with(49)Fvzt=−a2mvaw¨xft,t+a1mvaw¨xrt,tIvaw¨xft,t−Ivaw¨xrt,tmfw¨xft,tmrw¨xrt,t−a2mvar¨xft,t+a1mvar¨xrt,tIvar¨xft,t−Ivar¨xrt,tmfr¨xft,tmrr¨xrt,t.The first term is incorporated into the left-hand side of the coupled system, leaving the second term on the right-hand side. Define the second term asFvzr, which is written as(50)Fvzrt=−a2mvar¨xft,t+a1mvar¨xrt,tIvar¨xft,t−Ivar¨xrt,tmfr¨xft,tmrr¨xrt,t.Similarly, when the bridge model equation of motion is described using the relative displacement variables, only the excitation forceFbt is updated as(51)Fbt=00ktfz1+ctfz˙1−mf+a2mvagktrz2+ctrz˙2−mr+a1mvag,where w¨xft,t and w¨xrt,t are acceleration at the lumped mass of the beam at the front wheel and the rear wheel, respectively, locating at xft and xrt; r¨xrt,t and r¨xrt,t are second derivatives, with respect to time t, of the road surface profile at the contacting points at xft and xrt, which is derived according to (32), as(52)r¨x=−2πv2∑i=1m4GdniΔmni2cos2πnix+θi.Combining vehicle model in (47) and bridge model in (38), the equation of motion for the coupled vehicle-bridge interaction system can be concisely written as(53)MvzM4mn0Mbz¨w¨+Cvz0Cmn4Cbz˙w˙+Kvz0Kmn4Kbzw=FvzrtFbt,where M4mn is a sparse matrix of size 4×mn with six nonzeros elements. The two nonzero elements of the first row are a2mv/a and a1mv/a, corresponding to wxft,t and wxrt,t, respectively; the two nonzero elements of the second row are Iv/a and −Iv/a, pairing with wxft,t and wxrt,t separately; each nonzero entry for the third and fourth row are mf and mr, matching with wxft,t and wxrt,t, respectively. Cmn4 is a sparse matrix of size mn×4 with two nonzero entries. The two nonzero entries are ctf and ctr, corresponding to wxft,t and wxrt,t, in column 3 and column 4, respectively. Similarly, the two nonzero entries for Kmn4 are −ktf and −ktr, occupying the same location as in Cmn4.Using Newmark’s method, (53) can be solved and time responses of both the bridge model and vehicle model are obtained. With the known relative displacement from (53), the axle loads at both of the wheels are hence achieved through (51). Given the parameters of the models, Mb, Cb, and Kbfor the bridge model, Mvz, Cvz. and Kvz for the vehicle model, the speed of the vehicle, the pavement surface roughness, and the dynamic responses of vehicle-bridge interaction are readily solvable.


## 11. Case Study

In this section, we use a single-span simply supported beam bridge for a case study. The fundamental frequency of the 30 m single-span bridge is 3.78 Hz. The parameters of the vehicle-bridge interaction system are presented in Table4. Parameters of the 4-DOF vehicle model are adopted from [102], which are measurements from the actual vehicle and have been utilized by [77].Table 4
Parameters for the vehicle-bridge interaction system.
Bridge parametersVehicle parametersL=30mmv=17735kgmf=1500kgmr=1000kgEI=2.5×1010N⋅m2Iv=1.47×105kg⋅m2csf=3.00×104NM/Scsr=4.00×104N/m/sρA=5.0×103kg/ma1=2.2mctf=3.90×103N/m/sctr=4.30×103N/m/sξ=0.02a2=1.1mksf=2.47×106N/mksr=4.23×106N/mktf=3.74×106N/mktr=4.60×106N/m


### 11.1. Influence of Number of Elements on Accuracy

In this section, pavement roughness is not considered. In order to select a proper number of elements to represent the bridge to achieve an acceptable level of accuracy, the displacement and the acceleration of the bridge at typical locations of the bridge are calculated with different numbers of elements, while the vehicle velocity remains a constant. It can be seen from Figures13 and 14 that when the element number is above 300, the influence on it is insignificant.Figure 13
Bridge displacementx bridge location (m) and y displacement (m). (a) Front wheel at L/3. (b) Front wheel at L/2. (c) Front wheel at 2L/3.
(a)(b)(c)Figure 14
Bridge acceleration historyx time (s) and y acceleration (m2/s)). (a) At L/3. (b) At L/2. (c) At 2L/3.
(a)(b)(c)


### 11.2. Dynamic Force

In this section, pavement roughness is not considered. The dynamic responses of the front wheel and the rear wheel described by relative displacement with respect to the road surface are shown in Figure15, which incorporate responses at different velocities. It can be observed that when speed increases, the magnitude of relative displacement goes up accordingly. Displacement of the middle span of the bridge when the vehicle traverses the bridge at a constant speed is plotted in Figure 16. The maximum displacement at the middle span of the bridge does not vary significantly. Displacement of the bridge when the front wheel is at L/2 is shown in Figure 17. It shows that, at an increasing vehicle velocity, the displacement of the bridge does not necessarily increase accordingly. The interaction force between the vehicle wheel and the bridge is presented in Figure 18. It is observed that peak interaction force grows when vehicle speed increases.Figure 15
The relative displacement of the front wheel and the rear wheel.Figure 16
Dynamic response of bridge at the middle span.Figure 17
Displacement of the bridge with different vehicle velocities.Figure 18
Moving force at different vehicle velocities.


### 11.3. Dynamic Load Coefficient (DLC)

In this section, dynamic load coefficient (DLC) [101] is used to analyze the influence of vehicle speed and pavement roughness. DLC is defined as the ratio between maximum dynamic axle loads and static vehicle load (Fstatic) as follows:(54)DLC=FbmaxFstatic.For pavement surface roughness levels A, B, and C, the simulated moving force of the front and the rear wheel is shown in Figure19, and the DLC is shown in Figure 20. DLC is usually in the range of [1.0, 2.0], as shown in Figures 19 and 20. It is noted that, for the dynamic load, coefficient of the front wheel is commonly higher than that of the rear wheel.Figure 19
Moving force at different velocities and different pavement roughness.
(a)(b)(c)(d)(e)(f)(g)(h)Figure 20
Peak dynamic load coefficient of the front wheel and rear wheel.
(a)(b)(c)(d)(e)For a smooth pavement, the DLC is in the magnitude of 1.01. Clearly, pavement roughness has a significant influence on DLC. At the same velocity, moving force and DLC generally grow with pavement roughness, while for the same pavement roughness, moving force and DLC generally grow with vehicle velocity.


### 11.4. DLC Standard

Different countries have adopted different criteria for allowable DLC [101]. American Association of Highways and Transportation Officials [103] technical specification uses 0.75 for deck joints, 0.15 for fatigue and fracture limit state, and 0.33 for all other limit states for all other components except joints. Australian technical specification [104] recommends a factor of 0.4 based on static wheel and axle load. Switzerland technical specification [105] adopts a factor of 0.8 for the bridge with a fundamental frequency between 2 and 4 Hz and 0.4 for frequency above 5 Hz. Chinese technical specification [85] recommends the following impact factor μ (i.e., allowable DLC):(55)μ=0.05,f<1.5Hz,0.1767lnf−0.0157,1.5Hz≤f≤14Hz,0.45,f>14Hz.Here,f is the fundamental frequency of the bridge. For pavement roughness in classes A and B, this study conforms with specifications of ASSHTO, Australia standard [104] and Switzerland standard [105]. However, the allowable DLC recommended by Chinese technical specification [85] is a little smaller for a bridge with a fundamental frequency of around 4 Hz. For a 30m long single-span bridge, Table 5 shows the recommended DLC based on technical specifications from different countries.Table 5
Recommended DLC based on technical specifications from different countries.
CodeSpecifications for considering dynamic effectDLCASSHTO [103]0.75 for deck joints, 0.15 for fatigue and fracture limit state, and 0.33 for all other limit states, for all other components except joints1.33Chinese code [85]μ=0.05,f<1.5Hz,0.1767lnf−0.0157,1.5Hz≤f≤14Hz,0.45,f>14Hz1.22Australian standard [104]Recommends a factor of 0.4 based on static wheel and axle load1.4Switzerland standard [105]A factor of 0.8 is adopted for a bridge with a fundamental frequency between 2 and 4 Hz and 0.4 for frequency above 5 Hz1.8EN 1991-2 [106]For carefully maintained trackΦ3=1.44/LΦ−0.2+0.82,1.0≤Φ3≤1.67 for a track with standard maintenance Φ3=2.16/LΦ−0.2+0.731.0≤Φ3≤2.0 is determinant length depending on the structural element type1.14 (standard maintenance


## 11.1. Influence of Number of Elements on Accuracy

In this section, pavement roughness is not considered. In order to select a proper number of elements to represent the bridge to achieve an acceptable level of accuracy, the displacement and the acceleration of the bridge at typical locations of the bridge are calculated with different numbers of elements, while the vehicle velocity remains a constant. It can be seen from Figures13 and 14 that when the element number is above 300, the influence on it is insignificant.Figure 13
Bridge displacementx bridge location (m) and y displacement (m). (a) Front wheel at L/3. (b) Front wheel at L/2. (c) Front wheel at 2L/3.
(a)(b)(c)Figure 14
Bridge acceleration historyx time (s) and y acceleration (m2/s)). (a) At L/3. (b) At L/2. (c) At 2L/3.
(a)(b)(c)


## 11.2. Dynamic Force

In this section, pavement roughness is not considered. The dynamic responses of the front wheel and the rear wheel described by relative displacement with respect to the road surface are shown in Figure15, which incorporate responses at different velocities. It can be observed that when speed increases, the magnitude of relative displacement goes up accordingly. Displacement of the middle span of the bridge when the vehicle traverses the bridge at a constant speed is plotted in Figure 16. The maximum displacement at the middle span of the bridge does not vary significantly. Displacement of the bridge when the front wheel is at L/2 is shown in Figure 17. It shows that, at an increasing vehicle velocity, the displacement of the bridge does not necessarily increase accordingly. The interaction force between the vehicle wheel and the bridge is presented in Figure 18. It is observed that peak interaction force grows when vehicle speed increases.Figure 15
The relative displacement of the front wheel and the rear wheel.Figure 16
Dynamic response of bridge at the middle span.Figure 17
Displacement of the bridge with different vehicle velocities.Figure 18
Moving force at different vehicle velocities.


## 11.3. Dynamic Load Coefficient (DLC)

In this section, dynamic load coefficient (DLC) [101] is used to analyze the influence of vehicle speed and pavement roughness. DLC is defined as the ratio between maximum dynamic axle loads and static vehicle load (Fstatic) as follows:(54)DLC=FbmaxFstatic.For pavement surface roughness levels A, B, and C, the simulated moving force of the front and the rear wheel is shown in Figure19, and the DLC is shown in Figure 20. DLC is usually in the range of [1.0, 2.0], as shown in Figures 19 and 20. It is noted that, for the dynamic load, coefficient of the front wheel is commonly higher than that of the rear wheel.Figure 19
Moving force at different velocities and different pavement roughness.
(a)(b)(c)(d)(e)(f)(g)(h)Figure 20
Peak dynamic load coefficient of the front wheel and rear wheel.
(a)(b)(c)(d)(e)For a smooth pavement, the DLC is in the magnitude of 1.01. Clearly, pavement roughness has a significant influence on DLC. At the same velocity, moving force and DLC generally grow with pavement roughness, while for the same pavement roughness, moving force and DLC generally grow with vehicle velocity.


## 11.4. DLC Standard

Different countries have adopted different criteria for allowable DLC [101]. American Association of Highways and Transportation Officials [103] technical specification uses 0.75 for deck joints, 0.15 for fatigue and fracture limit state, and 0.33 for all other limit states for all other components except joints. Australian technical specification [104] recommends a factor of 0.4 based on static wheel and axle load. Switzerland technical specification [105] adopts a factor of 0.8 for the bridge with a fundamental frequency between 2 and 4 Hz and 0.4 for frequency above 5 Hz. Chinese technical specification [85] recommends the following impact factor μ (i.e., allowable DLC):(55)μ=0.05,f<1.5Hz,0.1767lnf−0.0157,1.5Hz≤f≤14Hz,0.45,f>14Hz.Here,f is the fundamental frequency of the bridge. For pavement roughness in classes A and B, this study conforms with specifications of ASSHTO, Australia standard [104] and Switzerland standard [105]. However, the allowable DLC recommended by Chinese technical specification [85] is a little smaller for a bridge with a fundamental frequency of around 4 Hz. For a 30m long single-span bridge, Table 5 shows the recommended DLC based on technical specifications from different countries.Table 5
Recommended DLC based on technical specifications from different countries.
CodeSpecifications for considering dynamic effectDLCASSHTO [103]0.75 for deck joints, 0.15 for fatigue and fracture limit state, and 0.33 for all other limit states, for all other components except joints1.33Chinese code [85]μ=0.05,f<1.5Hz,0.1767lnf−0.0157,1.5Hz≤f≤14Hz,0.45,f>14Hz1.22Australian standard [104]Recommends a factor of 0.4 based on static wheel and axle load1.4Switzerland standard [105]A factor of 0.8 is adopted for a bridge with a fundamental frequency between 2 and 4 Hz and 0.4 for frequency above 5 Hz1.8EN 1991-2 [106]For carefully maintained trackΦ3=1.44/LΦ−0.2+0.82,1.0≤Φ3≤1.67 for a track with standard maintenance Φ3=2.16/LΦ−0.2+0.731.0≤Φ3≤2.0 is determinant length depending on the structural element type1.14 (standard maintenance


## 12. Concluding Remarks

The accuracy and efficiency of discretizing the beam in space as lumped masses using the flexibility method and as finite elements using the stiffness method are studied and compared against the analytical solutions of a case study. In general, both methods result in a similar level of accuracy for the same number of elements or nodes as well as the same time step, though the former is slightly better than that of the latter.The time step plays a significant role in improving accuracy. Decreasing the time step improves significantly the accuracy for both methods at the cost of increased computation time.Increasing the number of elements or nodes not only slightly improves the accuracy of the flexibility method but also improves considerably the accuracy of the stiffness method.When the time step is small and the element number is over 50, the increase of the element number does not improve the accuracy considerably for both methods. When the time step is small enough, even a small number of elements could achieve satisfactory accuracy.In terms of computation efficiency, discretizing the beam as lumped masses will be in theory eight times or in reality two to six times faster than discretizing the beam as finite elements. The former is, even more, faster than the latter when the number of elements or nodes increases. This conclusion paves the way for developing a high efficiency algorithm for the dynamics of vehicle-bridge interaction.This high efficiency algorithm and corresponding discretization schema are applied to study vehicle-bridge interaction systems, in which a 4-DOFs vehicle model, a simply supported single-span bridge, and ISO pavement roughness model are used. The coupled equation of the vehicle-bridge interaction system is solved directly.When pavement is smooth, peak interaction force of the front wheel and the rear wheel in general increases when vehicle speed increases, and the DLC is typically below 1.01 for vehicle speed in the range of [0, 30 m/s].Pavement roughness significantly influences DLC. Therefore, preventative maintenance or timely maintenance of the pavement surface on a bridge to reduce pavement roughness while the bridge is still in a good condition is of significant importance for the bridge’s longevity and life cycle cost benefit.For class A and B level pavement roughness, DLC is simulated within 2.0, compatible with specifications of ASSHTO, Australian standard [104], and Switzerland standard [105]. However, the Chinese code [85] underestimates the DLC for a bridge with a fundamental frequency of around 4 Hz.



---

*Source: 1964200-2021-08-11.xml*