## Abstract

Part I of this paper develops a framework that is an extension of the Darboux–Keller shock dynamics towards frictionless multiple impacts between bodies composed of rate-independent materials. A numerical algorithm is proposed in this paper, in which the impulsive differential equation is discretized with respect to the primary impulse corresponding to the contact with the highest potential energy, and the integration step size is estimated at the momentum level. This algorithm respects the energetic constraints and avoids the stiff ordinary differential equation problem arising by directly using the compliant model. The well-known example of Newton's cradle, as well as Bernoulli's system, is used to illustrate the developments.

## 1. Introduction

In this paper, an algorithm associated with the theory developed in part I (Liu *et al*. 2008) for frictionless collisions between bodies composed of rate-independent materials is proposed. Several simple systems such as Newton's cradle and Bernoulli's problem are used to illustrate and validate the scheme. Let us first introduce a brief summary of the developments of part I.1

We consider a Lagrangian system with *n* degrees of freedom and *s* unilateral constraints. The kinematic state of a contact is determined by the distance *δ*_{i}(** q**) between the contact bodies, and the relative velocity of the contact point

*i*is expressed as , where

*w*_{i}(

**) is related to the Jacobian matrix of the contact point. The directions of the relative normal velocities depend on the contact kinematics as illustrated in Pfeiffer (1999), and we always define that for the colliding bodies approaching (compression phase), while for separation (expansion phase). We add compliances with the mono-stiffness model (2.5)**

*q*_{I}or the bi-stiffness model, (3.1)

_{I}and (3.2)

_{I}, into rigid body models, and use energetic coefficients

*e*

_{i}at the contact points to constrain local dissipation of energy at various contacts. By transferring the time scale into the impulse scale, and assuming that the exponent

*η*in the force–indentation relationship takes equal values at all contacts, and the same for compression and expansion, the Darboux–Keller shock dynamics that are used to compute the post-impact velocities can be expressed as follows.

Contact parameters:

*γ*_{ij};*e*_{j}; 1≤*i*≤*s*; 1≤*j*≤*s*; and*E*_{0,j}.Dynamical equation:(1.1)with(1.2)and(1.3)where

is the transpose of a Jacobian matrix related to the contact point velocities.*W*Mono-stiffness model:(1.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,(1.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:(1.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}).

In (1.1)–(1.3), *γ*_{ji} and *E*_{ji} are the ratios of the contact stiffnesses and of the potential energies *E*_{i} and *E*_{j} between points *i* and *j*, respectively; *e*_{j}∈[0,1] is the energetic coefficient of restitution at the contact point *j*; *η* characterizes the elasticity law (*η*=1 (linear stiffness) or =3/2 (Hertz contact), or other suitable values); is the work of the contact forces during the compression phases of the last cycle; and *E*_{0,j} is the initial precompression energy at the contact point *j*. Let us recall that the position ** q** is assumed to be constant, and that the inertia matrix may depend on

**in (1.1). The**

*q**distributing rule*in (1.2) implies that the increase of normal impulse at the contact point

*i*(d

*P*

_{i}) will drive the normal impulses at other contacts (d

*P*

_{j}) continuously changing in terms of the relative contact stiffnesses and the accumulation of energy between contacts

*i*and

*j*. Once a normal impulse at the contact with the highest potential energy is selected as an independent variable (the

*primary normal impulse*), the distribution of all normal impulses at a certain time in space can be determined according to (1.2). This relationship presents a clear picture for the evolution of multiple impacts. When the normal impulse is used as a ‘time-like’ independent variable for a set of first-order velocity/impulse dynamical equations using Darboux–Keller's approach (Darboux 1880; Keller 1986), some problems will arise in the implementation of the algorithm, such as the numerical singularity due to the variable structure of the equations generated by the repeating impacts, or the multi-compression processes, the switch of the primary colliding point due to the change of the potential energy at contacts, as well as the selection of the step size for the primary impulse that is the independent variable. The robustness of the algorithm should also be investigated since different energetic constraints are respectively applied into the mono-stiffness and bi-stiffness compliant contact models. The objective of this paper is to present a numerical procedure and numerical results to illustrate the developments. More results can be found in Liu

*et al*. (2008).

The organization of this paper is as follows. Section 2 provides a procedure for the time discretization of (1.1)–(1.6). Several chains of balls are presented in §3 to test the method by comparing its outcomes with the solutions obtained from compliant methods and the experimental results found in Ceanga & Hurmuzlu (2001). The investigation of the method applied into Bernoulli's example of impacts in the plane is carried out in §4. This paper ends with conclusions in §5.

## 2. The numerical algorithm

In this section, the discretization of the Darboux–Keller shock dynamics (1.1)–(1.6) is presented. Let *N* denote the number of steps and adopt a step size Δ*P* for the primary impulse. Obviously, the accuracy and efficiency of the numerical results depend on the selection of *N* and Δ*P*. In order to choose reasonable values of Δ*P* and *N*, we can define the following norm to estimate the possible value of the normal impulse at the momentum level:(2.1)where *m*_{ij} is the *ij*th element of the inertia matrix ** M** in (1.1), and is the initial generalized velocity of the system with

*s*frictionless contacts and

*n*degrees of freedom. The step size can be set as(2.2)At the beginning of the impact, the primary colliding point

*i*(see §2

*c*of part I) can be selected as the point that takes the maximum value of the relative velocity(2.3)and is the vector of initial relative velocities that can be obtained using equation (2.3)

_{I}. Then, the increment of normal impulse at the contact point

*i*can be set as , and the increments at the other contact points can be determined according to (2.16)

_{I}, in which(2.4)In most cases, we can take a constant value Δ

*P*as the step size. However, we may change the step size in order to find some critical points. Let us denote Δ

*P*

^{(l)}as the step size at step

*l*and use Euler's explicit difference scheme to discretize the impulsive differential equations (2.4)

_{I}. The quantities with the superscript (

*l*) represent their initial values at step

*l*, and the ones with the superscript (

*l*+1) represent the terminal values. The difference formulations for (2.4)

_{I}can be expressed as(2.5)where

*Γ*^{(l)}is a matrix related to the ratios of contact stiffness and those of the potential energy at the instant

*P*

^{(l)}, i.e.(2.6)Let us assume that

*i*is the primary colliding point that takes the normal impulse with the maximum potential energy, then

*Γ*^{(l)}in (2.5) can be expressed as(2.7)where

*γ*

_{(.),i}is the relative stiffness at the contact point (.) compared with the primary colliding point

*i*, and is the relative potential energy of the contact point (.) at the impulse step

*l*, and is defined by (2.3). Once is obtained, we can use the expression (2.3)

_{I}to obtain the relative velocity is set to zero and

*W*^{T}is a constant matrix) at the instant

*P*

^{(l+1)},(2.8)If , the contact point

*j*is located in the compression phase. According to the expression (2.20)

_{I}for the mono-stiffness model and the expression (3.7)

_{I}for the bi-stiffness model, the potential energy of the contact

*j*at the instant

*P*

^{(l+1)}can be calculated by(2.9)If , the contact point

*j*is located in the expansion phase. Therefore, the work done by the normal impulse during the interval Δ

*P*

^{(l)}is negative, and part of the potential energy accumulated in the compression phase will be released. For the mono-stiffness model where the compression and expansion phases take the same force–indentation relationship , the negative work done by the normal impulse during the interval Δ

*P*

^{(l)}of the expansion phase equals the potential energy released (the local dissipation in this model will be taken into account by the energetic constraint added at the termination of the impact). Therefore, we can directly use (2.9) to calculate the potential energy residing in the contact point and do not need to identify whether the contact point

*j*is located in the compression or expansion phase. For the bi-stiffness model, however, the work done by the normal impulse during the expansion phase needs more potential energy to be released due to the local dissipation. According to the definition of energetic coefficient and using the expression in (3.12)

_{I}, the residual potential energy at step (

*l*+1) for the contact point

*j*is(2.10)Since the potential energy at all contacts will vary continuously during the impacts, we should compare the values of the potential energy at the instant

*P*

^{(l+1)}to reselect the primary colliding point

*i*in order that(2.11)It is noteworthy that the stiff ordinary differential equation problem arising by directly using the compliant model (Hairer & Wanner 1996, §VII.7) can be avoided by this method, since we can set d

*P*

_{j}=0 if the ratio of the potential energies

*E*

_{ji}approaches zero. Moreover, the multiple compression phenomenon does not require any special treatment in the calculation of the potential energy, since its evolution is continuous during the impacts.

The termination of the impact can be identified according to the values of . Since the energetic constraint has been applied in the bi-stiffness compliant contact model, the terminating condition for that can be expressed as(2.12)For the mono-stiffness model, which is elastic, the impact at the contact point *j* will be terminated at the instant of , if no energy dissipation occurs. However, in the case of contact points with certain local dissipation, the energetic constraint defined by the coefficient *e* should be applied into the contact process. In this case, part of the potential energy should be discarded before the expansion phase finishes, and the impact at the contact point *j* will be terminated if(2.13)where *W*_{c,j} is the work done by the compression force, which is also related to the maximum value of *E*_{j} for an impact with a single compression–expansion cycle (if the impact has a multi-compression phase, the value of *W*_{c,j} will correspond to the last compression phase).

In order to improve the accuracy of the numerical results, a variable step size can be applied for searching the critical points, such that and *E*_{j}=0. If the energetic constraint is satisfied at the contact point *j*, the contact point *j* will remain open on some non-zero interval of time, and the accumulated normal impulse, as well as the relative velocities, and , will not change if separation at the contact point *j* is kept ( is always less than zero). These values can be thought of as the outcomes of the contact point *j* after multiple impacts.

However, the relative normal velocity will be influenced by the motion of the adjacent colliding bodies that still participate in the impacts at other contact points. Thus, the contact point *j* may participate in the impact again, and the value of the may change from a negative to a positive value. Based on the assumption of impacts within infinitesimal time intervals, the injected velocity for the new impact can be assigned as the value of when the separation occurs at the instant of *P*^{(l+1)}.

Let us denote the impulse instant at which the new impact appears by *P*^{(m)}, in which the contact *i* is the primary colliding point with the potential energy and the step size at this instant is Δ*P*^{(m)}. Similar to the situation at the beginning of the impact, we have to determine the possible increments of the normal impulse at the contact point *j* when it participates in the impact again. The work done by the compression force at the contact point *j* during the time interval [0,*t*] (respectively Δ*P*^{(m)}) can be approximated as(2.14)where is the increment of the normal impulse at the contact point *j* when the normal impulse at point *i* increases by Δ*P*^{(m)} during a time interval [0,*t*]. Based on the distributing law in (1.2) (or in (3.27)_{I}), and using (2.14), we have(2.15)Thus, the possible increment of the normal impulse at the beginning of the new impact can be expressed as(2.16)Once is obtained, we can modify the matrix *Γ*^{(m)} and use the difference equation (2.5) to calculate the generalized velocities , and thus the quantities related to the next moment can be obtained.

In summary, the numerical procedure can be described as follows:

use the expressions (2.1) and (2.2) to estimate the step size,

determine the initial primary colliding point

*i*and the ratios of the increments of normal impulses based on (2.3) and (2.4),begin the numerical simulation and use expressions (2.5), (2.8) and (2.9) or (2.10) to obtain the quantities of velocity and potential energy,

determine which contact point is selected as the primary point based on (2.11), and use (2.8) to determine other increments of normal impulses,

the impact at the contact point

*j*will be terminated if the expressions (2.13) or (2.12) are satisfied, and the repeating impact can be solved using expression (2.16), andthe outcomes of multiple impacts correspond to the instant when all the contact points separate from each other.

The integration is carried out with the primary impulse as the independent variable and integration step size Δ*P* estimated initially by momentum. The ‘real’ time step *h* is given by *h*=Δ*P*/*λ*, where *λ* is the force. If *λ* is large, then *h* becomes small. It is noteworthy that this has no influence on the calculation process because one does not need to come back to the *h*-integration to compute the post-impact velocities. Consequently, the algorithm can be seen as a black box representing the collision mapping when inserted in a code for the simulation of a multibody system.

## 3. The problem of Newton's cradle

The well-known problem of Newton's cradle is a typical system with multiple impacts. This system has been studied by many authors as a benchmark for multiple impacts (Brogliato 1999) and also plays an important role for many applications, such as granular material (Luding *et al*. 1994; Falcon *et al*. 1998; Nesterenko 2001; Doney & Sen 2005; Daraio *et al*. 2006; Ma *et al*. 2006). In this section, this problem is investigated by comparing the numerical results obtained from the compliant contact models (for linear springs; Wei & Liu 2006) and the experimental results presented in Ceanga & Hurmuzlu (2001) with the results obtained from the integration of the impulsive dynamics studied in §2.

### (a) Problem description

Consider a chain of balls that consists of *n* aligned balls B_{i}, so that *s*=*n*−1. Initially, all the balls except B_{1} are at rest, in contact and unstressed. The first ball B_{1} with mass *m*_{1} collides with this chain with an initial velocity *v*_{0}. Let us denote *x*_{i} as the displacement of the mass centre of the ball B_{i}, with mass *m*_{i}. When a compliant contact model is used, the contact force *F*_{i} between the balls B_{i} and B_{i+1} is(3.1)where *K*_{i} and *δ*_{i} are the contact stiffness and the elastic deformation at the contact point *i*. The exponent *η* determines the kind of contacts between the balls. The kinematic state of a contact *i* is expressed as(3.2)If *δ*_{i}<0, the balls B_{i} and B_{i+1} are separated, while *δ*_{i}>0 corresponds to the closed contact situation between B_{i} and B_{i+1}. The dynamics can be expressed as(3.3)with the complementarity conditions(3.4)Since the initial impact occurs between the balls B_{1} and B_{2}, the impulse d*P*_{1} can be first specified as the independent variable for the impact differential equations. Therefore, the distributing rule for the changes of normal impulses at the other contact points is(3.5)Initially, the distribution of the change of normal impulses with respect to Δ*P*_{1} is(3.6)In the following, the values of the velocities after the impact will be denoted with a superscript ‘+’ (as ).

### (b) Numerical results

#### (i) Case 1. Three-ball chain

Ceanga & Hurmuzlu (2001) presented an impulsive correlation ratio (ICR) for a triplet of identical balls. With the help of the energetic coefficient (denoted here as *e* since *e*_{1}=*e*_{2}), they showed that the post-impact velocities of the balls can be numerically approximated. In their experiments, four types of balls, designated as A, B, C and D with masses *m*_{A}=45, *m*_{B}=53, *m*_{C}=53 and *m*_{D}=166 g, respectively, are used. The experimental results of the ICR and *e* for the system with identical balls are presented in table 1 (see tables 1 and 2 in Ceanga & Hurmuzlu 2001).

Let us assume that the initial velocity of B_{1} is . According to the values of the ICR in table 1 and using the expressions (31)–(33) of Ceanga & Hurmuzlu (2001), we can compute the post-impact velocities of the three balls, which are thought of as the experimental outcomes reported in table 2. Let us set *η*=3/2 for the Hertzian contact and *γ*_{12}=1 due to the identical balls. A comparison between the theoretical prediction according to the method described in this paper for the mono-stiffness model and these outcomes from the ICR is presented in table 2.

Observation of table 2 shows that the theoretical results coincide fairly well with the experimental results presented in Ceanga & Hurmuzlu (2001). Clearly, the qualitative properties are comparable, since the wave effects are present in both sets of results in table 2. The discrepancies between experiments and predictions for the outcomes of post-impact velocities and the values of are very small, especially for the cases of types A and D whose coefficients of restitution take a relatively higher value.

Figure 1 shows the evolution of the velocities accompanying the normal impulses for the ball of type A, calculated with the algorithm of §2. At the beginning of the impact, d*P*_{1} is selected as the independent variable since contact point 1 contains more potential energy than contact point 2. After the end of the compression phase at contact point 1, its potential energy will decrease, while the potential energy at contact 2 will continue to increase. When the potential energy at contact point 2 is greater than the one resided in contact point 1, the independent variable d*P*_{1} for the impulsive differential equations will be changed into d*P*_{2}. In this situation, d*P*_{1} will depend on d*P*_{2}. Clearly, the ratio between d*P*_{1} and d*P*_{2} is not linear, but depends on the potential energy stored at the two contact points. Before the compression phase at contact point 2 finishes, the ball B_{1} will lose contact from the chain. Finally, all balls will separate with different post-impact velocities.

Let us investigate the difference between the theoretical predictions obtained from the bi-stiffness (see (3.1)_{I} and (3.2)_{I}) and mono-stiffness (see (2.5)_{I}) models, in which both models are set with the same exponent *η*=3/2 for the force–indentation relationships, and *e*_{1}=*e*_{2}. In the mono-stiffness model, part of the potential energy is discarded using (2.13) to confine the local energy loss, while, in the bi-stiffness model, the local energy loss is considered using different force–indentation relationships for the compression and expansion phases. As shown in table 3, the theoretical predictions obtained from different compliant contact models with the same coefficients of restitution have only small discrepancies for the outcomes of the post-impact velocities. In other words, the method developed in this paper can provide relatively precise information for an impulsive process, even though the constitutive relationship associated with the contacts is ambiguous, a common situation in practice.

#### (ii) Case 2. Five-ball chain

Complex behaviour will appear when multiple impacts occur between different types of balls. In particular, the mass and stiffness ratios will significantly influence the process. In the following, a particle chain involving five balls will be investigated using the impulsive differential equations in (1.1)–(1.4) and the compliant contact model.

Except for the imparting ball, the balls in the chain are assumed to be identical. The mass of the first ball is set as *m*_{1}=1 kg. The mass and stiffness ratios between the first ball and other balls are defined byThe contact forces are assumed to be linear elastic. This permits *e*=1 and the exponent *η*=1. The initial velocity of the imparting ball is set as *v*_{0}=1 m s^{−1}. For the case of *α*=2 and *γ*=1, figure 2 shows the evolutions of the relative velocities during impacts, using the impulsive model of §3 in part I.

Observation of figure 2 clearly shows that the compression phase between B_{1} and B_{2} (the contact between them is still kept) finishes at point a because the velocities of these two balls are the same at a, while after a. At this instant, the left force from B_{1} acting on B_{2} reaches the maximum and is greater than the right force from B_{3}, so B_{2} will continue with positive acceleration. After time a, the expansion phase between B_{1} and B_{2} will lessen the contact force on the left-hand side of B_{2}, while the compression phase between B_{2} and B_{3} will increase the contact force on the right-hand side of B_{3}. When the potential energy between B_{1} and B_{2} is completely released, separation between them will occur. Before separation, however, balls B_{1} and B_{2} will again obtain the same velocities when the normal impulse reaches the point b; a new compression process occurs between B_{1} and B_{2} after b, since one has on the right-hand side of b. This process will change the velocities of B_{1} and B_{2} and make them obtain the same value at point c. After that, the potential energy will be completely released and the contact between B_{1} and B_{2} will be finally open. An analogous process can also be found between B_{2} and B_{3}, where the instants for appear at points 1, 2 and 3, respectively.

These multiple compression phenomena can also be observed using the compliant contact model and the algorithm used in Wei & Liu (2006; shown in figure 3). The results of figure 3 were obtained with a stiffness of 1 N m^{−1}, which is not realistic; however, the outcome is independent of the absolute values of the stiffnesses, but depends only on their ratios in this example (Acary & Brogliato 2003). This implicitly means that the small deformation at contact points has little influence on the multi-impact process, confirming the validity of assumptions used in the impulsive dynamics of impact.

The phenomenon that a repeating impact occurs at the same contact point can be observed by setting *α*=1 and *γ*=2, as shown in figure 4. Clearly, ball B_{1} and ball B_{2} will collide again at point 1 on the curves. The outcomes of post-impact velocities will change due to the second impact.

It may be helpful to compare the final outcomes of the post-impact velocity obtained from the impulsive method and the one obtained from the compliant model (ii) (see part I of this paper; Liu *et al*. 2008). By changing the values of the stiffness and mass ratios, table 4 presents the final velocities of the system using the impulsive method. The results obtained from the compliant contact model are presented in table 5. Clearly, the results obtained from the two different methods agree very well. Further numerical results on *n*-ball chains colliding *m*-ball chains, and chains impacting a rigid wall, together with accurate comparisons with experimental results in Falcon *et al*. (1998) may be found in Liu *et al*. (2008).

## 4. The Bernoulli problem

The more general situation that the impacts between the balls occur in the plane is also attractive and crucial for multi-impact problems in granular material. In this section, the well-known Bernoulli example (Bernoulli 1969–1993; Brogliato 1999) is analysed.

### (a) Problem description

Two geometrically identical balls B_{1} and B_{3} are stationarily placed on a smooth plane. Let a ball B_{2} with mass *m*_{2} and initial velocity collide with the two balls along their symmetric line. The masses for B_{1} and B_{3} are *m*_{1} and *m*_{3}, respectively. The frictionless contact between the balls is assumed to be satisfied with the Hertz model.

Figure 5 depicts the system for Bernoulli's problem, in which *α* is the angle formed by the mass centres of the balls at the instant of collision. Let us set the *x*-axis along the symmetric line. *x*_{i} and *y*_{i} represent the components of the position of the mass centre for the ball B_{i}, *i*=1, 2, 3. *F*_{1} and *F*_{2} are the normal contact forces between the balls B_{1} and B_{2}, and the ones between B_{2} and B_{3}, respectively. The rotation effects of the balls are not considered. Thus, as long as the contact forces are functions of time, the governing equations for the motion of the mass centres of the three balls can be written as(4.1)where the contact forces satisfy a complementarity condition similar to (3.4). If the three balls are identical and the dissipated energy during impacts is assumed to be zero, Bernoulli (1969–1993) postulated that the impact outcome should be symmetric, and presented a theoretical solution using momentum and energy conservations. This solution reads as(4.2)Clearly, this solution significantly depends on the symmetry conditions and on the assumption of elastic collisions at both contacts. Any condition that destroys the symmetry will make the real post-impact outcomes much different from the results in (4.2). In §4*b*, the impulsive method based on the method developed in Liu *et al*. (2008) and the algorithm proposed in this paper are applied to Bernoulli's problem.

### (b) Numerical results for Bernoulli's problem

The impulsive differential equations for Bernoulli's problem can be easily obtained from (4.1)(4.3)The relative velocity along the normal direction between the balls B_{1} and B_{2} and the one between the balls B_{2} and B_{3} can be easily expressed as(4.4)The Hertz model at the contact points and the identical materials for the balls yield the following relationship between d*P*_{1} and d*P*_{2} (see (1.2)):(4.5)During the numerical simulation, either d*P*_{1} or d*P*_{2} is selected as the independent variable for equation (4.3) based on the potential energy accumulated in the contact points.

Let us set the three balls with identical mass *m*=1 kg and the coefficients of restitution in both contact points *e*_{1}=*e*_{2}=1. The initial velocity of B_{2} is set as . This is a symmetric situation where a theoretical solution could be obtained. Figure 6 presents the evolution of the velocities during impacts for the configuration with *α*=*π*/6. Clearly, the contact between B_{1} and B_{2} and the one between B_{2} and B_{3} will reach the compression point and separate simultaneously, such that the outcomes for the post-impact velocities are symmetric.

Table 6 presents the numerical results obtained from equation (4.3) and the analytical results obtained from equation (4.2) of the post-impact velocities for the system in symmetric situations with various angles *α*. The results coincide very well.

Obviously, the symmetry of the system can be preserved if both collisions between the balls are plastic impacts. The symmetric solution for the post-impact velocities should also be anticipated. This situation is illustrated in figure 7 by setting *α*=*π*/6 while the impacts at the two contact points are fully plastic (*e*_{1}=*e*_{2}=0).2 After impacts, the solution of the system is still symmetric, while the kinetic energy preserved in the system becomes 0.185 J due to local energy dissipation (its initial value is 0.5 J).

The difference of the dissipated energy at the two contact points destroys the symmetry of the multiple impact. Let us consider a limit situation, in which the collision between B_{1} and B_{2} is assumed to be elastic (i.e. *e*_{1}=1), while the collision between B_{2} and B_{3} is plastic (i.e. *e*_{2}=0, there is only a compression phase and the potential energy at the contact is dissipated completely). Figure 8 shows the evolution of the velocities during the impact for the system with *α*=*π*/3. During the compression phase, the evolution of the velocities related to B_{1} and B_{3} is symmetric, and the velocity of ball B_{2} is located in the symmetric line. However, after the compression phase finishes, the motion of ball B_{2} will diverge from the symmetric line and approaches the side of B_{3}. Before the impact between B_{1} and B_{2} finishes, B_{3} will participate in the impact again to make B_{2} move towards B_{1}, such that a little change for the velocity of B_{2} is generated. Clearly, the outcomes for Bernoulli's problem are sensitive to the configuration, the initial conditions and the properties of the contact points (Ivanov 1995).

### (c) A planar four-ball system

It is crucial for an algorithm to keep the intrinsic properties of the systems in the numerical results. In terms of Bernoulli's system, we may speculate that any system with symmetric configurations and symmetric shock conditions should provide symmetric post-impact solutions. In this section, a system composed of four identical balls is investigated using the method proposed in this paper, and the symmetry of the solutions is emphasized. Moreover, how the configuration of the system influences the transfer of momentum through a ‘distance’ impact is also discussed according to the numerical results.

Figure 9 shows the possible configurations of four identical balls stationarily placed on a smooth plane. According to the geometrical relationship, the half-angle between the lines connecting the centres of the adjacent balls should be limited to the scope of *π*/6≤*α*≤*π*/3. For case A, shown in figure 9, four contact points (*s*=4) exist in the system, while in case B, as shown in figure 10, two balls will have three contact points and the total number of the contact points is *s*=5 when *α*=*π*/6 or *π*/3.

We omit the explanation for the nomenclatures shown in figures 9 and 10, which is similar to the illustration in Bernoulli's system. The impulsive differential equations for case A can be expressed as(4.6)The relative velocities along the normal direction between the contact balls *i* and *j* are(4.7)Let us set all balls with identical mass *m*=1 kg and assume that all contacts are elastic (i.e. *e*_{i}=1 for all *i*) and satisfy the Hertzian relationship (*η*=3/2). For this symmetric configuration, we can postulate that the impact outcomes should be symmetric with respect to the *x*-axis when only B_{2} takes an initial velocity along the *x* direction. In other words, the post-impact velocities of B_{2} and B_{4} will be kept in the *x* direction, while those of B_{1} and B_{3} will be symmetric with respect to the *x*-axis. Figure 11 shows the evolution of the velocities during impacts for the configuration with *α*=*π*/5 by an impacting ball B_{2} that has an initial velocity along the *x* direction with the value of . During the impact, the velocity of the ball B_{2} changes from 1 to −0.1363, and that of ball B_{4} changes from 0 to 0.259 along the *x*-axis. The other two balls B_{1} and B_{3} have the post-impact velocities (0.4206, 0.52) and (0.4206, −0.52), respectively. It is clear that the symmetry of the solution is well preserved in the numerical results.

Based on the configuration of the four-ball system, it is obvious that the post-impact velocity of the ball B_{4} will be influenced by the angle *α*. Large values of *α* will make B_{4} obtain little kinetic energy after the system is impacted by B_{2}. For the case of *π*/4≤*α*<*π*/3, B_{4} will receive no impulse since, as a result of the integration of the Darboux–Keller shock dynamics, the relative velocities and take non-positive values during the impact, so that B_{1} and B_{3} either separate from or tangentially slide on ball B_{4}. Thus, no kinetic energy can be transferred into ball B_{4} after impacts. Figure 12 shows the numerical results for the four-ball system with *α*=2*π*/7 impacted by B_{2}, in which the ball B_{4} keeps static after impacts, and the four-ball system is consequently similar to the three-ball system of Bernoulli's problem.

If the balls B_{2} and B_{4} have the same initial velocity from two inverse directions along the *x*-axis, the four-ball system has an initial condition symmetric in both the *x* and *y* directions. So, the impact outcomes should also be symmetric in both directions. In this initial condition, figure 13 shows that balls B_{1} and B_{3} will have the post-impact velocities 0.975 and −0.975 m s^{−1}, respectively, along the *y* direction, and balls B_{2} and B_{4} the post-impact velocities 0.2224 and −0.2224 m s^{−1}, respectively, along the *x* direction.

When the angle *α*=*π*/3 or *π*/6, the four-ball system will have five contact points, as shown in figure 10. Let us consider the situation of *α*=*π*/3, in which the corresponding impulsive differential equations and the relative velocities at contacts for this configuration become(4.8)and(4.9)In comparison with the case *α*=2*π*/7, the ball B_{4} initially collides with the ball B_{2}, so that it can gain more kinetic energy. However, the post-impact velocities of B_{1} and B_{3} will decrease. Figure 14 shows the numerical results, in which the velocity of ball B_{2} before and after the impact changes from 1 to −0.197 m s^{−1} along the *x* direction, and the post-impact velocity of ball B_{4} is 0.84 m s^{−1} in the *x* direction. The post-impact velocities of balls B_{1} and B_{3} are symmetric with respect to the *x*-axis, with the values (0.178, 0.308) and (0.178, −0.308) m s^{−1}, respectively.

It is imaginable that the post-impact velocities of the balls B_{1} and B_{3} will be along the *y* direction when the four-ball system is impacted by B_{2} and B_{4} with equal inverse velocities along the *x* direction. In comparison with the configuration of *α*=2*π*/7, most of the momentum will be exchanged between B_{2} and B_{4} due to the contact point added between B_{2} and B_{4}. Figure 15 shows that the velocity of B_{2} (respectively B_{4}) changes from 1 m s^{−1} (respectively −1 m s^{−1}) to 0.964 m s^{−1} (respectively −0.964 m s^{−1}), while the post-impact velocities along the *y*-axis of B_{1} and B_{3} are 0.266 and −0.266 m s^{−1}, respectively.

## 5. Conclusions

Based on the theory of part I, a numerical algorithm is developed in this paper. The advantages of this method are

the distribution of the normal impulses can be inserted in the rigid body model so that wave effects are modelled,

the kinetic energy that is dissipated locally at the contact points can be confined by energetic coefficients, such that an energetically consistent solution is calculable,

integrating stiff compliant problems is avoided since small displacements and large contact forces are not needed,

the calculated outcomes are robust with respect to uncertainties in the parameters (stiffnesses ratios and energetic restitution coefficients), and

the algorithm is easy to implement since the integral process respects the physical meaning of the multiple impacts dominated by the primary colliding point.

Comparisons with experimental results published elsewhere and the numerical results obtained with the compliant models validate the proposed algorithm. The examples of Newton's cradle and the planar Bernoulli ball system illustrate the developments.

## Acknowledgments

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

## Footnotes

↵Equations from part I are referred to with the subscript I, such as (2.16)

_{I}for (2.16) of the first part and so on.↵The lines in figures 6–8 are parametrized with

*P*_{i}, i.e. the start corresponds to*P*_{i}(0)=0 and the end corresponds to*P*_{i}(*t*_{f}).- Received February 25, 2008.
- Accepted July 16, 2008.

- © 2008 The Royal Society