# Earthquake Response Spectra Analysis of Bridges considering Pounding at Bilateral Beam Ends Based on an Improved Precise Pounding Algorithm

Authors: Ruijie Zhang; Lei Yan; Kefeng Yue; Junhong Yin; Kang An
Journal: Shock and Vibration
 (2021)

---

## Abstract

Asynchronous vibration was generated between the main bridge and approach spans or abutments due to differences in stiffness and mass during an earthquake, thus further leading to pounding at the bilateral beam ends. By taking a T-shaped rigid frame bridge as an example, the bilateral pounding model was abstracted, and the earthquake response spectra considering pounding at the bilateral beam ends were studied, including the maximum displacement spectrum, the acceleration dynamic coefficient spectrum, the pounding force response spectrum, and the response spectrum for the number of pounding events. An improved precise pounding algorithm was proposed to solve the dynamic equation of the bilateral pounding model. This algorithm is based on the precise integration method for solving the second-order dynamic differential equation and reduces the order thereof by introducing a new velocity vector and uses the series method to find the nonhomogeneous term. The system matrix is simpler, and the inversion of the system matrix can be avoided. On this basis, a multipoint earthquake-induced pounding response spectrum program was developed. A total of 18 seismic waves from Class II sites were selected, and the response spectra of 18 waves were analyzed using this new program. Furthermore, the effects of structural stiffness, mass, stiffness of contact element, pounding recovery coefficient, and peak ground acceleration (PGA) on the earthquake response spectrum were studied. Through the analysis of earthquake response spectra and a parametric study, the phenomenon of earthquake-induced pounding of bridges was clarified to the benefit of the analysis and engineering control of earthquake-induced pounding of bridges.


---

## Body

## 1. Introduction

The pounding that may arise between adjacent structures causes structural damage during earthquakes. Some scholars have studied earthquake-induced pounding through the analysis of earthquake response spectra considering pounding between adjacent structures. Kawashima and Sato [1] studied the maximum displacement response spectrum of adjacent structures; Ruangrassameeand and Kawashima [2] presented the relative displacement response spectrum under the impact of pounding; Jankowski [3] studied the structural pounding response spectrum under excitation by the El-Centro wave; and Yaghmaei-Sabegh and Jalali-Milani [4] studied the pounding response spectrum under both near-field and far-field earthquakes. In recent years, Chinese scholars have studied the pounding response spectrum from different perspectives, and the pounding response spectrum based on the precise integration method has been studied [5, 6]. Furthermore, the acceleration dynamic coefficient spectra under different field wave excitations were examined [7].The single-point pounding based on two SDOF structures has been widely explored in previous studies, while multipoint pounding events are more commonplace in bridge structures. The pounding between the main bridge of common continuous beam, continuous rigid frame bridge, or T-shaped rigid frame bridge, and the approach bridge or abutment with their inherent differences in mass and stiffness can be produced due to asynchronous vibration under seismic excitation. The pounding between the main bridge and the structure on both sides will affect the vibration response of the main bridge structure, which in turn will affect the pounding force and the number of pounding events during an earthquake; therefore, a model containing three SDOF structures and two contact elements is abstracted to simulate pounding events of the main bridge with bilateral abutments or approach spans. At present, most of the research on pounding spectra focuses on the pounding force spectrum, dynamic coefficient spectrum, or displacement spectrum. The number of pounding events is also taken as an important evaluative indicator for earthquake-induced pounding analysis of engineering structures, based on which the concept of number of pounding events as a response spectrum was proposed in the present research.The time-history analysis of second-order dynamic differential equation with pounding force is the basis of reaction spectrum calculation. Traditional solutions to second-order dynamic differential equations include the explicit central difference method and the implicit Newmark-β method [8, 9]. These methods are widely used in commercial software packages but have only second-order numerical accuracy. The precise integration method proposed by Zhong [10, 11] reduces the order of the second-order differential equation to that of a first-order equation by introducing dual variables and uses the 2N class algorithm to obtain the high-precision solution of the exponential matrix. This method has absolute numerical stability, zero amplitude attenuation rate, zero period extension rate, and no transcendence. The precise integration method has been applied to many structural dynamics problems [12, 13], and the structures involved have also been extended from linear elastic structures to nonlinear structures [14, 15]. The earthquake-induced pounding phenomenon of a bridge structure is a nonlinear problem of state, which is more complicated than the general structural nonlinear state. Zhang et al. [16, 17] extended the precise integration method to solve the earthquake-induced pounding problems of adjacent structures. Through combining it with a contact element model, a precise pounding algorithm using the strategy of different lengths of time step was proposed for use in different structural states, and the precise pounding algorithm is applied to analyse multipoint pounding of elastic-plastic structures. For the pounding dynamic equation [18, 19], the dual variable is introduced to reduce the order of the equation, and the analytical method based on the linear variation of load assumption is adopted to find the nonhomogeneous term. After order reduction, the order of the system matrix is doubled, and the sparsity characteristic of the stiffness matrix is lost. In addition, it is necessary to invert the system matrix when solving the nonhomogeneous term, which is onerous in terms of computational burden, and may also give rise to the situation where the inverse matrix does not exist.To avoid inversion of the system matrix and accelerate the calculation process, researchers have improved the precise integration method [18–23]. Gu et al. [19] proposed an augmented dimensional precise integration method. By expanding the scale of the system matrix, the nonhomogeneous terms were transformed into homogeneous terms to avoid matrix inversion; however, this method requires more data storage and is computationally more onerous. The direct integration of nonhomogeneous terms can also avoid inversion of the system matrix. Lin et al. [20] proposed the method of Fourier expansion for nonhomogeneous terms. Wang and Au [21] developed a Gaussian quadrature method and piece wise interpolation polynomial method [22]. In addition, Simpson, Romberg, Cots, and other direct integration methods were used to find the nonhomogeneous terms [23, 24]. For order reduction in systems of second-order dynamic differential equations, others [25] have combined the Newmark method therein for order reduction and adopted series solution methods to deal with nonhomogeneous terms. However, the calculation accuracy of this method is related to the step size, so it loses the superiority of the precise integration method of “increasing dimension and order reduction.” Ding [14] reduced the order of the equation by introducing a new velocity vector, thus simplifying the calculation of the exponential matrix.In time-history analysis considering pounding, contact element was used to calculate the pounding force, and a variety of contact element models were developed [26–30]. Based on the bilateral pounding model proposed herein, the dynamic differential equation considering the bilateral pounding was given by combining the Hertz-damping contact element [29]. The dynamic differential equation was reduced by introducing a new velocity vector, and the nonhomogeneous term was found by the use of a series method, which simplifies the system matrix and avoids inversion of the system matrix. This improvement can not only avoid the difficulty of solving the system matrix in special cases, but also greatly improve the calculation speed, especially in the case of response spectrum analysis, which requires a lot of repeated time-history calculation. The improved precise pounding algorithm was verified, and then a program for calculating multipoint pounding response spectra was developed.The maximum displacement spectrum, acceleration dynamic coefficient spectrum, pounding force spectrum, and the number of pounding event spectrum of selected seismic waves were analyzed, and the influence of abutment stiffness, abutment mass, contact element stiffness, impact recovery coefficient, and peak ground acceleration on displacement, dynamic amplification coefficient, pounding force, and the number of pounding events were studied.


## 2. Definition of Pounding Spectrum

The pounding spectrum is a curve or surface relationship that considers the variation in earthquake-induced response with the natural vibration period of the SDOF involved in pounding. In previous studies [3–6], the pounding reaction involves structural displacement, a pounding force, or dynamic amplification coefficient, and the object of interest has two elastic SDOF with only one pounding point. In fact, the number of pounding points (or events) is also important; for a T-shaped rigid frame bridge, in general, both sides of the abutment have larger stiffness and smaller mass, while the middle pier has a smaller antipush stiffness and larger pier beam mass. Under longitudinal earthquake excitation, the main girder may collide with the abutment structure on both sides (Figure 1), and in this simplified model, m1, k1, and ξ1 are the simplified mass, stiffness, and damping ratio of the SDOF structural abutment on both sides; m2, k2, and ξ2 are the simplified mass, stiffness, and damping ratio of the SDOF structure of the pier and main girder; and kc, e, and gp are, respectively, the stiffness parameters, recovery coefficient, and clearance of the contact element.Figure 1
Bilateral pounding model.Therefore, the pounding displacement spectrum is depicted as representing the variation of the pounding particle displacement with the natural vibration period of the SDOF considering the specific structural parameters, pounding parameters, and ground motion input. The structural parameters include the mass, stiffness, and damping ratio of the SDOF, and the pounding parameters include the contact element model type, pounding stiffness, pounding recovery coefficient, and gap width. The dynamic acceleration amplification coefficient is the ratio of the peak acceleration to the peak seismic wave. The spectrum of the dynamic acceleration amplification coefficient,m2, represents the variation in the dynamic acceleration amplification coefficient of a particle with its natural vibration period. The pounding force spectrum represents the variation in the maximum pounding forces at two possible pounding points with the natural vibration period of the SDOF. The number of pounding event response spectrum represents the variation in the maximum number of pounding events arising from the presence of two possible pounding points with the natural vibration period of the SDOF.


## 3. Pounding Dynamic Differential Equation and Improved Precise Pounding Algorithm

### 3.1. Pounding Dynamic Differential Equation

The beam ends pounding model shown in Figure1 contains three SDOF particles, in which the mass, damping, and stiffness of the left and right particles are the same. The pounding force time-history between the middle particle and the left particle is assumed to be f1t, and that between the right particle and the middle particle is assumed to be f2t. The dynamic balance equation of the system is as follows:(1)MX¨+CX˙+KX+Ft=−Mgt,where M=m1000m2000m1, C=c1000c2000c1, K=k1000k2000k1, X=x1x2x3, and Ft=f1t−f1t+f2t−f2t.The structural damping coefficient is given by(2)ci=2ξimiki,i=1,2.The Hertz damping contact element model is used for the pounding between particles, and the viscous damper is added on the basis of the Hertz model to allow energy loss.The spring stiffness and damping coefficient are variables. In the contact stage, spring stiffnessβ and damping coefficient ηare related to the penetration displacement of the two particles, as given by the following equations:(3)βit=kcxi−xi+1−gp,(4)ηit=ζxi−xi+1−gp3/2,(5)ζ=3kc1−e24u˙i−u˙i+1,where kc, e, and gp are, respectively, the spring constant, recovery coefficient, and initial gap of the contact element. xi and xi+1 are the displacements of two adjacent particles at time t, and u˙iu˙i+1 are the velocities of two adjacent particles before contact.The pounding force on a contact element is obtained as follows:(6)fit=βitxi−xi+1−gp+ηitx˙i−x˙i+1.In general, the pounding force can be treated as the hysteresis force, that is, the displacement and velocity calculated in the previous time step are used to calculate the pounding force with equation (6), which is incorporated into the right-hand side of the equation as the known load; therefore, the calculated pounding force always lags by one time step. To overcome this problem, the pounding force column vector of the contact element in the model can be written as follows:(7)Ft=βtX+ηtX˙+γtgp,where βtX+γtgp is the elastic restoring force vector and ηtX˙ is the damping force vector. In detail, βt=β1t−β1t0−β1tβ1t+β2t−β2t0−β2tβ2t refers to the stiffness matrix of contact elements at time t. ηt=η1t−η1t0−η1tη1t+η2t−η2t0−η2tη2t refers to the damping coefficient matrix of contact elements at time t, and γt=−β1tβ1t+β2t−β2t is the column vector to calculate elastic restoring force related to the initial gap.The decomposed pounding force column vector is introduced into the pounding dynamic differential equation, which can be rewritten as follows:(8)MX¨+C+ηtX˙+K+βtX=−Mgt−γtgp.


### 3.2. Improved Precise Integration Method

To reduce the order of the dynamic balance equation of the system, we introduce vectorp. Let p=MX˙, then X˙=M−1p in the dynamic equilibrium equation, which can be written as follows:(9)p˙=−K+βtX−C+ηtM−1p−γtgp−Mgt.Assumingv=Xp, we have(10)v˙=Hv+ft,where the system matrix H=0M−1−K+βt−C+ηtM−1 and nonhomogeneous term ft=0−γtgp−Mgt.Assumming that the solution of the reduced first-order equation of state at timetk is vtk, the system matrix H is time-variant, so tk is denoted asHtk, and the general solution at time tk+1is as follows:(11)vtk+1=expHtk⋅Δtvtk+∫tktk+1expHtktk+1−τfτdτ.The exponential matrix of the first term on the right of the above equation is represented as follows:(12)T=expHtk⋅Δt=expHtk⋅Δt2N2N=expHtk⋅s2N=I+Ta2N.The matrixTa is expanded by Taylor series as follows:(13)Ta=∑l=1∞Htk⋅sll!.In equation (13), the truncation order of the series expansion is five, and N is 20, so a more accurate solution can be obtained. After matrix Ta is found, matrix T is found by the addition principle as applied to an exponential function. Specifically, the matrix is iterated N times the following equations:(14)Ta=2Ta+Ta×Ta,(15)T=I+Ta.The series method is used to solve the nonhomogeneous integral, supposes=τ−tk, then(16)∫tktk+1expHtktk+1−τfτdτ=∫0ΔtexpHtkΔt−sftk+sds=T∫0Δtexp−Htksftk+sds.Suppose the nonhomogeneous term varies linearly over timetiti+1, that is, ftk+s=ftk+Δf⋅s. Then, the solution of the integral term of equation (16) is as follows:(17)∫0Δt∑l=0∞−Htksll!ftk+Δf⋅sds=∑l=0∞−Htkll!ftk⋅Δtl+1l+1+Δf⋅Δtl+2l+2.Suppose the truncation order of the series expansionl in equation (17) is five, then, the precise recursion format of the pounding dynamic differential equation is as follows:(18)vtk+1=Tvtk+∑l=05−Htkll!ftk⋅Δtl+1l+1+Δf⋅Δtl+2l+2.The elements in the first half of the vectorvtk+1 constitute the displacement X at time tk+1, then the rest elements constitute the vector p. The velocity at time tk+1 can be calculated by X˙=M−1p, and substitution into the dynamic balance equation yields the acceleration at time tk+1.The above derivation process is derived from the dynamic balance equation of the full state. For problems considering structural nonlinearity, the secant stiffness matrix is time-varying, so it needs to be solved by iteration. In some cases, the secant stiffness will be very large or very small, resulting in significant numerical error in the calculation; therefore, the incremental dynamic balance equation is usually adopted after the structure is deemed elastic-plastic. For the incremental dynamic balance equation, the precise integration method is also applicable after slight modification. The recursive solution requires that the increment from the previous step is deferred, and the increment of the next step is deferred. The state of the structure at a given time must be obtained by incremental superposition. In the literature [17], a precise integration method for solving the incremental dynamic balance equation of structural nonlinear problems was examined: in the present work, the relevant results were cited along with the description of the program design.


### 3.3. Precise Pounding Algorithm

The pounding force column vector is introduced when the dynamic differential equation of the pounding is established. It should be noted that the state of the structure in the earthquake-induced response can be divided into separate, and contact, states. The current state of the structure can be obtained by judging the displacement difference between the ends of the two contact elements. In the contact state, because the stiffness and damping coefficients of the pounding spring are time-variant, the system matrix is also time-variant. Upon separation, the elements inβt, ηt, and γt are all zero, the system matrix of elements is the time-varying stiffness matrix related to the linear elastic stage or incremental dynamic balance equation of the tangent stiffness constant and plastic phase, the elements in the system matrix are constant, as are the changes in secant stiffness or tangent stiffness in the plastic region, meaning that the system matrix is still time-variant.In the process of the step-by-step detailed recursive solution of the pounding dynamic equation, the separation-contact state transition should also be solved. By dividing the step sizeΔt into 2N substeps, the precise integration method can provide a high-precision solution. In the state of separation, the integral step size Δt can be larger, which is conducive to reducing the computational burden. In the contact state, due to the short duration of the process itself, a smaller integral step size should be adopted. In the program design, the integral step Δt1 of the separation state is set to the sampling interval of seismic wave, and the integral step of the pounding stage is set to 1% to 10% of Δt1.In the recursion process used in the precise integration method, the pounding phenomenon is judged by the value ofxi−xi+1−gp. Assume that time tk is xi−xi+1−gp<0, and time tk+1 is xi−xi+1−gp≥0, indicating the exact time of pounding. The displacement, velocity, and acceleration at time tk and time tk+1 have been deduced by the precise integration method of separate states. According to the assumption that the acceleration changes linearly from tk to tk+1, the structural displacement at time tk+τ can be given by the following equation:(19)xitk−xi+1tk+x˙itk−x˙i+1tkτ+x˙itk−x˙i+1tkτ22+x¨itk+1−x˙itk−x¨i+1tk+1−x˙i+1tkΔtτ36−gp=0.The equation is a cubic equation containing one variable, and the solution in the interval0Δt can be derived by using the dichotomy method, and time tk+τ is the initial time of pounding. The structural state at the exact start time tk of pounding can be obtained by changing the step size to τ, and then the integral step size in the contact state is changed to Δt2, and the system matrix H is set at each step. The structural displacement, velocity, and acceleration under pounding action can be solved by the iterative method until xi−xi+1−gp≤0. At this point, the structure changes from a state of contact to one of separation. Similarly, equation (19) can be used to find the contact separation time. After entering the separation state, the integral step size is changed again toΔt1, and the recursive solution of the separation state is started. The above processes can be conducted alternately until the input seismic wave analysis is completed, and the time-history response of the structure is obtained.


### 3.4. Algorithm Verification

The bilateral pounding model is abstracted from a 2 × 85 m T-shaped rigid frame bridge. The elevation of the bridge is shown in Figure2. In the abstract model, the mass of the main bridge is concentrated at the centroid of the main girder on the pier top, m2=8.12×106 kg, the height of the main pier is 52 m, the calculated antipush stiffness k2=3.51×107 N/m, the damping ratio of the structure is 5%, and the basic natural vibration period of the structure is 3.02 s. The abutment is built on a piled foundation and buried (its stiffness is thus large). Considering the spring support of the soil layer on the pile foundation of the abutment and the soil filling behind the abutment, k1=5.40×108. The mass of the abutment and the slabs is estimated as m1=3.16×105 kg. The damping ratio of the abutment is 5%, and its period of natural vibration is 0.15 s. The width of the expansion joint is 0.1 m. The Hertz-damping contact element model is used between the main beam and the abutment. The spring stiffness constant kc=1.25×109 N/m3/2, and the recovery coefficient of the concrete structure is 0.65. Plastic hinges may develop at the bottom of the main pier under strong seismic action. In the elastic-plastic analysis, the elastic-plastic model of the main pier is considered to conform to the assumption of the Clough model. The longitudinal yield displacement of the main pier is 0.16 m, and the postyield stiffness is 0.1 times the elastic stiffness, that is k¯2=3.51×106 N/m.Figure 2
Elevation view of the T-shaped rigid frame bridge (all dimensions: cm).Since there is no Hertz-damping contact element model available in commercial software, the improved precise pounding algorithm proposed herein is difficult to be verified using commercial software. The calculation results of the example were compared with those calculated using the traditional precise integration method [17] to indirectly verify the proposed method.The bridge site is a Class II site with a site characteristic period of 0.35 s and a peak basic seismic acceleration of0.1g. The time-history analysis of rare earthquakes was carried out according to the specifications for seismic design of highway bridges. The input seismic wave was the EW component of the El-centro seismic wave. The duration analyzed was 53.44 s, and the peak acceleration of the seismic wave was adjusted to 0.17g. Two types of precise integration methods (PIM) were used for analysis, and the step size of separation state was set to 0.02 s, and that in the contact state was set to 0.002 s. The time histories of displacement, velocity, and acceleration of the main pier top calculated by the two methods are shown in Figures 3∼5 . The pounding of the structure mainly occurred in the first 21 s. The history of the pounding force at the expansion joints on both sides is shown in Figure 6, and the relationship between the pounding force and penetration displacement is illustrated in Figure 7. In addition, the restorative force-displacement hysteresis curve of the main pier is plotted in Figure 8.Figure 3
Time-history curve of pier top displacement under rare earthquake excitation.Figure 4
Time-history curve of pier top velocity under rare earthquake excitation.Figure 5
Acceleration time-history curve of pier top under rare earthquake excitation.Figure 6
Pounding force time-history curve under rare earthquake excitation.Figure 7
Pounding force-penetration displacement curve under rare earthquake excitation.Figure 8
Restoring force-displacement curve of main pier under rare earthquake excitation.By comparing the results of the two time-history analysis methods, the time-history curves of displacement, velocity, and acceleration of the main pier top and the time-history curves of pounding forces on both sides are almost identical, which indicates that the precision of the improved PIM is not reduced compared with that of the traditional PIM. From Figure7, the pounding force-penetration displacement curve of the two contact points conforms to the hypothesis of the Hertz-damping contact element model. Under the action of a rare earthquake, the main pier enters the plastic stage, and the hysteresis model of restorative force-displacement conforms to the assumption of the Clough model. It is shown that the two methods can be used to analyse the elastic-plastic earthquake-induced response of structures and simulate earthquake-induced processes accurately. On the same computing platform, the program based on the traditional precise integration method and the program based on the improved precise integration method run three times. The average calculation time of the traditional precise integration method is 6.49 s, while the average calculation time of the improved precise integration method is 3.94 s, which is 39.3% lower than the traditional method, which means that the improved method has higher calculation efficiency. This is mainly because the improved precise integration method avoids the inverse operation of the system matrix, and the new order reduction method makes the matrix multiplication less. It will save a lot of time by using the improved precise integration method in the following response spectrum analysis.


## 3.1. Pounding Dynamic Differential Equation

The beam ends pounding model shown in Figure1 contains three SDOF particles, in which the mass, damping, and stiffness of the left and right particles are the same. The pounding force time-history between the middle particle and the left particle is assumed to be f1t, and that between the right particle and the middle particle is assumed to be f2t. The dynamic balance equation of the system is as follows:(1)MX¨+CX˙+KX+Ft=−Mgt,where M=m1000m2000m1, C=c1000c2000c1, K=k1000k2000k1, X=x1x2x3, and Ft=f1t−f1t+f2t−f2t.The structural damping coefficient is given by(2)ci=2ξimiki,i=1,2.The Hertz damping contact element model is used for the pounding between particles, and the viscous damper is added on the basis of the Hertz model to allow energy loss.The spring stiffness and damping coefficient are variables. In the contact stage, spring stiffnessβ and damping coefficient ηare related to the penetration displacement of the two particles, as given by the following equations:(3)βit=kcxi−xi+1−gp,(4)ηit=ζxi−xi+1−gp3/2,(5)ζ=3kc1−e24u˙i−u˙i+1,where kc, e, and gp are, respectively, the spring constant, recovery coefficient, and initial gap of the contact element. xi and xi+1 are the displacements of two adjacent particles at time t, and u˙iu˙i+1 are the velocities of two adjacent particles before contact.The pounding force on a contact element is obtained as follows:(6)fit=βitxi−xi+1−gp+ηitx˙i−x˙i+1.In general, the pounding force can be treated as the hysteresis force, that is, the displacement and velocity calculated in the previous time step are used to calculate the pounding force with equation (6), which is incorporated into the right-hand side of the equation as the known load; therefore, the calculated pounding force always lags by one time step. To overcome this problem, the pounding force column vector of the contact element in the model can be written as follows:(7)Ft=βtX+ηtX˙+γtgp,where βtX+γtgp is the elastic restoring force vector and ηtX˙ is the damping force vector. In detail, βt=β1t−β1t0−β1tβ1t+β2t−β2t0−β2tβ2t refers to the stiffness matrix of contact elements at time t. ηt=η1t−η1t0−η1tη1t+η2t−η2t0−η2tη2t refers to the damping coefficient matrix of contact elements at time t, and γt=−β1tβ1t+β2t−β2t is the column vector to calculate elastic restoring force related to the initial gap.The decomposed pounding force column vector is introduced into the pounding dynamic differential equation, which can be rewritten as follows:(8)MX¨+C+ηtX˙+K+βtX=−Mgt−γtgp.


## 3.2. Improved Precise Integration Method

To reduce the order of the dynamic balance equation of the system, we introduce vectorp. Let p=MX˙, then X˙=M−1p in the dynamic equilibrium equation, which can be written as follows:(9)p˙=−K+βtX−C+ηtM−1p−γtgp−Mgt.Assumingv=Xp, we have(10)v˙=Hv+ft,where the system matrix H=0M−1−K+βt−C+ηtM−1 and nonhomogeneous term ft=0−γtgp−Mgt.Assumming that the solution of the reduced first-order equation of state at timetk is vtk, the system matrix H is time-variant, so tk is denoted asHtk, and the general solution at time tk+1is as follows:(11)vtk+1=expHtk⋅Δtvtk+∫tktk+1expHtktk+1−τfτdτ.The exponential matrix of the first term on the right of the above equation is represented as follows:(12)T=expHtk⋅Δt=expHtk⋅Δt2N2N=expHtk⋅s2N=I+Ta2N.The matrixTa is expanded by Taylor series as follows:(13)Ta=∑l=1∞Htk⋅sll!.In equation (13), the truncation order of the series expansion is five, and N is 20, so a more accurate solution can be obtained. After matrix Ta is found, matrix T is found by the addition principle as applied to an exponential function. Specifically, the matrix is iterated N times the following equations:(14)Ta=2Ta+Ta×Ta,(15)T=I+Ta.The series method is used to solve the nonhomogeneous integral, supposes=τ−tk, then(16)∫tktk+1expHtktk+1−τfτdτ=∫0ΔtexpHtkΔt−sftk+sds=T∫0Δtexp−Htksftk+sds.Suppose the nonhomogeneous term varies linearly over timetiti+1, that is, ftk+s=ftk+Δf⋅s. Then, the solution of the integral term of equation (16) is as follows:(17)∫0Δt∑l=0∞−Htksll!ftk+Δf⋅sds=∑l=0∞−Htkll!ftk⋅Δtl+1l+1+Δf⋅Δtl+2l+2.Suppose the truncation order of the series expansionl in equation (17) is five, then, the precise recursion format of the pounding dynamic differential equation is as follows:(18)vtk+1=Tvtk+∑l=05−Htkll!ftk⋅Δtl+1l+1+Δf⋅Δtl+2l+2.The elements in the first half of the vectorvtk+1 constitute the displacement X at time tk+1, then the rest elements constitute the vector p. The velocity at time tk+1 can be calculated by X˙=M−1p, and substitution into the dynamic balance equation yields the acceleration at time tk+1.The above derivation process is derived from the dynamic balance equation of the full state. For problems considering structural nonlinearity, the secant stiffness matrix is time-varying, so it needs to be solved by iteration. In some cases, the secant stiffness will be very large or very small, resulting in significant numerical error in the calculation; therefore, the incremental dynamic balance equation is usually adopted after the structure is deemed elastic-plastic. For the incremental dynamic balance equation, the precise integration method is also applicable after slight modification. The recursive solution requires that the increment from the previous step is deferred, and the increment of the next step is deferred. The state of the structure at a given time must be obtained by incremental superposition. In the literature [17], a precise integration method for solving the incremental dynamic balance equation of structural nonlinear problems was examined: in the present work, the relevant results were cited along with the description of the program design.


## 3.3. Precise Pounding Algorithm

The pounding force column vector is introduced when the dynamic differential equation of the pounding is established. It should be noted that the state of the structure in the earthquake-induced response can be divided into separate, and contact, states. The current state of the structure can be obtained by judging the displacement difference between the ends of the two contact elements. In the contact state, because the stiffness and damping coefficients of the pounding spring are time-variant, the system matrix is also time-variant. Upon separation, the elements inβt, ηt, and γt are all zero, the system matrix of elements is the time-varying stiffness matrix related to the linear elastic stage or incremental dynamic balance equation of the tangent stiffness constant and plastic phase, the elements in the system matrix are constant, as are the changes in secant stiffness or tangent stiffness in the plastic region, meaning that the system matrix is still time-variant.In the process of the step-by-step detailed recursive solution of the pounding dynamic equation, the separation-contact state transition should also be solved. By dividing the step sizeΔt into 2N substeps, the precise integration method can provide a high-precision solution. In the state of separation, the integral step size Δt can be larger, which is conducive to reducing the computational burden. In the contact state, due to the short duration of the process itself, a smaller integral step size should be adopted. In the program design, the integral step Δt1 of the separation state is set to the sampling interval of seismic wave, and the integral step of the pounding stage is set to 1% to 10% of Δt1.In the recursion process used in the precise integration method, the pounding phenomenon is judged by the value ofxi−xi+1−gp. Assume that time tk is xi−xi+1−gp<0, and time tk+1 is xi−xi+1−gp≥0, indicating the exact time of pounding. The displacement, velocity, and acceleration at time tk and time tk+1 have been deduced by the precise integration method of separate states. According to the assumption that the acceleration changes linearly from tk to tk+1, the structural displacement at time tk+τ can be given by the following equation:(19)xitk−xi+1tk+x˙itk−x˙i+1tkτ+x˙itk−x˙i+1tkτ22+x¨itk+1−x˙itk−x¨i+1tk+1−x˙i+1tkΔtτ36−gp=0.The equation is a cubic equation containing one variable, and the solution in the interval0Δt can be derived by using the dichotomy method, and time tk+τ is the initial time of pounding. The structural state at the exact start time tk of pounding can be obtained by changing the step size to τ, and then the integral step size in the contact state is changed to Δt2, and the system matrix H is set at each step. The structural displacement, velocity, and acceleration under pounding action can be solved by the iterative method until xi−xi+1−gp≤0. At this point, the structure changes from a state of contact to one of separation. Similarly, equation (19) can be used to find the contact separation time. After entering the separation state, the integral step size is changed again toΔt1, and the recursive solution of the separation state is started. The above processes can be conducted alternately until the input seismic wave analysis is completed, and the time-history response of the structure is obtained.


## 3.4. Algorithm Verification

The bilateral pounding model is abstracted from a 2 × 85 m T-shaped rigid frame bridge. The elevation of the bridge is shown in Figure2. In the abstract model, the mass of the main bridge is concentrated at the centroid of the main girder on the pier top, m2=8.12×106 kg, the height of the main pier is 52 m, the calculated antipush stiffness k2=3.51×107 N/m, the damping ratio of the structure is 5%, and the basic natural vibration period of the structure is 3.02 s. The abutment is built on a piled foundation and buried (its stiffness is thus large). Considering the spring support of the soil layer on the pile foundation of the abutment and the soil filling behind the abutment, k1=5.40×108. The mass of the abutment and the slabs is estimated as m1=3.16×105 kg. The damping ratio of the abutment is 5%, and its period of natural vibration is 0.15 s. The width of the expansion joint is 0.1 m. The Hertz-damping contact element model is used between the main beam and the abutment. The spring stiffness constant kc=1.25×109 N/m3/2, and the recovery coefficient of the concrete structure is 0.65. Plastic hinges may develop at the bottom of the main pier under strong seismic action. In the elastic-plastic analysis, the elastic-plastic model of the main pier is considered to conform to the assumption of the Clough model. The longitudinal yield displacement of the main pier is 0.16 m, and the postyield stiffness is 0.1 times the elastic stiffness, that is k¯2=3.51×106 N/m.Figure 2
Elevation view of the T-shaped rigid frame bridge (all dimensions: cm).Since there is no Hertz-damping contact element model available in commercial software, the improved precise pounding algorithm proposed herein is difficult to be verified using commercial software. The calculation results of the example were compared with those calculated using the traditional precise integration method [17] to indirectly verify the proposed method.The bridge site is a Class II site with a site characteristic period of 0.35 s and a peak basic seismic acceleration of0.1g. The time-history analysis of rare earthquakes was carried out according to the specifications for seismic design of highway bridges. The input seismic wave was the EW component of the El-centro seismic wave. The duration analyzed was 53.44 s, and the peak acceleration of the seismic wave was adjusted to 0.17g. Two types of precise integration methods (PIM) were used for analysis, and the step size of separation state was set to 0.02 s, and that in the contact state was set to 0.002 s. The time histories of displacement, velocity, and acceleration of the main pier top calculated by the two methods are shown in Figures 3∼5 . The pounding of the structure mainly occurred in the first 21 s. The history of the pounding force at the expansion joints on both sides is shown in Figure 6, and the relationship between the pounding force and penetration displacement is illustrated in Figure 7. In addition, the restorative force-displacement hysteresis curve of the main pier is plotted in Figure 8.Figure 3
Time-history curve of pier top displacement under rare earthquake excitation.Figure 4
Time-history curve of pier top velocity under rare earthquake excitation.Figure 5
Acceleration time-history curve of pier top under rare earthquake excitation.Figure 6
Pounding force time-history curve under rare earthquake excitation.Figure 7
Pounding force-penetration displacement curve under rare earthquake excitation.Figure 8
Restoring force-displacement curve of main pier under rare earthquake excitation.By comparing the results of the two time-history analysis methods, the time-history curves of displacement, velocity, and acceleration of the main pier top and the time-history curves of pounding forces on both sides are almost identical, which indicates that the precision of the improved PIM is not reduced compared with that of the traditional PIM. From Figure7, the pounding force-penetration displacement curve of the two contact points conforms to the hypothesis of the Hertz-damping contact element model. Under the action of a rare earthquake, the main pier enters the plastic stage, and the hysteresis model of restorative force-displacement conforms to the assumption of the Clough model. It is shown that the two methods can be used to analyse the elastic-plastic earthquake-induced response of structures and simulate earthquake-induced processes accurately. On the same computing platform, the program based on the traditional precise integration method and the program based on the improved precise integration method run three times. The average calculation time of the traditional precise integration method is 6.49 s, while the average calculation time of the improved precise integration method is 3.94 s, which is 39.3% lower than the traditional method, which means that the improved method has higher calculation efficiency. This is mainly because the improved precise integration method avoids the inverse operation of the system matrix, and the new order reduction method makes the matrix multiplication less. It will save a lot of time by using the improved precise integration method in the following response spectrum analysis.


## 4. Earthquake-Induced Spectrum Analysis

### 4.1. Selection of Seismic Waves

Since the first earthquake wave (the El-centro event) was recorded in 1940, many earthquake records have been acquired. In general time-history analysis, a small number of applicable seismic waves must be selected from a large number of seismic records for analysis. According to the site type, aiming at designing acceleration response spectrum, the wave selection method adopts synthetic seismic waves compatible with it, or natural seismic waves with better compatibility in a statistical sense. A total of eighteen seismic waves were selected from the Peer database (Table1). On this basis, the peak value of the selected wave is adjusted, and the peak value of seismic wave is 0.17g except during PGA variable parameter analysis. The period of the reaction spectrum ranges from 0.1 s to 6 s, the step size is 0.1 s, and the corresponding pier stiffness is 8.90×106 and 3.21×1010 N/m.Table 1
Selected seismic waves.
No.EarthquakeObservation stationPGA (g)Magnitude1Anza (Horse Cany), 1980/02/255047 Rancho de Anza0.092M (4.9)2Cape Mendocino, 1992/04/2589486 Fortuna0.116M (7.1)3Chi-Chi, Taiwan, 1999/09/20CHY0220.065M (7.6)4Chi-Chi, Taiwan, 1999/09/20CHY0290.277M (7.6)5Coyote Lake, 1979/08/0657383 Gilroy Array #60.434M (5.7)6Coyote Lake, 1979/08/061492 SJB Overpass, Bent 3 g.l.0.124M (5.7)7Duzce, Turkey, 1999/11/12Sakarya0.023M (7.1)8Imperial Valley, 1979/10/155051 Parachute Test Site0.204M (6.5)9Kern County, 1952/07/211095 Taft Lincoln School0.178M (7.4)10Kobe, 19950 KJMA0.821M (6.9)11Kocaeli, Turkey, 1999/08/17Arcelik0.149M (7.4)12Landers, 1992/06/2823 Coolwater0.283M (7.3)13Livermore, 1980/01/2758219 APEEL 3E Hayward CSUH0.028M (5.4)14Loma Prieta, 1989/10/181652 Anderson Dam (Downstream)0.244M (6.9)15Loma Prieta, 1989/10/1858373 APEEL 10–Skyline0.103M (6.9)16Lytle Creek, 1970/09/12290 Wrightwood–6074 Park Dr.0.162M (5.4)17Northridge, 1994/01/1724278 Castaic–Old Ridge Route0.662M (7.1)18Imperial Valley, 1940/05/18El-Centro (EW)0.348M (7.1)


### 4.2. Pounding Spectrum Analysis of Selected Waves

Based on the improved PIM, a multipoint earthquake-induced response spectrum analysis program was developed on the MATLAB™ platform. In the seismic impact response spectrum analysis, the longitudinal natural vibration period of the main bridge is the basic variable. In the case of the fixed mass of the main bridge, the stiffness of the pier k2 is calculated according to the period of natural vibration. The values of other basic analysis parameters have been given in the example. Although the program can perform elastic-plastic analysis, it ignores the influence of elastic-plastic behaviour of the bridge pier in response spectrum analysis with considering pounding or not. First, the program was used to calculate the displacement spectrum and dynamic amplification coefficient spectrum of the simplified single-degree-of freedom structure of main pier and main beam under eighteen seismic waves without considering pounding, as shown in Figure 9. Then, the displacement spectrum, dynamic amplification coefficient spectrum, pounding force spectrum, and the number of pounding event response spectrum of the simplified SDOF of main pier and main beam were calculated, as shown in Figure 10.Figure 9
Analytical results without considering pounding. (a) Displacement spectrum. (b) Dynamic amplification factor spectrum.
(a)(b)Figure 10
Analytical results considering pounding. (a) Displacement spectrum. (b) Dynamic amplification factor spectrum. (c) Pounding force spectrum. (d) Number of pounding event response spectrum.
(a)(b)(c)(d)The standard dynamic amplification spectrum is shown in Figures9(b) and 10(b). It is calculated using the following formula, which was cited from the specifications for seismic design of highway bridges:(20)At=Amax0.6TT0,T≤T0,Amax,T0<T≤Tg,AmaxTgT,Tg<T≤10,where T0=0.1s, Tg=0.35s refers to the site characteristic period, and Amax=2.5CiCsCd refers to the maximum value of dynamic amplification factor, in which Ci=1.0, Cs=1.7, Cd=1.0, respectively.The structural displacement spectrum without considering pounding shows that, at the same PGA level, the displacement response of the structure is strongly related to the input seismic waveform, and the displacement spectra corresponding to the seven selected seismic waves never exceed 0.1 m, while the displacements associated with the Duzce and Kocaeli seismic waves increase rapidly with the displacement extreme value exceeding 1 m. The acceleration dynamic amplification factor changes significantly in a short period of less than 3 s, and the peak value of the spectral curve appears in a period of 0.2 s to 1 s, and tends to be smooth during long-period motions of more than 3 s, and the spectral value of the acceleration dynamic coefficient does not exceed 1.5.Considering the pounding between the beam ends on both sides of the intermediate structure and the abutment, pounding will occur when the difference in displacement between the intermediate structure and the abutment exceeds the reserved expansion joint gap. The longitudinal displacement of the intermediate structure is limited by the abutments with high stiffness on both sides, so the longitudinal displacement is greatly reduced, which is reflected in the displacement spectrum curve, and the long-period spectrum curve tends to be smooth. It can be seen from the spectrum curve of acceleration dynamic amplification coefficient that the acceleration dynamic amplification coefficient does not change in the short period when the natural vibration period is less than 1.5 s, but the spectrum curve of dynamic amplification coefficient of six seismic waves increases significantly when the natural vibration period exceeds 1.5 s. The design acceleration response spectrum given in the design code does not consider the influence of pounding, so the response spectrum analysis method should be avoided when assessing the earthquake-induced response of structures prone to pounding with a long period. Compared with the pounding force spectrum and the number of pounding events as a response spectrum, the structural pounding mainly occurs with a period of more than 1.5 s, which shows that the effect of pounding force on the acceleration is huge. From the analysis of the pounding force spectrum and the number of pounding event response spectrum curve, only eleven of eighteen seismic waves induce pounding, and the pounding initiation period, pounding force, and the number of pounding events as response spectra associated with eleven seismic waves are quite different. This shows that the pounding response of structures is closely related to the spectral characteristics of the input seismic wave. By analyzing the pounding force spectrum curve and the number of pounding events as a response spectrum curve under the same input excitation, it is found that the pounding force and the number of pounding event response spectrum increase rapidly at first and then do so more slowly with increasing period from the onset of pounding, but this does not preclude abrupt changes within individual periods, which shows that the pounding force and the number of pounding events (as response spectra) are sensitive to the natural period of vibration of the structure.


## 4.1. Selection of Seismic Waves

Since the first earthquake wave (the El-centro event) was recorded in 1940, many earthquake records have been acquired. In general time-history analysis, a small number of applicable seismic waves must be selected from a large number of seismic records for analysis. According to the site type, aiming at designing acceleration response spectrum, the wave selection method adopts synthetic seismic waves compatible with it, or natural seismic waves with better compatibility in a statistical sense. A total of eighteen seismic waves were selected from the Peer database (Table1). On this basis, the peak value of the selected wave is adjusted, and the peak value of seismic wave is 0.17g except during PGA variable parameter analysis. The period of the reaction spectrum ranges from 0.1 s to 6 s, the step size is 0.1 s, and the corresponding pier stiffness is 8.90×106 and 3.21×1010 N/m.Table 1
Selected seismic waves.
No.EarthquakeObservation stationPGA (g)Magnitude1Anza (Horse Cany), 1980/02/255047 Rancho de Anza0.092M (4.9)2Cape Mendocino, 1992/04/2589486 Fortuna0.116M (7.1)3Chi-Chi, Taiwan, 1999/09/20CHY0220.065M (7.6)4Chi-Chi, Taiwan, 1999/09/20CHY0290.277M (7.6)5Coyote Lake, 1979/08/0657383 Gilroy Array #60.434M (5.7)6Coyote Lake, 1979/08/061492 SJB Overpass, Bent 3 g.l.0.124M (5.7)7Duzce, Turkey, 1999/11/12Sakarya0.023M (7.1)8Imperial Valley, 1979/10/155051 Parachute Test Site0.204M (6.5)9Kern County, 1952/07/211095 Taft Lincoln School0.178M (7.4)10Kobe, 19950 KJMA0.821M (6.9)11Kocaeli, Turkey, 1999/08/17Arcelik0.149M (7.4)12Landers, 1992/06/2823 Coolwater0.283M (7.3)13Livermore, 1980/01/2758219 APEEL 3E Hayward CSUH0.028M (5.4)14Loma Prieta, 1989/10/181652 Anderson Dam (Downstream)0.244M (6.9)15Loma Prieta, 1989/10/1858373 APEEL 10–Skyline0.103M (6.9)16Lytle Creek, 1970/09/12290 Wrightwood–6074 Park Dr.0.162M (5.4)17Northridge, 1994/01/1724278 Castaic–Old Ridge Route0.662M (7.1)18Imperial Valley, 1940/05/18El-Centro (EW)0.348M (7.1)


## 4.2. Pounding Spectrum Analysis of Selected Waves

Based on the improved PIM, a multipoint earthquake-induced response spectrum analysis program was developed on the MATLAB™ platform. In the seismic impact response spectrum analysis, the longitudinal natural vibration period of the main bridge is the basic variable. In the case of the fixed mass of the main bridge, the stiffness of the pier k2 is calculated according to the period of natural vibration. The values of other basic analysis parameters have been given in the example. Although the program can perform elastic-plastic analysis, it ignores the influence of elastic-plastic behaviour of the bridge pier in response spectrum analysis with considering pounding or not. First, the program was used to calculate the displacement spectrum and dynamic amplification coefficient spectrum of the simplified single-degree-of freedom structure of main pier and main beam under eighteen seismic waves without considering pounding, as shown in Figure 9. Then, the displacement spectrum, dynamic amplification coefficient spectrum, pounding force spectrum, and the number of pounding event response spectrum of the simplified SDOF of main pier and main beam were calculated, as shown in Figure 10.Figure 9
Analytical results without considering pounding. (a) Displacement spectrum. (b) Dynamic amplification factor spectrum.
(a)(b)Figure 10
Analytical results considering pounding. (a) Displacement spectrum. (b) Dynamic amplification factor spectrum. (c) Pounding force spectrum. (d) Number of pounding event response spectrum.
(a)(b)(c)(d)The standard dynamic amplification spectrum is shown in Figures9(b) and 10(b). It is calculated using the following formula, which was cited from the specifications for seismic design of highway bridges:(20)At=Amax0.6TT0,T≤T0,Amax,T0<T≤Tg,AmaxTgT,Tg<T≤10,where T0=0.1s, Tg=0.35s refers to the site characteristic period, and Amax=2.5CiCsCd refers to the maximum value of dynamic amplification factor, in which Ci=1.0, Cs=1.7, Cd=1.0, respectively.The structural displacement spectrum without considering pounding shows that, at the same PGA level, the displacement response of the structure is strongly related to the input seismic waveform, and the displacement spectra corresponding to the seven selected seismic waves never exceed 0.1 m, while the displacements associated with the Duzce and Kocaeli seismic waves increase rapidly with the displacement extreme value exceeding 1 m. The acceleration dynamic amplification factor changes significantly in a short period of less than 3 s, and the peak value of the spectral curve appears in a period of 0.2 s to 1 s, and tends to be smooth during long-period motions of more than 3 s, and the spectral value of the acceleration dynamic coefficient does not exceed 1.5.Considering the pounding between the beam ends on both sides of the intermediate structure and the abutment, pounding will occur when the difference in displacement between the intermediate structure and the abutment exceeds the reserved expansion joint gap. The longitudinal displacement of the intermediate structure is limited by the abutments with high stiffness on both sides, so the longitudinal displacement is greatly reduced, which is reflected in the displacement spectrum curve, and the long-period spectrum curve tends to be smooth. It can be seen from the spectrum curve of acceleration dynamic amplification coefficient that the acceleration dynamic amplification coefficient does not change in the short period when the natural vibration period is less than 1.5 s, but the spectrum curve of dynamic amplification coefficient of six seismic waves increases significantly when the natural vibration period exceeds 1.5 s. The design acceleration response spectrum given in the design code does not consider the influence of pounding, so the response spectrum analysis method should be avoided when assessing the earthquake-induced response of structures prone to pounding with a long period. Compared with the pounding force spectrum and the number of pounding events as a response spectrum, the structural pounding mainly occurs with a period of more than 1.5 s, which shows that the effect of pounding force on the acceleration is huge. From the analysis of the pounding force spectrum and the number of pounding event response spectrum curve, only eleven of eighteen seismic waves induce pounding, and the pounding initiation period, pounding force, and the number of pounding events as response spectra associated with eleven seismic waves are quite different. This shows that the pounding response of structures is closely related to the spectral characteristics of the input seismic wave. By analyzing the pounding force spectrum curve and the number of pounding events as a response spectrum curve under the same input excitation, it is found that the pounding force and the number of pounding event response spectrum increase rapidly at first and then do so more slowly with increasing period from the onset of pounding, but this does not preclude abrupt changes within individual periods, which shows that the pounding force and the number of pounding events (as response spectra) are sensitive to the natural period of vibration of the structure.


## 5. Impact Analysis of Pounding Spectrum Parameters

The analysis of results from inputting each of eighteen seismic waves in the previous section show that the response of the structure is closely related to the selection of seismic waveform under the same peak ground motion; because of the randomness of earthquakes, it is difficult to characterize the structural response spectrum with a single seismic wave; therefore, the average spectrum of multiple seismic waves is used in parameter influence analysis, that is, the average value of the selected seismic wave response spectrum curves is used to characterize the spectra of this group of seismic waves. Seven seismic waves without pounding were eliminated. Due to the influence of the soil filling behind the abutment, it is difficult to determine the stiffness and mass of the abutment when it is simplified as an SDOF structure: to discuss the influence of the two parameters on the structural response, a variable parameter analysis is conducted. Among the contact element parameters, besides the structural clearance, the values of the pounding stiffness parameter and the recovery coefficient are also controversial, and the variable parameter analysis is also aimed at assessing how the two parameters affect the response spectrum. In addition, PGA is an important factor affecting the earthquake-induced response of structures, and variable parameter analysis of PGA is also conducted. It is also worth mentioning that the material nonlinearity of piers and abutments is not included in this part of parametric analysis. It is mainly considered that the collision between the main beam and the abutments on both sides limits the excessive displacement of the pier top, which limits the generation of the plastic hinge at the bottom of the pier, and this consideration can reduce the complexity of the analysis problem.


### 5.1. Stiffness of the Abutment

The initial stiffnessk1=5.40×108 N/m of the abutment is about 15.4 times that of the main pier. The stiffness of the abutment depends on the geological conditions, the geometric dimensions of the abutment, especially, the height of the abutment, and the support stiffness of the soil filling behind the abutment. Considering the changes of the above factors, it is possible for the stiffness of the abutment to vary from 0.1 to 10 times of the initial stiffness. Thus, the abutment stiffness is set to each of k1/10, k1/4, k1/2, 2k1, 4k1, and 10k1 with the other structural parameters and pounding parameters remaining unchanged. The pounding spectrum analysis of seven groups of abutment stiffness parameters is conducted. The average displacement spectrum, dynamic coefficient spectrum, pounding force spectrum, and the number of pounding events as a response spectrum were calculated for different abutment stiffnesses as shown in Figures 11(a)–11(d).Figure 11
Average spectral analysis results with different stiffnesses of abutment. (a) Displacement spectrum. (b) Dynamic amplification factor spectrum. (c) Pounding force spectrum. (d) Number of pounding event response spectrum.
(a)(b)(c)(d)Results show that the abutment stiffness has no effect on the displacement spectrum and dynamic amplification coefficient spectrum of short-period structures with a period of 1.5 s as the boundary because such short-period structures do not suffer structural pounding. For long-period structures, due to the pounding, the stiffness of the abutment affects the structural displacement, dynamic amplification factor, pounding force, and the number of pounding events (as response spectra). With the decrease in the abutment stiffness, the longitudinal displacement of the main girder increases, while the dynamic amplification coefficient decreases, and the pounding force also shows a decreasing trend, while the change in the number of pounding event response spectrum is irregular. From the variation of abutment stiffness from2k1 to 10k1, the curves of displacement spectrum, dynamic amplification coefficient spectrum, and pounding force spectrum are similar, suggesting that, when the stiffness of the abutment reaches a certain value, the difference in structural earthquake-induced response arising therefrom can be ignored.


### 5.2. Abutment Quality

The concrete quality in the abutment and slabs is considered when the abutment mass is calculated, and the initial mass ism1=3.16×105 kg. The quality of the abutment depends on the size of the abutment and the scope of the soil filling behind the abutment. The abutment mass varies from 0.1 to 10 times the initial value is possibly in considering various limit cases. Pounding spectrum analysis of seven groups of abutment mass parameters is performed while other structural parameters and pounding parameters are kept unchanged: the abutment mass is set to m1/10, m1/4, m1/2, 2 m1, 4 m1, and 10 m1. The average displacement spectrum, dynamic coefficient spectrum, pounding force spectrum, and the number of pounding events (as response spectra) arising from provision of different abutment mass are shown in Figures 12(a)–12(d).Figure 12
Average spectrum analysis results with different abutment mass. (a) Displacement spectrum. (b) Dynamic amplification factor spectrum. (c) Pounding force spectrum. (d) Number of pounding event response spectrum.
(a)(b)(c)(d)Comparing the displacement spectrum, dynamic coefficient spectrum, pounding force spectrum, and the number of pounding events (as a response spectrum) for different abutment mass, it is found that although the abutment mass changes significantly, the response spectrum curve changes little and no trend emerges therein. The influence of abutment mass variation on the earthquake-induced response of structures can be ignored, so only the mass of abutment concrete need be included in the earthquake-induced pounding analysis.


### 5.3. Stiffness of the Contact Element

The calculation of stiffness parameters of the Hertz-damping contact element model for simulating pounding between the main girder and abutment is inconclusive. Goldsmith [31] equated the impactor to a sphere, which is suggested to be derived from the formula related to the impactor material and the radius of the equivalent sphere. Xu et al. [32] derived a stiffness calculation formula for the Hertz-damping model based on the straight-rod coaxial pounding model, but it is only applicable to the simulation of pounding between beams. The present research does not consider the calculation of the contact element stiffness: the stiffness parameter is approximately calculated to be 0.2 times the longitudinal stiffness of the half-span main girder. The initial stiffness constant is kc=1.25×109 N/m3/2, and it is assumed to vary from 0.1 times to 10 times the initial stiffness constant. Keeping other structural parameters and pounding parameters unchanged, and changing the values of the pounding stiffness constant kc/10, kc/4, kc/2, 2kc, 4kc, and 10kc, the average displacement spectrum, dynamic coefficient spectrum, pounding force spectrum, and the number of pounding events (as a response spectrum) calculated using pounding stiffness constants of different contact elements are shown in Figures 13(a)–13(d).Figure 13
Average spectrum analysis results for different pounding stiffness constants. (a) Displacement spectrum. (b) Dynamic amplification factor spectrum. (c) Pounding force spectrum. (d) Number of pounding event response spectrum.
(a)(b)(c)(d)Analysis of Figures13(a)–13(d) shows that, once contact pounding occurs, the stiffness of the contact element exerts a significant influence on the earthquake-induced response of the structure. The greater the stiffness of the contact element, the smaller the structural displacement response, and the greater the acceleration dynamic amplification coefficient, the greater the pounding force. In long-period events (more than 3 s), the increased pounding stiffness will also tend to increase the number of pounding events, but the regularity in this trend is poor. At present, the stiffness of contact pounding between adjacent beams is reduced by setting rubber pads at adjacent beam ends, thus reducing the pounding force.


### 5.4. Pounding Recovery Coefficient

Jankowski [33] points out that the coefficient of restitution is not a fixed value, and the coefficient of restitution of concrete is recommended to be 0.43∼0.78. The larger the pounding recovery coefficient, the more it resembles elastic pounding, and the smaller it is, the more it resembles plastic pounding. In this part, the initial pounding recovery coefficient of the Hertz-damping model is such that e=0.65. While keeping other parameters unchanged, the pounding spectrum was calculated under four sets of parameters with pounding recovery coefficients of 0.45, 0.55, 0.75, and 0.85, and compared with the calculated initial parameters, as shown in Figures 14(a)–14(d).Figure 14
Average spectrum analysis results of different pounding recovery coefficients. (a) Displacement spectrum. (b) Dynamic amplification factor spectrum. (c) Pounding force spectrum. (d) Number of pounding event response spectrum.
(a)(b)(c)(d)Comparing the spectra under different pounding recovery coefficients, it is found that the pounding recovery coefficient mainly affects long-period structures with a period of natural vibration exceeding 2 s. The increase in the recovery coefficient means that the energy dissipated in pounding decreases, and the structural displacement and dynamic coefficient of pounding force increase (in particular, the pounding force of long-period structures becomes more significant), but the number of pounding events (as a response spectrum) shows no clear trend. Increasing the pounding energy dissipation helps to attenuate the earthquake-induced response of the structure.


### 5.5. Peak Ground Acceleration

The peak acceleration induced by a given input seismic wave is taken as an index with which the intensity of that seismic input is characterized. Due to the randomness of the earthquake source location and intensity and the variability of seismic wave transmission, there is significant uncertainty in the peak ground motion acceleration at a given site. Although it is impossible to control the intensity of seismic waves, revealing the influence of PGA on the earthquake-induced response of structures is helpful to understand the earthquake-induced phenomenon as they affect such structures; therefore, a PGA based parametric study was conducted in which structural parameters and pounding parameters are kept unchanged. The peak value of input seismic wave acceleration changed from0.09g to 0.29g in increments of 0.04g, and six groups of spectral analysis were run. The average displacement spectrum, dynamic coefficient spectrum, pounding force spectrum, and the number of pounding events (as a response spectrum) calculated using pounding stiffness constants of different contact elements are shown in Figures 15(a)–15(d).Figure 15
Average spectral analysis results using different PGA parameters. (a) Displacement spectrum. (b) Dynamic amplification factor spectrum. (c) Pounding force spectrum. (d) Number of pounding event response spectrum.
(a)(b)(c)(d)Comparing the calculated results pertaining to six groups of PGA (Figures15(a)–15(d)), it is found that the PGA affects the structural response. With increasing PGA, the dynamic amplification coefficients of displacement and acceleration of the structure increase, the initial pounding period of the structure decreases, and the pounding force and the number of pounding events (as a response spectrum) after pounding also increase to a significant extent. The differences in earthquake-induced response caused by changes to the PGA are the most significant among the five factors assessed in this parametric study, making PGA an important factor in determining earthquake-induced structural responses.


## 5.1. Stiffness of the Abutment

The initial stiffnessk1=5.40×108 N/m of the abutment is about 15.4 times that of the main pier. The stiffness of the abutment depends on the geological conditions, the geometric dimensions of the abutment, especially, the height of the abutment, and the support stiffness of the soil filling behind the abutment. Considering the changes of the above factors, it is possible for the stiffness of the abutment to vary from 0.1 to 10 times of the initial stiffness. Thus, the abutment stiffness is set to each of k1/10, k1/4, k1/2, 2k1, 4k1, and 10k1 with the other structural parameters and pounding parameters remaining unchanged. The pounding spectrum analysis of seven groups of abutment stiffness parameters is conducted. The average displacement spectrum, dynamic coefficient spectrum, pounding force spectrum, and the number of pounding events as a response spectrum were calculated for different abutment stiffnesses as shown in Figures 11(a)–11(d).Figure 11
Average spectral analysis results with different stiffnesses of abutment. (a) Displacement spectrum. (b) Dynamic amplification factor spectrum. (c) Pounding force spectrum. (d) Number of pounding event response spectrum.
(a)(b)(c)(d)Results show that the abutment stiffness has no effect on the displacement spectrum and dynamic amplification coefficient spectrum of short-period structures with a period of 1.5 s as the boundary because such short-period structures do not suffer structural pounding. For long-period structures, due to the pounding, the stiffness of the abutment affects the structural displacement, dynamic amplification factor, pounding force, and the number of pounding events (as response spectra). With the decrease in the abutment stiffness, the longitudinal displacement of the main girder increases, while the dynamic amplification coefficient decreases, and the pounding force also shows a decreasing trend, while the change in the number of pounding event response spectrum is irregular. From the variation of abutment stiffness from2k1 to 10k1, the curves of displacement spectrum, dynamic amplification coefficient spectrum, and pounding force spectrum are similar, suggesting that, when the stiffness of the abutment reaches a certain value, the difference in structural earthquake-induced response arising therefrom can be ignored.


## 5.2. Abutment Quality

The concrete quality in the abutment and slabs is considered when the abutment mass is calculated, and the initial mass ism1=3.16×105 kg. The quality of the abutment depends on the size of the abutment and the scope of the soil filling behind the abutment. The abutment mass varies from 0.1 to 10 times the initial value is possibly in considering various limit cases. Pounding spectrum analysis of seven groups of abutment mass parameters is performed while other structural parameters and pounding parameters are kept unchanged: the abutment mass is set to m1/10, m1/4, m1/2, 2 m1, 4 m1, and 10 m1. The average displacement spectrum, dynamic coefficient spectrum, pounding force spectrum, and the number of pounding events (as response spectra) arising from provision of different abutment mass are shown in Figures 12(a)–12(d).Figure 12
Average spectrum analysis results with different abutment mass. (a) Displacement spectrum. (b) Dynamic amplification factor spectrum. (c) Pounding force spectrum. (d) Number of pounding event response spectrum.
(a)(b)(c)(d)Comparing the displacement spectrum, dynamic coefficient spectrum, pounding force spectrum, and the number of pounding events (as a response spectrum) for different abutment mass, it is found that although the abutment mass changes significantly, the response spectrum curve changes little and no trend emerges therein. The influence of abutment mass variation on the earthquake-induced response of structures can be ignored, so only the mass of abutment concrete need be included in the earthquake-induced pounding analysis.


## 5.3. Stiffness of the Contact Element

The calculation of stiffness parameters of the Hertz-damping contact element model for simulating pounding between the main girder and abutment is inconclusive. Goldsmith [31] equated the impactor to a sphere, which is suggested to be derived from the formula related to the impactor material and the radius of the equivalent sphere. Xu et al. [32] derived a stiffness calculation formula for the Hertz-damping model based on the straight-rod coaxial pounding model, but it is only applicable to the simulation of pounding between beams. The present research does not consider the calculation of the contact element stiffness: the stiffness parameter is approximately calculated to be 0.2 times the longitudinal stiffness of the half-span main girder. The initial stiffness constant is kc=1.25×109 N/m3/2, and it is assumed to vary from 0.1 times to 10 times the initial stiffness constant. Keeping other structural parameters and pounding parameters unchanged, and changing the values of the pounding stiffness constant kc/10, kc/4, kc/2, 2kc, 4kc, and 10kc, the average displacement spectrum, dynamic coefficient spectrum, pounding force spectrum, and the number of pounding events (as a response spectrum) calculated using pounding stiffness constants of different contact elements are shown in Figures 13(a)–13(d).Figure 13
Average spectrum analysis results for different pounding stiffness constants. (a) Displacement spectrum. (b) Dynamic amplification factor spectrum. (c) Pounding force spectrum. (d) Number of pounding event response spectrum.
(a)(b)(c)(d)Analysis of Figures13(a)–13(d) shows that, once contact pounding occurs, the stiffness of the contact element exerts a significant influence on the earthquake-induced response of the structure. The greater the stiffness of the contact element, the smaller the structural displacement response, and the greater the acceleration dynamic amplification coefficient, the greater the pounding force. In long-period events (more than 3 s), the increased pounding stiffness will also tend to increase the number of pounding events, but the regularity in this trend is poor. At present, the stiffness of contact pounding between adjacent beams is reduced by setting rubber pads at adjacent beam ends, thus reducing the pounding force.


## 5.4. Pounding Recovery Coefficient

Jankowski [33] points out that the coefficient of restitution is not a fixed value, and the coefficient of restitution of concrete is recommended to be 0.43∼0.78. The larger the pounding recovery coefficient, the more it resembles elastic pounding, and the smaller it is, the more it resembles plastic pounding. In this part, the initial pounding recovery coefficient of the Hertz-damping model is such that e=0.65. While keeping other parameters unchanged, the pounding spectrum was calculated under four sets of parameters with pounding recovery coefficients of 0.45, 0.55, 0.75, and 0.85, and compared with the calculated initial parameters, as shown in Figures 14(a)–14(d).Figure 14
Average spectrum analysis results of different pounding recovery coefficients. (a) Displacement spectrum. (b) Dynamic amplification factor spectrum. (c) Pounding force spectrum. (d) Number of pounding event response spectrum.
(a)(b)(c)(d)Comparing the spectra under different pounding recovery coefficients, it is found that the pounding recovery coefficient mainly affects long-period structures with a period of natural vibration exceeding 2 s. The increase in the recovery coefficient means that the energy dissipated in pounding decreases, and the structural displacement and dynamic coefficient of pounding force increase (in particular, the pounding force of long-period structures becomes more significant), but the number of pounding events (as a response spectrum) shows no clear trend. Increasing the pounding energy dissipation helps to attenuate the earthquake-induced response of the structure.


## 5.5. Peak Ground Acceleration

The peak acceleration induced by a given input seismic wave is taken as an index with which the intensity of that seismic input is characterized. Due to the randomness of the earthquake source location and intensity and the variability of seismic wave transmission, there is significant uncertainty in the peak ground motion acceleration at a given site. Although it is impossible to control the intensity of seismic waves, revealing the influence of PGA on the earthquake-induced response of structures is helpful to understand the earthquake-induced phenomenon as they affect such structures; therefore, a PGA based parametric study was conducted in which structural parameters and pounding parameters are kept unchanged. The peak value of input seismic wave acceleration changed from0.09g to 0.29g in increments of 0.04g, and six groups of spectral analysis were run. The average displacement spectrum, dynamic coefficient spectrum, pounding force spectrum, and the number of pounding events (as a response spectrum) calculated using pounding stiffness constants of different contact elements are shown in Figures 15(a)–15(d).Figure 15
Average spectral analysis results using different PGA parameters. (a) Displacement spectrum. (b) Dynamic amplification factor spectrum. (c) Pounding force spectrum. (d) Number of pounding event response spectrum.
(a)(b)(c)(d)Comparing the calculated results pertaining to six groups of PGA (Figures15(a)–15(d)), it is found that the PGA affects the structural response. With increasing PGA, the dynamic amplification coefficients of displacement and acceleration of the structure increase, the initial pounding period of the structure decreases, and the pounding force and the number of pounding events (as a response spectrum) after pounding also increase to a significant extent. The differences in earthquake-induced response caused by changes to the PGA are the most significant among the five factors assessed in this parametric study, making PGA an important factor in determining earthquake-induced structural responses.


## 6. Conclusion

Taking a T-shaped rigid frame bridge as the prototype, the earthquake-induced pounding model on both sides was abstracted, the precise integration method for solving the earthquake-induced pounding equation of the system was improved, and the precise pounding algorithm and program were developed. The analysis of an example shows that the calculation results based on the improved precise pounding algorithm are consistent with those based on the traditional precise integration method, and the calculation efficiency of the improved algorithm is improved. Then, based on the improved precise pounding algorithm, earthquake-induced pounding response spectrum analysis and a parametric study were conducted. The conclusions are as follows:(1)
The structural response (when considering pounding) is closely related to the nature of the input seismic wave, and the earthquake-induced response of different seismic waves is quite different under the same peak acceleration. With the increase in the peak acceleration of a given seismic wave, the displacement, acceleration dynamic amplification coefficient, pounding force, and the number of pounding events (as a response spectrum) of the structure increase significantly.(2)
The dynamic amplification coefficient spectrum of acceleration after considering pounding is much greater than that when ignoring the effects of pounding in long-period events; however, the effects of pounding are not considered in the current design acceleration response spectrum, so the response spectrum analysis method should be avoided for the earthquake-induced response of structures prone to long-period pounding phenomena.(3)
The stiffness of the abutment affects the earthquake-induced pounding response, but the quality thereof has little influence on the structural response. It is suggested that the stiffness of the abutment and the filling behind the abutment should be considered in any analysis, whereas the quality of the backfill can be ignored.(4)
The stiffness of the contact element affects the earthquake-induced response of the structure; reducing the stiffness of the contact element can greatly reduce the pounding force and dynamic amplification factor, increasing the energy dissipated upon each impact in the pounding process, while also reducing the pounding response (albeit to a limited extent). From the perspective of seismic pounding control in such structures, it remains necessary to design devices that can reduce the pounding stiffness and increase the energy dissipated therein.Although this study of pounding response spectrum take traditional response spectrum as references, the limitation of this study is obvious. This study only draws some qualitative conclusions through limited parametric analysis, and it is difficult to apply this to the general seismic response spectrum analysis.



---

*Source: 9469405-2021-09-02.xml*