Numerical modeling method for 3D printing of micro-structure sneaker sole

Abstract

The invention relates to the technical field of 3D printing, in particular to a numerical modeling method for 3D printing of a micro-structure sports shoe sole. The invention relates to a numerical modeling method for a 3D printing micro-structure sports shoe sole, and the method comprises the following steps: (a) building a three-dimensional model according to the overall structure of the 3D printing shoe sole and the size of a micro-structure; (b) carrying out unit grid division according to the 3D printing sole integral structure and the three-dimensional model of the microstructure, and assembling the model; a display dynamics method is added on ABAQUS software through a time step module, and the method has the advantages that test analysis of the 3D printing shoe sole microstructure under the complex working condition is achieved, and the structure under the specific working condition can be rapidly and accurately analyzed through a numerical method after verification of the simple working condition; secondly, a numerical simulation method is adopted for achieving the mechanical property of the 3D printing shoe sole microstructure, results of different parameters are provided in batches, more references are provided, and a large amount of cost is saved.

Description

Technical Field

[0001] The present invention relates to the technical field of 3D printing, and particularly to a numerical modeling method for 3D printing the microstructure of sports shoe soles. Background Art

[0002] At present, the technology of 3D printing shoe soles has been widely applied. However, the existing analysis methods for the microstructures of 3D printed shoe soles often require preparing 3D printed structures through equipment and designing corresponding tests, such as compression tests, linear friction tests, etc., to conduct corresponding test measurements. The tests for analyzing the complex working conditions of 3D printed shoe sole microstructures cannot be achieved, which will lead to certain limitations in the analysis of 3D printed shoe sole microstructures under complex conditions.

[0003] Secondly, when conducting experimental analysis on the microstructures of 3D printed shoe soles, it is often necessary to design tests to specifically analyze the performance of the test structures. When analyzing the specific parameters of 3D printed shoe sole microstructures, a large number of experimental results are often required for support, resulting in excessively high actual costs and making it difficult to meet the parameter optimization and improved design requirements in the design and analysis work of 3D printed shoe sole microstructures.

[0004] Therefore, in view of the many drawbacks of the existing 3D printed shoe sole microstructure analysis technologies, such as limitations in working condition analysis and batch analysis limitations, etc., the present invention proposes a numerical modeling method for 3D printed microstructure sports shoe soles. Summary of the Invention

[0005] To solve the above problems, the present invention provides a numerical modeling method for 3D printed microstructure sports shoe soles, and the specific scheme is as follows: A numerical modeling method for 3D printed microstructure sports shoe soles includes the following steps: (a) Establish a three-dimensional model according to the overall structure of the 3D printed shoe sole and the dimensions of the microstructure; (b) Perform unit mesh division and assemble the model according to the overall structure of the 3D printed shoe sole and the three-dimensional model of the microstructure; (c) On the ABAQUS software, that is, the Abaqus software, add the explicit dynamics method through the time step module, apply corresponding boundary conditions to the three-dimensional model established in (a), and assign material parameters based on experimental results and reference documents; add the contact conditions between models; the reference documents at least include the test data obtained from uniaxial tensile tests, uniaxial compression tests, equibiaxial tensile tests, pure shear tests, and volume deformation tests on the test materials, as well as the parameters obtained by fitting the test data with the Neo-Hookean, Mooney-Rivlin, Ogden, and Yeoh theoretical models; (d) Construct the dynamic equilibrium equation according to the explicit dynamics method; (e) Solve the dynamic equilibrium equation in (d) using the central difference solution method; (f) Perform contact determination based on the added contact boundaries; (g) Calculate the relevant physical quantities of the 3D printed sole microstructure, where the relevant physical quantities at least include the displacements and velocities of the nodes in the sole microstructure; (h) Conduct output analysis on the physical quantities calculated in (g).

[0006] By combining 3D printing technology to precisely design the sole microstructure and establish a three-dimensional model, during the model establishment process, establishing a corresponding mechanical model according to the specific actual working conditions is a modeling factor that must be considered. Under working conditions such as walking and standing, the sole is subjected to vertical impact forces and lateral drag forces, resulting in compression and shear deformations of the sole structure, causing extrusion damage, slipping, and wear of the sole under walking conditions; at the same time, material parameters, structural shapes, and contact surface properties, etc., will significantly affect the mechanical properties of the sole structure and further affect the service performance of the sole.

[0007] Secondly, analyze and optimize the size and layout of the microstructure to maximize the mechanical properties of the sole. Then, by collecting experimental data and combining relevant theories such as elasticity, hyperelasticity, and tribology, evaluate and analyze the mechanical property data of the 3D printing material, such as elastic modulus, Poisson's ratio, friction coefficient, etc. Further, through numerical analysis methods, finally establish a mechanical model of the 3D printed sole microstructure.

[0008] The Neo-Hookean, Mooney-Rivlin, Ogden, and Yeoh theoretical models are all mechanical models describing hyperelastic materials, and their Chinese translations are: Neo-Hookean theoretical model, Mooney-Rivlin theoretical model, Ogden theoretical model, and Young's theoretical model respectively.

[0009] Further, it also includes discretizing the unit meshes divided in step (b) into uniform structured meshes by adopting FEM element types; after mesh division, enter the assembly module to complete the assembly of all models; the FEM element types include S3, S4R, C3D4, and C3D8R.

[0010] Among them, FEM is the abbreviation of the finite element method, and S3, S4R, C3D4, and C3D8R are several typical finite element cells; and the FEM cells are the cell types carried by the ABAQUS software itself. Among them, the S3 cell represents a three-node shell cell, the S4R cell represents a four-node shell cell, the C3D4 cell represents a four-node tetrahedral cell, and the C3D8R cell represents an eight-node hexahedral cell.

[0011] Furthermore, in step (c), the applied boundary conditions include: applying a gravity load to the sole model, applying a vertical load and a horizontal tangential load to the shoe sole model, and applying a fixed constraint to the ground model to ensure that the ground remains stationary.

[0012] Apply the corresponding boundary conditions according to the actual working condition model in (a): for example, apply gravity to the sole model, apply vertical and horizontal loads to the shoe sole model respectively, and apply a fixed constraint to the ground model.

[0013] According to the material experiment results and the literature material parameters: for example, use the incompressible material model of Young's model Yeoh for the rubber compound model in the literature, and set the material model according to the relevant parameters; according to the stress-strain curve of the material obtained from the experiment, import it into Abaqus, that is, in the ABAQUS software, fit the corresponding material model and material parameters, and set the material model according to the relevant parameters.

[0014] Furthermore, it also includes that after discretizing the structure, the dynamic equilibrium equations of each node in the motion state are constructed as shown in formula (1): (1) Among them, {} represents the vector of each parameter, represents the inertial force, with the unit of N; represents the damping force, with the unit of N; represents the external dynamic load; with the unit of N; represents the structural elastic force, with the unit of N; respectively represent their vectors, and the same applies hereinafter; Among them, the structural elastic force vector is represented by the node displacement vector and the stiffness matrix as ; According to D'Alembert's principle, the inertial force vector is represented by the node acceleration vector and the mass matrix vector where , then it is ; The damping force vector is represented by the node velocity vector and the damping matrix as , where ; The external dynamic load vector is represented by the external force acting acceleration vector and the mass matrix as , therefore, formula (1) is sorted as follows: (2) {} represents the vector of each parameter, which is represented as the nodal displacement vector, with the unit of m, which is represented as the nodal velocity vector, unit: m / s, which is represented as the nodal acceleration vector, unit: m / s², which is represented as the acceleration vector of the external force, unit: m / s².

[0015] Furthermore, the formula (2) is calculated using the central difference solution method, and the calculation method is as follows: (s1) Define any time period as Δ t , and adopt an equal time step Δ t i =Δ t Then t i the nodal velocity vector and the nodal acceleration vector at time (3) (4) (s2) Substitute the above formula into formula (2) to obtain t i the dynamic equilibrium equation at time t i+1 where the position of the node at time (5) In the formula, is t i the stiffness at time , with the unit of N / m, t i is the vector of the external load at time , with the unit of N, t i+1 is the vector of the nodal displacement at time where: (6) (7) where the mass matrix M represents the relationship between the mass of each object and its physical quantities of position and velocity, unit: kg; the damping matrix C represents the damping effect caused by the friction, viscosity or other forms of energy dissipation of the system during the dynamic response of the structure and mechanical system objects, unit: kg / s.

[0016] Furthermore, the contact conditions described in step (c) include normal contact and frictional contact.

[0017] Among them, the normal contact formula needs to be defined by a distance function, and the frictional contact needs to be defined by the relative tangential displacement; the additional contact conditions need to set contact options in the contact module of the ABAQUS software. First, set the properties of the normal contact and frictional tangential contact, and then select the master and slave contact surfaces. The judgment and calculation are carried out through the built-in functions of the software; the contact conditions are specifically the judgment functions used when independent models interact, ensuring the transfer of parameters such as force and displacement between models.

[0018] Furthermore, step (f) includes the following situations: When making contact with the physical model, it is necessary to judge the contact boundary. After contact, relevant parameters are analyzed. The relevant parameters include at least contact force, contact normal force, contact tangential force, contact stress, and contact displacement. Then, the judgment of the conditions is further calculated through equations. The judgment method is as follows: Step 1: According to the set master-slave surface relationship, call the function to find the node closest to the master-slave surface. The contact point on the slave surface is known X The corresponding unknown contact point on the master surface for contact point , as shown in the following formula: (8) Among them, X 1 represents the contact point on the slave surface, represents the unknown contact point on the master surface, represents the set range of the master surface; represents one of the unknown contact points on the master surface and is the closest to X 1, represents X 1 and take the modulus, and it is the shortest distance between the two points; Step 2: Calculate the distance between the contact pairs: By calculating the X 1 and X 2 The distance between the two points is used to judge whether contact occurs, so as to apply contact constraints. The distance between the two points is defined by the following formula: (9) Among them: is the distance between the two points, with the unit of m; is the unknown contact point on the master surface 's tangential vector; It can be obtained from the above formula that when >0, there is no contact; =0, perfect contact; <0, penetration, that is, contact is activated; Step 3: Obtain the equation of the governing equation through the contact surface equilibrium equation: (10) Wherein, is the integral term of the internal force term, is the integral term of the contact stress term, is the integral term of the boundary condition term, is the integral term of the body force term; Wherein, In the formula represents the inner product of stress and displacement gradient, with the unit of N / m 2 , represents the infinitesimal volume element, with the unit of m 3 ; In the formula represents the contact stress, with the unit of N / m 2 , represents the contact displacement increment, with the unit of m, represents the infinitesimal contact surface area, with the unit of m 2 ; In the formula represents the force acting on the boundary, with the unit of N / m 2 , represents the displacement increment, with the unit of m, represents the acting boundary area, with the unit of m 2 ; In the formula represents the body force per unit volume, which is the external force acting on the structure by the environment or external conditions, with the unit of N / m 3 , represents the displacement increment, with the unit of m; As can be seen from formula (2), the contact term in the governing equation is always non-positive. The penalty function algorithm is adopted, and this function satisfies: a. When the distance between the contacting bodies is greater than zero, the function value is 0; b. After contact occurs, the function value increases monotonically with the contact penetration distance. At this time, the contact stress integral term needs to be rewritten. After rewriting, it is as follows: (11) Wherein, k is the penalty factor in the penalty function, with the unit: N / m³, which controls the magnitude of the contact stiffness; d is the contact gap on the contact surface, which is the distance between two contact points, with the unit: m; When the objects are in contact, , when there is no contact, .

[0019] Further, the relevant physical quantities described in step (g) include the nodal displacement vector , the node velocity vector , the node acceleration vector For these three parameters, and then through the above three parameters combined with the reference motion equation formula (2) and equation (10), further derivation is carried out to obtain the required physical quantities, which include node strain, node stress, and node force.

[0020] Furthermore, the calculation of the required physical quantities of node strain, node stress, and node force is specifically as follows: a) The calculation method of node strain is as follows: A displacement function is established for the element nodes to obtain the shape function reflecting the element displacement form. Combining the shape function, the strain components in the directions of the element coordinate axes can be expressed as follows: (12) In the formula, the [B] matrix is determined by the element area or element volume of the mesh division and the coefficient constant; it is used to relate the displacement gradient to the strain; represents the vector of the strain at the node, with the unit of m / m; b) The calculation method of node stress is as follows: According to Hooke's law, a stress-strain relationship is established for each component direction of the node as follows: (13) In the formula, is the elastic matrix, with the unit of N / m³, represents the vector of the stress at the node, with the unit of N / m 2 ; is determined by the elastic constant E and Poisson's ratio μ of the material at the node; represents the vector of the strain at the node, with the unit of m / m; c) The calculation method of node force is as follows: The boundary force is equivalently concentrated on the node to ensure that the direction is consistent with the node displacement direction. According to the principle of static equivalence, a relationship between the element node force and the node displacement is established as follows: (14) In the formula, K is the stiffness matrix, which is composed of the [B] of the node and the elastic matrix , { F} represents the vector of the node force, with the unit of N / m.

[0021] Furthermore, it also includes that after obtaining the relevant physical quantities for the solution, the relevant data is processed through the post-processing function of ABAQUS, and the relevant parameters are output and the stress nephogram, contact force nephogram, and contact area nephogram are drawn.

[0022] After obtaining the relevant physical quantities to be solved, such as displacement, nodal force, etc., the relevant data is processed through the post-processing function of ABAQUS to output the relevant parameters such as stress, contact force, contact area, etc. of the 3D printed sole structure, and draw stress nephograms, contact force nephograms, contact surface nephograms, etc.

[0023] The advantages of the present invention are as follows: This method realizes the test analysis of the 3D printed sole microstructure under complex working conditions. After verification through simple working conditions, the structure under specific working conditions can be analyzed quickly and accurately by numerical methods, providing a fast and highly accurate analysis method compared with traditional experimental design test methods; Secondly, the numerical simulation method is used to realize the mechanical properties of the 3D printed sole microstructure, and the results of a batch of different parameters can be provided, providing more references for the parameter optimization analysis and improvement design of the structure, saving a large amount of costs compared with traditional methods. Description of the Drawings

[0024] Figure 1 is the method flow chart of the present invention; Figure 2 is the three-dimensional model schematic diagram of the sole - shoe sole - ground in Embodiment 1; (a) Schematic diagram of the combined model of the sole, shoe sole and ground; (b) Schematic diagram of the midsole model; (c) Schematic diagram of the outsole model; Figure 3 is the friction force and contact area curve diagram during rubber friction slip under different parameters in Embodiment 1; (a) Verification results of simulation and experiment: (b) Friction force - displacement curve diagram under different pattern conditions: (c) Relative contact area - displacement curve diagram under different pattern conditions: (d) Friction force - displacement curve diagram under different working conditions: (e) Relative contact area - displacement curve diagram under different working conditions: Figure 4 is the stress, strain, contact normal force, frictional shear force, relative slip and contact area nephogram during rubber friction slip with rectangular patterns in Embodiment 1; Figure 5 is the deformation nephogram of the sole - shoe sole - ground model under the full - foot landing working condition. Detailed Embodiment

[0025] Embodiment 1 A numerical modeling method for a 3D printed microstructure sports shoe sole includes the following steps: (a) Establish a three - dimensional model according to the overall structure of the 3D printed sole and the dimensions of the microstructure; (b) Perform element mesh division and assemble the model according to the overall structure of the 3D printed sole and the three - dimensional model of the microstructure; (c) Add the explicit dynamics method through the time step module in the ABAQUS software, apply corresponding boundary conditions to the three-dimensional model established in (a), and assign material parameters based on the experimental results and references; add the contact conditions between the models; the references at least include the test data obtained from uniaxial tensile tests, uniaxial compression tests, equibiaxial tensile tests, pure shear tests, and volume deformation tests on the test materials, as well as the parameters obtained by fitting the test data with the Neo-Hookean, Mooney-Rivlin, Ogden, and Yeoh theoretical models; (d) Construct the dynamic equilibrium equation according to the explicit dynamics method; (e) Solve the dynamic equilibrium equation in (d) using the central difference solution method; (f) Perform contact determination based on the added contact boundaries; (g) Calculate the relevant physical quantities of the 3D printed sole microstructure, where the relevant physical quantities at least include the displacement and velocity of the nodes in the sole microstructure; (h) Perform output analysis on the physical quantities calculated in (g).

[0026] Next, the numerical modeling method of the 3D printed microstructure sports shoe sole of the present invention is verified: Figure 2 It is a schematic diagram of the three-dimensional model of the foot-sole-ground. This model is composed of the foot, the midsole of the sole, the outsole of the sole, and the ground model. The foot, the sole, and the ground are arranged and combined from top to bottom. The midsole model uses four-node linear tetrahedral elements, namely C3D4 meshes. The basic size of the mesh is 0.2 mm, the structural size is 250 mm long × 100 mm wide, the minimum height is 7 mm, and the maximum height is 27 mm. The material used is rubber material; The outsole model uses C3D8R, that is, a hexahedron linear reduced integration element mesh. The minimum size of this mesh is 0.25 mm, its size is 250 mm long × 100 mm wide × 3 mm high. Among them, the size of the outsole pattern unit cell model is 6 mm long × 5 mm wide × 3 mm high, arranged at intervals of 2 mm. The material used is polymer composite resin material.

[0027] The foot model uses S4R, that is, a 4-node quadrilateral linear shell element mesh, with a basic size of 3 mm. The ground model uses C3D8R, that is, a hexahedron linear reduced integration element mesh, with a basic size of 5 mm. The foot model is set as a rigid body and a gravity load of 700 N is added; the ground model is set as a rigid body, and the boundary condition of fixed constraint is set. The outer side of the midsole model of the shoe is subjected to displacement constraint, and an external load of 200 mm / s is applied to the outsole model.

[0028] By setting the explicit dynamics method in the ABAQUS software, the software can automatically construct the dynamic equilibrium equation and use the central difference method to solve the equation.

[0029] At the same time, the general contact is set. By setting the explicit dynamics method, the software can judge the contact boundary of the added general contact. According to the adopted standard Coulomb friction model, the software sets the contact condition as the penalty function model for calculation.

[0030] Figure 3 It shows the friction force and contact area curve during the frictional slip of rubber under different parameter conditions. It can be seen from the figure that the frictional simulation results of the rubber block can accurately verify the test results and ensure the accuracy of the simulation analysis; under the pattern condition, the rubber block structure is deformed by friction, and the increase of the extrusion edge provides higher friction force, and the increase of the non-extrusion edge provides more relative contact area.

[0031] Under the working condition, the increase of the actual contact area of the outsole rubber structure provides more friction force, and the contact pressure of the actual contact surface decreases, reducing the relative contact area.

[0032] Figure 4 It shows the stress, contact normal force, frictional shear force, relative slip and contact surface nephogram during the frictional slip of rubber under the rectangular pattern condition.

[0033] It can be seen from the figure that the structure extrusion edge is bent and deformed due to friction induction, and the stress and strain are concentrated on the contact edge of the friction-induced bend; the contact normal force and frictional shear force are also concentrated on the contact surface between the friction-induced bend edge and the ground; during frictional sliding, the extrusion contact edge is bent and twisted, resulting in the front vertical surface continuing to contact the ground and the rear non-bent contact surface having relative sliding; due to the bending deformation of the extrusion contact edge, the contact surfaces at the extrusion end and the non-extrusion end are continuously reduced until the front vertical surface continuously contacts the ground.

[0034] Figure 5 It shows the deformation nephogram of the sole-shoe outsole-ground model under the full-foot contact working condition. The full-foot contact working condition mainly simulates the external load on the shoe outsole when a person is walking; it can be seen from Figure (a) that the midsole structure of the shoe is significantly deformed under the sole gravity load, and the stress and strain are mainly concentrated at the midfoot arch position in the middle of the midsole; it can be seen from Figure (b) that the outsole structure of the shoe is subjected to a velocity external load in the horizontal direction, and the stress, strain, contact normal force and frictional tangential force are concentrated on the extrusion contact edge; during frictional sliding, it causes the shear deformation of the pattern structure, and local contact surface detachment occurs at the non-extrusion contact edge, while the extrusion end remains in contact and has relative sliding.

Claims

1. A numerical modeling method for 3D printed microstructure sports shoe soles, characterized in that: The following steps are included: (a) Establish a three-dimensional model based on the overall structure and microstructure dimensions of the 3D printed sole; (b) Divide the unit grid and assemble the model according to the three-dimensional model of the overall structure and microstructure of the 3D printed sole; (c) Adding a display dynamics method through the time step module on the ABAQUS software, applying corresponding boundary conditions to the three-dimensional model established in (a), and assigning material parameters based on experimental results and references; adding contact conditions between models; the references at least include test data obtained by performing uniaxial tension tests, uniaxial compression tests, equibiaxial tension tests, pure shear tests, and volume deformation tests on the test materials, and parameters obtained by fitting the test data curves using Neo-Hookean, Mooney-Rivlin, Ogden, and Yeoh theoretical models; (d) construct the dynamic equilibrium equation according to the explicit dynamics method; (e) Solve the dynamic equilibrium equation in (d) using the central difference solution method; (f) Contact determination based on the added contact boundaries; (g) calculating relevant physical quantities of the 3D printed sole microstructure, wherein the relevant physical quantities at least include displacement and velocity of nodes in the sole microstructure; (h) Output analysis of the physical quantities calculated in (g).

2. The numerical modeling method of a 3D printed microstructure sports shoe sole according to claim 1, characterized in that: The method also includes discretizing the unit grid divided in step (b) into a uniform structured grid by adopting a FEM finite element method unit type; After meshing, enter the assembly module to complete the assembly of all models; the FEM finite element method unit types include S3, S4R, C3D4, and C3D8R.

3. The numerical modeling method of a 3D printed microstructure sports shoe sole according to claim 1, characterized in that: In step (c), the boundary conditions applied include: applying a gravity load to the sole model, applying a vertical load and a horizontal tangential load to the sole model, and applying a fixed constraint to the ground model to ensure that the ground is stationary.

4. The numerical modeling method of a 3D printed microstructure sports shoe sole according to claim 2, characterized in that: It also includes: after discretizing the structure, constructing the dynamic balance equation of each node in the motion state as shown in formula (1): （1） Among them, {} represents the vector of each parameter, Represents inertial force, unit is N; Indicates the damping force, in N; Indicates external dynamic load; the unit is N; Represents the structural elastic force, in N; Then they represent their vectors respectively, the same below; Among them, the structural elastic force vector Use node displacement vector and the stiffness matrix vector Indicates that ; According to D'Alembert's principle, the inertial force vector Using the node acceleration vector and the mass matrix Indicates that , then ; Damping force vector Using the node velocity vector and the damping matrix Indicates that ,in ; External dynamic load vector Using external force to apply acceleration vector and the mass matrix vector Indicates that , therefore, formula (1) is organized as follows: (2) Among them, {} represents the vector of each parameter, Expressed as a node displacement vector in m, Expressed as a node velocity vector, unit: m / s, Expressed as a node acceleration vector, unit: m / s², Expressed as the acceleration vector of the external force, unit: m / s².

5. The numerical modeling method of a 3D printed microstructure sports shoe sole according to claim 4, characterized in that: The central difference solution method is used to calculate formula (2), and the calculation method is as follows: (s1) Define any time period as Δ t , using equal time step Δ t i =Δ t but t i The vector of node velocity at time and the vector of the node acceleration The central difference of is approximately: （3） （4） (s2) Substitute the above formula into formula (2) and we get t i The dynamic balance equation at the moment t i+1 The position of the node at a given moment is derived by: （5） Among them, [] represents the matrix of each parameter, for t i The matrix of moment stiffness, unit: N / m, for t i The vector of the external load at the moment, unit: N, for t i+1 The vector of node displacement at the moment, unit: m; in: （6） （7） Among them, the mass matrix [ M ]: represents the relationship between the mass of each object and its position and velocity, unit: kg; damping matrix [ C ] represents the damping effect caused by friction, viscosity or other forms of energy dissipation in the dynamic response of a structure or mechanical system object, unit: kg / s.

6. The numerical modeling method of a 3D printed microstructure sports shoe sole according to claim 1, characterized in that: The contact conditions described in step (c) include normal contact and friction contact.

7. The numerical modeling method of a 3D printed microstructure sports shoe sole according to claim 1, characterized in that: Step (f) includes the following situations: when the physical model is in contact, the contact boundary needs to be determined, and relevant parameters are obtained by analysis after contact. The relevant parameters include at least contact force, contact normal force, contact tangential force, contact stress, and contact displacement. The determination of the conditions is further calculated through equations, and the determination method is as follows: Step 1: Based on the master-slave relationship that has been set, call the function to find the node closest to the master and slave surfaces. The contact points on the slave surface are known. X 1 corresponds to the unknown contact point on the main surface , as shown below: （8） in, X 1 represents the contact point on the surface, Indicates unknown contact points on the main surface. Indicates the set main surface range; Indicated as an unknown contact point belonging to the main surface One of the contact points of X 1 is the closest distance, express X 1 and Take the modulus, and it is the shortest distance between two points; Step 2: Calculate the distance between contact pairs: By calculating X 1 and X 2 The distance between two points is used to determine whether contact occurs, thereby applying contact constraints, and the distance between two points is defined by the following formula: （9） in: is the distance between two points, in meters; Unknown contact point on main surface The tangent vector of ; From the above formula, we can get >0, no contact; =0, perfect contact; <0, penetration, i.e. contact activation; Step 3: Obtain the equation of the governing equation through the contact surface equilibrium equation: （10） in, is the integral term of the internal force, is the integral term of the contact stress term, is the boundary condition integral term, It is the integral item of physical strength; in, in the formula Expressed as the inner product of stress and displacement gradient, in N / m 2 , Expressed as a small volume unit, the unit is m 3 ; in the formula Indicates contact stress in N / m 2 , Represents the contact displacement increment in m. Expressed as the small area of ​​the contact surface, unit is m 2 ; in the formula Represents the force acting on the boundary, in N / m 2 , Represents displacement increment, unit is m, Expressed as the boundary area of ​​action, in m 2 ; in the formula It represents the body force per unit volume. It is the external force acting on the structure by the environment or external conditions. The unit is N / m 3 ; From formula (2), we can see that the contact term in the control equation is always not positive. The penalty function algorithm is used. This function satisfies the following conditions at the same time: a. When the contact body distance is greater than zero, the function value is 0; b. After contact occurs, the function value increases monotonically with the contact penetration distance. At this time, the contact stress term integral needs to be rewritten: , rewritten as follows: （11） in, k is the penalty factor in the penalty function, unit: N / m³, which controls the size of the contact stiffness; d is the contact gap on the contact surface, which is the distance between two contact points, in meters; when the objects are in contact, , when there is no contact, .

8. The numerical modeling method for 3D printed microstructure sports shoe soles according to claim 1, characterized in that: The relevant quantities described in step (g) include the vectors of node displacements , the vector of node velocity , the vector of the node acceleration These three parameters are then further deduced through the above three parameters combined with the reference motion equation formula (2) and equation (10) to obtain the required physical quantities, which include node strain, node stress, and node force.

9. The numerical modeling method of a 3D printed microstructure sports shoe sole according to claim 8, characterized in that: The calculation of the required physical quantities, nodal strain, nodal stress, and nodal force are described as follows: a): The calculation method for node strain is as follows: establish a displacement function for the unit node, and obtain a shape function that reflects the unit displacement shape. Combined with the shape function, the strain component in the unit coordinate axis direction can be expressed as follows: （12） Where [B] is the matrix determined by the unit area or unit volume of the grid and the coefficient constant; it is used to relate the displacement gradient to the strain; It is expressed as the vector of strain at the node, in m / m; b): The calculation method for node stress is as follows: According to Hooke's law, the stress-strain relationship is established for the components in each direction of the node, as follows: （13） In the formula, is the elastic matrix, in N / m³, Expressed as a vector of stress at the node in N / m 2 ; Determined by the elastic constant E and Poisson's ratio μ of the material at the node; It is expressed as the vector of strain at the node, in m / m; c): The calculation method for the node force is as follows: the boundary force is equivalently concentrated on the node to ensure that the direction is consistent with the node displacement direction. According to the static equivalence principle, the relationship between the unit node force and the node displacement is established as follows: （14） In the formula, [ K ] is the stiffness matrix, which is composed of the [B] matrix and the elastic matrix composition,{ F } is represented as the vector of nodal force, with unit of N.

10. The numerical modeling method of a 3D printed microstructure sports shoe sole according to claim 9, characterized in that: It also includes, after obtaining the relevant physical quantities to be solved, processing the relevant data through the ABAQUS post-processing function, outputting the relevant parameters and drawing the stress cloud map, contact force cloud map, and contact area cloud map.

Images