Surface material point embedding method for multi-body contact problem

Abstract

The invention discloses a surface material point embedding method for a multi-body contact problem, and belongs to the technical field of computational mechanics and numerical simulation. The method comprises the following steps: 1, performing discretization and parameter setting on a physical model of a multi-body contact system; 2, calculating a contact force; 3, transmitting the physical quantity of the material points of the contact body to background grid nodes to obtain node mass, momentum and force; 4, updating the grids; 5, updating the speed, the displacement, the stress and the density of the material point; 6, resetting the physical quantity of the background grid; seventhly, repeating the first step, the fourth step, the fifth step and the sixth step, sequentially updating the physical quantities of the material points of other contact bodies and the positions and the speeds of the slave nodes, and distributing background grid node memories of the areas where the corresponding contact bodies are located; and 8, after all contact bodies are updated, updating the current time to obtain a dynamic simulation result. According to the method, the boundary of the object can be accurately represented, complicated problems such as impact and multi-body contact can be effectively simulated, and meanwhile, the method has good stability and relatively low memory consumption.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of computational mechanics and numerical simulation, and relates to a surface-embedded material point method for multi-body contact problems, in particular to a high-precision numerical simulation method for large deformation multi-body contact problems based on a material point method. BACKGROUND

[0002] In the engineering field, multi-body contact phenomena (such as collision, impact, fluid-structure coupling, etc.) are widely present, and numerical simulation thereof is crucial for engineering design and analysis. In traditional numerical methods, the finite element method (FEM) can accurately describe the contact interface through the grid boundary, but in the extreme deformation scenario, grid distortion easily occurs, leading to unstable internal force updating, reduced time step, and serious impact on calculation efficiency and accuracy.

[0003] The material point method (MPM), as a typical meshless method, combines the advantages of Lagrangian description and Eulerian description, and performs well in large deformation problems (such as metal plastic forming, geotechnical dynamic response, and super-high-speed collision), but traditional MPM has significant defects in contact problems: first, the interface representation is ambiguous, and contact detection relies on background grids, which easily leads to early contact or non-physical penetration due to the mismatch between physical boundaries and grid boundaries; second, the contact normal vector is approximated by the grid mass gradient, which has poor accuracy when the grid is rough or the material points are unevenly distributed; third, when there are multiple bodies in contact, background grid memory needs to be allocated for all contact bodies, leading to large memory consumption; fourth, the contact calculation based on background grids requires the background grid size and node position of each contact body to be consistent, which cannot set the appropriate background grid size for each contact body.

[0004] In existing improved methods, although the background grid contact algorithm based on physical distance auxiliary detection can alleviate the problem of premature contact, it still cannot solve the non-physical penetration of the contact interface; the implicit interface method based on level set is complex to implement in a three-dimensional multi-body system and has high computational cost; in addition, optimization schemes for contact detection logic also have limitations: for example, a kind of MPM contact algorithm disclosed in Chinese invention patent (publication number CN115410663A) directly calculates the contact force according to the position of the material point, which is equivalent to migrating the "pinball method" to the MPM framework, and realizing contact judgment by setting the virtual collision ball radius between the material points, effectively avoiding the false contact problem caused by the dependence of the traditional grid, but the method has low accuracy in representing complex topological contact interfaces (such as components with grooves and sharp corners), and the setting of the virtual ball radius easily leads to uneven distribution of contact force; another typical improvement scheme, such as the application of the NTS (Node-To-Segment) point-face contact algorithm in MPM proposed in Chinese invention patent (publication number CN114357717A), optimizes the contact detection accuracy by establishing the mapping relationship between the material points and the surface segments of the contact body, but it uses a fixed penalty coefficient to solve the contact constraint, which is difficult to match the dynamic contact stiffness requirement in extreme working conditions such as high-speed impact and high strain rate, and is easy to cause non-physical penetration of the contact interface.

[0005] In summary, the existing MPM contact framework is limited by the accuracy of complex interfaces, or is difficult to adapt to extreme conditions, or has a contradiction between implementation complexity and computational cost, and there is an urgent need for a multi-body contact numerical method that can accurately represent boundaries, efficiently handle contact constraints, and has low memory consumption. SUMMARY

[0006] In view of the problems of interface representation ambiguity, inaccurate contact detection, easy non-physical penetration, grid distortion affecting stability, and large memory consumption in traditional MPM and existing improved methods in multi-body contact problems, the present application provides a surface-embedded material point method (SEMPM) for multi-body contact problems, which realizes high-precision and high-stability numerical simulation of multi-body contact scenes.

[0007] To achieve the above purpose, the present application adopts the following technical solutions:

[0008] A surface-embedded material point method for multi-body contact problems, comprising the following steps:

[0009] Step 1: initialization. Discretize and set parameters for the physical model of the multi-body contact system to provide basic data for subsequent contact calculation and object update. Specifically as follows:

[0010] Step 1.1, Discrete treatment. The inner domain of each contactor is discretized into material points, which are used to carry physical information such as mass, stress, etc. Here, it is defined that contactors A, B, C, D and other contactors are in contact. The boundary of each contactor is discretized into massless master surface patches, which are only used for geometric boundary representation and do not carry physical quantities, and slave nodes matching the master surface patches are defined as reference points for contact detection. The calculation domain space is discretized into a regular structured background mesh, and the background mesh nodes in the area where contactor A is located are allocated memory. The size of the background mesh is taken as twice the average gap between the discretized material points.

[0011] Step 1.2, Parameter setting. Define the material parameters, contact parameters and calculation parameters of each contactor. The material parameters include elastic modulus, Poisson's ratio, density, yield stress, etc. The contact parameters include contact gap , master surface patch area , etc. The calculation parameters include time step , total simulation time, etc.

[0012] Second step, contact force calculation. At the beginning of each time step, identify the real contact pair through contact search, and calculate the contact force using the variable coefficient penalty stiffness mechanism for the real contact pair to ensure the reliability of the contact constraint. Specifically:

[0013] Step 2.1, contact search. Geometric detection is performed on the master surface patch and the slave node to determine whether it constitutes a real contact pair. Specifically: calculate the minimum distance from the slave node to the master surface patch. If the minimum distance is less than the set contact gap , then the slave node-master surface patch is determined to be a real contact pair. If the detected slave node-master surface patch contact pair is not a real contact pair, then the tangential cumulative displacement of the contact pair is set to zero.

[0014] Further, the minimum distance from the slave node to the master surface patch is calculated as follows: assume that point is the point on the master surface patch that is closest to the slave node position , at which time the perpendicular relationship should be satisfied:

[0015] (1)

[0016] In the formula, is the natural coordinate of a point on the master surface patch; represents the master surface patch function; represents the contact point position.

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

[0018] (2)

[0019] where, represents the main surface patch node number; represents the spatial position of the main surface patch node ; represents the finite element shape function for the main surface patch node .

[0020] According to formula (2), the contact point position can be solved. When the contact point position , it means that the contact point is outside the main surface patch, and at this time, the minimum distance between the main surface patch and the slave node is theoretically on the edge of the main surface patch, that is, the following correction needs to be made:

[0021] (3)

[0022] where, and represent the natural coordinates of the contact point on the main surface patch;

[0023] After obtaining the contact point position , the minimum distance from the slave node to the main surface patch is calculated, and if the minimum distance satisfies , it is considered that the contact pair is a real contact pair. The contact gap is taken as 1 / 10 of the average side length of the main surface patch.

[0024] Step 2.2, contact pair penetration amount calculation. For a real contact pair, the penetration amount of the slave node relative to the main surface patch is calculated, where ;

[0025] Step 2.3, contact force solving. For a real contact pair, the contact force of the main surface patch node and the contact force of the slave node are calculated based on the variable coefficient penalty method. Specifically:

[0026] First, the normal contact force of the contact pair is calculated, and the formula is as follows:

[0027] (4)

[0028] where, is the area of the main surface patch; is the bulk modulus of the material point of the main contact body; is the initial volume of the material point of the main contact body; is the normal of the main surface patch; is the contact gap.

[0029] Second, the normal contact force Then, the contact pair friction force is calculated according to the Coulomb friction model , which is as follows:

[0030] (5)

[0031] wherein, is the cumulative tangential displacement of the contact pair; is the Coulomb friction coefficient; is the tangential force direction vector, which is obtained according to . In the formula, and are the velocity vectors of the slave node and the master patch contact point, respectively, wherein the velocity vector of the master patch contact point , is obtained according to the finite element interpolation.

[0032] Finally, the contact force of the master patch node is obtained from the slave node contact force .

[0033] Step 2.4, contact pair cumulative tangential displacement update. For the real contact pair, the contact tangential displacement of this step is accumulated as , which is as follows:

[0034] (6)

[0035] wherein, is the time step.

[0036] Thirdly, the physical quantities of the material points of a contact body are transferred to the background grid nodes to obtain the mass, momentum and force of the background grid nodes. Specifically:

[0037] Step 3.1, mass transfer. The mass of the master contact body material point is transferred to the background grid node, which is as follows:

[0038] (7)

[0039] wherein, is the mass of the background grid node I at time n, is the weight of the material point position on the background grid node I at time n; the subscript p represents the material point serial number.

[0040] Step 3.2, force transfer. The external force and stress of each material point are applied to the background grid node, and the master patch and slave node contact forces are transferred to the background grid node, which is as follows:

[0041] ​ (8)

[0042] where, is the grid node force vector of background grid node I; is the external force of material point; is the volume of material point; is the Cauchy stress tensor of material point; is the gradient of interpolation function; and are the contact force and position of master facet node k, respectively; is the position of material point at time n is the weight of background grid node I.

[0043] Step 3.3, velocity transfer. The velocity of each material point is transferred to the background grid node, the formula is as follows:

[0044] (9)

[0045] where, is the velocity vector of background grid node I at time n; is the velocity vector of material point at time n; is the mass of material point p.

[0046] Fourth step, grid update. Based on the force and mass of the background grid node, the node velocity is updated, and the dynamics calculation of the grid level is completed. Specifically:

[0047] Step 4.1, velocity update. According to Newton's second law, the node velocity increment is calculated by the grid node force and the node mass , and the node velocity is updated, the formula is as follows:

[0048] (10)

[0049] where, is the velocity vector of background grid node I at time n+1; is the velocity of background grid node I at time n; is the time step.

[0050] Fifth step, according to the velocity of the background grid node, the velocity, displacement, stress and density of the material point are updated. Specifically:

[0051] Step 5.1, material point position update, the formula is as follows:

[0052] (11)

[0053] where, is the position of material point at time n+1.

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

[0055] (12)

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

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

[0058] (13)

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

[0060] Step 5.4, velocity update of surface nodes, the formula is as follows:

[0061] (14)

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

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

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

[0065] (15)

[0066] (16)

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

[0068] Then calculate the density of the material point at time n+1 , Cauchy stress :

[0069] (17)

[0070] (18)

[0071] Where, represents the density of the material point p at time n+1; represents the density of the material point p at time n; represents the spherical strain increment of the material point p; represents the Cauchy stress of the material point p at time n+1; represents the Cauchy stress of the material point p at time n; represents the Cauchy stress increment of the material point p;

[0072] The specific calculation method of the material point Cauchy stress depends on the material constitutive adopted.

[0073] Step 6, reset the background grid physical quantity. Set all background grid node mass, velocity, force to zero.

[0074] Step 7, repeat step 3, step 4, step 5, step 6, update the material point physical quantity of other contact bodies and the slave node position and velocity in turn, allocate the background grid node memory of the region where the corresponding contact body is located, and reuse the background grid node memory already allocated in the previous contact body update.

[0075] Step 8, after all contact bodies are updated, update the current time: , wherein is the time at time n+1. If is less than the end time of the simulation, return to step 2, otherwise exit, complete the solution, and obtain the dynamic simulation result.

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

[0077] (1) The contact boundary representation in the calculation process of the present application is accurate (corresponding to step 1): the object boundary is discretized by massless master surface patches to form a Lagrangian surface grid that fits the object, solving the problem of mismatch between traditional MPM physical boundaries and grid boundaries;

[0078] (2) The contact constraint in the calculation process of the present application is reliable (corresponding to step 2): 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 requirement under extreme working conditions;

[0079] (3) The memory consumption in the calculation process of the present application is low (corresponding to step 7): the sequential update and grid reuse strategy is adopted, and the background grid only needs to allocate "the node number required by the largest contact body", instead of the "sum of all contact body node numbers" in traditional MPM;

[0080] (4) The grid consistency requirement in the calculation process of the present application is low (corresponding to steps 3-7): the contact calculation is based on the geometric relationship between the master surface patch and the slave node, which is decoupled from the background grid, allowing each contact body to adopt a grid size suitable for itself. BRIEF DESCRIPTION OF DRAWINGS

[0081] Figure 1is a flow chart of the algorithm of the present application.

[0082] Figure 2 is a schematic diagram of a geometry model of a multi-body large deformation contact example.

[0083] Figure 3 is a comparison diagram of contact processes of various numerical methods of a multi-body large deformation contact example. Figure 3 (a) in is a structure deformation diagram of the present application at 0.5 seconds; Figure 3 (b) in is a structure deformation diagram of the present application at 0.75 seconds; Figure 3 (c) in is a structure deformation diagram of the present application at 1.0 seconds; Figure 3 (d) in is a structure deformation diagram of the finite element at 0.5 seconds; Figure 3 (e) in is a structure deformation diagram of the finite element at 0.75 seconds; Figure 3 (f) in is a structure deformation diagram of the finite element at 1.0 seconds; Figure 3 (g) in is a structure deformation diagram of the traditional material point method at 0.5 seconds; Figure 3 (h) in is a structure deformation diagram of the traditional material point method at 0.75 seconds; Figure 3 (i) in is a structure deformation diagram of the traditional material point method at 1.0 seconds. DETAILED DESCRIPTION

[0084] The present application will be further described below in combination with a specific implementation case.

[0085] The following tests a rigid body extrusion of multiple elastic-plastic bodies, wherein the specific geometry model and example setting are as follows: Figure 2 In the present embodiment, among the contact bodies A, B, C and D, the contact body A is a rigid disc, and the contact bodies B, C and D are elastic-plastic bodies; the mass of the rigid disc is 30000kg, the diameter is 4m, and the rigid disc freely falls under the action of gravity g=9.81m / s 2 , the initial distance between B, C and D is 0.1m, and D is placed on a smooth plane. The height of B, C and D is 2m, and the width is 6m, 8m and 10m respectively. The elastic-plastic body adopts a linear hardening constitutive model, the Young's modulus is 10Mpa, the Poisson's ratio is 0.3, the density , the yield stress is 0.157MPa, and the strain hardening coefficient is 0.425MPa.

[0086] A surface-embedded material point method for multi-body contact problems, as shown in Figure 1 , comprises the following steps:

[0087] Step 1, initialization. Discretize the physical model of the multi-body contact system and set the parameters, to provide basic data for subsequent contact calculation and object update.

[0088] Step 1.1, Discrete treatment. Discretize each contactor interior into 13326 material points with average material point spacing of 0.1m and background mesh spacing of 0.2m; discretize the contactor geometric surface into 626 quadrilateral unstructured meshes as contact surface (define all surface meshes as master patch and all surface mesh nodes as slave node).

[0089] Step 1.2, Parameter setting. Set the material parameters of each contactor material point, contact gap =0.01m, time step , total simulation time 1s according to the material of the example.

[0090] Second step, contact force calculation (input master patch node and slave node position, get master patch node and slave node contact force). At the beginning of each time step, identify the real contact pair through contact search, and calculate the contact force based on the variable coefficient penalty stiffness mechanism for the real contact pair to ensure the reliability of the contact constraint.

[0091] Step 2.1, contact search. Geometric detection is performed on the master patch and slave node to determine whether a real contact pair is formed; specifically: calculate the minimum distance from the slave node to the master patch, if the minimum distance is less than the set contact gap , then determine that the "slave node-master patch" is a real contact pair;

[0092] In this embodiment, the minimum distance from the slave node to the master patch is calculated as follows: assuming that the point is the point on the master patch closest to the slave node position , at this time the perpendicular relationship should be satisfied:

[0093] (1)

[0094] In the formula, is the natural coordinate of a point on the master patch. The master patch is discretized based on the finite element, which can be written as the following form:

[0095] (2)

[0096] In the formula, is the spatial position of the master patch node . is the finite element shape function for the master patch node .

[0097] According to formula (2), the contact point position can be solved. When the natural coordinate of the contact point is solved, it means that the contact point is outside the patch, at this time the minimum distance between the master patch and the slave node is theoretically on the edge of the master patch, that is, the following correction needs to be made:

[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 surface. 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 Contact force from node .

[0106] Step 2.4, contact pair cumulative tangential displacement update. For real contact pair, accumulate this step contact tangential displacement, formula as follows:

[0107] (6)

[0108] In the formula is the time step.

[0109] Third, the contact body A of the material point physical quantity is passed to the background grid node, and the mass, momentum and force of the background grid node are obtained. Specifically:

[0110] Step 3.1, mass transfer. The mass of the main contact body material point is transferred to the background grid node, and the formula is as follows:

[0111] (7)

[0112] In the formula is the mass of the background grid node I at time n, is the position of the material point at time n The weight of the background grid node I.

[0113] Step 3.2, force transfer. The external force and stress of each material point are applied to the background grid node, and the main face sheet and slave node contact force are transferred to the background grid node, and the formula is as follows:

[0114] (8)

[0115] In the formula, is the node force vector of the background grid node I; is the external force of the material point, is the volume of the material point, is the Cauchy stress tensor of the material point, is the interpolation function gradient; and are the contact force and position of the main face sheet node k respectively.

[0116] Step 3.3, velocity transfer. The velocity of each material point is transferred to the background grid node, and the formula is as follows:

[0117] (9)

[0118] In the formula is the velocity vector of the background grid node I at time n, is the velocity vector of the material point at time n.

[0119] Step 4, Mesh update. Based on the background mesh node force and mass, update the node velocity, complete the dynamics calculation at the mesh level.

[0120] Step 4.1, Velocity update. According to Newton's second law, the node velocity increment is calculated by the mesh node force and the node mass The node velocity is updated, and the formula is as follows:

[0121] (10)

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

[0123] Step 5, according to the background mesh node velocity, update the velocity, displacement, stress and density of the material point. Specifically:

[0124] Step 5.1, material point position update, the formula is as follows:

[0125] (11)

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

[0127] Step 5.2, surface node (main and slave nodes) position update, the formula is as follows:

[0128] (12)

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

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

[0131] (13)

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

[0133] Step 5.4, surface node (main and slave nodes) node velocity update, the formula is as follows:

[0134] (14)

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

[0136] Step 5.5, material point stress, density update. Update material point density, and update material point stress according to material constitutive, the specific formula as follows:

[0137] First, calculate the material point strain increment

[0138] (15)

[0139] (16)

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

[0141] Then, calculate the density of the material point at time n+1 , Cauchy stress :

[0142] (17)

[0143] (18)

[0144] wherein, the specific calculation method of the material point Cauchy stress depends on the adopted material constitutive.

[0145] Step 6, reset the background grid physical quantity. Set all background grid node mass, velocity, force to zero.

[0146] Step 7, repeat step 3, step 4, step 5, step 6, update the material point physical quantity of the contact body B, C, D and the node position and velocity in turn.

[0147] Step 8, after all the contact bodies are updated, update the current time: , wherein is the time at time n+1. If is less than the end time of simulation, return to step 2, otherwise exit the program, complete the solution.

[0148] After the solution is completed, the results of SEMPM (the present application), FEM and traditional MPM are compared as Figure 3 (a) in Figure 3 ​The results show that the traditional MPM (comparative example) has obvious gaps at the contact interface due to the premature contact problem; both the SEMPM and the FEM achieve accurate characterization of the contact interface without gaps or penetration. In terms of computational efficiency, the FEM needs to use a smaller time step to ensure stability due to element distortion, and the CPU time is 4.6 times that of the MPM method (see Table 1 for specific time steps and CPU time). The point-surface contact algorithm used by the SEMPM takes 1.4 times the time of the background grid contact algorithm used by the traditional MPM, but this is negligible compared to the total solving time. Table 2 compares the theoretical background grid 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 at the same time, V is the number of double-precision floating-point numbers required for each node, and D is the number of double-precision floating-point numbers for a single node). The results show that SEMPM uses sequential updating and memory reuse strategies, and the theoretical background grid memory consumption is only about 1 / 3 of that of the traditional MPM.

[0149] Table 1 Comparison of calculation time and time step of multi-body contact example

[0150]

[0151] Table 2 Comparison of theoretical memory consumption of multi-body contact example

[0152]

[0153] The verification of this embodiment shows that SEMPM can accurately characterize the contact interface and avoid non-physical penetration in the large deformation multi-body contact scenario, and can significantly reduce the background grid memory consumption through the memory reuse strategy while maintaining high time efficiency, verifying its stability and practicality in multi-body contact problems.

[0154] The above-described embodiments only express the implementation of the present application, but should not be construed as limiting the scope of the present application. It should be noted that for those skilled in the art, without departing from the concept of the present application, a number of modifications and improvements can be made, which are all within the scope of the present application.

Claims

1. A surface-embedded material point method for multi-body contact problems, characterized in that, The surface-embedded material point method comprises the following steps: The first step is initialization; the physical model of the multi-body contact system is discretized and parameters are set; The second step is contact force calculation; at the beginning of each time step, the real contact pairs are identified through contact search, and the contact force is calculated for the real contact pairs using a variable coefficient penalty stiffness mechanism; The third step is to transfer the material point physical quantities of one of the contact bodies to the background grid nodes to obtain the mass, momentum and force of the background grid nodes; The fourth step is grid update; based on the force and mass of the background grid nodes, the node velocity is updated to complete the dynamics calculation at the grid level; The fifth step is to update the velocity, displacement, stress and density of the material points according to the velocity of the background grid nodes; The sixth step is to reset the background grid physical quantities; the mass, velocity and force of all background grid nodes are set to zero; The seventh step is to repeat the third step, the fourth step, the fifth step and the sixth step to update the material point physical quantities and the node position and velocity of the other contact bodies in turn, allocate the background grid node memory of the region where the corresponding contact body is located, and reuse the background grid node memory allocated in the previous contact body update; Step 8, after all contacts are updated, update the current time: , where is the time at n+1; if is less than the end time of the simulation, go back to step 2, otherwise exit, and the solution is obtained, and the dynamic simulation result is obtained.

2. The surface-embedded material point method for multi-body contact problems of claim 1, wherein, The first step specifically comprises: Step 1.1, discretization; the internal domain of each contact body is discretized into material points, and the contact between contact bodies A, B, C, D and other contact bodies is defined; the boundary of each contact body is discretized into massless main face pieces, and slave nodes matched with the main face pieces are defined as reference points for contact detection; the calculation domain space is discretized into a regular structured background grid, and the background grid node memory of the region where the contact body A is located is allocated; Step 1.2, parameter setting; the material parameters, contact parameters and calculation parameters of each contact body are defined.

3. The surface-embedded SPH method for multi-body contact problems of claim 2, wherein, In the first step: In step 1.1, the size of the background grid is taken as twice the average gap between the material points after discretization; In step 1.2, the material parameters include elastic modulus, Poisson's ratio, density, yield stress, and the contact parameters include contact gap , the area of the main surface sheet , the calculation parameters include time step , total simulation time.

4. The surface-embedded material point method for multi-body contact problems of claim 3, wherein, The second step specifically comprises: Step 2.1, contact search; geometric detection of the master surface patch and the slave node to determine whether a real contact pair is formed; specifically: calculate the minimum distance from the slave node to the master surface patch, if the minimum distance is less than the set contact gap , then determine that the slave node-master surface patch is a real contact pair; if the detected slave node-master surface patch contact pair is not a real contact pair, then the tangential cumulative displacement of the contact pair is zeroed. Step 2.2, contact pair penetration amount Compute; for real contact pairs, compute penetration amount from node to master face where ; Step 2.3, contact force solution; For real contact pairs, the contact force at the master surface patch node is calculated based on a variable coefficient penalty method and the slave node contact force ;​ Step 2.4, contact pair cumulative tangential displacement update; for real contact pairs, accumulate contact tangential displacement The formula is as follows: （6） In the formulae, is the time step.

5. The surface-embedded material point method for multi-body contact problems of claim 4, wherein, The minimum distance from the slave node to the master patch in step 2.1 is calculated as follows: assume the point is the point on the master patch with the minimum distance to the slave node position satisfying the perpendicular relationship: （1） wherein is the natural coordinate of a certain point on the main face sheet; denotes the main face sheet surface function; denotes the contact point position; For a four-node main face piece, based on finite element discretization, it can be written as follows: （2） wherein denotes the main face sheet node number; denotes the spatial position of the main face sheet node ; denotes the finite element shape function for the main face sheet node ; According to equation (2), the contact point position can be solved When the contact point position When the contact point is outside the main face sheet, the minimum distance between the main face sheet and the slave node is at the edge of the main face sheet, and the following correction needs to be made: （3） wherein and denotes the natural coordinates of the contact point on the master surface patch. Obtaining the contact point position Then, the minimum distance from the node to the master patch is calculated If the minimum distance meets The contact pair is considered as a real contact pair; the contact gap The contact gap is 1 / 10 of the average edge length of the master patch.

6. The surface-embedded material point method for multi-body contact problems of claim 5, wherein, Step 3.2, force transmission; the external force and stress of each material point are applied to the background grid nodes, and the contact force of the main face piece and the slave node is transmitted to the background grid node, and the formula is as follows: First, the normal contact force of the contact pair is calculated and is given by the formula （4） wherein is the area of the master facet; is the bulk modulus of the master contact; is the initial volume of the master contact; is the normal to the master facet; is the contact gap; Second, the normal contact force of the contact pair is obtained Then, the friction force of the contact pair is calculated according to the Coulomb friction model The formula is as follows: （5） wherein is the contact pair cumulative tangential displacement; is the Coulomb friction coefficient; is the tangential force direction vector, obtained from ; wherein and are the slave node and master patch contact point velocity vectors, respectively, wherein the master patch contact point velocity vector , is the velocity vector of the master patch node ; Finally, the contact force of the main face sheet node from the node contact force .​ 7. A surface-embedded material point method for multi-body contact problems according to claim 6, wherein, Step 3.3, velocity transmission; the velocity of each material point is transmitted to the background grid node, and the formula is as follows: Step 3.1, Mass Transfer; the mass of the primary contactor material point is transferred to the background grid nodes, as follows: Step 3.1, Mass Transfer; the mass of the primary contactor material point is transferred to the background grid nodes, as follows: （7） wherein mass of background grid node I at time n, position of particle at time n weight of background grid node I; The fourth step specifically comprises: The fifth step specifically comprises: （8） where, is the grid node force vector for background grid node I; is the external force on the material point; is the volume of the material point; is the Cauchy stress tensor of the material point; is the gradient of the interpolation function; and are the contact force and position of the master facet node k, respectively; is the material point position at time n is the weight of background grid node I. Step 5.1, material point position update, the formula is as follows: （9） wherein is the background grid node I velocity vector at time n; is the material point n velocity vector; is the mass of material point p.

8. The surface-embedded material point method for multi-body contact problems of claim 7, wherein, Step 5.2, position update of surface nodes, wherein the surface nodes refer to main face piece nodes and slave nodes, and the formula is as follows: By grid node force With node mass Calculate node velocity increment, update node velocity, formula as follows: （10） wherein is the velocity vector of the background grid node I at time step n+1; is the velocity of the background grid node I at time step n; is the time step.

9. The surface-embedded material point method for multi-body contact problems of claim 8, wherein, Step 5.3, material point velocity update, the formula is as follows: Step 5.4, velocity update of surface nodes, the formula is as follows: （11） wherein is the position of the material point at time n+1. Step 5.5, material point stress and density update; the material point density is updated, and the material point stress is updated according to the material constitutive relation, and the specific formula is as follows: （12） In the formula, is the position of the surface node at time n+1. ​ （13） wherein is the velocity of the material point at time n+1. ​ （14） wherein Vn+1is the velocity of the surface node at time n+1; ​ Computing material point strain increment : （15） （16） wherein, denotes the velocity gradient of the material point n+1 at time t; denotes the gradient of the interpolation function of the material point p to the background grid node I; calculating the density of the material point at the n+1 time , cauchy stress : （17） （18） wherein, represents the material point p density at time n + 1 ; represents the material point p density at time n ; represents the spherical strain increment of the material point p; represents the Cauchy stress of the material point p at time n + 1 ; represents the Cauchy stress of the material point p at time n ; represents the Cauchy stress increment of the material point p.

Images