INTRODUCTION

The nonlinear model of impacts based on the Hertz contact law has been used in numerous studies of earthquake-induced 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 Kelvin-Voigt 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 second-order 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 Runge-Kutta (IRK) methods to solve the system of first order ordinary differential equations.

HERTZ CONTACT LAW FORCE-BASED 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).

(1)

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

(2)

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

(3)

where, h₁ and h₂ are the material parameters defined by the formula (Goldsmith, 1960):

(4)

Here, γ_i and E_i are Poisson's ratio and Young's modulus, respectively. When R₂→4, i.e., the second structure becomes a massive plane surface, the nonlinear spring stiffness k_h is defined as (Goldsmith, 1960):

(5)

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

(6)

where, e is the coefficient of restitution (Goldsmith, 1960), is the relative velocity during contact and v₁–v₂ is the relative approaching velocity prior to contact which can be expressed in terms of nonlinear spring and maximum deformation (Lankarani and Nikravesh, 1990):

(7)

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:

(8)

where, is the impact stiffness parameter and

(9)

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

(10)

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

(11)

where, F is the response time history obtained experimentally, is the response time history obtained numerically and is the Euclidean norm. and in case of the time histories given in a discrete form can be calculated as:

(12)

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.

Steel-to-steel 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 force-time 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₁ dropping onto a beam can be written by drawing the free body diagram (Jankowski, 2005b) as shown in Fig. 1:

(13)

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¹⁰ 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¹⁰ 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.

Concrete-to-concrete 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₁ and a prestressed fixed pile can be written by drawing the free body diagram as shown in Fig. 3 (Jankowski, 2005b).

(14)

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⁹ N/m^3/2 (Jankowski, 2005b). For the Hertzdmap model, k_h = 1.54x10⁹ 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 stand-of 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

(15)

where, 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).

(16)

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_gt) as an excitation acting on SDOF 1 and SDOF 2 with the following dynamic characteristics m₁ = 98.0 kg, f₁ = 5.04 Hz, ξ₁ = 7.2%, m₂ = 146.4 kg, f₂ = 2.76 Hz, ξ₂ = 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

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

As a second input shaking table, we consider the harmonic excitation 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. 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.

Table 1:	Properties of SDOF systems used for the impact model

Fig. 8:	Mean displacement amplification for T1/T2 = 0.3

Fig. 9:	Mean displacement amplification for T1/T2 = 0.5

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. 8-13.

Fig. 10:	Mean displacement amplification for T1/T2 = 0.7

Fig. 11:	Mean acceleration amplification for T1/T2 = 0.3

Fig. 12:	Mean acceleration amplification for T1/T2 = 0.5

Fig. 13:	Mean acceleration amplification for T1/T2 = 0.7

CONCLUSIONS

The comparison between two nonlinear models for earthquake-induced 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 earthquake-induced 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 RUNGE-KUTTA (IRK) METHODS

The differential Eq. 13-15 can be written in a uniform version

(17)

where, u ∈ Rⁿ, p : Rⁿ→ Rⁿ is a continuously differentiable function and f : Rⁿ → Rⁿ 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 Runge-Kutta (IRK) method (Jay, 2000; Chen and Mahmoud, 2008; Mahmoud and Chen, 2008). One step of an s-stage Implicit Runge-Kutta (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₀ = u₀. For k≥0:

Step 1: Solve the s x n-dimensional system of nonlinear equations

(18)

to get a solution

Step 2: Solve the n-dimensional system of nonlinear equations

(19)

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 two-stages 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

(20)

and

(21)

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

(22)

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 and for every U∈u_k, G (ψ_k(U, h), U) = 0. Therefore, we can write

(23)

It is easy to verify that φ_k (·;h): Rⁿ → Rⁿ is a locally Lipschitz continuous function.

Theorem 1: Chen and Mahmoud (2008) suppose that there are positive constants K₁ and K₂ such that:

(24)

and

(25)

Then, there is a constant α > 0 such that

(26)

Figure 14 shows the convergence of IRK methods with different coefficient versus the Explicit Runge-Kutta (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.