INTRODUCTION
The nonlinear model of impacts based on the Hertz contact law has been used in numerous studies of earthquakeinduced structural pounding (e.g., Jing and Young, 1991; Pantelides and Ma, 1998; Chau and Wei, 2001). It simulates the relation between the pounding force and deformation during impact more realistically than the linear models (Goldsmith, 1960). However, the application of the Hertz elastic spring alone prevents us from simulation of energy dissipation during impact.
Lankarani and Nikravesh (1990) presented an improved version of the Hertz contact model based on the Hertz contact law, where a nonlinear damper is used along with the nonlinear spring, to represent the dissipated energy in impact.
Marhefka and Orin (1999) developed a compliant contact model with a nonlinear spring in parallel with a nonlinear damping, in order to overcome problems in the use of rigid body models with Coulomb friction and to eliminate the tension forces arising in the KelvinVoigt model.
Muthukumar and DesRoches (2006) and Muthukumar (2003) used the Hertzdamp model, which is equivalent to the improved version of the Hertz contact law model (Lankarani and Nikravesh, 1990; Marhefka and Orin, 1999), to study pounding simulation in structural engineering. Two Single Degree Of Freedom (SDOF) systems with different period ratios and a suit of 27 ground motion records of different Peak Ground Acceleration (PGA) levels were used to compare the Hertzdamp, linear spring, Kelvin, Hertz and stereomechanical models. Numerical results indicate that the Hertz model provides adequate results at low PGA levels, while the Hertzdamp model is recommended at moderate and high PGA levels (Muthukumar and DesRoches, 2006).
Jankowski (2005b) proposed a nonlinear viscoelastic model for more precise simulation of structural pounding based on the Hertz law of contact. This model gave the smallest simulation errors in the response time histories of the structural pounding experiments considered in comparison with the Hertz model and the linear viscoelastic model introduced by Anagnostopoulos (1988, 2004). Additionally, this model and the linear viscoelastic model gave the smallest simulation errors in the pounding force time histories during contact.
The nonlinear viscoelastic model was used to investigate earthquake induced
pounding between two buildings (Jankowski, 2008) as well as to determine the
pounding force response spectra under ground motion excitation (Jankowski, 2005a,
2006b). Numerical results confirmed the ability of the nonlinear viscoelastic
model to accurately simulate seismic pounding.
Two different models have been considered as improved versions of the Hertz model by adding a nonlinear damper. However, no comprehensive comparison has been given so far. The aim of this study is to compare the accuracy of the nonlinear viscoelastic model and the Hertzdamp model in seismic pounding analysis. The first comparison uses the results of two impact experiments as well as the results of shaking table experiments on pounding between two steel towers excited by harmonic waves. The second comparison applies a suite of thirty ground motion records from thirteen different earthquakes applied to pounding between two SDOF systems with varying period ratios.
Fung (2001, 2002) used the differential quadrature method to solve both dynamic problems governed by secondorder ordinary differential equations in time and first order initial value problems. Civalek (2007) used the method of Harmonic Differential Quadrature (HDQ) for the nonlinear dynamic response of Multi Degree Of Freedom (MDOF) systems. In this study, we use Implicit RungeKutta (IRK) methods to solve the system of first order ordinary differential equations.
HERTZ CONTACT LAW FORCEBASED MODELS
Hertz model: A nonlinear spring stiffness, depending on elastic properties
of the colliding structures, is used in the Hertz model to simulate structural
pounding. The impact force between the two structures of masses m_{i}
(i = 1,2) follows the relation (Goldsmith, 1960; Lankarani and Nikravesh, 1990,
1994).
where, δ (t) is the relative displacement. Assuming that the colliding
structures are spherical of density ρ and the radius R_{i} estimation
can be calculated through the following relation (Goldsmith, 1960):
The nonlinear spring stiffness k_{h} is linked to the material
properties and the radii of the colliding structures as stated through the following
formula (Goldsmith, 1960):
where, h_{1} and h_{2} are the material parameters defined
by the formula (Goldsmith, 1960):
Here, γ_{i} and E_{i} are Poisson's ratio and Young's
modulus, respectively. When R_{2}→∞, i.e., the second structure
becomes a massive plane surface, the nonlinear spring stiffness k_{h}
is defined as (Goldsmith, 1960):
Since the Hertz model is fully elastic (Eq. 1), it does not
allow us to consider the energy dissipation during the collision.
Hertzdamp model: Hertzdamp model has been considered by Muthukumar and
DesRoches (2006) and Muthukumar (2003) to study the pounding phenomenon in the
field of structural engineering. The energy loss during impact has been taken
into account by adding nonlinear damping to the Hertz model. The pounding force
is written as (Lankarani and Nikravesh, 1990):
where, e is the coefficient of restitution (Goldsmith, 1960),
is the relative velocity during contact and v_{1}–v_{2}
is the relative approaching velocity prior to contact which can be expressed
in terms of nonlinear spring and maximum deformation (Lankarani and Nikravesh,
1990):
Nonlinear viscoelastic model: Another improved version of the Hertz
model has been introduced by Jankowski (2005b) by connecting a nonlinear damper
in unison with the nonlinear spring. The contact force for this model is expressed
as:
where, is the impact stiffness parameter and
is the impact element's damping. Here is
an impact damping ratio corresponding to a coefficient of restitution e which
can be defined as Jankowski (2006a):
MATERIALS AND METHODS
To investigate the performance of the Hertzdamp and the nonlinear viscoelastic
models in capturing pounding, three different procedures of comparison are used.
The first one is based on two impact experiments conducted for different types
of structural members with various materials and contact surface geometries.
The second procedure uses the shaking table experiments on pounding between
two steel towers. The accuracy of each model in the first and second procedure
is assessed by calculating the Normalized Error (NE) to indicate the difference
between the experimental and numerical results (Jankowski, 2005b).
where, F is the response time history obtained experimentally, is
the response time history obtained numerically and is
the Euclidean norm. in
case of the time histories given in a discrete form can be calculated as:
where, n is a number of values in the time history record.
The third procedure is based on simulation using thirty ground motion records
from thirteen different earthquakes of Peak Ground Acceleration (PGA) levels
ranging from 0.1 to 1 (Muthukumar and DesRoches, 2006; Muthukumar, 2003). The
ground motion records were selected from PEER Strong Motion Database (http://www.peer.berkely.edu/smcat/).
Two Single Degree Of Freedom (SDOF) systems of structures with equal masses
and three different period ratios are subjected to the ground motion records.
Two cases are considered. In the first case, the initial separation distance
is chosen such that no pounding takes place. In the second case, the initial
separation distance is chosen such that pounding occurs. The performance of
the models is evaluated by comparing the numerically obtained displacement and
acceleration amplifications to each other.
COMPARISON BASED ON CONDUCTED IMPACT EXPERIMENTS
The results of two impact experiments are considered in this study.
Steeltosteel impact: Goland et al. (1955) carried out an experiment
to measure load time histories and strain propagation in a square beam of different
dimensions subjected to sharp lateral impacts letting a steel ball with diameter
ranging from to inch, drops onto the top of the beam from a specified height.
A force gauge was used to measure the forcetime history exerted by the ball
on the beam. The strain gauges were placed at different locations on the beam
to record the strain time histories at those locations. The dynamic equation
of motion for pounding between a ball of mass m_{1} dropping
onto a beam can be written by drawing the free body diagram (Jankowski, 2005b)
as shown in Fig. 1:
where, g
and F(t) denote the acceleration, gravitational acceleration and the pounding
force, respectively. The pounding force follows the relations (6) or (8). In
the experiment the maximum pounding force was 80.7 N when a ball of diameter inch fell from a height of 2 inches (Jankowski, 2005b). As for the stiffness
parameter in the nonlinear viscoelastic model, we set =
1.03 10^{10} N/m^{3/2} which was determined through an
iterative procedure in order to keep the maximum pounding force in the numerical
analysis and experiment to be the same (Jankowski, 2005b).

Fig. 1: 
Free body diagram of a ball dropping onto a beam 

Fig. 2: 
Pounding force time histories during impact between falling
steel ball and a steel hemisphere mounted on a beam 
The same procedure is used to obtain the stiffness parameter of the Hertzdamp
model, which is found to be k_{h} = 0.60 10^{10}
N/m^{3/2}. For the two models, the coefficient of restitution e = 0.6
is used. The results from the numerical analysis and the experiment are shown
in Fig. 2. Using Eq. 11 the normalized errors for pounding
force histories are found to be equal to 21.45% for the nonlinear viscoelastic
model and 24.18% for the Hertzdamp model.
Concretetoconcrete impact: Van Mier et al. (1991) carried out
an experiment on collisions between a prestressed concrete pile and a concrete
striker. The dynamic equation of motion for pounding between a striker of mass
m_{1} and a prestressed fixed pile can be written by drawing the free
body diagram as shown in Fig. 3 (Jankowski, 2005b).
The pounding force F (t) follows the relation (6) for the Hertzdamp model and
relation (8) for the nonlinear viscoelastic model. The stiffness parameter used
for the nonlinear viscoelastic model is
= 2.75x10^{9} N/m^{3/2 }(Jankowski, 2005b). For the Hertzdmap
model, k_{h} = 1.54x10^{9} N/m^{3/2} was found through
an iterative procedure in order to keep the maximum pounding force in the numerical
analysis equals the maximum pounding force of the experiment as 102.5 N. For
the two models, the coefficient of restitution is e = 0.6. The results
from the numerical computations and the experiment are shown in Fig.
4. It is found that the normalized error equals 25.04% for the nonlinear
viscoelastic model and 54.01% for the Hertzdamp model.

Fig. 3: 
Free body diagram for pounding between a concrete pendulum
striker and prestressed concrete pile 

Fig. 4: 
Pounding force time histories during impact between a concrete
pendulum striker and prestressed concrete pile 
COMPARISON BASED ON SHAKING TABLE EXPERIMENTS
Chau et al. (2003) carried out shaking table tests to investigate the
pounding phenomenon between two steel towers of different natural frequencies
and damping ratios subject to different combinations of standof distance and
seismic excitations. The two Single Degree Of Freedom (SDOF) building systems
shown in Fig. 5 were used in the numerical simulation as a
representative for the two towers. For i =1, 2, let m_{i} be the masses,
c_{i} be the viscous damping coefficients and k_{i} be the stiffness
for SDOF 1 and SDOF 2, accordingly. The coupling equation of motion for two
adjacent buildings subjected to horizontal ground motion has
the following form:

Fig. 5: 
Model idealization of adjacent structures 
where, and represent
the displacement, velocity and acceleration of the system, respectively. The
pounding force F (t) follows the relation (6) for Hertzdamp model and relation
(8) for nonlinear viscoelastic model. The values of structural stiffness and
damping coefficients: k_{i}, c_{i} can be calculated from the
formulas (Harris and Piersol, 2002).
where, T_{i}, ξ_{i} (i = 1, 2) denote the natural
structural vibration period and structural damping, respectively. We solve the
equation of motion (15) for the harmonic waves as excitations and the two different
models of pounding force considered herein. Chau et al. (2003) performed
a series of shaking table experiments for n = 1 (1 pounding per cycle) and n
= 2 (1 pounding per 2 cycles) using various harmonic waves as excitations for
the two SDOF systems shown in Fig. 5. Experimental, numerical
and analytical results of the relative impact velocity versus the excitation
frequency had been observed and computed for the comparison purposes between
experiments and theories. In this study, the results of two different sinusoidal
excitations as input are presented. We consider first the input shaking table ü_{g}(t) = 2.6 (2πf_{g}t) as an excitation acting
on SDOF 1 and SDOF 2 with the following dynamic characteristics m_{1 }
= 98.0 kg, f_{1 } = 5.04 Hz, ξ_{1} = 7.2%,
m_{2} = 146.4 kg, f_{2 }= 2.76 Hz, ξ_{2
}= 1.5% and a separation distance of 19.8 mm (Chau et al.,
2003). The obtained numerical solutions for the steady state relative impact
velocity using both the nonlinear viscoelastic model and the Hertzdamp model
are compared with the experimental results which are the average of 10 cycles
in the steady state (Chau et al., 2003). The results from the numerical
analysis and the experiment are shown in Fig. 6. The simulation
errors for the obtained numerical impact velocity compared with the experimental
results for n = 1 and n = 2 are 25.46 and 47.25% for the Hertzdamp
model and 31.66 and 51.26% for the nonlinear viscoelastic model.

Fig. 6: 
The steady state relative impact velocity versus the excitation
frequency. Numerical results obtained by nonlinear viscoelastic and Hertzdamp
models. The input shaking table is ü_{g}(t) = 2.6sin(2πf_{g}t) 

Fig. 7: 
The steady state relative impact velocity versus the excitation
frequency. Numerical results obtained by nonlinear viscoelastic and Hertz
damp models. The input shaking table is ü_{g}(t) = 1.9sin (2πf_{g}t) 
As a second input shaking table, we consider the harmonic excitation ü_{g}(t) = 1.9sin(2πf_{g}t)applied
to the SDOF 1 and SDOF 2 with the same dynamic characteristics and separation
distance as used before. The results from the numerical analysis and the experiment
are shown in Fig. 7. The simulation errors for the obtained
numerical impact velocity compared with the experimental results for n = 1 and
n = 2 are 24.97 and 54.25% for the Hertzdamp model and 33.75
and 59.80% for the nonlinear viscoelastic model.
COMPARISON BASED ON GROUND MOTION RECORDS
Muthukumar and DesRoches (2006) and Muthukumar (2003) assessed the performance
of the Hertzdamp model by comparing it to the linear spring, Kelvin, Hertz and
stereomechanical models using two Single Degree Of Freedom (SDOF) building systems
shown in Fig. 5 with the parameters shown in Table
1. The systems were subjected to a suite of thirty ground motion records
(Muthukumar, 2003) from thirteen different earthquakes with Peak Ground Accelerations
(PGA) levels varying from 0.1 to 1 which were carefully chosen to fall within
zone I.
Table 1: 
Properties of SDOF systems used for the impact model 


Fig. 8: 
Mean displacement amplification for T_{1}/T_{2}
= 0.3 

Fig. 9: 
Mean displacement amplification for T_{1}/T_{2}
= 0.5 
Three ground motion records were used at each PGA level. In the present study, the performance of the two models is investigated
in the process of studying buildings pounding. As suggested by Muthukumar and
DesRoches (2006), two cases are considered. First is a no pounding case where
the gap distance d is set to be large enough to avoid pounding. In the second
case, the gap distance d is chosen to be small enough to induce pounding. The
ratio of the maximum responses obtained in second case and the first case, called
the amplification response, is used in the analysis. We solve the equation of
motion (15) for the thirty ground motion records and the three different period
ratios. The obtained amplification factor for displacement and acceleration
of both SDOF 1 and SDOF 2 indicates that the nonlinear viscoelastic model provides
the smallest displacement and acceleration amplifications for the whole PGA
levels and the different period ratios considered herein as shown in Fig.
813.

Fig. 10: 
Mean displacement amplification for T_{1}/T_{2}
= 0.7 

Fig. 11: 
Mean acceleration amplification for T_{1}/T_{2}
= 0.3 

Fig. 12: 
Mean acceleration amplification for T_{1}/T_{2}
= 0.5 

Fig. 13: 
Mean acceleration amplification for T_{1}/T_{2}
= 0.7 
CONCLUSIONS
The comparison between two nonlinear models for earthquakeinduced structural pounding, i.e., the Hertzdamp model and the nonlinear viscoelastic model, has been conducted in this research. Both models are based on the Hertz contact law and incorporate additional damping in order to simulate the energy dissipation during impact. First, the numerical results obtained for each model have been compared with the results of two impact experiments. Then, the comparison has been carried out based on the results from the shaking table experiments. Finally, the investigation has been conduced for pounding between two single degree of freedom systems with different period ratios subjected to thirty ground motion records of different PGA levels. The effectiveness of two models of pounding has been assessed by comparing the displacement and acceleration response amplifications.
The results of the study demonstrate that the use of the nonlinear viscoelastic
model leads to smaller normalized errors in the impact force time histories
comparing to the Hertzdamp model. This is mainly due to the fact, that the application
of the Hertzdamp model results in longer time of contact as can be shown in
Fig. 2 and 4. On the other hand, the impact
force time histories obtained for the nonlinear viscoelastic model show the
change of the curvature after passing the peak force value and this disadvantage
is related to the disengagement of the damping term in the restitution period
of impact (Eq. 8).
Studies with the use of the results from the shaking table experiments showed that the Hertzdamp model gives smaller normalized error than the nonlinear viscoelastic model in simulation of impact velocity for pounding between two steel towers of different dynamics characteristics subjected to sinusoidal excitations.
The results of further analysis indicate that the nonlinear viscoelastic model provides smaller displacement and acceleration amplifications of the response of colliding single degree of freedom systems for three different period ratios and PGA levels of ground motions. This makes the nonlinear viscoelastic model more convenient in the simulation of earthquakeinduced structural pounding.
Both impact force models considered in this paper have been found to have some advantages and disadvantages when used for modeling of structural pounding. The results of the study indicate that the accuracy of each of the models depends on the type of analysis conducted.
REMARK 1: IMPLICIT RUNGEKUTTA (IRK) METHODS
The differential Eq. 1315 can be written
in a uniform version
where, u ∈ R^{n}, p : R^{n}→ R^{n} is a continuously differentiable function and f : R^{n} → R^{n} is a continuous function but not necessarily differentiable. Problem (17) is called a system of nonsmooth ordinary differential equations of the first order. In order to solve the system (17) efficiently we use the Implicit RungeKutta (IRK) method (Jay, 2000; Chen and Mahmoud, 2008; Mahmoud and Chen, 2008). One step of an sstage Implicit RungeKutta (IRK) method for solving (17) has the following version (Jay, 2000):
Given a step size h, a coefficient matrix A∈R^{sxs} and a weight vector b∈R^{s}. Let U_{0} = u_{0}. For k≥0:
Step 1: Solve the s x ndimensional system of nonlinear equations
to get a solution .
Step 2: Solve the ndimensional system of nonlinear equations
and let the solution be U_{k+1}.
A practical IRK method can be defined, by choosing appropriate matrix A and vector b, such as coefficients of Gauss, Radau IA and IIA, Lobatto IIIA, Burrage, etc. (Jay, 2000; Burrage, 1982). Moreover, we use the slanting Newton method (Chen et al., 2000) to solve the system of nonsmooth equations in each iteration of the IRK method. Numerical results show that the IRK method with the slanting Newton method is efficient. Numerical results reported in this study were obtained by using the twostages Burrage IRK method (Burrage, 1982) which has stable property (Ferracina and Spliker, 2005). No numerical problems were encountered.
REMARK 2: CONVERGENCE ORDER OF IRK METHODS
Chen and Mahmoud (2008) have studied convergence order of IRK methods and tested
various IRK methods with the slanting Newton method on numerous problems in
structural oscillation and pounding. In theory, it can be shown that the order
of convergence is at least one for the Lipschitz continuous ODEs under mild
conditions. Consider the solution u(t) to (17) in a fixed interval [0, T] with
the number n of step chosen such that t_{n} = nh = T. Let
and
For simplicity, we consider p (u) = u. Assume that for any U_{k} there
are x^{k}, U_{k+1} such that H(x^{k}) = 0 and .=
0 Let
Since f is Lipschitz continuous, the function G is Lipschitz continuous. By
the Rademacher theorem, G is differentiable almost everywhere. Hence, we can
define the Clarke generalized Jacobian as αG (x,U) (Clarke, 1983). By
the implicit function theorem for Lipschitz continuous function (Clarke, 1983),
for small h, there exist a neighborhood u_{k} of U_{k}
and a Lispchitz function ψ_{k}(.;h) :u_{k}→ R^{s×n} such that x^{k} ψ_{k}(U_{k}; h) and for every U∈u_{k}, G (ψ_{k}(U, h), U) = 0. Therefore,
we can write
It is easy to verify that φ_{k} (.;h): R^{n} → R^{n} is a locally Lipschitz continuous function.
Theorem 1: Chen and Mahmoud (2008) suppose that there are positive constants
K_{1} and K_{2} such that:
and
Then, there is a constant α > 0 such that
Figure 14 shows the convergence of IRK methods with different
coefficient versus the Explicit RungeKutta (ERK) methods of order four in computation
of the response for the collapse of the Tacoma Narrows suspension bridge.

Fig. 14: 
Convergence of ERK method of fourth order and Burrage
and Radau coefficients of the IRK method with various constant step sizes
hl 
The numerical experiments were performed using MATLAB 7.0 on a Dell PC with
2MB memory and 800 MHZ.