Abstract

Spherical clearance joints are essential for the successful deployment of space structures. When the clearance is small enough, the contacts will be considered conformal contact, which probably leads to inaccuracies in existing contact force models. To address the limitation, this paper proposes a novel hyperbolic contact surface Winkler model. First, a new fundamental formula incorporating a modified variable exponent is presented. Based on the surrogate modeling method, an optimized surrogate function for the variable coefficient is developed. In the optimization process, the finite element and response surface methods (RSMs) are introduced to improve the precision and reliability of the model. Compared with previous models, this paper organizes a detailed discussion and evaluation to validate the accuracy and application of the new proposed model, after which a dynamic example demonstrates the model's effectiveness. The results highlight the model's accuracy and practical efficacy, showing a strong correlation and minimal margin of error, especially when compared to finite element method (FEM) results. This improvement is attributed to the refined variable exponent, which accurately characterizes the relationship between contact force and penetration depth, and the optimized variable coefficient, which fine-tunes the contact force magnitude. Additionally, the model's versatility extends beyond the geometric properties of the contact bodies, offering a broad application scope. As a foundation of precise impact modeling, it is crucial to address the structural dynamic challenges inherent in high-precision space structures.

1 Introduction

Large deployable structures are extensively utilized in aerospace missions for their lightweight and small size [14]. Numerous studies have focused on the dynamic modeling and nonlinear analysis of these structures, most of which assume the joints are ideal and overlook the impact of joint clearance [58]. Spherical joints are widely used in orbiting space deployable structures. The clearance in these joints is crucial for facilitating smooth deployment and successful assembly of structures in orbit [9]. Due to nonlinear factors such as impact, contact, and friction, joint clearance significantly affects the stability and precision of these structures [1014]. Therefore, accurately evaluating contact forces within spherical clearance joints is of great significance for the modeling and dynamic analysis of space deployable structures.

Many researchers have explored the influence of spherical clearance joints on the dynamics of mechanical systems and issues related to structural dynamics [1521]. Tian et al. [22] examined the effects of clearances and lubrication within the system's spherical joints using a contact law based on Hertzian principles. Erkaya [23] studied the influence of spherical clearance joints and flexible joints on the mechanism's vibration. He observed a chaotic phenomenon in the system dynamics, primarily resulting from spherical clearance joints. Research [2426] has introduced a methodology for predicting the dynamic characteristics of spatial parallel mechanisms that incorporate multiple spherical clearance joints. Wang and Wang [27] explored the influence of the characteristics of multiple flexible actuated rods on the dynamic responses of a parallel mechanism that includes a spherical clearance joint. Jing and Liu [28] examined the effects of various surface coatings, the absence of coating, and different clearance levels on the dynamic behavior of a mechanism with multiple spherical joint clearances. Ma and Li [29] analyzed the dynamic response of planar structures, considering both single and multiple revolute clearance joints under impact loads. However, most studies have relied on previous contact force models with specific limitations or inaccuracies in calculating contact forces.

The interaction within a spherical clearance joint is a typical example of an axisymmetric contact issue. Hertz [30] is the first to examine the stress distribution characteristics in the region where two elastic bodies come into contact. According to Hertz's theory, the contact region between these elastic bodies ideally takes an elliptical shape, while the axisymmetric contact issue between spheres can be regarded as an unusual case with a circular contact region. Based on Hertz's theory, Johnson [30] established a clear and quantifiable expression of the relationship between contact force and penetration depth. Lankarani and Nikravesh [31] introduced a simplified formula for contact force using an explicit indentation expression. Hunt and Crossley [32] refined the exponent of indentation within the contact force function through optimization. Goodman and Keer [33] proposed an approximate solution to the contact stress problem involving an elastic sphere penetrating an elastic cavity, assuming a confined contact region. Steuermann [34] further explored the pressure distribution within an elastic body exhibiting an exponential polynomial configuration. However, this approach requires prior knowledge of the elastic body's shape, making it not universally applicable for multiple contact problems. Liu et al. [35] developed an approximate contact model for spherical joints with clearances, using the distributed elastic forces to model the compliance of surfaces in contact. Fang et al. [36] proposed a theoretical formulation to estimate the contact pressure distribution for conformal contact in spherical plain bearings. The Winkler elastic foundation model is often used to evaluate contact forces [30,3739]. The Winkler model treats the contact between two elastic bodies as a rigid body with an equivalent curvature pressed against a series of longitudinal springs. The contact force is the combined elastic forces caused by independent spring deformation. However, the difficulty in determining the foundation's elastic modulus results in the Winkler model being unable to calculate contact forces accurately. It means that the contact forces can accurately calculate if the foundation's properties can be appropriately selected. Therefore, we decide to develop a new contact force model with excellent adaptability in spherical clearance joints on the basis of the Winkler model.

The finite element method (FEM) serves as a highly effective numerical approach for resolving contact forces and has been widely applied in the development of novel contact models [4043] and the analysis of complex contact-related issues [4446]. Jia and Chen [47] conceived a novel conformal contact force model for spherical clearance joints by integrating principles of elastic foundation theory, FEM, and least squares methods. Wang et al. [48] introduced an expanded Hertz model to predict the mechanical responses of an incompressible Mooney–Rivlin half-space subjected to finite spherical indentation, validating the model's constitutive parameters through finite element simulations. In a separate study, Ding et al. [49] proposed an approach for determining the normal contact forces of viscoelastic particles, using FE simulations to validate the model and define viscoelastic material properties. Yaylaci et al. [50,51] explored the frictionless contact issue in a functionally graded layer, utilizing analytical techniques, the finite element method, and a multilayer perceptron. Lan et al. [52] advanced an elastic–plastic contact analysis model for irregular ellipsoid surfaces based on the principles of elastic–plastic fractal theory and investigated fluctuations in contact loads associated with elastic ellipsoid bodies. Despite the widespread use of the FEM in examining contact-related challenges, it remains limited in addressing the complexities of mechanism dynamics and structural dynamics issues featuring clearance joints, as its computational efficiency may not meet the practical requirements of engineering practice [53]. Therefore, precise contact force models are indispensable in engineering practice for such challenges. In this paper, the FEM emerges as a crucial tool, contributing a lot to the accuracy and efficiency of contact force computations.

This paper proposes a hyperbolic contact surface Winkler model for calculating contact forces in spherical clearance joints. The model features a novel formula incorporating a modified variable exponent based on the circumferential Winkler model. The variable coefficient is fine-tuned by both the FEM and the response surface method (RSM). The remainder of this paper is organized as follows: Sec. 2 presents the physical model characterizing spherical clearance joints and reviews previous contact force models, highlighting their advantages and disadvantages. Section 3 introduces a novel formula for a dimensionless contact force model, including an optimization surrogate function applied to the variable coefficient, which is correlated with the profiles of the contacting bodies. Section 4 extensively compares the proposed model to previous models, validating its precision. A detailed analysis evaluates its applicability, interrelation with other variables, and margin of error. Finally, Sec. 5 concludes the study, offering insights into potential future research directions.

2 Problem Description

Figure 1 shows the physical model of a spherical clearance joint. Joint clearance ΔR is defined as
(1)

where R2 and R1 are the socket's and ball's radii, respectively.

Fig. 1
Physical model of a spherical clearance joint
Fig. 1
Physical model of a spherical clearance joint
Close modal
When the socket and ball come into contact, localized deformation occurs, generating a contact force FN. This analysis model is depicted in Fig. 2, showcasing the extent of penetration depth δ and the resulting contact force FN. The penetration depth δ can be obtained by
(2)
where e denotes the eccentricity, determined by the relative position of the socket's center in relation to the ball's center. The eccentricity e can be expressed as
(3)
Fig. 2
Penetration depth and contact force
Fig. 2
Penetration depth and contact force
Close modal
The essence of contact problems lies in developing a contact force model capable of numerically representing the relationship between the extent of penetration depth δ and the resultant contact force P. According to Hertz's theory, when the socket and the ball are in contact, the contact region forms a circle with a radius a, as illustrated in Fig. 3. Within this contact region, the pressure is distributed according to the function p(r) as
(4)
where p0 represents the maximum pressure at the center of the contact region, and r denotes the distance between the contact point and the center of the contact area. Hertz assumed that the approach distance between two objects equals the penetration depth δ, and the relationship satisfies
(5)
where E* is the composite modulus related to Young's modulus E and Poisson's ratio ν of the contact bodies' material properties: 1/E*=(1ν12)/E1+(1ν22)/E2. R is the equivalent radius of curvature on the contact region: 1/R=1/R1+1/R2. The total load P on two contact bodies equals the contact force between the two objects and can be expressed by the penetration depth δ as
(6)
Fig. 3
Pressure distribution of Hertz theory
Fig. 3
Pressure distribution of Hertz theory
Close modal
Steuermann [34] assumed that the contact bodies' profile is smooth and continuous, and the pressure distribution pn(r) having the form of Anr2n, n is the exponent of the polynomial to describe the contact bodies' profile. The force–displacement relationship for the spherical clearance joint can be expressed as Eq. (7). It should be noted that in Steuermann's theory, the case n=1 corresponds to Hertz's theory
(7)
Liu [35] proposed an approximate contact force model for the spherical joint with clearance. This contact model is founded upon the elementary Winkler elastic foundation model and integrates the pressure distribution principles outlined in Hertz's theory, incorporating the semi-angle ε of contact defined in the Steuermann theory. Liu induced a modifying coefficient, denoted as k, to preserve geometrical congruence. The determination of this coefficient was achieved through analytical outcomes from the FEM. According to the FEM results, Liu's model can be expressed as
(8)
Previous contact force models offer practical approximate approaches for calculating contact forces, many of which rely on certain assumptions, leading to significant discrepancies between numerical and experimental contact force calculations. For instance, Hertz's model assumes that the dimension of the contact region is far smaller than the contact body's size and the relative curvature radius of the contact surface. This assumption suggests nonconformal contact [54], making Hertz's theory suitable for scenarios with substantial clearances and minimal contact forces [30,5557]. In Steuermann's theory, an unreasonable exponent n representing the profiles in the polynomial will result in considerable errors in calculating contact force, which is also presented in Liu's model.

Furthermore, for deployable structures in orbit, the clearance of spherical joints tends to be minimal. In this condition, the conformal contact between the socket and the ball leads to significant numerical disparities when applying conventional contact force models. Therefore, it is urgent to develop a more accurate contact force model that is specifically tailored for joints with small clearance. The newly developed model enables precise evaluation of the contact force between the socket and the ball.

3 Model Establishment

3.1 Basic Formulation.

The Winkler elastic foundation model offers a simplified approach to calculating the contact force between contact bodies of any profile [30]. Figure 4 illustrates a basic Winkler elastic foundation utilized to simulate the compressed contact body, with a height of h and elastic modulus of K. Another contact body is regarded as a rigid indenter with the profile function z(x,y). As shown in Fig. 4, longitudinal deformation occurs due to compression, forming the contact region when the rigid indenter compresses the elastic foundation downward. Within the contact region, the deformation of the elastic foundation is nonzero. According to the definition, the elastic foundation will generate a longitudinal pressure distribution p(x,y) on the rigid indenter. The rigid indenter's profile function z(x,y) is taken as the sum of the profile functions of the two contact bodies, represented as
(9)
Fig. 4
Winkler elastic foundation model
Fig. 4
Winkler elastic foundation model
Close modal
Assuming no interaction among the springs within the elastic foundation, the distributed force in the elastic foundation is solely associated with the normal elastic displacement of the foundation, denoted as u¯z(x,y), which can be represented as
(10)
where δ is the center displacement of the indenter, which can be regarded as the penetration depth. The pressure distribution p(x,y) can be expressed as
(11)
where p(x,y) is the pressure distribution function, representing the magnitude of pressure per unit area, with units of (N·m−2). K is the elastic modulus of the foundation, with units of (N·m−2). h and u¯z(x,y), respectively, represent the height and deformation of the foundation, with units of (m).
For two spherical contact bodies with curvature radii R1 and R2, the maximum and minimum main curvature radii of each contact body's sections are the same, denoted as R1=R1=R1 and R2=R2=R2, respectively. In the contact region, the longitudinal deformation of the elastic foundation can be expressed as
(12)

where R and R denote the maximum and minimum relative main curvature radii and satisfy 1/R=1/R1+1/R2 and 1/R=1/R1+1/R2.

The pressure distribution p(x,y) in the contact region can be expressed as
(13)
According to Hertz's theory, the contact region is an ellipse with half axes a=(2δR)1/2 and b=(2δR)1/2, respectively. For the axisymmetric case, there is a=b=(2δR)1/2. Equation (13) becomes
(14)
Subsequently, the total load P of the indenter can be calculated and expressed as
(15)

The reasonable selection of elastic modulus and foundation height makes it easy to calculate the accurate contact forces. However, we still face the challenge of adopting the Winkler model, which has limitations when utilizing it. As depicted in Fig. 4, the total load equals the integral of the vertically distributed force in the contact area. Initially, the profiles of the two contact bodies expand along with the curvature radius of the lower contact body, with the vertically distributed force provided by the elastic foundation pointing toward the center of the socket.

When the curvature radii of the contact bodies vary significantly, and the penetration depth is small, the contact arc on the contact surface is also small, allowing for a reasonable approximation of the contact force within the Winkler model. However, as the penetration depth increases, the contact region proportionally expands, leading to a larger contact arc and its associated central angle, as shown in Fig. 5. It becomes essential to rigorously distinguish the distribution pressure directed toward the center of the contact body and the distribution pressure in the vertical direction, especially when the curvature radii of the two contact bodies differ only slightly.

Fig. 5
Contact area in the socket
Fig. 5
Contact area in the socket
Close modal

This paper introduces a hyperbolic contact surface Winkler model, an extension of the classical Winkler elastic foundation model. As depicted in Fig. 6, the elastic foundation transfers from its original flat state to a state characterized by a radius of curvature. The pressure field, distributed vertically, plays a significant role in determining the contact force between the socket and the ball, which contrasts with the force in the traditional Winkler model that is directed toward the socket's center and perpendicular to the elastic foundation. In this study, it is assumed that the radius of the contact area in the classical Winkler model equals the semi-arc length of the contact area as defined in the circumferential Winkler model.

Fig. 6
Distributed pressure in the circumferential Winkler model
Fig. 6
Distributed pressure in the circumferential Winkler model
Close modal
The pressure distribution of Eq. (13) can be translated to Eq. (15) based on the geometric relationships: a=R2α,x=R2θ
(16)
The total load P can be calculated by the following Eq. (16) as
(17)
where D denotes the contact region: D={0θα,0φ2π}, dσ denotes the area infinitesimal, and it can be obtained using Eq. (17)
(18)

where R2 is the radius of the socket, φ and θ, respectively, denote the azimuthal angle and the polar angle of the area infinitesimal dσ, and dφ and dθ denote the infinitesimal increment of the azimuthal angle and the polar angle of the area infinitesimal dσ.

Substituting Eq. (18) into Eq. (17), Eq. (17) becomes
(19)
where α is the semicentral angle of the contact region in the longitudinal direction. When the clearance of the spherical joint ΔR=R2R1 is small, the socket's and the ball's radius are close: R1R2. Hence, the total load P in the vertical direction comes to
(20)
According to the estimation of penetration depth δ by the Winkler model, the semicentral angle of the contact area α can be expressed as
(21)

In Eq. (19), the total load P is (N). The equation's right side comprises three terms. The first term encapsulates the elastic foundation's characteristics K and h, with the parameters selected based on the specific contact problem. The second term signifies the geometric attributes of the contact bodies. The third term explicitly represents the penetration depth δ, which is a crucial element in the contact force model's formulation.

It is assumed that the elastic foundation modulus K is equivalent to the composite modulus E* of the contacting bodies. Following this, we define the dimensionless load P* and the dimensionless penetration depth δ* accordingly
(22)
Equation (20) can be transformed into a dimensionless expression as
(23)

To establish an accurate contact force model, it is also necessary to appropriately define the elastic foundation's height h. This is a challenge that we often encountered when applying the traditional Winkler model. According to Eq. (23), both the socket radius and the height of the elastic foundation are parameters associated with the geometric properties of the contact bodies. Consequently, the height of the elastic foundation h should also be determined based on the profiles of the contact bodies.

Furthermore, the numerical value of the dimensionless contact force is equivalent to the dimensionless load P*. In this paper, Eq. (23) is transformed into the following form:
(24)

where f2 is a function of dimensionless penetration depth δ*, which provides a basic formula for the contact force model, including a specific variable exponent. f1 is a function to be optimized about the geometric properties of the contact bodies, which is assumed to be related to the dimensionless radius of the socket R2*=R2/ΔR. f1 provides an appropriate variable coefficient for the contact force model.

The key of the proposed model is to obtain a suitable functional relationship to describe the effect of dimensionless socket radius R2* on the variable coefficient.

3.2 Parameter Determination.

In this section, the FEM and the RSM are utilized to derive the surrogate function for the variable coefficient f1. Figure 7 outlines the step-by-step process of constructing a surrogate function for the variable coefficient f1. The process begins by investigating the relationship between the geometric properties of the contact bodies and the variable coefficient within the contact model. Among these geometric properties, the dimensionless socket radius R2* is the primary descriptor of the contact bodies' profiles, aside from the joint clearance ΔR.

Fig. 7
Solution procedure of the surrogate function of the variable coefficient
Fig. 7
Solution procedure of the surrogate function of the variable coefficient
Close modal

In the design of experiments, a series of simulations is carefully planned and executed by varying the value of the dimensionless socket radius R2*. Each corresponding variable coefficient, obtained from the FEM analysis, is optimized by a least square fitting algorithm, after which a surrogate function that includes the prespecified parameters is developed. This function is further refined by applying the least square fitting algorithm. Finally, the effectiveness of the surrogate function for the variable coefficient is evaluated, leading to the development of a surrogate model for the dimensionless contact force P*.

Figure 8 presents the FEM model for a spherical clearance joint. As depicted in Fig. 8(a), a ball with a radius R1 is situated within a socket characterized by a radius R2. Given that a sphere is a typical rotating body, an axisymmetric modeling method is utilized to model the spherical clearance joint, as shown in Fig. 8(b).

Fig. 8
Modeling of the FEM analysis: (a) geometric parameters and (b) FEM model
Fig. 8
Modeling of the FEM analysis: (a) geometric parameters and (b) FEM model
Close modal

A static analysis is performed to measure the contact force inherent to the spherical clearance joint to minimize the impact of inertial force on the calculation results. abaqus/standard is employed as the solver for the static analysis. Here, a displacement u=ΔR+δ is intentionally imparted at the center of the spherical ball. The entire external surface of the socket, labeled as AB, is completely fixed. In the setting of contact behavior in the “interaction” option, the normal contact behavior between the ball and socket is defined as “hard contact,” while the tangential behavior is set as “frictionless.” The structure is meshed using the linear standard triangle element CAX3, a three-node linear axisymmetric triangle stress element. The mesh is densified around the contact center to ensure accuracy.

In space deployable structures, the joint clearance ΔR typically ranges from tens to hundreds of micrometers, while the joint size varies from a few to tens of millimeters. Interestingly, the issues of “large penetration depth under large clearance” and “small penetration depth under small clearance” are physically identical. For the simulation, there are 15 cases when we set the joint clearance ΔR to 0.1 mm and the dimensionless socket radius R2* to values between 20 and 300. Based on the FEM results and research on clearance joints [15,35,39,42], the maximum penetration depth δm* for each case is set to 0.5, since this range is sufficient to address most engineering problems. The Young's modulus and Poisson's ratio values are 2×1011Pa and 0.3, respectively.

According to the FEM results, the surrogate function of the variable coefficient f1 can be expressed as
(25)

Figure 9 provides a comparative analysis between the surrogate function of the variable coefficient and the optimal fitting derived from the FEM results. The correlation between the proposed model and the FEM outcomes is evaluated and presented in Table 1. Key metrics in this evaluation include the multicorrelation coefficient R2, the adjusted multicorrelation coefficient Radj2, the root-mean-squared error RMSE, and the relative error RE.

Fig. 9
Comparison of the surrogate function of the variable coefficient and the best fitting of the FEM
Fig. 9
Comparison of the surrogate function of the variable coefficient and the best fitting of the FEM
Close modal
Table 1

Correlation evaluation of the surrogate function and FEM results

Evaluation coefficientsValue
R21.0000
Radj21.0000
RMSE0.0396
RE−0.014% to 0.022%
Evaluation coefficientsValue
R21.0000
Radj21.0000
RMSE0.0396
RE−0.014% to 0.022%

The higher the values of the multicorrelation coefficient R2 and the adjusted multicorrelation coefficient Radj2 (closer to 1) are, the stronger the alignment between the proposed model and the FEM results is. On the other hand, lower values of RMSE and RE indicate a lower level of prediction error. The results presented in Table 1 suggest that the surrogate function for the variable coefficient aligns closely with the optimal fitting obtained through FEM. This finding confirms the near-optimal alignment between the predictions of the surrogate function for the variable coefficient and the FEM results.

Substituting Eq. (25) into Eq. (24), the dimensionless contact force P* can be expressed as
(26)
The proposed contact force model can be written as
(27)

4 Model Analysis

4.1 Model Comparison.

In this section, we conduct a comparative analysis between the proposed model and various contact force models, including the Hertz, Steuermann, Liu model, classical Winkler model, and the FEM. Through comparisons from different perspectives, we get a comprehensive understanding of performance and characteristics. In this comparative study, we set the clearance size at 0.1 mm and standardized the dimensionless socket radius at 200. The contacting bodies' material properties are consistent with those outlined in Sec. 3.2 to ensure a fair and accurate comparison among all models.

4.1.1 Comparison of the Relationship Between the Dimensionless Contact Force and Penetration Depth.

Figure 10 illustrates the relationships between the dimensionless contact force and the dimensionless penetration depth, examined across various contact force models or methodologies. Figure 10(a) presents a comprehensive view of the relationship curve, including the maximum dimensionless penetration depth within a range of 0.5. Additionally, Figs. 10(b)10(d) display the relationship curves within specified intervals of the dimensionless penetration depth.

Fig. 10
Comparison of the relationship between the dimensionless contact force and dimensionless penetration depth obtained from different contact force models: (a) maximum dimensionless penetration depth range, 0<δ*<0.5, (b) local details, 0<δ*<0.1, (c) local details, 0<δ*<0.2, and (d) local details, 0.35<δ*<0.5
Fig. 10
Comparison of the relationship between the dimensionless contact force and dimensionless penetration depth obtained from different contact force models: (a) maximum dimensionless penetration depth range, 0<δ*<0.5, (b) local details, 0<δ*<0.1, (c) local details, 0<δ*<0.2, and (d) local details, 0.35<δ*<0.5
Close modal

Figure 10(a) shows that the proposed model's correlation between dimensionless contact force and dimensionless penetration depth aligns more closely with FEM results. Notably, the curve trend generated by the proposed model closely matches the trend of FEM results in contrast to other models. Figure 10(b) indicates that at the beginning of contact, the results of the Hertz and Steuermann models closely approximate that of FEM at the beginning. However, it becomes apparent that only when the dimensionless penetration depth exceeds approximately 0.13 does the error in the dimensionless contact force obtained from the proposed model become smaller than that of previous models (Fig. 10(b)) across the entire range of dimensionless penetration depth. This discrepancy is attributed to the superior variable component, which accurately captures the correlation between dimensionless contact force and dimensionless penetration depth. Conversely, a notable difference exists between the variable exponent in the Hertz and Steuermann models and the variable exponent in the FEM, resulting in a more precise approximation of the dimensionless contact force only within a limited range of dimensionless penetration depth (Fig. 10(c)). With the increase of dimensionless penetration depth, the error in dimensionless contact force significantly magnifies, ultimately leading to inadequate reference values for the results.

In Fig. 10(a), the curve trend observed in the Liu model closely resembles that of FEM; however, there is a notable disparity in the magnitude of the dimensionless contact force compared to the FEM results. This discrepancy in the results of dimensionless contact forces can be directly attributed to the trend of the curve, including both the amplitude of dimensionless contact force, which is inherently linked to variable exponent and coefficient. Reasonable selection of the variable coefficient can effectively adjust the magnitude of dimensionless contact force to ensure a better alignment with FEM results. Moreover, an optimized Winkler model was proposed for comparison, wherein its variable exponent is derived from the classical Winkler model, and its variable coefficient is optimized using the method adopted in this paper. The expression of the dimensionless contact force of the optimized Winkler model is shown as
(28)

As depicted in Fig. 10(a), the hyperbolic contact surface Winkler model optimizes the variable exponent based on the classical Winkler model. Throughout the entire range of dimensionless penetration depth, the error of the classical Winkler model exceeds that of the model proposed in this paper, which exactly proves the rationality and applicability of the proposed hyperbolic contact surface Winkler model in this study.

4.1.2 Comparison of the Correlation.

This paper not only compares the relationship curves between dimensionless contact force and dimensionless penetration depth but also employs correlation evaluation as a quantitative method to assess the proposed model. The results of this evaluation are presented in Table 2, with a visual interpretation provided in Fig. 11. Both Table 2 and Fig. 11 display the absolute values of the multicorrelation coefficient and the adjusted multicorrelation coefficient. It means that the correlation within the proposed model is improved when these coefficients approach a value of 1. Also, the magnitude of the root-mean-square error (RMSE) indicates the model's predictive accuracy.

Fig. 11
Comparison of the correlation evaluation among different contact force models
Fig. 11
Comparison of the correlation evaluation among different contact force models
Close modal
Table 2

Correlation evaluation of different contact force models

ModelR2Radj2RMSE
Hertz7.28397.744153.2215
Steuermann5.31795.668948.3240
Liu2.34292.528632.3994
Winkler0.91930.91487.7896
Proposed0.96200.95995.2469
ModelR2Radj2RMSE
Hertz7.28397.744153.2215
Steuermann5.31795.668948.3240
Liu2.34292.528632.3994
Winkler0.91930.91487.7896
Proposed0.96200.95995.2469

Table 2 and Fig. 11 reveal that the Hertz and Steuermann models display multicorrelation coefficients and adjusted multicorrelation coefficients significantly surpassing 1. Concurrently, the RMSE values are notably high. As depicted in Fig. 10, there is a variable exponent akin to FEM in Liu's model that enhances its correlation with FEM beyond the first two models. However, the suboptimal variable coefficient results in a larger overall error for the dimensionless contact force, as shown in Fig. 11. The classical Winkler model, with an optimized variable coefficient, exhibits a better correlation with FEM than previous contact force models, and its RMSE is considerably lower. The model proposed in this paper markedly ameliorates the correlation and overall error compared with the classical Winkler model. The results indicate that this paper's circumferential transformation approach in this paper can effectively refine the classic Winkler model and derive a more reasonable variable exponent. The dimensionless contact force can be computed more precisely when combined with the optimized variable coefficient.

The results demonstrate the superior congruence between the proposed contact force model and the FEM. This superiority is manifested in the variable exponent, the variable coefficient, or both. There is an enhanced variable exponent in the hyperbolic contact surface Winkler model proposed in this paper, which suggests that the circumferential transformation approach of the classical Winkler model is both practical and reasonable. Moreover, the results from the proposed model exhibit fewer discrepancies than other models, which can be attributed to the optimized variable coefficient. The variable coefficient, which is sensitive to the profiles of the contact bodies, adjusts the contact force accordingly. Therefore, it can be concluded that the hyperbolic contact surface Winkler model proposed in this paper offers a more accurate and reasonable method for modeling contact force.

4.2 Parameter Analysis.

This section delves into the effects of two crucial geometric properties, penetration depth and socket radius, on the proposed contact force model. A comprehensive investigation is conducted to evaluate the feasibility and versatility of the proposed contact force model.

4.2.1 Application Range of Penetration Depth.

For small-clearance spherical joints, the penetration depth increases rapidly during the initial contact stage, leading to a swift expansion of the contact region radius and a proportional enlargement of the central angle of the contact arc. However, as the penetration depth continues to escalate, the expansion rate for both the contact region and the contact arc begins to decrease, which makes Eq. (21) inappropriate for this scenario. Therefore, 0.5 is set as the reference range of the dimensionless penetration depth. This study extends the maximum dimensionless penetration depth to 1 and investigates its influence on the proposed contact force model.

Figure 12 visually depicts the comparison between the dimensionless contact force results obtained from the proposed model and those derived from the FEM. We can get the conclusion from Fig. 12 that the results of the proposed model are in satisfactory agreement with the FEM results, considering that the surrogated function was parameterized with a different penetration depth range.

Fig. 12
Comparison of the results obtained by the proposed model and FEM: (a) relationship of the dimensionless contact force and dimensionless penetration depth and (b) relative error to the dimensionless penetration depth
Fig. 12
Comparison of the results obtained by the proposed model and FEM: (a) relationship of the dimensionless contact force and dimensionless penetration depth and (b) relative error to the dimensionless penetration depth
Close modal

As the dimensionless penetration depth progresses beyond this range, the relative error between these two approaches gradually increases, as shown in Fig. 12(b). Based on the analysis of various contact force models presented in Sec. 4.1, it can be reasonably inferred that as the dimensionless penetration depth continues to increase, the discrepancy between the dimensionless contact force calculated by the proposed model and that derived from the FEM will also increase. However, it will still be smaller than other previous models, as demonstrated in Fig. 10(a), even though it is not presented.

This can be attributed to the variable exponent of the proposed contact force model. Despite the fact that the proposed model is designed for a specific range of dimensionless penetration depth, the similarity in variable exponent ensures that the model exhibits a trend of closely mirroring the FEM results and maintains a relatively low level of relative error.

According to Eq. (24), the contact force model proposed in this study adopts a dimensionless expression. This formulation simplifies the structure, reduces the relevant geometric parameters from 3 to 2, and also establishes a close correlation between the application range of the penetration depth and the clearance size within the proposed contact force model. Figure 13 illustrates the reasonable application range of penetration depth for the proposed contact force model, as well as providing evaluation results of the maximum relative error.

Fig. 13
Application range of penetration depth
Fig. 13
Application range of penetration depth
Close modal

Based on the above results, considering the geometric characteristics of the contact bodies in real-world engineering scenarios, the proposed contact force model maintains a reasonable application range for the penetration depth. Moreover, it exhibits greater accuracy compared with previous models.

4.2.2 Application of Socket Radius.

The socket radius is a crucial geometric parameter in contact problems. This study employs a dimensionless expression to establish the contact force model, substituting the socket radius with the dimensionless socket radius. This method simplifies both the contact force model and the contact problem. Theoretically, no matter how the clearance size changes, the corresponding contact problem remains the same when the dimensionless socket radius and dimensionless penetration depth remain constant. In this paper, we explore the applicability range of the proposed model, taking into account variations in the dimensionless socket radius.

Figure 14 illustrates the correlation assessment of the proposed model with the FEM across varying dimensionless socket radius. It is clear that within the defined range of the dimensionless socket radius, the proposed model maintains a significant level of correlation. Both the multicorrelation coefficient and the adjusted multicorrelation coefficient consistently exceed 0.96. As the dimensionless socket radius extends, there is a noticeable decrease in correlation, but it still surpasses 0.95 and tends to stabilize. The relative root-mean-squared error RMSE/P¯* is the metric for assessing the model's error. As depicted in Fig. 14, as the dimensionless socket radius increases, the relative root-mean-squared error exhibits a corresponding increase, converging to a specific value under 6.8%. The proposed model demonstrates its accuracy and effectiveness as the dimensionless socket radius increases.

Fig. 14
Correlation evaluation between the proposed model and FEM with different dimensionless socket radius
Fig. 14
Correlation evaluation between the proposed model and FEM with different dimensionless socket radius
Close modal

The dimensionless nature of the proposed model requires a nonzero joint clearance, which implies that the model is unsuitable for entirely conformal contact problems. However, the problem can be considered a quasi-conformal contact problem when the socket radius significantly exceeds the joint clearance. Figure 15 shows the feasible application range of the socket radius in the proposed model, from which we can observe that the proposed model provides an excellent fit and minor error, regardless of the magnitude of the socket radius. These results underscore the wide applicability of the proposed model.

Fig. 15
Application range of socket radius
Fig. 15
Application range of socket radius
Close modal

4.3 Numerical Example

4.3.1 Model Configuration.

The crank–slider mechanism is a prevalent mechanical transmission device that converts rotary motion into reciprocating linear motion. As depicted in Fig. 16, the crank–slider mechanism comprises a crank, a coupler rod, and a slider. The crank is a rod affixed to the shaft, with one end connected to the coupler rod and the other linked to the slider. The crankshaft is attached to the fixed shaft and the connecting rod, and the coupler rod is connected to the slider via a joint. The crank–slider mechanism is frequently employed to validate the accuracy and effectiveness of the impact contact force model and investigate the influence of joint clearance on the dynamic characteristics of the mechanism [17,18].

Fig. 16
Crank–slider mechanism with clearance joint
Fig. 16
Crank–slider mechanism with clearance joint
Close modal

The geometric parameters and material properties of the crank–slider mechanism are presented in Table 3, while the dynamic simulation parameters of the system are listed in Table 4. At the onset of the dynamic simulation, both the crank and the connecting rod are maintained horizontally. A constant moment of force is applied at the fixed axis of the crank to drive the mechanism.

Table 3

Geometric parameters and material characteristics of the crank–slider mechanism

MemberLength (m)Mass (kg)
Crank0.050.3
Coupler0.120.21
Slider0.14
Clearance0.00005
MemberLength (m)Mass (kg)
Crank0.050.3
Coupler0.120.21
Slider0.14
Clearance0.00005
Table 4

Simulation parameters of the dynamic analysis

ParameterValue
Young's modulus E (GPa)207
Poisson ratio υ0.3
Coefficient of restitution cr0.5
Sliding friction coefficient μk0.05
Tolerance velocity v0 (m/s)0.0001
Tolerance velocity v1 (m/s)0.01
Driving moment (N·m)0.5
ParameterValue
Young's modulus E (GPa)207
Poisson ratio υ0.3
Coefficient of restitution cr0.5
Sliding friction coefficient μk0.05
Tolerance velocity v0 (m/s)0.0001
Tolerance velocity v1 (m/s)0.01
Driving moment (N·m)0.5

4.3.2 Model Establishment

4.3.2.1 Finite particle method.
The finite particle method (FPM) is a method of structural dynamics analysis that originates from the vector intrinsic finite element method. During modeling, the continuous structure is discretized into finite particles that move under internal and external forces. The relative position of discrete particles represents the structure's shape, and the motion trajectory of discrete particles depicts the changes of the structure's shape. Each discrete particle's motion adheres to Newton's laws of motion. Therefore, for a given particle J, its control equation can be expressed as follows:
(29)

where mJ and IJ represent the equivalent mass and equivalent moment of inertia of discrete particles, respectively. FJint and MJint are the reaction forces on the particle caused by structural deformation, respectively. Their detailed derivation process can be referred to Ref. [58].

As described by Eq. (29), the motion equation of a discrete particle is an ordinary differential equation, which allows the use of the explicit integration method, specifically the central difference method, for computation. In the application of contact mechanics models, explicit algorithm indeed has their advantages in that they simplify the contact detection process and allow for real-time updates of contact forces, which makes them highly efficient for such applications. This approach increases the accuracy and effectiveness of the model, particularly in dynamic simulations involving contact forces. Among them, the contact force can be expressed as
(30)

where eN and eT denote the unit vector of the normal force and tangential force during the contact process, and FN and FT denote the normal contact force and tangential friction force.

4.3.2.2 Normal contact force.
The impact behavior of clearance joints is a dynamic process that necessitates the consideration of energy dissipation due to material damping. Lankarani and Nikravesh introduced the concept of the recovery coefficient and derived a specific expression for the damping factor [59]. This study integrates the static contact force model proposed in Sec. 3.2 with the Lankarani–Nikravesh model to derive a hybrid impact force model
(31)

where P denotes the elastic contact force, which can be obtained by Eq. (27). cr is the coefficient of restitution, δ˙ denotes the instantaneous normal impact velocity, and δ˙() is the initial normal impact velocity.

4.3.2.3 Tangential friction force.
The tangential force is mainly friction force. In this paper, a modified Coulomb's friction model suggested by Ambrósio [60] is employed and can be expressed as
(32)
where μk is the sliding friction coefficient, FN is the normal contact force, vT is the relative tangential velocity, and cd is the dynamic correction coefficient given as
(33)

where v0 and v1 are the tolerance velocities of the friction model.

4.3.3 Results and Discussion.

In order to validate the accuracy of the modeling method and showcase the effectiveness of the proposed contact force model, this study compares the dynamic response of the mechanism under both ideal joint and clearance joint conditions. This comparison is executed by the FEM and FPM, respectively. The dynamic response results are represented in Figs. 1721.

Fig. 17
Comparisons of the slider's displacement: (a) ideal joint and (b) clearance joint
Fig. 17
Comparisons of the slider's displacement: (a) ideal joint and (b) clearance joint
Close modal
Fig. 18
Comparisons of the slider's velocity: (a) ideal joint and (b) clearance joint
Fig. 18
Comparisons of the slider's velocity: (a) ideal joint and (b) clearance joint
Close modal
Fig. 19
Comparisons of the slider's acceleration: (a) ideal joint and (b) clearance joint
Fig. 19
Comparisons of the slider's acceleration: (a) ideal joint and (b) clearance joint
Close modal
Fig. 20
Comparisons of the clearance joint's contact force: (a) ideal joint and (b) clearance joint
Fig. 20
Comparisons of the clearance joint's contact force: (a) ideal joint and (b) clearance joint
Close modal
Fig. 21
Center's trajectory of clearance joint
Fig. 21
Center's trajectory of clearance joint
Close modal

Figures 1719, respectively, depict the slider's displacement, velocity, and acceleration responses. The dynamic response results of the slider reveal that the outcomes of the FEM and the FPM are mainly consistent for the ideal joint situation, thereby the accuracy of the FPM is validated. The slider's displacement, velocity, and acceleration responses are slightly advanced when considering the clearance joint. Particularly in terms of the slider's acceleration response, compared with the ideal joint, the acceleration of the slider exhibits significant and high-frequency impact when considering the joint clearance. These results illustrate the effectiveness of the dynamic modeling method for structures with clearances adopted in this study.

Figure 20 presents the contact force curve of the clearance joint. In the case of an ideal joint, the internal interaction forces are relatively smooth, while if a clearance exists, the contact force inside the joint exhibits a large amplitude in the early stage of motion, accompanied by significant fluctuations. This suggests the presence of high-frequency impact behavior at the joint during the initial phase of the mechanism's movement. This impact behavior persists throughout most of the mechanism's movement, resulting in a larger amplitude of contact force and an intense oscillation of an extended period. Furthermore, when a clearance is present in the joint, the amplitude of the contact force inside the joint is significantly larger than that in an ideal joint.

Figure 21 illustrates the relative motion trajectory within the clearance joint. The trajectory effectively portrays the relative motion within the joint during the initial stage of the mechanism's motion, with the joint center moving downwards relative to the shaft sleeve, accompanied by noticeable impact behavior, which is corroborated by the contact force response results depicted in Fig. 21. Throughout the mechanism's movement, the clearance joints remain in continuous contact for an extended period, in which the circumferential contact trajectories generate. Upon further observation of the contact trajectory, it is evident that the contact area of the clearance joint exhibits a pattern characterized by a large depth and prolonged duration in specific areas.

This study conducts a dynamic response analysis on the crank–slider mechanism with clearance joints. The FEM and the FPM are used to simulate and analyze the crank–slider mechanism under ideal joint and clearance joint conditions, respectively. It is indicated that there is a great agreement in the simulation outcomes of the FPM and FEM, which perfectly validates the effectiveness of the dynamic modeling method with clearances proposed in this study for the dynamic simulation analysis of such structures. If we take the joint clearance into consideration, the dynamic response of the structure significantly deviates from the response under ideal joint conditions. In particular, the structure's acceleration generates high-frequency and large-amplitude oscillations. The contact within the clearance joints also exhibits large depth and long-duration characteristics. These results fully prove the effectiveness of the proposed clearance joint contact force model in this study.

5 Conclusion

This study proposes a hyperbolic contact surface Winkler model for spherical clearance joints. The proposed hybrid theoretical and numerical model incorporates a novel basic formula and an optimized variable coefficient. The new basic formula, derived from the hyperbolic Winkler model, features a superior variable exponent. Meanwhile, the variable coefficient is optimized using a surrogate function developed from the response surface method. Comprehensive comparison and dynamic analyses are conducted to evaluate the effectiveness of the proposed contact force model.

The main conclusions can be summarized as follows:

  1. The proposed contact force model exhibits enhanced precision. It is characterized by a superior variable exponent that captures the trend relationship between the contact force and the penetration depth, and an optimized variable coefficient that effectively adjusts the contact force precisely.

  2. The dimensionless expression of the proposed model is reasonable. It simplifies the form of the contact force model, reduces the number of parameters to be considered, and elucidates the essence of contact problems with small-clearance sizes.

  3. The proposed contact force model has a wide range of applications. The geometric properties of the contact bodies do not limit its application. The proposed model demonstrates significant applicability when the socket radius is substantially larger than the clearance size, which is regarded as a quasi-conformal contact problem.

This study proposes a novel contact force model for spherical clearance joints and validates its accuracy and effectiveness through multiple analyses. Future research will take the energy dissipation during the impact process of the clearance joint into consideration and establish a continuous contact force model based on the model proposed in this study. In addition, in large deployable structures working in orbit, many spherical joints are used as deployment joints, where clearance is inevitable. The nonlinear factors such as contact, impact, and friction caused by the clearance affect the dynamic characteristics and response of the structure. Therefore, establishing an accurate nonlinear impact contact force model will be helpful in establishing a more accurate structural dynamics model; thus, a more precise analysis and prediction of the structure's dynamic response will become a reality.

Author Contribution Statement

All authors contributed to the study's conception and design. Huaibo Yao, Lei Liang, and Huibo Zhang performed material preparation, data collection, and analysis. Huaibo Yao wrote the first draft of the paper and all authors commented on previous versions of the paper. All authors read and approved the final paper.

Funding Data

  • National Natural Science Foundation of China (Grant No. U21B2075; Funder ID: 10.13039/501100001809).

Conflict of Interest

The authors have no competing interests to declare that are relevant to the content of this article.

Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

Ethics Approval and Consent to Participate

This work does not involve human participants, human material, or human data.

Nomenclature

a =

radius of contact region (m)

cd =

dynamic correction coefficient

cr =

coefficient of restitution

e =

eccentricity vector (m)

E =

Young's modulus of contact bodies (Pa)

E* =

composite modulus of contact bodies (Pa)

FN =

normal contact force (N)

FT =

tangential friction force (N)

Fcon =

particle's total contact force vector (N)

Fext =

particle's external force vector (N)

Fint =

particle's internal force vector (N)

h =

foundation's height (m)

I =

inertia matrix of particle (kg m2)

K =

foundation's elastic modulus (Pa)

m =

particle mass (kg)

Mext =

particle's external moment vector (N·m)

Mint =

particle's internal moment vector (N·m)

p =

contact pressure distribution function (N)

P =

total contact force (N)

P* =

dimensionless contact force

p0 =

contact center's maximum pressure (N)

r =

distance between the contact point and the contact center (m)

R =

contact region's equivalent radius (m)

RB =

ball's outer radius (m)

RS =

socket's inner radius (m)

R2* =

dimensionless socket's inner radius

R2 =

multicorrelation coefficient

Radj2 =

adjusted multicorrelation coefficient

RMSE =

root-mean-squared error

u¯z =

normal elastic displacement in the contact region (m)

v0, v1 =

tolerance velocity of friction (m/s)

vT =

relative tangential velocity (m/s)

xB =

position vector of ball's center (m)

xS =

position vector of socket's center (m)

z =

contact bodies' profile function (m)

α =

contact region's semi-angle (rad)

δ =

penetration depth of contact (m)

δ* =

dimensionless penetration depth

δ˙ =

instantaneous normal impact velocity (m/s)

δ˙() =

initial normal impact velocity (m/s)

ε =

relative error

μk =

sliding friction coefficient

υ =

Poisson's ratio of contact bodies

References

1.
Thornton
,
E. A.
, and
Paul
,
D. B.
,
1985
, “
Thermal-Structural Analysis of Large Space Structures—An Assessment of Recent Advances
,”
J. Spacecr. Rockets
,
22
(
4
), pp.
385
393
.10.2514/3.25761
2.
Li
,
T.
,
Jiang
,
J.
,
Deng
,
H.
,
Lin
,
Z.
, and
Wang
,
Z.
,
2013
, “
Form-Finding Methods for Deployable Mesh Reflector Antennas
,”
Chin. J. Aeronaut.
,
26
(
5
), pp.
1276
1282
.10.1016/j.cja.2013.04.062
3.
Puig
,
L.
,
Barton
,
A.
, and
Rando
,
N.
,
2010
, “
A Review on Large Deployable Structures for Astrophysics Missions
,”
Acta Astronaut.
,
67
(
1–2
), pp.
12
26
.10.1016/j.actaastro.2010.02.021
4.
Gruber
,
P.
,
Häuplik
,
S.
,
Imhof
,
B.
,
Özdemir
,
K.
,
Waclavicek
,
R.
, and
Perino
,
M. A.
,
2007
, “
Deployable Structures for a Human Lunar Base
,”
Acta Astronaut.
,
61
(
1–6
), pp.
484
495
.10.1016/j.actaastro.2007.01.055
5.
Xu
,
Y.
,
Guan
,
F.
,
Chen
,
J.
, and
Zheng
,
Y.
,
2012
, “
Structural Design and Static Analysis of a Double-Ring Deployable Truss for Mesh Antennas
,”
Acta Astronaut.
,
81
(
2
), pp.
545
554
.10.1016/j.actaastro.2012.09.004
6.
Liu
,
M.
,
Cao
,
D.
, and
Zhu
,
D.
,
2020
, “
Equivalent Dynamic Model of the Space Antenna Truss With Initial Stress
,”
AIAA J.
,
58
(
4
), pp.
1851
1863
.10.2514/1.J058647
7.
Zhang
,
J.
,
Deng
,
Z.
,
Guo
,
H.
, and
Liu
,
R.
,
2014
, “
Equivalence and Dynamic Analysis for Jointed Trusses Based on Improved Finite Elements
,”
Proc. Inst. Mech. Eng., Part K: J. Multi-Body Dyn.
,
228
(
1
), pp.
47
61
.10.1177/1464419313512469
8.
Siriguleng
,
B.
,
Zhang
,
W.
,
Liu
,
T.
, and
Liu
,
Y. Z.
,
2020
, “
Vibration Modal Experiments and Modal Interactions of a Large Space Deployable Antenna With Carbon Fiber Material and Ring-Truss Structure
,”
Eng. Struct.
,
207
, p.
109932
.10.1016/j.engstruct.2019.109932
9.
Bai
,
Z.
,
Ning
,
Z.
, and
Zhou
,
J.
,
2022
, “
Study on Wear Characteristics of Revolute Clearance Joints in Mechanical Systems
,”
Micromachines
,
13
(
7
), p.
1018
.10.3390/mi13071018
10.
Guo
,
H.
,
Zhang
,
J.
,
Liu
,
R.
, and
Deng
,
Z.
,
2013
, “
Effects of Joint on Dynamics of Space Deployable Structure
,”
Chin. J. Mech. Eng.
,
26
(
5
), pp.
861
872
.10.3901/CJME.2013.05.861
11.
Quinn
,
D. D.
,
2012
, “
Modal Analysis of Jointed Structures
,”
J. Sound Vib.
,
331
(
1
), pp.
81
93
.10.1016/j.jsv.2011.08.017
12.
Yoshida
,
T.
,
2006
, “
Dynamic Characteristic Formulations for Jointed Space Structures
,”
J. Spacecr. Rockets
,
43
(
4
), pp.
771
779
.10.2514/1.15537
13.
Qi
,
Z.
,
Xu
,
Y.
,
Luo
,
X.
, and
Yao
,
S.
,
2010
, “
Recursive Formulations for Multibody Systems With Frictional Joints Based on the Interaction Between Bodies
,”
Multibody Syst. Dyn.
,
24
(
2
), pp.
133
166
.10.1007/s11044-010-9213-z
14.
Yao
,
H.
,
Huang
,
Y.
,
Ma
,
W.
,
Liang
,
L.
, and
Zhao
,
Y.
,
2023
, “
Dynamic Analysis of a Large Deployable Space Truss Structure Considering Semi-Rigid Joints
,”
Aerospace
,
10
(
9
), p.
821
.10.3390/aerospace10090821
15.
Flores
,
P.
,
Ambrósio
,
J.
,
Claro
,
J. C. P.
,
Lankarani
,
H. M.
, and
Koshy
,
C. S.
,
2006
, “
A Study on Dynamics of Mechanical Systems Including Joints With Clearance and Lubrication
,”
Mech. Mach. Theory
,
41
(
3
), pp.
247
261
.10.1016/j.mechmachtheory.2005.10.002
16.
Flores
,
P.
,
Ambrósio
,
J.
,
Claro
,
J. P.
, and
Lankarani
,
H. M.
,
2007
, “
Dynamic Behaviour of Planar Rigid Multi-Body Systems Including Revolute Joints With Clearance
,”
Proc. Inst. Mech. Eng., Part K: J. Multi-Body Dyn.
,
221
(
2
), pp.
161
174
.10.1243/14644193JMBD96
17.
Flores
,
P.
,
2010
, “
A Parametric Study on the Dynamic Response of Planar Multibody Systems With Multiple Clearance Joints
,”
Nonlinear Dyn.
,
61
(
4
), pp.
633
653
.10.1007/s11071-010-9676-8
18.
Flores
,
P.
,
Leine
,
R.
, and
Glocker
,
C.
,
2011
, “
Modeling and Analysis of Rigid Multibody Systems With Translational Clearance Joints Based on the Nonsmooth Dynamics Approach
,”
Multibody Dynamics: Computational Methods and Applications
, Springer, Berlin, Germany, pp.
107
130
.10.1007/978-90-481-9971-6
19.
Zhao
,
Y.
, and
Bai
,
Z. F.
,
2011
, “
Dynamics Analysis of Space Robot Manipulator With Joint Clearance
,”
Acta Astronaut.
,
68
(
7–8
), pp.
1147
1155
.10.1016/j.actaastro.2010.10.004
20.
Marques
,
F.
,
Isaac
,
F.
,
Dourado
,
N.
,
Souto
,
A. P.
,
Flores
,
P.
, and
Lankarani
,
H. M.
,
2017
, “
A Study on the Dynamics of Spatial Mechanisms With Frictional Spherical Clearance Joints
,”
ASME J. Comput. Nonlinear Dyn.
,
12
(
5
), p.
051013
.10.1115/1.4036480
21.
Hou
,
J.
,
Yao
,
G.
, and
Huang
,
H.
,
2018
, “
Dynamic Analysis of a Spatial Mechanism Including Frictionless Spherical Clearance Joint With Flexible Socket
,”
ASME J. Comput. Nonlinear Dyn.
,
13
(
3
), p.
031002
.10.1115/1.4038508
22.
Tian
,
Q.
,
Zhang
,
Y.
,
Chen
,
L.
, and
Flores
,
P.
,
2009
, “
Dynamics of Spatial Flexible Multibody Systems With Clearance and Lubricated Spherical Joints
,”
Comput. Struct.
,
87
(
13–14
), pp.
913
929
.10.1016/j.compstruc.2009.03.006
23.
Erkaya
,
S.
,
2018
, “
Experimental Investigation of Flexible Connection and Clearance Joint Effects on the Vibration Responses of Mechanisms
,”
Mech. Mach. Theory
,
121
, pp.
515
529
.10.1016/j.mechmachtheory.2017.11.014
24.
Chen
,
X.
, and
Sun
,
C.
,
2019
, “
Dynamic Modeling of Spatial Parallel Mechanism With Multi-Spherical Joint Clearances
,”
Int. J. Adv. Rob. Syst.
,
16
(
5
), p.
172988141987591
.10.1177/1729881419875910
25.
Chen
,
X.
,
Li
,
Y.
, and
Jia
,
Y.
,
2019
, “
Dynamic Response and Nonlinear Characteristics of Spatial Parallel Mechanism With Spherical Clearance Joint
,”
ASME J. Comput. Nonlinear Dyn.
,
14
(
4
), p.
041010
.10.1115/1.4042636
26.
Chen
,
X.
,
Wu
,
R.
, and
Jia
,
Y.
,
2022
, “
Dynamic Response and Chaotic Characteristics Analysis of Spatial Parallel Mechanisms With Multiple Spherical Clearance Joints
,”
J. Vib. Eng. Technol.
, 10, pp.
1
27
.10.1007/s42417-021-00387-7
27.
Wang
,
G.
, and
Wang
,
L.
,
2019
, “
Dynamics Investigation of Spatial Parallel Mechanism Considering Rod Flexibility and Spherical Joint Clearance
,”
Mech. Mach. Theory
,
137
, pp.
83
107
.10.1016/j.mechmachtheory.2019.03.017
28.
Jing
,
Q.
, and
Liu
,
H.
,
2023
, “
Dynamics and Wear Prediction of Mechanisms Considering Multiple Clearances and Coatings
,”
Lubricants
,
11
(
7
), p.
310
.10.3390/lubricants11070310
29.
Ma
,
X. F.
, and
Li
,
T. J.
,
2018
, “
Wave Analysis of Planar Deployable Structures With Revolute Clearance Joints Based on Spectral Element Method
,”
Int. J. Appl. Mech.
,
10
(
8
), p.
1850090
.10.1142/S1758825118500904
30.
Johnson
,
K. L.
,
1987
,
Contact Mechanics
,
Cambridge University Press
,
London
.
31.
Lankarani
,
H. M.
, and
Nikravesh
,
P. E.
,
1994
, “
Continuous Contact Force Models for Impact Analysis in Multibody Systems
,”
Nonlinear Dyn.
,
5
(
2
), pp.
193
207
.10.1007/BF00045676
32.
Hunt
,
K. H.
, and
Crossley
,
F. R. E.
,
1975
, “
Coefficient of Restitution Interpreted as Damping in Vibroimpact
,”
ASME J. Appl. Mech.
,
42
(
2
), pp.
440
445
.10.1115/1.3423596
33.
Goodman
,
L. E.
, and
Keer
,
L.
,
1965
, “
The Contact Stress Problem for an Elastic Sphere Indenting an Elastic Cavity
,”
Int. J. Solids Struct.
,
1
(
4
), pp.
407
415
.10.1016/0020-7683(65)90005-3
34.
Noble
,
B.
, and
Hussain
,
M. A.
,
1969
, “
Exact Solution of Certain Dual Series for Indentation and Inclusion Problems
,”
Int. J. Eng. Sci.
,
7
(
11
), pp.
1149
1161
.10.1016/0020-7225(69)90081-0
35.
Liu
,
C. S.
,
Zhang
,
K.
, and
Yang
,
L.
,
2006
, “
Normal Force-Displacement Relationship of Spherical Joints With Clearances
,”
ASME J. Comput. Nonlinear Dyn.
,
1
(
2
), pp.
160
167
.10.1115/1.2162872
36.
Fang
,
X.
,
Zhang
,
C.
,
Chen
,
X.
,
Wang
,
Y.
, and
Tan
,
Y.
,
2015
, “
A New Universal Approximate Model for Conformal Contact and Non-Conformal Contact of Spherical Surfaces
,”
Acta Mech.
,
226
(
6
), pp.
1657
1672
.10.1007/s00707-014-1277-z
37.
Ascione
,
L.
, and
Grimaldi
,
A.
,
1984
, “
Unilateral Contact Between a Plate and an Elastic Foundation
,”
Meccanica
,
19
(
3
), pp.
223
233
.10.1007/BF01743736
38.
Ascione
,
L.
, and
Olivito
,
R. S.
,
1985
, “
Unbonded Contact of a Mindlin Plate on an Elastic Half-Space
,”
Meccanica
,
20
(
1
), pp.
49
58
.10.1007/BF02337062
39.
Zhang
,
J.
,
Guo
,
H.
,
Liu
,
R.
, and
Deng
,
Z.
,
2015
, “
Nonlinear Characteristic of Spherical Joints With Clearance
,”
J. Aerosp. Technol. Manage.
,
7
(
2
), pp.
179
184
.10.5028/jatm.v7i2.464
40.
Sassi
,
M.
, and
Desvignes
,
M.
,
1998
, “
A Seminumerical Method for Three-Dimensional Frictionless Contact Problems
,”
Math. Comput. Modell.
,
28
(
4–8
), pp.
413
425
.10.1016/S0895-7177(98)00131-9
41.
Rokach
,
I. V.
,
2003
, “
On the Accurate Determination of Contact Compliance for Impact Test Modelling
,”
Int. J. Solids Struct.
,
40
(
11
), pp.
2715
2729
.10.1016/S0020-7683(03)00088-X
42.
Pereira
,
C.
,
Ramalho
,
A.
, and
Ambrosio
,
J.
,
2015
, “
An Enhanced Cylindrical Contact Force Model
,”
Multibody Syst. Dyn.
,
35
(
3
), pp.
277
298
.10.1007/s11044-015-9463-x
43.
Yao
,
H.
,
Liang
,
L.
,
Ma
,
W.
,
Zhang
,
H.
, and
Zhao
,
Y.
,
2024
, “
An Enhanced Circumferential Winkler Contact Force Model for Cylindrical Clearance Joints
,”
Int. J. Appl. Mech.
,
16
(
2
), p.
2450019
.10.1142/S1758825124500194
44.
Liu
,
X.
,
Zhao
,
L.
,
Mao
,
J.
, and
Li
,
T.
,
2020
, “
Tangential Force Model for the Combined Finite-Discrete Element Method
,”
Int. J. Comput. Methods
,
17
(
9
), p.
1950068
.10.1142/S0219876219500683
45.
Bai
,
X.
,
Dong
,
Q.
,
Zheng
,
H.
, and
Zhou
,
K.
,
2021
, “
A Finite Element Model for Non-Newtonian Starved Thermal-Elastohydrodynamic Lubrication of 3D Line Contact
,”
Int. J. Appl. Mech.
,
13
(
9
), p.
2150107
.10.1142/S1758825121501076
46.
Ding
,
Y.
,
Liang
,
X.
, and
Wang
,
G.
,
2022
, “
An Incremental Contact Model for Rough Surfaces of Strain-Hardening Solids
,”
Int. J. Appl. Mech.
,
14
(
8
), p.
2250088
.10.1142/S1758825122500880
47.
Jia
,
Y.
, and
Chen
,
X.
,
2022
, “
Application of a New Conformal Contact Force Model to Nonlinear Dynamic Behavior Analysis of Parallel Robot With Spherical Clearance Joints
,”
Nonlinear Dyn.
,
108
(
3
), pp.
2161
2191
.10.1007/s11071-022-07344-3
48.
Wang
,
B.
,
Lyu
,
Q.
,
Jiang
,
L.
,
Chen
,
Y.
, and
Guo
,
Z.
,
2022
, “
An Extended Hertz Model for Incompressible Mooney–Rivlin Half-Space Under Finite Spherical Indentation
,”
Int. J. Appl. Mech.
,
14
(
10
), p.
2250103
.10.1142/S1758825122501034
49.
Ding
,
S.
,
Jian
,
B.
,
Zhang
,
Y.
,
Hu
,
Y.
,
Xia
,
R.
, and
Hu
,
G.
,
2023
, “
A Finite Element Solution to Normal Contact Forces of Viscoelastic Particles
,”
Int. J. Appl. Mech.
,
15
(
1
), p.
2350003
.10.1142/S1758825123500035
50.
Yaylacı
,
M.
,
Abanoz
,
M.
,
Yaylacı
,
E. U.
,
Ölmez
,
H.
,
Sekban
,
D. M.
, and
Birinci
,
A.
,
2022
, “
Evaluation of the Contact Problem of Functionally Graded Layer Resting on Rigid Foundation Pressed Via Rigid Punch by Analytical and Numerical (FEM and MLP) Methods
,”
Arch. Appl. Mech.
,
92
(
6
), pp.
1953
1971
.10.1007/s00419-022-02159-5
51.
Yaylacı
,
M.
,
Şabano
,
B. Ş.
,
Özdemir
,
M. E.
, and
Birinci
,
A.
,
2022
, “
Solving the Contact Problem of Functionally Graded Layers Resting on a HP and Pressed With a Uniformly Distributed Load by Analytical and Numerical Methods
,”
Struct. Eng. Mech., Int. J.
,
82
(
3
), pp.
401
416
.10.12989/sem.2022.82.3.401
52.
Lan
,
W.
,
Fan
,
S.
, and
Fan
,
S.
,
2021
, “
A Fractal Model of Elastic–Plastic Contact Between Rough Surfaces for a Low-Velocity Impact Process
,”
Int. J. Comput. Methods
,
18
(
9
), p.
2150039
.10.1142/S0219876221500390
53.
Zheng
,
Y.
,
Yang
,
C.
,
Wan
,
H. P.
,
Luo
,
Y.
,
Li
,
Y.
, and
Yu
,
Y.
,
2020
, “
Dynamics Analysis of Spatial Mechanisms With Dry Spherical Joints With Clearance Using Finite Particle Method
,”
Int. J. Struct. Stab. Dyn.
,
20
(
3
), p.
2050035
.10.1142/S0219455420500352
54.
Alves
,
J.
,
Peixinho
,
N.
,
da Silva
,
M. T.
,
Flores
,
P.
, and
Lankarani
,
H. M.
,
2015
, “
A Comparative Study of the Viscoelastic Constitutive Models for Frictionless Contact Interfaces in Solids
,”
Mech. Mach. Theory
,
85
, pp.
172
188
.10.1016/j.mechmachtheory.2014.11.020
55.
Pereira
,
C. M.
,
Ramalho
,
A. L.
, and
Ambrósio
,
J. A.
,
2011
, “
A Critical Overview of Internal and External Cylinder Contact Force Models
,”
Nonlinear Dyn.
,
63
(
4
), pp.
681
697
.10.1007/s11071-010-9830-3
56.
Dubowsky
,
S.
, and
Freudenstein
,
F.
,
1971
, “
Dynamic Analysis of Mechanical Systems With Clearances—Part 1: Formation of Dynamic Model
,”
ASME J. Manuf. Sci. Eng.
,
93
(
1
), pp.
305
309
.10.1115/1.3427895
57.
Bauchau
,
O. A.
, and
Rodriguez
,
J.
,
2002
, “
Modeling of Joints With Clearance in Flexible Multibody Systems
,”
Int. J. Solids Struct.
,
39
(
1
), pp.
41
63
.10.1016/S0020-7683(01)00186-X
58.
Ting
,
E. C.
,
Shih
,
C.
, and
Wang
,
Y. K.
,
2004
, “
Fundamentals of a Vector Form Intrinsic Finite Element: Part I. Basic Procedure and a Plane Frame Element
,”
J. Mech.
,
20
(
2
), pp.
113
122
.10.1017/S1727719100003336
59.
Lankarani
,
H. M.
, and
Nikravesh
,
P. E.
,
1990
, “
A Contact Force Model With Hysteresis Damping for Impact Analysis of Multibody Systems
,”
ASME J. Mech. Des.
,
112
(
3
), pp.
369
376
.10.1115/1.2912617
60.
Ambrósio
,
J. A. C.
,
2003
, “
Impact of Rigid and Flexible Multibody Systems: Deformation Description and Contact Models
,”
Virtual Nonlinear Multibody Systems
, Springer, Dordrecht, The Netherlands, pp.
57
81
.10.1007/978-94-010-0203-5