FFT-Based Numerical Method for Nonlinear Elastic Contact

In theoretical research pertaining to sealing, a contact model must be used to obtain the leakage channel. However, for elastoplastic contact, current numerical methods require a long calculation time. Hyperelastic contact is typically simplified to a linear elastic contact problem, which must be improved in terms of calculation accuracy. Based on the fast Fourier transform, a numerical method suitable for elastoplastic and hyperelastic frictionless contact that can be used for solving two-dimensional and three-dimensional (3D) contact problems is proposed herein. The nonlinear elastic contact problem is converted into a linear elastic contact problem considering residual deformation (or the equivalent residual deformation). Results from numerical simulations for elastic, elastoplastic, and hyperelastic contact between a hemisphere and a rigid plane are compared with those obtained using the finite element method to verify the accuracy of the numerical method. Compared with the existing elastoplastic contact numerical methods, the proposed method achieves a higher calculation efficiency while ensuring a certain calculation accuracy (i.e., the pressure error does not exceed 15%, whereas the calculation time does not exceed 10 min in a 64 × 64 grid). For hyperelastic contact, the proposed method reduces the dependence of the approximation result on the load, as in a linear elastic approximation. Finally, using the sealing application as an example, the contact and leakage rates between complicated 3D rough surfaces are calculated. Despite a certain error, the simplified numerical method yields a better approximation result than the linear elastic contact approximation. Additionally, the result can be used as fast solutions in engineering applications.

In theoretical research pertaining to sealing, the contact pressure and the area between the asperities of a sealing pair must be calculated to obtain the leakage rate, which is the most effective parameter for characterizing the sealing performance of sealing systems [1][2][3][4][5].For axisymmetric sealing structures, the contact between sealing pairs is typically simplified to a two-dimensional (2D) problem.However, in specific cases, such as those pertaining to the anisotropy of the sealing surface texture, the contact problem is difficult to simplify and must solved using a three-dimensional (3D) model.Metal/rubber and metal/plastic pairs are typically used in sealing systems.The deformation of rubber or plastic is primarily considered because of their smaller elastic modulus than that of metal.Unlike linear elastic contact, rubber is a hyperelastic material whose uniaxial compressive stressstrain curve is nonlinear.Plastic deformation occurs when the contact pressure is extremely high.Therefore, it is necessary to consider whether the contact problem is a 2D or 3D problem, and whether the contact involved is linear or nonlinear elastic contact for different sealing systems; subsequently, an appropriate contact model should be selected to perform calculation.
Existing contact models can be classified into three categories: statistical models (such as the GW model [6] and improved models based on the GW model [7,8]), fractal models [9,10] and deterministic models.For statistical and fractal models, the description of a rough surface is simplified using statistical and fractal parameters, respectively.In addition, the calculation for both models are performed by assuming a simplified interaction between asperities, or by disregarding it altogether.Compared with the first two models mentioned above, the deterministic model is based on the deterministic description of a rough surface.More accurate pressure distributions and deformation results can be obtained using this model.Deterministic models can be classified into analytical and numerical models.For the analytical model, the Hertz theory can be used to solve the contact area and average contact pressure of the ball-plane elastic contact.Westerggard [11] derived a complete solution in a closed form for the elastic contact of a one-dimensional (1D) sinusoidal surface with a flat surface.Subsequently, Johnson [12] extended Westerggard's results to a 2D sinusoidal surface with a flat surface.Zhu [13] proposed an elastic-plastic model for line contact structures based on the yield mechanism.In terms of numerical simulation, the finite element method (FEM) can be used to analyze the contact behavior between solids of different materials [14][15][16].However, for complicated rough surface contact, particularly in 3D problems, a refined mesh is required for the FEM; this would result in a long calculation time or non-convergence.Another approach is the semi-analytical method.Polonsky [17] proposed a numerical method to solve rough surface contact problems using multilevel multi-summation and conjugate gradient techniques.Stanley [18] solved the elastic contact between a rigid plane and a half-space using the fast Fourier transform (FFT) and analytical solutions from Westerggard [11] and Johnson [12]; it is noteworthy that the interaction between asperities was considered in the abovementioned method.Similarly, the FFT was used in the numerical methods proposed by Johnson [12] and Wang [19].In addition to the linear elastic contact calculation, Jacq [20], Wang [21], Chen [22], and other researchers [23,24] proposed numerical methods for elastoplastic contact.Moreover, other contact models that consider friction (tangential deformation) [25,26] and heat [27,28] have been developed.
For elastoplastic contact, the current semi-analytical methods [20][21][22] are primarily based on purely elastic contact considering plastic deformation, in which plastic strain is calculated using the von Mises yield criterion and a specific hardening law.This method involves the calculation of the stress and strain of the bulk element, which can be used to analyze contact fatigue, damage, etc.However, for complicated 3D rough surface contact, a large number of 3D grid elements would result in a long calculation time.Therefore, a simplified method adapted to engineering applications must be devised.In addition, the hyperelastic contact problem is typically simplified to a linear elastic contact problem, which must be improved in terms of calculation accuracy.
Based on the study of Stanley et al. [18], a numerical method suitable for elastoplastic and hyperelastic frictionless contact that can be used to solve 2D and 3D contact problems is proposed herein.The compressive stress-strain curve of a material with linear hardening or a hyperelastic material such as rubber is simplified to a combination of two linear segments with different slopes, and the solution of the residual deformation is simplified as well.The nonlinear elastic contact problem is converted into a linear elastic contact problem considering residual deformation (or the equivalent residual deformation).Results from numerical simulations for elastic, elastoplastic, and hyperelastic contact between a hemisphere and a rigid plane are compared with those obtained using the FEM to verify the accuracy of the numerical method.For nonlinear elastic contact, the numerical method yields a better approximation result than the linear elastic contact approximation.Finally, a numerical method is used to calculate the contact between complicated 3D rough surfaces.

Existing Numerical Method for Linear Elastic Contact
Based on a study by Stanley [18], an improved method to calculate elastoplastic and hyperelastic contact is proposed herein.To facilitate the understanding of the following sections, the main idea of Stanley's study is briefly introduced in this section.Readers can perform detailed calculations by referring to the original literature.
Based on the results of Westerggard [11] and Johnson [12], when a 1D (or 2D) sinusoidal surface and a rigid plane are in complete contact, a linear relationship exists between the sinusoidal pressure variation and displacement.For any dimensionless pressure distribution p * =p/E * (where p is the actual pressure, and E * is the equivalent elastic modulus), p * can be expressed as a superposition of the trigonometric series.Therefore, the expression between deformation u and p * for a discrete representation can be expressed using the FFT, as follows: (1) u p * = FFT −1 w. × FFT p * .
For 2D contact problems, 1D FFT and 1D inverse fast Fourier transform (FFT −1 ) are used; both p * and the numerical factor w are vectors.For 3D contact problems, 2D FFT and 2D FFT −1 are used, and both p * and w are matrices; w can be derived from the analytical formulations proposed by Westerggard and Johnson.The final problem is to obtain p * and u when the two surfaces are pressed together.Based on minimizing the total complementary energy, the final pressure distribution p * and u must satisfy Eq. ( 2) in the constraint region expressed in Eq. ( 3), where f is the total complementary energy, g i the gap between the rigid plane and the undeformed surface, n the number of discrete points, and p target the average pressure.Based on the relationship between u and p * shown in Eq. ( 1), the problem involving two variables shown in Eq. ( 2) is converted into a problem with only one variable: In Stanley's study, the gradient descent method was used to iteratively solve the distributions of p * and u.The details are available in Ref. [18].

Solution of Residual Strain for Nonlinear Elastic Contact
Figure 1 shows the uniaxial compressive stress-strain curve of a material with linear hardening.The slopes of the two linear segments are E and E' (E > E').For a unit in a uniaxial compression state, when the stress is greater than the initial yield stress (σ n ), the unit undergoes plastic deformation.The corresponding strains of the two linear (2) segments are ε and ε', respectively.The elastic strain is ε e , and the plastic strain is ε r (residual strain).Therefore, For units with multidirectional stresses, σ refers to the von Mises stress of the unit.As shown in Eq. ( 4), for a unit that undergoes plastic deformation (σ > σ n , p > p n , p n is the critical yield pressure), the strain solution under stress σ is equivalent to applying residual strain ε r to the unit, and the linear elastic contact (elastic modulus E) is used to solve the strain under stress σ.For the elastoplastic contact of the entire rough surface, the solutions for the pressure and deformation can be equivalent to the solution for a linear elastic contact problem considering residual deformation.
The uniaxial compressive stress-strain curve of a hyperelastic material such as rubber can be simplified to a combination of two linear segments with different slopes, as shown in Figure 2 (E < E').The corresponding strains of the two linear segments are ε and ε', respectively.Similar to the elastoplastic contact problem, the strain of the unit is equal to the sum of the equivalent elastic strain ε e and the equivalent residual strain ε er (ε er is negative), as shown in Eq. ( 5).For the hyperelastic contact of the entire rough surface, the solutions for the pressure and deformation can be equivalent to the solution for the linear elastic contact problem (elastic modulus E) considering the equivalent residual deformation.

Proposed Numerical Method
Based on the analysis above, it is clear that in terms of the solution for residual deformation, if the interaction between units is disregarded, then all units under pressure can be simplified to a uniaxial compression state. (4) (5) ε e +ε er = ε + ε ′ .Therefore, the residual deformation solution is simplified.Both the elastoplastic contact and hyperelastic contact problems can be equivalent to the linear elastic contact problem considering residual deformation (or the equivalent residual deformation), i.e., the piecewise linear elastic contact problem; as such, the solution algorithms are similar, as shown in Figure 3.For the solution of the elastoplastic contact, the pressure distribution p and deformation u e can be obtained by calculating the linear elastic contact using the elastic modulus E. According to the critical pressure p n , the p at each unit is decomposed into p e and p p , as shown in Eq. ( 6).For p e and p p , linear elastic contact with elastic moduli E and E' are calculated, respectively.Subsequently, the corresponding deformations u and u' and the residual deformation u r2 are calculated, as shown in Eq. ( 7).Next, it is determined whether the residual deformation converges with max|u r2 −u r1 | / max|u r2 | ≤ er (where er is the specification error) or u r2 =0 (no plastic deformation in this case).If convergence does not occur, then the residual deformation u r1 is modified to u r1 = u r1 + k 0 (u r2 − u r1 ) [20], and the surface topography z is updated.The process is repeated until the residual deformation converges.The solution process for the hyperelastic contact is similar.( 6) )

Comparison with Existing Method
Using Jacq's study [20] as an example of existing numerical methods, the simplified method for elastoplastic contact proposed herein is compared with the existing method, as shown in Figure 4.The main difference between the two methods is the method of solving the residual displacement (step 3).In the existing method, the contact pressure distribution is first solved via a linear elastic contact calculation (step 1).Subsequently, the plastic strain distribution of the bulk elements is solved using the contact pressure distribution (step 2).Finally, the residual deformation is obtained using the plastic strain distribution (step 3).The solution for the plastic strain (step 2) involves an iterative process: the pressure stress distribution of bulk elements is solved based on the contact pressure distribution; subsequently, the von Mises stress distribution is obtained.The plastic strain distribution of bulk elements is calculated using the linear hardening law, and the residual stress distribution is obtained.A new stress distribution for the bulk elements is obtained by the superposition of the pressure and residual stresses.Subsequently, the von Mises stresses and plastic strains are solved again.The loop (P-loop) is repeated until the plastic strain converges.For the simplified method, the solution for the residual deformation can be simplified (steps 6 to 3) using Eqs.( 6) and (7) when the interaction between units is disregarded.
A comparison of the calculation efficiency is shown in Tables 1, 2. For the simplified method, the average loop number of step 1 (E-loop) was approximately 10, which is similar to Jacq's method.Because the iterative solution of plastic strain was not involved, it can be assumed that the number of P-loops was 1. Furthermore, the total P-loop number was only 10 for each loading step, which is significantly less than 300 in Jacq's method.In addition, as shown in Table 2, the effort to perform one E-loop of the two methods was similar, where m is the iteration number for one E-loop, and N s is the number of discrete points on the rough surface.For the simplified method, the effort to perform one P-loop increases as 2N s lnN s .However, Jacq's method requires the calculations of pressure and residual stresses; therefore, the effort to perform one P-loop increases as N vz (N s +N vs )ln(N s +N vs ) and N vz 2 N vs ln N vs (where N vz is the number of points in the plastic volume at a certain depth, and N vs is the number of points in the plastic volume per depth), which increases the calculation time.In summary, both the number of iterations and the calculation effort of each iteration of the simplified method are significantly less than those of Jacq's method; therefore, the calculation efficiency can be improved significantly.
However, because the calculation of the residual displacement is based on a certain simplification, the calculation accuracy is reduced to a certain extent, whereas the calculation efficiency is improved.For an initial average pressure distribution p 0 , the pressure distribution p can be calculated via linear elastic contact calculation, and the actual stress state of any element is shown in Figure 5(a); the von Mises stress of the element is σ VM1 .Based on Eq. ( 6), p is decomposed into p e and p p .The corresponding stress states of the element under p e and p p are shown in Figure 5(b) and (c), respectively; the corresponding von Mises stresses are σ VM2 and σ VM3 , respectively.In fact, if the interaction between the units is considered, the pressure decomposition based on Eq. ( 6) will result in a von Mises stress error Δe1 and a residual displacement error Δe2, as shown in Eqs. ( 8) and ( 9), respectively.The value of the error is associated with the actual stress state of the elements.

Verification
In this study, the contact pressure and contact area of the elastic contact, elastoplastic contact, and hyperelastic contact between the hemisphere and rigid plane were calculated using the numerical method.The results were compared with those obtained using the FEM and the Hertz analytical solutions to verify the accuracy of the algorithm.The radius (R) of the hemisphere was 0.5 mm.

Linear Elastic Contact
According to the Hertz theory, when a load F = 4p 0 ER 2 = p 0 EA 0 is applied between a ball of radius R and a rigid plane, the equivalent elastic modulus E * , maximum deformation of the ball ω, contact area A, and maximum contact pressure P max can be expressed as shown in Eq. ( 10) [29]: Figure 6(a) shows the results of elastic contact between the hemisphere (E = 210 GPa; Poisson's ratio ν = 0.3) and the rigid plane calculated using the Hertz theory and the FEM (using ABAQUS software), where P max /E--FEM, P max /E--Hertz, CA--FEM, and CA--Hertz represent P max /E and CA calculated using FEM and Hertz theory, respectively; error--P max and error--CA are relative error of two methods.As shown, as the load increased, the dimensionless contact area CA = A/A 0 , and the maximum contact pressure P max increased gradually.When the load was small, the result calculated using the Hertz theory was similar to that obtained using the FEM, thereby verifying the accuracy of the FEM calculation (8) .  within the range of small loads.Moreover, the relative error between the two methods increased gradually as the load increased.The contact area calculated using the Hertz theory was larger than that obtained using the FEM, whereas the maximum contact pressure result was the opposite.Based on the derivation process of the Hertz theory, it was discovered that the calculated contact area was larger than the actual value, whereas the maximum contact pressure was smaller than actual value, which is consistent with the results above.In addition, the Hertz theory is only suitable for small loads and deformations.When the deformation is large, the calculation error becomes extremely large and hence the Hertz theory cannot be applied.For example, when the dimensionless load p 0 = F/A 0 E = 0.0003, the maximum deformation of the ball is ω≈0.0042mm = 0.0084R, the radius of contact area is approximately 0.0465 mm = 0.093R, and ω 2 /a 2 ≈ 0.008 ≪ 1; therefore, the calculation error is negligible, and the Hertz theory is applicable.When p 0 = 0.1, ω ≈ 0.16 mm = 0.32R, a ≈ 0.3 mm = 0.6R, and ω 2 /a 2 ≈ 0.28, the calculation error is not negligible.To further verify the accuracy of the FEM under a large load (p 0 = 0.1), the mesh was refined from the original 17908 elements to 60282.It was discovered that the maximum contact pressure increased from 0.5298E to 0.5319E, and the contact area increased from 0.2884A 0 to 0.2896A 0 , and the relative errors for both were less than 0.5%.Therefore, the calculated values using the FEM were regarded as accurate for the entire load range.Because Stanley's method [18] is applicable to complete contact, w must be modified to w' = k × w, based on the FEM results.The modified results are presented in Figure 6(b), where k is correction coefficient; error--P max and error--CA are relative errors of numerical method and FEM, respectively.As shown, the value of k remained stable over the entire load range of approximately 0.135.Furthermore, elastic contact with an elastic modulus of 0.5E and 2E was obtained using the FEM.The results show that, regardless of the elastic modulus, the calculation results under the same p 0 were the same, indicating that the coefficient k is applicable to the linear elastic contact calculation of materials with different elastic moduli.The CA and P max were calculated using the numerical method, and their relative errors compared with the FEM results under different loads p 0 were calculated as well, as shown in Figure 6(b).In the entire load range, the relative error of P max increased with the load, and all error values were less than 10%.The relative error of the CA first decreased and then increased as the load increased.For the numerical simulation, the CA is the ratio of the number of elements in the contact area (approximately circular) to the total number of elements.Compared with a large load, when the load is smaller, the same error for the number of elements in the contact area will result in a larger relative error for the CA.In addition, the total number of elements (64 × 64 in this section) and the presence of a planar area around the hemisphere (as shown in the inset of Figure 6(b)) will affect the relative error of the CA.

Elastoplastic Contact
In this study, the numerical method was used to calculate the elastoplastic contact between a hemisphere and a rigid plane.The results were compared with the linear elastic and elastoplastic results calculated using the FEM, as shown in Figure 7.In Figure 7, P max,E,FEM , P max,EP,FEM , and P max,EP,S represent P max /E of elastic contact by FEM, elastoplastic contact by FEM, and elastoplastic contact by numerical method.Representation for CA is denoted similarly.P max,S,FEM , CA S,FEM are relative error of elastoplastic contact for numerical method and FEM, respectively.The uniaxial compressive stress-strain curve of the material exhibits linear hardening, an elastic modulus E = 210 GPa, E′ = 0.5E = 105 GPa, and an initial yield stress σ n = 235 MPa = 0.0011E.The yield stress σ Y and Figure 6 Results of elastic contact plastic strain ε p of the material can be expressed using Eq. ( 11) [21], and the material properties for the FEM can be defined accordingly.Based on the CEB model [29], an initial yield occurs when p > KH = K•a•σ n = 0.6 × 2.8 × 235 MPa = 394.8MPa = 0.0019E.Therefore, Eq. ( 12) was used to define the critical pressure p n in the numerical method.If the von Mises stress is greater than the yield stress σ n in the FEM, or the maximum contact pressure P max is greater than the critical yield pressure p n in the numerical method, then the material will yield: As shown in Figure 7, if plastic deformation is considered, then P max decreases, and CA increases under the same load.Compared with the linear elastic contact results, the elastoplastic contact results of P max and CA calculated using the numerical method were similar to those of the FEM results.Yielding occurred over the entire load range, and the plastic deformation area increased with the load.The error of the elastoplastic contact caused by the numerical method was due to the calculation of the linear elastic contact and the solution for the residual deformation.When the load p 0 ranged between 0.0003 and 0.1, the relative error of the contact area was small, and the relative error of the maximum contact pressure was less than 15%.
In the numerical calculation for the elastoplastic contact, the total number of elements was 64 × 64 = 4096, and the specification error er (as shown in Figure 3) was 0.03.To understand the effects of the number of grids and the value of er on the calculation results, the calculation time and calculation error of P max for different er values (0.01, 0.03, 0.05, 0.1) and different grid numbers (11) (16 × 16, 32 × 32, 64 × 64, 128 × 128) were analyzed, as presented in Figure 8.As shown in Figure 8(a), within the load range, the calculation time was 1-9 min, and the calculation error ranged between 9.5% and 14.5%, and both decreased as the load increased.At the same grid number and same load, the calculation time decreased slightly, whereas the calculation error increased slightly as er increased.In terms of the effect of the grid number, as shown in Figure 8(b), at the same er and the same load, the calculation time increased with the grid number, and the time span was large.The calculation time did not exceed 1 min when the grid number was 16 × 16.When the grid number was 32 × 32, the calculation time was 0-2.5 min.When the grid number was 64 × 64, the calculation time was 2-8 min.When the grid number was 128 × 128, the calculation time was 8-47 min.However, except for the larger calculation error with a grid number of 16 × 16, the calculation errors corresponding to the other three sets of grid numbers were similar.Therefore, considering the calculation time and calculation error, a grid number of 32 × 32 or 64 × 64 appeared to be a better option.

Hyperelastic Contact
The contact behavior between the hyperelastic hemisphere and the rigid plane was calculated.In the FEM, the constitutive model of rubber must be defined based on the compressive stress-strain curve of the material, as presented in Figure 9(a).As shown, among the five fitted constitutive models, the Ogden 3 model yielded the best fitting result; therefore, it was selected to define the material properties in the FEM.For numerical simulation, the original stress-strain curve must be simplified to a combination comprising two linear segments.As shown in Figure 9(b), three simplified results were obtained: FC1, FC2, and FC3.The parameters of each fitted curve are shown in Table .As shown, four values were obtained for elastic moduli E and E′, whereas three values were obtained for the critical stress.
The hyperelastic contact calculation under different loads (0.001-50 N) was performed using the FEM.As shown in Figure 9(b), two elastic moduli, E 1 = 3.1251 MPa and E 2 = 32.4054MPa, were selected as the elastic modulus of the linear elastic contact to approximate the hyperelastic contact (ν = 0.49), as presented in Figure 10.In Figure 10, Elastic E 1 and Elastic E 2 represent linear elastic contact results using FEM with elastic modulus E 1 and E 2 , respectively.HE-FEM represents results of hyperelastic contact calculation using FEM.As shown, the approximation results can be artificially segmented into three regions.When the load was small, the linear elastic contact result with a lower elastic modulus E 1 was similar to the hyperelastic contact result (area A).When the load was large, the linear elastic contact calculation with a higher elastic modulus E 2 yielded a better approximation result (area C).Between the two regions (area B), the approximation errors of the linear elastic contact with both elastic moduli were relatively large.For the approximation via linear elastic contact calculation, the selection of the elastic modulus E and the approximation result are associated with the load and material properties, which are yet to be elucidated currently, particularly for complicated rough surfaces.Therefore, the numerical method used in this study was utilized to approximate the hyperelastic contact.
In reference to Figure 10, load values of 0.01, 1, and 5 N were selected for areas A, B, and C, respectively.For the three fitted curves FC1, FC2, and FC3, nonlinear elastic contact calculation was performed using the numerical method under the selected three loads (the critical pressure p n is calculated using Eq. ( 12)).The results of the pressure distribution along the radius of the contact area are shown in Figure 11.In Figure 11, "Hyperelastic--FEM" represents result of hyperelastic contact calculation in FEM."FC1, " "FC2, " "FC3" represent results of non-linear elastic contact calculation with fitted curves FC1, FC2, and FC3 using numerical method, respectively."Elastic-E1-FEM" and "Elastic-E2-FEM" represent results of elastic contact calculation with elastic moduli E 1 and E 2 using FEM, respectively.When the load was 0.01 N (Figure 11(a)), the maximum contact pressures of the three fitted curves calculated using the numerical method were smaller than the corresponding critical pressures; therefore, they were in the state of linear elastic contact.Curve 2 corresponds to a linear elastic contact approximation with an elastic modulus of 0.6027 MPa, whereas curves 3, 4, and 5 were equivalent to a linear Figure 9 Fitting results of compressive stress-strain curve of material elastic contact approximation with an elastic modulus of 3.1251 MPa; curve 6 corresponds to a linear elastic contact approximation with an elastic modulus of 32.4054 MPa.It was observed that under a small load, linear elastic contact with a lower elastic modulus E 1 and nonlinear elastic contact of the three fitted curves yielded good approximation results.When the load was 1 N (Figure 11(b)) or 5 N (Figure 11(c)), the three fitted curves were in the state of nonlinear elastic contact.When the load was 1 N, compared with the large approximation error of linear elastic contact with either a lower elastic modulus E 1 or a higher elastic modulus E 2 , the three nonlinear elastic contact calculations yielded better approximation results.When the load was 5 N (5 N was extremely large for linear elastic contact with elastic modulus E 1 to converge in the FEM; therefore, it was not considered), the two nonlinear elastic contacts of fitted curves FC2 and FC3 fitted better.In conclusion, compared with the approximation by linear elastic contact with a specific elastic modulus, the nonlinear elastic contact approximation offers calculation results of hyperelastic contact that can be better fitted within the load range in this study.Considering the calculation efficiency and approximation results, FC3 was the best among the three fitted curves.The calculation error of this numerical method originated primarily from the fitting results of the stress-strain curve, linear elastic contact calculation, and the solution of the equivalent residual deformation.In theory, to better fit the calculation results of hyperelastic contact, the original stress-strain curve can be simplified to a combination of multiple linear segments with a shape that resembles the original curve.

Machined Rough Surface
For the contact between the machined rough surface and the rigid plane (Figure 12(a)), the typical form of the cross section of the leak path formed at the interface contact was discovered to be a triangle with α≈4° [30].As shown in Figure 12(b), the height of the initial undeformed cross section of the leak path was h 0 .Under an applied load F = p 0 E 0 × l 0 , the contact area between the rough surface and the rigid plane was A c .The cross section of the leak path was reduced proportionally, and the deformed height was reduced to h.Therefore, the relationship between h and A c can be expressed as shown in Eq. ( 13), where CA denotes the dimensionless contact area.And the topography and pressure distribution of deformed surface can be obtained by contact calculations, shown in Figure 11 and Figure 13.Using the laminar flow of an incompressible fluid as an example, the viscosity of the fluid is η.When the fluid pressure difference across the leak path is Δp, the leak rate of a single deformed leak path can be calculated using Eq. ( 15) [31], where q V is the volume flow rate: (13)  As shown in Eqs. ( 9) and (10), to obtain the volume flow rate, the contact area A c or dimensionless contact area CA must be solved based on the applied load.Therefore, the numerical method was used to calculate the linear elastic, elastoplastic, and hyperelastic contact between the machined rough surface (the relevant parameters are shown in Table 3) and the rigid plane.In this example, the height of the initial undeformed cross section of leak path h 0 was 2 μm, the total number of elements was 64 × 64 = 4096, and the specification error er was 0.03.The calculation results are presented in Figure 14.As shown, as the load increased, the contact area increased gradually.For the elastoplastic contact, under the same load, the contact area was larger than that of the elastic contact and increased as E′ decreased.For the hyperelastic contact, compared with the elastic contact, the contact area reduced and increased as E' decreased.
As presented in Eq. ( 15), q V is proportional to h 4 when Δp, l 0 , and η remain constant; therefore, h 4 was calculated to compare the relative value of the leakage rate, as shown in Figure 15.Contrary to the calculation result of the contact area, q V reduced, and the difference in q V for different materials increased with the load.For the elastoplastic contact, under the same load, the flow rate was less than that of the elastic contact and decreased with E'.For the hyperelastic contact, compared with the elastic contact, the flow rate increased and then decreased as E' decreased.

Random Rough Surface
For an isotropic complicated rough surface, if the surface topography is known, then the leakage channel can be determined using the lattice leakage model [32], and the leakage rate can be calculated using Eq. ( 14).The topography of the undeformed rough surface is shown in Figure 16, where the length and width of the rough surface (14) are Δx and Δy, respectively.The undeformed profile is expressed as shown in Eq. ( 16), and the units of z was micrometers.
The contact area and leakage rate were calculated, and the results are shown in Figures 17, 18, respectively.In this example, the total number of elements was 64 × 64 = 4096, and the specification error er was 0.03.Similar to the calculation results presented in Section 3.2.1, for the same material, the contact area increased, and the leakage rate decreased as the load increased.The effects of different materials on the calculation results was the same as that presented in (16) z = cos( πx 8 ) cos( πy 8 ).   b3)), blue indicates the contact field, green the through leak channel, and white the leak channel with no fluid passing through.When the applied load p 0 was in the range of 0.5-2, the dimensionless contact area was less than 0.35, and a through leakage channel was present.When p 0 exceeded 3, the contact area was extremely large, which hindered the formation of a through leakage channel; therefore, the leakage rate was 0.

Conclusions
(1) Based on the FFT, a numerical method suitable for elastoplastic and hyperelastic frictionless contact that can be used to solve 2D and 3D contact problems was proposed herein.(2) For elastoplastic contact, the calculation efficiency improved at the expense of calculation accuracy when using the proposed method.(3) For hyperelastic contact, the proposed method reduced the dependence of the approximation result on the load used in the linear elastic approximation.Therefore, the elastic modulus need not be selected based on the load, material properties, and other factors, as in the linear elastic contact approximation.(4) In some engineering applications, contact calculation between complicated rough surfaces and fluid-solid coupling is necessitated.In most processing methods, the linear elastic contact result is used as the approximation result to reduce the calculation time, which must be improved to achieve better calculation accuracy.The proposed method provides a balance for such problems; therefore, it is expected to be advantageous in engineering applications.(5) To better fit the calculation results of the hyperelastic contact, the original stress-strain curve can be simplified to a combination of multiple linear segments, and the approximation results can be further discussed.

Figure 1 Figure 2
Figure 1 Uniaxial compressive stress-strain curve with linear hardening

Figure 3 Figure 4
Figure 3 Algorithm of nonlinear elastic contact

Figure 7 Figure 8
Figure 7Results of elastoplastic contact obtained using numerical method

Figure 10 Figure 11
Figure 10Results of hyperelastic contact and approximation via linear elastic contact calculation using FEM

Figure 12 Figure 13
Figure 12 Contact between machined rough surface and rigid plane

Figure 14 Figure 15
Figure 14 Dimensionless contact area CA for different materials

Figure 16 Figure 17 Figure 18
Figure 16 Topography of undeformed surface

Table 1
Iteration number comparison between Jacq's method and simplified method

Table 2
Calculation effort comparison between Jacq's method and simplified method One E-loop (step 1) One P-loop (step 2) vz (N s +N vs )ln(N s +N vs ) N vz 2 N vs ln N vs Simplified method mN s lnN s 2N s lnN s (step 6) Figure 5 Stress state of element under different pressure distributions

Table 3
Relevant parameters of different materials