Further investigation on improved viscoelastic contact force model extended based on hertz's law in multibody system

Highlights

• The shortcoming of the existing improved contact stiffness coefficients is analyzed systematically.

• A new contact force model suited for the internal contact between the cylindrical bodies is proposed.

• The other contact force model suited for the external contact between the cylindrical bodies is proposed.

• The correctness of the external contact force model is validated by the bouncing example, the reasonability of the internal contact force model is proved by the experimental data.

ABSTRACT

Since the Hertzian contact stiffness coefficient is closely related to the power exponent, the modification of the Hertzian contact stiffness coefficient leads to some problems: (i) the power exponent must be changed with the improved contact stiffness coefficient, otherwise, the unit of the contact force is not newton; (ii) since the power exponent is one of the factors to determine the magnitude of the contact force, the inappropriate power exponent makes the multibody system hard to convergence or lose physical meaning. The main goal of this investigation is to reveal the shortcoming of the improved contact stiffness coefficient, and develop two different contact force models suited for the internal and external contact form associated with the hysteresis damping factor from Lankarani-Nikravesh contact force model (L-N model). The effectiveness of the external contact force model is validated by the bouncing example. The correctness of the internal contact force model used in the multibody system with a clearance revolute joint is validated by the experimental data.

Introduction

The contact-impact event is a complicated physical phenomenon that happens in a very short duration, which owns a lot of conspicuous behaviors including high force level, rapid energy dissipation, and jumps in the velocities, etc. [1]. In order to depict these behaviors during the contact process, the normal contact force model plays a crucial role in the modeling of the contact phenomenon [2]. In general, the contact forces generated by the contact-impact event can be expressed as a function of the penetration depth by smoothing the discontinuity of the impact force. Most of the existing contact force models are developed based on Hertz's law [1,3]. There are two primary approaches that can be used to develop the contact forces in the multibody system [4]: (i) non-smooth dynamics formulation based on the geometric constraints [5]; (ii) regularized approach belongs to the continuous methods [2,4,[6], [7], [8].

The non-smooth system usually displays nonlinearities or discontinuities, such that the contact bodies are treated as an ideal rigid body that cannot deform in the contact process [9], or the contact areas cannot be penetrated into each other. The complementarity approach [9], [10], [11], [12] is one of the most important methods to resolve the contact dynamics problem using unilateral constraints. However, this formulation has some challenges which are hard to be fixed. Namely, it is inconvenient to cooperate with the friction force model to calculate the contact forces that occurred between the contact bodies; secondly, this approach is prone to violate the energy conservation during contact with friction [13]. Recently, Hu et al. proposed the structure-preserving approach [14,15] to reproduce the nonlinear characteristics for the non-smooth infinite-dimensional dynamic systems [16], [17], [18], which owns excellent long-time numerical stability and high-accurcy in the nuermical simulation.

As for the regularized approach, it can be treated as a series of spring-damper elements covering on the surfaces of the contact bodies, therefore the normal contact force incorporates two terms which are the elastic and damping force. Thereby, the approach is referred to as the penalty or compliant method [19]. Being a sharp contrast to the complementarity approach, the normal contact force can be described as a continuous function of the penetration depth between the contact bodies [20], that is, the magnitude of the contact force depends on the detection and estimation of the penetration depth. This characteristic brings some challenges to the system; the high-frequency components are introduced since the contact happens in a very short period [21], which forces the integration algorithm to take a very small time step for the sake of capturing the high-frequency components [22] and in consequence cause low computational efficiency. Moreover, the other shortcoming of this approach is hard to choose the contact parameters, including the stiffness coefficient, the coefficient of restitution, and the degree of nonlinearity of the penetration depth, etc. [23]. Considering this approach can accurately depict the energy dissipation using the hysteresis damping factor in the contact process, therefore, this method is widely used in some commercial software to tackle the contact problems.

As we know, the Hertz's contact law is obtained through the experimental contact test of two spheres [24], and the general formulation is $F_n=K{\delta }^n$, when the power exponent is equal to 1, this formulation is similar to the Hooke's law, but the essence between them is significantly different since the Hertzian contact stiffness includes the material and geometric properties of contact bodies. However, the stiffness coefficient of the Hooke's law is a measurement of the spring, which has completely different physical meaning from the Hertzian contact stiffness coefficient. Some scholars consider that the Hertzian contact law is Hooke's law when the power exponent is equal to 1 [25], which makes a serious problem. Because the power exponent is determined by the contact material, contact form, and distribution status of the contact stress, the value of the power exponent can be regarded as 1 in a causal manner [22]. Hence, it must be clear that understanding the relationship between the Hertzian stiffness coefficient and power exponent is very important for the improvement of the Hertzian stiffness coefficient [26,27].

Considering that the original Hertzian contact force model is a pure elastic model that cannot estimate the energy dissipation in the contact process [3], in order to constantly approach the fidelity of the contact event occurred in the multibody system, the original Hertzian contact force model is improved by imposing the different hysteresis damping factors [28], [29], [30], or by modifying the Hertzian contact stiffness coefficient. Therefore, these improved methods can be classified into two principal ideas: (i) enhancing the accuracy of the damping coefficient D for the sake of describing a more precise energy dissipation [31]. (ii) improving the physical property of the contact stiffness coefficient K in order to generate the relationship with penetration depthδ, because the original Hertzian contact stiffness coefficient is a constant value [22]. Nonetheless, most investigations adopted the first method to improve the normal contact force model by using different ways to approximate the energy dissipation during the contact process [32]. In order to get a better understanding of the normal contact force model, it is necessary to introduce the physical meaning of each contact parameter in detail. To this end, the general form of the normal contact force F_n can be expressed as $F_n=K{\delta }^n+D\dot{\delta }\left(D=\chi {\delta }^n\right)$. The first term on the right-hand side represents the elastic force components; the second term explains the energy dissipation. $\dot{\delta }$ denotes the normal contact velocity obtained by differentiating the constraint equation of the potential contact points between the contact bodies, n is the power exponent, χis the hysteresis damping factor. First of all, regarding the second term, the damping coefficient is a function of the penetration depth δand hysteresis damping factor χ. The research progress of the improved normal contact force model is the development process of the hysteresis damping factor that includes the coefficient of restitution c_r, stiffness coefficient K, and initial contact velocity${\dot{\delta }}^{(-)}$. Lankarani and Nikravesh indicated [33] that the initial contact velocities should smaller than the propagation velocity of elastic waves across the two colliding bodies when the normal contact force model extended from Hertzian contact law are valid, that is, ${\dot{\delta }}^{(-)}\le {10}^{-5}\sqrt{E/\rho }$, where E is Young's modulus and ρ is the density of the contact bodies. However, when the initial contact velocity is larger than ${10}^{-5}\sqrt{E/\rho }$ [22], there are some permanent penetrations left behind on the surfaces of the contact bodies, which accounts for the energy loss in contact [34]. The coefficient of restitution c_r is a controller for the energy dissipation, which is a constant value that represents the proportionality between the pre-impact velocity and the post-impact velocity using newton's law of restitution [1,35]. During the impact process, the initial kinetic energy from the contact bodies has three different destinations [24,36]: (i) one part is the kinetic energy of the system at the end of the compression phase or the beginning of the restitution phase; (2) the second part is the elastic strain energy stored in the contact bodies; (3) the last part is dissipated by the vibration, sound, and heat, etc. When the coefficient of restitution c_r is equal to one, the hysteresis damping factor is equal to zero [35]. All the normal contact force models are back to the original Hertzian contact force model because there is no energy to be dissipated. However, when c_r is equal to zero, the hysteresis damping factor is infinite that represents the full plastic contact [24]. Namely, two contact bodies become one body when the contact occurred in a short time. Therefore, as for the usual oblique collision, the value of the coefficient of restitution should reside in the range from zero to one. For the high value of the coefficient of restitution, most of the energy is stored as the elastic strain energy. However, for the low value, most of the energy is dissipated [37].

Since the Hertzian stiffness coefficient is a constant value, which is independent with the penetration depth. In order to introduce the penetration depth into the contact stiffness coefficient, there are three different methods to enhance the accuracy of the first term in this equation $F_n=K{\delta }^n+D\dot{\delta }$: (i) differentiating the normal contact force model with respect to the penetration depth [22,27,38,39], or modify the Hertzian stiffness coefficient based on the contact mechanics. (ii) in order to describe contact status between cylindrical contact bodies, the length parameter of the cylindrical bodies is imposed on the stiffness coefficient to achieve a more accurate contact force model [40]. (iii) in order to depict the energy dissipation in a more precise way, the power exponent is modified by exponent fit function [41]. Unfortunately, these investigations ignored a very important implicit relationship between the power exponent and contact stiffness coefficient, in other words, the power exponent must keep consistent with the unit of the contact stiffness coefficient for the sake of ensuring that the product of K and δ_n has a force unit. If the Hertzian contact stiffness coefficient keeps its original formulation, the value of the power exponent cannot be fitted in any approach. Otherwise, the improved contact force model by fitting the power exponent violates the physical meaning of the contact mechanics because the unit of the contact force is not newton. In general, the power exponent is equal to 1.5 in the Hertzian law. As for the power exponent, it is necessary to emphasize at least three points: (i) when the power exponent is equal to 1.5 for metallic materials, it indicates the contact pressure is distributed as a parabolic shape [42]; (ii) when the contact material is change to glass or polymers, etc. The distribution of the contact pressure is not parabolic, and consequently the power exponent is not equal to 1.5 [20]; (iii) since the original Hertzian contact law is suitable for the point contact [43], [44], [45], [46], [47], such that the power exponent should be modified accordingly when the bodies contact in a line or surface [48], [49], [50], [51]. However, most of the existing publications on the revolute joint with clearance [50,[52], [53], [54], [55], or cylindrical joint with clearance [56], [57], [58] do not follow this law, but they can still obtain a reasonable solution in the multibody system [23,46,47,52,[58], [59], [60], [61].

In light of the investigations above, the contact force model with an improved contact stiffness coefficient has some limitations as follows:

• From the physical point of view, when remaining the Hertzian contact stiffness coefficient, the power exponent can not be changed, otherwise, the unit of the first term in the contact force model is not newton. Simultaneously, the unit of the damping force is not newton as well since the damping coefficient is a function of the contact stiffness coefficient and then penetration depth associated with the power exponent. The normal contact force model loses its physical meaning and violates the foundation theory of the mechanics.

• From the numerical analysis point of view, once the power exponent is changed in order to match the unit of the improved contact stiffness coefficient [26], the magnitude of the contact force becomes an abnormal value since the term δⁿplays a decisive role despite the penetration depth is reasonable. These abnormalities lead to a meaningless contact force. This situation not only causes the integration process is easy to divergence but also makes the computation inefficient. When the contact force is large, the integration algorithm either reduces the time step until capturing the large amplitude of the vibration or directly go to divergence [20]. When the contact force is small, the contact deformation will not be detected, which makes a non-ideal situation tend to the ideal one.

Considering the large number of publications related to the contact force model [62], [63], [64], [65], [66], [67], [68], [69], [70],22,42,[77], [78], [79],54,63,[71], [72], [73], [74], [75], [76], in order to make the future investigation move to the correct direction, this investigation should attract sufficient attention so as to figure out how to improve Hertzian contact stiffness coefficient in a rational manner, and propose two contact force models suited for the line contact form. The main contributions of this investigation can be summarized as follows:

• The relationship between the Hertzian contact stiffness coefficient and the power exponent is systematically studied by a comparative analysis. The drawbacks of the improved contact force model with an improved stiffness coefficient are indicated by a series of numerical simulations. The simulation results account for that the improved contact stiffness coefficient must keep consistent with the power exponent for the sake of ensuring the unit of the improved contact force model is the newton. Even in some cases, the improved contact force models can be used to obtain the accepted solution, but these models lose the original physical meaning of the contact mechanics.

• Two contact force models with energy dissipation suited for internal and external contact are tailored to the impact between the cylindrical bodies. That is mainly because the two new contact force models not only introduce the length of the joint but also make the power exponent equal to 1. More importantly, two contact force models can describe all existing contact force models when their hysteresis damping factor is replaced by the corresponding contact force model.

The structure of this investigation can be organized as follows: In Section 2, the contact force models with energy dissipation are simply introduced. In Section 3, the contact force models are studied when the power exponent is not equal to 1.5. In Section 4, some contact force models with the improved contact stiffness coefficient are introduced in detail. In Section 5, two contact force models suited for internal and external contact are proposed, the correctness of the external contact is verified using the bouncing example. In Section 6, the slider-crank mechanism with a clearance revolute joint is taken as an example to account for the shortcomings of the improved contact force models and prove the reasonability of the internal contact force model. The simulation results are validated by the experimental data. The main conclusions are summarized in Section 7.