A Joint Surface Contact Stiffness Model Considering Micro-Asperity Interaction

: Mechanical joint interfaces are widely found in mechanical equipment, and their contact stiffness directly affects the overall performance of the mechanical system. Based on the fractal and elastoplastic contact mechanics theories, the K-E elastoplastic contact model is introduced to establish the contact stiffness model for mechanical joint interfaces. This model considers the interaction effects between micro-asperities in the fully deformed state, including elasticity, first elastoplasticity, second elastoplasticity, and complete plastic deformation state. Based on this model, the effects of fractal parameters on normal contact stiffness and contact load are analyzed. It can be found that the larger fractal dimension D or smaller characteristic scale coefficient G will weaken the interaction between micro-asperities. The smoother processing surfaces lead to higher contact stiffness in mechanical joint interfaces. The applicability and effectiveness of the proposed model are verified by comparing it with the traditional contact model calculation results. Under the same load, the interaction between micro-rough surfaces leads to an increase in both overall deformation and contact stiffness. The accuracy of the predicted contact stiffness model is also validated by comparing it with experimental results.


Introduction
Contact interfaces are indispensable in the mechanical structure of aero-engines, including bearings and bolted joints, which ensure tight connections and stable operation between different components [1].The contact problem of the rough surface is a fundamental barrier for accurately describing the interface characteristics and analysis of the static and dynamic behavior of the mechanical system [2].Considering the normal contact stiffness in the joint interface is a challenge when studying the mechanical characteristics of aero-engines.
In practice, engineering surfaces are inherently rough.At the microscopic scale, the contact interface is formed by the mutual contact between micro-asperities of different sizes [3,4].In practical engineering research, the contact of asperities will affect the mechanical properties of the joint structure.The total stiffness of these contacting micro-asperities constitutes the contact stiffness of the joint surface [2].In practice, the contact stiffness at the component interface typically accounts for 60% to 80% of the structure's overall stiffness [5,6].The interactions among micro-asperities cannot be ignored, and establishing a high-precision contact stiffness model is of significant importance in studying the overall stability and dynamic characteristics of aircraft engines [2][3][4][5][6][7].
Rough surface contact models can be classified into deterministic, statistical, and fractal contact models.Statistical contact models are widely employed in contact analysis because statistical models can effectively capture the randomness and quantify the uncertainty of contacting surfaces [8][9][10][11][12].However, the model is characterized by statistical parameters, which are inevitably limited by the resolution of measurement instruments and sampling length [13,14].Statistical contact models require high-resolution instruments and substantial computational resources [15,16].Deterministic models are also affected by sampling intervals [17].The fractal contact model is described using fractal parameters, which solves the model scale dependence and ensures the uniqueness of the surface contact analysis results [18].
Scholars have conducted extensive research on contact interface modeling.Zhang et al. [19] established the normal contact stiffness model of the mechanical joint based on fractal theory.Wen et al. [20] proposed the fractal contact stiffness model considering the effect of the domain expansion factor in the area distribution function of micro-asperity.Yang et al. [21] established the micro asperity contact stiffness model covering regimes from elasticity to elastoplasticity and then to total plastic deformation, and derived the relationship between contact load and contact stiffness for different plastic indices.Fractal models for normal and tangential stiffness considering the friction impact have been developed in the literature [22,23].In summary, the aforementioned fractal models neglect the influence of interactions between micro-asperities, resulting in limited reference values of the established models.Wang et al. [24] proposed the fractal stiffness model considering the interactions among micro-asperities.Xiao et al. [25] developed an approximate expression for the normal contact stiffness in terms of normal load based on the K-E statistical contact model and then analyzed the influence of surface roughness parameters on the normal contact stiffness.Kang et al. [26] considered the coupling effect of macroscopic profile deviation and mesoscale roughness on the contact behavior, resulting in the proposal of an improved rough surface contact model to describe the normal and tangential contact behavior of rough interfaces with multi-scale morphology.Considering the actual contact area of the single asperity, Li et al. [27] proposed a multi-scale contact mechanics modeling method.Later, based on the fractal theory, the influence of contact load, fractal parameters and friction factor on the real contact area was evaluated.Similarly, Sun et al. [28] established the thermal contact fractal model with the consideration of friction coefficient.Pan et al. [29] used a novel truncated area size distribution function to construct a multiasperity fractal contact model.To construct a more comprehensive contact mechanics model, Yu et al. [30] constructed a multi-stage contact model of the fractal curved surface based on the attitude angle of the asperities.Li et al. [31] proposed a multi-physics electrical contact model considering the interaction of micro-asperities and investigated the influence of surface fractal parameters on electrical contact behavior.Shen et al. [32] studied the size distribution of truncated spots on the fractal surfaces by employing the three-dimensional Weierstrass-Mandelbrot function to characterize the fractal surface.However, the influence of micro-asperity interactions in this model is only reflected in the elastic deformation stage, and the model only includes two parts: elastic and fully plastic deformation stages, which have certain limitations.Up to now, contact mechanics modeling that considers the complex interaction mechanism between micro-asperities has not been fully researched.Establishing a more precise rough surface contact mechanics model is crucial for investigating the contact behavior and predicting the contact stiffness of fixed joint interfaces.
In this study, a new contact stiffness modeling method is proposed based on the influence of asperity interaction to enhance the predictive accuracy of contact models.The proposed model comprehensively considers the elastic-plastic deformation mechanism of micro-asperities and their interactions.Based on the fractal theory and elastic-plastic contact mechanics theory, the K-E elastic-plastic contact theory is introduced.The contact stiffness model for mechanical joint surfaces is established, which encompasses the complete deformation regimes from elasticity, first elastoplasticity, second elastoplasticity, to total plastic deformation, and considers the effects of micro-asperity interactions.Subsequently, the influence of fractal parameters on the normal contact stiffness of the joint interface is determined.

Single Micro-Asperity Body Deformation Mechanism
The following assumptions are proposed to simplify calculations and analysis: a. Hardening of micro-asperities during the deformation process is not considered; b.The contact between two rough surfaces can be equivalently represented as the contact between a smooth plane and a rough surface; c.The distribution of micro-asperities on the rough surface follows the Gaussian distribution.
The contact between the single micro-asperity and the rigid contact plane is illustrated in Figure 1.h denotes the height of the micro-asperity, R represents the curvature radius of the single micro-asperity, l denotes the base length of the micro-asperity, l r and l t represent the actual contact length and the truncation length, respectively.Figure 1 illustrates that when micro-asperities are subjected to load from the rigid plane, they will also experience interactions with adjacent micro-asperities.These interactions not only result in the deformation of the micro-asperities but also transfer the load to the substrate, causing additional deformation in the normal direction of the micro-asperities' bottom surfaces [33], and this additional deformation affects the actual contact area and contact load of the micro-asperity.In this study, the model focuses on variations in the normal direction of the micro-asperities and establishes the fractal contact model that considers interactions between micro-asperities.

Single Micro-Asperity Body Deformation Mechanism
The following assumptions are proposed to simplify calculations and analysis: a. Hardening of micro-asperities during the deformation process is not considered; b.The contact between two rough surfaces can be equivalently represented as the contact between a smooth plane and a rough surface; c.The distribution of micro-asperities on the rough surface follows the Gaussian distribution.
The contact between the single micro-asperity and the rigid contact plane is illustrated in Figure 1.h denotes the height of the micro-asperity, R represents the curva- ture radius of the single micro-asperity, l denotes the base length of the micro-asperity, r l and t l represent the actual contact length and the truncation length, respectively.Figure 1 illustrates that when micro-asperities are subjected to load from the rigid plane, they will also experience interactions with adjacent micro-asperities.These interactions not only result in the deformation of the micro-asperities but also transfer the load to the substrate, causing additional deformation in the normal direction of the micro-asperities' bottom surfaces [33], and this additional deformation affects the actual contact area and contact load of the micro-asperity.In this study, the model focuses on variations in the normal direction of the micro-asperities and establishes the fractal contact model that considers interactions between micro-asperities.Based on the W-M function, the surface microtopography height function before deformation can be defined as follows [34].Based on the fractal theory, the size of asperities is related to fractal parameters.From Equation (1), when the contact area a and the base length l satisfy a specific relation- ship, 2 a l = .The asperity curvature radius R can be described as: Based on the W-M function, the surface microtopography height function before deformation can be defined as follows [34].
where l = 1 γ n denotes the base length, n represents the frequency coefficient of microasperities, and γ n is the spatial frequency of the surface profile.
Based on the fractal theory, the size of asperities is related to fractal parameters.From Equation (1), when the contact area a and the base length l satisfy a specific relationship, a = l 2 .The asperity curvature radius R can be described as: The contact deformation process sequentially undergoes three stages: elastic deformation, elastoplastic deformation, and plastic deformation [35].Furthermore, the elastoplastic deformation regime can be further subdivided into the first and second elastoplastic deformation.In this study, the subscripts "e", "ep1", "ep2", and "p" indicate the different deformation stages.Additionally, "ec", "epc", and "pc" are used to represent parameters related to critical elasticity, first critical elastoplasticity, and second critical elastoplasticity, respectively.

Elastic Deformation Regime
Based on the Hertzian contact theory [36], a e and f e are the contact area and contact load of the single asperity with a e = πRω (3) where E is the equivalent elastic modulus.By combining Equations ( 2)-( 4), the relationship between the contact area and contact load in the elastic deformation stage can be rewritten as The critical deformation ω ec could be expressed as K represents the hardness-related coefficient, which is related to the material's Poisson's ratio and satisfies the relationship K = 0.454 + 0.41v, and H is the hardness coefficient.
The critical contact area can be calculated by

Elastoplastic Deformation Regime
Kogut and Etsion put forward the empirical expressions between contact load and deformation for the single asperity and gave the contact force relationships for the single asperity contact in the different deformation stages [37].The contact area and contact force for different deformation stages can be determined as [37] a ep1 = 0.93πRω ec ω ω ec 1.136 ( 8) The relationship between contact area and contact load in the different deformation stages can be described by The first critical elastic-plastic deformation ω epc and the second critical elastic-plastic deformation ω pc can be expressed as ) The critical contact area can be described as

Plastic Deformation Stage
In the fully plastic stage, the contact load and contact area of the single micro-asperity can be represented as [37] a p = 2πRω ( 18)

Micro-Asperity Deformation
The rough surface is formed by the accumulation of many micro-asperities.As shown in Figure 2, the contact between two rough surfaces is simplified to a rough surface in contact with a smooth surface [38].The increase in the normal deformation of the joint will intensify the relative interaction between the asperities.The deformation caused by load on the micro-asperity is given as where ω indicates the local deformation of the asperities under the contact load, σ de- notes the displacement generated on the average plane of the rough surface, ξ repre- sents the total deformation of the rough surface [39].The total deformation of asperities ξ can be represented by the amplitude of the asperity profile function [39].
The additional amount of deformation can be represented by [40] m 1.12 where f represents the contact load, and m P denotes the contact pressure between the The deformation caused by load on the micro-asperity is given as where ω indicates the local deformation of the asperities under the contact load, σ denotes the displacement generated on the average plane of the rough surface, ξ represents the total deformation of the rough surface [39].The total deformation of asperities ξ can be represented by the amplitude of the asperity profile function [39].
The additional amount of deformation can be represented by [40] σ = 1.12 f P m E where f represents the contact load, and P m denotes the contact pressure between the contact interface.
Combining Equations ( 20)-( 22), the elastic deformation of contact interfaces is as follows: √ Combining Equations ( 12) and ( 23), the deformation in the first elastoplastic stage can be further described as Similarly, based on Equations ( 13) and ( 23), the deformation in the second elastoplastic stage is rewritten as Substituting Equation (19) into Equation ( 23), the deformation in the entire plastic deformation stage is rewritten as

Micro-Asperity Contact Stiffness
Substituting Equation (2) into Equation ( 4), the contact stiffness of the micro-asperity in the elastic deformation stage can be written as where da dω can be obtained via Equation (23).Similarly, substituting Equation (2) into Equation ( 9), the contact stiffness of the microasperity in the first elastoplastic deformation stage can be calculated by  where da dω can be obtained via Equation (24).Combining Equations ( 2) and ( 11), the contact stiffness in the second elastoplastic deformation stage can be described as  where da dω can be determined via Equation (25).

The Contact Characteristics of the Interface 4.1. Normal Contact Load on the Contact Interface
The three-dimensional size-distribution function of the rough surface can be written as [34] where a l represents the largest micro-asperity contact area.Substituting Equation (30) into Equation ( 4), the total elastic contact load can be expressed as Similarly, based on Equations ( 9) and ( 30), the total contact load of micro-asperities in the first elastoplastic deformation stage (ω ec < ω ≤ 6ω ec ) can be calculated by Substituting Equation (30) into Equation ( 11), the total contact load of micro-asperities in the second elastoplastic deformation stage (6ω ec < ω ≤ 110ω ec ) can be calculated as The total contact load in the plastic deformation stage is denoted by To sum up, the total contact load of the rough interface can be described as

The Normal Contact Stiffness of the Contact Interface
Based on Equations ( 27) and ( 30), the normal contact stiffness of the elastic deformation stage can be expressed as Combining Equations ( 28)-( 30), the contact stiffness from the different elastoplastic deformation regimes can be described as It is worth noting that normal contact stiffness only exists when the micro-asperities are in the elastic and elastoplastic deformation stages.Considering the asperity interaction, the contact stiffness of the joint interface can be expressed as [41] K t = K e +K ep1 +K ep2 (39)

Model Simulation
The material and fractal parameters of the rough surface were D = 1.4,1.45, 1.5, 1.55, 1.6 and the nominal contact area was Aa = 1 × 10 −4 m 2 .Considering that the final total contact load is challenging to calculate directly, a new iterative method is employed to solve the proposed model [42].
Figure 3 illustrates the variation trend between the maximum micro-asperity contact area a l and the contact load.As the contact load F t increases, the maximum micro-asperity contact area a l also increases.

Model Simulation
The material and fractal parameters of the rough surface were 1.4 1.45 1.5 1.55 1.6 , and the nominal contact area was .Considering that the final total contact load is challenging to calculate directly, a new iterative method is employed to solve the proposed model [42].
Figure 3 illustrates the variation trend between the maximum micro-asperity contact area l a and the contact load.As the contact load t F increases, the maximum micro- asperity contact area l a also increases.To enhance the accuracy of calculations and eliminate dimensional influences, the following dimensionless topographical parameters are introduced.Figure 5 shows the impact of different characteristic scale coefficients G on the contact stiffness when  = 1.4.It can be seen from Figure 4 that there is a positive correlation between contact stiffness and contact load.As the contact load increases, the contact stiffness gradually rises, which can be attributed to an increase in contact load leading to more micro-asperities making contact on the joint surfaces.Therefore, more micro-asperities share the load, leading to a gradual decrease in the overall deformation of the contact surface.To enhance the accuracy of calculations and eliminate dimensional influences, the following dimensionless topographical parameters are introduced.
Figure 4 illustrates the influence of different fractal dimensions D on the contact stiffness G = 10 −8 .Figure 5 shows the impact of different characteristic scale coefficients G on the contact stiffness when D = 1.4.It can be seen from Figure 4 that there is a positive correlation between contact stiffness and contact load.As the contact load increases, the contact stiffness gradually rises, which can be attributed to an increase in contact load leading to more micro-asperities making contact on the joint surfaces.Therefore, more micro-asperities share the load, leading to a gradual decrease in the overall deformation of the contact surface.Moreover, this study emphasizes the influence of various fractal parameters on the stiffness model.The fractal dimension D describes the complexity and irregularity of rough contact surface profiles across all scales.As shown in Figure 4, with an increase in the fractal dimension D, the contact stiffness also increases because a higher fractal dimension D shows a smoother contact surface.Under the same load, more micro-asperities are in contact, leading to a larger total contact area and reduced sensitivity of the rough surface deformation.The reduced deformation of micro-asperities leads to an enhanced stiffness of the contact surface.Moreover, this study emphasizes the influence of various fractal parameters on the stiffness model.The fractal dimension D describes the complexity and irregularity of rough contact surface profiles across all scales.As shown in Figure 4, with an increase in the fractal dimension D, the contact stiffness also increases because a higher fractal dimension D shows a smoother contact surface.Under the same load, more micro-asperities are in contact, leading to a larger total contact area and reduced sensitivity of the rough surface deformation.The reduced deformation of micro-asperities leads to an enhanced stiffness of the contact surface.
The characteristic scale factor G reflects the size of the micro-asperities and contact surfaces in the vertical direction.The more considerable value indicates greater heights of rough micro-asperities, resulting in a rougher surface.As shown in Figure 5, it is evident that a more minor scale factor corresponds to greater contact stiffness.Moreover, smoother surfaces significantly impact the contact stiffness of rough surfaces (smaller G values leading to more extensive slopes), indicating a more pronounced interaction among micro-asperities.
When the characteristic scale factor G increases, the contact surface becomes increasingly rough, making micro-asperities more sensitive to elastic and elastoplastic deformation and increasing the total deformation.An increase in deformation indicates a decrease in contact stiffness and reduces the ability of micro-asperities to balance normal loads.Therefore, contact stiffness can be increased by reducing surface roughness and improving smoothness.

Model Validation
In order to evaluate the applicability and effectiveness of the proposed model, a comparison is made with a stiffness fractal model that does not consider the interaction between micro-asperities in this study.The fractal parameters are D = 1.6, G = 1 × 10 −9 .As shown in Figure 6, when the load is relatively small, the difference between these two models is not significant.However, when the contact load exceeds 10 −5 , with the contact load gradually increasing, the contact stiffness of the model without considering the interaction between asperities becomes significantly lower than the one considering interaction.
Aerospace 2024, 11, x FOR PEER REVIEW 13 of 17 In order to evaluate the applicability and effectiveness of the proposed model, a comparison is made with a stiffness fractal model that does not consider the interaction between micro-asperities in this study.The fractal parameters are .As shown in Figure 6, when the load is relatively small, the difference between these two models is not significant.However, when the contact load exceeds 5 10 − , with the contact load gradually increasing, the contact stiffness of the model without considering the interaction between asperities becomes significantly lower than the one considering interaction.
The reason is that considering interactions among micro-asperities leads to an increase in the total contact area of the contact surface, resulting in increasing overall deformation and reducing contact stiffness.Additionally, the applicability of the proposed model was verified by comparing it with experimental data [43].The equivalent parameters under grinding machining were adopted, and the relevant mechanical parameters are shown in Table 1.Compared to the JZZ model [43], the error between the calculation results of the pro- The reason is that considering interactions among micro-asperities leads to an increase in the total contact area of the contact surface, resulting in increasing overall deformation and reducing contact stiffness.
Additionally, the applicability of the proposed model was verified by comparing it with experimental data [43].The equivalent parameters under grinding machining were adopted, and the relevant mechanical parameters are shown in Table 1.Compared to the JZZ model [43], the error between the calculation results of the proposed model and the experimental data in this research significantly decreased, as shown in Figure 7. Additionally, the proposed model exhibited consistent trends with the experimental data, indicating that the model can reasonably evaluate the contact stiffness of the joint surfaces.Figure 8 depicts the correlation between the dimensionless normal stiffness of the joint interface and the dimensionless contact load, as determined by different contact mechanics models.The comparison results show that the proposed model is in basic agreement with the classic contact mechanics model, further illustrating the applicability and effectiveness of the model in this research.
The model proposed in this research has two main limitations: the area distribution function of the surface roughness elements was not corrected, and the influence of frictional factors on the contact surface was not considered.These factors may diminish the accuracy of the model; further research is required to address these limitations.

Conclusions
Based on the M-B fractal contact mechanics theory, the normal contact stiffness fractal model that considers the interaction between micro-asperities is established in this study.The proposed fractal contact model will provide a theoretical reference for accurately pre-

Conclusions
Based on the M-B fractal contact mechanics theory, the normal contact stiffness fractal model that considers the interaction between micro-asperities is established in this study.The proposed fractal contact model will provide a theoretical reference for accurately predicting the mechanical behavior of fixed joint interfaces.The main conclusions are as follows: (1) The introduction of the K-E elastoplastic contact theory has led to the development of a comprehensive contact stiffness model, considering the impact of micro-asperity interactions throughout the entire deformation process, including elasticity, first elastoplasticity, second elastoplasticity, and complete plastic deformation.(2) The influence of fractal parameters on contact stiffness was analyzed.Specifically, the increase in fractal dimension D and decrease in scale parameter G will weaken the interaction among asperities, reducing deformation and increasing contact stiffness.
In other words, smoother machined surfaces exhibit higher contact stiffness.(3) Considering the interaction of micro-asperities will lead to an overall increase in the total deformation, resulting in lower stiffness values calculated using the model proposed in this paper.(4) With the consideration of asperity interaction, the established fractal contact model is suitable for fixed joint interfaces, such as bolted joint structures, but not for moveable interfaces.
However, this study focused on the relative interaction between the asperities without considering the area distribution function of the contact asperities and the influence of the frictional factor.Some meaningful work can be carried out in future research, such as constructing the fractal mechanical model with the corrected area distribution function of the rough interfaces and analyzing the effects of frictional factors on the mechanical model that considers micro-asperity interactions.

Figure 1 .
Figure 1.Deformation diagram of the single asperity.
length, n represents the frequency coefficient of mi- cro-asperities, and n γ is the spatial frequency of the surface profile.

Figure 1 .
Figure 1.Deformation diagram of the single asperity.

Figure 2 .
Figure 2. Schematic diagram of the rough surface contact.

Figure 2 .
Figure 2. Schematic diagram of the rough surface contact. 2544 1021

Aerospace 2024 , 17 5.
11, x FOR PEER REVIEW 9 of Simulation and Result Analysis of the Normal Stiffness Model for the Contact Interface

Figure 3 .
Figure 3.The relationship between the contact load and the contact area of the largest asperity.

Figure 4
Figure 4 illustrates the influence of different fractal dimensions D on the contact stiffness 8 10 G − = .Figure5shows the impact of different characteristic scale coefficients G on the contact stiffness when  = 1.4.It can be seen from Figure4that there is a positive correlation between contact stiffness and contact load.As the contact load increases, the contact stiffness gradually rises, which can be attributed to an increase in contact load leading to more micro-asperities making contact on the joint surfaces.Therefore, more micro-asperities share the load, leading to a gradual decrease in the overall deformation of the contact surface.

Figure 3 .
Figure 3.The relationship between the contact load and the contact area of the largest asperity.

Figure 4 .
Figure 4.The influence of different fractal dimensions on contact stiffness.Figure 4. The influence of different fractal dimensions on contact stiffness.

Figure 4 .
Figure 4.The influence of different fractal dimensions on contact stiffness.Figure 4. The influence of different fractal dimensions on contact stiffness.

6 Figure 5 .
Figure 5.The influence of different characteristic scale coefficients on contact stiffness.

Figure 5 .
Figure 5.The influence of different characteristic scale coefficients on contact stiffness.

Figure 6 .
Figure 6.Comparison of contact stiffness calculated by different models.

Figure 6 .
Figure 6.Comparison of contact stiffness calculated by different models.
Aerospace 2024, 11, x FOR PEER REVIEW 14 of 17 frictional factors on the contact surface was not considered.These factors may diminish the accuracy of the model; further research is required to address these limitations.

Figure 7 .
Figure 7.Comparison of theoretical model and experimental data.Figure 7. Comparison of theoretical model and experimental data.

Figure 7 . 17 Figure 8 .
Figure 7.Comparison of theoretical model and experimental data.Figure 7. Comparison of theoretical model and experimental data.

Figure 8 .
Figure 8.Comparison of theoretical model and existing models.

Table 1 .
Parameter values of rough surfaces.

Table 1 .
Parameter values of rough surfaces.