Modelling multiple-simultaneous impact problems with a nonlinear smooth approach: pool/billiard application

Smooth approaches are able to model reasonably well contact/impact events between two bodies, showing some peculiarities when dealing with certain geometries and arising certain issues with the detection of the initial instant of contact. The characterization of multiple-simultaneous interaction systems, considering (or not) energy dissipation phenomena (mainly friction), is always an interesting research topic, addressed from different perspectives. In the present work, the process of design, optimization and verification of a multiple-impact, day-to-day multibody novel model is shown. Specifically, we have decided to focus on a pool/billiard game due to its geometry simplicity. The model involves several balls moving freely and rolling, suffering different kinds of contacts/impacts among them and against the cushions and the cloth. In this system, the proper modelling of both contact and friction forces in the multiple, simultaneous contacts and impacts events is critical to obtain consistent results. In addition, these forces are complicated to model because of its nonlinear behaviour. The different existing approaches when dealing with multiple-contact events are briefly described, along with their most distinctive features. Then, the interactions identified on the model are implemented using several nonlinear contact-force models, following a smooth-based approach and considering friction phenomena, aiming at determining the most suitable set of both contact and friction force models for each of these implemented interactions, which take place simultaneously, thus resulting in a complex system with multiple impacts. Subsequently, the solving method that provides the most accurate results at the minimum computational cost is determined by testing a simple shot. Finally, the different interactions on the model are verified using experimental results and previous works. One of the main goals of this work is to show the some of the issues that arise when dealing with multiple-simultaneous impact multibody systems from a smooth-contact approach, and how researchers can deal with them.

Abstract Smooth approaches are able to model reasonably well contact/impact events between two bodies, showing some peculiarities when dealing with certain geometries and arising certain issues with the detection of the initial instant of contact. The characterization of multiple-simultaneous interaction systems, considering (or not) energy dissipation phenomena (mainly friction), is always an interesting research topic, addressed from different perspectives. In the present work, the process of design, optimization and verification of a multiple-impact, day-to-day multibody novel model is shown. Specifically, we have decided to focus on a pool/billiard game due to its geometry simplicity. The model involves several balls moving freely and rolling, suffering different kinds of contacts/impacts among them and against the cushions and the cloth. In this system, the proper modelling of both contact and friction forces in the multiple, simultaneous contacts and impacts events is critical to obtain consistent results. In addition, these forces are complicated to model because of its nonlinear behaviour. The different existing approaches when dealing with multiple-contact events are briefly described, along with their most distinctive features. Then, the interactions identified on the model are implemented using several nonlinear contact-force models, following a smooth-based approach and considering friction phenomena, aiming at determining the most suitable set of both contact and friction force models for each of these implemented interactions, which take place simultaneously, thus resulting in a complex system with multiple impacts. Subsequently, the solving method that provides the most accurate results at the minimum computational cost is determined by testing a simple shot. Finally, the different interactions on the model are verified using experimental results and previous works. One of the main goals of this work is to show the some of the issues that arise when dealing with multiple-simultaneous impact multibody systems from a smoothcontact approach, and how researchers can deal with them.
Contact forces and impact events, which are present in almost all branches of engineering [17][18][19], are responsible for the appearance of harmful phenomena in mechanical systems such as vibrations [20], wave propagation [21], fatigue [22], wear [23] and crack [24]. In impact events, sudden changes occur, in which the conditions of the mechanical system vary in very short times. This implies the appearance of great magnitude forces, energy dissipation processes and discontinuities of velocities and accelerations, among other issues. Contact events are difficult to model and pose a challenge for the engineers due to the large number of variables that must be considered: geometry of the contacting surfaces, material properties, inclusion of friction phenomena, multiple-simultaneous impacts, … Two main approaches are identified when modelling contact events: the non-smooth approach and the models based on contact forces [25,26].
One of the main purposes of this work is to develop a realistic multibody model of a pool/billiard game that takes into account all contact interactions of the system, thus properly predicting the behaviour of the different bodies (namely the balls, the cloth, the cushions and the pockets), given a certain shot. Unlike most of the researches developed so far [27][28][29], a smooth approach has been considered to define interactions between bodies [30]. In this approach, normal contact forces are expressed as continuous functions of the relative indentation between the contacting bodies, as well as their geometric and material characteristics [31]. Parameters from previous experiments have been used to calculate the values of the contact forces, for the different interactions identified.
The equations of motion have been implemented in a Matlab code designed to perform and solve forward dynamic analysis, which allows multiple customization choices as well as a good control of the integration process and provides accurate results at a reasonable computational cost.
Several simulation scenarios corresponding to different shots have been considered. In the following sections, the chosen contact force models for the different interactions defined in the model (namely ball-cloth, ball-ball and ball-cushion) are discussed, considering some features such as the computational efficiency and the contact detection. The values of the force model parameters are remarked too [32]. Friction forces are also implemented using different models with the purpose to appraise the most relevant and suitable options [33][34][35].
2 Formulation used to perform dynamic analysis on multibody systems The formulation of the equations of motion for constrained multibody adopted in this work follows the nomenclature defined by Flores and Machado in [30,36,37], respectively. It follows closely the work developed by Nikravesh, who used generalized Cartesian coordinates to describe the configuration of the system [2].

Position and orientation characterization
The position of the mass centre of a body i is located by vector r, which consists of three independent variables associated to Cartesian coordinates xyz Any point of body i can be described from the origin of the local reference frame (usually located in its mass centre) by its vector s P i , so its location with respect to the global system can be expressed as Vector s P i , defined in global-system coordinates, can also be described with respect to the local reference system, being then denoted as s 0P Being A i the rotation matrix that describes the orientation of the local coordinate frame with respect to the global system. For a 3D system, A is defined as [36] A ¼ where w, h and r are the so-called Euler angles. Letters s and c refer to sine and cosine trigonometric functions, respectively. Euler angles represent three independent, consecutive rotations that define the orientation of a given body-fixed coordinate system with respect to a global one. As matrix multiplication is not commutative, a convention when rotating the successive intermediate reference systems must be considered, being this usually zxz [36]. However, this formulation presents some singularities that make it impossible to implement it in an engineering software. For example, when sin h ¼ 0, intermediate systems defined by angles w and r are indistinguishable. For this reason, and based on Euler's theorem on finite rotation, the fourth coordinate is added, so a rotation in the three-dimensional space can always be described by a rotation along a certain axis over a certain angle /. These fourth magnitudes, denoted by e 0 , e 1 , e 2 and e 3 are known as Euler's parameters, and must meet the following conditions: So, these four parameters, along with the three Cartesian coordinates that locate the position of the mass centre, describe uniquely the position and orientation of a certain body in a multibody system: The kinematic joints of a constrained multibody system, usually considered as holonomic, are defined through a set of algebraic equations which are functions of the generalized coordinates: being t the time variable. The first derivative of Eq. (2.9) with respect to time yields the velocity constraint equations, whereas the second one leads to the acceleration constraint equations: where D is the Jacobian matrix of dimension k Á n, being k the total number of constraints of the system and n the number of coordinates. Term À _ Dv is usually placed on the right-hand side of the acceleration equations, thus being denoted as c and agglutinating those functions explicitly dependent on positions, velocities and time: constraints represented in Eq. (2.9) are nonlinear in terms of q and are, usually, solved iteratively by employing the Newton-Raphson method.
On the other hand, Newton-Euler equations of motion of an unconstrained system are defined as [36]: where € q is the acceleration vector and M is the system mass matrix, written as ð2:20Þ This system is solved for € q and k to obtain the dynamic evolution of the system. In each integration time step, the accelerations vector € q together with the velocities vector _ q are integrated to obtain the positions and velocities for the next time step. This algorithm is repeated until the final analysis time reached.

Kinematic joints and Newton-Euler equations
of motion: constraint violation However, the system described in Eq. (2.20), known as the standard Lagrange multipliers method, does not consider explicitly the position and velocity equations related to the kinematic constraints. This leads to a violation of the constraint equations when dealing with moderately long simulations, due to the error accumulated over the integration process and/or to inaccurate initial conditions. In order to avoid this error propagation and eliminate any possible slip in the position and velocity equations throughout the analysis, or at least keep such errors within a defined margin, different methods have been defined. In this work, three different error control algorithms have been considered: the Augmented Lagrangian formulation [39], the stabilization method defined by Baumgarte [40] and the direct correction method proposed by Marqués et al. [41]. As will be seen in the following sections, the proper choice of a solving method has a huge impact on the results obtained. The augmented Lagrangian formulation consists of solving the system equations of motion by an iterative process with the goal of penalizing constraint violations. Being i the i th iteration of the integration process [39], the iterative scheme evaluate the system accelerations using the expression [36] M where a is a large real value (penalty number) and l and x are parameters related to the natural frequency and the damping ratio of the penalty system, respectively (mass, dashpot and spring). _ U represents the velocity constraint equations. The iterative process continues until the following condition is met being e a specified tolerance defined by the user. Even when the system is close to a singular position or when dealing with redundant constraints the system of equations can still be solved. Baumgarte proposed, in a similar way as the augmented Lagrangian, an alternative expression for Eq. (2.10) [40] Equation (2.23) is a differential equation that leads to a closed-loop system in terms of kinematic constraint equations. The two added terms, 2 Á a Á _ U and b 2 Á U, operate as control terms to damp the violation of acceleration constraint by inserting the data of both position and velocity constraint violations. If both parameters a and b are chosen as positive constants, the stability of the system is assured. Baumgarte proved that values of these parameters smaller than 10 are advantageous. Nonetheless, he stated that the suitable choice of the constant should be obtained by numerical experiments. A study carried out in [42] and based on the stability analysis procedure used in control theory shows that, for a fixed time step, the values a = b = 5 lead to a convergence of the integration process without oscillation, being this process even faster if a = b = 20. However, this implied a stiffer system. For certain problems, small parameters are preferred, whereas for others, large values provide better results [43], these being always in the interval 1-20 [39]. A bigger increase in b compared to one of a produces oscillations of the results and, ultimately, the divergence of the integration process. If both values are equal, the critical damping is reached. For the model developed in the following sections the values a = b = 5 were chosen. Unlike Lagrangian formulation, Baumgarte's method uses to fail when dealing with kinematic singular configurations or with redundant constraints [39].
Marqués et al. proposed a direct correction approach to get rid of constraints violation by considering the positions vector at a ith iteration as where q c represents the corrected positions of the bodies of the system at that time step, q u denotes the uncorrected positions and dq represents the set of corrections that eliminates the constraints violation. Thus, Eq. (2.9) became where the last term, dU, denotes the variation of the constraint equations Combining Eqs. (2.25) and (2.26) and considering the concept of Moore-Penrose inverse matrix, dq can be rewritten and a new definition of q c is obtained In a similar way, the vector of generalized velocities can be corrected as With these two expressions, positions are corrected at each time step until deviation is contained. Then, by using Eq. (2.29) with the new values of q c , the velocity constraints violation is corrected.

Definition of the application model: pool/billiard table
In this section, a practical application of the equations developed above will be described. The authors decided to choose a pool/billiard multibody due to the different, multiple-contact interactions that take place simultaneously (see Fig. 1). In addition, it was not excessively difficult to model it, as geometries are relatively simple, and there are experimental and theoretical data available to verify the model and draw some worthwhile conclusions about the issues arisen during the characterization and verification stages. As will be seen, the proper definition of both, contact and friction forces, will have a critical impact on the results.

Contact modelling approaches and contact interactions considered in the proposed model
In general, two different approaches are identified when modelling contact/impact phenomena in multibody systems: the non-smooth or impulse-momentum based approach and the models based on contact forces or force based models [30,32,44,45]. The former assumes that the deformations experienced by the bodies involved in the process are small, compared to their geometry (so they can be considered as rigid solids), and that the contact process is almost instantaneous (so the forces involved are modelled through the concept of impulse) [5,46]. Some authors subdivide this first group into two classes [47,48]: algebraic methods that relate post and pre-impact velocities through a certain function [49,50], which can be explicitly or implicitly defined; and first-order dynamics models that follows the Darboux-Keller approach [51,52]. The first group consider the contact/ impact event as instantaneous, whereas the second one establishes a small-time interval. Emphasis must be given in this second group to the LZB model [48], an extension of the aforementioned DK approach applied to single impacts, which considers plastic deformation through a bi-stiffness contact model and energetic coefficients of restitution [53,54]. A comparative analysis between these two groups can be found in [55]. Non-smooth models allow a simple, computationally efficient modelling of the contact event. However, some drawbacks are the requirement of different numerical strategies for certain contact scenarios, such as a permanent contact or an intermittent one [5,30] or the complexity of some of the algorithmic procedures used to model the impact event [37,56,57]. Some studies regarding the specific contact scenarios issues can be found in [58,59].
On the other hand, models based on contact forces, also known as penalty models or compliant models [32], are characterized by their simplicity, computational efficiency and implementation easiness [37,60]. These approaches define contact forces as continuous functions of the relative penetration (and its temporal derivative) of the contacting bodies [61,62], which are supposed to be deformable [30]. The contact regions of each body are modelled as a set of spring-damper elements scattered over their surfaces. Nevertheless, these approached present some drawbacks, being the main one the proper choice of the parameters associated with the definitions of the forces (for example, the equivalent stiffness or the degree of nonlinearity of the indentation) [32,63]. In these models, the right detection of the initial instant of contact is critical to obtain accurate results [64].
A quick search on the subject shows that most of the works made so far about pool ball-ball interaction are based on non-smooth approaches [27-29, 65, 66]. The main aim of this research is to develop a complete model penalty-based that considers not only ball-ball interactions, but also the ball-cloth and the ballcushion ones, so multiple, different contact events take place simultaneously. Within the spirit of this methodology, special attention is paid to the contact detection itself, in terms of both the consistency of the results and computational efficiency [64].
There are multiple commercial software in the market focused on multibody system modelling, both commercial [67,68] and non-commercial [69,70] However, in this project equations will be implemented and solved in a Matlab-based code designed to perform forward dynamic simulations for spatial multibody systems [36]. Bodies will be defined by their relative coordinates, using the Euler parameters to describe the rotation of their respective local reference systems with respect to the global one. The translational motion will be described in terms of the Cartesian coordinates.

Material characteristics of the model and geometric modelling of the contact interactions
An Eight-ball table, the most popular pool billiard game, was considered to develop the model. It consists of 15 coloured balls (seven solid-coloured balls following the sequence yellow-blue--red-purple-orange-green-brown, seven striped balls with the same sequence and the black 8 ball) and a white ball, which is hit by a cue. According to the regulations established by the WPA (World Pool Association), a 9-Foot table was chosen, with dimensions of 2.54 9 1.27 m.
Although the cushions should be triangular-shaped bodies, the ball-cushion interaction was modelled as a sphere-plane contact, so they would be considered as plane bodies. The balls were also defined according to the regulations (57.15 mm diameter, 170 g mass, solid bodies made of Phenol-formaldehyde resin). Table 1 shows the values of the elastic properties of the respective materials. Three contact interactions are identified in the proposed model: ball-ball, ball-cloth, and ball-cushion. Some geometric data of the bodies must be provided to model the impact events and to calculate the contact detection condition, the contact direction and the value of the indentation [77]. For the ball-ball interaction, the radii of the balls must be defined, the contact condition being: where d is a vector set between the centres of the balls i and j (see Fig. 2a). The contact between the balls and the cloth is defined as a ball-plane interaction, for which the radius of the ball and a normal vector of the plane must be provided. The contact condition is calculated as: ð3:2Þ being dist ij the minimum distance from the centre of the ball to the plane (Fig. 2b). The interactions between the balls and the cushions are defined as sphere-parallelepiped interactions, which are specific cases derived from the sphere-plane general contact case. Each cushion is divided into three half-bodies, as shown in Fig. 3, thus being a total of 18 planes, as well as the cloth. Therefore, each cushion consists of a central, rectangular prism to which two triangular prisms are attached at its bases. The corners of the cushions are critical contact points, where three different situations can take place (see Fig. 4).
Depending on the specific contact point, the ball will rebound on the direction normal to the longest side of the cushion (Fig. 4a), the original direction (Fig. 4b) or the normal to the shortest side of the cushion (Fig. 4c).
To define accurately the boundaries of each plane, the dimensions of the prisms are implemented through their inertia matrixes, their terms being functions of these, as well as their mass [78]. If any of the coordinates of the projection of the sphere is out of the boundaries of the plane (s C ! jÀsph;x AE R and s C ! jÀsph;y AE R in Fig. 5), contact will not happen.

Contact and friction forces
An appropriate set of contact and friction force models had to be chosen for each of the contact cases described above, considering the properties of the contacting bodies and the contact event itself, as well as looking for a balance between consistency of results and computational efficiency. For this reason, a research on the modelling of the different pool/billiard interactions was carried out [34,79].
The model defined by Lankarani and Nikravesh was selected to calculate the normal force between the balls [80]. This approach stands out by its reasonably good results when dealing with general mechanical contacts, especially when dealing with cases in which the energy dissipated during the contact is relatively small compared to the maximum absorbed elastic where K denotes the contact stiffness parameter, which, for a contact between two spheres, i and j, with radii R i and R j , respectively [32], is where R eff is the equivalent radius. r i and r j are material parameters given by being v l and E l the Poisson's ratio and the Young's modulus of each sphere, showing the influence of the material properties in the contact event. For the contact between a sphere i and a planar surface j [32] v denotes the hysteresis factor that quantifies energy dissipation, which for LN model is defined as [32]: being c r the restitution coefficient and _ d À ð Þ the initial impact velocity. Regarding the friction force, two models were compared: Threlfall's model, chosen due to its simplicity and the data availability, (see Fig. 6a) and an experimental expression defined by Alciatore (see Fig. 6c) [34], based on Marlow's data [35], which is dependent on the relative tangential contact velocity. The latter is defined by Eq. (3.8) l v ð Þ ¼ 9:951 Â 10 À3 þ 0:108 Á e À1:088Áv ð3:8Þ On the other hand, Hertzian contact model provided the best results to characterize the ball-cloth interaction, due to its simplicity and computational efficiency [31,32] As there is an almost-permanent contact between the balls and the cloth, for the shots tested in this work, energy dissipation due to elongation can be omitted, thus making this approach a perfect choice. However, some considerations regarding this interaction will be commented in the following section. The model proposed by Ambrósio was used to define the friction phenomena associated to this contact event (see Fig. 6b) [33]. This model, as well as Threlfall's, set up a null friction force for very low tangential velocities. As it will be commented later, balls will move at much higher speeds than the limits defined, so they will tend to slip rather than roll. Therefore, stickslip motion can be neglected in this model: when the cue ball is hit, the friction between the balls and the cloth is so low that they start to slip immediately. Although the models described in Fig. 6 are 1D approximations, the friction forces will be 3D vectorial quantities, as they are obtained from the normal forces whose direction is defined by the normal vector calculated through the algorithms described above. Lastly, the model defined by Lankarani and Nikravesh model was chosen again to characterize the interaction between the balls and the cushions, as cushions dissipated most of the impact energy through elastic deformation, being the impact almost fully elastic. Friction was considered as minimal in this contact event at first, so it was neglected. However, Mathavan et al. calculated, based on experimental results, some specific values of restitution and friction that were later considered to validate the ball-cushion contact interaction, being these 0.98 and 0.14, respectively [79]. The friction force was implemented then using again the model defined by Ambrósio [33].
The rest of values of both friction and restitution coefficients tested, for the three interactions described above, were extracted from the researches carried out by Alciatore [34] and are shown in Table 2. For the contact force expressions, a value of the nonlinear power exponent n = 1.5 was adopted when computing the contact forces values.
For both the Ambrósio and the Threfall models, some values of tolerance velocity (v 0 and v 1 , see    Fig. 6) had to be defined [33]. In this work, a value of v 0 = 0.01 m/s was taken for the Threlfall model, while for the Ambrósio model values of v 0 = 0.0001 m/s and v 1 = 0.001 m/s were adopted. Both models show similar behaviours to Coulomb's law and provide computationally efficient results. In addition, for tangential velocity values close to zero, the friction force will always be low, independently of the displacement.

Initial conditions for the verification experiments
Three experiments were carried out with the proposed model: firstly, a fine-tune process was made, along with the comparison of four solving methods of the Newton-Euler equations [36]; secondly, ball-ball and ball-cloth interactions were verified simultaneously using Alciatore's experimental results and theoretical work [34]; finally, the interaction between the balls and the cushions was checked through a straight shot compared with both experimental results and a slowmotion video of a real shot. For all cases, the inertia of the balls was J = 0.000055 kg m 2 . All bodies had their local coordinate system aligned with the global frame, so their Euler parameters were (1, 0, 0, 0). The centre of the cloth was set at the position (0, 0, 0) m of the global reference system.
Four balls were placed on the cloth for the first experiment: the cue ball and the first three balls of the rack (see Fig. 7a). The shot was implemented as an initial velocity of the white ball of value 10.729 m/s (24 mph) in the positive direction of x-axis. This value corresponded to a professional break shot [81]. All balls were able to move freely. The initial position of the centre of the white ball was (-0.635, 0, 0.028576) m, while the centre of the first ball of the rack was located on (0.635, 0, 0.028575) m. Some comments will be made on the effects of the separation of the static balls on the results, and how the different contact events that happen simultaneously affected each other.
For the second experiment, the verification of ballcloth and ball-ball interactions, two balls were arranged on the cloth: the cue ball and a blue solidcoloured one. Two cases were considered: a medium speed shot and a fast speed shot, with values of the initial velocity of the cue ball of 3 and 7 mph, respectively [34] (see Fig. 7b, c). For the first shot, the centre of the cue ball was located on the global position (0, 0, 0.028576) m, while the local system of the blue ball was on x = 0.635 m. On the other hand, for the second case the cue ball was located on (-0.635, 0, 0.028576) m and the blue ball on x = 0 m, preventing it from colliding with the cushion on x = 1.27 in excessively early stages of the simulation.
Regarding the third experiment, the two sets of contact and friction force described above for the ballcushion interaction were tested and verified. First, a simple, straight shot in the positive direction of the yaxis was compared with the experimental results from real slow-motion footage. The cue ball had an initial velocity value of 2.724 m/s, and its initial position in global-system coordinates was on (1, 0.55, 0.028576) m. Then, the values proposed by Mathavan et al. were implemented and verified, testing different impact velocities on x-axis direction (see Fig. 6d, e).
It must be highlighted that no complex plays were considered in the early design and verification stages. Professional pool/billiard players consider additional effects such as Coriolis effect in order to make unbelievable moves that allow getting multiple balls into the pockets with a single shot [82]. More information regarding the so-called English effect in pool/billiard can be found in [34] (see Fig. 8).

Numerical results
The dynamic analysis was carried out using four different solving methods: the standard Lagrange multipliers [2], the Augmented Lagrangian formulation [39], the stabilization method defined by Baumgarte [40] and the direct correction method proposed by Marqués et al. [41], all of them run by the ordinary differential equations (ODE) integrator of Matlab. For the first experiment, analysis times of 3 s and 4 s were set, with different values of default time step. Conversely, for the second test an analysis time of 2 s was defined, with a 0.0001 s time differential by default. Finally, for the third experiment, an analysis of 0.2 s was carried out to recreate the real shot, with a default time step value of 0.0001 s, whereas for the experimental values verification, the analysis time was set at 1 s, with a time step of 0.001 s. These values of time steps would be reduced automatically by ODE integrator at the critical moments, mainly when contact starts, in order to perform a proper and accurate detection of the initial instant of contact and the value of indentation then [64].

First experiment: model optimization and solving methods comparison
In the early stages of the development of the model, balls static at the outset (coloured spheres in Fig. 7a) were arranged with a certain distance between their surfaces, avoiding the system to stall in excessively early times of the simulation. The coloured balls were always arranged so their mass centres formed an equilateral triangle. The distance between the mass centres of the coloured balls were defined by the expression d ¼ 2 Á Radius þ a, being a a value set by the user, as shown in Fig. 9. The positions tested are resumed in Table 3. The last value corresponds to a random arrangement of the balls, whereas the first four values were set to determine the minimum precision of the model. These values proved to have a decisive impact on the results obtained. For instance, with lower values of a, some of the balls (mainly the cue and the yellow ones) tended to increase their position in global z-axis as the analysis progresses and they rolled on the cloth (see Fig. 10). As a was increased, the maximum values reached by these balls reduced, providing a more reasonable behaviour of the model, i.e. balls do not tend to gain energy and increase their z-axis position. For the random value set, cue and yellow balls hardly move in z-axis. This is also an example of how a certain configuration of one of the contact interactions affects to the rest of contact events that happen simultaneously. For this reason, in the first simulations performed with this model, when it was being tuned, both normal force in the ball-cloth contact and gravity force were neglected, therefore keeping constant z-axis component. This resulted in a substantial reduction of the This modification of the model had some consequences: friction forces normal to the cloth had to be neglected to avoid the balls to either fly over the cloth or to sink into it. The value of the normal force was calculated as the product of the ball weight by the value of acceleration of gravity, 9.81 m/s 2 . In Fig. 11 a comparison of the evolution of the position of the yellow ball along x-axis is shown, for both the ''first'' version of the model and this second, ''fine-tuned'' one. The first model was simulated using Baumgarte's stabilization method, whereas the second one was tested with three of the methods introduced above: Baumgarte's, standard Lagrange multipliers and Augmented Lagrangian formulation. As can be observed, for the lower values of a, the model provided consistent results, as the ball reached the boundaries faster since there was no energy dissipated in z-axis movement. However, for the greatest values of a, 10 -6 and 2.1325735 9 10 -4 m, different results from the previous ones were obtained. In the case of a = 10 -6 m, the reason is that the yellow ball impacted with both the red and the blue ones at t * 0.7 s. On the other hand, for the random value of a, the yellow ball did not move in positive x axis direction after the first impact with the cue ball, but in the negative one. This was later checked in a nonsmooth model and proved to be wrong. In consequence, the value of a that should be used to obtain proper results should be in the interval 10 -9 -10 -7 .
Once the model was tuned, the different solving methods were compared in order to define which one of them provided the most accurate results at the minimum computational cost. For this experiment, the random value of a was chosen, so the computing times were not too much long. Furthermore, gravity and forces normal to the cloth were considered. The first simulation was carried out with a default time step of 0.01 s. As can be seen in Fig. 12, both standard and Baumgarte methods provided logical results, as the cue ball remained along the global x-axis throughout the simulation (i.e. the y-position component kept null) for an initial shot with only not-null x-component, while the Augmented Lagrangian and the direct correction methods showed that the cue ball moved on the global y-axis. This was confirmed by the values of y-axis linear acceleration in Fig. 13.
One of the reasons of these incongruous results could be the large size of the default time step, which is reduced by the ODE integrator at critical instants (for example, when a contact starts), but that it is kept constant as long as the error is within the allowed margin for the rest of the analysis. For this reason, the test was repeated, this time defining a default time differential value of 0.00001 s. The new graphs of both y-axis position and y-axis linear acceleration of the cue ball are shown in Figs. 14 and 15, respectively. On this occasion, the deviations were way lower, and, for the augmented Lagrangian method, the symmetry of the problem was kept. Nevertheless, the direct correction method still showed a slight deviation of the cue ball on y-axis, although it could be considered negligible. Going over the computing times shown in Fig. 16, it can be said that the Augmented Lagrangian formulation and the direct correction method are not computationally efficient for the proposed model. This a perfect example of how the choice of the solving

Ball-cloth and ball-ball interactions verification
In order to verify completely the proposed model, the results obtained by Alciatore, mentioned above, were used [34]. Alciatore defined the ''308 Rule'', which stated that for a rolling-CB shot, over a wide range of cut angles, between a 1/4-ball and 3/4-ball hit, the cue ball will deflect by very close to 30°from its original direction after hitting the objective ball (see Figs. 17 and 18 for further clarification). He proposed the following expression to calculate the deflection angle h as a function of the fraction hit f :  [33,34]. In addition, two different shots were tested: a medium speed shot (3 mph) and a fast speed shot (7 mph), both of them with only not-null x-axis direction component [34]. The results obtained for both cases are shown in Fig. 19. It is observed that, for a medium speed shot, with the same cue ball speed as that used by Alciatore in his experimental tests (3 mph), the deflection angle values provided by the model were quite similar, for both friction force models. However, with the velocity-dependent model a slight deviation was noted, for low ball hit fraction values, with respect to the experimental results, except for the augmented Lagrangian method, which had been discarded before due to its lack of computational efficiency. On the other hand, it is observed that, as the cue ball speed increased, the deflection angle also did, for the same ball hit fraction value, which is consistent from a physical point of view, as kinetic energy of the balls was greater. For both cases, the maximum values of deflection angle were obtained for values of f between 0.7 and 0.9.

Ball-cushion interaction verification
With the aim of completing the verification process of the proposed model, two experiments were carried out to verify the remaining interaction: the ball-cushion contact. Both tests provided quite promising results and allowed the model consistency verification process to be completed. Firstly, some footage of a billiard ball impacting against the cushion were analysed in order to replicate the shot shown in the video in the designed model. According to the author of the video, the shot was recorded with a high-speed camera, with a frequency of 2000 frames/s. Calculating the initial speed of the ball using two alternate frames gave a value of  Fig. 20 were obtained. The proposed model provided results almost identical to the real shot from the video, thus verifying this third interaction.
After this test, the authors found the work developed by Mathavan et al. [79], in which they obtained the values of coefficient of restitution and ball-cushion sliding friction coefficient that provided the minimum error for the ball-cushion interaction. These values were c r = 0.98 and l = 0.14, for a shot in which the cue ball moved normal to the cushion, for different initial speeds. These shots were replicated in the proposed model, implementing the values of restitution and friction mentioned. The friction force was simulated using again the model defined by Ambrósio [33], obtaining the results shown in Fig. 21. As can be observed, for low incident velocities, both models complied reasonably well with the experimental results, whereas for greater values, specially from 2.5 onwards, there was a slight deviation with respect to the experimental results. Nevertheless, it can be said that the proposed model provides once again consistent, realistic results.

Concluding remarks
In the present paper a multibody model of a pool billiard game with multiple impacts has been designed, developed and verified. Taking an innovative contact force approach with respect to previous works, the different contact and friction interactions identified (namely the contact events between the balls, between the balls and the cloth and ball-cushion contact) have been implemented. The suitable set of contact and friction models, for each one of these interactions, has been defined, looking for a balance  Several experiments have been carried out with the proposed model. In the first test, a process of optimization has been made to reduce the computing time and to fine-tune the results provided by the model. Then, different solving methods of the equations of motion have been tested, identifying those that provide consistent results at a minimal computational cost. Then, the ball-ball and ball-cloth interactions have been verified simultaneously using both theoretical and experimental results from previous works. Finally, the interaction between the balls and the cushions has been validated using both footage of a real pool billiard shot and experimental results. Both verification processes have delivered consistent, promising results.
As can be seen, the study of the phenomenon of contact is a research field that covers a wide range of applications. Considering the most characteristic (a) Constant friction coefficient (0.06), implemented through a Threlfall model [33] (b) Velocity dependent friction coefficient defined on Equation (3.6) [34] Fig. 19 Verification results for the ball-ball interaction features of each interaction (elasticity, energy dissipation mechanisms and friction phenomena), the most suitable set of contact force and friction force models has been defined for each interaction. For the interaction between balls and between the balls and the cushions, a contact force model conceived to characterize almost elastic impacts as that defined by Lankarani and Nikravesh has been chosen. Conversely, in order to get as much computational efficiency as possible, a simple Hertzian model has proved to be the best choice to define the permanent contact between the balls and the cushion, for the tested shots. It should be noted that, in complex multibody models with different and simultaneous contact interactions as the one developed in this work, multiple factors must be taken in order to obtain consistent results. The arrangement of the bodies, distances between them and contact and friction  Comparison of the simulation progress for the ball-cushion interaction experiment models, as well as time step differential have a critical impact in the results. Besides that, penalty-based approaches have proved to work reasonably well with simultaneous-multiple-impact multibody systems.
Acknowledgements The authors would like to thank the Spanish Government through the MCYT Project ''RETOS2015: sistema de monitorización integral de conjuntos mecánicos críticos para la mejora del mantenimiento en el transportemaqstatus''. The authors would also like to acknowledge the support received by the Community of Madrid through its multiyear agreement with University Carlos III focused on its policy ''Excelencia para el Profesorado Universitario''.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.
Data availability The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

Declarations
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.