## Abstract

A new method is proposed that can deal with multi-impact problems and produce energetically consistent and unique post-impact velocities. A distributing law related to the energy dispersion is discovered by mapping the time scale into the impulsive scale for bodies composed of rate-independent materials. It indicates that the evolution of the kinetic energy during the impacts is closely associated with the *relative contact stiffness* and the *relative potential energy* stored at the contact points. This distributing law is combined with the Darboux–Keller method of taking the normal impulse as an independent ‘time-like’ variable, which obeys a guideline for the selection of an independent normal impulse. Local energy losses are modelled with energetic coefficients of restitution at each contact point. Theoretical developments are presented in the first part in this paper. The second part is dedicated to numerical simulations where numerous and accurate results prove the validity of the approach.

## 1. Introduction

Multiple impacts occur in many mechanical systems (granular matter and kinematic chains), and are a phenomenon of major importance. Systems with multiple impacts have a time-varying structure due to the unilateral features of the constraints (Delassus 1920; Bernoulli 1969–1993; Moreau 1994; Pfeiffer & Glocker 1996; Brogliato 1999; Stronge 2000*a*; Ceanga & Hurmuzlu 2001; Frémond 2002; Glocker 2004). Moreover, the configurations and the possible flexibilities at the contact points significantly influence the outcomes. Consequently, these phenomena require a special treatment for analytical and numerical solution methods. There are basically four ways of modelling a contact–impact phenomenon in multibody systems: (i) a finite-element method to directly discretize the contact bodies (e.g. Seifried *et al*. 2003; Schiehlen *et al*. 2006; Liu *et al*. 2007*b*), (ii) compliant models (spring–dashpot, with linear or nonlinear stiffness) and second-order dynamics (Khulief & Shabana 1987; Yigit *et al*. 1990), (iii) a purely algebraic impact law (or mapping) of the form , where *p* is a set of parameters (Han & Gilmore 1993; Moreau 1994; Pfeiffer & Glocker 1996; Frémond 2002; Leine & van de Wouw 2008), and (iv) the first-order Darboux–Keller impact dynamics (Darboux 1880; Keller 1986; Stronge 1994, 2000*a*; Batlle & Cardona 1998; Brogliato 1999; Ceanga & Hurmuzlu 2001; Liu *et al.* 2007*a*).

The approach (iv) is chosen in this paper for reasons that will be made clear below. Let us consider a mechanical Lagrangian system with an *n*-dimensional generalized coordinate ** q** and configuration space . This system is supposed to have a mass matrix , and is subjected to several unilateral constraints

*h*

_{i}(

**)≥0, 1≤**

*q**i*≤

*s*, where the functions

*h*

_{i}(.) are supposed differentiable and the set has an interior with positive volume. We define the

*s*hypersurfaces . A multiple impact is the one that occurs at the intersection of two (or more than two) surfaces

*Σ*

_{i}. Physically, this corresponds to having several points at which contact is established simultaneously.

There are many properties and characteristics that are associated with a multiple impact, independently of the chosen class of models (i)–(iv). Let us provide a rapid summary as follows:

the kinetic angle between the surfaces

*Σ*_{i}involved in the impact,the continuity of the solutions with respect to the initial data,

the kinetic energy behaviour at the impact,

the wave effects due to the coupling between various contacts,

the local energy loss during impacts,

the ability of the impact rule to span the whole admissible post-impact velocities domain,

the ability of the parameters defining the impact rule to be identified from experiments,

the (in)dependence of these parameters on the initial data,

the physical meaning of the parameters of the impact rule,

the ability of the impact rule to provide post-impact velocities in agreement with the experimental results,

the well-posedness of the non-smooth dynamics when the impact rule is incorporated in it,

the law should be applicable (or easily extendable) to general mechanical systems,

the determination of the impact termination, and

the impact law has to be numerically tractable.

Items (1), (3) and (6) enter into the geometrical considerations in Glocker (2004) to develop a framework for multiple impacts, extending Moreau's sweeping process rule (Brogliato 1999; Leine & van de Wouw 2008). Item (2) is important as it implies that sequences of single impacts may completely fail to deduce the multiple impact law (moreover, in general, impacts are not propagated in a sequential way through the contacts (Acary & Brogliato 2003; Wei & Liu 2006)). It is also related to (11) and (1). Indeed, under some conditions on the kinetic angles, continuity of the solutions (2) may hold. These cases may be considered as isolated cases, however. Item (4) implies that the coupling between various contacts will generate possible ‘distance’ effects, i.e. a collision may be propagated through the system and induce a detachment at another contact point (this is also called the *dispersion*, or the *scattering* effects). This implies that other relationships more than the usual restitution coefficients at each contact have to be discovered. For instance, Ceanga & Hurmuzlu (2001) were the first ones to quantitatively study the coupling effects by introducing an impulse correlation ratio (ICR) together with energetic coefficients of restitutions. Their algorithm is tested on chains of balls and on the rocking block (Hurmuzlu & Ceanga 2001). Another approach using the ICR and a global energetic coefficient is proposed by Acary & Brogliato (2003). Both Ceanga & Hurmuzlu (2001), Hurmuzlu & Ceanga (2001) and Acary & Brogliato (2003) belong to class (iv) of the models. Other examples with different ‘coupling’ coefficients may be found in Moreau (1994) and Frémond (2002); however, they may fail to satisfy item (9), rendering (7) tricky. Item (5) is related to the meaning of the restitution of coefficient. Newton's coefficient reflects the local energy loss on the velocity level, while Poisson's coefficient puts the energy loss on the impulse level. Obviously, these variables must change due to item (4) and may violate the basic law of energy conservation (Kane & Levinson 1985; Leine & van de Wouw 2008). Items (7), (8), (9) and (10) are related to the choice of the parameters that enter the impact law. They are fundamental from the point of view of the practical usefulness of the law. For instance, one may want that some parameters are to be identified through simple experiments before being injected in a rule for a more complex collision. It is noteworthy that the main drawback of models (ii) is the necessity to identify stiffness and damping parameters, which is usually not easy when there are many contacts with various materials. Item (12) is closely linked to (4), since the wave effects may differ a lot depending on the bodies' shapes, materials, etc. Item (4) also indicates that a guideline should be established for item (13), since the separations of contacts do not occur simultaneously. The last item is not the least one. It has led several researchers to look for multiple impact laws that lend themselves to a linear complementarity problem formulation (Moreau 1994; Pfeiffer & Glocker 1996; Frémond 2002; Glocker 2004), and are extensions of the Newton or Poisson restitution models. Pioneering works related to the multi-impact problems can also be found in several papers (see, e.g. Delassus 1920, Bernoulli 1969–1993, Han & Gilmore 1993 and Brogliato 1999 for references).

The contribution of this paper is to extend the Darboux–Keller approach (iv) to multiple impacts. This allows us to discover a distributing rule that links the ratios of infinitesimal impulses to the relative contact stiffnesses and the relative potential energies. This, together with the energetic coefficients of restitution and a suitable compliant contact model, results in a set of first-order nonlinear differential equations using normal impulse at the primary contact as an independent variable. These equations are non-stiff and depend on mechanical parameters that satisfy items (7), (8) and (9). In the second part of the paper (see also the report Liu *et al*. 2008), it is shown, through comparisons between numerical and experimental results, that the proposed model encapsulates very well the main features of multiple impacts in chains of balls.

The paper is organized as follows. In §2, a general method to solve multi-impact problems by using impulsive dynamics is developed, in which the ratios for the distribution of the increments of normal impulses at all contact points are established based on a mono-stiffness contact model. Section 3 presents a similar relationship for the distributing law, which is deduced from a bi-stiffness model that satisfies the energetic constraint. A compact formulation for the analysis procedure of multi-impact problems is presented in §4. Conclusions are given in §5.

## 2. The Darboux–Keller multi-impact dynamics

### (a) The impulsive dynamics

Let us consider a multibody system with *s* frictionless contacts. The maximum number of degrees of freedom *n* is obtained when none of the contacts is closed. In this state, the system may be described by a set of generalized coordinates ** q**∈

^{n}, and the equations of motion take the form(2.1)where and

**(.,.,.) contain the inertial and applied forces, respectively.**

*h***=[**

*Λ**λ*

_{1},

*λ*

_{2}, …,

*λ*

_{s}]

^{T},

*λ*

_{i}is the scalar value of the normal contact forces along the common normal to the surfaces of the contact bodies for contact point

*i*. The connection between the normal contact forces and the generalized forces is defined as , which is related to the Jacobian matrices of the contact points (Pfeiffer & Glocker 1996).

The kinematic state of a contact is determined by the distance *δ*_{i}(** q**,

*t*) between the contact bodies. Clearly,

*δ*

_{i}(

**,**

*q**t*) represents the normal deformation when contact is established. The relative velocity of the contact points is expressed as(2.2)where is a nonlinear term related to time, and is often zero if all constraints, excluding contact, are ideal and time invariant. The directions of the relative normal velocities depend on the contact kinematics, and we always define that for the colliding bodies approaching (compression phase), while for separation (expansion phase). In matrix notation, equation (2.2) becomes (we drop the arguments)(2.3)Let [

*t*

_{0},

*t*

_{f}] denote the time interval of the impact, which can be further divided into much smaller intervals [

*t*

_{i},

*t*

_{i+1}]. According to the Darboux–Keller model (Brogliato 1999, §4.2.5), an integration over [

*t*

_{i},

*t*

_{i+1}] has to be carried out in order to achieve a representation of the equations of motion at the impulse level. Thus, a set of differential equations with respect to the normal impulses can be obtained,(2.4)The terms

**and**

*M***remain unchanged during the integration because of the assumption of constant configuration**

*W***on [**

*q**t*

_{0},

*t*

_{f}]. The vector

**consists of finite, non-impulsive terms and therefore vanishes by the integration (Brogliato 1999, ch. 1). The quantities and d**

*h***are the changes of generalized velocities and normal impulses during [**

*P**t*

_{i},

*t*

_{i+1}], respectively.

### (b) The distributing rule for the normal impulses in a mono-stiffness model

Based on experiments or some approximate theory, different types of local compliance models could be found in the existing literature. If no dissipation occurs at contact points, usually the force–indentation mapping at the contact point *i* could take the following form:(2.5)where *K*_{i} is the contact stiffness and the exponent *η*_{i} determines the kind of contacts between bodies (*η*_{i}=(3/2) is for a Hertz contact, *η*_{i}=1 is linear elasticity). The variable *δ*_{i} is the normal deformation, which is assumed to be only a function related to the generalized coordinates ** q**. Therefore, the term in equals zero in (2.2).

This force–indentation mapping is denoted as a *mono-stiffness model*, since it is identical for the compression and expansion phases. Let *P*_{i}(*t*) denote the total normal impulse accumulated during the time interval [0,*t*], . So, and(2.6)In terms of the compliant model expressed by (2.5), we have(2.7)Note that *δ*_{i} can always be expressed as(2.8)Substituting (2.8) and (2.7) into (2.6) leads to(2.9)The initial value of the normal impulse can be set to *P*_{i}(0)=0, and the static contact force before impact is *λ*_{i}(0)=0 for the case without initial precompression energy. The integration of equation (2.9) leads to(2.10)Noting that *λ*_{i}=d*P*_{i}/d*t* and considering only the variation in space, the ratio of the changes of normal impulses at the contact points *i* and *j* can therefore be expressed as(2.11)We should note that d*P*_{i} and d*P*_{j} are the variations of the normal impulses of the two contact points during the same time interval. Meanwhile, *P*_{i}(*t*) and *P*_{j}(*t*) are the accumulated normal impulses during the same time interval [0,*t*]. The work functions are defined as(2.12)It is easy to find that *E*_{i} and *E*_{j} are just the work done by the normal contact forces at contact points *i* and *j* from the beginning of impacts to the time when the respective impulses are *P*_{i} and *P*_{j}.

These terms can also be thought of as the potential energy stored in the springs at the contact points *i* and *j*. During the compression phase, the spring will transfer the kinetic energy of the contact point into potential energy, such that *E*_{(.)} increases with *P*_{(.)}. Once the potential energy at the contact point (.) is saturated due to , the expansion phase of the spring will begin to release the potential energy that has been stored, such that *E*_{(.)} will decrease with *P*_{(.)} and the kinetic energy is transferred. Clearly, the exchange process between kinetic and potential energies will finish at the instant when *E*_{(.)}=0.

Expression (2.11) can be further simplified if *η*_{i}=*η*_{j}=*η*. In this case, we can introduce the ratios of contact stiffnesses *γ*_{ji}=*K*_{j}/*K*_{i}, and define(2.13)to represent a ratio of the energies stored at any time *t* at the contact points *i* and *j*. This ratio varies during the impacts to make the contact bodies change their motions. At the same time, the relationship between d*P*_{i} and the changes of normal impulses at other contact points can be expressed using (2.11) and (2.13) as(2.14)Obviously, these expressions reflect the wave behaviours generated in multiple impacts and depend only on the properties of the contact constraints: *the relative stiffness and the relative potential energies accumulated in the contact points*.

### (c) The selection for the independent variable

Combining the distributing relationships (2.14) with the Darboux–Keller model (2.4), the impulsive equations (2.4) can be expressed as a set of first-order differential equations with respect to a single integral variable related to the normal impulse *P*_{i}, in which the new integration dumb variable is *τ*=*P*_{i}(*t*).

A problem arising when solving the first-order differential equations is which normal impulse among all the contact points can be selected as a ‘time-like’ independent variable. Since all normal impulses are initially monotonously increasing for the points keeping contact, in principle, any one among all the normal impulses can be taken as the independent variable. However, inappropriate selection for the independent variable will result in certain numerical difficulties. For instance, if the normal impulse at a contact point with little potential energy is selected as the independent variable, it is clear from (2.13) and (2.14) that a little change of the independent variable will make the normal impulses at other contact points vary abruptly. In order to avoid the numerical difficulties and respect the physical meaning that the multiple impacts should be dominated by the contact point that takes the maximum value of potential energy, we present a guideline for the selection of the independent variable. If the energy at the contact point *i* satisfiesthen the normal impulse corresponding to this contact point can be selected as the independent variable. We denote it as the *primary colliding point*.1 Once the primary colliding point is determined, the impulse related to this point can be considered to increase monotonously like the time variable. The normal impulses at other contact points will increase according to the relationship (2.14).

At the beginning of the impact, a problem will arise for the selection of the independent variable since no energy is stored at any contact point (recall that we assume there is no precompression between the bodies, i.e. *δ*_{i}(*t*_{0})=0 for all *i*≥2 and *λ*_{i}(0)=0 in (2.10)). This corresponds to a singular point in the simulation. Let us denote Δ*P*_{i}(1) and Δ*P*_{j}(1) the possible increments of the normal impulses. and are the initial relative velocities at the contact points *i* and *j*. For simplicity, the exponents at various contact points are assumed to be equal (*η*_{i}=*η*, *i*=1,2, …, *s*). In terms of (2.11) and (2.2), we have(2.15)This form can be further simplified to(2.16)Obviously, the impact behaviours will be firstly dominated by the contact point at which the relative velocity takes the maximum value. Therefore, the normal impulse at this point can be selected as the initial independent variable. The normal impulses at other contact points will vary with the independent variable in terms of the expression (2.16).

### (d) The energetic constraint for the local energy loss

The strong interactions between contact points usually dissipate a part of the energy that cannot be recovered by the expansion phase. For the colliding bodies made of materials without viscosity (i.e. rate-independent materials), the dissipation is mainly due to the plastic deformation at the local contact region (Johnson 1992). So, the energetic coefficient can always be expressed as a function with respect to the work done by the contact force during the compression phase (the potential energy stored in the elastic deflection; Stronge 2000*a*). Then, we can use *e*_{j} as an energetic constraint to define the local energy loss at contact *j*. According to the definition given by Stronge (2000*a*) and his predecessors such as Routh and Boulanger (see Brogliato 1999, p. 147), the energetic constraint *e*_{j} at contact *j* is given by(2.17)where *W*_{c,j}≥0 and *W*_{r,j}≤0 are the work done by the normal contact force at point *j* during the compression phase [0,*t*_{c}] and the expansion phase [*t*_{c},*t*_{f}], respectively. Obviously, *W*_{c,j} also corresponds to the potential energy accumulated during the compression phase, and can be obtained by summing the scalar product of d*P*_{j} and from the beginning of the impact to the instant of *P*_{j}(*t*_{c}) that corresponds to . Thus,(2.18)and(2.19)If the contact point *j* is the primary colliding point, d*P*_{j} will define the size of the numerical step. Otherwise, d*P*_{j} is calculated by the distributing rule in (2.14). The time *t*_{c} can be thought of as an instant when the potential energy is saturated. After that, the potential energy will be transferred into kinetic energy through the expansion phase. Since *e*_{j} can always be identified as off-line (as a function of *W*_{c,j}), the total recovered energy *W*_{r,j} can be determined by *e*_{j} and *W*_{c,j}. For instance, *W*_{r,j} equals *W*_{c,j} for a fully elastic impact and *W*_{r,j}=0 for a plastic impact. At any instant, during the compression phase *P*_{j}(*t*), the residual energy equals(2.20)When *P*_{j}(*t*) is located in the expansion phase, *E*_{j} is(2.21)Since is continuous and equal to , the expressions (2.20) and (2.21) can always take the same form as in (2.12), and represent the potential energy stored at the contact point *j*. If the normal impulse reaches the instant of *P*_{j}(*t*_{f}) satisfying(2.22)the process of energy transfer at the contact point *j* will finish, as the residual potential energy *E*_{j}(*P*_{j}(*t*_{f})) will be dissipated based on the energetic constraint expressed by (2.17). At this time, the outcome of the post-impact velocities at this contact point can be obtained if it does not again participate in impacts, and all the ones related to the termination of multiple impacts can be determined when all the accumulated energy is completely released.

Some complex situations may appear in the process. For example, the contact point may experience multiple compression–expansion phases due to the interactions between contact points,2 i.e. the relative velocity may change from , related to an expansion phase, to , which corresponds to a new compression phase. We denote this as the *multiple compression phenomenon*. In this case, before all the potential energy *W*_{c,j} is transferred into the kinetic energy (the kinetic energy that can be transformed is equal to *W*_{r,j}), the contact point may experience a new compression phase to absorb the external kinetic energy. Let us denote the beginning instant of the new expansion phase as , in which(2.23)During the new compression phase, *E*_{j} will increase again and obtain a new value at the instant , in which . We can put the energetic constraint *e*_{j} on the potential energy to define the kinetic energy that will be recovered during a new expansion phase. In other words, if(2.24)the impact at the contact point *j* finishes at the instant .

When the force–indentation relationship in (2.5) is used, it is obvious that *e*_{j} cannot constrain the energy loss appearing in the compression phases before the last one, which is detected when equality (2.24) is satisfied. This flaw can be overcome by using the bi-stiffness compliant model in §3.

Another interesting phenomenon can also be found in the process of multiple impacts. At the same contact point, bodies will participate again in impacts after separation. We denote this case as a *repeating impact*. Obviously, a new process of energy accumulation and transfer will appear at that point, and the energetic constraints can be applied on each isolated impact based on the potential energy obtained from the compression phases. It is possible to introduce a variation of the restitution coefficient as a function of the number of impacts, based, for example, on experimental results (Schiehlen *et al*. 2006). If the plastic deformation remains small enough, we may, however, assume in a first study that the coefficients *e*_{j} are constant, even when repeated impacts occur.

## 3. The distributing rule for the compliant bi-stiffness contact model

The distributing rule in (2.14) is based on the assumption that the compression and expansion phases satisfy the same relationship between the contact forces and the indentation. The energetic coefficient is applied as a constraint to take into account the local dissipation of energy at the last compression phase. We should, however, use two different modes of the contact forces to represent the compression and expansion phases, respectively, such that the dissipated energy at each loading–unloading cycle can be involved in the compliant contact model (Schiehlen *et al*. 2006). In this section, we will adopt a compliant contact model that satisfies the definition in (2.17) to discuss the distributing rule.

### (a) The bi-stiffness compliant contact model

A compliant contact model has been proposed in Lankarani & Nikravesh (1994) to represent the relationship between the contact force and the indentation for the compression and expansion phases, respectively. This model lacks fundamental mechanical meaning as no material seems to satisfy such laws; however, it has been proved to be quite useful in various contexts of impact dynamics (Sadd *et al*. 1993; Stronge 2000*a*,*b*, 2003). This bilinear model is the simplest means of representing the plastic deformation (Schiehlen *et al*. 2006). It has certain drawbacks such as not limiting maximal force. Nevertheless, it gives a correct form of dissipation for rate-independent materials. It is not supposed to accurately represent any real material, and may be seen as an approximate representation of the plastic indentation effect in elastic solids (Johnson 1992, §6.4). Figure 1 shows the bi-stiffness compliant contact model by setting different force–indentation relationships for the expansion and compression phases. The relationship for the compression phase at the contact point *j* is expressed as(3.1)and the one for the expansion phase is(3.2)where *δ*_{r,j} is the plastic deformation, and *λ*_{M,j} and *δ*_{M,j} correspond to the maxima of the normal contact force and normal deformation at the end of the compression phase, which corresponds to the values when . Clearly, the dissipated energy is just the area enclosed by the compression and expansion curves. For simplicity, we omit the subscript *j* in the following expressions (and also in figure 1).

The scalar *δ*_{r} represents the permanent plastic deformation generated when the compression phase finishes (the elastic part of the deformation generated in the compression phase will be recovered in the expansion phase). Its value should depend on the dissipated energy at the contact point, such that the energetic constraint should be applied. At any contact point *j*, work done during the compression and restitution phases is obtained from the integration of the expressions (3.1) and (3.2)(3.3)Based on the assumption in (2.17), *W*_{c} and *W*_{r} can be connected by a function *e*=*f*(*W*_{c}). In the case of a constant *e*, we have from (3.3) and (2.17)(3.4)so that *e*=1 implies *δ*_{r}=0 (no plastic deformation occurs). The indentation *δ*_{M} corresponds to the first instant *t*_{c} such that , so it is not a parameter of the impact dynamics. From (3.4), neither is *δ*_{r}. When the energetic coefficient takes a constant value, the energetic definition presents a linear relationship between the local plastic deformation *δ*_{r} and the maximum indentation *δ*_{M}, such that the dissipation of energy is reflected in an average level.

### (b) The potential energy at a contact point

The potential energy at a contact point will be accumulated during the compression phase and released in the expansion phase. Let us pick a point p in the curve of the compression phase. The potential energy *E*_{p} at p is equal to the work done by the contact force along the path ,(3.5)Let us assign a time factor on the compression process. *δ*^{p} denotes the deformation when the contact point moves along the path to the point p within a time interval [0,*t*]. This means that, at any instant *τ* during the interval [0,*t*], *λ*_{c}(*τ*) can always be expressed as *λ*_{c}(*τ*)=(d*P*(*τ*)/d*τ*), in which *P*(*τ*) is the normal impulse accumulated in the time interval [0,*τ*]⊂[0,*t*] corresponding to a deformation *δ*(*τ*). Thus, we have(3.6)Since *P*(*τ*) and *δ*(*τ*) can be connected by a one-to-one mapping during the compression phases, we can use the variable *P*(*τ*) to replace *δ*(*τ*) as the integral variable. Thus,(3.7)where *P*(*t*) is the normal impulse that is needed to make the indentation change from zero to *δ*^{p} by obeying the relationship (3.1).

The accumulation of energy at the contact point *j* will be ended when the compression process finishes. After that, the potential energy will be released through an expansion phase. Let us select a point R in the curve of the expansion phase , as shown in figure 2. The work done by the contact force along the compression path can be expressed as(3.8)This term is related to the potential energy accumulated in the compression phase, and corresponds to the area enclosed by the curve . The total energy that can be recovered through the expansion phase is *e*^{2}*W*_{c}, based on the energetic constraint and is associated with the area enclosed by the curve .

When the contact point *j* moves from M to R along the expansion curve, the recovered energy is , where *δ*_{R} is the deformation related to the contact force *λ*_{R},(3.9)The second term in expression (3.9) equals the area enclosed by the curve . Let us determine the value of the potential energy *E*_{R} when the contact point is located at R, which should make the expansion force generate the work with a value equal to the area enclosed by the curve , such that the contact point can move from R to along the expansion curve. The following result is based on the properties of the contact force and the energetic constraint and deals with the residual energy.

*When the contact point is expanded from* M *to reach the position* R*, the residual potential energy E*_{R} *at this point is equal to the work done by the contact force along the compression path moving from* O *to the position* R′*,**in which the contact force takes the same value λ*_{R} *as the one in position* R.

The proof is omitted for the sake of briefness of the paper. It can be found in Liu *et al*. (2008). ▪

Theorem 3.1 shows that a compression–expansion cycle (curves 1 and 3) is equivalent, from the energetic point of view, to a cycle (curves 2 and 4), where the compression would end at R′ (respectively at R). The work done by the expansion force from R to along curve 4 can be expressed as(3.10)Substituting (3.10) into (3.9) leads to(3.11)Similar to the compression phase, the residual potential energy for a dynamical process reaching the position R within a time interval [0,*t*] can be expressed as(3.12)where *P*_{c} is the normal impulse when the compression phase finishes (so *δ*=*δ*_{M}), and *P*(*τ*) is an integral variable that is related to the normal impulse of the contact force experiencing a time interval [0,*τ*] by obeying the energetic relationship defined by the compliant contact model.

### (c) The energetic constraint for a complex impact

In some cases, the contact point may take some initial energy due to the precompression or the experience of a complex process with multiple compression phases owing to the coupling between the contact points. In a physical sense, it is obvious that the dissipated energy should depend on the energy stored at the contact points. Here, we will extend the energetic coefficient *e* into the situations where the contact point has some initial energy or is experiencing multiple compression phases (the multiple compression phenomenon).

Let us assume that the contact point *j* has an initial pressure *λ*_{0} that makes the contact point with an initial energy *E*_{0} and an initial deformation *δ*_{0}. Furthermore, we suppose that *λ*_{0} and *δ*_{0} satisfy the compression relationship expressed in (3.1). Then, integration yields(3.13)Under the initial pressure, the energetic constraint is defined as follows:(3.14)where is the work done by the expansion force; is the work done by the compression force; and *δ*_{r} is in (3.4) and figure 1.

Involving the initial potential into the energetic constraint permits the use of the coefficient *e* as the index to describe the local dissipated energy for the contact point with multi-compression phases.

Let us present a clear scenario by conducting an example of the contact point *j* with two compression phases. Figure 3 shows that the contact point *j* will first experience a compression phase along the curve , then begin an expansion process from M_{1} to R. Before the total potential energy is released, a new compression phase starts at R, and the contact point will begin a second compression phase along the curve . This second compression phase stops at M_{2} and then the accumulated potential will be completely released through an expansion phase along the curve .

Let us assume that the relationship between the contact force and the indentation is not changed for the second compression–expansion cycle. Therefore, the residual potential *E*_{R} at R can be thought of as an initial energy for a new compression–expansion cycle. So, from (3.14) we get(3.15)The energy *E*_{R} can be obtained by using (3.11) or (3.12),(3.16)Let us set a time interval [0,*t*] on the micro-movement of the contact point and replace *t* by . The potential energy at the point Q of the force–indentation mapping in different phases can be expressed as(3.17)where *E*_{(.)} and *P*_{(.)} are the residual potential energy and the normal impulse at the location (.)=M_{1}, R and M_{2}, respectively.

### (d) The distributing rule for the bi-stiffness compliant contact model

Since the compression and expansion forces adopt different relationships, we should separately analyse the evolution of energy during the compression and expansion phases. Let us suppose that the initial pressure at the contact point *j* is *λ*_{0,j}, such that, at the beginning of the impact process, the point will take the initial energy *E*_{0,j}, given by (3.16). Based on (2.13) and setting *λ*_{j}(0)=*λ*_{0,j}, we can obtain(3.18)where *E*(*P*_{j}(*t*)) takes the form of the first term in (3.17). In terms of the compliant model expressed in (3.2) for the expansion force, we can deduce(3.19)Since and(3.20)the expression (3.19) can be further simplified to(3.21)At the end of the compression phase, the maximum compression force *λ*_{m,j} can always be expressed as(3.22)So,(3.23)The initial value of the normal impulse at the beginning of the expansion phase is set as *P*_{c,j}, which is related to . The contact force at this instant, i.e. *λ*(*P*_{c,j}), can be obtained using (3.18)(3.24)The integration of (3.23) and using (3.24) leads to(3.25)Thus, the contact force at the impulse instant *P*_{j}(*t*) is(3.26)Since the contact forces can always be expressed as the differential of the normal impulse with respect to time, from (3.26) and assuming the exponent *η*_{j} in the force–indentation relationship takes equal values *η* at all contacts, the ratios between the increments of normal impulses among various contact points *j* can be expressed as in (4.2), i.e.(3.27)The distribution of the increments of normal impulses depends on *the relative stiffness and the relative potential energy* among various contact points.

The distributing law in (3.27) takes the same form as the one expressed in (2.14), even though we adopt two different kinds of constitutive relationships for the compliances at the contacts. One may thus expect that an impulsive process with any kind of compliance can be dominated by the underlying law with the form expressed in (3.27). In other words, the evolution of motion in an impulsive process just depends on the relativity of the contact stiffness and of the potential energy resided in the system. This may extremely facilitate the understanding of the energy transmitted through a network of contacts.

## 4. The compact formulation for the analysis procedure of multi-impact problems

Based on the distributing law proposed in this paper for mono- and bi-stiffness models, multi-impact problems can be solved using the Darboux–Keller impulsive equations plus the energetic constraint to take into account the local energy loss. Let us provide a compact formulation for the analysis procedure of multi-impact problems (the exponents *η*_{i} at various contacts take the same value *η*).

Contact parameters:

*γ*_{ij};*e*_{j}; 1≤*i*≤*s*; 1≤*j*≤*s*(*η*=1(linear stiffness) or=3/2(Hertz contact), or other suitable values);*E*_{0,j}; and*λ*_{0,j}for the initial precompression cases.Dynamical equation:(4.1)with(4.2)and(4.3)The calculation for the potential energy

*E*_{j}(*P*_{j}) and the termination of each contact depend on the compliance model added at the contact points.Mono-stiffness model:(4.4)where the time

*t*_{**}at the contact*j*is calculated from for the last loading–unloading cycle, while*t*_{*}is the instant when the compression phase in the last cycle begins. For the contact with a single loading–unloading cycle,(4.5)The impact at contact*j*will be terminated at the instant*t*_{f}when . Part of the potential energy in the expansion phase of the last cycle is discarded in order to take into account the local energy loss.Bi-stiffness model:(4.6)Tra is a parameter to transfer the work done by the normal impulse into the potential energy, in which Tra=1 if and if . When and , the contact

*j*will be open at time*P*_{j}(*t*_{f}).

The solution of the set of ordinary differential equations (4.1)–(4.6) is the function on the impulsive intervals, where *P*_{i} is the dominant impulse. It is denoted as the derivative of the position for obvious reasons; however, it is to be considered as the solution of this particular set of differential equations since the energy is mapped into the velocity–impulse level. The assumption that the position is constant is imposed for equation (4.1), and the infinitesimal time interval for the impact duration is adopted for the calculation of the post-impact velocities.

## 5. Conclusions

This paper presents a method to solve impact problems with multiple frictionless points of contact between bodies composed of rate-independent materials. The theoretical analysis indicates that, except for the local energy loss, the dispersion of energy during impacts must be considered by adding compliances into the rigid body model. The local energy loss is confined by using an energetic constraint that can be defined by energetic coefficients, while the dispersion of energy during impacts (or the wave effects) is well represented by mapping the impact process into a velocity–impulse level, hence extending the Darboux–Keller approach to multiple shocks. A distributing law that seems to be independent of the compliant contact model is developed by using mono-stiffness and bi-stiffness models. It shows that the distribution of the increments of normal impulses in space is associated with the relative potential energy and the relative stiffness between various contact points. Combining this distributing law with the impulsive dynamics, one can deduce a set of first-order nonlinear differential equations with respect to an independent time-like normal impulse that corresponds to a so-called primary colliding point, in which the potential energy takes the maximum value. This method clearly separates dissipation and wave effects during multiple impacts and gives a legible picture for the dynamics of multiple impacts. The validity and the advantages of this method are illustrated in part II (Liu *et al*. in press) and Liu *et al*. (2008), where numerous numerical results are presented and carefully compared with the experimental results.

## Acknowledgments

The support of the National Science Foundation of China (10772002 and 10642001) is gratefully acknowledged.

## Footnotes

↵The primary colliding point may change from one contact point to another during the impact due to the exchange and transformation of the energy. In some cases, there may be several possible candidates for the primary colliding point (when equality holds). There is no strict limitation for the selection of the primary impulse when

*E*_{i}and*E*_{j}are near the same value for some*i*and*j*.↵Interestingly enough, such phenomena have also been noted for single impacts with Coulomb friction (Batlle & Cardona 1998).

- Received February 25, 2008.
- Accepted July 16, 2008.

- © 2008 The Royal Society