Modeling and verification of an improved contact force in multibody systems

Abstract

In many engineering applications, contact and collision phenomena are commonly observed in multibody systems, which can lead to issues such as dynamic output oscillations, reduced motion accuracy, decreased reliability and lifespan, and even functional failures in mechanical systems. To more accurately describe common collision phenomena and their impact on the dynamic characteristics of multibody systems, this paper introduces a contact force model with improved nonlinear stiffness and damping coefficients. This model is based on Hertz theory and considering the relationship between indentation and velocity during the collision process by incorporating a nonlinear index factor, m. Additionally, the Newton restitution coefficient is used as an evaluation criterion to verify the effectiveness of the improved model. Numerical analyses were conducted on single collisions between joints at different initial collision velocities and restitution coefficients, using both the proposed improved contact force model and existing classical models. Further numerical simulations of steel ball collisions were performed using the improved model, and the results were compared with experimental test data. Ultimately, the study results demonstrate that the improved contact force model is more accurate and effective than existing classical models under various restitution coefficients and collision velocity conditions, and has a wider range of applicability.

Introduction

Contact collision phenomena are commonly present in multibody systems and can lead to issues such as vibration, wear, and fatigue, which disrupt the continuity of the system structure and affect dynamic accuracy, reliability, and lifespan. Dynamic parameters during collisions, such as contact force, velocity, and indentation, are crucial for accurate and effective modeling of contact forces. These parameters are foundational for modeling multisystem dynamics, studying dynamic impact patterns, and evaluating the performance of systems.^1–3 Therefore, contact modeling and analysis hold significant theoretical reference and engineering application value in many fields, including granular materials, ⁴ mechanical systems, ⁵ collision analysis, ⁶ biomechanics, ⁷ robotics, ⁸ vibrational impact drilling, ⁹ and others. ¹⁰

Typical application scenarios include the landing and docking processes of the Chang’e-5 lunar mission spacecraft, collision tests in vehicle development, and frictional contact between tires and the ground. With technological advancements, various systems are evolving toward lightweight, precise, high-speed, and high-load capabilities. ¹¹ In collision events, instantaneous impact forces and elastoplastic deformation can significantly affect the stability or structural integrity of multibody systems, thereby posing potential threats to the safe use of high-precision equipment. ¹²

Therefore, for the safe and stable operation of engineering mechanical systems, it is essential to study the physical mechanisms of impact and vibration in multibody systems triggered by collision events, and to understand the coupled relationships between collisions and the overall dynamic performance of the system. ¹³ However, research into collision behaviors in multibody systems is closely intertwined with the development of contact force models. ¹⁴ Collision force models serve as critical criteria for evaluating whether collision behaviors could lead to structural damage in multibody systems, underscoring the crucial need for their accuracy. ¹⁵ The dynamics of multibody systems involve the study of motion constraints and interaction forces within mechanical systems. Contact, friction, and impact events are among the complex and highly significant phenomena in multibody systems. However, accurately assessing these phenomena hinges on the geometric shapes of contacting surfaces, material properties, and the constitutive models that simulate contact responses. ¹⁶ Contact dynamics modeling primarily involves identifying potential contact points between colliding objects within multibody systems, assessing contact-impact forces, and constructing transitions between different states of contact and non-contact. ¹⁷ Analyzing contact and collision motions in multibody systems remains one of the most challenging and complex problems in the fields of science and engineering. ¹⁸

Since the 1970s, many researchers around the world have conducted extensive studies on the dynamics of multibody systems with contact joints, and the core issues being the modeling and solving of system dynamics with collision. ¹⁹ The aims of this paper are to review and compare existing classical contact force models from three aspects: surface geometry, material properties, and simulated contact responses. By introducing a nonlinear index factor denoted as m, this paper established an improved external contact force model. Through numerous numerical simulations of two balls, based on the new model and six existing classical contact models, and by comparing the results with the collision experiments of steel balls in the reference, the effectiveness and accuracy of the improved contact force model presented in the paper are verified.

The structure of the remaining sections of this paper is outlined as follows. Section “Geometry of contact in multibody dynamics” provides a brief overview and model construction of intra-joint and inter-joint contact collisions in multibody systems. In Section “Contact force models of joints,” subsection “Classical contact force models” first discusses the development of classical contact force models, ranging from the non-dissipative Hertz contact force model to Hunt and Crossley’s nonlinear spring-damping model, and includes improved models adaptable to various conditions, such as those by Flores et al. Section “Improved contact force model” builds upon the analysis and synthesis of these models to establish an equivalent contact model, deriving and constructing the improved model used in this study. In Section “Comparison and analysis of models,” Subsection “Application 1” uses a set of joint models as an example, employing the coefficient of restitution as an evaluation criterion. It initially validates the effectiveness of the improved contact force model by varying collision velocities. Subsequently, different restitution coefficients are set to analyze the relative errors in actual collision restitution coefficients between the improved model and six other models (Hertz, Liu, H-C, L-N, Gonthier, and Flores), thereby confirming the accuracy of the improved model. Section “Application 2” introduces the classical pendulum steel ball collision test model, comparing and analyzing the restitution coefficients of the improved model proposed in this paper with those obtained experimentally, thereby revalidating the accuracy of the improved contact force model. Section “Conclusion” discusses the main conclusions of this study, establishing a theoretical foundation for joint contact issues in multibody systems under varying restitution coefficients and loading conditions.

Geometry of contact in multibody dynamics

The modeling accuracy and efficiency of contact problems in multibody dynamics primarily depend on the complexity of the contact surfaces, ²⁰ the number of potential colliding bodies, ²¹ and the kinematics of the objects. ²² The contact geometry in multibody dynamics consists of three fundamental aspects ²³ : (i) the geometric description of the contact surfaces; (ii) the identification of potential contact points; and (iii) the evaluation of contact kinematics. Therefore, accurately describing contact collision phenomena is crucial for contact force modeling. The main focus of this study is to model the collision force between two contacting bodies in a multibody system. Figure 1 illustrates the collision model of two objects.

Figure 1.

Collision model: (a) δ > 0, (b) δ = 0, and (c) δ < 0.

Assume that two objects undergo deformation δ during the collision process. δ is the relative penetration depth, typically representing the amount of deformation, which is related to factors such as material properties and collision velocity. Figures 1(a) to (c) show three different scenarios involving two generic contacting bodies, where δ takes on three distinct values, ²⁴ allowing for the identification of contact and non-contact states between the two objects. Figure 1(a)δ > 0, refers to the state where the two colliding bodies are in a non-contact situation; (b) δ = 0, represents the instant of contact initiation; (c) δ < 0, indicates that the two colliding bodies are in contact. The expression for the contact state is as follows:

$$f\left(x\right)=\{\begin{array}{c}\delta >0,\mathrm{Not}\mathrm{contact}\hfill \\ \delta =0,\mathrm{Just}\mathrm{contact}\hfill \\ \delta <0,\mathrm{Already}\mathrm{contact}\hfill \end{array}$$

(1)

Contact force models of joints

Classical contact force models

As the pioneer of contact force modeling, Hertz ²⁵ proposed a fully elastic penetration contact force model suitable for non-conforming contact surfaces, which is outlined as follows

$$F=K{\delta }^n$$

(2)

In the equation, F represents the normal contact force, and K represents the equivalent contact stiffness. The generalized stiffness depends on the material properties and the geometry of the contacting bodies. For two spheres in contact, the generalized stiffness coefficient is a function of the radii of the spheres i and j, as well as the material properties as ²⁶

$$K=\frac{4}{3\pi ({\sigma }_i+{\sigma }_j)}\sqrt{\frac{R_i{R}_j}{R_i+R_j}}$$

(3)

Here, n represents the nonlinear power exponent (n ≥ 1), which is determined by the material and geometric properties of the local area of the contacting bodies. Researchers like Zhang Jie have applied the law of energy conservation to develop a continuous contact force model that accommodates arbitrary values of penetration depth and velocity exponent. Their study results indicate that by choosing appropriate deformation and velocity exponents, the performance of the continuous contact force model can be enhanced. ²⁷ To more accurately simulate the diversity of different materials, this paper proposes a variable nonlinear index, denoted as m, and establishes an improved contact force model incorporating this variable nonlinear index.

In the equation (3), the stress ${\sigma }_k$ can be expressed as

$${\sigma }_k=\frac{1-v_k^2}{\pi E_k}k=i,j$$

(4)

where $E_k$ and $v_k$ are the elastic modulus and Poisson’s ratio of the material.

In order to overcome the limitations of the Hertz model, which is only applicable to large clearances, small loads, and low collision speeds, Liu et al.^28,29 proposed a normal contact force model suitable for conformal contact surfaces. This model describes the relationship between the external load P and the elastic indentation δ at any contact surface, represented as

$$P=\frac{1}{2}\pi E^{*}\delta \sqrt{\frac{\delta }{2(\mathrm{\Delta }R+\delta )}}$$

(5)

where $E^{*}$ is the equivalent elastic modulus, and ΔR is expressed as the difference in radius between the two colliders.

All contact force models above are non-damping elastic penetration models. However, energy dissipation occurs in actual contact collisions. Therefore, Kelvin-Voigt ³⁰ and others proposed a contact force model considering damping consumption, as follows

$$F_N=K\delta +D\stackrel{·}{\delta }$$

(6)

where $\stackrel{·}{\delta }$ represents the normal relative velocity during the collision process, D is the damping coefficient, which is independent of the restitution coefficient. Since D is a nonzero constant, at the instant of collision initiation, the indentation is zero while the contact force is greater than zero. The Kelvin-Voigt contact force model is generally not suitable for high-speed collisions, as it primarily describes elastic and viscous behavior at low speeds. In high-speed collisions, the inertial effects and nonlinear characteristics of the materials become significant. Additionally, throughout the entire collision and recovery process, the damping coefficient remains constant, which is inconsistent with actual working conditions.

To overcome the limitations of the Kelvin-Voigt model, Hunt and Crossley (H-C) ³¹ proposed a nonlinear spring-damping model, expressed as

$$F_N=K{\delta }^n[1+\frac{3(1-C_r)\stackrel{·}{\delta }}{2{\stackrel{·}{\delta }}^{(-)}}]$$

(7)

where ${\stackrel{·}{\delta }}^{(-)}$ represents the initial collision velocity, and $C_r$ is the restitution coefficient.

Based on the H-C model, Lankarani and Nikravesh ³² considered factors such as material damping, geometric shape, elastic penetration, and collision velocity, and established a contact force model applicable to a contact area that is small relative to the contacting object, denoted as

$$F_N=K{\delta }^n[1+\frac{3(1-{C_r}^2)\stackrel{·}{\delta }}{4{\stackrel{·}{\delta }}^{(-)}}]$$

(8)

Gonthier et al. ³³ proposed an improved contact force model based on the H-C model, expressed as

$$F_N=K{\delta }^n[1+\frac{(1-{C_r}^2)\stackrel{·}{\delta }}{C_r{\stackrel{·}{\delta }}^{(-)}}]$$

(9)

Flores et al. ³⁴ regarded the contact collision process as a single-degree-of-freedom dynamic system. Assuming that internal damping of the material causes energy dissipation during the collision process, the established model can be represented as

$$F_N=K{\delta }^n[1+\frac{8(1-C_r)\stackrel{·}{\delta }}{5C_r{\stackrel{·}{\delta }}^{(-)}}]$$

(10)

In addition to the classical contact force models mentioned above, numerous scholars have proposed various collision contact force models applicable to different working conditions.^35,36

Improved contact force model

According to Hertz theory, when collisions occur between joint surfaces with different radii, the indentation area can be approximated as two spherical caps with equal base areas but different heights. As shown in Figure 2(a), the shaded part represents the collision contact volume.

Figure 2.

External contact collision model and its equivalent model: (a) external contact model, (b) the equivalent external contact model, (c) actual contact indentation area, (d) equivalent 2D diagram of the contact area, and (e) equivalent elliptic cylinder model.

In order to accurately simulate the collision forces between two contacting bodies in a multi-body system, this study proposes an improved contact force model.

Based on the collision theory of spheres with arbitrary radii from Hertz, ³⁷ an equivalent contact model using two spheres of the same radius is proposed, as shown in Figure 2(b). Where the computational model of equivalent radius $\overline{R}$ can be expressed as

$$\overline{R}=R_1{R}_2/(R_1+R_2)$$

(11)

The indentation on both sides of the sphere in the equivalent model is the same, with a height of $\overline{h}$ . The mathematical model is as follows

$$\mathrm{OC}=\overline{h}=\delta /2$$

(12)

The volume of the actual contact penetration zone in Figure 2(c) is the sum of the volumes of two spherical caps with radius $\overline{R}$ and height $\overline{h}$ . The volume $V_{\mathrm{Lack}}$ can be presented as

$$V_{\mathrm{Lack}}=\frac{2}{3}\pi {\overline{h}}^2(3\overline{R}-\overline{h})$$

(13)

The equivalent two-dimensional front view of the contact penetration zone interface is shown in Figure 2(d). The cross-section of the contact penetration zone is equivalent to an ellipse with a short axis OC and a long axis OA. The calculation model for OA can be indicated as

$$\mathrm{OA}=\sqrt{{\overline{R}}^2-{(\overline{R}-\delta /2)}^2}=\sqrt{\delta (4\overline{R}-\delta )}/2$$

(14)

As shown in Figure 2(c) and (e), based on the equivalence of the actual contact penetration zone model to an elliptic cylinder model, this study proposes an improved contact force model. Here, L represents the height of the equivalent elliptic cylinder.

Combining equations (12) and (14), the volume $V_{\mathrm{Equal}}$ of the equivalent contact indentation elliptic cylinder can be obtained as

$$V_{\mathrm{Equal}}=\pi ·\mathrm{OA}·\mathrm{OC}·L=\frac{1}{2}\pi ·\sqrt{\delta (4\overline{R}-\delta )}·\overline{h}·L$$

(15)

As shown in Figure 2(c) and (e), $V_{\mathrm{Lack}}$ is defined as the sum of the volumes of the two spherical caps and represents the volume of the actual contact penetration zone. $V_{\mathrm{Equal}}$ represents the volume of the equivalent deformed elliptic cylinder. From this, an equilibrium equation (16) can be derived.

$$V_{\mathrm{Lack}}=V_{\mathrm{Equal}}$$

(16)

By combining equations (13), (15), and (16), the equivalent height of the elliptic cylinder L can be expressed as

$$L=\frac{6\overline{R}-\delta }{3}\sqrt{\frac{\delta }{4\overline{R}-\delta }}$$

(17)

Without considering damping, by combining formula (5) and (17), the equivalent contact force $F_e$ between colliding bodies is obtained as follows

$$F_e=\underset{0}{\overset{L}{\int }}\mathrm{PdL}=\frac{1}{6}\pi E^{*}{\delta }^2\frac{6\overline{R}-\delta }{\sqrt{2(4\overline{R}-\delta )(\mathrm{\Delta }R+\delta )}}$$

(18)

In order to describe the nonlinear properties of materials more accurately, the nonlinear index factor m, which is related to the penetration of different materials in contact collision, is introduced in this paper, which is different from the fixed nonlinear index n in the existing studies. Based on Hertz theory, the force balance formula is obtained

$$\frac{1}{6}\pi E^{*}{\delta }^m\frac{6\overline{R}-\delta }{\sqrt{2(4\overline{R}-\delta )(\mathrm{\Delta }R+\delta )}}=K{\delta }^n$$

(19)

where n is the general material coefficient, and the improved nonlinear exponential factor m is

$$m=n+{\log }_{\delta }\frac{6K}{\pi E^{*}}\sqrt{\frac{2(4\overline{R}-\delta )(\mathrm{\Delta }R+\delta )}{{(6\overline{R}-\delta )}^2}}$$

(20)

Without considering damping, the improved nonlinear stiffness coefficient $K_i$ can be expressed as

$$K_i=\frac{F_e}{{\delta }^n}=\frac{1}{6}\pi E^{*}{\delta }^{m-n}\sqrt{\frac{{(6\overline{R}-\delta )}^2}{2(4\overline{R}-\delta )(\mathrm{\Delta }R+\delta )}}$$

(21)

Then, based on the nonlinear damping coefficient formula proposed in Wang et al., ³⁸ the stiffness coefficient is replaced with the variable stiffness coefficient improved in this study, as shown in equation (21), resulting in the enhanced nonlinear damping factor $D_i$ . It can be expressed as

$$D_i=\frac{3(1-{C_r}^2)e^{2(1-C_r)}{\delta }^n}{4{\stackrel{·}{\delta }}^{(-)}}{K}_i$$

(22)

By combining equations (21), (22), and equation (6) containing nonlinear indices, an improved contact force model for this study is obtained, which is expressed as

$$\begin{array}{c}\multicolumn{1}{c}{F_N=K_i{\delta }^n+D_i\mathit{\delta }=\frac{1}{6}\pi E^{*}{\delta }^m\sqrt{\frac{{(6\overline{R}-\delta )}^2}{2(4\overline{R}-\delta )(\mathrm{\Delta }R+\delta )}}}\\ +\frac{3(1-{C_r}^2)e^{2(1-C_r)}{\delta }^n}{4{\stackrel{·}{\delta }}^{(-)}}\frac{1}{6}\pi E^{*}{\delta }^{m-n}\sqrt{\frac{{(6\overline{R}-\delta )}^2}{2(4\overline{R}-\delta )(\mathrm{\Delta }R+\delta )}}\stackrel{·}{\delta }\hfill \\ =\frac{1}{6}\pi E^{*}{\delta }^m\sqrt{\frac{{(6\overline{R}-\delta )}^2}{2(4\overline{R}-\delta )(\mathrm{\Delta }R+\delta )}}1+\frac{31-{C_r}^2{e}^{2(1-C_r)}}{4{\stackrel{·}{\delta }}^{-}}\stackrel{·}{\delta }\hfill \end{array}$$

(23)

Comparison and analysis of models

The restitution coefficient is an important parameter for evaluating the accuracy of contact force modeling. ³⁹ It is a global metric used to measure energy loss during collision and is widely applied in engineering and theoretical research. ⁴⁰ Common types of restitution coefficients include kinematic, dynamic, and energy restitution coefficients.^41–43 These can be further categorized into four types: Newton, Poisson, Stronge, and Hedrih restitution coefficients. Among these, the Newton restitution coefficient is suitable for various collision scenarios, including elastic, inelastic, and mixed-type collisions. Consequently, it has been extensively applied in fields such as engineering, physics, and kinematics. The Newton restitution coefficient is both intuitive and practical, in actual engineering and physical problems, it simplifies the calculation process and improve the solution efficiency. ⁴⁴ According to the law of momentum conservation, the Newton restitution coefficient is defined as the ratio of the relative velocities of objects before and after a collision, where ${\stackrel{·}{\delta }}^{(-)}$ represents the initial velocity before collision, and ${\stackrel{·}{\delta }}^{(+)}$ represents the velocity after collision. The Newton restitution coefficient can be expressed as

$$C_r\prime =\frac{{\stackrel{·}{\delta }}^{(+)}}{{\stackrel{·}{\delta }}^{(-)}},0\le C_r\prime \le 1$$

(24)

Application 1

To better compare and analyze various models, this study uses a set of sphere collision models as an example, as shown in Figure 3. Numerical simulations are conducted on parameters such as contact force, penetration depth, and collision velocity at the contact surface.

Figure 3.

Sphere collision model.

Assuming the material of the joints is metal and the nonlinear index factor n is 1.5.^31,44 The object 2 is fixed while a positive collision occurs between the object 1 and object 2. The mass of each object is 0.04 kg. The radius $R_1$ of object 1 is 9.9 mm, the radius $R_2$ of object 2 is 10 mm, and the radial difference c between the two spheres is $0.1\mathrm{mm}0.1\mathrm{mm}$ . Additionally, the Poisson’s ratio σ of the material is 0.3, the elastic modulus E ₁ and E ₂ are both 207 GPa.

To initially verify the effectiveness of the improved model, the initial velocities of the object 1 are assumed to be 0.5, 5, and 50 m/s, with restitution coefficients of 0.3, 0.5, and 0.8. Table 1 presents the actual restitution coefficients $C_r\prime$ and the relative error results of the improved model.

Table 1.

Actual restitution coefficient and relative error of the improved model at different initial velocities and restitution coefficients.

$C_r$	$C_r{\stackrel{·}{\delta }}^{(-)}(\mathrm{m}/\mathrm{s})$	0.5	5	50
0.3	${\stackrel{·}{\delta }}^{(+)}/(\mathrm{m}/\mathrm{s})$	0.1633	1.6335	16.3347
	$C_r\prime$	0.3266	0.327	0.3267
	Error/(%) (%)	8.6667	9	8.9
0.5	${\stackrel{·}{\delta }}^{(+)}/(\mathrm{m}/\mathrm{s})$	0.2413	2.4130 m/s.	24.1301
	$C_r\prime$	0.4826	0.4826	0.4826
	Error (%).	3.48	3.48%	3.48%
0.8	${\stackrel{·}{\delta }}^{(+)}/(\mathrm{m}/\mathrm{s})$	0.3937	3.9367	39.3685
	$C_r\prime$	0.7874	0.7873	0.7874
	Error (%).	1.575	1.5875	1.575

According to Table 1, it can be concluded clearly that: (1) With the same restitution coefficient, as the initial collision velocity increases by factors of 10 and 100 times from 0.5 m/s, the separation velocity of the improved model increases by factors of 10 and 100, respectively, while the actual restitution coefficient and relative error remain relatively stable. (2) Across different restitution coefficients, the relative error of the improved model gradually decreases as the restitution coefficient increasing. When the restitution coefficient is 0.3, the relative error is approximately 8.9%; when it’s 0.5, the relative error is approximately 3.48%; and when it reaches 0.8, the error is only 1.575%.

The results indicate that the error of the improved model is almost unaffected by collision speed under different impact velocities. Under various restitution coefficient conditions, the error of the improved model remains below 10% and gradually decreases as the restitution coefficient increases. In conditions with a high restitution coefficient, the improved model demonstrates very high accuracy, with an error of only 1.575%.

To further analyze the differences between the improved model and existing contact force models, this paper assumes an initial velocity of 0.5 m/s and restitution coefficients is 0.3, 0.5, and 0.8, respectively. The relationship between indentation and contact force for seven different contact force models is obtained through a great number of numerical simulations, as shown in Figure 4(a) to (c). Specific data on the maximum contact force and indentation are listed in Tables 2 and 3. It should be noted that the models of Hertz and Liu do not have a restitution coefficient, so in the subsequent analysis regarding the restitution coefficient, these two models are not mentioned, and their restitution coefficient errors are defined as 0.

Figure 4.

Indentation versus contact force relationship graph: (a) C _r  = 0.3, (b) C _r  = 0.5, and (c) C _r  = 0.8.

Table 2.

The maximum contact force of different contact t force models at three restitution coefficient values.

Contact force model		Hertz	Liu	H-C	L-N	Gonthier	Flores	Author
Maximum contact force/(N)	$C_r=0.3$	740.9571	811.4666	674.8040	669.4002	768.2244	803.5917	754.9717
	$C_r=0.5$	740.9571	811.4666	669.3942	670.9966	691.5379	696.360	693.006
	$C_r=0.8$	740.9571	811.4666	685.9138	689.0489	675.0241	677.8196	677.6132

Table 3.

The maximum indentation of different contact force models at three restitution coefficient values.

Contact force model		Hertz	Liu	H-C	L-N	Gonthier	Flores	Author
Maximum indentation/(×10−5 m)	$C_r=0.3$	1.6867	1.54	1.3774	1.4595	1.1162	1.0594	1.1400
	$C_r=0.5$	1.6867	1.54	1.4430	1.4910	1.2973	1.2819	1.2928
	$C_r=0.8$	1.6867	1.54	1.5703	1.5086	1.5230	1.5382	1.5374

As shown in Figure 4, the contact force models proposed by Hertz and Li do not include a damping term because they do not account for energy dissipation. This results in the absence of a closed-loop region for energy loss during the collision process. However, the H-C, L-N, Gonthier, Flores, and the Author model presented in this paper introduce a damping term, leading to the existence of a closed-loop region during the collision penetration and recovery process.

When the restitution coefficient is low, there are significant differences in contact force and penetration among the different contact force models. As the restitution coefficient increases, the differences in penetration and recovery processes among the models decrease, and the region of energy dissipation gradually diminishes.

When $C_r=0.3$ , the force-indentation curves of Gonthier, Flores, and the Author model proposed in this paper are similar, with the Author model being closer to the Gonthier model. When $C_r=0.5$ , the force-indentation curves of Gonthier, Flores, and the Author model largely overlap, while those of the H-C and L-N models are close, indicating similar energy dissipation tendencies. When $C_r=0.8$ , except for the H-C model, the contact forces and indentations of the L-N, Gonthier, Flores, and Author models are very similar, with the penetration curves of the Author and Flores models almost identical. In summary, these results validate the effectiveness of the improved contact force model.

Figure 5 illustrates the relationship between penetration and velocity under different restitution coefficients. Tables 4 to 6 list the actual restitution coefficients and relative errors of different models for three restitution coefficients. It can be found that with the increase of restitution coefficient, the variation trend between penetration and velocity of each model becomes more consistent. Combining with the tables, it is evident that as the restitution coefficient approaches 1, the relative errors in the restitution coefficients of each model also decrease gradually, all of which are less than 10%.

Figure 5.

Relationship between indentation and velocity: (a) C _r  = 0.3, (b) C _r  = 0.5, and (c) C _r  = 0.8.

Table 4.

Actual restitution coefficients of $C_r\prime$ and relative errors of different contact force models at $C_r=0.3$ .

Contact force model	Hertz	Liu	H-C	L-N	Gonthier	Flores	Author
${\stackrel{·}{\delta }}^{(+)}/(\mathrm{m}/\mathrm{s})$	0.5	0.5	0.5	0.5	0.5	0.5	0.5
${\stackrel{·}{\delta }}^{(-)}(\mathrm{m}/\mathrm{s})$	0.5	0.5	−0.2907	−0.3421	−0.1521	−0.1281	−0.1633
$C_r\prime$	—	—	0.5814	0.6842	0.3042	0.2562	0.3266
Relative error/(%)	—	—	93.7333	128.0667	1.4	14.6	8.8667

Table 5.

Actual restitution coefficients and relative errors of different contact force models at $C_r=0.5$ .

Contact force model	Hertz	Liu	H-C	L-N	Gonthier	Flores	Author
${\stackrel{·}{\delta }}^{(+)}(\mathrm{m}/\mathrm{s})$	0.5	0.5	0.5	0.5	0.5	0.5	0.5
${\stackrel{·}{\delta }}^{(-)}(\mathrm{m}/\mathrm{s})$	0.5	0.5	−0.3315	−0.3626	−0.2439	−0.2352	−0.2413
$C_r\prime$	—	—	0.663	0.7252	0.4878	0.4704	0.4826
Relative error/(%).	—	—	32.6	45.04	2.44	5.92	3.48

Table 6.

Actual restitution coefficients and relative errors of different contact force models at $C_r=0.8$ .

Contact force model	Hertz	Liu	H-C	L-N	Gonthier	Flores	Author
${\stackrel{·}{\delta }}^{(+)}(\mathrm{m}/\mathrm{s})$	0.5	0.5	0.5	0.5	0.5	0.5	0.5
${\stackrel{·}{\delta }}^{(-)}(\mathrm{m}/\mathrm{s})$	0.5	0.5	−0.4164	−0.4236	−0.3840	−0.3943	−0.3937
$C_r\prime$	—	—	0.8328	0.8472	0.768	0.7886	0.7874
Relative error/(%).	—	—	4.1	5.9	4	1.425	1.5

When $C_r$ is 0.3, the relative error of the actual restitution coefficient of the Gonthier model is 1.4%, while the error of the improved model proposed in this paper is 8.8667%. The results show that, although the error in the improve d model is less than 10% under low restitution coefficients, there is still a small error. Additionally, the error values of other models are all greater than 10%, with the L-N model having the highest error value at 128.0667%, followed by the H-C and Flores models with errors of 93.7333% and 14.65%, respectively. Hence, except for Gonthier and the improved Author model, the other mentioned contact force models are not suitable for situations with low restitution coefficients.

With $C_r=0.5$ , the relative errors of the actual restitution coefficients of H-C and L-N models are greater than 10%, specifically 32.6% and 45.04%, respectively. The error values of the other models are all less than 10%, with the relative errors of Gonthier, Author, and Flores models being 2.44%, 3.48%, and 5.92%, respectively. Although the Gonthier model has the smallest relative error, it increases as the restitution coefficient increases, while the errors of other models show an opposite trend. The error of the improved model under moderate restitution coefficients has also decreased from 8.8667% at low coefficients restitution to 3.48%. Therefore, it can be inferred that as the restitution coefficient increases to 0.8, the error of the improved model will likely decrease further.

When $C_r=0.8$ , while the relative error of the Gonthier model’s restitution coefficient increases to 4%, the actual relative errors of the restitution coefficient for other models gradually decrease as the restitution coefficient increases, all falling below 10%. Among them, the L-N model has the maximum error of 5.9%, while the minimum error is 1.425% for the Flores model, followed by the Author model with 1.5%.

Figure 6 compares the relative errors of the actual restitution coefficient for different contact force models. It is clear that, excluding the Hertz and Liu models, which do not include the restitution coefficient term, only the Gonthier model and the improved model presented in this paper have errors less than 10% under various conditions with different restitution coefficients. Specifically, when the restitution coefficient is 0.3, the relative error of Gonthier model is 1.4%, while that of the improved model is 8.8667%. However, when the restitution coefficient is 0.8, the relative error of the Gonthier model increases to 4%, while the relative error of the improved model decreases to 1.5%.

Figure 6.

Comparison of actual restitution coefficient errors for different contact force models.

Additionally, the error in the Gonthier model gradually increases with the restitution coefficient, while the error in the improved model and other models decreases as the restitution coefficient increases. Figure 6 shows that the error of the improved model remains below 10%, indicating a wide range of applicability.

To further compare the Gonthier model with the improved model, the initial collision velocity, material properties, and other parameters were set the same as in the previous section, with only the restitution coefficient adjusted to 0.99. The post-collision velocities of the two models were obtained, and the specific results are shown in Table 7. From the analysis, it can be concluded that when the restitution coefficient increases to 0.99, the error of the Gonthier model is 0.3273%, while the improved model demonstrates extremely high accuracy, with an error of only 0.0040%, which is nearly zero. This indicates that the improved model has very high precision under high restitution coefficient conditions.

Table 7.

Actual restitution coefficients and relative errors of Gonthier and Author models at $C_r=0.99$ .

Contact force model	${\stackrel{·}{\delta }}^{(+)}(\mathrm{m}/\mathrm{s})$	${\stackrel{·}{\delta }}^{(-)}(\mathrm{m}/\mathrm{s})$	$C_r\prime$	Relative error (%).
Gonthier	0.5	0.49338	0.98676	0.3273
Author	0.5	0.49498	0.98996	0.0040

In summary, under different restitution coefficients and using the Newton restitution coefficient as the evaluation standard, the proposed improved contact force Author model consistently exhibits an error of less than 10%. Furthermore, as the restitution coefficient increases, the relative error of the improved model gradually decreases, corresponding to an improvement in simulation accuracy.

Application 2

Utilizing the collision experiment depicted in Figure 7 of Minamoto and Kawamura, ⁴⁵ further verification was conducted on the improved contact force model presented in this paper.

Figure 7.

Steel ball collision simulation experimental model.

The experiment involved two steel balls with identical parameters: elastic modulus E ₁  = E ₂  = 210 GPa, Poisson’s ratio σ = 0.3, and yield strength of 518.4 MPa. In the initial state of the experiment, steel ball 2 is suspended vertically, while steel ball 1 is released from a specified swing angle θ ₀ and falls under gravity at a certain initial velocity to collide with steel ball 2. The angle θ ₀ is used to determine the impact velocity of the spheres. The angles θ ₁ and θ ₂ represent the maximum angles of each pendulum after the impact, which are measured from videos recorded by a digital camera at 30 frames per second. ⁴⁵ Ignoring air resistance, the coefficient restitution is then calculated.

Under different initial velocity conditions, Tables 8 to 10 and Figure 8 show the actual restitution coefficients obtained from the experiments, as well as the restitution coefficient values of the Gonthier, Flores, and the improved model in this paper.

Table 8.

Comparison and error data between experimental testing ⁴⁵ and numerical simulations of restitution coefficients from the Gonthier model.

Initial velocity/(m/s)		0.158	0.372	0.763	1.158	1.526	1.895	2.237	2.579	2.895
$C_r$	Experimental data	0.9	0.82	0.754	0.71	0.688	0.666	0.644	0.639	0.622
	Gonthier model data	0.8766	0.7887	0.7218	0.6791	0.6583	0.6377	0.6173	0.6128	0.5972
Relative error/(%).		2.6	3.817	4.27	4.35	4.32	4.25	4.14	4.11	3.98

Table 9.

Comparison and error data between experimental testing ⁴⁵ and numerical simulations of restitution coefficients from the Flores model.

Initial velocity/(m/s)		0.158	0.372	0.763	1.158	1.526	1.895	2.237	2.579	2.895
$C_r$	Experimental data	0.9	0.82	0.754	0.71	0.688	0.666	0.644	0.639	0.622
	Flores model data	0.8937	0.8097	0.7401	0.6937	0.6705	0.6596	0.6240	0.6187	0.6007
Relative error/(%).		0.7	1.26	1.84	2.30	2.54	2.81	3.11	3.18	3.43

Table 10.

Comparison and error data between experimental testing ⁴⁵ and numerical simulations of restitution coefficients from the improved model.

Initial velocity/(m/s)		0.158	0.372	0.763	1.158	1.526	1.895	2.237	2.579	2.895
$C_r$	Experimental data	0.9	0.82	0.754	0.71	0.688	0.666	0.644	0.639	0.622
	Improve model data	0.8962	0.8091	0.7374	0.6902	0.6669	0.6439	0.6212	0.6161	0.5989
Relative error/(%).		0.42	1.33	2.2	2.87	3.07	3.32	3.54	3.58	3.71

Figure 8.

Comparison and error plot between experimental restitution coefficients ⁴⁵ and numerical simulations from the improved model.

The analysis shows that among the three models, the Gonthier model exhibits a large error compared to the experimental data. The Flores model has the smallest error. However, as discussed earlier, the Flores model exhibits errors exceeding 10% under low restitution coefficients, making it unsuitable for conditions with low restitution coefficients. The improved model falls in between the two. As the restitution coefficient decreases, the error of the improved model gradually increases, but the maximum error is only 3.71%. In summary, the improved model has a wide range of applications, and under conditions with high restitution coefficients, the model exhibits very high accuracy.

Conclusion

Based on Hertz theory and classical contact force model, an improved contact force model at joint contact surface is constructed by introducing a nonlinear index factor m. Through steel ball collision experiments, numerous numerical simulations of the improved model were conducted by varying parameters such as collision speed and restitution coefficient. The results were then compared and analyzed with the simulated values of existing contact force models and the experimental test results from the literature, verifying the validity and accuracy of the improved model.

Key findings include: (1) Through numerical simulations with different recovery coefficients and collision velocities, it can be concluded that the relative error of the actual recovery coefficient in the improved model is less than 10%. As the restitution coefficient increases, the error of the improved model gradually decreases. This indicates that the improved model has a broad range of applications and high calculation accuracy. (2) By comparing the results with experimental data, it can be observed that the improved model aligns well with the experimental data at different collision velocities, with a maximum error of only 3.71%, further validating the effectiveness of the improved model.