A Hyperbolic Contact Surface Winkler Contact Force Model for Spherical Clearance Joints

Large deployable structures are extensively utilized in aerospace missions for their lightweight and small size [1–4]. 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 [5–8]. 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 [10–14]. 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 [15–21]. 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 [24–26] 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,37–39]. 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 [40–43] and the analysis of complex contact-related issues [44–46]. 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.

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

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.

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 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.

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.

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σ$⁠.

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.

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.

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.

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$⁠.

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).

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×1011 Pa$ and $0.3$⁠, respectively.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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].

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.

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].

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.

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.

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

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. 17–21.

Figures 17–19, 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.

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.

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.

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

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

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

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

$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 m²)

$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