A method for calculating a smooth contact potential of an arbitrary shape of a continuous-discontinuous model

Abstract

The application discloses a method for calculating a smooth contact potential of an arbitrary shape of a continuous-discontinuous model, comprising the following steps: S1, constructing a continuous-discontinuous coupling simulation model, initializing model parameters, discriminating contact pairs, identifying and storing element pairs with possible contact relations, and establishing a data list of contact element pairs; S2, for each contact element pair, assigning a master element and a slave element, and judging whether the contact element pair is overlapped; S3, for the contact element pair with overlap, recording the vertex coordinates of the master element embedded in the slave element, the vertex is called an embedded vertex, and a weight function model of the embedded vertex is established; S4, based on the weight function model established in step S3, calculating the smooth potential function of the embedded vertex in the slave element; and S5, according to the smooth potential function calculated in step S4, calculating the normal contact force generated by the master element of the slave element pair, and further distributing the calculated contact force into equivalent node forces of the master element and the slave element according to a shape function.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of geotechnical engineering numerical simulation, and particularly relates to a calculation method of a smooth contact potential of an arbitrary shape of a continuous-discontinuous model. BACKGROUND

[0002] In current engineering practice, numerical calculation methods such as finite element method (FEM), boundary element method (BEM), finite difference method (FDM) and discrete element method (DEM) have been widely applied in the field of geotechnical engineering stability analysis. Through the construction of geomechanical models, the stress-strain response of rock mass under different load conditions can be predicted. These methods provide important numerical analysis tools for engineering structure optimization design, construction safety evaluation and geological disaster warning.

[0003] However, the prior art has the following significant defects: the traditional continuous medium mechanics method is difficult to effectively simulate the dynamic evolution of cracks, the penetration of fracture network and the motion behavior of blocks, and cannot completely characterize the whole process of the failure mechanism of rock mass from continuous to discontinuous state; the typical discontinuous medium analysis method has inherent contradictions between the prediction accuracy of crack propagation path and the calculation efficiency, especially under large-scale complex conditions, the calculation cost is high, the convergence is poor and the numerical oscillation problem is particularly prominent; the existing algorithm has limitations in characterizing the multi-scale structural characteristics of rock mass, and cannot effectively couple the mesoscopic damage evolution and macroscopic structural instability process, resulting in insufficient reliability of the prediction of complex jointed rock mass failure mode.

[0004] This technical bottleneck is due to the multi-phase heterogeneous characteristics, anisotropic mechanical response and complex joint network distribution characteristics of rock mass material itself. Its failure process essentially shows cross-scale progressive failure under the coupling of multiple physical fields, involving nonlinear and discontinuous mechanical behaviors such as micro-crack initiation, propagation, penetration, and block sliding and separation.

[0005] As a crucial part of the discrete block simulation process, the contact potential usually occupies a large amount of calculation resources in the overall calculation of the continuous-discontinuous coupled model. There are three key problems in the traditional contact potential calculation method in practical application: first, under the same embedding depth condition, the contact force calculation results have significant differences, which leads to discontinuous mechanical response, affecting the stability and accuracy of the simulation; second, in the modeling process using quadrilateral elements and higher-order variable elements, the conventional contact potential function has contact force singularity or unreasonable oscillation in some characteristic areas, which seriously affects the continuity and rationality of the crack propagation path; third, the existing contact potential model is only for two-dimensional triangular grid and three-dimensional tetrahedron.

[0006] Therefore, it is necessary to design a new smooth contact potential calculation method for arbitrary element shape to overcome the above problems. SUMMARY

[0007] To achieve the above object, the technical scheme adopted by the present application is:

[0008] A method for calculating a smooth contact potential of an arbitrary shape of a continuous-discontinuous model, comprising the following steps:

[0009] S1, a continuous-discontinuous coupled simulation model is constructed, model parameters are initialized, contact pair discrimination is performed, and element pairs with possible contact relationships are identified and stored, referred to as contact element pairs, and a data list of the contact element pairs is established;

[0010] S2, for each contact element pair, main and slave element allocation is performed, and fine judgment is performed on whether the contact element pair overlaps;

[0011] S3, for the contact element pair that overlaps, the vertex coordinates of the main element embedded in the slave element are recorded, the vertex is referred to as an embedded vertex, and a weight function model of the embedded vertex is established;

[0012] S4, based on the weight function model established in step S3, the smooth potential function of the embedded vertex in the slave element is directly calculated;

[0013] S5, according to the smooth potential function calculated in step S4, the normal contact force generated by the main element of the slave element pair is calculated, and the calculated contact force is further distributed to the equivalent node forces of the main and slave elements according to the shape function.

[0014] Further, in step S1, the continuous-discontinuous coupled simulation model constructed is a two-dimensional model or a three-dimensional model of an arbitrary element shape.

[0015] Further, in step S2, specifically:

[0016] For each contact element pair, the main element e m and a slave element e s are selected from the element ID arrangement;

[0017] A coordinate list of each vertex of the main and slave elements is established, if the continuous-discontinuous coupled simulation model is a two-dimensional model, an outer normal vector list of each edge of the main and slave elements is established, and if the continuous-discontinuous coupled simulation model is a three-dimensional model, an outer normal vector list of each plane of the main and slave elements is established;

[0018] Based on the separation axis theorem, it is judged whether the main and slave elements overlap, if any vertex x p of the main element, any vertex v j of the slave element, and the outer normal vector f of each edge or each plane of the main and slave elements satisfy:

[0019] max(x p ·f)-min(x j• f) > 0

[0020] or

[0021] max(x j • f) - min(x p • f) > 0

[0022] The master and slave units are considered to overlap, where x p is any vertex in the list of vertex coordinates of the master unit, v j represents any vertex in the list of vertex coordinates of the slave unit, and f represents the outward normal vector of any edge or any plane of the master and slave units.

[0023] Further, in step S3, specifically:

[0024] For the contact unit pair that overlaps, the vertex coordinate x p of the master unit embedded in the slave unit is recorded, which is referred to as the embedded vertex;

[0025] Assuming that there are n vertices in the list of coordinates of the vertices of the slave unit, when establishing the weight function model of the embedded vertex x p , the weight ratio of the point to each vertex of the slave unit is calculated:

[0026]

[0027] where N i represents the weight ratio of the i-th vertex v i of the slave unit, w i (x p ) is the weight of the vertex v i , the denominator is the sum of the weights of all vertices, and i = 1, 2, …, n;

[0028] Assuming that the vertex v i is on l planes or edges, then w i (x p ) is calculated by the following formula:

[0029]

[0030] where:

[0031]

[0032]

[0033] where, is the shortest distance from the embedded vertex x p to the a-th plane or edge where the vertex v i is located, is the weight of the vertex v iThe normal vector of the a-th plane or edge in which it is located, where a = 1, 2, ..., l.

[0034] Furthermore, in step S3, if the continuous-discontinuous coupling simulation model is a two-dimensional model, then w i (x p Specifically, it is calculated using the following formula:

[0035]

[0036] Where A i It is the area of ​​the triangle formed by the i-th vertex of the unit and its adjacent vertices. The formula for calculation is as follows:

[0037]

[0038] Where v i-1 v i and v i+1 These are from the list of unit vertex coordinates V s The coordinates of the (i-1), i, and i+1th vertices in the array;

[0039] Indicates that the main unit is embedded at vertex x p The area of ​​the triangle constructed from the (i-1)th and ith vertices in the cell, This is the main unit embedded in vertex x p The area of ​​the triangle constructed from the i-th and i+1-th vertices in the unit is calculated using the following formula:

[0040]

[0041]

[0042] w is obtained using the above method. i (x p After that, the weight ratio N is obtained by further normalization calculation. i .

[0043] Furthermore, in step S3, if the continuous-discontinuous coupled simulation model is a three-dimensional model, then w i (x p Specifically, it is calculated using the following formula:

[0044]

[0045] From the i-th vertex v of the unit i Adjacent vertices are called adjacent vertices, m represents the number of adjacent vertices, and V is the vertex v. i The volume of the polyhedron formed by m adjacent vertices is decomposed into several tetrahedra for solution, with v as the unit of volume. iFor the vertex, the triangle f is the single tetrahedron volume of the bottom surface The calculation is performed by the following formula:

[0046]

[0047] Wherein the subscript f1, f2, f3 is the vertex number of the triangle f arranged counterclockwise, the superscript (1), (2), (3) represents the x, y, z coordinate, and then the sum of the several individual tetrahedron volumes is obtained V i ;

[0048] The jth adjacent vertex is represented by v j , j = 1, 2, …, m, The main unit embedded vertex x p , the vertex v i , and the polyhedral volume formed by the other m-1 adjacent vertices except the adjacent vertex v j ;

[0049] After w i (x p ) is obtained by the above method, the weight ratio N i is further normalized to obtain the weight ratio.

[0050] Further, in step S4, when calculating the normalized potential function of the embedded vertex in the unit, the calculation formula is as follows:

[0051]

[0052] Compared with the prior art, the continuous-discontinuous model smooth contact potential calculation method for arbitrary shape provided by the application solves the problems of non-smoothness, high calculation complexity, non-conservative force, and inability to adapt to various problems of arbitrary unit shape in the prior art, so that the calculation result is closer to the theoretical value, which is beneficial to improve the stability and precision of simulation, and the application range is improved, and the method is suitable for contact potential calculation of two-dimensional and three-dimensional grid of arbitrary geometric shape. BRIEF DESCRIPTION OF DRAWINGS

[0053] Figure 1 The flowchart of the continuous-discontinuous model smooth contact potential calculation method provided by the application is shown in the figure;

[0054] Figure 2 The schematic diagram of the smooth potential distribution of the triangular unit is shown in the figure;

[0055] Figure 3 The schematic diagram of the smooth potential distribution of the quadrilateral unit is shown in the figure;

[0056] Figure 4 The schematic diagram of the smooth potential distribution of the multi-modal unit is shown in the figure;

[0057] Figure 5 Modeling diagrams of the methods provided by the conventional fillet method, the Munjiza method and the method provided by the present application;

[0058] Figure 6 A diagram of the friction force calculated by the conventional fillet method, the Munjiza method and the method provided by the present application changing with the sliding distance. DETAILED DESCRIPTION

[0059] In order to make the technical means, creative features, purposes and effects of the present application easy to understand, the following further describes how the present application is implemented in combination with specific embodiments.

[0060] Referring to Figure 1 , the present application provides a smooth contact potential calculation method for a continuous-discontinuous model with an arbitrary shape, which comprises the following steps:

[0061] S1, a continuous-discontinuous coupled simulation model is constructed, model parameters are initialized, contact pair discrimination is performed, and element pairs with possible contact relationships are identified and stored, which are referred to as contact element pairs, and a data list of the contact element pairs is established;

[0062] S2, for each contact element pair, master and slave element allocation is performed, and whether the contact element pair overlaps is finely judged;

[0063] S3, for the contact element pair that overlaps, the vertex coordinates of the master element embedded in the slave element are recorded, which is referred to as an embedded vertex, and a weight function model of the embedded vertex is established;

[0064] S4, based on the weight function model established in step S3, the smooth potential function of the embedded vertex in the slave element is directly calculated;

[0065] S5, according to the smooth potential function calculated in step S4, the normal contact force generated by the master element of the slave element pair is calculated, and the calculated contact force is further distributed to the equivalent node forces of the master and slave elements according to the shape function.

[0066] In the first specific embodiment:

[0067] In step S1, a continuous-discontinuous coupled simulation model is constructed, model parameters are initialized, contact pair discrimination is performed, and element pairs with possible contact relationships are identified and stored, which are referred to as contact element pairs, and a data list of the contact element pairs is established. The constructed continuous-discontinuous coupled simulation model is a two-dimensional model or a three-dimensional model.

[0068] In step S2, for each contact element pair, the master element e m and a slave element e s are selected from the element ID arrangement; Figure 2 andFigure 3 The slave unit can be a triangle, a quadrilateral or other shape, and the master unit is a two-dimensional planar figure;

[0069] A coordinate list of each vertex of the master and slave units is established, and a list of outer normal vectors of each edge of the master and slave units is established;

[0070] Based on the Separating Axis Theorem, it is determined whether the master and slave units overlap. If any vertex x p of the master unit satisfies the following conditions with respect to any vertex v j of the slave unit and the outer normal vector f of each edge of the master and slave units:

[0071] max(x p ·f)-min(v j ·f)≥0

[0072] or

[0073] max(v j ·f)-min(x p ·f)≥0

[0074] then it is considered that the master and slave units overlap, wherein x p is any vertex in the coordinate list of the vertices of the master unit, v j represents any vertex in the coordinate list of the vertices of the slave unit, and f represents the outer normal vector of any edge of the master and slave units.

[0075] In step S3, for the contact unit pair that overlaps, the vertex coordinate x p of the master unit embedded in the slave unit is recorded, and this vertex is referred to as an embedded vertex;

[0076] Assuming that there are n vertices in the coordinate list of the vertices of the slave unit, when establishing a weight function model of the embedded vertex x p , the weight ratio of this point with respect to each vertex of the slave unit is calculated:

[0077]

[0078] wherein N i represents the weight ratio of the i-th vertex v i of the slave unit, w i (x p ) is the weight of the vertex v i , the denominator is the sum of the weights of all vertices, and i=1, 2, …, n;

[0079] Assuming that the vertex v i is on l edges, then w i (x p ) is calculated by the following formula:

[0080]

[0081] in:

[0082]

[0083]

[0084] in, It is an embedded vertex x p To vertex v i The shortest distance to the a-th edge. For vertex v i The normal vector of the a-th edge, where a = 1, 2, ..., l.

[0085] Since the continuous-discontinuous coupling simulation model is a two-dimensional model, then w i (x p Specifically, it can be calculated using the following formula:

[0086]

[0087] Where A i It is the area of ​​the triangle formed by the i-th vertex of the unit and its adjacent vertices. The formula for calculation is as follows:

[0088]

[0089] Where v i-1 v i and v i+1 These are from the list of unit vertex coordinates V s The coordinates of the (i-1), i, and i+1th vertices in the array;

[0090] Indicates that the main unit is embedded at vertex x p The area of ​​the triangle constructed from the (i-1)th and ith vertices in the cell, This is the main unit embedded in vertex x p The area of ​​the triangle constructed from the i-th and i+1-th vertices in the unit is calculated using the following formula:

[0091]

[0092]

[0093] w is obtained using the above method. i (x p After that, the weight ratio N is obtained by further normalization calculation. i .

[0094] In step S4, it is understandable that, due to the weight ratio N i satisfy:

[0095]

[0096] And for any vertex of the cell, N i >0, therefore the embedding vertex x of the principal unit p The interpolation function within the element can be expressed as:

[0097]

[0098] Therefore N i That is, the main unit embeds vertex x p For the i-th vertex v in the cell i The shape function, when x p =v i There are N times i =1, therefore the shape function N i satisfy:

[0099]

[0100] Its value is maximized at the center, that is As can be seen from the above equation, the shape function gradually decreases as the embedding point approaches the boundary, and then decreases to zero at the boundary.

[0101] When calculating the potential function of the embedded vertex in the cell, a smooth potential function is constructed by multiplying the areas of multiple adjacent triangles of the polygon. The specific calculation formula is as follows:

[0102]

[0103] In step S5, based on the potential function obtained in step S4, the normal contact force generated by the slave element on the master element is calculated, and the calculated contact force is further distributed to the equivalent nodal forces of the master and slave elements according to the weighting function model (shape function). This completes the calculation of the contact force.

[0104] Steps S2-S5 can then be repeated until the initialization, boundary loading, stress equilibrium, and simulation completion of the continuous-discontinuous coupled simulation model are achieved. Furthermore, stress changes in the model should be monitored throughout the process to ensure that no significant oscillations occur during the model calculation.

[0105] Reference Figure 2 and Figure 3 The figures show the calculation principle and smooth potential distribution diagrams for triangular and quadrilateral elements, respectively. It can be understood that when the elements are other two-dimensional shapes, the contact force calculation can also be performed using the above method to obtain the smooth potential distribution.

[0106] In the second specific embodiment:

[0107] In step S1, a continuous-discontinuous coupling simulation model is constructed, model parameters are initialized, contact pair discrimination is performed, and element pairs in possible contact relationship are identified and stored, referred to as contact element pairs, and a data list of contact element pairs is established. The continuous-discontinuous coupling simulation model constructed is a three-dimensional model.

[0108] In step S2, for each contact element pair, a master element e m and a slave element e s are selected from the element ID arrangement according to the element ID arrangement. Figure 4 The slave element can be a polyhedron as shown in the figure, and the master element is a three-dimensional planar figure.

[0109] A coordinate list of each vertex of the master and slave elements is established, and a list of outer normal vectors of each plane of the master and slave elements is established.

[0110] Based on the separation axis theorem, it is determined whether the master and slave elements overlap. If any vertex x p of the master element and any vertex v j of the slave element and the outer normal vector f of each plane of the master and slave elements satisfy:

[0111] max(x p ·f)-min(v j ·f)≥0

[0112] or

[0113] max(v j ·f)-min(x p ·f)≥0

[0114] then it is considered that the master and slave elements overlap, where x p is any vertex in the coordinate list of the vertices of the master element, v j represents any vertex in the coordinate list of the vertices of the slave element, and f represents the outer normal vector of any plane of the master and slave elements.

[0115] In step S3, for the contact element pair that overlaps, the vertex coordinate x p of the master element embedded in the slave element is recorded, and this vertex is referred to as an embedded vertex.

[0116] Assuming that there are n vertices in the coordinate list of the vertices of the slave element, when establishing a weight function model of the embedded vertex x p , the weight ratio of this point to each vertex of the slave element is calculated:

[0117]

[0118] where N i represents the i-th vertex v iThe corresponding weight ratio, w i (x p ) is the vertex v i The corresponding weight, the denominator is the sum of the weights of all vertices, i = 1, 2, …, n;

[0119] Assuming that the vertex v i On the l plane, then w i (x p ) is calculated by the following formula:

[0120]

[0121] Wherein:

[0122]

[0123]

[0124] Wherein, is the shortest distance from the embedding vertex x p to the a-th plane where the vertex v i is located, is the normal vector of the a-th plane where the vertex v i is located, a = 1, 2, …, l.

[0125] Since the continuous-discontinuous coupling simulation model is a three-dimensional model, w i (x p ) is specifically calculated by the following formula:

[0126]

[0127] From the i-th vertex v i of the element, the adjacent element vertices are called adjacent vertices, and m represents the number of adjacent vertices, V i is the polyhedral volume formed by the vertex v i and the adjacent m vertices, which is decomposed into several tetrahedral volumes, with v i as the vertex and triangle f as the bottom surface The following formula is used to calculate:

[0128]

[0129] Wherein, the subscripts f1, f2, f3 are the vertex numbers of the triangle f arranged counterclockwise, and the superscripts (1), (2), (3) represent the x, y, z coordinates, and then the sum of the decomposed several tetrahedral volumes is obtained V i ;

[0130] The j-th adjacent vertex is represented by v j , j = 1, 2, …, m, Main unit embedded vertex x p Vertex v i And except for adjacent vertices v j The volume of the polyhedron formed by the other m-1 adjacent vertices is calculated using a method similar to the above formula;

[0131] w is obtained using the above method. i (x p After that, the weight ratio N is further calculated. i .

[0132] In step S4, it is understandable that, due to the weight ratio N i satisfy:

[0133]

[0134] And for any vertex of the cell, N i >0, therefore the embedding vertex x of the principal unit p The interpolation function within the element can be expressed as:

[0135]

[0136] Therefore N i That is, the main unit embeds vertex x p For the i-th vertex v in the cell i The shape function, when x p =v i There are N times i =1, therefore the shape function N i satisfy:

[0137]

[0138] Its value is maximized at the center of the polyhedron, that is... As can be seen from the above equation, the shape function gradually decreases as the embedding point approaches the boundary, and then decreases to zero at the boundary.

[0139] When calculating the potential function of embedded vertices in the element, a smooth potential function is constructed by multiplying the volumes of multiple adjacent polyhedra. The specific calculation formula is as follows:

[0140]

[0141] In step S5, based on the potential function obtained in step S4, the normal contact force generated by the slave element on the master element is calculated, and the calculated contact force is further distributed to the equivalent nodal forces of the master and slave elements according to the weighting function model (shape function). This completes the calculation of the contact force.

[0142] The steps S2-S5 can be repeated subsequently until the initialization, boundary loading, stress balancing and simulation completion of the continuous-discontinuous coupling simulation model are completed. In addition, the stress change of the whole process model is monitored to ensure that there is no obvious oscillation in the model calculation process.

[0143] Referring to Figure 4 , the calculation principle and smooth potential distribution diagram when the unit is a polyhedron are shown.

[0144] In addition, referring to Figure 5 , the model diagram of the conventional round corner method, the Munjiza method and the method provided by the present application (referred to as the smooth potential method) in the prior art are shown. The results of the friction force changing with the sliding distance calculated by using the above methods are shown in Figure 6 . It can be seen that the calculation result of the method provided by the present application is obviously closer to the theoretical value, thereby verifying the effectiveness of the present application.

[0145] In summary, the smooth contact potential calculation method of the continuous-discontinuous model of any shape provided by the present application solves various problems existing in the prior art algorithms, makes the calculation result closer to the theoretical value, is beneficial to improve the stability and precision of simulation, and improves the application range, and is suitable for the contact potential calculation of two-dimensional, three-dimensional and any geometric shape grid.

[0146] Finally, it should be pointed out that the above embodiments are only used to illustrate the technical solutions of the present application but not limit the present application. Although the present application has been described in detail with reference to the preferred embodiments, it should be understood by those skilled in the art that the technical solutions of the present application can be modified or replaced equivalently without departing from the purpose and scope of the present application, and they should be covered in the scope of the claims of the present application.

Claims

1. A method for calculating a smooth contact potential of an arbitrary shape of a continuous-discontinuous model, characterized by, The method comprises the following steps: S1, constructing a continuous-discontinuous coupling simulation model, initializing model parameters, performing contact pair discrimination, identifying and storing element pairs in possible contact relationship, referred to as contact element pairs, and establishing a data list of contact element pairs; S2, for each contact element pair, performing master and slave element allocation, and finely judging whether the contact element pair overlaps; S3, for the contact element pair that overlaps, recording the vertex coordinates of the master element embedded in the slave element, the vertex being referred to as an embedded vertex, and establishing a weight function model of the embedded vertex; S4, based on the weight function model established in step S3, directly calculating the smooth potential function of the embedded vertex in the slave element; S5, according to the smooth potential function calculated in step S4, calculating the normal contact force generated by the master element of the slave element pair, and further distributing the calculated contact force to the equivalent node forces of the master and slave elements according to the shape function; In step S3, specifically: For contact element pairs that occur in overlap, the master element records the vertex coordinates of the slave element that it is embedded in This vertex is called the embedded vertex; Assuming there are vertices in total from the coordinate list of the element vertexes, when building the weight function model of the embedded vertex , the weight ratio of the point to each vertex corresponding to the element vertex is calculated: ； wherein, represents the weight ratio of the first vertex in the cell to the corresponding weight, is the weight of the vertex in the denominator is the sum of the weights of all vertices, ; Assuming vertex In One plane or edge, then Is calculated by the following formula: ； Wherein: ； ； wherein, is the shortest distance from the vertex to the plane or edge on which the vertex is located, is the normal vector to the plane or edge on which the vertex is located, ; In step S4, when calculating the smooth potential function of the embedded vertex in the slave element, the calculation formula is as follows: 。 2. The method of claim 1, wherein, In step S1, the continuous-discontinuous coupling simulation model constructed is a two-dimensional model or a three-dimensional model of any element shape.

3. The method of claim 2, wherein the method is a method of calculating a smooth contact potential of an arbitrary shape of a continuous-discontinuous model. In step S2, specifically: For each pair of contact units, the master unit is selected from among the units arranged according to the unit ID and a slave unit ; Establishing a coordinate list of each vertex of the master and slave elements, if the continuous-discontinuous coupling simulation model is a two-dimensional model, then establishing an outer normal vector list of each edge of the master and slave elements, and if the continuous-discontinuous coupling simulation model is a three-dimensional model, then establishing an outer normal vector list of each plane of the master and slave elements; Based on the separating axis theorem, it is determined whether the master and slave units overlap, and if any vertex of the master unit , relative to any vertex of the slave unit , and the outer normal vectors of each edge or each plane of the master and slave units satisfy: ； Or ； then the master and slave cells are considered to overlap, where is any vertex in the list of master cell vertex coordinates, represents any vertex in the list of slave cell vertex coordinates, represents an outward normal to any edge or any plane of the master and slave cells.

4. The method of claim 3, wherein, In step S3, if the continuous-discontinuous coupling simulation model is a two-dimensional model, then Specifically, it is calculated by the following formula: ； wherein is the area of the triangle formed by the i-th vertex and its adjacent vertices, calculated as follows: is the area of the triangle formed by the i-th vertex and its adjacent vertices, calculated as follows: ； wherein , and are the coordinates of the , , and vertices from the list of unit vertex coordinates is the area of the triangle constructed with the first and the second vertex of the master cell and the first vertex of the slave cell is the area of the triangle constructed with the first and the second vertex of the master cell and the first vertex of the slave cell ； ； The weight ratio is further normalized by the following method After that, the weight ratio is further normalized by the following method .

5. The method of claim 3, wherein the method is a method of calculating a smooth contact potential of an arbitrary shape of a continuous-discontinuous model, characterized by, In step S3, if the continuous-discontinuous coupling simulation model is a three-dimensional model, then Specifically, it is calculated by the following formula: ； From the unit A vertex Adjacent from the unit vertex is called adjacent vertex, Indicates the number of adjacent vertices, For vertex With adjacent Polyhedral volume formed by adjacent vertex, decomposed into several tetrahedrons, wherein Vertex, triangle Single tetrahedral volume with bottom surface Calculated by the following formula: ； where the subscript , , is the counterclockwise vertex order of the triangle , and the superscript , , denotes , , the coordinates, and the sum over the decomposed tetrahedra volumes is ; with representing the first adjacent vertex, , the main unit embedding vertex , the vertex and the polyhedral volume formed by the other -1 adjacent vertices except the adjacent vertex The weight ratio is further normalized by the following method After that, the weight ratio is further normalized by the following method .

Images