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.

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 [132], bridge dynamics and vibration control [3350], bridge structural health monitoring, bridge ambient testing, bridge damage identification [5164], and high-speed railway [29, 6568].

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, 6976].

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, 7684]. 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 load q(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 length L subjected to an arbitrarily distributed force , as shown in Figure 1. The parameters of the beam are which denotes Young’s modulus, which denotes the cross-section moment of inertia, which denotes the density, and which denotes the area of cross section. It is assumed that the beam bridge is uniform in cross section.

Taking out an infinitesimal segment from the beam, with which the force applied on this segment is shown in Figure 2. Here, denotes the moment applied on left cross section, denotes the shear force applied on left cross section, denotes the moment applied on right cross section, denotes the shear force applied on right cross section, denotes the distributed force applied, denotes the mass of the infinitesimal segment volume, denotes the acceleration of the infinitesimal segment volume, , where is the deflection of the beam and is time variable, is the unit damping force caused by the ambiance and material viscosity, which is proportional to the velocity, , where is the damping coefficient of the unit length.

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:or equivalently

Take moment with respect to point O. The moment balance equation can be obtained as follows:

Substituting and into (2) leads to

Divide both sides by and take the limit of , obtaining

Substitution of into (1b) leads to

According to Euler–Bernoulli beam theory,

The equation of motion for the beam bridge is obtained as follows:where and could be functions of both position x and time t, that is, , . This generalization leads to

The coefficient of the third term in the continuous equation of motion (8) could be set to 1 by dividing into both sides of (8), resulting in

For a uniform beam with a constant , (7) can be further simplified aswhere , , and .

When neglecting the damping, i.e., , (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 , and is the length of the bridge. The coefficients of , , and 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 , as shown in Figure 3.

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]:where is the concentrated force acting on the beam, are defined the same way as before, a is the distance from the origin, and .

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 [84]. The mass matrix is obtained by constructing a diagonal matrix of size , with the element of value . The force vector applied on the lumped mass could be directly represented by .

4. Discretization of Beam as Finite Elements

The beam bridge could also be discretized as n finite elements, with 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 .

Multiplying on both sides of equation (1a), discretizing the system into

Equation (12) can be written in a matrix form as

Here, is equivalent to , is equivalent to , is equivalent to , is equivalent to , mass matrix of the discretized system, is the damping coefficient matrix, is the stiffness matrix, is the excitation force that is applied on each lumped mass, and is the deflection vector of the beam.

In Figure 5, 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 . For a distributed load, as shown in Figure 6, the equivalent node force vector in the local coordinates is , Here, and are respective the concentrated force and the distributed force applied locally on the element.

The local stiffness matrix and mass matrix can be written as

Besides, the damping matrix is constructed in a similar way as the mass matrix in (17):

The exciting force can be written accordingly as follows:

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:

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 Table 1. Here, is the deflection at lumped mass when the unit force is applied to the beam at lumped mass .

Equation (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, using the lumped masses and the finite elements, respectivelyStep 2: compute the excitation force vector Step 3: solve equation (19) using Newmark’s method

The size of the coefficient matrix for the discretization of the beam as lumped masses is , whereas the size of the coefficient matrix for the discretization of the beam as finite elements is 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 Table 2, the computational operations (addition/subtraction and multiplication/division) for the algorithm are approximately and , 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.

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 and are constant. The parameters of the bridge are adopted from Law et al. [77]; the total length of the beam bridge: , , and :

The boundary conditions and the initial conditions are given in equations (21) and (22):

The theoretical solution for this beam bridge is provided in (23), which can be verified by directly substituting (23) into (19):

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 Figures 4 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 is defined aswhere the subscript of stands for the flexibility method, and stands for the stiffness method.

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

According to Figures 7 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.

Figure 9 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 10 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.

7. Vehicle Model

Figure 11 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.

The kinetic energy equation of the vehicle model is

The potential energy equation of the vehicle model is

The energy dissipation equation is

Substituting into Lagrange formation.equation of motion of the vehicle model is obtained, presenting in the matrix form,where

Here, and 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, (x, t) and r(x, t) stand for bridge deflection and pavement roughness, respectively. The values of and are dependent on the location on the bridge during the motion of the vehicle. For a moving vehicle at a constant speed, ; therefore, and 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 into 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 , excluding boundary conditions.

The mass matrix for the bridge model can be conveniently obtained as a diagonal matrix. For the stiffness matrix , 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 is obtained.

In establishing the damping matrix, the classical Rayleigh damping for structure dynamics is adopted. The damping matrix is , which is a combination of both mass matrix and stiffness matrix :where , , is the damping ratio and is used in this study, and and are the first and second fundamental frequency, respectively.

The equation of motion of the beam can be derived by Newton’s second law and readswhere is the mass matrix of the discretized system, which is a diagonal matrix, is damping coefficient matrix, is the stiffness matrix and , is excitation force, concentrated on each lumped mass, which is induced by the motion of the vehicle, and is the deflection vector at each of the lumped mass of the beam, the size of which is .

The excitation force applied on the bridge model at each of the vehicle wheels is consisted of static and dynamic interaction load, which can be written as

9. Pavement Roughness Model

Pavement roughness is an important parameter affecting the dynamic response of the vehicle-bridge interaction system [9597]. 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

Here, is the pavement surface profile function, is the location indicator, which is time dependent, is the random phase angle uniformly distributed from 0 to , is the wavenumber (cycle/m), , is the PSD function for the pavement surface profile (), is the interval of successive points of the pavement surface profile, is the total number of sample points of the surface profile, and is the PSD function.

In this paper, the PSD function recommended by ISO-8608 [100] is employed:where is the reference spatial frequency (= 0.1 cycles/m), is the roughness coefficient (), are upper and lower cut-off frequencies, is the exponent of the fitted PSD, and 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 . 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.

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, , 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 Figure 12. The relation between the variables is as follows:where are the absolute displacement of the bodies and are shown in Figure 1, and .

The conversion is conducted through the following transformation matrix:

After transformation, the equation of motion is written aswhere

By introducing relative displacement variables between the bodies and the corresponding bridge contact points, the equation of relative motion is acquired. The relative displacement variables are defined aswhere is the relative displacement with respect to the absolute displacement of the two contact points and ; are the absolute displacement of the bodies.

Substituting (46) into (44), the relative equation of motion is written aswhere , , , and ,

Applying (35) and (36), (49) is equivalent with

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 as , which is written as

Similarly, when the bridge model equation of motion is described using the relative displacement variables, only the excitation force is updated aswhere and are acceleration at the lumped mass of the beam at the front wheel and the rear wheel, respectively, locating at and ; and are second derivatives, with respect to time , of the road surface profile at the contacting points at and , which is derived according to (32), as

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 aswhere is a sparse matrix of size with six nonzeros elements. The two nonzero elements of the first row are and , corresponding to and , respectively; the two nonzero elements of the second row are and , pairing with and separately; each nonzero entry for the third and fourth row are and , matching with and , respectively. is a sparse matrix of size with two nonzero entries. The two nonzero entries are and , corresponding to and , in column 3 and column 4, respectively. Similarly, the two nonzero entries for are and , occupying the same location as in .

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, , , and , , . and 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 Table 4. Parameters of the 4-DOF vehicle model are adopted from [102], which are measurements from the actual vehicle and have been utilized by [77].

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 Figures 13 and 14 that when the element number is above 300, the influence on it is insignificant.

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 Figure 15, 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 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.

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 () as follows:

For pavement surface roughness levels A, B, and C, the simulated moving force of the front and the rear wheel is shown in Figure 19, 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.

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):

Here, 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.

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.

Data Availability

The data, models, or code generated or used during the study are proprietary or confidential in nature and can be obtained from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This study was sponsored in part by the National Science Foundation, under Grant CMMI-0644552, by the National Natural Science Foundation of China, under grant U1134206, by Yunnan Department of Transportation (7921000079), and by Shangdong Expressway Co. Ltd, under Grant 2017-BLKY4, to which the authors are very grateful.