A surface-embedded material point method for multi-body contact problems

Abstract

The application discloses a surface-embedded material point method for a multi-body contact problem, and belongs to the technical field of computational mechanics and numerical simulation. In a first step, a physical model of a multi-body contact system is discretized and parameters are set. In a second step, contact forces are calculated. In a third step, physical quantities of material points of contact bodies are transferred to background grid nodes to obtain node mass, momentum and force. In a fourth step, a grid is updated. In a fifth step, velocities, displacements, stresses and densities of the material points are updated. In a sixth step, background grid physical quantities are reset. In a seventh step, the first, fourth, fifth and sixth steps are repeated to update physical quantities of material points of other contact bodies and to allocate corresponding contact body background grid nodes in a memory. In an eighth step, after all contact bodies are updated, a current time is updated to obtain a dynamic simulation result. The application can accurately represent object boundaries, effectively simulate complex problems such as impact and multi-body contact, and has good stability and low memory consumption.

Description

Technical Field

[0001] This invention belongs to the field of computational mechanics and numerical simulation technology, and relates to a surface-embedded material point method for multibody contact problems, and more particularly to a high-precision numerical simulation method for large deformation multibody contact problems based on the material point method. Background Technology

[0002] In the field of engineering, multibody contact phenomena (such as collisions, impacts, and fluid-structure interactions) are widespread, and their numerical simulation is crucial for engineering design and analysis. Among traditional numerical methods, the finite element method (FEM) can accurately describe the contact interface through mesh boundaries, but it is prone to mesh distortion under extreme deformation scenarios, leading to unstable internal force updates and reduced time steps, which seriously affects computational efficiency and accuracy.

[0003] The Material Point Method (MPM), as a typical meshless method, combines the advantages of Lagrangian and Eulerian descriptions and performs well in large deformation problems (such as metal plastic forming, geotechnical dynamic response, and ultra-high-speed collision). However, traditional MPM has significant drawbacks in contact problems: First, the interface representation is fuzzy, relying on the background mesh for contact detection. Due to the mismatch between the physical boundary and the mesh boundary, premature contact or non-physical penetration is prone to occur. Second, the contact normal vector is approximated by the mesh quality gradient, which has poor accuracy when the mesh is coarse or the material points are unevenly distributed. Third, in multi-body contact, background mesh memory needs to be allocated for all contact bodies, resulting in high memory consumption. Fourth, contact calculation based on the background mesh requires that the background mesh size and node position of each contact body be consistent, making it impossible to set a suitable background mesh size for each contact body.

[0004] Among existing improved methods, background mesh contact algorithms based on physical distance-assisted detection can alleviate the problem of premature contact, but they still cannot solve the problem of non-physical penetration at the contact interface. Implicit interface methods based on level sets are complex to implement and computationally expensive in 3D multi-body systems. In addition, optimization schemes for contact detection logic also have limitations: for example, an MPM contact algorithm disclosed in Chinese invention patent (publication number CN115410663A) directly calculates the contact force based on the position of the material points, which is equivalent to migrating the "pinball method" to the MPM framework. It achieves contact judgment by setting the radius of the virtual collision ball between material points, effectively avoiding the problems caused by traditional mesh dependence. While this method addresses the issue of spurious contact, it suffers from low accuracy in representing complex topological contact interfaces (such as components with grooves or sharp corners). The setting of the virtual sphere radius can easily lead to uneven distribution of contact force. Another typical improvement scheme is the application of the NTS (Node-To-Segment) point-to-surface contact algorithm in MPM proposed in Chinese invention patent (publication number CN114357717A). This algorithm optimizes contact detection accuracy by establishing a mapping relationship between material points and line segments on the contact body surface. However, it uses a fixed penalty coefficient to solve the contact constraint. Under extreme conditions such as high-speed impact and high strain rate, the fixed penalty coefficient is difficult to match the dynamic contact stiffness requirements, which can easily lead to non-physical penetration of the contact interface.

[0005] In summary, existing MPM contact frameworks are either limited by the accuracy of complex interfaces, difficult to adapt to extreme working conditions, or have a contradiction between implementation complexity and computational cost. There is an urgent need for a multibody contact numerical method that can accurately characterize boundaries, efficiently handle contact constraints, and consume little memory. Summary of the Invention

[0006] To address the problems of fuzzy interface representation, inaccurate contact detection, easy non-physical penetration, mesh distortion affecting stability, and high memory consumption in traditional MPM and existing improved methods for multi-body contact problems, this invention provides a surface embedded material point method (SEMPM) for multi-body contact problems, which can achieve high-precision and high-stability numerical simulation of multi-body contact scenarios.

[0007] To achieve the above objectives, the present invention adopts the following technical solution:

[0008] A method for embedding material points on a surface for multi-body contact problems includes the following steps:

[0009] The first step is initialization. The physical model of the multi-body contact system is discretized and its parameters are set to provide the foundational data for subsequent contact calculations and object updates. Specifically:

[0010] Step 1.1, Discretization. The internal domain of each contact body is discretized into material points to carry physical information such as mass and stress. Here, it is defined that there are contact bodies A, B, C, D, and other contact bodies in contact. The boundaries of each contact body are discretized into massless master patches, which are only used for geometric boundary representation and do not carry physical quantities. Slave nodes that match the master patches are defined as reference points for contact detection. The computational domain space is discretized into a regular structured background mesh, and memory is allocated to the background mesh nodes in the region where contact body A is located. The background mesh size is twice the average particle spacing of the discretized material points.

[0011] Step 1.2, Parameter Setting. Define the material parameters, contact parameters, and calculation parameters for each contact element. The material parameters include elastic modulus, Poisson's ratio, density, yield stress, etc., and the contact parameters include contact gap. Main surface area Parameters to be calculated include time step. Total simulation time, etc.

[0012] The second step is contact force calculation. At the beginning of each time step, actual contact pairs are identified through contact search. For these pairs, a variable coefficient penalty stiffness mechanism is used to calculate the contact force, ensuring the reliability of the contact constraints. Specifically:

[0013] Step 2.1, Contact Search. Perform geometric detection on the master face and slave node to determine if they form a true contact pair; specifically: calculate the minimum distance from the slave node to the master face. If the minimum distance is less than the set contact gap... If the detected slave node-master facet contact pair is not a true contact pair, then the cumulative tangential displacement of the contact pair is determined. Set to zero.

[0014] Furthermore, the method for calculating the minimum distance from a node to a main face is as follows: Assuming point... It is the main face. Upper and lower node positions The point with the smallest distance should satisfy a perpendicular relationship:

[0015] (1)

[0016] In the formula, These are the natural coordinates of a point on the main facet. Represents the surface function of the principal patch; Indicates the location of the contact point.

[0017] For a four-node master face, based on finite element discretization, it can be written in the following form:

[0018] (2)

[0019] In the formula, Indicates the node number of the main facet; Represents the main face node Spatial location; This indicates that for the main face node Finite element shape functions.

[0020] The location of the contact point can be obtained from formula (2). When the contact point is located When the contact point is outside the main facet, theoretically the minimum distance between the main facet and the slave node is at the edge of the main facet, requiring the following correction:

[0021] (3)

[0022] in, and Indicates the natural coordinates of the contact point on the main facet;

[0023] Get the location of the contact point Then, calculate the minimum distance from the node to the main face. If the minimum distance satisfies This means that the contact pair is considered a true contact pair. The aforementioned contact gap... Take 1 / 10 of the average side length of the main face.

[0024] Step 2.2, Contact Penetration Calculation. For actual contact pairs, calculate the penetration amount from the node relative to the master face. ,in ;

[0025] Step 2.3, Solving for contact forces. For the actual contact pairs, calculate the main face nodes based on the variable coefficient penalty method. Contact force and contact force at the node Specifically:

[0026] First, calculate the normal contact force between the contact pairs. The formula is as follows:

[0027] (4)

[0028] In the formula, The area of ​​the main face; The bulk modulus of the material point of the main contact body; The initial volume of the material point of the main contact body; Normal to the main facet; This refers to the contact gap.

[0029] Secondly, the contact force relative to the normal direction is obtained. Then, the contact friction force was calculated based on the Coulomb friction model. The formula is as follows:

[0030] (5)

[0031] In the formula, For the cumulative tangential displacement of the contact pair; Coulomb friction coefficient; Let be the tangential force direction vector, according to We obtain the formula. and These are the velocity vectors from the node and the contact point of the main facet, respectively, where, according to finite element interpolation, the velocity vector of the contact point of the main facet is... , Main face node The velocity vector.

[0032] Finally, the main face node Contact force From the nodal contact force .

[0033] Step 2.4, Update the cumulative tangential displacement of the contact pair. For the actual contact pair, sum the contact tangential displacements from this step. The formula is as follows:

[0034] (6)

[0035] In the formula, For time step.

[0036] The third step involves transferring the physical quantities of a contacting body's material point to the background mesh nodes, obtaining the mass, momentum, and force of the background mesh nodes. Specifically:

[0037] Step 3.1, Mass Transfer. The mass of the main contact body material point is transferred. The formula for passing this to the background mesh node is as follows:

[0038] (7)

[0039] In the formula, Let I be the quality of the background mesh node I at time n. The position of the matter point at time n The weight of background grid node I; the subscript p indicates the material point number.

[0040] Step 3.2, Force Transfer. The external forces and stresses at each material point are applied to the background mesh nodes, and the contact forces between the master and slave nodes are transferred to the background mesh nodes, as shown in the following formula:

[0041] (8)

[0042] In the formula, The force vector of the mesh node I in the background mesh; External force at a point on matter; Let V be the volume of a point mass. Let Cauchy stress tensor be the material point; The gradient of the interpolation function; and These represent the contact force and position of node k on the main facet, respectively. The position of the matter point at time n Weights for background grid node I.

[0043] Step 3.3, Velocity Transfer. The velocity of each material point is transferred to the background mesh nodes, using the following formula:

[0044] (9)

[0045] In the formula, Let I be the velocity vector of the background mesh node I at time n; Let n be the velocity vector of the material point at time n; Let p be the mass of the substance point.

[0046] Step 4: Mesh Update. Based on the forces and masses of the background mesh nodes, update the nodal velocities to complete the dynamic calculations at the mesh level. Specifically:

[0047] Step 4.1, speed update. According to Newton's second law, the force at the mesh nodes... With node quality Calculate the node speed increment and update the node speed using the following formula:

[0048] (10)

[0049] In the formula, Let I be the velocity vector of background grid node I at time n+1; Let I be the velocity of background grid node I at time n; For time step.

[0050] Step 5: Update the velocity, displacement, stress, and density of the material points based on the velocities of the background mesh nodes. Specifically:

[0051] Step 5.1, update the position of the substance point, using the following formula:

[0052] (11)

[0053] In the formula, Let n be the position of the material point at time n+1.

[0054] Step 5.2, update the position of surface nodes, where surface nodes refer to master patch nodes and slave nodes, as shown in the following formula:

[0055] (12)

[0056] In the formula, Let n be the position of the surface node at time n+1.

[0057] Step 5.3, update the velocity of the matter point, the formula is as follows:

[0058] (13)

[0059] In the formula, Let be the velocity of the material point at time n+1.

[0060] Step 5.4, update the velocity of the surface nodes, using the following formula:

[0061] (14)

[0062] In the formula, Let n be the velocity of the surface node at time n+1.

[0063] Step 5.5, Update Material Point Stress and Density. Update the material point density and update the material point stress according to the material constitutive model. The specific formula is as follows:

[0064] First, calculate the strain increment at the material point. :

[0065] (15)

[0066] (16)

[0067] in, This represents the velocity gradient of the material point at time n+1; This represents the gradient of the interpolation function of the material point p with respect to the background grid node I;

[0068] Calculate the density of the material point at time n+1. Cauchy stress :

[0069] (17)

[0070] (18)

[0071] in, Let p represent the density of the matter point at time n+1; Let p represent the density of the matter point at time n; The spherical strain increment represents the change in material point p. This represents the Cauchy stress at material point p at time n+1; This represents the Cauchy stress at point p of matter at time n; This represents the Cauchy stress increment at material point p;

[0072] The specific calculation method for the Cauchy stress at the material point depends on the constitutive model of the material used.

[0073] Step 6: Reset the physical quantities of the background mesh. Set the mass, velocity, and force of all background mesh nodes to zero.

[0074] Step 7: Repeat steps 3, 4, 5, and 6 to update the physical quantities of the material points of other contact bodies, as well as the positions and velocities of the slave nodes. Allocate background grid node memory for the corresponding contact body region and reuse the background grid node memory allocated in the previous contact body update.

[0075] Step 8: After all contacts have completed the update, update the current time: In the formula Let n be the time at time n+1. If the simulation ends before the timeout, return to step two; otherwise, exit the simulation, complete the solution, and obtain the dynamic simulation results.

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

[0077] (1) The contact boundary is accurately represented in the calculation process of the present invention (corresponding to the first step): the object boundary is discretized by massless master facets to form a Lagrange surface mesh that fits the object, thus solving the problem of mismatch between the physical boundary and the mesh boundary in traditional MPM.

[0078] (2) The contact constraint is reliable in the calculation process of the present invention (corresponding to the second step): a variable coefficient penalty stiffness is introduced, and the contact stiffness increases with the increase of the penetration amount, which can dynamically match the contact stiffness requirements under extreme working conditions.

[0079] (3) The memory consumption during the calculation process of the present invention is low (corresponding to step 7): the sequential update and mesh reuse strategy is adopted, and the background mesh only needs to allocate memory for the "number of nodes required for the maximum contact body", rather than the "total number of nodes of all contact bodies" of the traditional MPM.

[0080] (4) The mesh consistency requirement is low during the calculation process of the present invention (corresponding to steps 3-7): the contact calculation is based on the geometric relationship between the master facet and the slave node, decoupled from the background mesh, allowing each contact body to adopt a mesh size suitable for itself. Attached Figure Description

[0081] Figure 1This is the flowchart of the algorithm of this invention.

[0082] Figure 2 This is a schematic diagram of the geometric model of a multibody large deformation contact example.

[0083] Figure 3 This is a comparison diagram of the contact process of various numerical methods in a multibody large deformation contact example; Figure 3 (a) in the figure is a structural deformation diagram at 0.5 seconds in this invention; Figure 3 (b) in the figure is the structural deformation diagram at 0.75 seconds of the present invention; Figure 3 (c) in the figure is the structural deformation diagram at 1.0 second of the present invention; Figure 3 (d) in the figure represents the structural deformation at 0.5 seconds using the finite element method. Figure 3 (e) in the figure is the structural deformation diagram at 0.75 seconds in the finite element method. Figure 3 (f) in the figure is the structural deformation diagram at 1.0 second using the finite element method; Figure 3 (g) in the diagram represents the structural deformation at 0.5 seconds using the traditional material point method. Figure 3 (h) in the diagram represents the structural deformation at 0.75 seconds using the traditional material point method. Figure 3 (i) in the diagram represents the structural deformation at 1.0 seconds using the traditional material point method. Detailed Implementation

[0084] The present invention will be further illustrated below with reference to specific implementation examples.

[0085] The following tests several elasto-plastic extrusion examples of rigid bodies, with the specific geometric models and example settings as follows: Figure 2 In this embodiment, among contact bodies A, B, C, and D, contact body A is a rigid disk, while contact bodies B, C, and D are all elasto-plastic bodies. The rigid disk has a mass of 30,000 kg and a diameter of 4 m, and operates under gravity g = 9.81 m / s². 2 Under the influence of force, the elastoplastic body is allowed to fall freely. Initially, the distance between points B, C, and D is 0.1m, with point D placed on a smooth plane. Points B, C, and D are each 2m high, and their widths are 6m, 8m, and 10m, respectively. The elastoplastic body adopts a linearly hardened constitutive model, with a Young's modulus of 10 MPa, a Poisson's ratio of 0.3, and a density of... The yield stress is 0.157 MPa and the strain hardening coefficient is 0.425 MPa.

[0086] A method for surface-embedded material points for multi-body contact problems, such as Figure 1 As shown, it includes the following steps:

[0087] The first step is initialization. The physical model of the multi-body contact system is discretized and its parameters are set to provide basic data for subsequent contact calculations and object updates.

[0088] Step 1.1, Discretization. The interior of each contact body is discretized into 13326 material points with an average material point spacing of 0.1m and a background mesh spacing of 0.2m. The geometric surface of the contact body is discretized into 626 quadrilateral unstructured meshes as the contact surface (all surface meshes are defined as master patches and all surface mesh nodes are defined as slave nodes).

[0089] Step 1.2, Parameter Setting. Initialize the material parameters and contact gaps for each contact point based on the calculation example's material settings. =0.01m, time step The total simulation time is 1 second.

[0090] The second step is contact force calculation (inputting the positions of the master and slave nodes to obtain the contact forces of the master and slave nodes). At the beginning of each time step, real contact pairs are identified through contact search. For real contact pairs, the contact force is calculated based on a variable coefficient penalty stiffness mechanism to ensure reliable contact constraints.

[0091] Step 2.1, Contact Search. Perform geometric detection on the master face and slave node to determine if they form a true contact pair; specifically: calculate the minimum distance from the slave node to the master face. If the minimum distance is less than the set contact gap... If so, then the "from node-main face" is determined to be a real contact pair;

[0092] In this embodiment, the minimum distance from a node to the main face is calculated as follows: Assuming the point... It is the main face. Upper and lower node positions The point with the smallest distance should satisfy a perpendicular relationship:

[0093] (1)

[0094] In the formula These are the natural coordinates of a point on the main facet. The main facet, based on finite element discretization, can be written in the following form:

[0095] (2)

[0096] In the formula It is the main face node Spatial location. This refers to the main face node. Finite element shape functions.

[0097] The location of the contact point can be obtained from formula (2). For the obtained natural coordinates of the contact point When the contact point is outside the facet, theoretically the minimum distance between the main facet and the slave node is at the edge of the main facet, requiring the following correction:

[0098] (3)

[0099] Get the location of the contact point Then, calculate the minimum distance from the node to the main face. If the minimum distance satisfies This means that the contact pair is considered a real contact pair in this embodiment. The value is set to 0.01m. If the detected "slave node-master face" contact pair is not a true contact pair, then the cumulative tangential displacement of the contact pair is... Set to zero.

[0100] Step 2.2, Contact Penetration Calculation. For actual contact pairs, calculate the penetration amount from the node relative to the master face. ( );

[0101] Step 2.3, Solving for contact forces. For the actual contact pairs, calculate the principal face nodes based on the variable coefficient penalty method mechanism. Contact force and contact force at the node First, calculate the normal contact force between the contact pairs. The formula is as follows:

[0102] (4)

[0103] In the formula The area of ​​the main face. The bulk modulus of the material point of the main contact body. The initial volume of the material point of the main contact body. The normal direction of the master surface is used to obtain the normal contact force. Then, the contact friction force was calculated based on the Coulomb friction model. The formula is as follows:

[0104] (5)

[0105] In the formula To contact the cumulative tangential displacement, Let be the Coulomb friction coefficient. Let be the tangential force direction vector, according to We obtain the formula. and These are the velocity vectors from the node and the contact point of the main facet, respectively, where, according to finite element interpolation, the velocity vector of the contact point of the main facet is... , Main face node The velocity vector. Finally, the main facet node. Contact force From the nodal contact force .

[0106] Step 2.4, Update the cumulative tangential displacement of the contact pair. For the actual contact pair, sum the contact tangential displacements from this step, using the following formula:

[0107] (6)

[0108] In the formula For time step.

[0109] The third step involves transferring the physical quantities of the material points of contacting body A to the background mesh nodes, obtaining the mass, momentum, and force of the background mesh nodes. Specifically:

[0110] Step 3.1, Mass Transfer. The mass of the main contact body material point is transferred. The formula for passing this to the background mesh node is as follows:

[0111] (7)

[0112] In the formula Let I be the quality of the background mesh node I at time n. The position of the matter point at time n Weights for background grid node I.

[0113] Step 3.2, Force Transfer. The external forces and stresses at each material point are applied to the background mesh points, and the contact forces between the master facet and slave nodes are transferred to the background mesh points, as shown in the following formula:

[0114] (8)

[0115] In the formula, The nodal force vector of node I in the background mesh; For a point of external force on matter, Let V be the volume of a point mass. Let Cauchy stress tensor be a point mass. The gradient of the interpolation function; and These represent the contact force and position of node k on the main facet, respectively.

[0116] Step 3.3, Velocity Transfer. The velocity of each material point is transferred to the background mesh nodes, using the following formula:

[0117] (9)

[0118] In the formula Let I be the velocity vector of the background mesh node I at time n. Let n be the velocity vector of the material point at time n.

[0119] The fourth step is mesh update. Based on the forces and masses of the background mesh nodes, the node velocities are updated to complete the dynamic calculations at the mesh level.

[0120] Step 4.1, speed update. According to Newton's second law, the force at the mesh nodes... With node quality Calculate the node speed increment and update the node speed using the following formula:

[0121] (10)

[0122] In the formula Let I be the velocity vector of the background mesh node I at time n+1. For time step.

[0123] Step 5: Update the velocity, displacement, stress, and density of the material points based on the velocities of the background mesh nodes. Specifically:

[0124] Step 5.1, update the position of the substance point, using the following formula:

[0125] (11)

[0126] In the formula Let n be the position of the material point at time n+1.

[0127] Step 5.2, update the positions of surface nodes (master patch nodes and slave nodes), using the following formula:

[0128] (12)

[0129] In the formula Let n be the position of the surface node at time n+1.

[0130] Step 5.3, update the velocity of the matter point, the formula is as follows:

[0131] (13)

[0132] In the formula Let be the velocity of the material point at time n+1.

[0133] Step 5.4, update the node speed of surface nodes (master patch nodes and slave nodes), using the following formula:

[0134] (14)

[0135] In the formula Let n be the velocity of the surface node at time n+1.

[0136] Step 5.5, Update Material Point Stress and Density. Update the material point density and update the material point stress according to the material constitutive model. The specific formula is as follows:

[0137] First, calculate the strain increment at the material point. :

[0138] (15)

[0139] (16)

[0140] in, This represents the velocity gradient of the material point at time n+1; This represents the gradient of the interpolation function of the material point p with respect to the background grid node I;

[0141] Then the density of the matter point at time n+1 is calculated. Cauchy stress :

[0142] (17)

[0143] (18)

[0144] The specific calculation method for Cauchy stress at a material point depends on the constitutive model of the material used.

[0145] Step 6: Reset the physical quantities of the background mesh. Set the mass, velocity, and force of all background mesh nodes to zero.

[0146] Step 7: Repeat steps 3, 4, 5, and 6 to update the physical quantities of the material points of contact bodies B, C, and D, as well as the node positions and velocities.

[0147] Step 8: After all contacts have completed the update, update the current time: In the formula Let n be the time at time n+1. If the simulation ends before the timeout, return to step two; otherwise, exit the program and complete the solution.

[0148] After completing the solution, the results of SEMPM (this invention), FEM, and traditional MPM are compared as follows: Figure 3 (a) ~ Figure 3As shown in (i) in the table. The results show that the traditional MPM (comparative example) has obvious gaps at the contact interface due to premature contact; SEMPM and FEM both achieve accurate characterization of the contact interface without gaps or penetration. In terms of computational efficiency, FEM requires a smaller time step to ensure stability due to element distortion, and its CPU time is 4.6 times that of MPM-like methods (see Table 1 for specific time steps and CPU time); the point-to-surface contact algorithm used by SEMPM takes 1.4 times longer than the background mesh contact algorithm used by traditional MPM, but this is negligible compared to the total solution time. Table 2 compares the theoretical background mesh memory consumption of SEMPM and traditional MPM at the initial time step. Where N is the number of grid nodes that need to be allocated memory simultaneously, V is the number of double-precision floating-point numbers required for each node, and D is the number of bytes in a single double-precision floating-point number. The results show that SEMPM, through sequential update and memory reuse strategies, theoretically consumes only about 1 / 3 of the background grid memory compared to traditional MPM.

[0149] Table 1 Comparison of calculation time and time step for multi-body contact examples.

[0150]

[0151] Table 2 Comparison of theoretical memory overhead for multi-body contact examples

[0152]

[0153] This embodiment demonstrates that SEMPM can accurately characterize contact interfaces and avoid non-physical penetration in large deformation multi-body contact scenarios. It can also significantly reduce background mesh memory consumption through memory reuse strategies while maintaining high time efficiency, thus verifying its stability and practicality in multi-body contact problems.

[0154] The above embodiments are merely illustrative of the implementation methods of the present invention, but should not be construed as limiting the scope of the present invention. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of the present invention, and these modifications and improvements all fall within the protection scope of the present invention.

Claims

1. A method for embedding material points on a surface for multi-body contact problems, characterized in that, The surface embedding material point method includes the following steps: The first step is initialization; the physical model of the multi-body contact system is discretized and its parameters are set. The second step is contact force calculation. At the beginning of each time step, real contact pairs are identified through contact search, and the contact force is calculated for real contact pairs using a variable coefficient penalty stiffness mechanism. The second step is specifically as follows: Step 2.1, Contact Search; Perform geometric detection on the main face and slave node to determine if they form a true contact pair; specifically: calculate the minimum distance from the slave node to the main face. If the minimum distance is less than the set contact gap... If the detected slave node-master facet contact pair is not a true contact pair, then the cumulative tangential displacement of the contact pair is determined. Set to zero; Step 2.2, Contact Penetration Calculate the penetration amount from the node relative to the master face for a real contact pair. ,in In the formula For the finite element shape function of node J of the main facet, and Indicates the natural coordinates of the contact points on the main facet. Indicates the spatial location of the main patch node J. The spatial location of the node; Step 2.3, Solving for contact forces; for the actual contact pairs, calculate the main face nodes based on the variable coefficient penalty method. Contact force and contact force from the node ; Step 2.4, update the cumulative tangential displacement of the contact pair; for the actual contact pair, accumulate the contact tangential displacement. The formula is as follows: （6）； In the formula, For time step, and These are the velocity vectors from the node and the contact point of the main surface piece, respectively. Normal to the main facet; The third step is to transfer the physical quantities of the material point of one of the contact bodies to the background mesh node to obtain the mass, momentum and force of the background mesh node. The fourth step is mesh update; based on the forces and masses of the background mesh nodes, the node velocities are updated to complete the dynamic calculations at the mesh level. Fifth step: Update the velocity, displacement, stress, and density of the material points based on the velocities of the background mesh nodes; Step 6: Reset the physical quantities of the background mesh; set the mass, velocity, and force of all background mesh nodes to zero. Step 7: Repeat steps 3, 4, 5, and 6 to update the physical quantities of the material points of other contact bodies, as well as the positions and velocities of the nodes. Allocate the background mesh node memory for the corresponding contact body region and reuse the background mesh node memory allocated in the previous contact body update. Step 8: After all contacts have completed the update, update the current time: In the formula Let n+1 be the time; if If the simulation ends before the timeout, return to step two; otherwise, exit the simulation, complete the solution, and obtain the dynamic simulation results.

2. The surface-embedded material point method for multi-body contact problems according to claim 1, characterized in that, The first step is specifically as follows: Step 1.1, Discretization: Discretize the internal domain of each contact body into material points, define the existence of contact bodies A, B, C, D and other contact bodies in contact; Discretize the boundary of each contact body into massless master facets, define slave nodes matching the master facets as reference points for contact detection; Discretize the computational domain space into a regular structured background mesh, and allocate memory for the background mesh nodes in the region where contact body A is located. Step 1.2, Parameter Settings: Define the material parameters, contact parameters, and calculation parameters for each contact body.

3. A method for embedding material points on a surface for multi-body contact problems according to claim 2, characterized in that, In the first step described above: In step 1.1, the background grid size is taken as twice the average spacing between discrete material particles; In step 1.2, the material parameters include elastic modulus, Poisson's ratio, density, and yield stress, and the contact parameters include contact gap. Main surface area The calculation parameters include the time step. Total simulation time.

4. A method for embedding material points on a surface for multi-body contact problems according to claim 1, characterized in that, In step 2.1, the method for calculating the minimum distance from the node to the main face is as follows: Assuming the point... It is the main face. Upper and lower node positions The points with the shortest distance satisfy a perpendicular relationship: （1）； In the formula, These are the natural coordinates of a point on the main facet. Represents the surface function of the principal patch; Indicates the location of the contact point; For a four-node master face, based on finite element discretization, it can be written in the following form: （2）； In the formula, Indicates the node number of the main facet; Represents the main face node Spatial location; This indicates that for the main face node Finite element shape functions; The location of the contact point can be obtained from formula (2). When the contact point is located When the contact point is outside the main facet, the minimum distance between the main facet and the slave node is at the edge of the main facet, and the following corrections are required: （3）； in, and Indicates the natural coordinates of the contact point on the main facet; Get the location of the contact point Then, calculate the minimum distance from the node to the main face. If the minimum distance satisfies This means that the contact pair is considered a true contact pair; the contact gap mentioned above. Take 1 / 10 of the average side length of the main face.

5. A method for embedding material points on a surface for multi-body contact problems according to claim 4, characterized in that, Step 2.3 specifically involves: First, calculate the normal contact force between the contact pairs. The formula is as follows: （4）； In the formula, The area of ​​the main face; The bulk modulus of the material point of the main contact body; The initial volume of the material point of the main contact body; Normal to the main facet; For contact gap; Secondly, the contact force relative to the normal direction is obtained. Then, the contact friction force was calculated based on the Coulomb friction model. The formula is as follows: （5）； In the formula, For the cumulative tangential displacement of the contact pair; Coulomb friction coefficient; Let be the tangential force direction vector, according to We obtain; in the formula and These are the velocity vectors from the node and the contact point of the main facet, respectively, where, according to finite element interpolation, the velocity vector of the contact point of the main facet is... , Main face node The velocity vector; Finally, the main face node Contact force From the nodal contact force .

6. A method for embedding material points on a surface for multi-body contact problems according to claim 5, characterized in that, The third step is specifically as follows: Step 3.1, mass transfer; transfer the mass of the main contact body material point. The formula for passing this to the background mesh node is as follows: （7）； In the formula, Let I be the quality of the background mesh node I at time n. The position of the matter point at time n Weights for background grid node I; The subscript p indicates the substance point number; Step 3.2, Force Transfer: Apply the external forces and stresses of each material point to the background mesh nodes, and transfer the contact forces between the master facet and slave nodes to the background mesh nodes, as shown in the following formula: （8）； In the formula, The force vector of the mesh node I in the background mesh; External force at a point on matter; Let V be the volume of a point mass. Let Cauchy stress tensor be the material point; The gradient of the interpolation function; and These represent the contact force and position of node k on the main facet, respectively. The position of the matter point at time n Weights for background grid node I; Step 3.3, Velocity Transfer; The velocity of each material point is transferred to the background mesh nodes, using the following formula: （9）； In the formula, Let I be the velocity vector of the background mesh node I at time n; Let n be the velocity vector of the material point at time n; Let p be the mass of the substance point.

7. A method for embedding material points on a surface for multi-body contact problems according to claim 6, characterized in that, The fourth step is specifically as follows: Forces at mesh nodes With node quality Calculate the node speed increment and update the node speed using the following formula: （10）； In the formula, Let I be the velocity vector of background grid node I at time n+1; Let I be the velocity of background grid node I at time n; For time step.

8. A method for embedding material points on a surface for multi-body contact problems according to claim 7, characterized in that, The fifth step is specifically as follows: Step 5.1, update the position of the substance point, using the following formula: （11）； In the formula, The position of the material point at time n+1; Step 5.2, update the position of surface nodes, where surface nodes refer to master patch nodes and slave nodes, as shown in the following formula: （12）； In the formula, The position of the surface node at time n+1; Step 5.3, update the velocity of the matter point, the formula is as follows: （13）； In the formula, Let be the velocity of the material point at time n+1; Step 5.4, update the velocity of the surface nodes, using the following formula: （14）； In the formula, The velocity of the surface node at time n+1; Step 5.5, Update material point stress and density; Update the material point density and update the material point stress according to the material constitutive model, as shown in the following formula: Calculate the strain increment at a material point : （15）； （16）； in, This represents the velocity gradient of the material point at time n+1; This represents the gradient of the interpolation function of the material point p with respect to the background grid node I; Calculate the density of the material point at time n+1. Cauchy stress : （17）； （18）； in, Let p represent the density of the matter point at time n+1; Let p represent the density of the matter point at time n; The spherical strain increment represents the change in material point p. This represents the Cauchy stress at material point p at time n+1; This represents the Cauchy stress at point p of matter at time n; This represents the Cauchy stress increment at material point p.

Images