A Microscale Framework for Seismic Stability Analysis of Bridge Pier Rocking Isolation Using the Discrete Element Method

: This study presents a microscale framework for investigating the seismic stability of bridge-pier structures using the discrete element method (DEM), with a focus on rocking isolation mechanisms. The piers and the deck are modeled as rigid blocks that follow rigid body dynamics. The rigid block is modeled as a collection of glued particles with geometrical arrangement and physical properties that mimic an actual block. To facilitate numerical contact points between the base of the block and the flat base wall, smaller particle sizes were introduced at the base of the block. A Hertz contact model was employed to model the interaction between contacting entities for better estimation of the contact constitutive parameters. Validation was performed using well-documented experimental data featuring the free-rocking of a granite stone block as well as existing analytical techniques. DEM simulations were performed on single blocks as well as on a bridge deck-pier system subjected to dynamic and seismic loadings. The study shows the effectiveness of rocking isolation through a comparative analysis of acceleration and angular velocity under varying seismic intensities, with acceleration reduction up to 70% for piers and 60% for the deck in a high-intensity scenario, affirming the potential of rocking isolation as a viable seismic mitigation strategy. The study monitors the structural response, contact mechanics, and energy dissipation of the pier–deck system. The application of the DEM model advances the analysis of bridge pier and deck interactions under seismic loads, providing new insights into the detailed behavior of rocking bridge piers and their potential for seismic isolation.


Introduction
The ability of bridge pier constructions to withstand seismic occurrences is a crucial topic in civil engineering, highlighted by the destructive effects of recent earthquakes on infrastructure globally [1].The nonlinear nature of rocking has drawn the attention of many researchers, using a number of analytical and numerical models, since its original formulation by Housner [2].Traditional methodologies such as static analysis and finite element analysis (FEA) have provided foundational insights into the seismic behavior of these structures.However, the evolution of seismic analysis techniques has been driven by the relentless pursuit of more accurate, efficient, and predictive models that can encapsulate the complex dynamics of earthquake-structure interactions.Modern seismic designs of rigid structures based on foundation systems, designed to prevent uplifting and sliding, have taken center stage, where two distinct designs have emerged: seismic base isolation and rocking isolation.In seismic base isolation designs, additional soft layers are interleaved between the foundation and the superstructure to help contain the seismic displacement demands, thereby ensuring safer structures.Rocking isolation designs, on the other hand, limits the transmission of forces to the structure, as the rocking interfaces help to dissipate the input ground motion energy.
Makris and Rousso [3-5] built upon Housner's [2] groundbreaking work, which demonstrated that a smaller block may overturn due to the high-frequency fluctuations that override the long-duration pulse, whereas a larger block will overturn due to the long-duration pulse.Makris [6] investigated the transient rocking response of a rigid block subjected to trigonometric pulses and near-source ground motions.They discovered that Housner's calculation of the minimum acceleration amplitude of a half-sine pulse required to topple a rigid block was inaccurate, claiming that, contrary to what Housner believed, a block subjected to a half-sine pulse actually overturns throughout its free vibration regime.The study also shows that the normalized pulse acceleration, frequency ratio, and slenderness are the three parameters that control overturning under horizontal excitation, with smaller blocks being more sensitive to peak ground acceleration and larger blocks being more sensitive to incremental ground velocity.A 3D formulation of rocking dynamics of a rigid rectangular prism on a deformable base was studied by the development of two models, the Concentrated Springs Model and the Winkler Model [7].By numerical analysis, sliding and uplift are both considered and represented by fully nonlinear equations of the problem.The behavior of bodies subjected to general ground excitations are then examined using the developed models by Chatzis et al. [7].The rocking motion of a rigid rectangular prism on a moving base is a complex three-dimensional phenomenon.
The concept of rocking isolation in bridge piers represents a paradigm shift in seismic design philosophy.Unlike traditional fixed-base piers, rocking piers are designed to lift off their foundations during an earthquake, allowing the structure to sway and then re-center itself, thus minimizing permanent deformation.This technique is underpinned by the principle of energy dissipation through controlled motion, a concept extensively explored in the work of Mergos et al. [8], where the feasibility and effectiveness of pre-stressed rocking piers in dissipating seismic energy were demonstrated.Makris et al. [6,9,10] examined the seismic behavior of slender, free-standing structures, concentrating on the phenomenon of rocking isolation.They explored historical contexts, such as ancient temple design and Housner's important study in 1963, which provided the groundwork for future research into rocking structures.The study underlined the benefits of rocking isolation as a feasible seismic protection alternative, pointing out that seismic resistance is predominantly caused by member rotational inertia.It also discusses how heavy decks atop slender columns might improve the stability of rocking systems.
Rele et al. [11] explored a controlled-rocking isolation technique on an overpass bridge highway with the use of elastomeric pads, stoppers, and external restrainers to isolate the base, restrict horizontal sliding, and reduce horizontal displacement, respectively.The findings indicated that the technique reduces the forces on the pier, minimizes residual drifts, and enhances post-earthquake serviceability compared to conventional systems.Hussam et al. tested bridge pier models for their acceleration response, seismic displacement, settlement, and failure mechanisms [12].Thomaidis et al. [13,14] concludes that the seismic response of rocking bridges with different levels of pier asymmetry can resist high seismic excitations.This emphasized that reducing the height of one pier in an asymmetric configuration significantly increases its rotational demand during rocking, affecting the deck's rotation and uplift but not the longitudinal displacement demand of the bridge.Skinner et al. and Solberg et al. [15,16] investigated bridge pier design utilizing earthquake records from Loma Prieta, Northridge, and Kobe to ensure post-shock serviceability and minimize financial losses.Astaneh-Asl et al. [17] presented rocking foundations as an alternative to fixed-base systems, which enhanced seismic stability.
Gajan et al. [18][19][20] analyzed the response of rocking and fixed-base foundations to a realistic earthquake using numerical simulation.Deng et al. [21] conducted a study in which they utilized a damped elastic column, a rocking base, and a deck mass to compute the linear damping and period.Kourkoulis et al. [22] investigated rocking isolation as an alternative to over-designing frame footings based on an emerging seismic design approach.Gelagoti et al. [23] presented a simple yet realistic reinforced concrete moment-resisting frame with one bay and two stories to compare the two methods.Scattarreggia et al. [24][25][26] numerically examined damage propagation and progressive failure up to full collapse using the applied element method.Galvez et al. [27] found that the rigid-body rocking motion induced by non-reinforced masonry parapets and facades exceeded the static stability threshold before collapse.
Rocking isolation has been found to be a highly effective method for mitigating inertial forces.Anastasopoulos et al. and Loli et al. [28,29] investigated the rocking isolation design philosophy using appropriately scaled models of contemporary reinforced concrete bridge piers.Antonellis et al. [30] reported on the same tests carried out by the University of California at San Diego's Network for Earthquake Engineering Simulation using its large outdoor high-performance shake table.Tsatsis et al. [31] examined a square foundation for both heavily loaded and lightly loaded structures.Liu et al. and Madaschi et al. [32,33] provided test findings from cyclic loading tests as well as experiments for the dynamic phase of the test.Ko et al. [34] proposed unconnected piles beneath a foundation as a preferable solution to reduce persistent deformation caused by shaking.The models underwent an earthquake series on a dry sandbed.Pelekis et al. [35] conducted an experiment to identify the required range of rocking amplitudes for base isolation.
Rocking motion is increasingly being utilized to isolate structures from large stresses caused by earthquakes [36][37][38].Several authors have utilized the discrete element method (DEM), first introduced by Cundall et al. [39], to examine the rocking motion of rigid bodies [40][41][42].Recent studies with the application of the DEM technique include the prediction of the seismic behavior of multi-block tower structures [43].These studies have demonstrated that the DEM is a useful tool for studying a rigid structure in freestanding motion or subjected to base motion.However, in most studies, the parameters determining the model's features (primarily stiffness and mechanical damping) were determined by fitting numerical findings to actual tests.One of the common assumptions in employed DEM models of rigid blocks in the literature is that the contact between the block and any entity is defined at only one point (typically the corner point).To overcome these challenges, this study utilizes a DEM technique based on modeling the rigid block as a group of glued (clumped) particles that follow rigid body dynamics.The proposed model was used to simulate the physical experimental results of Pena et al. [44,45], where a series of experiments with accurate measurements for the freestanding motions of four rectangular blocks with different geometrical characteristics was reported.The advantage of the proposed model is that it inherently accounts for nonlinear behavior, whereby the clumped particles behave as a rigid body that will not break apart, regardless of the forces acting upon it.Additionally, the base of the block was modeled using smaller particle sizes to provide numerous contact points between the block and the base wall.In some sense, the proposed model is similar to having a series of nonlinear Winkler springs (with stiffness predicted using the Hertz contact law) at the contact points between the block particles and the base wall.Further, this study pioneers the use of the DEM to analyze bridge pier and deck interactions under seismic loads, building upon the validated single block model response.

Computational Model
The methodology utilized for the analytical approach employs a classical nonlinear formulation based on the fundamental principles of energy and momentum conservation.It involves deriving the equations of motion for a rigid body to characterize its response under dynamic loading conditions.These equations characteristically describe the planar rocking response of rigid bodies, as depicted in the validation plots by Zhang et al. [46] and Shenton III et al. [47].This process requires accounting for the geometry of the structure, material properties, and the critical angles at which instabilities such as overturning may occur.The analytical models provide insight into the system's natural frequencies and the coefficient of restitution, which are essential for understanding the dynamics of the structure during seismic events.Validation of the freestanding response also includes Pena's experimental data.The analytical method provides a comparable baseline for the test, allowing us to calibrate our DEM models, and then analyze the seismic resistance in the components of the structures.
The DEM models the rigid block as a collection of distinct spherical particles, as shown in Figure 1.The displacement of each particle is independent of the other particles, and interactions between particles only occur at contacts, allowing for the behavior of the system to be described in terms of the movement of each particle and the inter-particle forces acting at each contact point.The motions of each particle are calculated through an explicit time integration scheme, which operates under the principle that during a small timestep the change in velocities and accelerations is so small that they can be assumed to be constant within that timestep.
In the DEM, particles are assumed to be rigid but can overlap.The amount of overlap is dictated by a contact force model.A contact force, f c , between two particles consists of normal, f n , and shear, f s , components.The normal component was idealized using a linear spring model connected in parallel with a viscous dashpot: where u n and v n are the relative displacement and velocity at the contact along the line connecting the spherical particle centers, n is the unit normal vector at the contact, c n is the normal viscous damping coefficient, and k n is the normal contact stiffness.The incremental shear contact force, df s , was modeled using an elastic spring in series with a frictional slider: where du s and dv s are the incremental shear displacement and velocity vectors at the contact, c s is the shear viscous damping coefficient, and k s is the shear contact stiffness.The shear and normal forces are related by a slip Coulomb, as explained in the model by Itasca et al. [48], such that F s ≤ µF n ; where F s is the resultant shear force at the contact, µ is the Coulomb friction coefficient, and F n is the resultant normal force at the contact.The Hertz contact model is an advanced nonlinear contact force model that takes into account the deformation characteristics of elastic bodies in contact.It extends the classical Hertzian contact theory, which describes the contact between smooth, elastic spheres, to include frictional effects, based on the work of Mindlin et al. [49].In DEM simulations, this model is crucial for capturing force-displacement behavior during collisions and is implemented at both ball-ball and ball-facet contacts.The Hertz model, as implemented in the DEM simulations, enables the calculation of normal and shear forces by assuming a finite-sized contact area, which is negligible compared to the contacting bodies' radii.While the contact model produces no resistance to relative rotation, thus setting the contact moment to zero, it defines a force-displacement law that includes both Hertzian elastic force and viscous damping components, enhancing the simulation's realism in dynamic impact scenarios.
The relationship between force and displacement within the Hertzian contact model is articulated through the equation: where F contact and M contact denote the overall contact force and moment, respectively, F Hertz represents the force attributed to Hertzian mechanics, and F damp is the damping force, accounting for viscous damping.This model differentiates between engaged and disengaged contacts by evaluating the separation between surfaces, ensuring that force calculations are only performed when surfaces are in direct contact [48].
For the Hertzian force, the normal and tangential components are given by the following: with the shear force within the contact plane described by the following: The equation above defines the shear force, F H shear , in a three-dimensional context.The terms F H ss , F H st , and F H su correspond to the components of the shear force along the ŝ, t, and û directions, respectively.The vector û is orthogonal to both ŝ and t, thereby completing the representation of the force in three-dimensional space.This Hertzian contact framework is pivotal for simulating the interactions and mechanical responses in granular materials, collision phenomena, and particulate system behavior under various stress conditions.The implementation steps for the force-displacement law within the Hertz model are as follows.
Step 1: Normal Force Update.The normal component of the Hertzian force is updated with Here, h n is influenced by the geometric and material properties of the contact, defined as where G is the effective shear modulus, ν is Poisson's ratio, and R is the mean radius of the contact, computed from the radii R (1) and R (2) of the contacting bodies.
Step 2: Shear Force Calculation.The trial shear force in a three-dimensional context is computed as: where F H shear,0 represents the vector of the initial Hertz shear force components in 3D space and ∆θ θ θ s denotes the vector of the shear displacement increments along all three dimensions.The scalar k s symbolizes the initial tangent shear stiffness, which is evaluated based on the prevailing normal force.The updating of the shear force considers the shear strength, F µ s , to comply with the Mohr-Coulomb failure criterion, adapted for the multidimensional analysis.
Step 3: Dashpot Normal Force Update.The normal dashpot force is updated according to with F * defined as where k n is the normal tangent stiffness and δn is the relative normal velocity.
Step 4: Dashpot Shear Force Update.The dashpot shear force, F D shear , is updated based on the dashpot behavior mode, distinguishing between the full shear and slip-cut scenarios, with the force set to zero in the latter case when sliding occurs.

Model of a Rigid Body Formulation
The rigid body formulation provides a means to create and modify groups of slaved particles or rigid bodies.The model behaves as a rigid body (i.e., the particles comprising the rigid body remain at a fixed distance from each other).Contacts internal to the rigid body are skipped during the calculation cycle, resulting in a saving of computational time compared to a similar calculation in which all contacts are active.However, contacts with particles external to the rigid body are not affected (i.e., such contacts will develop when the particles comprising the boundary of a rigid body come into contact with other particles).Particles within a rigid body may overlap to any extent; contact forces are not generated between these particles, but any contact forces that exist when the rigid body is created or when a particle is added to the rigid body will be preserved unchanged during cycling.This model acts as a rigid body (with a deformable boundary) that will not break apart, regardless of the forces acting upon it.In this sense, the model differs from a group of particles that are bonded to one another [48].The basic mass properties of a rigid body are its total mass, m, location of the center of mass, x G , and moments and products of inertia, I ii and I ij .For a general rigid body comprised of N p balls, each of which has mass m [p] , radius R [p] , and centroid location x i , the mass properties are defined by the following equations:

Full Equations of Motion for a Rigid Body
The motion of a rigid body is determined by the resultant force and moment vectors acting upon it.Because the model is treated as a rigid body, its motion can be described in terms of the translational motion of a point in the rigid body and the rotational motion of the entire model.The translational motion of the center of mass is described in terms of its position, x i , velocity, ẋi , and acceleration, ẍi ; the rotational motion of the rigid body is described in terms of its angular velocity, ω i , and angular acceleration, ωi .The equations of motion can be expressed as two vector equations, one of which relates the resultant force to the translational motion, while the other relates the resultant moment to the rotational motion.The equation for translational motion can be written in the vector form: in which F i is the resultant force, the sum of all externally applied forces acting on the rigid body, m is the total mass of the rigid body, and g i is the body force acceleration vector arising from gravity loading.The resultant force is computed by: i is the externally applied force acting on particle (p) and is the force acting on particle (p) at contact (c).The equation for rotational motion can be written in the vector form by Ginsberg and Genin [50] : in which M i is the resultant moment about the center of mass and H i is the time rate-ofchange of the angular momentum of the rigid body.The resultant moment in Equation ( 17) is computed by i is the externally applied moment acting on particle (p), k is the resultant force acting on particle (p) at its centroid, and is the force acting on particle (p) at contact (c).
Equation ( 17) is referred to a local coordinate system that is attached to the rigid body at its center of mass.For such a system, the time rate-of-change of the angular momentum can be written as where The equations of motion, given by Equations ( 15) and ( 17), are integrated using a centered finite-difference procedure involving a timestep of ∆t.The quantities ẋi and ω i are computed at the mid-intervals of t + ∆t/2, while the quantities x i , ωi , F i , and M i are computed at the primary intervals of t + n∆t.The following finite-difference expressions describe the translational and rotational accelerations at time t in terms of the velocity values at mid-intervals.The accelerations are calculated as: Inserting the first of these expressions into Equation ( 15), and solving for the velocity at time (t + ∆t/2), results in Obtaining the angular velocity at time (t + ∆t/2) requires solving a nonlinear set of equations using an iterative procedure during each timestep.Substituting Equation (19) into Equation ( 17) allows one to write which can be expressed in matrix form as with: Equation (24) provides three equations for the six unknowns, ω i and α i .These six unknowns are found by using the following iterative procedure.

Set ω
[0] i equal to the initial angular velocity (i.e., before the motion computation).

3.
Solve Equation ( 24) for α Determine a new (provisional) angular velocity, assuming no damping: 5. Revise the estimate of ω i as Set n := n + 1 and go to Step 3.
Numerical experimentation has shown that this solution converges after four iterations.
Thus, after four iterations have been performed, we set ω i = ω [4] i .Note that the above iterative procedure breaks if one or more rotational degrees-of-freedom are fixed.For such cases, the approximate solution method is employed.The same local nonviscous damping formulation used to dissipate energy from each individual particle is also used to dissipate energy from each rigid body.The damping coefficient for a rigid body is the arithmetic mean of the coefficients of its component particles, and the angular velocity used in the damping formulation is taken to be ω[∆t] i .After determining the translational and rotational velocities of the rigid body, the position of the rigid body model center of mass is updated by and the velocity of each ball in the rigid body is set to The position of each ball in the rigid body is then updated based on the ball velocity according to Equation (27).Equation (28) uses position vectors of the particles and the rigid body that exist half a timestep behind the calculated velocity vector.This causes a gradual drift in the relative locations of particles within a rigid body: the error is related to the magnitude of the spin components.For rapidly spinning rigid bodies, it is possible to obtain a more accurate result by using time-centered position vectors in the calculation of particle velocities.We replace Equation (28) with the following expression. where: and x [p] i denotes the position vector of particle [p] before the update.Expanding Equation ( 29) and substituting from Equation (30), we obtain: Equations ( 32)-( 34) can be expressed in matrix form as where: Equation ( 35) is solved by LU decomposition for the velocity vector of each particle, ẋ[p] i , when the full mode of calculation is selected in three dimensions.In two dimensions, the time-centered scheme is always used, and Equation (35) reduces to Equations ( 36) and (37), as follows: where: .

Validation
Pena's experimental work serves as a critical benchmark for validating the contact models employed in the DEM simulations.This work is particularly vital for establishing a baseline for the accuracy of microscale interactions within our proposed framework.The configuration of the experimental work was mirrored in the performed DEM simulations, leveraging the physical parameters outlined in Table 1 to validate the DEM model against Pena's experimental findings and replicate its dynamic behavior.This focused approach allows for an in-depth validation of the DEM simulation properties by comparing them with the experimental results obtained from the shake table tests.The model was configured to reflect the base loading scenarios encountered in the experimental setup, thus facilitating a direct correlation between the numerical simulations and physical experiments.The construction of the base wall in the DEM model, representing the shake table from Pena's experiments, involved defining its properties, such as stiffness, friction, and density, to ensure fidelity of the experimental conditions (as shown in Table 2).Additionally, the built-in energy-tracking function of PFC was utilized to monitor energy loss during the simulations, a critical aspect of dynamic systems.To emulate the physical experiment, the DEM model of the block was discretized into three sections with varying ball radii.This granularity was necessary to accurately replicate the contact dynamics at the base, facilitate computational efficiency, and maintain symmetry to reliably reflect the block's mass and moment of inertia.A parametric study was conducted to optimize the ball radius, position, and volume, ensuring sufficient overlap and contact at the base, as well as to clump the particles effectively into a rigid body (Figure 2).The mechanical damping in the model was calibrated to simulate the energy dissipation observed in the physical experiments.Comparative plots of the DEM simulation results and experimental data were generated using MATLAB, demonstrating the ability of the DEM to replicate the physical experiments.The stress evolution at the base, as the block underwent its rocking response, was captured and analyzed to provide insights into the contact stress dynamics.
Validation of the free-standing rigid block model is a critical step in ensuring the accuracy of the DEM simulations.Figure 3 shows a comparison of the rocking response over time for Pena's experimental block.The plot constitutes the Hertz contact model, analytical prediction [46], and Pena's experimental data; the plot illustrates the efficacy of each approach in capturing the dynamic behavior of the block.The analytical model shows a slight phase difference, indicating the potential omission of intricate interactions that are otherwise accounted for in the DEM simulations and physical experiments.The DEM simulation results exhibit a high degree of agreement with Pena's experimental outcomes, confirming the model's capacity to capture the complex dynamics involved in freestanding motion.The close agreement between the DEM results, analytical response, and the experimental data strengthens confidence in the DEM as a robust tool for seismic analysis-a pivotal step in the validation process for the microscale framework developed in this research.

Contact Mechanics in a Free-Standing Rigid Block
Figure 4 presents the evolution of the total impact contact force ratio over a time period of 6 s.This ratio is defined as the ratio of the total contact force to the static weight of the block while it undergoes free-rocking motion.The force ratio starts with an initial high ratio of approximately 150, indicating that the first impact force at contact is two orders of magnitude higher than the static weight of the block, suggesting a high impact at the beginning of the rocking motion.It also shows a steady decline in the force ratio from the initial peak, decreasing to lower values as time progresses over 6 seconds.The decreasing magnitude of these peaks over time suggests that the block is gradually losing energy, leading to smaller oscillations until the rigid block settles at a minimal force ratio level toward the end of the observed period.
Figures 5 and 6 depict the contact force ratio at the left and right corners, respectively, of the rigid block undergoing its free-rocking motion, from the DEM simulation.Figure 5 (left corner) shows a series of spikes in the force ratio, beginning sharply at approximately 2.8 and decreasing progressively with each oscillation until the rigid block stabilizes around the 6-s mark.The peaks are sharp and closely spaced initially but spread out and diminish in magnitude over time.Figure 6, similar to Figure 5, also displays a series of diminishing force ratio peaks over time.The pattern starts with a peak at approximately 3, decreasing progressively.The force ratio behavior at the right corner mirrors that of the left, indicating the symmetrical rocking behavior of the block.
Figure 7 presents the temporal evolution of force ratio distribution along the width (x) of the block at peak rocking angles (θ y ).Each sub-figure corresponds to different rocking angles ranging from θ y = 2.9 • to θ y = 1.5 • .The sub-figures correspond to the instants of peaks in the rocking angle.At these instants, the corner edge force is provided by two or three contact points whose sum is equal to one.This observation verifies that static equilibrium is maintained at those instants.

Microscale Energy Monitoring and Dissipation Mechanism
Individual energy components were continuously monitored during the simulation.In this work, the different energy components can be separated into two categories; an input energy and a dissipation energy.The main source of input energy is the body work, E b , resulting from the gravitational force of the block.This energy is then dissipated through inter-particle sliding (friction), E f , damping energy, E d , kinetic energy, E k , and elastic strain energy at the contacts, E s .Energy balance was verified for all simulations.That is, the input energy is equal to the dissipated energy: The figures presented (Figures 8 and 9) provide a quantitative analysis of the energy dynamics associated with the rigid block under free-rocking motion.Figure 8 illustrates the distribution of energy components over time: body work (potential), damping, kinetic, frictional, and strain energies.Notably, kinetic energy displays sharp peaks of approximately 27 J, aligning closely with peaks in potential energy.The peaks occur at the time of impact.The damping energy shows a step-like behavior with the steps corresponding to the sudden increase in damping energy upon impact.Frictional and strain energies remain minimal throughout the motion, with frictional energy reaching 3 J, indicating negligible losses due to friction.In Figure 9, the energy balance plot verifies the equity demonstrated by Equation (38).

Seismic Isolation of Rocking Bridge Pier Structures
The notion of 'rocking bridge piers' has emerged as an exceptional design philosophy, utilizing the dynamic interplay between structural inertia and base flexibility to mitigate seismic forces.This approach is characterized by allowing controlled rotational movements around the base of the pier, which can significantly reduce seismic demand on the superstructure.This section focus on the seismic isolation mechanisms inherent in rocking bridge pier designs, supported by a series of computational scenarios.Each scenario progressively builds on the preceding one, beginning with the fundamental single-pier model and culminating in a complex multi-pier system with enhanced mass and frictional properties.Figure 10 graphically encapsulates the comparative stability across the examined cases.The seismic performance of a pier or deck is evaluated in terms of tilt angle, which is derived from the lateral displacement, known as drift, between two points of the structure.The permissible tilt is derived from drift ratios stipulated in design codes.
Given the drift ratio Drift = Lateral displacement at the top of the pier Height of the Pier (39) Drift = Lateral displacement at one end of the deck Length of the Deck (40) the tilt angle, θ, can be inferred from the drift as The American Association of State Highway and Transportation Officials (AASHTO) [51] prescribes a drift limit of 2.5%, translating to a tilt angle of approximately 1.43 • , ensuring that piers can endure lateral shifts while maintaining stability.Deck structures are also subject to drift limitations, with AASHTO suggesting a more conservative tilt angle of approximately 0.86 • , acknowledging the different roles and behaviors of decks compared to piers in seismic response.

Seismic Response of a Single-Pier
A single pier-length, 5 m; width, 1.6 m; and height, 9.6 m-is positioned on a base wall and subjected to sinusoidal motion with a maximum amplitude of 0.5 g, as shown in Figure 11 and Equation (42).This motion, defined by a precise velocity profile, persists for 16 s and gradually diminishes to zero over a period from 16 to 20 s.For the model's comprehensive representation, the pier is segmented into three distinct zones-bottom, middle, and top (as seen in the case of the experimental work validation).Each exhibits unique geometric characteristics, with the top and bottom sections modeled with an 8 cm particle radius and the middle section with a larger radius of 20 cm.An overlapping ratio of 0.75 was used between the particles.The properties of the single-pier simulation are presented in Table 3.The pier was subjected to a seismic motion by shaking the base wall using the following velocity signal: where the coefficient A is defined as The parameters are defined as follows: L1 = 14.5 (end of constant amplitude motion), L2 = 16 (time when sinusoidal motion stops), f = 1 (frequency), g = 9.81 m/s 2 (acceleration due to gravity), T = 4.5 (duration of increasing amplitude), ω = 2π f (angular frequency), b = 0.5 (scaling factor for amplitude), and A = amplitude of the motion.Overturning Boundary, =9.46 °Figure 12. Analytical [46,47] and DEM validation of the seismic response for a single-pier structure.

Rocking Isolation in Pier-Deck Systems under Low-Intensity Seismic Excitation
This section examine the effectiveness of rocking isolation in reducing seismic responses when subjected to low-intensity seismic excitations, specifically a 0.1g input motion.This study includes extensive time histories of different structural response parameters, such as acceleration, angular velocity, rocking angle, displacements, and the resultant forces at the contact interface between each pier and the base, and between the piers and the deck.
The input motion is simulated with the response of the pier and deck system recorded as acceleration, angular velocity, and rocking angles.The system is subjected to a lower peak ground acceleration (PGA) of 0.1 g with a frequency of 1Hz.The input acceleration of 0.1 g is nearly matched by the response of the bridge piers and deck, with only a marginal reduction in peak acceleration values (see Figure 13).This suggests that, at a lower PGA of 0.1g, the rocking mechanism's effectiveness is less pronounced than at higher intensities.Figure 13 shows that the acceleration response for piers 1 through 4 closely tracks the input motion, with peak values at or slightly below 0.1 g.That is, the motion of all the system components (piers and deck) were in-phase.This indicates a reduction ranging from negligible to a maximum of 10%.The bridge deck shows a similar trend with a slight reduction in peak acceleration, which suggests that the rocking isolation system is less engaged at this lower level of input motion.
Given that the maximum amplitude for üg = 0.1g < g tan α = 0.167 g, the bridge piers experienced a less pronounce rocking motion, with the pier-deck system having the same input motion, supporting the horizontal acceleration in Figure 13.The vertical acceleration in Figure 15 shows different dynamic responses from each pier and the bridge deck under the same seismic conditions.The magnitude of the vertical acceleration, however, was negligible.There was no significant rocking of the piers under this low amplitude input motion, as shown by the amplitudes of the angular velocity (Figure 16) and rocking angle (Figure 17).

Contact Mechanics in Pier-Deck System under Low-Intensity Seismic Excitation
In order to investigate the contact behavior between the piers and base wall as well as the piers and the deck, the contact forces at the interfaces in the vertical direction and shaking lateral direction were monitored during the course of the simulation.These forces were used to compute the resultant horizontal force and vertical force at the contact between the piers and base wall as well as those between the deck bottom surface and the top of the piers.Figure 18 displays the horizontal resultant force at the interface between the piers and the base wall under a 0.1 g and 1 Hz seismic loading.The resultant forces for all the piers peak at approximately 120 kN.All piers display similar patterns to Pier 1 but with slightly different amplitudes, indicating that the response of the piers was perfectly in-phase.Figure 19 illustrates the normalized vertical resultant forces.In this Figure, the total reaction force at each pier-wall interface was normalized by the weight of the pier plus one quarter the weight of the deck.As shown in Figure 19, the normalized reaction values vary around baseline of 1.00, indicating slight deviations in force application during the seismic loading.
Figure 20 presents the normalized eccentricity at the interface between the piers and the base wall.The computed eccentricity for every pier was normalized by the pier width.The normalized eccentricity, therefore, measures the offset of the resultant force from the center of the pier-deck system to assess full contact between the piers and the base wall.The eccentricity boundary ranged from −0.17 to 0.17 (which corresponds to B/6), with all four piers exhibiting eccentricities within the boundary, and oscillating around the neutral axis, at a range from −0.1 to 0.1.This suggests that the contact surface between the piers and base wall remained fully effective and there was no uplift between the piers and the base wall since the resultant force always remained within the middle one third of the pier.Figure 21 shows the frictional ratio at the interface between the piers and base wall.The frictional ratio, indicating the ratio of lateral force resisted through friction to the normal force, remains significantly below the friction coefficient (µ) represented by the dashed line at 0.59.For all piers, the frictional ratio oscillates just above zero, peaking slightly but never approaching the static friction limit.Figure 22 displays the resultant horizontal contact force at the interface between the piers and deck under a 0.1 g and 1Hz seismic loading.The forces for all four piers fluctuate between approximately −30 kN and 30 kN, following a somewhat sinusoidal pattern with peak activity between 10 and 20 s. Figure 23 shows the normalized vertical resultant forces, which are clustered close to 1.00, with slight deviations up to approximately 1.025 or down to approximately 0.99.Note that the vertical reaction at each pier-deck interface was normalized by one quarter the weight of the deck.Figure 24 presents the normalized eccentricity of the vertical reaction at the interface between the pier and the deck, fluctuating at approximately −0.05 and 0.05, staying within the boundary of the eccentricity between −0.17 and 0.17.Again, this suggests that there was no uplift between the deck and the piers and the contact surface at the interface between the deck and the piers was fully effective throughout the simulation.
Figure 25 shows the frictional ratio remains well below the static friction coefficient (µ) indicated by the dashed line at 0.59.The frictional ratios are very close to zero, suggesting minimal relative movement between the deck and the piers.This mean that sliding potential at the interface remains consistently below 0.1 for all piers throughout the simulation.This is well below the threshold of 0.59, suggesting that, under the current seismic conditions, there is minimal risk of sliding at the interface, indicating good frictional resistance.

Rocking Isolation in Pier-Deck Systems under High-Intensity Excitation
The same bridge pier-deck system shown in Figure 14 was subjected to a more intense input motion applied to the base wall.The input motion had a maximum amplitude of 0.5 g and a frequency of 1 Hz.A significantly different response was observed in this case compared to the low-intensity seismic loading discussed in the previous section.The horizontal acceleration transmitted from the base wall to the pier-deck system was smaller in amplitude than that of the input motion (Figure 26).The peak accelerations of the piers and deck are lower than the input motion, with a reduction of approximately 60-70% relative to the input motion.In addition, there was a prolonged period of free vibration post shaking that lasted for approximately 20 s after the base wall stopped shaking the system.There was significant rocking in this case, as shown by the angular velocity time histories in Figure 27 that continued post shaking.Figure 28 shows the rocking responses of the four bridge piers under sinusoidal seismic excitation.Close inspection of the response of the piers shows a slight phase shift in their motion.The figures also show the rocking motion of pier 1 was accumulating opposite of pier 4 and pier 2 rocking opposite of piers 3. Note that the interface surface between the piers is not flat but rather corrugated due to the nature of the pier and deck model representation by a collection of particles.The implication of having a non-flat surface at the interface is that the contact points are not perfectly similar between the top of the piers and the bottom deck.This dissimilarity becomes more pronounced once the piers start to tilt, and new contacts are formed and other contacts are deleted.In addition, the deck merely sits on top of the piers, with only frictional resistance acting at the interface between the deck and the piers.This means that the deck is free to translate and rotate at the interface between the bottom of the deck and the top of the piers.The response of the rocking frame obtained by the DEM simulation was compared against the analytical approach presented by Makris and Vassiliou [52], which is shown in Figure 29.It can be seen that the obtained amplitudes of angular velocities (Figure 27) and rocking angles (Figure 28) for the piers using DEM are within a reasonable agreement with the analytical approach.Figure 30 presents the lateral displacement of the piers and the deck.The cyclic amplitude of the deck displacement is smaller than that of the piers (in agreement with the lateral acceleration time history).The lateral displacements of the piers and deck show that sliding took place during the simulation.By the end of the simulation, the deck displaced about 6 cm.The piers displaced at a different magnitude and direction.Piers 1 and 4 experienced an absolute final displacement of about 6 cm, while piers 2 and 3 experienced an absolute displacement of about 1 cm.The vertical displacement in Figure 31 shows that the deck experienced a peak vertical displacement of about 40 mm during shaking, which dropped to about 15 mm once shaking stopped.The vertical displacement of the deck stabilized to zero after about 60 s.

Contact Mechanics in Pier-Deck System under High-Intensity Excitation
Figure 32 shows significant fluctuations in horizontal resultant forces between the piers and the base wall, reaching up to an absolute magnitude of about 1000 kN for all piers under the 0.5g and 1Hz seismic loading.The forces fluctuate widely after 15 s of shaking and during the free vibration phase of the response.The normalized vertical resultants at the interface between the piers and the base wall are shown in Figure 33.The highest magnitude in oscillation was experienced by piers 2 and 3, with ratios reaching values as high as 2.2 and as low as 0.5.Figure 34 shows the horizontal resultant forces between the piers and the deck.Piers 2 and 3 experienced the largest magnitudes of horizontal force post shaking.The normalized resultant forces at the interface between the piers and the deck fluctuated around the baseline values and show multiple peaks exceeding a factor of 1, with values spiking up to approximately 3.2 (Figure 35).The pattern of the spikes is not uniform among the piers, indicating diverse responses to the seismic load.The normalized resultant represent the vertical load relative to the static weight of the pier-deck system, highlighting the dynamic loading effects.Peaks above 1 indicate that the dynamic loads exceed the normal static loads and separation between the deck and piers occurred at points of zero normalized resultant forces.

Conclusions
This study presented a comprehensive analyses on the dynamics of rigid bodies under seismic loading, focusing on their rocking response and energy dissipation mechanisms.The study introduced an innovative microscale framework for seismic stability analysis of bridge pier rocking isolation, employing the DEM along with the Hertz contact model formulation.The advantage of the presented DEM model is that it accounts for several contact points at the impact interface, with the stiffness of those contact points being nonlinear elastic springs.The comparative analysis between DEM simulation results and Pena's experimental data exhibited a high degree of agreement, establishing the DEM's reliability in predicting the dynamics involved in the seismic behavior of free-standing rigid blocks.The presented validated model offers new insights into the microscale energy interactions and stress distributions, which are important for designing structures capable of withstanding seismic forces.
Results from the performed simulations of the seismic response of the pier-deck systems showed that, under low excitation input motion, there was no significant rocking and therefore there was a reduction in the magnitude of input acceleration transmitted to the piers and deck.In this case, there was no uplift nor sliding between the elements of the pier-deck system.The piers maintained full contact with the base wall as well as with the deck, with the vertical resultant forces falling within the middle third of the pier width.The vertical resultant forces of each pier remained virtually equal and within the magnitude of the initial static resultant force.Upon subjecting the pier-deck system to a strong ground motion, rocking of the piers was observed.The work further confirms the rocking isolation hypothesis as an effective seismic protection mechanism for bridge structures, with results showing a substantial reduction of seismic loads across a range of loading conditions, effectively reducing input acceleration by 70% in the piers and 60% in the deck during higher intensity events.However, uplift and sliding was observed during the course of the high intensity ground motion simulation.The vertical resultant force of each pier varied significantly and instantaneously vanished due to separation between the deck and the piers at contact points.

Figure 1 .
Figure 1.Schematic of the constitutive law of the normal and shear contact forces between two particles.

Figure 2 .
Figure 2. The DEM model setup for Pena's experimental block (blue) with base wall (gray).

Figure 3 .
Figure3.Validation of the DEM simulation, analytical approach[46], and experimental data of a rigid block in freestanding rocking motion.

Figure 4 .Figure 5 .
Figure 4. Total impact contact force ratio from DEM simulation of the experimental rigid block in free-standing rocking motion.

Figure 6 .Figure 7 .
Figure 6.Contact force ratio at the right corner of the experimental rigid block undergoing free-rocking motion from DEM simulation.

Figure 8 .Figure 9 .
Figure 8. Energy component for the rocking motion of Block 1

Figure 11 .
Figure 11.DEM model setup for single pier with base wall.

Figure 12 illustrates
Figure12illustrates the dynamic behavior of a single-pier structure under sinusoidal ground motion, juxtaposing theoretical instability predictions from Zhang et al.[46] and Shenton III et al.[47] with DEM simulations.Both the analytical and DEM approaches closely track each other, affirming the validity of the proposed DEM model in capturing the seismic response.The rocking angle is normalized with the single pier slenderness, α, where α = arctan b h .

Figure 13 .
Figure 13.Horizontal acceleration plot of multi-pier and deck system under 0.1 g, 1 Hz.

Figure 14 .
Figure 14.DEM model of a bridge pier-deck system.

Figure 18 .Figure 19 . 6 Figure 20 .Figure 21 .
Figure 18.Horizontal resultant force at the interface between the piers and base wall under 0.1 g and 1 Hz.

Figure 22 .Figure 23 . 6 Figure 24 .Figure 25 .
Figure 22.Resultant horizontal contact force at the interface between the piers and deck under 0.1 g and 1 Hz.

Figure 26 .
Figure 26.Horizontal acceleration plot of multi-pier and deck system under 0.5 g and 1 Hz.

Figure 27 .
Figure 27.Normalized angular velocity plot of multi-pier and deck system under 0.5 g and 1 Hz.

Figure 29 .
Figure 29.Analytical response of rocking pier-deck system under 0.5 g and 1 Hz.

Figure 30 .
Figure 30.Lateral displacement of the piers and deck under 0.5 g and 1 Hz.

Figure 31 .
Figure 31.Vertical displacement on the bridge deck in pier-deck system under 0.5 g and 1 Hz.

Figure 32 .
Figure 32.Horizontal contact force at the interface between the piers and base wall under scenario 1 (0.5 g, 1 Hz).

Figure 33 .
Figure 33.Normalized vertical reaction at the interface between the piers and base wall under 0.5 g and 1 Hz.

Figure 34 .Figure 35 .
Figure 34.Horizontal contact force at the interface between the piers and deck under 0.5 g and 1 Hz.

Table 1 .
Physical parameters and geometric size for Pena experiments.

Table 2 .
DEM simulation properties for Pena's experiments.

Table 3 .
DEM simulation properties for the single-pier configuration.