Research on Elastic–Plastic Contact Behavior of Hemisphere Flattened by a Rigid Flat

The contact behavior of a hemisphere pressed by a rigid plane is of great significance to the study of friction, wear, and conduction between two rough surfaces. A flattening contact behavior of an elastic–perfectly plastic hemisphere pressed by a rigid flat is researched by using the finite element method in this paper. This behavior, influenced by different elastic moduli, Poisson’s ratios, and yield strengths, is compared and analyzed in a large range of interference values, which have not been considered by previous models. The boundaries of purely elastic, elastic–plastic, and fully plastic deformation regions are given according to the interference, maximum mean contact pressure, Poisson’s ratio, and elastic modulus to yield strength ratio. Then, a new elastic–plastic constitutive model is proposed to predict the contact area and load in the elastic–plastic range. Compared with previous models and experiments, the rationality of the present model is verified. The study can be applied directly to the contact between a single sphere and a plane. In addition, the sphere contact can also be used to simulate the contact of single asperity on rough surfaces, so the present proposed model can be used to further study the contact characteristics of rough surfaces.


Introduction
The elastic-plastic contact behavior between hemisphere and rigid plane is one of the fundamental problems in particle mechanics [1,2], contact mechanics [3][4][5][6][7][8][9], and biomechanics [10]. Mechanical surfaces are microscopically rough. Contact between two rough surfaces can be equivalent to contact between a hemisphere and a rigid flat [11][12][13]. The study of elastic-plastic contact behavior between hemisphere and a rigid flat is of great significance to the analysis of electrical contact [14][15][16][17], friction [18,19], and wear [20,21] between two rough surfaces. The early contact model only predicted the elastic contact behavior between the hemisphere and the rigid flat. Greenwood and Williamson (GW model) [22] analyzed the pure elastic contact behavior between hemisphere and rigid flat. With the increase of interference, the hemispheres yield, and the contact state changes from pure elastic to elastic-plastic. Based on the Hertz solution [23], the GW model presented the critical normal interference (ω c ), critical contact load (F c ), and critical contact area (A c ) formulas for the initial yield of the hemisphere. However, the Poisson's ratio (ν) effect is ignored in their predicted formulas for these critical contact parameters. Later, Lin and Lin (LL model) [24] and Chang et al. (CEB model) [25,26] improved the GW model and further considered the effect of different ν on these critical contact parameters. In the elastic-plastic contact stage, the contact behavior is highly nonlinear due to the existence of plasticity. Some scholars have predicted the elastic-plastic contact behavior between hemispheres and rigid planes through experiments [27,28]. However, the experimental results only apply to specific A flattening model of a deformable hemisphere pressed by a rigid flat before and after loading is shown in Figure 1. A uniform downward load (F) is applied to the top of the hemispheres to simulate loading. At a small interference (ω), the contact state is purely elastic. With an increase of interference, the initial yield occurs at the contact subsurface depth z of the hemisphere, marking the end of purely elastic deformation or the inception of elastic-plastic deformation.
Due to the complexity of elastic-plastic problems, there is still no closed solution. Although some hemispheric contact models exist, there is still a lack of a model to compare and analyze hemispheric elastic-plastic contact behavior over different E, ν, Y, and large ω/ωc values. In this paper, the elastic-plastic contact behavior of the hemisphere pressed by the rigid flat by E, ν, and Y was analyzed based on FEM. According to ω/ωc, maximum p/Y, ν, and E/Y, the boundaries of elastic, elastic-plastic, and fully plastic deformation regimes were given. A new elastic-plastic constitutive model was proposed to predict the A/Ac and F/Fc in the elastic-plastic range. Compared with the results predicted by previous models and experiments, the rationality of the present model is verified. In Section 2, the critical contact parameters for the initial yield of hemispheres are analyzed. The FE model is presented in Section 3. In Section 4, the empirical formula of hemispheric elastic-plastic contact is demonstrated by curve fitting, and the results are compared with previous models. The conclusions are presented in the last section.

Critical Formula
A flattening model of a deformable hemisphere pressed by a rigid flat before and after loading is shown in Figure 1. A uniform downward load (F) is applied to the top of the hemispheres to simulate loading. At a small interference (ω), the contact state is purely elastic. With an increase of interference, the initial yield occurs at the contact subsurface depth z of the hemisphere, marking the end of purely elastic deformation or the inception of elastic-plastic deformation. The interference (ω) at the initial yield is called the critical interference, ωc, which is calculated according to the formula provided by Johnson [3] and given by  The interference (ω) at the initial yield is called the critical interference, ω c , which is calculated according to the formula provided by Johnson [3] and given by where E is the equivalent elastic modulus and R is the equivalent radius. ν 1 , ν 2 , and E 1 , E 2 are the Poisson's ratios and elastic moduli of the two materials in contact, respectively. R 1 and R 2 are the radii of the contact pairs. p 0 is maximum Hertzian pressure, which is listed in Table 1. As shown in Table 1, p 0 is related to material hardness (H), where H is equal to 2.8 Y predicted by Tabor [34] in the GW model [22], CEB model [26], and LL model [24]. However, this relationship between H and Y, recently pointed out by the JG model [33] and SM model [38], was not the constant value equal to 2.8. In the Green model [32], BK model [37], and JG model, p 0 is predicted by using the von Mises yield criterion and the stress field of Johnson, which seems more reasonable than the GW model, CEB model, and LL model. The yield strength coefficient (C) is a function of Poisson's ratio (ν). It is shown in Figure 2. The critical results of all models are the same at ν = 0. 35 where E is the equivalent elastic modulus and R is the equivalent radius. ν1, ν2, and E1, are the Poisson's ratios and elastic moduli of the two materials in contact, respectively. and R2 are the radii of the contact pairs. p0 is maximum Hertzian pressure, which is list in Table 1.
As shown in Table 1, p0 is related to material hardness (H), where H is equal to 2.8 predicted by Tabor [34] in the GW model [22], CEB model [26], and LL model [24]. Ho ever, this relationship between H and Y, recently pointed out by the JG model [33] a SM model [38], was not the constant value equal to 2.8. In the Green model [32], BK mod [37], and JG model, p0 is predicted by using the von Mises yield criterion and the stre field of Johnson, which seems more reasonable than the GW model, CEB model, and model. The yield strength coefficient (C) is a function of Poisson's ratio (ν). It is shown Figure 2. The critical results of all models are the same at ν = 0.35. For the common ma rials (0.2 ≤ ν ≤ 0.45), except for the GW model, the critical interference of the JG mod CEB model, Green model, LL model, and BK model are similar.  Considering the above discussion and the meaning of p0, the critical interference chosen in this study by the JG model and given by The contact state is purely elastic before the initial yielding of the hemispherical co tact (ω < ωc). The contact load (Fe) and the contact area (Ae) can be expressed as Considering the above discussion and the meaning of p 0 , the critical interference is chosen in this study by the JG model and given by The contact state is purely elastic before the initial yielding of the hemispherical contact (ω < ω c ). The contact load (F e ) and the contact area (A e ) can be expressed as when the interference equals the critical interference (ω = ω c ), their critical values can be expressed as These critical values predict the particular contact parameters at the inception of elastic-plastic deformation. Therefore, they are selected to nondimensionalize the results in all models.

Finite Element Model
The finite element model presently used in this paper is similar to the finite element model described by the JG model [33], SM model [38,39], and SC model [41]. The commercial program ANSYS18.2-Workbench was used for modeling and analyzing the contact behavior of a deformable hemisphere pressed by a rigid flat. Due to its axial symmetry, a 2-D model was selected. The hemisphere was modeled by a quarter of a circle (R = 1 mm), while a half-space represented the rigid flat surface. The boundary load was applied, as shown in Figure 3. The nodes of a half-space were fixed in all directions to model the rigid flat. The nodes on the symmetry axis of the hemisphere were fixed in the radial direction to model a half-sphere. The tangent modulus and the friction coefficient were assumed as zero within ANSYS. The von Mises yield criterion was used to define material yield. The present analysis covers a wide range of material properties [48]. As shown in Tables 2 and 3, the analyzed material properties were divided into two groups. Specifically, the property ranges of the 79.4 ≤ E/Y ≤ 800 and 0.2 ≤ ν ≤ 0.45 were considered. The half-sphere and half-space by PLANE 183 triangle elements were discretized in such a way that there were many elements near the contact edge, as shown in Figure 3. The mesh far from the contact edge became coarser to improve computing speed. The half-sphere and half-space total number of elements was 37,391 and 15,028, respectively. Additionally, their total number of nodes were 75,662, and 33,164, respectively. For the half-sphere, the JG model used a constant mesh of 11,101 elements, the KE model used a maximum of 2944 nodes in total, and the SM model used 9933 elements in total for their analysis. By comparison, the present model mesh was better. There are 2696 two-dimensional three-node surface contact elements, designated as CONTA172 and TARGE169 in ANSYS, to detect the contact behavior of a deformable hemisphere pressed by a rigid flat. The uniform displacement was applied to the bottom section of the hemisphere, and then the contact force was obtained by extracting the reaction force of the node at its bottom. The contact area was calculated by the contact radius, which was obtained by finding the contacting edge. The results compare well with the Hertz elastic solution at ω ≤ ω c . The error between them is less than 2%. In addition, to ensure mesh convergence, the mesh density was iterated twice until the result of each iteration changed by no more than 1%. Due to the material and geometric nonlinearity, the Augmented Lagrange algorithm was adopted, and a large deflection was activated. Additionally, multi-step loading was used to guarantee the solution convergence. The maximum number of substeps in the calculation for very large interference was 15,000. In order to obtain a generalized model, all contact parameters were in dimensionless form, i.e., F/F c , ω/ω c , p/Y, and A/A c .

Contact Load
The effects of E and ν shown in Table 2 on F/Fc as a function of ω/ωc were analyzed Here Y was a constant value equal to 0.25 GPa, while ν were 0.2, 0.3, 0.4, and 0.45. With E ranging from 45 to 200 GPa, corresponding the E/Y ratios ranging from 180 to 800, respec tively. Here only the results of ν = 0.2 and 0.45 for E/Y = 180, 280, 400, 600, and 800 were selected and plotted in Figure 4. F/Fc increases with an increasing ω/ωc. For 1< ω/ωc ≤ 110 the results indicate the similarity with the KE model, F/Fc is independent of E and ν of the material. For 110 < ω/ωc ≤ 500, F/Fc decreased with a decreasing E. Moreover, it is clear from Figure 4 that the SM model underestimates F/Fc. The reason of this discrepancy is that it is valid only for the E/Y ratio of 79.4. It is clear from Figure 5 that with the increase of ω/ωc, the present result is lower than F/Fc predicted by the SC model. This may be the SC model analyzed the effect of strain hardening of materials on the elastic-plastic contac behavior, while the strain hardening effect is ignored in the present work. As ω/ωc in creases, as shown in Figure 4, F/Fc is also affected by ν. Additionally, for the fixed E/Y, the

Apply Displacement
Face Fixed Support Symmetry axis   Table 3. Material properties in the second group.

Contact Load
The effects of E and ν shown in Table 2 on F/F c as a function of ω/ω c were analyzed. Here Y was a constant value equal to 0.25 GPa, while ν were 0.2, 0.3, 0.4, and 0.45. With E ranging from 45 to 200 GPa, corresponding the E/Y ratios ranging from 180 to 800, respectively. Here only the results of ν = 0.2 and 0.45 for E/Y = 180, 280, 400, 600, and 800 were selected and plotted in Figure 4. F/F c increases with an increasing ω/ω c . For 1< ω/ω c ≤ 110, the results indicate the similarity with the KE model, F/F c is independent of E and ν of the material. For 110 < ω/ω c ≤ 500, F/F c decreased with a decreasing E. Moreover, it is clear from Figure 4 that the SM model underestimates F/F c . The reason of this discrepancy is that it is valid only for the E/Y ratio of 79.4. It is clear from Figure 5 that with the increase of ω/ω c , the present result is lower than F/F c predicted by the SC model. This may be the SC model analyzed the effect of strain hardening of materials on the elastic-plastic contact behavior, while the strain hardening effect is ignored in the present work. As ω/ ω c increases, as shown in Figure 4, F/F c is also affected by ν. Additionally, for the fixed E/Y, the lower the ν, the higher F/F c is.  The effects of Y and ν shown in Table 3 on F/Fc as a function of ω/ωc are analyzed. Here E is a constant value equal to 200 GPa, while ν ranges from 0.2 to 0.45. With Y ranging from 1.25 to 0.25 GPa, corresponding to the E/Y ratios ranging from 160 to 800, which are all greater than 133.3, respectively. Here only the results of ν = 0.2 and 0.45 for E/Y = 160, 200, 356.6, 571.4, and 800 are selected and plotted in Figure 5. In the present interference domain, the SM model distinctly underestimates F/Fc. For 1< ω/ωc ≤ 110, the results indicate that the similarity with the KE model, F/Fc is independent of the Y and ν of the material. For 110 < ω/ωc ≤ 500, the present results show that the value of F/Fc decreases with an increasing Y. In addition, as ω/ωc increases for the material with a lower ν, the higher the F/Fc is for the fixed E/Y. In this work, the curve fitting of simulation data is carried out. For 133.3 < E/Y ≤ 800, 0.2 ≤ ν ≤ 0.45, and 1 < ω/ωc ≤ 500, the empirical expressions of F/Fc as a function of ω/ωc is presented in this paper as follows: where m, n, and q are shown in Table 4. Table 4. Parameters in Equations (9) and (10). When E/Y is less than 133.3, F/Fc as a function of ω/ωc for different Y and ν values is plotted in Figure 6. As shown in Table 3  The effects of Y and ν shown in Table 3 on F/F c as a function of ω/ω c are analyzed. Here E is a constant value equal to 200 GPa, while ν ranges from 0.2 to 0.45. With Y ranging from 1.25 to 0.25 GPa, corresponding to the E/Y ratios ranging from 160 to 800, which are all greater than 133.3, respectively. Here only the results of ν = 0.2 and 0.45 for E/Y = 160, 200, 356.6, 571.4, and 800 are selected and plotted in Figure 5. In the present interference domain, the SM model distinctly underestimates F/F c . For 1< ω/ω c ≤ 110, the results indicate that the similarity with the KE model, F/F c is independent of the Y and ν of the material. For 110 < ω/ω c ≤ 500, the present results show that the value of F/F c decreases with an increasing Y. In addition, as ω/ω c increases for the material with a lower ν, the higher the F/F c is for the fixed E/Y. In this work, the curve fitting of simulation data is carried out. For 133.3 < E/Y ≤ 800, 0.2 ≤ ν ≤ 0.45, and 1 < ω/ω c ≤ 500, the empirical expressions of F/F c as a function of ω/ω c is presented in this paper as follows: where m, n, and q are shown in Table 4.  (9) and (10).
When E/Y is less than 133.3, F/F c as a function of ω/ω c for different Y and ν values is plotted in Figure 6. As shown in Table 3

Contact Area
The effects of E and ν shown in Table 2 on A/Ac as a function of ω/ωc are analyzed. Here Y is a constant value equal to 0.25 GPa, with ν ranging from 0.2 to 0.45. With E ranging from 45 to 200 GPa, corresponding the E/Y ratios ranging from 180 to 800, respectively. Here only the results of ν = 0.2, 0.3, 0.4, and 0.45 for E/Y = 180 and 800 are selected and plotted in Figure 7. For 1 < ω/ωc ≤ 110, the results indicate the similarity with the KE model, and A/Ac is independent of E and ν. For 110 < ω/ωc ≤ 500, A/Ac increases with an increasing ω/ωc. In addition, the results show that E and ν have little effect on A/Ac. It is clear from Figure 7 that the SM model underestimates A/Ac as ω/ωc increases, especially for large interference. The reason of this discrepancy is that the SM model is valid only for an E/Y ratio of 79.4. The present result is slightly higher than the A/Ac predicted by the SC model with an increasing ω/ωc. The SC model analyzed the effect of tangential modulus on the elastic-plastic contact behavior, while the strain hardening effect is ignored in the present study. Compared with the predicted results of the SC model, it can be observed that tangential modulus has a great effect on the elastic-plastic contact behavior. In view of the shortcomings of present study, we will consider the effect of tangential modulus on hemispheric elastic-plastic contact behavior under large interference in future work. It can be seen that F/F c is affected by ν and that the lower the ν, the higher the F/F c is with the increased ω/ω c . In this work, the curve fitting of simulation data is carried out. For 79.4 ≤ E/Y ≤ 133.3, 0.2 ≤ ν ≤ 0.45, and 1 < ω/ω c ≤ 500, the empirical expressions of F/F c as a function of ω/ω c are presented in this paper as follows: where m and n are shown in Table 5.

Contact Area
The effects of E and ν shown in Table 2 Figure 7. For 1 < ω/ω c ≤ 110, the results indicate the similarity with the KE model, and A/A c is independent of E and ν. For 110 < ω/ω c ≤ 500, A/A c increases with an increasing ω/ω c . In addition, the results show that E and ν have little effect on A/A c . It is clear from Figure 7 that the SM model underestimates A/A c as ω/ω c increases, especially for large interference. The reason of this discrepancy is that the SM model is valid only for an E/Y ratio of 79.4. The present result is slightly higher than the A/A c predicted by the SC model with an increasing ω/ ω c . The SC model analyzed the effect of tangential modulus on the elastic-plastic contact behavior, while the strain hardening effect is ignored in the present study. Compared with the predicted results of the SC model, it can be observed that tangential modulus has a great effect on the elastic-plastic contact behavior. In view of the shortcomings of present study, we will consider the effect of tangential modulus on hemispheric elastic-plastic contact behavior under large interference in future work.

Contact Area
The effects of E and ν shown in Table 2 on A/Ac as a function of ω/ωc are analyzed. Here Y is a constant value equal to 0.25 GPa, with ν ranging from 0.2 to 0.45. With E ranging from 45 to 200 GPa, corresponding the E/Y ratios ranging from 180 to 800, respectively. Here only the results of ν = 0.2, 0.3, 0.4, and 0.45 for E/Y = 180 and 800 are selected and plotted in Figure 7. For 1 < ω/ωc ≤ 110, the results indicate the similarity with the KE model, and A/Ac is independent of E and ν. For 110 < ω/ωc ≤ 500, A/Ac increases with an increasing ω/ωc. In addition, the results show that E and ν have little effect on A/Ac. It is clear from Figure 7 that the SM model underestimates A/Ac as ω/ωc increases, especially for large interference. The reason of this discrepancy is that the SM model is valid only for an E/Y ratio of 79.4. The present result is slightly higher than the A/Ac predicted by the SC model with an increasing ω/ωc. The SC model analyzed the effect of tangential modulus on the elastic-plastic contact behavior, while the strain hardening effect is ignored in the present study. Compared with the predicted results of the SC model, it can be observed that tangential modulus has a great effect on the elastic-plastic contact behavior. In view of the shortcomings of present study, we will consider the effect of tangential modulus on hemispheric elastic-plastic contact behavior under large interference in future work.  The effects of Y and ν, as shown in Table 3, on A/A c as a function of ω/ω c are analyzed. Here E is 200 GPa, with ν ranging from 0.2 to 0.45. With Y ranging from 1.25 to 0.25 GPa, corresponding the E/Y ratios ranging from 160 to 800, respectively. Here only the results of ν = 0.2, 0.3, 0.4, and 0.45 for E/Y = 160 and 800 are selected and plotted in Figure 8. A/A c increases with an increasing ω/ω c . For 1 < ω/ω c ≤ 110, the results indicate the similarity with the KE model, A/A c is independent of Y and ν. For 110 < ω/ω c ≤ 500, Y and ν have little effect on A/A c . Since the SM model is only suitable for E/Y = 79.4, A/A c is underestimated in the present interference domain.  When the ratio of E/Y is less than 133.3, A/Ac as a function of ω/ωc for different Y and ν values is plotted in Figure 9. Here E is a constant value equal to 200 GPa, with ν ranging from 0.2 to 0.45. With Y being 2.52, 2, and 1.5 GPa, corresponding the E/Y ratio is 79.4, 100, and 133.3, respectively. ω/ωc ranges from 1 to 500. At small interferences, the dependence of A/Ac on Y and ν is weak, which is similar to the prediction of KE model. However, with an increasing ω/ωc, the influence of Y on A/Ac is strong. It is clear from Figure 9 that the present A/Ac is consistent with that predicted by the SM model for E/Y = 79.4. However, A/Ac is underestimated by the SM model when the E/Y is greater than 79.4. It can also be seen that in the present ω/ωc domain, the JG model overestimates A/Ac. At a ω/ωc value of 500, A/Ac deviation between the present work and JG model is 70.73% when E/Y is 133.3 and 25.1% when E/Y is 79.4. When E/Y is greater than 133.3, the error of the JG model is In Figures 7 and 8, the effects of E, ν, and Y on A/A c are comprehensively considered under an E/Y greater than 133.3. The results show that the E/Y ratio affects A/A c . As the ω/ω c increases, the lower the E/Y, the higher the A/A c is. ν has little effect on A/A c . In this work, the curve fitting of simulation data is carried out. For 133.3 < E/Y ≤ 800, 0.2 ≤ ν ≤ 0.45, and 1 < ω/ω c ≤ 500, the empirical expressions of A/A c as a function of ω/ω c is presented in this paper as follows: where m and n are shown in Table 6. Table 6. Parameters in Equation (12).
1 < ω/ω c ≤ 10 10 < ω/ω c ≤ 500 When the ratio of E/Y is less than 133.3, A/A c as a function of ω/ω c for different Y and ν values is plotted in Figure 9. Here E is a constant value equal to 200 GPa, with ν ranging from 0.2 to 0.45. With Y being 2.52, 2, and 1.5 GPa, corresponding the E/Y ratio is 79.4, 100, and 133.3, respectively. ω/ω c ranges from 1 to 500. At small interferences, the dependence of A/A c on Y and ν is weak, which is similar to the prediction of KE model. However, with an increasing ω/ω c , the influence of Y on A/A c is strong. It is clear from Figure 9 that the present A/A c is consistent with that predicted by the SM model for E/Y = 79.4. However, A/A c is underestimated by the SM model when the E/Y is greater than 79.4. It can also be seen that in the present ω/ω c domain, the JG model overestimates A/A c . At a ω/ω c value of 500, A/A c deviation between the present work and JG model is 70.73% when E/Y is 133.3 and 25.1% when E/Y is 79.4. When E/Y is greater than 133.3, the error of the JG model is slightly less than that of other models compared with the present predicted results. To overcome the drawbacks of previous models, the curve fitting of the simulation data is carried out in this work. For 79.4 ≤ E/Y ≤ 133.3, 0.2 ≤ ν ≤ 0.45, and 1 < ω/ω c ≤ 500, the empirical formulas for predicting the contact area are presented in this paper as follows: where m and n are shown in Table 7.

Contact Pressure
The effects of E and ν shown in Table 2 on p/Y as a function of ω/ωc are analyzed. Here, Y is a constant value equal to 0.25 GPa, with ν ranging from 0.2 to 0.45. With E ranging from 45 to 200 GPa, corresponding the E/Y ratios ranging from 180 to 800, respectively. Here, only the results of ν = 0.2 and 0.45 for E/Y = 180, 400, and 800 are selected and plotted in Figure 10. The KE model considered that p/Y reaching the peak could be marked the contact state transition from the elastic-plastic to the fully plastic and that the corresponding ω/ωc is a constant value equal to 110. Through this, the elastic-plastic range (1 ≤ ω/ωc ≤ 110) is proposed by the KE model. However, this tendency is not observed in the present work. It is obvious that ω/ωc corresponding to the peak value of p/Y is not fixed but is dependent on E and ν in all cases. It is clear from Figure 10 that the peak value of p/Y also increases with an increasing E. Additionally, the value of ν increases, the peak value of p/Y also increases under constant E/Y. By comparison, the effects of ν are less than that of E on the peak value of p/Y.
After p/Y reaches its peak value, namely the fully plastic range, it is no longer affected by ν, as shown in Figure 10. In the fully plastic range, p/Y decreases with an increasing ω/ωc. Moreover, the results also show that this trend is more evident with the lower E value, which is the same as the observational results of the SC model. Interestingly, the results obtained show that before p/Y reaches its peak in all cases, namely the elasticplastic range, the higher ν, the higher p/Y is.

Contact Pressure
The effects of E and ν shown in Table 2 on p/Y as a function of ω/ω c are analyzed. Here, Y is a constant value equal to 0.25 GPa, with ν ranging from 0.2 to 0.45. With E ranging from 45 to 200 GPa, corresponding the E/Y ratios ranging from 180 to 800, respectively. Here, only the results of ν = 0.2 and 0.45 for E/Y = 180, 400, and 800 are selected and plotted in Figure 10. The KE model considered that p/Y reaching the peak could be marked the contact state transition from the elastic-plastic to the fully plastic and that the corresponding ω/ ω c is a constant value equal to 110. Through this, the elasticplastic range (1 ≤ ω/ ω c ≤ 110) is proposed by the KE model. However, this tendency is not observed in the present work. It is obvious that ω/ ω c corresponding to the peak value of p/Y is not fixed but is dependent on E and ν in all cases. It is clear from Figure 10 that the peak value of p/Y also increases with an increasing E. Additionally, the value of ν increases, the peak value of p/Y also increases under constant E/Y. By comparison, the effects of ν are less than that of E on the peak value of p/Y. The effects of Y and ν, as shown in Table 3, on p/Y as a function of ω/ωc are analyzed. Here E is a constant value equal to 200 GPa, with ν ranging from 0.2 to 0.45. With Y ranging from 2.52 to 0.25 GPa, corresponding the E/Y ratios ranging from 79.4 to 800, respec- After p/Y reaches its peak value, namely the fully plastic range, it is no longer affected by ν, as shown in Figure 10. In the fully plastic range, p/Y decreases with an increasing ω/ ω c . Moreover, the results also show that this trend is more evident with the lower E value, which is the same as the observational results of the SC model. Interestingly, the results obtained show that before p/Y reaches its peak in all cases, namely the elastic-plastic range, the higher ν, the higher p/Y is.
The effects of Y and ν, as shown in Table 3, on p/Y as a function of ω/ ω c are analyzed. Here E is a constant value equal to 200 GPa, with ν ranging from 0.2 to 0.45. With Y ranging from 2.52 to 0.25 GPa, corresponding the E/Y ratios ranging from 79.4 to 800, respectively. Here only the results of ν = 0.2 and 0.45 for E/Y = 79.4, 133.3, 200, 356.6, and 800 are selected and plotted in Figure 11. In the elastic-plastic range, p/Y increases with an increasing ω/ ω c . Additionally, the higher ν, the higher p/Y is. p/Y decreases with an increasing ω/ ω c in the plastic range. The results show that this trend is more evident with the higher Y value. Moreover, ν has no longer has an effect on p/Y in the fully plastic range. The effects of Y and ν, as shown in Table 3, on p/Y as a function of ω/ωc are analyzed. Here E is a constant value equal to 200 GPa, with ν ranging from 0.2 to 0.45. With Y ranging from 2.52 to 0.25 GPa, corresponding the E/Y ratios ranging from 79.4 to 800, respectively. Here only the results of ν = 0.2 and 0.45 for E/Y = 79.4, 133.3, 200, 356.6, and 800 are selected and plotted in Figure 11. In the elastic-plastic range, p/Y increases with an increasing ω/ωc. Additionally, the higher ν, the higher p/Y is. p/Y decreases with an increasing ω/ωc in the plastic range. The results show that this trend is more evident with the higher Y value. Moreover, ν has no longer has an effect on p/Y in the fully plastic range.
The present work has observed that in all cases, ω/ωc corresponding to the peak value of p/Y is not limited to a specific value of 110 predicted by the KE model but depends on Y and ν. It is clear from Figure 11 that the peak value of p/Y also increases when Y decreases. Moreover, as ν increases, the peak value of p/Y also increases under constant E/Y. By comparison, the effects of ν are less than Y for the peak value of p/Y.   The present work has observed that in all cases, ω/ ω c corresponding to the peak value of p/Y is not limited to a specific value of 110 predicted by the KE model but depends on Y and ν. It is clear from Figure 11 that the peak value of p/Y also increases when Y decreases. Moreover, as ν increases, the peak value of p/Y also increases under constant E/Y. By comparison, the effects of ν are less than Y for the peak value of p/Y. Figure 12 shows the end of the elastic-plastic deformation (ω * p ) as a function of E/Y. In addition, comparisons with the KE model are also depicted. The KE model suggested that p/Y reaches the peak and marks the transition from the elastic-plastic to the fully plastic deformation regime, and the corresponding ω * p is the constant value equal to 110. The method of using the maximum value of p/Y to find the end of the elastic-plastic deformation regime can also be found in Kogut et al. [31,36]. The present results show that the value of ω * p increases with an increasing of E/Y and is not the constant value equal to 110 predicted by the KE model. For all of E/Y, the higher the value of ν, the lower value of ω * p . The present results also show that the value of ω * p . is about 110 at E/Y = 200, which is consistent with the predicted result of the KE model. Moreover, when the E/Y is less than 200, the KE model overestimates the value of ω * p . Therefore, to estimate the value of ω * p and predict the end of the elastic-plastic range, the empirical formula is proposed in this work as follows:  As mentioned above, the corresponding ωc at the inception of the elastic-plastic deformation regime can be determined by using Equation (4). Additionally, the corresponding * p  at the end of elastic-plastic deformation regime can be determined by using Equation (14). Furthermore, to estimate the contact parameters of the elastic-plastic range, the empirical formula is proposed in this work as follows:

Elastic-Plastic Range
where m, n, q, and t are shown in Table 8.  Figure 12. ω * p as a function of E/Y.
As mentioned above, the corresponding ω c at the inception of the elastic-plastic deformation regime can be determined by using Equation (4). Additionally, the corresponding ω * p at the end of elastic-plastic deformation regime can be determined by using Equation (14). Furthermore, to estimate the contact parameters of the elastic-plastic range, the empirical formula is proposed in this work as follows: where m, n, q, and t are shown in Table 8. In order to more conveniently use the elastic-plastic constitutive relation presented in this paper, a flowchart for predicting the contact parameters in the purely elastic, elasticplastic, and fully plastic ranges is shown in Figure 13. First, for the purely elastic contact stage, the contact parameters are calculated according to Equations (5) and (6). Secondly, for the elastic-plastic range, the contact parameters are calculated according to Equations (15) and (16). Additionally, and lastly, for large deformations, it includes the elastic-plastic and fully plastic range. If 0.2 ≤ ν ≤ 0.45 and 79.4 ≤ E/Y ≤ 133.3, the contact parameters are calculated according to Equations (11) and (13), otherwise according to the Equations (9)- (12). tic-plastic, and fully plastic ranges is shown in Figure 13. First, for the purely elastic contact stage, the contact parameters are calculated according to Equations (5) and (6). Secondly, for the elastic-plastic range, the contact parameters are calculated according to Equations (15) and (16). Additionally, and lastly, for large deformations, it includes the elastic-plastic and fully plastic range. If 0.2 ≤ ν ≤ 0.45 and 79.4 ≤ E/Y ≤ 133.3, the contact parameters are calculated according to Equations (11) and (13), otherwise according to the Equations (9)-(12). Figure 13. Flowchart for predicting the contact parameters.

Comparison with Experimental Results
To verify the rationality of the present model, the present predicted contact load and area are compared with Ovcharenko et al. [28] and Jamari and Schipper [27] experimental data.
Jamari and Schipper [27] experimentally measured the relationship between the contact parameters of the copper sphere and the SiC ceramic flat. In their experiment, the material parameters of the copper sphere are E = 120 GPa, R = 1.5 mm, ν = 0.35, H = 1.2 GPa, and E/Y = 280, respectively. As their material parameter is E/Y = 280, Jamari and Schipper's experimental results were compared with the present prediction formula of the Equation (12) and the results are shown in Figure 14. The present prediction results are in good agreement with the experimental ones, which verifies the accuracy of the present work.

Comparison with Experimental Results
To verify the rationality of the present model, the present predicted contact load and area are compared with Ovcharenko et al. [28] and Jamari and Schipper [27] experimental data.
Jamari and Schipper [27] experimentally measured the relationship between the contact parameters of the copper sphere and the SiC ceramic flat. In their experiment, the material parameters of the copper sphere are E = 120 GPa, R = 1.5 mm, ν = 0.35, H = 1.2 GPa, and E/Y = 280, respectively. As their material parameter is E/Y = 280, Jamari and Schipper's experimental results were compared with the present prediction formula of the Equation (12) and the results are shown in Figure 14. The present prediction results are in good agreement with the experimental ones, which verifies the accuracy of the present work. Ovcharenko et al. [28] experimentally measured the relationship between the contact load and the area of the stainless-steel and copper spheres in contact with the Sapphire flat. As shown in the Figure 15  Ovcharenko et al. [28] experimentally measured the relationship between the contact load and the area of the stainless-steel and copper spheres in contact with the Sapphire flat. As shown in the Figure 15. the present predicted results are compared with the Ovcharenko et al. experimental results for the copper and stainless-steel. Compared with the Jamari and Schipper experimental study, Ovcharenko et al. only analyzed the contact area under small loads through their experiments. It is clear from Figure 15 that the present model is in good agreement with the experimental results and the error is less than 10%, indicating that the proposed model has a high accuracy. This discrepancy may be due to the fact that the SiC ceramic flat in the experiment is not as rigid as assumed in the theoretical model, leading to some error between the experimental and theoretical results. Ovcharenko et al. [28] experimentally measured the relationship between the contact load and the area of the stainless-steel and copper spheres in contact with the Sapphire flat. As shown in the Figure 15. the present predicted results are compared with the Ovcharenko et al. experimental results for the copper and stainless-steel. Compared with the Jamari and Schipper experimental study, Ovcharenko et al. only analyzed the contact area under small loads through their experiments. It is clear from Figure 15 that the present model is in good agreement with the experimental results and the error is less than 10%, indicating that the proposed model has a high accuracy. This discrepancy may be due to the fact that the SiC ceramic flat in the experiment is not as rigid as assumed in the theoretical model, leading to some error between the experimental and theoretical results.

Conclusions
A flattening contact behavior of an elastic-perfectly plastic hemisphere against a rigid flat was researched by using FEM. The present analysis addresses the effects of Poisson's ratio, elastic modulus, and yield strength on the contact behavior. The main conclusions can be drawn as follows: (1) The present study considers a large range of dimensionless interference from 1 to 500. A new elastic-plastic constitutive model is proposed to predict the contact area

Conclusions
A flattening contact behavior of an elastic-perfectly plastic hemisphere against a rigid flat was researched by using FEM. The present analysis addresses the effects of Poisson's ratio, elastic modulus, and yield strength on the contact behavior. The main conclusions can be drawn as follows: (1) The present study considers a large range of dimensionless interference from 1 to 500.
A new elastic-plastic constitutive model is proposed to predict the contact area and load based on the curve fitting the finite analyses results. Compared with previous models and experiments, the rationality of the present model is verified. (2) p/Y is mainly affected by the ν in the elastic-plastic range, and the higher the ν, the higher p/Y. However, this influence disappears in the fully plastic range. The maximum of p/Y is not a constant value for the different E, Y, and ν. The higher E/Y, the higher the maximum of p/Y is. Moreover, the higher ν, the higher the maximum value of p/Y is when E/Y is constant. However, the effects of ν are less than that of E and Y for the maximum value of p/Y. In the plastic range, p/Y decreases with increasing interference. Additionally, the lower E/Y, the more noticeable this trend is. (3) The boundaries between the elastic, elastic-plastic, and fully plastic deformation regimes are determined according to the interference, maximum mean contact pressure, Poisson's ratio, and the elastic modulus to yield strength ratio. When the interference is small, the contact state is purely elastic, and the contact parameters can be calculated according to the Hertz formula. When the interference increases, the contact state changes from the purely elastic to the elastic-plastic. The present work shows that the JG model can more reasonably determine the inception of elasticplastic deformation regime. The end of the elastic-plastic deformation regime is defined according to the interference corresponding to the maximum contact pressure. New dimensionless constitutive relationships are proposed to predict the contact parameters in the elastic-plastic range. The value of asperity 2 The value of rigid flat c Critical values at yielding inception p Critical values at the end of the elastic-plastic range * Dimensionless e Elastic range ep Elastic-plastic range