A method for calculating the contact area of ​​three-dimensional rough surfaces

Abstract

The invention discloses a method for calculating the contact area of ​​a three-dimensional rough surface. The method combines an analytical model of MATLAB with a simulation model of finite element software COMSOL. First, a three-dimensional anisotropic rough surface is simulated and generated in MATLAB software based on fractal theory, and contact behaviors of elastic-plastic rough surfaces with different material characteristic parameters and surface morphology characteristics are numerically simulated and calculated in COMSOL software. According to the material characteristic parameters and fractal parameters of the bonding surface, a mathematical model of the contact area of ​​the bonding surface is constructed in MATLAB, and a numerical solution of the contact pressure obtained by simulation is imported into the mathematical model to calculate an analytical solution of the contact area, thereby realizing analytical calculation and simulation verification of the contact area of ​​the rough contact body. The mathematical model can be further modified according to the simulation results of COMSOL, thereby providing a theoretical basis for subsequent research on contact resistance and contact thermal resistance under the influence of contact area.

Description

Technical Field

[0001] The invention belongs to the research field of microscopic contact mechanism of bonding surfaces, and specifically is a method for calculating the contact area of ​​three-dimensional rough surfaces based on COMSOL and MATLAB. Background Art

[0002] Contact resistance exists on the contact surface of electrical connectors in all electrical circuits, including power lines (strong electricity) and electronic circuits (weak electricity). It is an important indicator for evaluating the working stability and safety and reliability of electrical connectors, and the size of the contact area is an important factor affecting the contact resistance.

[0003] In practical contact models, the contacting surfaces are irregular with a certain degree of roughness. Due to this, actual contact between the two surfaces occurs at discrete contact pairs, each of which can be simplified as a pair of micro-asperities. The sum of the contact areas of all micro-asperities is the actual contact area. Therefore, studying the relationship between the actual contact area and contact load during contact on rough surfaces provides important theoretical guidance for further studying contact resistance and improving the load-carrying capacity and service life of electrical contact components.

[0004] Regarding the study of the actual contact area of ​​mating surfaces, some researchers have proposed a series of statistical microscopic contact models based on Hertz contact theory and from a statistical perspective. This type of model simplifies the actual rough surface into a surface composed of highly random or exponentially distributed spherical or other shaped micro-convex bodies. While speeding up the calculation, it also introduces errors into the model. The main reasons for the errors are: the rough surface relies heavily on the instrument's detection accuracy of the rough surface wave number range, and cannot accurately depict the rough surface morphology. On the other hand, with the advancement of computer technology and the development of numerical calculation methods, numerical calculation methods have gradually been applied to the field of contact mechanics and become a commonly used method for contact area research. Some scholars have studied the contact problem of rough surfaces with different parameters such as Poisson's ratio and surface roughness, and have also concluded that the actual contact area and contact pressure are approximately proportional.

[0005] However, neither of these two methods systematically studies the contact behavior of elastic-plastic rough surfaces. The evolution of the actual contact area of ​​these surfaces and the factors influencing it remain unclear, which in turn affects the accuracy of contact structure analysis and subsequent calculations of contact resistance and contact thermal resistance based on these structures. Furthermore, most models rely on two-dimensional contact curves when establishing contact models, which do not conform to actual three-dimensional contact characteristics. Summary of the Invention

[0006] In view of the shortcomings of the prior art, the technical problem to be solved by the present invention is to provide a method for calculating the contact area of ​​a three-dimensional rough surface.

[0007] The technical solution of the present invention to solve the above technical problem is to provide a method for calculating the contact area of ​​a three-dimensional rough surface, characterized in that the method comprises the following steps:

[0008] Step 1: Set the micromorphological parameters of the equivalent rough surface, generate point cloud data of the equivalent rough surface based on the micromorphological parameters in MATLAB software, import the point cloud data into the finite element software COMSOL, and use the parametric surface function in COMSOL to construct the equivalent rough surface;

[0009] Step 2: In COMSOL, the equivalent rough surface constructed in step 1 is used as the contact surface of a virtual contact to obtain a rough contact body; the contact surface of another virtual contact is set as a rigid smooth plane to obtain a rigid smooth contact body; the two virtual contacts together constitute an equivalent three-dimensional contact model;

[0010] Step 3. In COMSOL, determine the material characteristic parameters of the two actual contact parts, and then calculate the material characteristic parameters of the rough contact body; then import the calculated material characteristic parameters of the rough contact body into COMSOL, and then import the material hardening curve and geometric displacement curve; finally, import the material characteristic parameters of the rigid smooth contact body into COMSOL;

[0011] Step 4. Set the physical field and initial conditions in COMSOL:

[0012] Set the rigid smooth plane as a fixed constraint and the equivalent rough surface as a specified displacement. Input the displacement curve of the geometry imported in step 3 at the specified displacement.

[0013] The integral operator is added to the equivalent rough surface;

[0014] The equivalent rough surface and the rigid smooth plane are set as a contact pair;

[0015] Step 5. In COMSOL after the settings in steps 3 and 4, mesh the equivalent three-dimensional contact model obtained in step 2, and then perform simulation calculations; if the simulation results converge, the COMSOL simulation results are obtained; the COMSOL simulation results are composed of the numerical solutions of the contact pressure and the contact area; if the simulation results do not converge, return to step 3, adjust the slope of the geometric body displacement curve in step 3 and the mesh unit size of the equivalent rough surface and rough contact body in the meshing of step 5, until the simulation results converge;

[0016] Step 6. Construct a mathematical model of the contact area of ​​the bonding surface in MATLAB based on the material characteristic parameters of the two actual contact parts determined in step 3 and the micromorphology parameters of the equivalent rough surface set in step 1. Then import the numerical solution of the contact pressure obtained in step 5 into the mathematical model to calculate the MATLAB simulation result; the MATLAB simulation result is the analytical solution of the contact area.

[0017] Compared with the prior art, the present invention has the following beneficial effects:

[0018] (1) The present invention combines the analytical model of MATLAB with the simulation model of finite element software COMSOL. First, based on the fractal theory, a three-dimensional anisotropic rough surface is simulated and generated in MATLAB software, and the contact behavior of elastic-plastic rough surfaces with different material characteristic parameters and surface morphology characteristics is numerically simulated and calculated in COMSOL software. According to the material characteristic parameters and fractal parameters of the bonding surface, a mathematical model of the contact area of ​​the bonding surface is constructed in MATLAB, and the numerical solution of the contact pressure obtained by simulation is imported into the mathematical model to calculate the analytical solution of the contact area, thereby realizing the analytical calculation and simulation verification of the contact area of ​​the rough contact body. According to the simulation results of COMSOL, the mathematical model can be further modified, which provides a theoretical basis for the subsequent research on contact resistance and contact thermal resistance under the influence of contact area.

[0019] (2) The mathematical model of the present invention breaks away from the traditional method of using two-dimensional fractal theory to describe the three-dimensional contact process. It uses three-dimensional fractal functions to characterize the surface morphology of the bonding surface, which is closer to the actual situation and makes the prediction results more accurate. The mathematical model derived from the three-dimensional fractal theory is compared with the simulation results of COMSOL to verify the accuracy of the mathematical model of the present invention.

[0020] (3) The present invention uses a three-dimensional WM function to construct an equivalent rough surface, which breaks away from the previous heavy reliance on the instrument's detection accuracy of the wavenumber range of the rough surface and can uniquely characterize the rough surface morphology. BRIEF DESCRIPTION OF THE DRAWINGS

[0021] Figure 1 is a flow chart of the method of the present invention;

[0022] Figure 2 This is a surface topography image of the equivalent rough surface of Example 1 of the present invention;

[0023] Figure 3 This is an equivalent three-dimensional contact model diagram constructed in Example 1 of the present invention;

[0024] Figure 4 : is the distribution diagram of the contact pressure of the equivalent rough surface of Example 1 of the present invention;

[0025] Figure 5is the domain expansion factor of embodiment 1 of the present invention Relationship diagram with three-dimensional fractal dimension D;

[0026] Figure 6 The dimensionless actual contact area A of the COMSOL simulation results and the MATLAB simulation results of Example 1 of the present invention is r * The curve diagram of the change of normal total load F. DETAILED DESCRIPTION

[0027] The specific embodiments of the present invention are given below. The specific embodiments are only used to further illustrate the present invention and do not limit the scope of protection of the claims of this application.

[0028] The present invention provides a method for calculating the contact area of ​​a three-dimensional rough surface (hereinafter referred to as the method), characterized in that the method comprises the following steps:

[0029] Step 1: Set the micromorphological parameters of the equivalent rough surface, generate point cloud data of the equivalent rough surface based on the micromorphological parameters in MATLAB software, import the point cloud data into the finite element software COMSOL, and use the parametric surface function in COMSOL to construct the equivalent rough surface;

[0030] Preferably, in step 1, the microscopic morphology parameters include three-dimensional fractal dimension D, fractal roughness G, bonding surface side length L, minimum metric size δmin, fractal coefficient γ and ridge number M.

[0031] Preferably, in step 1, the three-dimensional WM function is used in MATLAB software to generate point cloud data of the equivalent rough surface according to the microscopic morphology parameters, and the formula is:

[0032]

[0033] In formula (1), z(x,y) represents the height of the point with coordinates x and y in the horizontal direction on the equivalent rough surface; L is the side length of the bonding surface of the equivalent rough surface; G is the fractal roughness of the equivalent rough surface; D is the three-dimensional fractal dimension of the equivalent rough surface; γ is the fractal coefficient, which is usually taken as 1.5; M is the number of ridges of the equivalent rough surface; δmin is the minimum metric size of the equivalent rough surface; φ m,n is a random phase; m is [1, M]; n is a frequency index.

[0034] Step 2. In COMSOL, the contact mechanics of the equivalent rough surface is equivalent to the contact between a rough surface and a rigid smooth plane: the equivalent rough surface constructed in step 1 is used as the contact surface of a virtual contact part to obtain a rough contact body; the contact surface of another virtual contact part is set as a rigid smooth plane to obtain a rigid smooth contact body; the two virtual contact parts (i.e., the rough contact body and the rigid smooth contact body) together constitute an equivalent three-dimensional contact model;

[0035] Step 3. In COMSOL, determine the respective materials and material characteristic parameters of the two actual contact parts, and then calculate the material characteristic parameters of the rough contact body; then import the calculated material characteristic parameters of the rough contact body into COMSOL as the material characteristic parameters of the rough contact body, and then import the material hardening curve and geometric body displacement curve; finally, import the material characteristic parameters of the rigid smooth contact body into COMSOL;

[0036] Preferably, in step 3, the material characteristic parameters include Young's modulus, Poisson's ratio and yield strength. The material characteristic parameters of the two actual contact parts are obtained from the COMSOL software library. The material characteristic parameters of the rough contact body include equivalent Young's modulus E, equivalent Poisson's ratio υ and equivalent yield strength σ; the calculation formula is υ=(υ1+υ2) / 2; where E1 and E2 are the Young's moduli of the two actual contact materials, and υ1 and υ2 are the Poisson's ratios of the two actual contact materials. The equivalent yield strength σ is the yield strength of the softer material of the two actual contact parts. A rigid smooth contact body is an ideal rigid body with infinite Young's modulus, zero Poisson's ratio, and infinite yield strength.

[0037] Preferably, in step 3, the material hardening curve is obtained through the material library inside COMSOL; the displacement curve of the geometric body is a piecewise linear function with respect to time, wherein the segmentation points are smoothed using continuous second-order inverses, and the slope of the piecewise linear function is large at first and small at the end, so that the calculation results of the COMSOL transient model converge.

[0038] Step 4. Set the physical field and initial conditions in COMSOL:

[0039] Set the rigid smooth plane as a fixed constraint and the equivalent rough surface as a specified displacement. Input the displacement curve of the geometry imported in step 3 at the specified displacement.

[0040] An integral operator is added to the equivalent rough surface, which is used for meshing in step 5 and then solving the calculation to obtain the numerical solution of the contact pressure of the entire equivalent rough surface;

[0041] The equivalent rough surface and the rigid smooth plane are set as a contact pair, and the contact method is penalty function;

[0042] Step 5. In COMSOL after the settings in steps 3 and 4, mesh the equivalent three-dimensional contact model obtained in step 2, and then perform simulation calculations in the transient solver; if the simulation results converge, the COMSOL simulation results are obtained; the COMSOL simulation results are composed of the numerical solutions of the contact pressure and the contact area; if the simulation results do not converge (i.e., the results cannot be calculated), return to step 3, adjust the slope of the geometric body displacement curve in step 3 and the mesh unit size of the equivalent rough surface and rough contact body in the meshing of step 5, until the simulation results converge;

[0043] Preferably, in step 5, meshing is performed by meshing the equivalent rough surface with a free triangle mesh, and then meshing the entire rough contact body with a free tetrahedron mesh; and meshing the rigid smooth contact body with a COMSOL sweep function to reduce the simulation calculation time.

[0044] Preferably, in step 5, the contact pressure numerical solution is: simulating and calculating the pressure distribution data of the equivalent rough surface under a specified displacement (i.e., the pressure value of each contact point on the equivalent rough surface), and then summing the pressure values ​​of each contact point to obtain the contact pressure numerical solution;

[0045] Preferably, in step 5, the numerical solution of the contact area is: first, the pressure of the equivalent rough surface contact area is set to 1, and the pressure of the non-contact area is set to 0, and then the contact area is integrated using an integral operator to obtain the numerical solution of the contact area.

[0046] Step 6: Construct a mathematical model of the contact area of ​​the bonding surface in MATLAB based on the material characteristic parameters of the two actual contact parts determined in Step 3 and the microscopic morphology parameters of the equivalent rough surface set in Step 1. Then, import the numerical solution of the contact pressure obtained in Step 5 into the mathematical model to calculate the MATLAB simulation result; the MATLAB simulation result is the analytical solution of the contact area;

[0047] The bonding surface is a pair of surfaces in contact state, that is, the equivalent rough surface and the rigid smooth plane in contact state constitute the bonding surface.

[0048] Preferably, in step 6, the specific steps of establishing the mathematical model are as follows:

[0049] (6.1) Assuming the number of ridges M = 1, the equivalent rough surface is cross-sectioned, and the cross-sectional profile curve Z(x) of the equivalent rough surface in two dimensions is obtained as follows:

[0050]

[0051] In formula (2), G is the fractal roughness; D is the three-dimensional fractal dimension; γ is the fractal coefficient; l is the size width of the micro-convex body;

[0052] Then, based on the cross-sectional profile curve, the deformation height δ of a single asperity on the equivalent rough surface is obtained as:

[0053]

[0054] The peak curvature radius R of the microconvex body is:

[0055]

[0056] In formula (4), δ min The smallest measurement size; is the peak curvature radius of the micro-convex body under the minimum measurement size;

[0057] (6.2) Calculate the critical asperity contact area a of a single asperity on the equivalent rough surface c ;

[0058] According to Hertz contact theory, under the action of plane pressure, the critical deformation δ of the micro-convex body undergoing elastic deformation is c for:

[0059]

[0060] In formula (5), σ is the equivalent yield strength of the rough contact body; K is the hardness coefficient, which is related to the equivalent Poisson's ratio υ of the rough contact body: K = 0.4645 + 0.3141υ + 0.1943υ 2 ; E is the equivalent Young's modulus of the rough contact body;

[0061] The relationship between the width l of the micro-convex body and the contact area a of the micro-convex body is l=a 0.5 , and the critical asperity contact area a is obtained from this c :

[0062]

[0063] (6.3) According to Hertz contact theory, when the peak curvature radius R of the micro-convex body before the contact point deformation is much larger than the deformation height δ, the elastic contact load F e The relationship between (a) and the contact area a is as follows:

[0064]

[0065] When the contact surface is locally plastically deformed, the plastic contact load F p (a) is:

[0066] F p (a)=Kσa (8)

[0067] (6.4) The distribution of the asperity contact area a follows the distribution of the Earth's island area, so the size distribution function n(a) of the asperity contact area is expressed as:

[0068]

[0069] In formula (9), a l is the maximum asperity contact area; is the domain expansion factor of the asperity contact area size distribution (referred to as the domain expansion factor), and is a function of the three-dimensional fractal dimension D, satisfying:

[0070]

[0071] Use MATLAB to solve equation (10) and get and the solution of D (i.e. and D in one-to-one correspondence); then Fit the solution of D to The implicit function of D (i.e., Equation 10) is fitted into an explicit function, and the fitting result is obtained.

[0072] Total contact area A of equivalent rough surface r for:

[0073]

[0074] In formula (11), a s is the minimum contact area of ​​the micro-convex body, which is 0 under the infinite subdivision scale of the fractal;

[0075] To simplify the expression, Substituting into formula (11), we get formula (12):

[0076]

[0077] In formula (12), A r1 is the domain expansion factor after fitting by formula (10) The contact area under the rough surface (i.e., the actual contact area of ​​the equivalent rough surface);

[0078] (6.5) According to equations (6), (7), (8), (9) and (12), the total normal load F and A are calculated. r1 The relationship is:

[0079]

[0080] In formula (12), F is the numerical solution of the contact pressure calculated in step 5, so according to formula (12) we can get A r1; For the rough contact model, the dimensionless actual contact area A is often used r * express, Among them A a =L 2 is the nominal contact area, so according to A r1 Calculated dimensionless actual contact area A r * (i.e., the analytical solution of the contact area) is the simulation result of MATLAB.

[0081] Preferably, the method also includes step 7: importing the numerical solution of contact area obtained in step 5 into MATLAB, and using MATLAB software to draw curves of the numerical solution of contact pressure and the analytical solution of contact area, and curves of the numerical solution of contact pressure and the numerical solution of contact area in the same coordinate system; and then using the relative error formula to calculate the relative error Δ of the two curves to verify the correctness of the simulation results and the accuracy of the mathematical model, so as to achieve the effect of mutual verification between the simulation results of COMSOL and the simulation results of MATLAB.

[0082] Preferably, in step 7, the relative error formula is:

[0083]

[0084] In formula (13), N is the number of discrete measurement points of the two curves, and Δ is the relative error.

[0085] Example 1

[0086] Step 1: Set the three-dimensional fractal dimension D = 2.4 and the fractal roughness G = 2.56×10 -9 m, the side length of the bonding surface L = 1.5 × 10 -4 m, minimum measurement size δ min =10 -8 m, fractal coefficient γ = 1.5, ridge number M = 1, generate the three-dimensional point cloud data of the equivalent rough surface in MATLAB software according to formula (1), then import the three-dimensional point cloud data into the finite element software COMSOL, and use the parametric surface function in COMSOL to construct the equivalent rough surface, as shown in Figure 2 As shown, in the COMSOL settings, the relative tolerance is set to 10 -6 , the maximum number of knots is set to 2000.

[0087] Step 2: Construct an equivalent three-dimensional contact model in COSMOL: Use two cuboids as the main structures of the rough contact body and the rigid smooth contact body, use the equivalent rough surface as the contact surface of the rough contact body, and set the contact surface of the other cuboid as a rigid smooth plane; the two cuboids with contact surface features together constitute an equivalent three-dimensional contact model, such as Figure 3As shown;

[0088] Step 3: In COMSOL, determine the materials of the two actual contact parts. In this embodiment, aluminum and copper are used as examples.

[0089] According to the COMSOL material library data, the Young's modulus of aluminum is E1 = 70 GPa, Poisson's ratio υ1 = 0.33, and yield strength σ1 = 220 MPa; the Young's modulus of copper is E2 = 126 GPa, Poisson's ratio υ2 = 0.34, and yield strength σ2 = 300 MPa. The equivalent Young's modulus of the rough contact body is E = 50.64 GPa, the equivalent Poisson's ratio υ = 0.34, and the equivalent yield strength σ is the yield strength of aluminum = 220 MPa;

[0090] The calculated equivalent Young's modulus E, equivalent Poisson's ratio υ, and equivalent yield strength σ are imported into COMSOL as the material characteristic parameters of the rough contact body, and then the material hardening curve and geometric body displacement curve are imported; finally, the material characteristic parameters of the rigid smooth contact body are imported into COMSOL, and its Young's modulus is infinite, its Poisson's ratio is 0, and its yield strength is infinite;

[0091] Step 4: Set the rigid smooth plane as a fixed constraint, the equivalent rough surface as a specified displacement, and input the displacement curve of the geometric body at the specified displacement; add the integral operator intop3 to the equivalent rough surface; set the equivalent rough surface and the rigid smooth plane as a contact pair, and use the penalty function as the contact method;

[0092] Step 5: In COMSOL after the settings in steps 3 and 4, mesh the equivalent three-dimensional contact model obtained in step 2, and then perform simulation calculations to obtain COMSOL simulation results (i.e., numerical solutions for contact pressure and contact area);

[0093] The meshing is as follows: the equivalent rough surface is meshed using a free triangle mesh with the element size set to normal; the entire rough contact body is meshed using a free tetrahedron mesh with the element size set to coarse; the rigid smooth contact body is meshed using the COMSOL sweep function to reduce the simulation calculation time;

[0094] Among them, the distribution diagram of contact pressure ( Figure 4 ) to obtain the distribution of contact pressure, and then obtain the contact area corresponding to each contact pressure, and then use the integral operator intop3 to integrate the area of ​​the contact area to obtain the total contact area (i.e., the numerical solution of the contact area). Figure 4 It can be seen that with the increase of contact depth, the actual contact area of ​​the equivalent rough surface also increases. At the same time, the contact stress cloud map tends to be strip-shaped as a whole, reflecting the anisotropic characteristics of the three-dimensional VM function.

[0095] The integral operator is intop3, and its discriminant is: intop3(if(solid.Tn>0,1,0)); where solid.Tn is the contact pressure distribution of the rigid smooth plane, and intop3 is the contact surface integral operator;

[0096] Step 6: Construct a mathematical model of the contact area of ​​the bonding surface in MATLAB based on the material characteristic parameters of the two actual contact parts determined in Step 3 and the microscopic morphology parameters of the equivalent rough surface set in Step 1. Then, import the numerical solution of the contact pressure obtained in Step 5 into the mathematical model to calculate the MATLAB simulation result; the MATLAB simulation result is the analytical solution of the contact area;

[0097] Among them, according to the COMSOL material library data, the following material characteristic parameters can be set in the mathematical model: the hardness coefficient of aluminum K = 0.5893, the hardness coefficient of copper K = 0.5934; after substituting the material characteristic parameters, l>δ can be obtained according to formula (6). min and l≤δ min Critical asperity contact area a under c are 1.1593×10 -9 m 2 and 1.8843×10 -20 m 2 ; Use MATLAB to solve equation (10) and use MATLAB toolbox cftool to fit the rational function. The domain expansion factor The curve relationship with the three-dimensional fractal dimension D is as follows Figure 5 As shown, when D = 2.4, we can get Finally, according to equations (6) and (12), the total contact area Ar of the equivalent rough surface can be calculated; the side length of the bonding surface L = 0.0014m, the nominal contact area Aa = 1.96×10 -6 m 2 , and then calculate the dimensionless actual contact area Ar * Relationship with the numerical solution F of contact pressure;

[0098] Step 7. Import the contact area numerical solution obtained in step 5 into MATLAB. Use MATLAB software to draw the curves of the contact pressure numerical solution and the contact area analytical solution, and the curves of the contact pressure numerical solution and the contact area numerical solution in the same coordinate system. The simulation results of COMSOL and MATLAB are as follows: Figure 6 As shown. Figure 6 It can be seen that the dimensionless actual contact area Ar * As the contact pressure increases, the numerical solution increases approximately linearly. The two results are consistent with each other, and the relative error Δ is within 5%, which is 3.4%, which verifies the accuracy of the MATLAB simulation results. Figure 6 In the simulation, COMSOL's simulation results will become more linear as the number of acquisition points increases, but this will also increase the amount of calculation and make convergence difficult.

[0099] Any matters not described in the present invention are applicable to the prior art.

Claims

1. A method for calculating the contact area of ​​a three-dimensional rough surface, characterized in that: The method comprises the following steps: Step 1: Set the micromorphological parameters of the equivalent rough surface, generate point cloud data of the equivalent rough surface based on the micromorphological parameters in MATLAB software, import the point cloud data into the finite element software COMSOL, and use the parametric surface function in COMSOL to construct the equivalent rough surface; Step 2: In COMSOL, the equivalent rough surface constructed in step 1 is used as the contact surface of a virtual contact to obtain a rough contact body; the contact surface of another virtual contact is set as a rigid smooth plane to obtain a rigid smooth contact body; the two virtual contacts together constitute an equivalent three-dimensional contact model; Step 3. In COMSOL, determine the material characteristic parameters of the two actual contact parts, and then calculate the material characteristic parameters of the rough contact body; then import the calculated material characteristic parameters of the rough contact body into COMSOL, and then import the material hardening curve and geometric displacement curve; finally, import the material characteristic parameters of the rigid smooth contact body into COMSOL; Step 4. Set the physical field and initial conditions in COMSOL: Set the rigid smooth plane as a fixed constraint and the equivalent rough surface as a specified displacement. Input the displacement curve of the geometry imported in step 3 at the specified displacement. The integral operator is added to the equivalent rough surface; The equivalent rough surface and the rigid smooth plane are set as a contact pair; Step 5. In COMSOL after the settings in steps 3 and 4, mesh the equivalent three-dimensional contact model obtained in step 2, and then perform simulation calculations; if the simulation results converge, the COMSOL simulation results are obtained; the COMSOL simulation results are composed of the numerical solutions of the contact pressure and the contact area; if the simulation results do not converge, return to step 3, adjust the slope of the geometric body displacement curve in step 3 and the mesh unit size of the equivalent rough surface and rough contact body in the meshing of step 5, until the simulation results converge; Step 6. Construct a mathematical model of the contact area of ​​the bonding surface in MATLAB based on the material characteristic parameters of the two actual contact parts determined in step 3 and the micromorphology parameters of the equivalent rough surface set in step 1. Then import the numerical solution of the contact pressure obtained in step 5 into the mathematical model to calculate the MATLAB simulation result; the MATLAB simulation result is the analytical solution of the contact area.

2. The method for calculating the contact area of ​​a three-dimensional rough surface according to claim 1, wherein: In step 1, the microscopic morphology parameters include three-dimensional fractal dimension D, fractal roughness G, bonding surface side length L, minimum metric size δmin, fractal coefficient γ and ridge number M; The three-dimensional WM function is used in MATLAB software to generate point cloud data of the equivalent rough surface according to the microscopic morphology parameters. The formula is: In formula (1), z(x,y) represents the height of the point with horizontal coordinates x and y on the equivalent rough surface; φ m,n is a random phase; m is [1, M]; n is a frequency index.

3. The method for calculating the contact area of ​​a three-dimensional rough surface according to claim 1, wherein: In step 3, the material characteristic parameters include Young's modulus, Poisson's ratio and yield strength; The material characteristic parameters of the two actual contact parts are obtained from the COMSOL software library; The material characteristic parameters of the rough contact body include equivalent Young's modulus E, equivalent Poisson's ratio υ and equivalent yield strength σ; the calculation formula is: υ=(υ1+υ2) / 2; where E1 and E2 are the Young's modulus of the two actual contact materials, υ1 and υ2 are the Poisson's ratios of the two actual contact materials, and the equivalent yield strength σ is the yield strength of the softer material of the two actual contact parts. The rigid smooth contact body is an ideal rigid body, whose Young's modulus is infinite, its Poisson's ratio is 0, and its yield strength is infinite.

4. The method for calculating the contact area of ​​a three-dimensional rough surface according to claim 1, wherein: In step 3, the material hardening curve is obtained from the material library inside COMSOL; the displacement curve of the geometric body is a piecewise linear function of time, in which the segment points are smoothed using continuous second-order inverses.

5. The method for calculating the contact area of ​​a three-dimensional rough surface according to claim 1, wherein: In step 5, the meshing is as follows: the equivalent rough surface is meshed with a free triangle mesh, and then the entire rough contact body is meshed with a free tetrahedron; the rigid smooth contact body is meshed using the COMSOL sweep function.

6. The method for calculating the contact area of ​​a three-dimensional rough surface according to claim 1, wherein: In step 5, the contact pressure numerical solution is: simulate and calculate the pressure value of each contact point on the equivalent rough surface under the specified displacement, and then sum the pressure values ​​of each contact point to obtain the contact pressure numerical solution.

7. The method for calculating the contact area of ​​a three-dimensional rough surface according to claim 1, wherein: In step 5, the numerical solution of the contact area is: first set the pressure of the equivalent rough surface contact area to 1 and the pressure of the non-contact area to 0 to obtain each contact area, and then use the integral operator to integrate the area of ​​the contact area to obtain the numerical solution of the contact area.

8. The method for calculating the contact area of ​​a three-dimensional rough surface according to claim 1, wherein: The specific steps of step 6 are as follows: (6.1) Assuming the number of ridges M = 1, the equivalent rough surface is cross-sectioned, and the cross-sectional profile curve Z(x) of the equivalent rough surface in two dimensions is obtained as follows: In formula (2), G is the fractal roughness; D is the three-dimensional fractal dimension; γ is the fractal coefficient; l is the size width of the micro-convex body; Then, based on the cross-sectional profile curve, the deformation height δ of a single asperity on the equivalent rough surface is obtained as: The peak curvature radius R of the microconvex body is: In formula (4), δ min The smallest measurement size; is the peak curvature radius of the micro-convex body under the minimum measurement size; (6.2) Calculate the critical asperity contact area a of a single asperity on the equivalent rough surface c ; According to Hertz contact theory, under the action of plane pressure, the critical deformation δ of the micro-convex body undergoing elastic deformation is c for: In formula (5), σ is the equivalent yield strength of the rough contact body; K is the hardness coefficient, which is related to the equivalent Poisson's ratio υ of the rough contact body: K = 0.4645 + 0.3141υ + 0.1943υ 2 ; E is the equivalent Young's modulus of the rough contact body; The relationship between the width l of the micro-convex body and the contact area a of the micro-convex body is l=a 0.5 , and the critical asperity contact area a is obtained from this c : (6.3) According to Hertz contact theory, when the peak curvature radius R of the micro-convex body before the contact point deformation is much larger than the deformation height δ, the elastic contact load F e The relationship between (a) and the contact area a is as follows: When the contact surface is locally plastically deformed, the plastic contact load F p (a) is: F p (a)=Kσa (8) (6.4) The distribution of the asperity contact area a follows the distribution of the Earth's island area, so the size distribution function n(a) of the asperity contact area is expressed as: In formula (9), a l is the maximum asperity contact area; is the domain expansion factor for the asperity contact area size distribution and is a function of the three-dimensional fractal dimension D, satisfying: Use MATLAB to solve equation (10) and get and the solution of D; then Fit with the solution of D to get the fitting result Total contact area A of equivalent rough surface r for: In formula (11), a s is the minimum contact area of ​​the micro-convex body, which is 0 under the infinite subdivision scale of the fractal; To simplify the expression, Substituting into formula (11), we get formula (12): In formula (12), A r1 Equivalent rough surface true contact area; (6.5) According to equations (6), (7), (8), (9) and (12), the total normal load F and A are calculated. r1 The relationship is: In formula (12), F is the numerical solution of the contact pressure calculated in step 5, so according to formula (12) we can get A r1 ; Dimensionless actual contact area Among them A a =L 2 is the nominal contact area, so the MATLAB simulation result is obtained; the MATLAB simulation result is the contact area analytical solution A r * .

9. The method for calculating the contact area of ​​a three-dimensional rough surface according to claim 1, wherein: The method also includes step 7: importing the numerical solution of contact area obtained in step 5 into MATLAB, using MATLAB software to draw curves of the numerical solution of contact pressure and the analytical solution of contact area, and curves of the numerical solution of contact pressure and the numerical solution of contact area in the same coordinate system; and then using the relative error formula to calculate the relative error Δ between the two curves to verify the correctness of the simulation results and the accuracy of the mathematical model.

10. The method for calculating the contact area of ​​a three-dimensional rough surface according to claim 9, wherein: In step 7, the relative error formula is: In formula (13), N is the number of discrete measurement points of the two curves, and Δ is the relative error.

Images