A cell-based peridynamic contact modeling method

Abstract

The present invention relates to the technical fields of numerical modeling and material damage, and provides a cell-based peridynamic contact modeling method. The core lies in the precise definition of contact elements, surfaces and judgment criteria, as well as the calculation of contact forces. By constructing contact cells and introducing a distance influence function, this method effectively solves the problems of abnormal tangential forces and unstable contact forces. Through numerical simulations of impact, friction and rolling contact, the accuracy of this method is verified by comparison with methods such as finite element analysis. Research has confirmed that this method can reliably simulate complex contact problems involving material damage.

Description

Technical Field

[0001] The present invention relates to the technical fields of numerical modeling and material damage, and more particularly, to a cell-based peridynamic contact modeling method. Background Art

[0002] Contact damage is a common engineering problem that can lead to the fracture and failure of engineering structures, such as rolling contact fatigue, impact fracture, etc. The finite element method based on classical continuum mechanics is a commonly used tool for contact modeling. However, classical continuum mechanics uses partial differential equations to describe mechanical problems and faces the problem of the failure of mathematical tools in discontinuous regions such as damage. Therefore, it is difficult for the finite element method to simulate contact problems involving material damage and fracture.

[0003] Nonlocal peridynamics is an emerging mechanical theory that uses integral equations instead of partial differential equations to describe mechanical problems, and there is no problem of the failure of mathematical tools at discontinuous points such as material damage. Peridynamics directly characterizes material damage through bond fracture and can spontaneously simulate damage propagation without additional criteria. Therefore, the peridynamics theory has great advantages in the study of damage problems.

[0004] Compared with the classical finite element analysis method, peridynamic contact modeling has not received wide attention and is still an open research field. The existing domestic and foreign peridynamic contact models generally have the following problems in the application process: abnormal tangential contact force, unstable contact force, and these models are usually only applicable to a single impact or sliding friction contact scenario. Therefore, there is an urgent need to develop a new peridynamic contact algorithm for studying contact problems involving material damage and fracture. Summary of the Invention

[0005] The content of the present invention is to provide a cell-based peridynamic contact modeling method for contact damage research, which overcomes the defects of false tangential force, unstable contact force, and only being applicable to a single contact scenario in the existing peridynamic contact modeling method.

[0006] A cell-based peridynamic contact modeling method according to the present invention includes the following steps:

[0007] S1. Set and initialize parameters and variables such as materials, geometries, simulations, etc.;

[0008] S2. Divide the two objects participating in the contact action into a slave object and a master object;

[0009] S3. Adopt a meshless method to discretize the slave object and the master object into peridynamic particles, define the particles that may participate in the contact interaction as contact particles, and call the contact particles on the slave object and the master object slave particles and master particles, and construct a set R of slave particles and master particless and R m ;

[0010] S4. For each contact particle on the main object, i.e., the main particle y i Construct a contact unit

[0011] S5. Set the main particle y i and the secondary particle y j The total normal contact force received and the total tangential contact force 、 to 0;

[0012] S6. Access all secondary particles y in the set j , and traverse each secondary particle y in the secondary object j ;

[0013] S7. Each secondary particle y j Access all main particles y in the set i , and traverse each main particle y in the main object i ;

[0014] S8. Determine the contact plane of the contact unit constructed by the main particle y i ;

[0015] S9. Judge whether the condition for contact occurrence is satisfied. If satisfied, continue with the subsequent operations; otherwise, terminate the contact calculation between particles y i and y j ;

[0016] S10. After contact occurs, calculate the normal overlap depth of the material Contact stiffness and the distance influence function ω ij ;

[0017] S11. Calculate the normal contact force between the secondary particle y j and the main particle y i

[0018] S12. Add the calculated normal contact force between the secondary particle y j and the main particle y i to to obtain the total normal contact force of y i 、y j ;

[0019] S13. If y j and y i are in contact for the first time, record the secondary particle y at the current moment​​​j With the main particle y i , for subsequent analysis; otherwise, recalculate and the normal contact force and accumulate it into the total normal contact force to update the contact state;

[0020] S14. Define the trial tangential contact force;

[0021] S15. Calculate the relative slip distance

[0022] S16. Define the tangential contact boundary according to the Coulomb friction criterion

[0023] S17. Determine the tangential contact force according to the tangential contact boundary;

[0024] S18. Sum up all the tangential contact forces generated by the slave particles and the main particle to obtain the total tangential contact force acting on y j 、y i ;

[0025] S19. Apply the calculated total tangential contact force to the contacting particles in the form of body force density, and update the particle displacement by the numerical scheme of PD;

[0026] S20. If the termination condition is satisfied, terminate the calculation; otherwise, return to S6.

[0027] Preferably, in S4, each contact unit includes the current contacting particle and its two adjacent contacting particles, that is, any contacting particle y i on the main object, and the contact unit is defined as the set of y i and its two adjacent contacting particles y i1 and y i2 :

[0028]

[0029] Based on the spatial position relationship, the contacting particles on the main object are divided into three categories: non-boundary particles, boundary particles, and particles adjacent to defects; for non-boundary particles, the adjacent particles of the contact unit are the two closest contacting particles to it; for the particles y b at the left and right boundaries, the adjacent particles y b2 inside the main object are determined by searching for the closest contacting particles, but there are no adjacent particles outside the main object. Therefore, a virtual adjacent particle y bv is introduced, and its position is determined by linear extrapolation from the internal adjacent particle:

[0030] y bv = 2y b - yb2

[0031] For the contact particle y close to the defect d , it is y d Introduce a virtual adjacent particle y dv , y dv The spatial position of is determined by linear extrapolation through another adjacent particle;

[0032] All contact particles should have at least one adjacent particle that does not require additional processing, that is, there is at least one adjacent particle at a distance of Δx from y i to avoid the problem of not being able to find an adjacent particle. i

[0033] Preferably, in S8, two adjacent particles in the contact unit are used to determine the contact surface direction, and the contact surface position is determined by three particles together;

[0034] First, the contact surface direction is defined by two adjacent particles, where the direction of the vector y i2 -y i1 is the tangential direction of the contact surface, and the main particle y i constructs a plane χ i2 -y i1 parallel to, which is the contact plane defined by the contact unit i , and the projection points of the adjacent particles y and y i1 and y i2 on χ i are the boundaries of the plane χ and ; i

[0035] The unit tangential vector of the contact surface is written as:

[0036]

[0037] is determined by the following formula and to determine the spatial positions of and:

[0038] k = 1 or 2

[0039] where k represents any one of the adjacent particles;

[0040] χ i is the contact plane defined by the contact unit and is expressed as:

[0041]

[0042] Define​​ is χ i of the unit normal vector, pointing outside the contact surface, the unit normal vector and the unit tangential vector are perpendicular, that is, there is the following relationship:

[0043]

[0044] From particle y j The projection point on χ i is Its position is determined by the following formula:

[0045]

[0046] Preferably, in S9, from the contact particle y j The projection point on χ i is When the projection point of the contact particle y j with respect to the contact surface χ i is inside χ i and the projection distance is less than the critical value, it is determined that the contact event occurs:

[0047]

[0048] Δx i represents the particle spacing of the main particle, and Δx j represents the particle spacing of the secondary particle; is the projection distance of the projection inside χ i ; is y j and χ i the critical contact distance between.

[0049] Preferably, in S10, the material normal overlap depth the contact stiffness and the distance influence function ω ij The calculation formulas are as follows:

[0050]

[0051] In the formula, p si is the control coefficient of the contact stiffness, K ij represents the bulk modulus, δ ij is the neighborhood radius, ΔV i and ΔV j are the bulk moduli of the material.

[0052] Preferably, in S11, the normal contact force The calculation formula is as follows:

[0053]

[0054] f cn (y i , y j ) represents the normal contact force between y j and y i .

[0055] Preferably, in S12, the total normal contact force of y i is as follows:

[0056]

[0057] f cn (y i , t) represents the total normal force on the main particle y i ;

[0058] According to Newton's third law, the contact force received by the particle y j is equal in magnitude and opposite in direction to that of y i , so the expression is:

[0059]

[0060] is the contact force received by the secondary particle, is the contact force received by the main particle;

[0061] Sum the reaction forces of the secondary particle y j and all the main particles y i to obtain the total normal contact force acting on y j as:

[0062]

[0063] where R m is the set of main particles.

[0064] Preferably, in S14, the trial tangential contact force is as follows:

[0065]

[0066] is the trial tangential contact force.

[0067] Preferably, in S15, the relative sliding distance is the difference between the tangential distances of two contacting particles from the initial contact time t c to the current time t, and its expression is as follows:

[0068]

[0069] In S16, the tangential contact boundary has the following expression:

[0070]

[0071] μ represents the friction coefficient;

[0072] In S17, the expression of the tangential contact force is as follows:

[0073]

[0074] Preferably, in S18, the expression of the total tangential contact force acting on y i is as follows:

[0075]

[0076] Similarly, according to Newton's third law, the expression of the total tangential contact force acting on y j is as follows:

[0077]

[0078] The present invention provides a cell-based peridynamic contact modeling method for contact damage research, which can solve the defects of false tangential forces, unstable contact forces, and only being applicable to a single contact scenario existing in the existing peridynamic contact modeling methods; the present invention is very suitable for contact damage problems and research on various contact scenarios. BRIEF DESCRIPTION OF THE DRAWINGS

[0079] Figure 1 is a flowchart of a cell-based peridynamic contact modeling method;

[0080] Figure 2 is a schematic diagram of contact area discretization and contact elements;

[0081] Figure 3 is a schematic diagram of the Kalthoff-Winkler experiment;

[0082] Figure 4(a) is a schematic diagram of the comparison of contact forces varying with time calculated by the EBCM and EMU program methods;

[0083] Figure 4(b) is a schematic diagram of the comparison of projectile velocities varying with time calculated by the EBCM and EMU program methods;

[0084] Figure 5 is a schematic diagram of the damage distribution at the termination moment calculated by the EBCM method and the EMU program;

[0085] Figure 6Schematic diagram of the frictional contact test for two flexible flat plates;

[0086] Figure 7(a) is a schematic diagram comparing the contact stress distributions calculated by the PD model and the FE model;

[0087] Figure 7(b) is a schematic diagram comparing the contact stress distributions calculated by the PD model and the FE model. Detailed implementation manners

[0088] To further understand the content of the present invention, the present invention will be described in detail with reference to the accompanying drawings and embodiments. It should be understood that the embodiments are only for explaining the present invention rather than limiting it.

[0089] Embodiment

[0090] As Figure 1 shown, this embodiment provides a cell-based peridynamic contact modeling method (EBCM method), which includes the following steps:

[0091] S1. Set and initialize parameters and variables such as materials, geometries, simulations, etc.

[0092] S2. Divide the two objects participating in the contact action into a slave object and a master object.

[0093] S3. Adopt a meshless method to discretize the slave object and the master object into peridynamic particles, define the particles that may participate in the contact interaction as contact particles, and call the contact particles on the slave object and the master object slave particles and master particles, and construct sets R s and R m .

[0094] S4. For each contact particle on the master object, that is, master particle y i construct a contact cell

[0095] In S4, each contact cell includes the current contact particle and its two adjacent contact particles, that is, any contact particle y i on the master object, and the contact cell is defined as the set of y i and its two adjacent contact particles y i1 and y i2 :

[0096]

[0097] Based on the spatial position relationship, the contact particles on the master object are divided into three categories: non-boundary particles, boundary particles, and particles adjacent to defects; for non-boundary particles, the adjacent particles of the contact cell are the two contact particles closest to it; for the particles y b, its neighboring particle y in the main object b2 It is determined by searching for the nearest contact particle, but the neighboring particles outside the main object do not exist. Therefore, a virtual neighboring particle y is introduced. bv , whose position is determined by linear extrapolation of internal neighboring particles:

[0098] y bv =2y b -y b2

[0099] For the contact particle y near the defect d , where the neighboring particle y on the left d1 On the other side of the defect, resulting in y d With y d1 There is a long distance between them, such as Figure 2 As shown, y d With y d1 The distance is 3Δx. The defect causes the contact surface to be discontinuous. If y is selected d If the particles on the left and right sides construct contact units, the contact surface will be forced to be continuous, which will also lead to inaccurate representation of the contact surface. d -y da ‖>1.5Δx d When da No longer selected as an adjacent particle. Similar to the particles at the boundary, d Introducing virtual neighboring particle y dv ,y dv The spatial position of is determined by linear extrapolation from another neighboring particle;

[0100] To ensure the accuracy and precision of contact calculation, all contact particles should have at least one neighboring particle that does not require additional processing, that is, there should be at least one particle that is connected to y i Distance Δx i adjacent particles to avoid the problem of not being able to find adjacent particles.

[0101] S5, the main particle y i and from particle y j Total normal contact force and the total tangential contact force , The value of is initialized to 0.

[0102] S6. Visit all particles y in the collection j , and traverse each particle y from the object j .

[0103] S7, by each particle y j Access all main particles y in the collection i, and traverse each primary particle y in the primary object i .

[0104] S8. Determine the primary particle y i The constructed contact element of the contact plane.

[0105] In S8, two adjacent particles in the contact element are used to determine the contact surface direction, and the contact surface position is jointly determined by three particles;

[0106] First, the contact surface direction is defined by two adjacent particles, where the vector y i2 -y i1 is the tangential direction of the contact surface, and the primary particle y i constructs a plane χ i2 -y i1 parallel to i , which is the contact plane defined by the contact element . The projection points i1 and i2 of the adjacent particles y i on χ and are the boundaries of the plane χ i ;

[0107] The unit tangential vector of the contact surface is written as:

[0108]

[0109] is determined by the following formula and to determine the spatial positions of:

[0110] k = 1 or 2

[0111] where k represents any one of the adjacent particles;

[0112] χ i is the contact plane defined by the contact element , expressed as:

[0113]

[0114] Define as the unit normal vector of χ i , pointing outside the contact surface, the unit normal vector and the unit tangential vector are perpendicular, that is, there is the following relationship:

[0115]

[0116] From particle y j At χ i The projection point is Its position is determined by the following formula:

[0117]

[0118] S9. Determine whether the condition for contact occurrence is satisfied. If satisfied, continue with subsequent operations; otherwise, terminate the contact calculation between particle y i and y j Inter - particle contact calculation

[0119] In S9, from the contact particle y j At χ i The projection point is When from the contact particle y j Regarding the contact surface χ i The projection point of At χ i Inside and the projection distance Is less than the critical value, determine the occurrence of the contact event:

[0120]

[0121] Δx i Represents the particle spacing of the main particle, Δx j Represents the particle spacing of the secondary particle; Is The projection distance projected inside χ i ; Is y j And χ i The critical contact distance between

[0122] S10. After contact occurs, calculate the material normal overlap depth Contact stiffness And the distance influence function ω ij .

[0123] In S10, the material normal overlap depth Contact stiffness And the distance influence function ω ij The calculation formulas are as follows:

[0124]

[0125] Where p si Is the control coefficient of the contact stiffness, K ij Represents the bulk modulus, δ ij Is the neighborhood radius, ΔV i And ΔV j Are the bulk moduli of the materials

[0126] S11. Calculate the normal contact force from particle y j and the main particle y i in between

[0127] In S11, the normal contact force is calculated as follows:

[0128]

[0129] f cn (y i , y j ) represents the normal contact force between y j and y i .

[0130] S12. Accumulate the calculated normal contact force from particle y j and the main particle y i in between to to obtain the total normal contact force of y i , y j .

[0131] In S12, the total normal contact force of y i is as follows:

[0132]

[0133] f cn (y i , t) represents the total normal force on the main particle y i ;

[0134] According to Newton's third law, the contact force received by particle y j is equal in magnitude and opposite in direction to that of y i , so the expression is:

[0135]

[0136] is the contact force received by the particle, is the contact force received by the main particle;

[0137] Sum up the reaction forces of particle y j and all main particles y i to obtain the total normal contact force acting on y j as:

[0138]

[0139] In the formula, R m is the set of main particles.

[0140] S13. If y j and y i come into contact for the first time, record the current moment of the slave particle y j and the master particle y i for subsequent analysis; otherwise, recalculate and the normal contact force and add it to the total normal contact force to update the contact state.

[0141] S14. Define the trial tangential contact force.

[0142] The tangential contact force is proportional to the relative sliding distance between particles. The tangential contact force is constrained by the Coulomb friction criterion. Therefore, first define a trial tangential contact force as follows:

[0143]

[0144] is the trial tangential contact force.

[0145] S15. Calculate the relative slip distance

[0146] In S15, the relative sliding distance is the difference between the tangential distances of two contacting particles from the initial contact moment t c to the current moment t, and its expression is as follows:

[0147]

[0148] S16. Define the tangential contact boundary according to the Coulomb friction criterion

[0149] In S16, the expression of the tangential contact boundary is as follows:

[0150]

[0151] μ represents the friction coefficient.

[0152] S17. Determine the tangential contact force according to the tangential contact boundary.

[0153] In S17, the expression of the tangential contact force is as follows:

[0154]

[0155] S18. Sum up all the tangential contact forces generated by the slave particles and the master particle to obtain the total tangential contact force acting on y j and y i .

[0156] In S18, the total tangential contact force acting on yi The expression for the total tangential contact force is as follows:

[0157]

[0158] Similarly, according to Newton's third law, the total tangential contact force acting on y j has the following expression:

[0159]

[0160] S19. Apply the calculated total tangential contact force to the contact particles in the form of body force density, and update the particle displacement by the numerical scheme of PD.

[0161] S20. If the termination condition (iteration time or iteration number) is satisfied, terminate the calculation; otherwise, return to S6.

[0162] This embodiment demonstrates the performance of the proposed EBCM (Element-based contact modeling, EBCM) method through two numerical examples.

[0163] The first example is the Kalthoff-Winkler experiment, which studies the fracture behavior of materials under dynamic loads and involves the impact contact between a rigid body and a flexible body. As Figure 3 shown, a 1.57 kg rigid object vertically impacts a flexible steel plate with two cracks. The steel plate is not displacement-constrained, and the two cracks are symmetrically distributed about the central axis. The material parameters of the steel plate include: Young's modulus of 191 GPa, Poisson's ratio of 0.25, density of 8000 kg / m 3 and a critical bond elongation rate of 0.01. The geometries of the impact object and the steel plate are shown in Figure 3 .

[0164] Figures 4(a) and 4(b) show the variation of the total contact force and the impact object velocity with time. The results show that at all impact velocities, the contact results calculated by EBCM are in good agreement with those of other EMU methods, verifying the effectiveness and accuracy of EBCM in the impact contact between a rigid body and a flexible body. In the first half of the simulation, due to the impact, the contact force and the compression wave increase significantly. In the second half, the contact force and the compression wave cause the impact object to separate from the steel plate, the contact force drops to zero, and the velocity of the impact object basically no longer changes.

[0165] Figure 5The crack propagation results under different impact velocities are presented. Consistent with the contact force results, the damage modes and crack propagation velocities calculated by the EBCM and EMU methods are highly consistent. Under all working conditions, due to the loading conditions and the initial crack being axisymmetric about the center, the two propagating cracks are also symmetrically distributed. In addition, the damage degree of the steel plate and the average crack propagation velocity increase with the increase of the impact velocity. When the crack length reaches a certain value, the crack propagation stops; with the increase of the impact velocity, the crack length also increases accordingly, and finally extends to the left and right boundaries of the steel plate. When the velocity further increases, additional damage appears at the bottom of the steel plate.

[0166] The second example is the friction contact test of two flexible plates, studying the influence of the friction contact and material strength of two elastic plates on the contact behavior under the condition of fixed compression displacement. As Figure 6 shown, displacement loading is applied to the top of the upper plate, while full constraints are applied to the bottom of the lower plate. The geometric shapes of the two elastic plates are exactly the same, with the same density and Poisson's ratio, but different Young's moduli: the density and Poisson's ratio of the two plates are 1 / 3 and 7850 kg / m 3 , respectively. The Young's modulus of the upper plate is set to 200 GPa, while two working conditions of 100 GPa and 200 GPa are considered for the lower plate. The friction coefficient between the two plates is set to 0.2. This embodiment does not involve material damage and can be verified by the FE (finite element) model. Due to the commercialization and standardization of the software used, when there is no involvement of damage and discontinuity, the accuracy and effectiveness of the peridynamics model can be verified by choosing to compare with the FE model. Therefore, an FE model with exactly the same parameters as the peridynamics model is also established for comparative analysis.

[0167] Figures 7(a) and 7(b) show the normal and tangential contact stress distributions calculated by the PD (the peridynamics) model and the FE model, and Table 1 records the contact force solutions under the two working conditions. The results show that the contact stress distributions calculated by the PD model and the FE model are highly consistent, and the error of the total contact force is less than 1.4%. With the increase of the material strength, the contact force also increases. When the strength increases from to, the normal contact force increases from 667.6 N to 1000.6 N, with an obvious increase. Due to the application of the fixed displacement loading condition, the deformation of the plate is relatively unchanged under the change of material strength. Under similar deformation conditions, it is reasonable that the material with greater strength generates a greater contact force.

[0168] Table 1 PD model solutions and FE model solutions of the total contact force components

[0169]

[0170] In summary, the element-based peridynamics contact modeling method proposed in this embodiment is applicable to the study of contact damage problems and various contact scenarios.

[0171] The above has schematically described the present invention and its embodiments. This description is not restrictive. What is shown in the drawings is only one of the embodiments of the present invention, and the actual structure is not limited thereto. Therefore, if those of ordinary skill in the art are inspired by this and, without departing from the gist of the present invention, design structures and embodiments similar to this technical solution without creative efforts, they shall fall within the protection scope of the present invention.

Claims

1. A cell-based peridynamic contact modeling method, characterized in that: It includes the following steps: S1. Set and initialize materials, geometries, simulation parameters and variables; S2. Divide the two objects participating in the contact into the slave object and the master object; S3. Adopt the meshless method to discretize the slave object and the master object into peridynamic particles, define the particles that may participate in contact interaction as contact particles, and the contact particles on the slave object and the master object are called slave particles and master particles, and construct the sets \(R\) of slave particles and master particles s and \(R\) m ; S4. For each contact particle on the main object, i.e., the main particle y i Construct a contact unit In S4, each contact unit includes the current contact particle and its two adjacent contact particles, i.e., any contact particle y on the main object i , the contact unit is defined as y i and its two adjacent contact particles y i1 and y i2 : S5. Initialize the numerical values of the total normal contact force i and the total tangential contact force j exerted on the primary particle y and the secondary particle y to 0; S6. Access all particles y in the set and traverse each particle y in the object j ; j ; S7. From each slave particle y j access all master particles y in the set i , and traverse each master particle y in the master object i ; S8. Determine the main particle y i The constructed contact unit Contact plane of S9. Determine whether the conditions for contact occurrence are met. If they are met, continue with subsequent operations; otherwise, terminate particle y i and y j for contact calculation S10. Calculate the normal overlapping depth of the material after contact occurs Contact stiffness and the distance influence function ω ij ; S11. Calculate the normal contact force between particle y j and the main particle y i ​ S12. The calculated normal contact force between particle y j and the main particle y i is accumulated to to obtain the total normal contact force of y i and j y; S13. If y j and y i come into contact for the first time, record the slave particle y j and the master particle y i at the current moment for subsequent analysis; otherwise, recalculate and the normal contact force and accumulate it into the total normal contact force to update the contact state; S14. Define the trial tangential contact force; S15. Calculate the relative slip distance S16. Define the tangential contact boundary according to the Coulomb friction criterion S17. Determine the tangential contact force according to the tangential contact boundary; S18. Sum up all the tangential contact forces generated by the slave particles and the master particle to obtain the total tangential contact force acting on y j y i ; S19. Apply the calculated total tangential contact force to the contact particles in the form of body force density, and update the particle displacement by the numerical scheme of PD; S20. If the termination condition is satisfied, terminate the calculation; otherwise, return to S6.

2. The unit-based peridynamic contact modeling method according to claim 1, wherein: In S4, based on the spatial position relationship, the contact particles on the main object are divided into three categories: non-boundary particles, boundary particles, and particles adjacent to defects; for non-boundary particles, the adjacent particles of the contact unit are the two closest contact particles to it; for the particles y at the left and right boundaries b , its adjacent particles y within the main object b2 are determined by searching for the closest contact particles to it, but there are no adjacent particles outside the main object. Therefore, virtual adjacent particles y bv are introduced, and their positions are determined by linear extrapolation from the internal adjacent particles: y bv = 2y b -y b2 For the contact particle y near the defect d , it is y d that a virtual adjacent particle y dv , y dv is introduced, and the spatial position of y is determined by linear extrapolation through another adjacent particle; All contacting particles should have at least one adjacent particle that requires no additional processing, i.e., there should be at least one adjacent particle that is Δx i away from y i to avoid the problem of not being able to find an adjacent particle.

3. A unit-based peridynamic contact modeling method according to claim 2, characterized in that: In S8, two adjacent particles in the contact unit are used to determine the contact surface direction, while the contact surface position is jointly determined by three particles; First, the contact surface direction is defined by two adjacent particles, where the vector y i2 -y i1 is in the tangential direction of the contact surface. The main particle y i constructs a plane χ i2 parallel to y i1 -y i , which is the contact plane defined by the contact element . The projection points of the adjacent particles y i1 and y i2 on χ i are and , which are the boundaries of the plane χ i ; Contact surface unit tangential vector Written as: Determined by the following formula and for the spatial position: In the formula, k represents any adjacent particle; χ i is the contact unit defines a contact plane, expressed as: Definition is χ i 's unit normal vector, pointing outside the contact surface, the unit normal vector and the unit tangential vector are perpendicular, that is, there is the following relationship: From particle y j At χ i The projection point is Its position is determined by the following formula:

4. The method for modeling the peridynamic contact based on units according to claim 3, wherein: In S9, from the contacting particle y j at χ i the projection point is When from the contacting particle y j with respect to the contact surface χ i the projection point is within χ i and the projection distance is less than the critical value, it is determined that a contact event has occurred: Δx i represents the particle spacing of the main particles, Δx j represents the particle spacing of the secondary particles; is the projection distance within the χ i projection; is j the critical contact distance between y i and χ.

5. A unit-based peridynamic contact modeling method according to claim 4, characterized in that: In S10, the material normal overlap depth Contact stiffness and the distance influence function ω ij are calculated as follows: where p si is the control coefficient of contact stiffness, K ij represents the bulk modulus, δ ij is the neighborhood radius, ΔV i and ΔV j are the bulk moduli of the material.

6. The unit-based peridynamic contact modeling method according to claim 5, wherein: In S11, the normal contact force is calculated as follows: f cn (y i ,y j ) represents the normal contact force between y j and y i and y 7. A unit-based peridynamic contact modeling method according to claim 6, characterized in that: In S12, y i The total normal contact force is as follows: f cn (y i , t) represents the total normal force on the primary particle y i ; According to Newton's third law, from particle y j the contact force received is equal in magnitude and opposite in direction to that of y i so the expression is: is the contact force received by the particle, is the contact force received by the main particle; will be from particle y j and sum the reaction forces with all the primary particles y i to obtain the total normal contact force acting on y j as follows: where R m is the set of primary particles.

8. A unit-based peridynamic contact modeling method according to claim 7, characterized in that: In S14, the trial tangential contact force is as follows: is the tangential contact force for trial cutting.

9. A cell-based peridynamic contact modeling method according to claim 8, characterized in that: In S15, the relative sliding distance is the difference in the tangential distances of the two contacting particles from the initial contact time t c to the current time t, and its expression is as follows: In S16, the tangential contact boundary has the following expression: μ represents the friction coefficient; In S17, the expression of the tangential contact force is as follows:

10. A unit-based peridynamic contact modeling method according to claim 9, characterized in that: In S18, the expression of the total tangential contact force acting on y i is as follows: Similarly, according to Newton's third law, the expression for the total tangential contact force acting on y j is as follows:

Images