Optimal damping coefficient for a class of continuous contact models

In this study, we develop an analytical formula to approximate the damping coefficient as a function of the coefficient of restitution for a class of continuous contact models. The contact force is generated by a logical point-to-point force element consisting of a linear damper connected in parallel to a spring with Hertz force–penetration characteristic, while the exponent of deformation of the Hertz spring can vary between one and two. In this nonlinear model, it is assumed that the bodies start to separate when the contact force becomes zero. After separation, either the restitution continues or a permanent penetration is achieved. Therefore, this model is capable of addressing a wide range of impact problems. Herein, we apply an optimization strategy on the solution of the equations governing the dynamics of the penetration, ensuring that the desired restitution is reproduced at the time of separation. Furthermore, based on the results of the optimization process along with analytical investigations, the resulting optimal damping coefficient is analytically expressed at the time of impact in terms of system properties such as the effective mass, penetration velocity just before the impact, coefficient of restitution, and the characteristics of the Hertz spring model.

. Different methods can be used to model the impact in a multibody system. These approaches have been reviewed in detail in [7,19].
Piecewise approach has largely been used to model the impact. This technique is mainly applicable when the duration of contact is significantly short. Based on this assumption, the configuration of the system including translational and angular coordinates are considered to remain unchanged during the period of contact. However, the system experiences a discontinuous change at the velocity level. In the piecewise method, the integration of the equations of motion is halted just before the impact occurs. Then impulse-momentum equations are formed. These equations relate the change in both linear and angular momenta of the system to impulsive forces and moments. Since applied loads (forces and moments) are bounded, their impulses during the very short period of contact are ignored in the momentum balance equations. These equations along with the elastic characteristic of the impact, which is reflected in the coefficient of restitution, are then used to compute the jumps in the velocities of the system [22,23]. These results in addition to the velocities just before the impact can now be used to find the velocities immediately after the impact. The numerical integrator uses the resulting velocities and the coordinates of the system just before the impact as a set of initial conditions to integrate the equations of motion of the system.
Another approach treats an impact as a continuous phenomenon. Allowing penetration between contacting bodies during the period of contact, this model results in a continuous change in the positions and velocities. Additionally, without interrupting the integration of equations of motion, contact forces are activated at the contact points and included in the dynamic equations of motion as long as the bodies remain in contact. In the Kelvin-Voigt viscoelastic model, the contact force is generated by a logical point-to-point linear spring-damper element [8]. However, this linear model is incapable of capturing the nonlinear nature of the impact. Furthermore, it generates a nonphysical tensile contact force if zero penetration is considered as the separation condition. To remove this shortcoming, it is assumed in [2] that the bodies separate at the moment when the contact force becomes zero, while allowing the deformation to continue after the separation of the bodies until the end of the complete restitution. This model has been further extended by using a fractional order damping element [3]. The Hertz model [10] can better capture nonlinear characteristics of the impact. In this continuous model, the force element consists of a spring with nonlinear force-penetration relationship which is obtained based on the local stress distribution. This method can only be used for impacts with elastic characteristics since it does not consider any energy dissipation. Adding a logical linear damper to the Hertz spring can further capture the energy dissipation during impact. However, the model generates a nonphysical tensile force before the complete restitution is achieved if the separation of the bodies is considered at the instant at which the penetration becomes zero [12]. To overcome this artifact, Hunt and Crossley introduced a hysteresis damping factor to modify the damping force [12]. In fact, the modification of the damping force in this model results in zero contact force and penetration at the beginning and the end of the period of contact. Although a multibody system experiences an impulsive behavior, because of the extremely short contact period, the dissipation of damping forces in the system except for the one related to the local deformation are mainly neglected. Therefore, due to the importance of local deformations at the contact point, different methods have been used to characterize the hysteresis damping factor based on the properties of the direct central impact of two spherical objects [6,11,15,17,27,29]. Similar estimation has been presented in [16,29] when permanent indentation after the separation is included in the model. Other estimations for the hysteresis damping factor have been developed when the force model changes during compression and restitution phases of contact [13,20]. These estimations have been expressed in terms of the coefficient of restitution and of the geometrical and mechanical properties of the contacting bodies, which are reflected in the parameters of the Hertz contact model. Since all of these methods use approximations, in some situations the continuous model with the estimated damping factor may not necessarily recover the desired restitution based on which the estimation has been generated.
The damping coefficient has been characterized as a function of the coefficient of restitution with a model containing a linear damper connected in parallel to a spring with Hertz force-deformation property of two interacting spheres [25], where the exponent of deformation is known to be n = 3/2. In this paper, we extend the characterization of the damping coefficient for a general class of Hertz spring models in which the exponent of deformation in the force-penetration relation can vary between one and two. Furthermore, it is assumed that the separation occurs when the contact force becomes zero, while either the complete restitution may have not yet been achieved or the permanent penetration may have been reached. To determine the damping coefficient, we apply an optimization process to ensure that the desired restitution is reproduced at the time of separation. In fact, this optimization is conducted on the solution of the equations governing the dynamics of the direct central impact. We further conduct a series of analytical studies and numerical experiments on the characteristics of the resulting optimal damping coefficients. It is shown that the optimal damping coefficient can analytically be expressed in terms of system properties at the time of impact such as the effective mass, penetration velocity just before the impact, coefficient of restitution, and the geometrical and mechanical properties of the contacting bodies which are reflected in the parameters of the Hertz contact model. Simulation results show that the contact force model with the proposed analytical approximation of the damping coefficient performs perfectly in recovering the desired restitution, while guarantying non-zero contact force at the time of separation.
In this paper, we first provide an overview of the continuous analysis of frictionless direct central impact in Sect. 2. This is followed by implementing the optimization-based approach to compute the damping coefficient for the zero contact force separation condition in Sect. 3. Characterization of the optimal damping coefficient through analytical investigations and numerical experiments are then conducted in Sect. 4. This is followed by numerical verification and simulation in Sects. 5 and 6, respectively. Finally, conclusions are drawn in Sect. 7.

Analysis of direct central impacts through continuous force models
In many applications of multibody dynamics, the impact is assumed to generate a local deformation. Therefore, the most common continuous contact models are constructed based on the characteristics of the direct central impact between two objects, mainly in spherical shapes [1,6,11,12,14,16]. Figure 1 shows the schematic of two colliding spheres of masses m i and m j . The unit vector n is the direction normal to the surface of contact. In the direct central impact it is assumed that the velocities of the colliding bodies just before the impact at − t , i.e. − V i and − V j , lie on the axis of vector n. Additionally, similar to the piecewise modeling of the direct central impact, in the continuous contact model, the effect of the applied loads is considered insignificant compared to the contact force due to the very short contact duration. Therefore, as shown in Fig. 1, during the contact period only contact force f c is continuously being applied to the spheres along the n axis. In this scheme, the contact cycle is divided into two phases: compression and restitution. The compression phase starts at − t when contact between the bodies occurs. The indentation (penetration) then increases until it reaches its maximum value, at which both bodies find the same velocity. At this moment, the restitution phase starts which is followed by a continuous decrease in the penetration until both bodies separate at + t with velocities + V i and + V j , respectively.
Following the free-body diagrams of Fig. 1, we can construct the equations of motion of each body during the period of contact along the n axis as where δ i and δ j denote the coordinates of the mass centers of spheres i and j along the n axis, respectively. Defining the penetration δ and effective mass m as we develop the following Ordinary Differential Equation (ODE) governing the dynamics of the penetration: The main challenge in the continuous contact model is to determine a mathematical expression for f c and tune its parameters based on known information about the system during the impact in order to better capture the physics of this phenomenon. Hertz model uses stress distribution and presents the contact force as the following nonlinear spring [10]: In this model, both k and n depend on the geometry and/or the material properties of the contacting bodies. For instance, for two spheres in contact or for a locally spherical surface and a plane n becomes 3/2, while in the contact between two locally plane surfaces n is established as one [15]. Since the Hertz model does not introduce any energy dissipation, it may not accurately present non-elastic impacts. To consider the energy dissipation during an impact, a linear damping force cδ is included, resulting in This yields the following ODE governing the dynamics of the penetration with the same initial conditions of Eq. (4): The solid curve in Fig. 2 shows the behavior of this contact force versus penetration. This figure shows a non-zero compressive contact force at point A which corresponds to − t , the start of the compression phase. The compression continues to point B where the maximum deformation δ m is reached. The restitution phase then starts at B and continues to point C where the spheres separate. The dotted curve in this figure indicates that using this model beyond point C generates a nonphysical tensile force. Therefore, in order to avoid any contact force at the time of separation, the bodies need to be released at point C which corresponds to + t . Due to computational limitations, as of several decades ago, it was not possible to predict the exact time of separation of the two bodies in forward dynamic integration of the equations of motion. Therefore, in order to overcome the artifact of tensile force of Eq. (8), Hunt and Crossley [12] proposed the damping force to be presented as μδδ n , where μ is called the "hysteresis damping factor". This description for damping transforms the continuous contact force to This results in the following ODE: with the same initial conditions of Eq. (4). The dashed curve in Fig. 2 shows the behavior of this contact force versus penetration. As it is observed, colliding bodies experience zero contact force and penetration at the beginning and at the end of contact. As a result, the model does not produce a tensile force. The hysteresis damping factor μ in Eq. (10) is a parameter representing the dissipation of energy in an impact. Therefore, researchers have attempted to find a relationship between μ and the corresponding coefficient of restitution since the introduction of this contact force model. The coefficient of restitution between the two colliding spheres of Fig. 1 is defined as The proposed relationships between μ and e have been derived using the following energy balance during the entire period of the direct central impact: The left hand side of this equation represents the drop in the kinetic energy of the system during the impact, and the integral on the right side is the work done by the contact force during the entire cycle of impact. Hunt and Crossley suggested the following expression approximating μ as a function of e [12]: Lankarani and Nikravesh [15,17] used the general trend of Hertz contact model and expressed an estimation of the damping factor as The estimations presented in Eqs. (13) and (14) can be used if the coefficient of restitution is close to unity [6,15,17]. This drawback was removed with the following estimation: proposed by Hu and Gua [11]. In another attempt, using fundamental characteristics of the Kelvin-Voigt model, Ye et al. [29] and Flores et al. [6] overcame the shortcoming of Eq. (14) by estimating the damping factor as The estimations presented in Eqs. (15) and (16) are applicable to a wider range of coefficients of restitution [6,11,29]. It should be noted that these expressions have been proposed based on different approximations of the solution of Eq. (10). As such, using the estimated μ of Eqs. (13)-(16) and solving Eq. (10) for the entire period of contact, the output coefficient of restitution may not necessarily be the same as the desired e based on which μ has been estimated. In other words, the estimated μ may not accurately reproduce the desired restitution.
As reported in [9], the exact analytical solution of ODE of Eq. (10) has been presented in [21,28]. Therefore, without any approximation, the damping factor based on which the desired restitution can be recovered is computed using the following expression [9]: where d ∈ [0, 1] is the solution of the following equation:

Hertz contact model with linear damping
The solid curve in Fig. 2 shows typical behavior of the Hertz contact model with linear damping, as expressed by Eq. (8). This curve clearly shows that, if the contact force is not removed from the model at point C when the two bodies separate, an undesirable tensile force will continue to be present in the equations of motion. Due to the computational deficiencies that existed several decades ago, the exact times of contact and separation such as points A and C in Fig. 2, respectively, could not be accurately detected during the integration of the equations of motion. In order to avoid such computational shortcomings, researchers adopted the contact model of Hunt and Crossley. However, as the dashed curve in Fig. 2 shows, this model artificially changes the characteristics of the original Hertz contact model with linear damping. With advances in numerical methods and computational capabilities, most well developed numerical integrators are now capable of detecting the occurrence of an event, such as the start or the end of a contact. In order to take advantage of the existing numerical capabilities, in the study that is presented in this paper, we return to the Hertz contact model with linear damping, as expressed by Eq. (8). Then we attempt to develop an analytical formula approximating the damping coefficient c as a function of a desired value for the coefficient of restitution.
In this model, to avoid any tensile contact force during the period of impact, we assume that in the restitution phase the bodies start to separate at t = s t provided that the contact force is zero, i.e.
As shown in Fig. 2, in this model point C corresponds to t = + t = s t at which the separation occurs. In this study, s δ could represent either a permanent penetration with no restitution after the separation or a temporary indentation which allows the continuation of the restitution after the bodies separate.
In order to find the damping coefficient c with the main objective of recovering the desired restitution e, the solution of the ODE presented in Eq.
Without using any approximation for the penetration velocity profile, in this work, we apply the following optimization process to obtain the value of c: This optimization in fact ensures that the computed damping coefficient can accurately recover the desired restitution. We note that, for given e, m, −δ , k, and n, the values of s δ and sδ at the time of separation + t = s t are dependent on c, which is unknown in advance.
To numerically solve the proposed optimization problem and compute the damping coefficient, any optimization technique could potentially be used. For instance, this problem has been addressed in [25] by using the method of bisection. Since most optimization methods use an iterative procedure, the optimal damping coefficient c opt is obtained when the following condition is satisfied: where is a small number and e out is defined in Eq. (22). Furthermore, in order to halt the integration process of the equations of motion at each iteration, we can use MATLAB built-in function "event" to detect the instant at which the separation condition of Eq. (20) is met. It should be noted that the result of the optimization problem automatically satisfies the energy balance of Eq. (12) since we have not used any approximation.

Analytical expression for optimal damping coefficient
We now conduct an analytical and a numerical investigation to characterize c opt as a function of system parameters m, −δ , k, n, and the desired e.

Analytical investigation to characterize c opt
In this section, we analytically study some characteristics of c opt . We start with the model of Eq. (8) in the absence of damping force which reduces to the original Hertz model with the following analytical solution: Settingδ = 0 yields the maximum penetration of the Hertz model as We note that the model of Eq. (8) contains a damping force. The dissipation due to the damping force does not allow the penetration at the separation instant to exceed the maximum indentation obtained in Eq. (27). This can mathematically be expressed as On the other hand, the damping coefficient can be obtained using the condition of Eq. (20): Comparing Eqs. (28) and (29) yields Imposing the desired restitution condition sδ = −e ( −δ ) on this inequality, we are allowed to replace c with c opt , resulting in Using the maximum penetration of Eq. (27) in the absence of damping, we can express Eq. (31) as This inequality can be further simplified as Introducing a slack variable λ ∈ [0, 1], this inequality can be transformed to the following equality: The natural logarithm of this equation can be expressed as ln c opt = ln λ − ln e + n n + 1 ln n + 1 2 The inequality of Eq. (32) can alternatively be developed by rescaling penetration and time byδ = δ δm andt = t δm/ −δ , respectively, and using them in the separation condition of Eq. (20) and enforcing sδ = −e ( −δ ) andδ ≤ 1.

Computational experiments to characterize c opt
Based on the analytical study of c opt in the previous section, we realize that the optimal damping can be expressed using Eq. (34) or (35). The term λ in these equations is not just a parameter; i.e., it could potentially be a function of system properties. In the following sections, we perform four simulation experiments to investigate the dependency of λ on m, −δ , k, n, and the desired e.
We should note that, for each of the following investigations, the value of in Eq. (24) is selected as 10 −8 . Furthermore, the curve fittings in these studies are performed using MATLAB curve-fitting toolbox.

Dependency of λ on m
In the first investigation, we numerically study the dependency of λ on m. In order to do so, we assume that −δ = 1 m/s, k = 1 N/m n , e = 0.7. Then we sweep the effective mass and n by choosing, respectively, 200 and 50 values for these quantities linearly separated in the following ranges: For each resulting system, we solve the optimization problem of Eq. (23) for c opt .
Based on the optimization results, we can show the behavior of ln c opt versus ln m in Fig. 3(a) for the representative values of n. It is observed that this relationship can be expressed by the following linear function: where α 1 and β 1 are auxiliary variables. The points shown by markers • in Fig. 3(b) indicate the behavior of the slope β 1 as a function of n as we solve the optimization problem. It is observed that the curve in the form of n n+1 can be fitted to β 1 (solid line in Fig. 3(b)), resulting in ln c opt = ln α 1 + n n + 1 ln m, We note that n n+1 , the slope of ln m which has been found through curve fitting on the solution of the optimization problem, is the same as the corresponding term in Eq. (35). Therefore, the linear relationship between ln c opt and ln m in the form shown in Eq. (37) for given −δ , k, and e indicates that λ in Eq. (35) is independent of the effective mass. Furthermore, the dependency of c opt on n confirms that λ should be a function of this parameter.

Dependency of λ on −δ
In this investigation, we study the dependency of λ on −δ . We first assume that m = 1 kg, k = 1 N/m n , e = 0.7. We then sweep −δ and n by selecting, respectively, 200 and 50 values for these quantities linearly separated in the following ranges: 0.01 ≤ −δ ≤ 50 m/s, For each resulting system, we solve the optimization problem of Eq. (23) for c opt . Figure 4(a) shows the behavior of ln c opt as a function of ln ( −δ ) for the selected values of n. It is observed that the following linear function can describe the relationship between ln c opt and ln ( −δ ): In this expression α 2 and β 2 are auxiliary variables. The slope β 2 is shown in Fig. 4(b) for each system with a specified n using • markers. As is shown in this figure, n−1 n+1 is the best fit to these points, yielding It is observed that the slope of ln ( −δ ) is the same as the corresponding term in Eq. (35). As such, based on the linear relationship between ln c opt and ln ( −δ ) in the form presented in Eq. (40) for given m, k, and e, we conclude that λ in Eq. (35) is independent of −δ .

Dependency of λ on k
In previous investigations, we showed that λ is independent of −δ and m. Therefore, in order to study the dependency of λ on k, we assume that m = 1 kg, −δ = 1 m/s, and e = 0.7 in Eq. (35). Then we sweep k and n by, respectively, choosing 200 and 50 values, selected in the following ranges: 10 5 ≤ k ≤ 10 10 N/m n , For each resulting system, we solve the optimization problem in Eq. (23) for c opt . The representative results in Fig. 5(a) shows that ln c opt varies linearly as a function of ln k as where α 3 and β 3 are auxiliary variables. The value of the slope β 3 for each system with a specified n is shown in Fig. 5(b) with • markers. The best fit of 1 n+1 for β 3 results in ln c opt = ln α 3 + 1 n + 1 ln k, or, equivalently, As it is observed, the slope of ln k is the same as the one in Eq. (35). Therefore, the linear relationship between ln c opt and ln k in the form presented in Eq. (43) for given m, −δ , and e indicates that λ in Eq. (35) is independent of k.

Dependency of λ on e and n
So far our investigations have shown that λ is independent of −δ , m, and k. In order to investigate the dependency of λ on e and n, we assume that m = 1 kg, −δ = 1 m/s, and k = 1 N/m n . We then pick 600 points for e and 100 points for n linearly separated in the following ranges: 0.45 ≤ e ≤ 1, 1 ≤ n ≤ 2.
For each system, we solve the optimization problem of Eq. (23) for c opt . As shown in Fig. 6 for some representative values of n, the optimization results indicate that the behavior of c opt versus e obeys approximately the following pattern: where α and β are only functions of n and computed using the following procedure. First for each individual value of n ∈ [1, 2], we plot c opt vs. e (similar to the graph shown in Fig. 6), and compute the corresponding values of α and β. Figure 7 shows the resulting values of α and β as we sweep n, which can then be curve fitted using the following expressions: (48)

Closed-form formula for c opt
In previous sections, we concluded that λ is independent of m, −δ , and k. Furthermore, we found that λ can be expressed in the form shown in Eq. (48). In order to present an expression for c opt as a function of system parameters, we replace the function obtained for λ from Eq. (48) into Eq. (34), resulting in the following closed-form solution: where α and β are computed using Eqs.
where c opt is computed using Eq. (49). The value of n in the literature is mainly assigned to 3/2 which represents the contact between two spheres or a sphere and plane surface. Therefore, we explicitly express the This expression has also been reported in [25].

Numerical verification
In this section, we perform several numerical experiments to verify the validity and accuracy of the formula of Eq. (49).

Randomly selected values for the system parameters
The closed-form formula of Eq. (49) was derived through a series of numerical investigations where arbitrary values for mass, stiffness, and initial penetration velocity were assigned to the system parameters. In order to verify that the accuracy of this formula does not depend on the assigned values, we perform the following numerical investigation. We randomly assign values to each of the system parameters in the following ranges: 10 −5 ≤ m ≤ 10 2 kg, 10 5 ≤ k ≤ 10 10 N/m n , We also assign a value to the desired coefficient of restitution in the range 10 −4 ≤ e ≤ 1.
Then we determine the optimal value for c opt using Eq. (49). We solve the following equations governing the dynamics of the continuous contact model: mδ + c optδ + k δ n = 0, and stop the simulation at s t when the separation condition of c opt sδ + k ( s δ) n = 0 is met. Based on the simulated value of the separation velocity, we determine e out using Eq. (22), then the absolute error of the simulation is computed as |e out − e|.
We repeat the above simulation 300 times for different sets of randomly selected values of system parameters m, k, −δ , and n. For each given value of e, we determine the maximum absolute error. Figure 8 shows this upper bound of the absolute error versus e. It is observed that the application of the formula of Eq. (49) can accurately reproduce the desired restitution.

Optimal damping for n = 1
For n = 1, the system is reduced to the linear model of Kelvin-Voigt [8] presented by the following ODE: By applying the separation condition of Eq. (20) to the exact analytical solution of this model, Butcher and Segalman [2] have proven that the exact value of the coefficient of restitution can be related to the damping ratio ξ as where It should be noted that special consideration should be made for the correct determination of the quadrant of the argument of tan −1 when ξ < √ 2/2. Using n = 1, we can express the damping coefficient of Eq. (49) as In order to study the accuracy of Eq. (57), we assume that m = 1 kg, −δ = 1 m/s, and k = 1 N/m. We then sweep e in the range 10 −4 ≤ e ≤ 1, and use Eq. (57) to find the corresponding optimal damping coefficient. We further calculate the damping ratio for each resulting c opt using Eq. (56) and substitute it in Eq. (55) to compute the exact value of the corresponding coefficient of restitution. Figure 9(a) shows the behavior of e exact versus e. To better compare these values, we plot the absolute error |e exact − e| in Fig. 9(b). It is observed that the formula presented in Eq. (57) can accurately capture the damping properties of the system. Schematics of a two-pendulum system: the initial configuration of the system is shown by solid bodies and the configuration of the system at the first impact is shown by dotted bodies

Numerical simulation
In this section, we use the presented contact force model for the impact simulation of two pendulums and compare the results with the traditional modeling techniques. As shown in Fig. 10, each pendulum is modeled as a rod with negligible thickness and uniform mass distribution. We assume the following system parameters: m 1 = 1 kg, m 2 = 0.2 kg, L 1 = 6 m, L 2 = 4 m, and a = 2 m. The simulation starts with the configuration shown in Fig. 10 for the rods with solid boundaries, representing the initial condition θ 1 = π/12 rad, θ 2 = −π/2 rad, andθ 1 =θ 2 = 0 rad/s. The impact between the bodies are modeled using different techniques. First, the piecewise method with the application of the impulse-momentum formulation is used to address the dynamics of the system during impact [23]. Then, the continuous contact force model of Eq. (50), developed in this paper, is used to express the dynamics of impact during the period of contact. Finally, the Hunt-Crossley continuous contact model of Eq. (9) with the damping factor of Eq. (16) is used to simulate the dynamics of the system during the period of contact. Figure 11 shows the behavior ofθ 1 andθ 2 for these simulations while e = 0.3, n = 3/2, and k = 10 8 N/m 1.5 are used for the continuous contact models. It is observed that the system experiences three impacts during an 8-second simulation. The configuration of the system for the first impact is shown by bodies with dotted boundaries in Fig. 10. Furthermore, Fig. 11 shows that the presented method in this paper can capture the

Conclusion
In this paper, we have presented an analytical expression to compute the optimal damping coefficient in a class of nonlinear continuous force models. The contact force model consists of a spring-damper force element in which the spring obeys force-penetration characteristics of the Hertz model. Furthermore, it has been assumed that the exponent of the deformation in the Hertz model can be assigned values between one and two. We have applied an optimization process to compute the damping coefficient of the linear damper in the contact model such that the desired restitution is reproduced at the time of separation. Furthermore, unlike most of other reported models in literature, in which the separation condition is considered at the instant at which both penetration and contact force are zero, in this work we have considered zero contact force as the only condition for separation. After separation, either the restitution continues or permanent indentation may have been reached. As such, the presented model covers a wider range of impact problems. We have verified the validity and accuracy of the proposed formula through conducting several numerical simulations on systems with randomly-assigned parameters. Additionally, in the case of the Kelvin-Voigt model, it has been shown that the proposed formula can successfully recover the exact solution of this model. Finally, we have applied the presented formulation to demonstrate the impact of a two-pendulum system and compared the results against those of the piecewise method and the Hunt-Crossley-based continuous model.