An Analytical Solution for Non-Linear Viscoelastic Impact

The paper presents an analytical solution for the centric viscoelastic impact of two smooth balls. The contact period has two phases, compression and restitution, delimited by the moment corresponding to maximum deformation. The motion of the system is described by a nonlinear Hunt–Crossley equation that, when compared to the linear model, presents the advantage of a hysteresis loop closing in origin. There is only a single available equation obtained from the theorem of momentum. In order to solve the problem, in the literature, there are accepted different supplementary hypotheses based on energy considerations. In the present paper, the differential equation is written under a convenient form; it is shown that it can be integrated and a first integral is found—this being the main asset of the work. Then, all impact parameters can be calculated. The effect of coefficient of restitution upon all collision characteristics is emphasized, presenting importance for the compliant materials, in the domain of small coefficients of restitution. The results (variations of approach, velocity, force vs. time and hysteresis loop) are compared to two models due to Lankarani and Flores. For quasi-elastic collisions, the results are practically the same for the three models. For smaller values of the coefficient of restitution, the results of the present paper are in good agreement only to the Flores model. The simplified algorithm for the calculus of viscoelastic impact parameters is also presented. This algorithm avoids the large calculus volume required by solving the transcendental equations and definite integrals present in the mathematical model. The method proposed, based on the viscoelastic model given by Hunt and Crossley, can be extended to the elasto–visco–plastic nonlinear impact model.


Introduction
The behavior of a mechanical system is described by using ordinary differential equations or differential equations which contain partial derivatives. However, there are few cases for which these equations are linear. Most of the cases involve nonlinear equations, and only the strict local study of the comportment of the system allows the linearization of the equations.
Concerning the domain of dynamics of mechanical systems, the main causes conducting to the occurrence of nonlinearities are as follows: - The impact phenomena, characterized by the sudden variation of the kinematical parameters of the system [1]; - The frictional phenomena: dry friction and rolling friction.
The requirement of precise modeling of the phenomenon of impact for dynamic systems imposed the use of modern mathematical techniques [2,3].
The kinematics discontinuities assume considering the time parameter in order to describe the evolution of the geometrical and kinematical parameters. The sudden variation in time of the geometrical and kinematical parameters of a mechanical system leads to the incidence of forces of considerable intensities. The phenomena characteristic to the sudden variation of geometrical and kinematical parameters have been acknowledged in technical literature as percussions.
Two bodies may act together via two manners: by the interaction of a field and by direct contact. The latter interaction method is especially accustomed in current activities and the good reason for this affirmation is immediate since we can encounter from closing doors and windows to landing of airplanes. However, bringing into contact two bodies intended to interact assumes the existence of relative motion between them. After the two bodies come into contact, the state of relative motion changes. The transition from a kinematical state to another is accomplished during a finite period of time, longer or shorter, and thus leads to the development of interaction forces having greater intensities with the shorter duration of the interaction.
So is explainable the destructive power of the firearms, power which results not from the value of the mass of the projectile, that can be of order of grams, but from the short period of time during which the motion state of the projectile is modified. Consequently, a well-defined category of mechanical phenomena, named percussion, collision or impact, having as common characteristic the development during a very short period of time is accepted. The study of the impact phenomena is particularly contemporary and the numerous treatises, monographs, journals dedicated to the subject sustain this affirmation.
There are two manners to approach the collision mechanics problem [4]. The first manner is by stereo dynamics [5][6][7][8][9][10][11][12][13]. According to this method, the bodies are considered perfectly rigid and the collision is instantaneous. Among the characteristic impact parameters, undoubted the most important is the coefficient of restitution, c r . This is defined as the ratio, with changed sign, between the projections on the normal direction, n, of the relative velocities post-impact, v , and pre-impact, v [14][15][16][17].
The requirement of obtaining an almost faithful model imposed the consideration of friction in the contact regions. Keller [18] proposes a model for the study of dry friction collisions but his model is difficult to apply. For the planar collisions with dry friction, the paper of Wang and Mason [10] represents a reference; for the study of dry friction impact the Wang and Mason use the Routh [19] percussion' plane. Using the Routh's percussions plane one can explain the paradox reported by Kane [20]. Kane considers a double pendulum subjected to impact with a rough horizontal plane.
To break the deadlock, the hypothesis of rigid colliding bodies hypothesis must be abandoned. Thus, the initial moment t i of the collision is the instant when the first points belonging to the two bodies come into contact. The moment when the distance between the initial points of contact reaches the maximum value is denoted t c . At this moment t c , the relative velocity between the bodies is zero. After the moment t c , the elastic properties of the two bodies manifest and the effect is a removal till the bodies separate completely. The collision is considered finished at the moment t f , when the last points from the surfaces of the two bodies detach from contact.
The time between the instants t i and t c is named compression. The time duration between the moments t c and t f is called restitution.
To explain the Kane's paradox, the kinematical (Newton) definition of the coefficient of restitution must be disregarded and replaced by a new definition, (Poisson) that employs the dynamic concept of percussion, defined as the integral from the interaction force over the period of action of the force [21,22]. The new definition of the coefficient of restitution states that it equals the ratio between the normal percussions corresponding to the restitution and compression phases, respectively [23][24][25].
Wang et al. [26] made a classification of percussions, and concluded that tangential or frictional percussions might exist, which are not caused by normal percussions, as it happens in the case of forces. One of the main limitations of the model with rigid elements consists in the impossibility of force estimation during the collision process.
A second method of approach considers that during the time of impact the impact forces and the kinematic and dynamic parameters present a continuous and finite variation. It is accepted as an additional hypothesis that, as long as the collision takes place, the configuration of the system remains unaltered.
The evolution of the impact models was a natural one, beginning with the simplest models, consisting of two bodies of simple shape (sphere, cylinder, rod and plate) and between which simple motions exist (mono-axial translation and planar motion), or from small [27] to large deformations [28].
Timoshenko [29] starts from the relation force-deformation for a punctual Hertz contact, Ref. [30] particularized for the contact of two metallic spheres and obtains the expressions for the maximum normal approach, impact time and maximum force. It was experimentally revealed that the Timoshenko model, which accepts the same kinetic energy of the system, even after impact, corresponds to a reduced number of actual cases [31,32].
In chronological order, the second impact model is the impact with internal friction. The damped impact described by using a Kelvin-Voigt model [33][34][35] presents the advantage of hysteresis loop presence, but the weak points are the open hysteresis loop and the fact predicts attraction between the two bodies at the end of the collision, a statement that is different from physical reality.
Dubowski and Freudenstein [36,37] and, later, Hunt and Crossley [38], proved that, in order to obtain a model with the loop closed in origin, the differential equation of the model must be changed, namely the expressions of the coefficients from the Kelvin-Voigt model should be variable and proportional to a power of deformation.
To be remarked that in order to characterize the frictional losses the coefficient of restitution is used.
Next, from the evolution of the damped centric collision models it should be mentioned a reference work due to Lankarani and Nikravesh [39]. Here the authors consider that the entire variation of kinetic energy of the system is retrieved as damping work-the area of the hysteresis loop. Additionally, the hypothesis that the damping work is the same for the two phases of collision is also introduced. Under these conditions, Lankarani and Nikravesh [39] obtained a nonlinear differential equation that after integration leads to the expression of the approach between the balls. The authors mention that the model is applicable only for quasi-elastic collisions (c r > 0.9). Applications of multi-body models are suitable for different domains of engineering, from classical parts of machine elements [40,41], to granular [34,35] and powder materials [42,43] or biomechanics and robotics [44][45][46].
The Lankarani-Nikravesh equation was used nearly for twenty years as the basis for solving numerous multibody problems. Recently, Flores [47] resumed the model proposed by Lankarani and Nikravesh with the adjustment introduced by the hypothesis that in the phase plane, during the two phases of the impact, the characteristic point moves on a quarter of ellipse arch, hypothesis that allows for finding the coefficient of the term characteristic to damping. The Flores equation is similar to the Lankarani-Nikravesh equation with the difference that the numeric coefficient of the damping term is dissimilar. The Flores model eliminates the drawback of the previous model and ensures quasi-equality between the initial coefficient of restitution and the one generated by the proposed model.
A major application in engineering concerns mechanical systems with joint clearances, consequently these were studied by numerous researchers [48][49][50][51]. The most recent research studies concerning the viscoelastic impact phenomenon employ recent results from the theory of differential equations [52,53] validated by experimental results. There are also new theoretical models that simulate the behavior of materials via numerical methods [54], or, some that use for validation experimental data from the classical literature [42,43].
To summarize the above review, in order to apply the Hunt and Crossley equation for the impact model, three coefficients must be known (reduced mass, damping coefficient and contact stiffness), but the viscous damping coefficient it is not known. The actual method is compared to the models proposed by Lankarani and Flores, respectively, who adopt supplementary hypothesis for finding the damping coefficient, but their results do not validate the experimental results for small values of the coefficient of restitution.
Here there are deduced expressions for the next impact parameters: the maximum impact force, the maximum normal deformation and the deformation corresponding to the maximum force that allow for finding the precise values of these; the previous models, for the calculus of these very important parameters in practical applications, give values dependent on the chosen damping characteristic.
The novelty introduced by the present work consists in finding a first integral of the equation that permits the exact finding of damping characteristic and the method is valid for any value of coefficient of restitution, presenting importance for the compliant materials, in the domain of small coefficients of restitution. The proposed model aims to be an instrument for validation of experimental data from impact tests. Additionally, the method proposed, that is based on the viscoelastic model given by Hunt and Crossley, can be extended to the elasto-visco-plastic nonlinear impact model. A simplified algorithm for the calculus of impact parameters is also proposed, based on the remark that for a certain value of the exponent of deformation from the expression of impact force, all normalized contact parameters are invariants with respect to initial impact velocity and stiffness.

Materials and Methods
Two manners of approach were highlighted in the study of impact phenomenon. First of these considers an instantaneous impact process between rigid bodies while the second assumes deformable bodies and a finite period of time for impact occurrence.
The first methodology uses simple mathematical tools but unlike the second method, presents the major drawback that the interaction forces during impact cannot be estimated. For the second approaching manner, even for the simplest situations, the mathematical apparatus is much more difficult.
For the simplest impact model, as in Figure 1, the collision initiates when the first points from the two balls come into contact. The bodies are deformable and the regions from vicinity of the contact points distort; the normal approach between the two bodies is defined as the distance between the initial contact points from the un-deformed bodies. To summarize the above review, in order to apply the Hunt and Crossley equation for the impact model, three coefficients must be known (reduced mass, damping coefficient and contact stiffness), but the viscous damping coefficient it is not known. The actual method is compared to the models proposed by Lankarani and Flores, respectively, who adopt supplementary hypothesis for finding the damping coefficient, but their results do not validate the experimental results for small values of the coefficient of restitution.
Here there are deduced expressions for the next impact parameters: the maximum impact force, the maximum normal deformation and the deformation corresponding to the maximum force that allow for finding the precise values of these; the previous models, for the calculus of these very important parameters in practical applications, give values dependent on the chosen damping characteristic.
The novelty introduced by the present work consists in finding a first integral of the equation that permits the exact finding of damping characteristic and the method is valid for any value of coefficient of restitution, presenting importance for the compliant materials, in the domain of small coefficients of restitution. The proposed model aims to be an instrument for validation of experimental data from impact tests. Additionally, the method proposed, that is based on the viscoelastic model given by Hunt and Crossley, can be extended to the elasto-visco-plastic nonlinear impact model. A simplified algorithm for the calculus of impact parameters is also proposed, based on the remark that for a certain value of the exponent of deformation from the expression of impact force, all normalized contact parameters are invariants with respect to initial impact velocity and stiffness.

Materials and Methods
Two manners of approach were highlighted in the study of impact phenomenon. First of these considers an instantaneous impact process between rigid bodies while the second assumes deformable bodies and a finite period of time for impact occurrence.
The first methodology uses simple mathematical tools but unlike the second method, presents the major drawback that the interaction forces during impact cannot be estimated. For the second approaching manner, even for the simplest situations, the mathematical apparatus is much more difficult.
For the simplest impact model, as in Figure 1, the collision initiates when the first points from the two balls come into contact. The bodies are deformable and the regions from vicinity of the contact points distort; the normal approach between the two bodies is defined as the distance between the initial contact points from the un-deformed bodies. The normal approach x increases from the instant of contact beginning until it reaches a maximum value and afterwards it decreases; the collision ends when the bodies do not have common points any more. The period while the normal approach increases is called the compression phase, and the time elapsed from the moment of maximum approach till the end of impact is called the restitution phase. A parameter that has critical significance The normal approach x increases from the instant of contact beginning until it reaches a maximum value and afterwards it decreases; the collision ends when the bodies do not have common points any more. The period while the normal approach increases is called the compression phase, and the time elapsed from the moment of maximum approach till the end of impact is called the restitution phase. A parameter that has critical significance for impact process is the coefficient of restitution. It can be defined in kinematical (Newton) or dynamical (Poisson) manner.
In the case of centric collision, the two definitions are equivalent. For this situation, the coefficient of restitution (c r ) is defined according to Newton [8] as the ratio with changed sign between the relative velocities corresponding to the end and beginning of collision: and according to Poisson [9], as the ratio between the percussions from restitution phase and compression phase: The velocities corresponding to the initiation instant of collision were denoted by ( ) and the ones corresponding to the ending moment, by ("). The values of coefficient of restitution are situated between two extreme values: c r = 1 is characteristic for the perfect elastic collision, and c r = 0 is the value related to perfect plastic collision. For the perfect elastic collision between two smooth spheres, Timoshenko accepts that, during impact phenomenon, the variation of impact force is described by the relation of Hertz [38]: where the K coefficient, is the generalized stiffness parameter that considers the elastic characteristics, E, ν, of the two bodies and the geometry of contact regions, by r 1,2 : where: The exponent of deformation x is 3/2 for the cases when the impacting spheres are made of elastic materials [55,56].
Assuming the above hypotheses, for an imposed value of initial relative velocity, v r = v 0 , Timoshenko [29] finds the following: the maximum normal approach between the bodies x H , the value of maximum impact force, F H , and values for compression and restitution times, t c and t r , proving that the two values are equal and that the restitution and compression phases are symmetrical with respect to the moment of reaching maximum normal approach.
Concerning the kinetical energy variation of the system for the case of sub-unitary value of coefficient of restitution, c r , a quantitative characterization a an extremely difficult task. Goldsmith [1] shows that the kinetical energy variation is regained as heat generated by internal friction in the materials of the bodies, work of plastic deformation and energy of the waves traveling through the two bodies.
Allowing for the entire kinetical energy variation to be retrieved as heat generated by internal friction, Hunt and Crossley [38] show that it is inadequate to describe the impact by using a second-order linear homogenous differential equation (the internal friction force being proportional to the deformation rate and to the elastic force), since it leads to a hysteresis loop open in origin, and additionally foreseeing that, at the end of collision, the bodies attract each other instead of rejecting. They prove that, in order to get a hysteresis loop closing in origin when the elastic force is of Hertzian type, it is necessary to adopt a nonlinear viscoelastic model, which ensures zero values of the damping force for the initial and final moments [38], and for which the coefficient of deformation rate should be proportional to deformation. The equation proposed by Hunt and Crossley has the following form: where α = 3/2, Refs. [1,30,55,56], c is a constant (hysteresis damping factor) and m is the reduced mass of the system: The relations corresponding to the elastic impact of two spheres can be particularized for the impact between a sphere and a half-space, considering that the mass and radius of one of the spheres tend to infinity [1,56].
By denoting with γ = c/m; the characteristic equation of the impact process is as follows: In order to process the Equation (9), the value of the γ coefficient should be found. The theorem of kinetic energy in finite form is not practicable, since it cannot be obtained an analytical expression of the work done by internal friction forces.
To overcome this difficulty, Lankarani [39] accepts that the work performed by dissipative forces has the same value both for compression and for restitution phase, and, thus, the coefficient γ takes the following expression: The limits of the model are underlined even by Lankarani, who highlights that the model is applicable only for quasi-elastic collisions with c r > 0.9. For a Hunt-Crossley model, with the damping coefficient given by Relation (11), Ref. [39], a value for the initial velocity v 0 and a value for the coefficient of restitution, c r_in , accepted a priori, the equation is integrated, and then the final value of relative velocity is calculated, and the coefficient of restitution, c r_out , can be found. It is noticed that the two values of the coefficient can be considered quasi-equal only for quasi-elastic collisions, c r ≥ 0.9.
Recently, Flores retakes the problem of finding a solution for Equation (9). He starts from the observation that, in the case of damped oscillatory motion described by a linear equation, the characteristic point depicts in the phase space an arch from an ellipse, and accepts that for Equation (9) too, in the phase plane, both for compression and for restitution phase the characteristic point describes different quarters of ellipse. Based on this hypothesis, Flores [47] establishes the following expression for γ coefficient: Using Relation (12), Flores finds an analytical expression for maximum approach (x m ) and shows that, in the case of damped motion, the maximum impact force is not attained at the end of compression phase, but earlier. The model proposed by Flores satisfies much better the condition of coincidence between the final relative velocity generated by the model and the one established by using the definition of coefficient of restitution.
One considers again Equation (9) in the attempt of applying a numerical integration method, but two impediments are obstructing: the damping coefficient is not known precisely, and the extent of collision time is also unknown. Using Expression (12) given by Flores for damping coefficient, it was noticed that the Runge-Kutta method applied for Equation (9) is inappropriate when, for the impact period, it is specified a greater value than the actual one. Trying to integrate numerically, using the Runge-Kutta IV algorithm, the precise value of the impact time (t f ) must be specified, but it is not known in the Flores model. Smaller values for t f conduct to incomplete solution and greater values for t f conduct to blocking of the program. To surmount this impediment, the method was applied by starting with small values assumed for the collision period, and following that, the assumed impact time was increased successively till the moment when the method become inappropriate. It was observed that, for small values of coefficient of restitution, c r < 0.3, the final value obtained via integration of the equation is greater than the one obtained by using the definition of coefficient of restitution. In Figure 2, we present the plot of impact velocity variation with dimensionless time for the collision between a ball, r = 0.015, colliding at a velocity 2.8 m/s an immobile plane, for the coefficient of restitution c r = 0.2. ated by the model and the one established by using the definition of coefficient of restitution.
One considers again Equation (9) in the attempt of applying a numerical integration method, but two impediments are obstructing: the damping coefficient is not known precisely, and the extent of collision time is also unknown. Using Expression (12) given by Flores for damping coefficient, it was noticed that the Runge-Kutta method applied for Equation (9) is inappropriate when, for the impact period, it is specified a greater value than the actual one. Trying to integrate numerically, using the Runge-Kutta IV algorithm, the precise value of the impact time ( ) must be specified, but it is not known in the Flores model. Smaller values for conduct to incomplete solution and greater values for conduct to blocking of the program. To surmount this impediment, the method was applied by starting with small values assumed for the collision period, and following that, the assumed impact time was increased successively till the moment when the method become inappropriate. It was observed that, for small values of coefficient of restitution, < 0.3, the final value obtained via integration of the equation is greater than the one obtained by using the definition of coefficient of restitution. In Figure 2, we present the plot of impact velocity variation with dimensionless time for the collision between a ball, = 0.015, colliding at a velocity 2.8 m/s an immobile plane, for the coefficient of restitution = 0.2.

Parameters for the Un-Damped Impact
The present paper proposes a method for finding an exact solution of Lankarani-Nikravesh Equation (9). One can notice that, when there is no damping, Equation (9) describes the perfect elastic collision. For exemplification, we consider the next data: the masses of the balls = 10 kg, = 0.11 kg; the balls are made of the same material,

Parameters for the Un-Damped Impact
The present paper proposes a method for finding an exact solution of Lankarani-Nikravesh Equation (9). One can notice that, when there is no damping, Equation (9) describes the perfect elastic collision. For exemplification, we consider the next data: the masses of the balls m 1 = 10 9 kg, m 2 = 0.11 kg; the balls are made of the same material, steel, with ρ = 7800 kg/m 3 , E = 2.1 · 10 11 Pa, ν = 0.3; the exponent α = 3/2 is valid for any pair of elastic materials [1]; the radii of the balls are r 1 = 31.28 m, r 2 = 0.015 m, r 2 /r 1 = 4.71 · 10 −4 and r 1 >> r 2 ; the sphere 1 can be approximated in the vicinity of the contact region to an elastic half-space; the initial velocity v 0 = 2.8 m/s; and the coefficient of restitution c r = 0.4. For this situation, the maximum approach found by Timoshenko by applying the energy theorem for compression period is as follows: and it permits us to obtain the maximum impact force, F H : Therefore, the work for the compression phase L H is as follows: At a certain time (t) when the approach between the bodies is x, the velocity has the following expression: The duration of the compression and restitution phases is the same, denoted as t H , and the velocity given by Relation (16) is used in finding it: In the above expression, Γ(x) is the Eulerian integral of second kind. The values of the parameters x H , F H , L H and t H are especially useful because, first, they allow for an estimation of the magnitude of the characteristic parameters of damped impact, and second, they allow for expressing the results by normalized expressions. Moreover, even if in the expression of t H , the denominator is zero for the value of the upper limit of the integral, the integral is convergent, and it is expressed by using the improper integral, Γ(x).

Finding a First Integral of Nonlinear ODE
The method is based on the remark that we can find a first integral of Equation (9). By multiplying Equation (9) by .
x, it can be expressed under the following form: and then it is rewritten as follows: By taking into account that v = dx dt (20) the following equation with separable variables is obtained: A first integral of this equation is as follows: The initial conditions are imposed in order to find the constant of integration: Moreover, a particular solution of Equation (9) is obtained as follows: Based on this first integral, all the characteristic parameters of the viscoelastic impact can be obtained, and from here, the importance of this expression arises.
To be remarked that the first integral of Timoshenko's model cannot be obtained directly from Relation (16) by letting γ = 0 as a particular case.
Equation (17) represents the trajectory of the characteristic point from phase plane. Different from Reference [57], where the velocity was expressed as a function of normal approach, requiring numerical solution of the equation with respect to v, the present work emphasizes the fact that the equation may be analytically solved with respect to x, and therefore substantially reducing the computing time.
By introducing Expression (25) into Equation (9), we obtain the dependency between relative acceleration (a) and velocity (v).
Or in explicit manner:

Finding the Maximum Approach x m
At the instant the normal approach has reached the maximum value, it results that the derivative dx/dt = v should be zero. By denoting with x m the value of maximum approach, it can be found from Equation (17) and has the expression:

Obtaining the Coefficient of the Damping Force γ
It results: Equation (30) allows us to find the exact value of the damping coefficient. The equation is a transcendental one, and solving it requires a numerical procedure.
Denoting by f (γ), the function from the left side of Equation (30), one can easily notice that function is defined on the closed interval [−χ/v 0 , χ/(c r v 0 )]. The derivative of the function from the left side of the equation has a root: Moreover, in this point, the derivative of the function presents an extremum, specifically as follows: It is seen that f (0) = 0. At once, it is shown that, for 0 < c r < 1, f (γ ) > 0. At the ends of the domain, the function has vertical asymptotes, since we have the following: From the abovementioned, we draw the conclusion that Equation (30) has one root to be found in the interval 1−c r c r Figure 3 presents the graph of the function f (γ) obtained for the data considered above, and the normalization was made by using the value χ/v 0 = 6.11364 N·s m α+1 ·kg .
Mathematics 2021, 9, x FOR PEER REVIEW 11 of 32 In Figure 3, the solution of Equation (19), was found by using the bi-partition methodology.

Finding the Maximum Impact force
The presence of damping leads to a delay between the times when the maximum impact force and maximum approach are reached. In order to obtain the approach corresponding to the maximum force, Relation (5) is used, namely considering the proportionality between the impact force and acceleration. The right member of Equation (27) is differentiated with respect to velocity and from the imposed condition that the expression is zero, the velocity corresponding to maximum force result as solution of the transcendental equation: The solution is denoted as , and, according the physical meanings, it must be contained by the interval [− , ]. At the end points of the interval, the values of the function ℎ( ) are as follows: and Based Relation (24), the logarithm from Relation (36) can be solved, and it results in the following: The function ℎ( ) is continuous on the closed interval [− , ], and the values with opposite sign for the limits of the interval ensure the existence of at least one root of Equation (34). The derivative of the function ℎ( ) is as follows: In Figure 3, the solution of Equation (19), γ 0 was found by using the bi-partition methodology.

Finding the Maximum Impact force
The presence of damping leads to a delay between the times when the maximum impact force and maximum approach are reached. In order to obtain the approach corresponding to the maximum force, Relation (5) is used, namely considering the proportionality between the impact force and acceleration. The right member of Equation (27) is differentiated with respect to velocity and from the imposed condition that the expression is zero, the velocity corresponding to maximum force result as solution of the transcendental equation: The solution is denoted as v Fmx , and, according the physical meanings, it must be contained by the interval [−c r v 0 , v 0 ]. At the end points of the interval, the values of the function h(v) are as follows: and Based Relation (24), the logarithm from Relation (36) can be solved, and it results in the following: The function h(v) is continuous on the closed interval [−c r v 0 , v 0 ], and the values with opposite sign for the limits of the interval ensure the existence of at least one root of Equation (34). The derivative of the function h(v) is as follows: It has the following roots: and It was shown that γ 0 < χ/(c r v 0 ), and then it results in the following: To conclude, in the interval [−c r v 0 , v 0 ], the derivative h (v) has maximum one root v 2 . If the root v 2 is outside the interval, h(v) is monotonous with opposite sign at the ends of the interval (and, then, with a single root inside the interval). If the root v 2 is inside the interval, the function h(v) has an extremum inside the interval, and the root v Fmx is between v 2 and the endpoint for which the sign of the function is opposite to the sign of extremum of function. Figure It was shown that < /( ), and then it results in the following: To conclude, in the interval [− , ], the derivative ℎ′( ) has maximum one root . If the root is outside the interval, ℎ( ) is monotonous with opposite sign at the ends of the interval (and, then, with a single root inside the interval). If the root is inside the interval, the function ℎ( ) has an extremum inside the interval, and the root is between and the endpoint for which the sign of the function is opposite to the sign of extremum of function. Figure 4 presents the plot of the function ℎ( ) on the interval [− , ] and the found root. The approach corresponding to this velocity is calculated with the following relation: The maximum force is obtained with the next relation:

Phase Plane and Hysteresis Loop
With the relations of normal approach and acceleration expressed as functions of velocity, the phase space ( Figure 5) and hysteresis loop ( Figure 6) were plotted.  The approach x Fmx corresponding to this velocity is calculated with the following relation: The maximum force is obtained with the next relation:

Phase Plane and Hysteresis Loop
With the relations of normal approach and acceleration expressed as functions of velocity, the phase space ( Figure 5) and hysteresis loop ( Figure 6) were plotted. Mathematics 2021, 9, x FOR PEER REVIEW 13 of 32 In Figure 5, there were also traces of the quarters of ellipses passing through the point that separates the phases and through the points corresponding to impact starting and impact ending. In Figure 6, alongside the hysteresis loop, we traced the loading-unloading curve for un-damped impact.

Relations for Finding the Characteristic Periods of Impact
To illustrate the variation with time of the parameters characteristic to collision, the definition relation of the velocity is considered: It results in the following: From the definition of acceleration, we get the following: = It results in the following: In Figure 5, there were also traces of the quarters of ellipses passing through the point that separates the phases and through the points corresponding to impact starting and impact ending. In Figure 6, alongside the hysteresis loop, we traced the loading-unloading curve for un-damped impact.

Relations for Finding the Characteristic Periods of Impact
To illustrate the variation with time of the parameters characteristic to collision, the definition relation of the velocity is considered: It results in the following: From the definition of acceleration, we get the following: = It results in the following: In Figure 5, there were also traces of the quarters of ellipses passing through the point that separates the phases and through the points corresponding to impact starting and impact ending. In Figure 6, alongside the hysteresis loop, we traced the loading-unloading curve for un-damped impact.

Impact Parameters as Time
It results in the following: From the definition of acceleration, we get the following: It results in the following: Relations (44)- (47) were written based on the relation x = x(v) that permits the explicit form of deformation versus velocity. For the general case, the velocity is expressed as follows: The following equation must be added to Relations (44) and (46): from which it results: Notice that, during the entire impact process, the velocity has a continuous and monotonous variation from v 0 to −c r v 0 , while the approach increases from 0 to x m during the compression period, and after that, it decreases to zero during the restitution period, and thus, for a given value of deformation, there are two possible values of velocity, v + for compression phase and v − for restitution phase. v Knowing the variation of normal approach as function of velocity, x = x(v) (relation 25) and respective a = a(v) (relation 27) and considering t = 0 when v = v 0 as starting time, the instant at which the velocity v end has an imposed value is given by the following: and: where dx(v)/dv is the derivative of Relation (25) with respect to velocity and we have the following: where v = v(x), respectively, and x end is the deformation corresponding to the instant when the velocity is v end . In order to compute the time of reaching the maximum force, t Fmx , the compression time (t c ) and the total time (t f ), the following relations are used in Equations (52) and (53), where v end takes the values v end = v Fmax , v end = 0 and v end = −c r v 0 , respectively. It results in the following: and The same parameters can be obtained by using Relation (54), where x end must be replaced by x Fmx , x c and 0.

Elimination of Singularities from the Expressions of Impact Periods
The complicate form of the integrands in Relations (55) and (56) cannot allow for closed forms for primitives, and the numerical methods are therefore demanded. Analyzing the integrands from Equations (55) and (56), we see that the integrands tend to infinity for v = v 0 and for v = −c r v 0 when the entire duration of collision is considered, as presented in Figures 7 and 8. ( ) ( ) ( )

Elimination of Singularities from the Expressions of Impac
The complicate form of the integrands in Relations (55) an closed forms for primitives, and the numerical methods are there ing the integrands from Equations (55) and (56), we see that the in for = and for = − when the entire duration of collis sented in Figures 7 and 8.
The same parameters can be obtained by using Relation (54), where must be replaced by , and 0.

Elimination of Singularities from the Expressions of Impact Periods
The complicate form of the integrands in Relations (55) and (56) cannot allow for closed forms for primitives, and the numerical methods are therefore demanded. Analyzing the integrands from Equations (55) and (56), we see that the integrands tend to infinity for = and for = − when the entire duration of collision is considered, as presented in Figures 7 and 8.   If the integration is performed with respect to deformation, the integrals from Relation (57) become divergent when the domain contains the time of maximum approach (corresponding to v = 0). A similar aspect is met for the elastic collision when calculating the approach time with Relation (17), with the integral being expressed with the gamma function Γ(x) and obtaining a finite approach time, t H , which is in agreement with the physical meaning. In order to obtain the value of time t end corresponding to x end and v end (for restitution phase), two values of velocity are considered in the interval of velocity variation [v 0 , v end ], denoted v c > 0 for compression and v r < 0 for restitution, as in Figure 9.
If the integration is performed with respect to deformation, the integrals from Relation (57) become divergent when the domain contains the time of maximum approach (corresponding to = 0). A similar aspect is met for the elastic collision when calculating the approach time with Relation (17), with the integral being expressed with the gamma function ( ) and obtaining a finite approach time, , which is in agreement with the physical meaning. In order to obtain the value of time corresponding to and (for restitution phase), two values of velocity are considered in the interval of velocity variation [ , ], denoted > 0 for compression and < 0 for restitution, as in The time is obtained with the following relation: All integrals from Expression (59) are convergent. The difficulty is in the fact that there are no closed forms for the dependencies = ( ) and = ( ), these being obtained by numerical solving of Equation (24), with appreciable calculus volume. Aiming to obtain the values of the characteristic impact parameters for different times, it is considered that the explicit dependency of approach versus velocity = ( ) is available by the first integral, Relation (25). Based on this fact, the interval [ , ] is divided into equal subintervals by + 1 points. A vector ( ) is generated with the following elements: Another vector ( ) is generated, with the following elements: Two values are chosen for the index, for the point ( ; ) from compression phase and for the point ( ; ) from restitution phase. The integrals from Relation (59) Figure 9. The subintervals from the phase portrait ensuring the convergence of the integrals from the expression of time.
The time t end is obtained with the following relation: All integrals from Expression (59) are convergent. The difficulty is in the fact that there are no closed forms for the dependencies v = v + (x) and v = v − (x), these being obtained by numerical solving of Equation (24), with appreciable calculus volume. Aiming to obtain the values of the characteristic impact parameters for different times, it is considered that the explicit dependency of approach versus velocity x = x(v) is available by the first integral, Relation (25). Based on this fact, the interval [v 0 , v end ] is divided into n p equal subintervals by n p + 1 points. A vector (v) is generated with the following elements: Another vector (x) is generated, with the following elements: Two values are chosen for the index, n + for the point (v c ;x c ) from compression phase and n − for the point (v r ;x r ) from restitution phase. The integrals from Relation (59) are replaced by the attached Darboux sums, with the value of the function calculated for the middle of the subinterval, and it is obtained as follows: Based on the relation above, a recurrence formula can be obtained for the moments (t k ), and the vector (t) is obtained, the vector containing the moments for which the approach and velocity are x k and v k .
The last value from the vector (t) corresponds to the moment t end . The Figure 10 presents the dependency of the values t Fmx , t C and t f found by the above algorithm, for a set of increasing values of the number of points, n p . It can be observed that the values of t Fmx , t C and t f stabilize after a relatively small number of points, n p = 100.
Based on the relation above, a recurrence formula can be obtained for the moments ( ), and the vector ( ) is obtained, the vector containing the moments for which the approach and velocity are and .
The last value from the vector (t) corresponds to the moment . The Figure 10 presents the dependency of the values , and found by the above algorithm, for a set of increasing values of the number of points, . It can be observed that the values of , and stabilize after a relatively small number of points, = 100. For testing the precision of the method, the results of the times and calculated with the above-described methodology are presented next. As shown in Johnson [30] and Popov [56], the Hertzian point contact between two spheres can be considered as an equivalent contact between a sphere and a half-space. Goldsmith [1] also considers the spherehalf-space contact as a particular situation of sphere-sphere contact. Therefore, here it is considered the quasi-elastic impact between a smooth steel ball, radius = 15 mm that hits with velocity = 2.8 m/s an immobile smooth steel prismatic block of infinite mass.
The coefficient of restitution is adopted, = 0.999. This value was chosen, since, for the perfect elastic impact of two spheres, the relation for total impact time given by Timoshenko is as follows:  For testing the precision of the method, the results of the times t C and t f calculated with the above-described methodology are presented next. As shown in Johnson [30] and Popov [56], the Hertzian point contact between two spheres can be considered as an equivalent contact between a sphere and a half-space. Goldsmith [1] also considers the sphere-half-space contact as a particular situation of sphere-sphere contact. Therefore, here it is considered the quasi-elastic impact between a smooth steel ball, radius r = 15 mm that hits with velocity v 0 = 2.8 m/s an immobile smooth steel prismatic block of infinite mass. The coefficient of restitution is adopted, c r = 0.999. This value was chosen, since, for the perfect elastic impact of two spheres, the relation for total impact time given by Timoshenko is as follows: The equation of motion was integrated by using γ 0 and t f obtained by the presented methodology for two values of the number of discretization points. The results for the approach versus time obtained both by the proposed method and by the Runge-Kutta method are presented comparatively in Figure 11a for n p = 20 and in Figure 11b for n p = 200. It can be observed that, for the Runge-Kutta algorithm, the points are equidistant with respect to time, but for the present method, the points are congested in the vicinity of the time corresponding to maximum approach, which is also the interval of maximum attention for engineering applications.
In Figure 12, we present the variation of compression and restitution times as functions of index of points. the approach versus time obtained both by the proposed method and by the Runge-Kutta method are presented comparatively in Figure 11a for = 20 and in Figure 11b for = 200. It can be observed that, for the Runge-Kutta algorithm, the points are equidistant with respect to time, but for the present method, the points are congested in the vicinity of the time corresponding to maximum approach, which is also the interval of maximum attention for engineering applications. In Figure 12, we present the variation of compression and restitution times as functions of index of points.

Graphical Representation of Impact Parameters with Time
At this stage, for a set of equidistant values of velocity, the vectors and , corresponding to deformations and times, respectively, are available. The time variations for approach and velocity, for = 0.4, are presented in Figures 13 and 14, respectively.

Graphical Representation of Impact Parameters with Time
At this stage, for a set of equidistant n p values of velocity, the vectors x k and t k , corresponding to deformations and times, respectively, are available. The time variations for approach and velocity, for c r = 0.4, are presented in Figures 13 and 14, respectively. One can observe that the presence of damping has as consequence a shorter duration of compression compared to restitution. Concerning the velocity, it is observed that, by the present method, the final velocity is precisely −c r v 0 . Figure 12. Dependence of times on the order number of the points.

Graphical Representation of Impact Parameters with Time
At this stage, for a set of equidistant values of velocity, the vecto corresponding to deformations and times, respectively, are available. The t for approach and velocity, for = 0.4, are presented in Figures 13 and 14 One can observe that the presence of damping has as consequence a short compression compared to restitution. Concerning the velocity, it is observe present method, the final velocity is precisely − . The time variations of forces are represented in Figure 15: the Hertzian (elastic) ; the damping force, ; and the total impact force, . The work produced by forces is represented in Figure 16. Moreover, during the restitution phase, the dam force is negative, and that is an attraction force.  The time variations of forces are represented in Figure 15: the Hertzian (elastic) force, F Hz ; the damping force, F a ; and the total impact force, F. The work produced by these forces is represented in Figure 16. Moreover, during the restitution phase, the damping force is negative, and that is an attraction force. The time variations of forces are represented in Figure 15: the Hertzian (elastic) forc ; the damping force, ; and the total impact force, . The work produced by thes forces is represented in Figure 16. Moreover, during the restitution phase, the dampin force is negative, and that is an attraction force.

Characteristic Curves
By applying the method detailed above, we estimated the total impact time for the models of Lankarani ( ) and Flores ( ) were estimated, and after that, using the values for the coefficients , and in Runge-Kutta methodology, the differential equations were integrated for the proposed model and for the Lankarani and Flores models.   By applying the method detailed above, we estimated the total impact time for the models of Lankarani (t L f ) and Flores (t F f ) were estimated, and after that, using the values for the coefficients γ L , γ F and γ 0 in Runge-Kutta methodology, the differential equations were integrated for the proposed model and for the Lankarani and Flores models. For two values of the coefficient of restitution, c r = 0.7 and c r = 0.2, the time dependence of approach, velocity and impact force are presented comparatively in Figures 17-19. For the same values of coefficient of restitution, the characteristic curves in the phase plane and the hysteresis loops are plotted in Figures 20 and 21, respectively. models of Lankarani ( ) and Flores ( ) were estimated, and after that, using the values for the coefficients , and in Runge-Kutta methodology, the differential equations were integrated for the proposed model and for the Lankarani and Flores models. For two values of the coefficient of restitution, = 0.7 and = 0.2, the time dependence of approach, velocity and impact force are presented comparatively in Figures 17-19. For the same values of coefficient of restitution, the characteristic curves in the phase plane and the hysteresis loops are plotted in Figures 20 and 21, respectively.     From the presented figures, one can conclude that, for values of coefficient of restitution, c r ≥ 0.7, all three models give similar results, while for smaller values, c r < 0.7, the current model and the Flores model present close results of the characteristic parameters for impact phenomenon, but the Lankarani model leads to significantly different results.

The Effect of Coefficient of Restitution
In order to describe quantitatively the effect of coefficient of restitution upon the solution of the differential equation of impact, the presented methodology was applied for studying the dependency of the coefficient (γ), maximum approach (x m ), maximum impact force (F max ) and total impact time (t f ). These variations are presented in Figures 22-25       Proposed model Flores model [47] Lankarani model [39] Coefficient of restitution, cr Dimensionless impact time tf/ tH

A Simplified Algorithm for Finding the Impact Parameters
As noticed from above, the proposed methodology requires an appreciabl of numerical calculus for solving the transcendental equations and for evaluatio inite integrals. This drawback is overcome after analyzing the effect of initial im locity ( ) and contact stiffness ( ) upon the characteristic parameters of impact enon, when the following remark arises: all normalized contact parameters of imp variants with respect to velocity ( ) and stiffness ( ). From here, an important prac sequence is that, for a given value of coefficient of restitution, it is sufficient to Lankarani model [39] Coefficient of restitution, cr Dimensionless maximum impact force Fmax/FH

A Simplified Algorithm for Finding the Impact Parameters
As noticed from above, the proposed methodology requires an appreciable volume of numerical calculus for solving the transcendental equations and for evaluations of definite integrals. This drawback is overcome after analyzing the effect of initial impact velocity (v 0 ) and contact stiffness (χ) upon the characteristic parameters of impact phenomenon, when the following remark arises: all normalized contact parameters of impact are invariants with respect to velocity (v 0 ) and stiffness (χ). From here, an important practical consequence is that, for a given value of coefficient of restitution, it is sufficient to solve the equation of the model for a velocity (v 0 ) and a contact stiffness (χ), and after that, for any other values of the two parameters, v 0 and χ, all the characteristic parameters of impact can be found by using a much simpler methodology.
For this purpose, for the values of coefficient of restitution are as follows: We evaluated the values of the main normalized parameters characteristic to impact: γ 0 , viscous damping coefficient; t C , compression time; t f , total impact time; v Fmx , velocity corresponding to the maximum impact force; x Fmx , approach corresponding to the maximum impact force; x m , maximum approach; F max , maximum impact force; F emax , elastic force corresponding to maximum approach, x m ; F amax , maximum damping force for compression phase (rejection force); F amin , maximum damping force for restitution phase (attraction force); L C , work done during the compression phase; L r , work done during the restitution phase; L e , work of elastic force during the compression phase; L t , total work.
The calculated values are presented in Table 1.  Based on the results presented in Table 1, the variations of dimensionless par versus coefficient of restitution are plotted in Figures 26-31.       For a given impact, there must be stipulated the radii of the balls, , ; the masses, , ; the elastic parameters, , and , ; the exponent, = 3/2 for Hertzian contact; the initial impact velocity, ; and the value of the coefficient of restitution.
Then, the next steps must be followed: 1. The elastic constant of the contact ( ) and the contact stiffness ( ) are determined; 2. The parameters characteristic to Hertzian un-damped impact, , , and , are For a given impact, there must be stipulated the radii of the balls, r 1,2 ; the masses, m 1,2 ; the elastic parameters, E 1,2 and ν 1,2 ; the exponent, α = 3/2 for Hertzian contact; the initial impact velocity, v 0 ; and the value of the coefficient of restitution.
Then, the next steps must be followed: 1.
The elastic constant of the contact (K) and the contact stiffness (χ) are determined; 2.
The parameters characteristic to Hertzian un-damped impact, t H , F H ,x H and L H , are found; 3.
The normalized parameters are found by interpolation; 4.
The actual values of the parameters are found by multiplication of the normalized parameters by the corresponding constants.

Conclusions
The main conclusions of the work are as follows: The values of viscous damping coefficient obtained by the actual method were used for integration of the equation of motion for the whole range of coefficients of restitution (0 < c r < 1). For the same values of coefficients of restitution, the equations of the models Lankarani and Flores were numerically integrated (Runge-Kutta method) and the dimensionless parameters of impact were compared in a graphical manner. The discrepancies between the proposed model and the previous ones are obvious in the range of 0 < c r < 0.4. We emphasize that the present model is applicable for any coefficient of restitution, being especially useful for soft materials. - The calculus of the impact parameters can be substantially simplified based on the observation that for a fixed value of α, the exponent of deformation, all dimensionless impact parameters are invariants with respect to initial velocity and stiffness. A table of dimensionless impact parameters is given for rapid interpolation, using the known c r and the Hertzian impact parameters.
The authors identify, for future research, the following directions: -The model can be extended to hyperelastic materials having as force-deformation law a most general expression. - The present model can also be employed for modeling the dashpots with large velocities (>50 m/s) when the damping force is proportional to the velocity squared. - The model can be extended to model the collision of elasto-visco-plastic nonlinear bodies-which are frequently encountered in concrete dynamical systems. This assertion is based on Johnson's remark [30] that, only for low-impact velocities, the collision can be regarded as an elastic phenomenon. As an example, Johnson shows that, for steel, the impact velocity must satisfy the condition v < 7.5 m/s This is an aimed future work that can be applied for any domain, from classical industrial dynamical systems to medical devices, food-processing industry, pharmaceutic industry and powder technology. - The combined forces (normal and tangential) frequently met in practical applications are an important subject, and, therefore, the issue of tangential impact forces is a task for future research, both for theoretical modeling and for experimental work. The friction forces, which are tangential forces, are a consequence of normal forces, and we meet them in modern but usual applications, for example, the final element of robots-the grasper, the medical prosthesis and the landing gear from aircrafts. Funding: This research received no external funding.