Method for calculating smooth contact potential of any shape of continuous and discontinuous model

Abstract

The invention discloses a method for calculating the smooth contact potential of any shape of a continuous and discontinuous model. The method comprises the following steps: S1, constructing a continuous-discontinuous coupling simulation model, initializing model parameters, discriminating contact pairs, identifying and storing unit pairs with possible contact relations, and establishing a data list of the contact unit pairs; s2, performing master unit and slave unit distribution on each contact unit pair, and judging whether the contact unit pairs are overlapped or not; s3, for the overlapped contact unit pairs, recording coordinates of vertexes of the master unit embedded in the slave unit, the vertexes being called embedded vertexes, and establishing a weight function model of the embedded vertexes; s4, on the basis of the weight function model established in the step S3, calculating a smooth potential function of the embedded vertex in the slave unit; and S5, according to the smooth potential function calculated in the step S4, calculating a normal contact force generated by the slave unit to the master unit, and further distributing the calculated contact force to equivalent node forces of the master unit and the slave unit according to a shape function.

Description

Technical Field

[0001] The present invention relates to the technical field of geotechnical engineering numerical simulation, and in particular to a method for calculating the smooth contact potential of a continuous or discontinuous model with arbitrary shapes. Background Art

[0002] In current engineering practice, numerical calculation methods such as the finite element method (FEM), boundary element method (BEM), finite difference method (FDM), and discrete element method (DEM) are widely used in geotechnical engineering stability analysis. By constructing geomechanical models, they predict the stress-strain response of rock masses under different loading conditions. These methods provide important numerical analysis tools for engineering structure optimization design, construction safety assessment, and geological disaster early warning.

[0003] However, the existing technology has the following significant defects: traditional continuum mechanics methods are difficult to effectively simulate the dynamic evolution of cracks, the connection of crack networks and the movement behavior of blocks, and cannot fully characterize the whole process of rock fracture mechanism from continuous to discontinuous state; typical discontinuous medium analysis methods have an inherent contradiction between the accuracy of crack propagation path prediction and computational efficiency, especially under large-scale and complex working conditions, the high computational cost, poor convergence and numerical oscillation problems are particularly prominent; existing algorithms have limitations in characterizing the multi-scale structural characteristics of rock masses, and fail to effectively couple the microscopic damage evolution of materials with the macroscopic structural instability process, resulting in insufficient reliability in predicting the failure mode of complex jointed rock masses.

[0004] This technical bottleneck stems from the inherent multiphase heterogeneity, anisotropic mechanical response, and complex joint network distribution of rock materials. The failure process is essentially a progressive, cross-scale failure driven by the coupling of multiple physical fields, involving nonlinear and discontinuous mechanical behaviors such as microcrack initiation, propagation, and penetration, as well as block slip and separation.

[0005] Contact potential, a crucial step in discrete block simulation, typically occupies a significant portion of computational resources in the overall calculation of continuous-discontinuous coupling models. Traditional contact potential calculation methods face three key issues in practical applications: First, under the same embedding depth, the contact force calculation results differ significantly, leading to discontinuous mechanical responses and affecting simulation stability and accuracy. Second, when modeling with quadrilateral elements and higher-order deformation elements, conventional contact potential functions exhibit contact force singularities or unreasonable oscillations in certain characteristic regions, severely impacting the continuity and rationality of the crack propagation path. Third, existing contact potential models are only applicable to two-dimensional triangular meshes and three-dimensional tetrahedrons.

[0006] Therefore, it is necessary to design a new method for calculating the smooth contact potential of arbitrary unit shapes to overcome the above problems. Summary of the Invention

[0007] To achieve the above object, the technical solution adopted by the present invention is:

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

[0009] S1. Build a continuous-discontinuous coupling simulation model, initialize model parameters, perform contact pair discrimination, identify and store unit pairs with possible contact relationships, called contact unit pairs, and establish a data list of contact unit pairs;

[0010] S2. For each contact unit pair, assign the master and slave units, and carefully determine whether the contact unit pairs overlap;

[0011] S3. For the overlapping contact unit pairs, record the vertex coordinates of the master unit embedded in the slave unit, which is called the embedded vertex, and establish a weight function model for the embedded vertex;

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

[0013] S5. Calculate the normal contact force generated by the slave unit on the master unit according to the smooth potential function calculated in step S4, and further distribute the calculated contact force as equivalent nodal forces of the master and slave units according to the shape function.

[0014] Furthermore, in step S1 , the constructed continuous-discontinuous coupling simulation model is a two-dimensional model or a three-dimensional model of arbitrary unit shapes.

[0015] Furthermore, in step S2, specifically:

[0016] For each contact element pair, select the main element e according to the element ID arrangement. m and a slave unit e s ;

[0017] Establish a coordinate list of each vertex of the master and slave units. If the continuous-discontinuous coupling simulation model is a two-dimensional model, establish a list of the external normal vectors of each edge of the master and slave units. If the continuous-discontinuous coupling simulation model is a three-dimensional model, establish a list of the external normal vectors of each plane of the master and slave units.

[0018] Based on the separating axis theorem, we can judge whether the master and slave units overlap. If any vertex x of the master unit exists p , relative to any vertex v on the cell j And the external normal vector f of each side or each plane of the master and slave units satisfies:

[0019] max(x p f)-min(x jf)≥0

[0020] or

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

[0022] It is considered that the master and slave units overlap, where x p is any vertex in the main unit vertex coordinate list, v j Represents any vertex in the slave unit vertex coordinate list, and f represents the external normal vector of any edge or plane of the master or slave unit.

[0023] Furthermore, in step S3, specifically:

[0024] For overlapping contact element pairs, record the vertex coordinates x where the master element is embedded in the slave element. p , this vertex is called an embedded vertex;

[0025] Assuming that there are n vertices in the coordinate list of the unit vertex, establish the embedded vertex x p When using the weight function model, calculate the weight ratio of the point to each vertex of the slave unit:

[0026]

[0027] Among them, N i Indicates the i-th vertex v in the cell i The 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;

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

[0029]

[0030] in:

[0031]

[0032]

[0033] in, is the embedded vertex x p To the vertex v i The shortest distance to the a-th plane or edge, For vertex v iThe normal vector of the a-th plane or edge, 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 ) is specifically calculated by the following formula:

[0035]

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

[0037]

[0038] where v i-1 、v i and v i+1 They are respectively from the unit vertex coordinate list V s The coordinates of the i-1, i and i+1 vertices in ;

[0039] Indicates that the main unit is embedded in the vertex x p and the area of ​​the triangle constructed from the i-1th and ith vertices in the cell, The main unit is embedded in the vertex x p The area of ​​the triangle constructed from the i-th and i+1-th vertices in the cell is calculated as follows:

[0040]

[0041]

[0042] By the above method, we can get w i (x p ), further normalization calculation is performed to obtain the weight ratio N i .

[0043] Furthermore, in step S3, if the continuous-discontinuous coupling simulation model is a three-dimensional model, then w i (x p ) is specifically calculated by the following formula:

[0044]

[0045] From the i-th vertex v i Adjacent vertices from a unit 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 the adjacent m vertices is decomposed into several tetrahedrons for solution, with v iThe volume of a single tetrahedron with vertices and triangle f as the base Calculated by the following formula:

[0046]

[0047] The subscripts f1, f2, and f3 are the vertex numbers of the triangle f arranged counterclockwise, and the superscripts (1), (2), and (3) represent the x, y, and z coordinates. Then, the sum of the volumes of the decomposed tetrahedrons is obtained to obtain V i ;

[0048] Use v j represents the jth adjacent vertex, j = 1, 2…, m, Embed vertex x into the main unit p , vertex v i and except for the adjacent vertex v j The volume of the polyhedron formed by the other m-1 adjacent vertices;

[0049] By the above method, we can get w i (x p ), further normalization calculation is performed to obtain the weight ratio N i .

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

[0051]

[0052] Compared with the existing technology, the method for calculating the smooth contact potential of arbitrary shapes of continuous and discontinuous models provided by the present invention solves the various problems of non-smoothness, high computational complexity, non-conservative force, and inability to adapt to arbitrary unit shapes in the existing algorithms, making the calculation results closer to the theoretical values, which is beneficial to improving the stability and accuracy of the simulation, and expanding the scope of application. It is suitable for contact potential calculation of two-dimensional and three-dimensional grids of arbitrary geometric shapes. BRIEF DESCRIPTION OF THE DRAWINGS

[0053] Figure 1 A schematic flow chart of the method for calculating the smooth contact potential of an arbitrary shape of a continuous or discontinuous model provided by the present invention;

[0054] Figure 2 Schematic diagram of the smooth potential distribution of triangles from the unit;

[0055] Figure 3 Schematic diagram of the smooth potential distribution of quadrilateral from the unit;

[0056] Figure 4 Schematic diagram of the multimodal slave unit smooth potential distribution;

[0057] Figure 5 Schematic diagram of the models of the conventional fillet method, the Munjiza method and the method provided by the present invention;

[0058] Figure 6 Schematic diagram of the friction force calculated by the conventional fillet method, Munjiza method and the method provided by the present invention as a function of sliding distance. DETAILED DESCRIPTION

[0059] In order to make the technical means, creative features, objectives and effects achieved by the present invention easier to understand, the following further describes how the present invention is implemented in conjunction with specific implementation methods.

[0060] Reference Figure 1 As shown, the present invention provides a method for calculating the smooth contact potential of an arbitrary shape of a continuous or discontinuous model, comprising the following steps:

[0061] S1. Build a continuous-discontinuous coupling simulation model, initialize model parameters, perform contact pair discrimination, identify and store unit pairs with possible contact relationships, called contact unit pairs, and establish a data list of contact unit pairs;

[0062] S2. For each contact unit pair, assign the master and slave units, and carefully determine whether the contact unit pairs overlap;

[0063] S3. For the overlapping contact unit pairs, record the vertex coordinates of the master unit embedded in the slave unit, which is called the embedded vertex, and establish a weight function model for the embedded vertex;

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

[0065] S5. Calculate the normal contact force generated by the slave unit on the master unit according to the smooth potential function calculated in step S4, and further distribute the calculated contact force as equivalent nodal forces of the master and slave units according to the shape function.

[0066] In the first specific embodiment:

[0067] In step S1, a continuous-discontinuous coupling simulation model is constructed, model parameters are initialized, contact pairs are identified, pairs of elements with potential contact relationships are identified and stored, and a data list of contact element pairs is created. The constructed continuous-discontinuous coupling simulation model is a two-dimensional model or a three-dimensional model.

[0068] In step S2, for each contact unit pair, the main unit e is selected according to the unit ID arrangement. m and a slave unit e s ;reference Figure 2 and Figure 3 As shown, the slave unit can be a triangle, quadrilateral or other shape, and the master unit is a two-dimensional plane figure;

[0069] Create a list of coordinates of each vertex of the master and slave units, and a list of external normal vectors of each edge of the master and slave units;

[0070] Based on the separating axis theorem, we can judge whether the master and slave units overlap. If any vertex x of the master unit exists p , relative to any vertex v on the cell j And the external normal vector f on each side of the master and slave units satisfies:

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

[0072] or

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

[0074] It is considered that the master and slave units overlap, where x p is any vertex in the main unit vertex coordinate list, v j Represents any vertex in the slave unit vertex coordinate list, and f represents the external normal vector of any side of the master or slave unit.

[0075] In step S3, for the overlapping contact unit pairs, record the vertex coordinates x of the master unit embedded in the slave unit. p , this vertex is called an embedded vertex;

[0076] Assuming that there are n vertices in the coordinate list of the unit vertex, establish the embedded vertex x p When using the weight function model, calculate the weight ratio of the point to each vertex of the slave unit:

[0077]

[0078] Among them, N i Indicates the i-th vertex v in the cell i The 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;

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

[0080]

[0081] in:

[0082]

[0083]

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

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

[0086]

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

[0088]

[0089] where v i-1 、v i and v i+1 They are respectively from the unit vertex coordinate list V s The coordinates of the i-1, i and i+1 vertices in ;

[0090] Indicates that the main unit is embedded in the vertex x p and the area of ​​the triangle constructed from the i-1th and ith vertices in the cell, The main unit is embedded in the vertex x p The area of ​​the triangle constructed from the i-th and i+1-th vertices in the cell is calculated as follows:

[0091]

[0092]

[0093] By the above method, we can get w i (x p ), further normalization calculation is performed to obtain the weight ratio N i .

[0094] In step S4, it can be understood that due to the weight ratio N i satisfy:

[0095]

[0096] And for any vertex from the unit, there are N i > 0, so the embedded vertex x of the main unit p The interpolation function within the slave unit can be expressed as:

[0097]

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

[0099]

[0100] Its value satisfies the maximum at the center, that is It can be seen from the above formula that as the embedded point approaches the boundary, the shape function gradually decreases and decreases to zero at the boundary.

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

[0102]

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

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

[0105] Reference Figure 2 and Figure 3 As shown in the figure, the calculation principle and smooth potential distribution diagram are respectively shown when the slave unit is a triangle and a quadrilateral. It can be understood that when the slave unit is a two-dimensional figure of other shapes, the contact force calculation can also be completed according to 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 pairs are identified, pairs of elements with potential contact relationships are identified and stored, and a data list of contact element pairs is created. The constructed continuous-discontinuous coupling simulation model is a three-dimensional model.

[0108] In step S2, for each contact unit pair, the main unit e is selected according to the unit ID arrangement. m and a slave unit e s ; The slave unit can be Figure 4 The polyhedron shown has a main unit that is a three-dimensional plane figure;

[0109] Establish a list of coordinates of each vertex of the master and slave units, and establish a list of external normal vectors of each plane of the master and slave units;

[0110] Based on the separating axis theorem, we can judge whether the master and slave units overlap. If any vertex x of the master unit exists p , relative to any vertex v on the cell j And the external normal vector f of each plane of the master and slave units satisfies:

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

[0112] or

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

[0114] It is considered that the master and slave units overlap, where x p is any vertex in the main unit vertex coordinate list, v j Represents any vertex in the slave unit vertex coordinate list, and f represents the external normal vector of any plane of the master or slave unit.

[0115] In step S3, for the overlapping contact unit pairs, record the vertex coordinates x of the master unit embedded in the slave unit. p , this vertex is called an embedded vertex;

[0116] Assuming that there are n vertices in the coordinate list of the unit vertex, establish the embedded vertex x p When using the weight function model, calculate the weight ratio of the point to each vertex of the slave unit:

[0117]

[0118] Among them, N i Indicates the i-th vertex v in the cell 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] Assume that vertex v i On l planes, then w i (x p ) is calculated by the following formula:

[0120]

[0121] in:

[0122]

[0123]

[0124] in, is the embedded vertex x p To the vertex v i The shortest distance to the a-th plane, For vertex v i The normal vector of the a-th plane, 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 Adjacent vertices from a unit are called adjacent vertices, m represents the number of adjacent vertices, V i For vertex v i The volume of the polyhedron formed by the adjacent m vertices is decomposed into several tetrahedrons for solution, with v i The volume of a single tetrahedron with vertices and triangle f as the base Calculated by the following formula:

[0128]

[0129] The subscripts f1, f2, and f3 are the vertex numbers of the triangle f arranged counterclockwise, and the superscripts (1), (2), and (3) represent the x, y, and z coordinates. Then, the sum of the volumes of the decomposed tetrahedrons is obtained to obtain V i ;

[0130] Use v j represents the jth adjacent vertex, j = 1, 2…, m, Embed vertex x into the main unit p , vertex v i and except for the adjacent vertex 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] By the above method, we can get w i (x p ), further calculation is performed to obtain the weight ratio N i .

[0132] In step S4, it can be understood that due to the weight ratio N i satisfy:

[0133]

[0134] And for any vertex from the unit, there are N i > 0, so the embedded vertex x of the main unit p The interpolation function within the slave unit can be expressed as:

[0135]

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

[0137]

[0138] Its value satisfies the maximum at the center of the polyhedron, that is, It can be seen from the above formula that as the embedded point approaches the boundary, the shape function gradually decreases and decreases to zero at the boundary.

[0139] When calculating the potential function of the embedded vertex in the slave unit, a smooth potential energy function is constructed by multiplying the volumes of multiple adjacent polyhedrons. The specific calculation formula is as follows:

[0140]

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

[0142] Steps S2-S5 can be repeated until the initialization, boundary loading, stress balance, and simulation of the continuous-discontinuous coupled simulation model are complete. Furthermore, the stress changes in the model are monitored throughout the entire process to ensure that there are no significant oscillations during the model calculation.

[0143] Reference Figure 4 As shown in the figure, it is a schematic diagram of the calculation principle and smooth potential distribution when the unit is a polyhedron.

[0144] In addition, refer to Figure 5 The following is a schematic diagram of the model of the conventional fillet method, Munjiza method and the method provided by the present invention (called smooth potential method) in the prior art. The results of the friction force calculated by the above methods as a function of sliding distance are shown as follows: Figure 6 As shown, it can be seen that the calculation results of the method provided by the present invention are obviously closer to the theoretical value, thereby verifying the effectiveness of the present invention.

[0145] In summary, the method for calculating the smooth contact potential of arbitrary shapes of continuous and discontinuous models provided by the present invention solves various problems existing in existing algorithms, makes the calculation results closer to the theoretical values, is conducive to improving the stability and accuracy of the simulation, and expands the scope of application. It is suitable for the contact potential calculation of two-dimensional and three-dimensional grids of arbitrary geometric shapes.

[0146] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not limiting. Although the present invention has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that the technical solutions of the present invention may be modified or replaced by equivalents without departing from the purpose and scope of the technical solutions of the present invention, which should all be covered by the scope of the claims of the present invention.

Claims

1. A method for calculating the smooth contact potential of an arbitrary shape of a continuous or discontinuous model, characterized in that: The steps include: S1. Build a continuous-discontinuous coupling simulation model, initialize model parameters, perform contact pair discrimination, identify and store unit pairs with possible contact relationships, called contact unit pairs, and establish a data list of contact unit pairs; S2. For each contact unit pair, assign the master and slave units, and carefully determine whether the contact unit pairs overlap; S3. For the overlapping contact unit pairs, record the vertex coordinates of the master unit embedded in the slave unit, which is called the embedded vertex, and establish a weight function model for 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 unit; S5. Calculate the normal contact force generated by the slave unit on the master unit according to the smooth potential function calculated in step S4, and further distribute the calculated contact force as equivalent nodal forces of the master and slave units according to the shape function.

2. The method for calculating the smooth contact potential of an arbitrary shape of a continuous or discontinuous model according to claim 1, characterized in that: In step S1 , the constructed continuous-discontinuous coupling simulation model is a two-dimensional model or a three-dimensional model of arbitrary unit shapes.

3. The method for calculating the smooth contact potential of an arbitrary shape of a continuous or discontinuous model according to claim 2, characterized in that: In step S2, specifically: For each contact element pair, select the main element e according to the element ID arrangement. m and a slave unit e s ; Establish a coordinate list of each vertex of the master and slave units. If the continuous-discontinuous coupling simulation model is a two-dimensional model, establish a list of the external normal vectors of each edge of the master and slave units. If the continuous-discontinuous coupling simulation model is a three-dimensional model, establish a list of the external normal vectors of each plane of the master and slave units. Based on the separating axis theorem, we can judge whether the master and slave units overlap. If any vertex x of the master unit exists p , relative to any vertex v on the cell j And the external normal vector f of each side or each plane of the master and slave units satisfies: max(x p ·f)-min(v j ·f)≥0 or max(v j ·f)-min(x p ·f)≥0 It is considered that the master and slave units overlap, where x p is any vertex in the main unit vertex coordinate list, v j Represents any vertex in the slave unit vertex coordinate list, and f represents the external normal vector of any edge or plane of the master or slave unit.

4. The method for calculating the smooth contact potential of an arbitrary shape of a continuous or discontinuous model according to claim 3, characterized in that: In step S3, specifically: For overlapping contact element pairs, record the vertex coordinates x where the master element is embedded in the slave element. p , this vertex is called an embedded vertex; Assuming that there are n vertices in the coordinate list of the unit vertex, establish the embedded vertex x p When using the weight function model, calculate the weight ratio of the point to each vertex of the slave unit: Among them, N i Indicates the i-th vertex v in the cell i The 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; Assume that vertex v i On l planes or edges, then w i (x p ) is calculated by the following formula: in: in, is the embedded vertex x p To the vertex v i The shortest distance to the a-th plane or edge, For vertex v i The normal vector of the a-th plane or edge, a=1,2…,l.

5. The method for calculating the smooth contact potential of an arbitrary shape of a continuous or discontinuous model according to claim 4, characterized in that: In step S3, if the continuous-discontinuous coupling simulation model is a two-dimensional model, then w i (x p ) is specifically calculated by the following formula: Among them A i It is the area of ​​the triangle formed by the i-th vertex of the unit and its adjacent vertices. The calculation formula is as follows: where v i-1 、v i and v i+1 They are respectively from the unit vertex coordinate list V s The coordinates of the i-1, i and i+1 vertices in ; Indicates that the main unit is embedded in the vertex x p and the area of ​​the triangle constructed from the i-1th and ith vertices in the cell, The main unit is embedded in the vertex x p The area of ​​the triangle constructed from the i-th and i+1-th vertices in the cell is calculated as follows: By the above method, we can get w i (x p ), further normalization calculation is performed to obtain the weight ratio N i .

6. The method for calculating the smooth contact potential of an arbitrary shape of a continuous or discontinuous model according to claim 4, characterized in that: In step S3, if the continuous-discontinuous coupling simulation model is a three-dimensional model, then w i (x p ) is specifically calculated by the following formula: From the i-th vertex v i Adjacent vertices from a unit are called adjacent vertices, m represents the number of adjacent vertices, V i For vertex v i The volume of the polyhedron formed by the adjacent m vertices is decomposed into several tetrahedrons for solution, where v i The volume of a single tetrahedron with vertices and triangle f as the base Calculated by the following formula: The subscripts f1, f2, and f3 are the vertex numbers of the triangle f arranged counterclockwise, and the superscripts (1), (2), and (3) represent the x, y, and z coordinates. Then, the sum of the volumes of the decomposed tetrahedrons is obtained to obtain v i ; Use v j represents the jth adjacent vertex, j = 1, 2…, m, Embed vertex x into the main unit p , vertex v i and except for the adjacent vertex v j The volume of the polyhedron formed by the other m-1 adjacent vertices; By the above method, we can get w i (x p ), further normalization calculation is performed to obtain the weight ratio N i .

7. The method for calculating the smooth contact potential of an arbitrary shape of a continuous or discontinuous model according to claim 4, characterized in that: In step S4, when calculating the normalized potential function of embedding the vertex in the slave unit, the calculation formula is as follows:

Images