A continuum phase field damage and fracture simulation method combining finite element method and material point method

Abstract

This invention discloses a method for simulating phase-field damage and fracture in continuums by combining finite element method (FEM) and material point method (MPD). A numerical model of the continuum to be simulated is established through 3D scanning. The fundamental physical and mechanical parameters of the continuum are measured, and an Eulerian background mesh and a finite element solid mesh are constructed. Based on the fundamental physical and mechanical parameters, the physical material parameters of the finite element mesh are defined, and the material parameters, phase-field parameters, and finite element Gaussian point parameters are initialized. Combining the finite element method, the MPD, and phase-field fracture theory, and employing the USF update step, the method accurately and effectively simulates the dynamic crack propagation process in continuums. This provides a new method and approach for simulating crack propagation in continuums, avoiding the limitations of the finite element method in strongly nonlinear problems, and overcoming the numerical fracture problem in the MPD calculation process. It eliminates the need for continuous crack propagation path tracking, effectively improving the efficiency of crack propagation simulation.

Description

Technical Field

[0001] This invention relates to the field of continuum mechanics, specifically a method for simulating continuum phase field damage and fracture by combining finite element method and material point method. Background Technology

[0002] Under external loads, continuums often undergo significant deformation processes. These processes not only exhibit complex stress and strain responses but can also lead to local instability, causing rapid crack propagation and ultimately resulting in the material losing its load-bearing capacity, significantly impacting the safety of engineering structures. Numerical simulation of fracture behavior in continuums, a highly discontinuous problem, is currently an important tool for studying its mechanisms.

[0003] Currently, continuum fracture simulation methods can be broadly categorized into continuous medium methods and discrete medium methods. The former has seen significant improvements and attempts in simulating discontinuous problems, but for complex rock fracture problems, the continuity assumption still struggles to truly address discontinuities. The latter, while overcoming the continuity limitation, suffers from unsatisfactory results in rigorous mechanical proofs and boundary condition handling. The material point method, based on the continuity assumption, avoids numerical dissipation during convection term averaging by employing particle integrals, giving it an advantage in simulating discontinuous problems. However, the material point method, essentially a particle method, still faces common problems inherent in particle methods, such as the application of boundary conditions, numerical fracture caused by particles crossing the grid, and Gaussian point migration. Furthermore, for fracture problems, most current algorithms are still based on traditional fracture mechanics theory, requiring crack propagation path tracking during simulation, resulting in low computational efficiency and an inability to simulate complex crack propagation problems. Summary of the Invention

[0004] The purpose of this invention is to provide a method for simulating continuum phase field damage and fracture by combining finite element method and material point method, comprising the following steps:

[0005] Step 1) Construct the Eulerian background mesh and finite element mesh using the material point method, and initialize the physical material parameters and Gaussian point parameters within the elements of the finite element mesh;

[0006] Step 2) Based on the physical material density and phase field viscous damping of the finite element mesh, calculate the nodal mass and phase field coefficient matrix of the finite element nodes using the Gaussian point parameters inside the element.

[0007] Step 3) Initialize the finite element node parameters and background mesh node parameters;

[0008] Step 4) Based on the current velocity of the finite element node, update the stress and strain of the Gaussian points inside the finite element. Using the current phase field of the finite element node, calculate the phase field and phase field gradient of the Gaussian points inside the element.

[0009] Step 5) Based on the updated Gaussian point stress and strain, calculate the internal forces and external forces of the finite element nodes. Using the phase field and phase field gradient of the Gaussian point inside the element at the current moment, calculate the geometric resistance of the phase field of the finite element nodes and the driving force of the phase field of the finite element nodes.

[0010] Step 6) Apply boundary conditions to the finite element nodes;

[0011] Step 7) Map the finite element node information to the background mesh node, and calculate the background mesh node mass, background mesh node momentum, background mesh node internal force, background mesh node external force, background mesh node phase field coefficient, background mesh node phase field driving force, and background mesh node phase field geometric resistance at the current moment.

[0012] Step 8) Based on the internal forces and external forces of the background mesh nodes at the current moment, calculate the momentum increment of the background mesh nodes and update the momentum of the background mesh nodes;

[0013] Step 9) Based on the internal forces, external forces, and mass of the background mesh nodes at the current moment, calculate the velocity increment of the background mesh nodes and update the velocity of the background mesh nodes;

[0014] Step 10) Calculate the phase field increment of the background grid nodes based on the phase field coefficients, phase field geometric resistance, and phase field driving force of the background grid nodes at the current moment;

[0015] Step 11) Update the velocity and position of the finite element nodes based on the updated background mesh node velocities and the background mesh node velocity increments;

[0016] Step 12) Apply boundary conditions to the finite element nodes again;

[0017] Step 13) Update the phase field of the finite element nodes based on the phase field increment of the background mesh nodes;

[0018] Step 14) Return to Step 2), substitute the updated parameters and continue iteratively solving until the crack propagates, the continuum breaks, the internal stress approaches 0, the fracture energy is exhausted, and the energy is completely released.

[0019] Furthermore, in step 1), the steps of constructing the material point method Eulerian background mesh and the finite element mesh include:

[0020] Step 1.1) Establish a numerical model of the continuum to be simulated through 3D scanning;

[0021] Step 1.2) Measure the basic physical and mechanical parameters of the continuum;

[0022] The basic physical and mechanical parameters of the continuum include the elastic modulus E and density. Poisson's ratio and critical energy release rate ;

[0023] Step 1.3) Based on the numerical model of the continuum to be simulated and the basic physical and mechanical parameters of the continuum, construct the Eulerian background mesh and the finite element mesh using the material point method;

[0024] The physical material parameters of the finite element mesh include elastic modulus E and density. Compared to Poisson ;

[0025] The Gaussian point parameters within the unit include position. Gaussian point stress Gaussian point strain Gaussian point phase field Critical energy release rate of material cracks Phase field viscous damping Phase field characteristic width Maximum crack driving force in historical phase field .

[0026] Furthermore, in step 2), the finite element node mass... As shown below:

[0027]

[0028] In the formula, e is the finite element number, g is the Gaussian point number inside the element, and K is the node number of the finite element. For the density of a continuous material, The total number of Gaussian points within a finite cell. This is the interpolation of the Gaussian point g at the node K of the finite element. The weights of the Gaussian points, Let the Jacobian matrix be a Gaussian point;

[0029] Finite element node phase field coefficient matrix As shown below:

[0030]

[0031] In the formula, It is a Gaussian point phase field viscous damping.

[0032] Furthermore, in step 3), the initialization of finite element node parameters includes the internal forces of the finite element nodes at the current moment. External forces at finite element nodes Finite element node phase field driving force Geometric resistance of phase field at finite element nodes ;

[0033] Initialize the background mesh node parameters, including the internal forces of the background mesh nodes at the current moment. External forces at background mesh nodes Background mesh node quality Momentum of background mesh nodes Background grid node phase field lumped coefficient matrix Background mesh node phase field driving force Geometric resistance of phase field with background mesh nodes .

[0034] Furthermore, in step 4), the strain update method for the Gaussian points inside the finite element is as follows:

[0035]

[0036] In the formula, Indicates the next moment, For the current moment, For time step, For the Gaussian point strain updated in the next moment, For the strain at the Gaussian point at the current moment, Let i be the total number of nodes in a single finite element, and j be the spatial dimensions. and These are the interpolation functions of the Gaussian point g at the current time at the finite element node K in the i and j directions, respectively. The derivative, and These are the velocities of the finite element node in the i-direction and j-direction at the current moment, respectively.

[0037] Gaussian point stress The result is obtained through constitutive equation updating, as shown below:

[0038]

[0039] In the formula, This is the constitutive relation function. For the updated deformation gradient, The deformation gradient rate of change It is an internal variable;

[0040] Phase field of Gaussian point inside the unit and phase field gradient As shown below:

[0041] ,

[0042] In the formula, Let K be the phase field at the current moment for a finite element node. The phase field value at the Gaussian point at the current moment, The gradient of the phase field at the Gaussian point at the current moment. Let g be the interpolation function of the Gaussian point g at the node K of the finite element at the current moment.

[0043] Furthermore, in step 5), the internal forces at the nodal points of the finite element... and external forces at finite element nodes As shown below:

[0044] ,

[0045] In the formula, The stress at the Gaussian point at the next moment. Let Jacobian matrix be the value of the Gaussian point at the current time. Gaussian point stamina;

[0046] Maximum tensile strain energy at the current moment at the Gaussian point As shown below:

[0047]

[0048]

[0049] In the formula, The tensile strain energy at the Gaussian point. and The Lamé coefficient is given. , and To adapt to changing circumstances ;

[0050] Phase field driving force of finite element nodes at the current moment Geometric resistance of phase field As shown below:

[0051]

[0052]

[0053] In the formula, It is a degenerate function; , The derivative of the degenerate function. The critical energy release rate of the Gaussian point crack. This represents the historical maximum tensile strain energy at the Gaussian point at the current moment.

[0054] Furthermore, in step 6), the finite element node boundary conditions refer to the conditions that satisfy the external forces at the finite element nodes at the current moment, for fixed boundary conditions. And the momentum of the nodal points of the finite element .

[0055] Furthermore, in step 7), the quality of the background mesh nodes at the current moment... Momentum of background mesh nodes As shown below:

[0056] , ,

[0057] In the formula, Let the momentum of the finite element node at the current moment be denoted as . Let K be the interpolation function of the finite element node K at the current moment on the background mesh node I, where I is the background mesh node number. This represents the total number of finite element nodes;

[0058] The internal forces at the nodes of the background mesh at the current moment and nodal external forces As shown below:

[0059] ,

[0060] Phase field coefficients of finite element nodes and the phase field driving force at the current moment. Geometric resistance of phase field Mapped onto the background mesh nodes, calculate the phase field coefficients of the background mesh nodes at the current time. Phase field driving force Geometric resistance of phase field As shown below:

[0061] , ,

[0062] In the formula, It is the interpolation function of the finite element node K at the current moment on the background mesh node I.

[0063] Furthermore, in step 8), the momentum increment of the background mesh nodes... and updated background mesh node momentum As shown below:

[0064] ,

[0065] In step 9), the velocity increment of the background mesh nodes and the updated background mesh node speed As shown below:

[0066] ,

[0067] In step 10), the phase field increment of the background mesh nodes. As shown below:

[0068]

[0069] In the formula, For time step.

[0070] Furthermore, in step 11), the speed of the updated background mesh nodes is... With speed increment Map back to finite element node and update the finite element node velocity. and location As shown below:

[0071] ,

[0072] In the formula, , These represent the velocity and position of the finite element node at the current moment. This represents the total number of nodes in a single background grid.

[0073] In step 12), the boundary conditions applied again to the finite element nodes are, for fixed boundaries: , ;

[0074] In step 13), update the phase field of the finite element nodes. As shown below:

[0075]

[0076] In the formula, Let K be the interpolation function of the finite element node K on the background mesh node I.

[0077] The technical effects of this invention are undeniable. This invention combines the finite element method and the material point method, supplemented by phase field fracture, and innovatively proposes a phase field-finite element-material point coupling algorithm for continuous body fracture simulation.

[0078] First, under the material point method calculation framework, this invention introduces the finite element concept, divides the continuum into finite elements, regards the element nodes as material points, and establishes a finite element-material point algorithm.

[0079] Secondly, compared with the traditional material point method, due to the existence of element nodes, the finite element material point method in this invention can effectively apply boundary conditions and avoid numerical fractures during the calculation process.

[0080] In addition, compared with traditional fracture mechanics, this invention proposes a phase-field method based on damage mechanics. This method describes strong discontinuous cracks as continuous diffuse cracks, with crack propagation driven by crack driving force. It does not require tracking the crack propagation path and can simulate complex crack propagation modes.

[0081] Finally, this invention adopts the USF (Update Stress First) time integration strategy, which updates the stress at the beginning of the time step when dealing with dynamic fracture problems. This allows the momentum equation and phase field evolution equation to be calculated at the same time level, which can effectively improve the accuracy of crack driving force calculation and improve the numerical stability of crack propagation simulation, making it more suitable for dynamic fracture simulation problems.

[0082] In summary, the method of this invention, by combining the finite element method, the material point method, and phase-field fracture theory, and based on the USF stress integral scheme, can accurately and effectively simulate the dynamic propagation process of continuum cracks. It provides a new method and approach for simulating continuum crack propagation, avoids the limitations of the finite element method in strongly nonlinear problems, and overcomes the numerical fracture problem in the material point method during the calculation process. It does not require continuous tracking of the crack propagation path, thus effectively improving the efficiency of crack propagation simulation. Attached Figure Description

[0083] Figure 1 The phase-field finite element material point method for mapping Gaussian points, finite element nodes, and background mesh nodes for the continuous dynamic fracture simulation of the present invention is described below.

[0084] Figure 2 The calculation process of the phase field finite element material point method for continuous dynamic fracture simulation of the present invention is as follows:

[0085] Figure 3 The model dimensions and boundary conditions for calculating a pre-cracked plate with one side under the velocity boundary conditions of this invention;

[0086] Figure 4This is the finite element discretization model for calculating a pre-cracked plate with one side under the velocity boundary conditions of the present invention.

[0087] Figure 5 The crack propagation at time 0.7 μs is calculated under the velocity boundary conditions of this invention;

[0088] Figure 6 The crack propagation at time 0.9 μs is calculated under the velocity boundary conditions of this invention.

[0089] Figure 7 The crack propagation at 1.5 μs is calculated under the velocity boundary conditions of this invention;

[0090] Figure 8 The model dimensions and boundary conditions for calculating a plate with a central precast crack under stress boundary conditions according to the present invention;

[0091] Figure 9 This is the finite element discretization model of the plate with a central precast crack under the stress boundary conditions of the present invention.

[0092] Figure 10 The crack propagation at 4.0 μs under the stress boundary conditions of this invention is shown.

[0093] Figure 11 The crack propagation at 6.0 μs under the stress boundary conditions of this invention is shown.

[0094] Figure 12 The crack propagation at 15.0 μs under the stress boundary conditions of this invention is shown.

[0095] Figure 13 The model dimensions and boundary conditions for calculating a pre-cracked plate with tilted prefabricated cracks under uniaxial compression according to the present invention;

[0096] Figure 14 This is the finite element discretization model for calculating a pre-cracked plate with tilted tilt under uniaxial compression according to the present invention.

[0097] Figure 15 This describes the crack propagation of the pre-fabricated crack plate at time 57.60 μs according to the present invention.

[0098] Figure 16 The crack propagation of the pre-fabricated crack plate at time 72.0 μs is shown in the present invention.

[0099] Figure 17 A comparison of the crack propagation of the precast crack plate at time 104.0 μs with the results of indoor tests. Detailed Implementation

[0100] The present invention will be further described below with reference to embodiments, but it should not be construed that the scope of the present invention is limited to the following embodiments. Various substitutions and modifications made based on ordinary technical knowledge and common practices in the art without departing from the above-described technical concept of the present invention should be included within the scope of protection of the present invention.

[0101] Example 1:

[0102] See Figures 1 to 17 A method for simulating continuum phase field damage and fracture combining finite element method and material point method includes the following steps:

[0103] Step 1) Construct the Eulerian background mesh and finite element mesh using the material point method, and initialize the physical material parameters and Gaussian point parameters within the elements of the finite element mesh;

[0104] Step 2) Based on the physical material density and phase field viscous damping of the finite element mesh, calculate the nodal mass and phase field coefficient matrix of the finite element nodes using the Gaussian point parameters inside the element.

[0105] Step 3) Initialize the finite element node parameters and background mesh node parameters;

[0106] Step 4) Based on the current velocity of the finite element node, update the stress and strain of the Gaussian points inside the finite element. Using the current phase field of the finite element node, calculate the phase field and phase field gradient of the Gaussian points inside the element.

[0107] Step 5) Based on the updated Gaussian point stress and strain, calculate the internal forces and external forces of the finite element nodes. Using the phase field and phase field gradient of the Gaussian point inside the element at the current moment, calculate the geometric resistance of the phase field of the finite element nodes and the driving force of the phase field of the finite element nodes.

[0108] Step 6) Apply boundary conditions to the finite element nodes;

[0109] Step 7) Map the finite element node information to the background mesh node, and calculate the background mesh node mass, background mesh node momentum, background mesh node internal force, background mesh node external force, background mesh node phase field coefficient, background mesh node phase field driving force, and background mesh node phase field geometric resistance at the current moment.

[0110] Step 8) Based on the internal forces and external forces of the background mesh nodes at the current moment, calculate the momentum increment of the background mesh nodes and update the momentum of the background mesh nodes;

[0111] Step 9) Based on the internal forces, external forces, and mass of the background mesh nodes at the current moment, calculate the velocity increment of the background mesh nodes and update the velocity of the background mesh nodes;

[0112] Step 10) Calculate the phase field increment of the background grid nodes based on the phase field coefficients, phase field geometric resistance, and phase field driving force of the background grid nodes at the current moment;

[0113] Step 11) Update the velocity and position of the finite element nodes based on the updated background mesh node velocities and the background mesh node velocity increments;

[0114] Step 12) Apply boundary conditions to the finite element nodes again;

[0115] Step 13) Update the phase field of the finite element nodes based on the phase field increment of the background mesh nodes;

[0116] Step 14) Return to Step 2), substitute the updated parameters and continue iteratively solving until the crack propagates, the continuum breaks, the internal stress approaches 0, the fracture energy is exhausted, and the energy is completely released.

[0117] Example 2:

[0118] The main structure of this embodiment is the same as that of embodiment 1. Further, in step 1), the steps of constructing the material point method Eulerian background mesh and the finite element mesh include:

[0119] Step 1.1) Establish a numerical model of the continuum to be simulated through 3D scanning;

[0120] Step 1.2) Measure the basic physical and mechanical parameters of the continuum;

[0121] The basic physical and mechanical parameters of the continuum include the elastic modulus E and density. Poisson's ratio and critical energy release rate ;

[0122] Step 1.3) Based on the numerical model of the continuum to be simulated and the basic physical and mechanical parameters of the continuum, construct the Eulerian background mesh and the finite element mesh using the material point method;

[0123] The physical material parameters of the finite element mesh include elastic modulus E and density. Compared to Poisson ;

[0124] The Gaussian point parameters within the unit include position. Gaussian point stress Gaussian point strain Gaussian point phase field Critical energy release rate of material cracks Phase field viscous damping Phase field characteristic width Maximum crack driving force in historical phase field .

[0125] Example 3:

[0126] The main structure of this embodiment is the same as any one of embodiments 1-2. Further, in step 2), the mass of the finite element nodes... As shown below:

[0127]

[0128] In the formula, e is the finite element number, g is the Gaussian point number inside the element, and K is the node number of the finite element. For the density of a continuous material, The total number of Gaussian points within a finite cell. This is the interpolation of the Gaussian point g at the node K of the finite element. The weights of the Gaussian points, Let the Jacobian matrix be a Gaussian point;

[0129] Finite element node phase field coefficient matrix As shown below:

[0130]

[0131] In the formula, It is a Gaussian point phase field viscous damping.

[0132] Example 4:

[0133] The main structure of this embodiment is the same as any one of embodiments 1 to 3. Further, in step 3), the finite element node parameters are initialized, including the internal forces of the finite element nodes at the current moment. External forces at finite element nodes Finite element node phase field driving force Geometric resistance of phase field at finite element nodes ;

[0134] Initialize the background mesh node parameters, including the internal forces of the background mesh nodes at the current moment. External forces at background mesh nodes Background mesh node quality Momentum of background mesh nodes Background grid node phase field lumped coefficient matrix Background mesh node phase field driving force Geometric resistance of phase field with background mesh nodes .

[0135] Example 5:

[0136] The main structure of this embodiment is the same as any one of embodiments 1-4. Further, in step 4), the strain update method of the Gaussian point inside the finite element is as follows:

[0137]

[0138] In the formula, Indicates the next moment, For the current moment, For time step, For the Gaussian point strain updated in the next moment, For the strain at the Gaussian point at the current moment, Let i be the total number of nodes in a single finite element, and j be the spatial dimensions. and These are the interpolation functions of the Gaussian point g at the current time at the finite element node K in the i and j directions, respectively. The derivative, and These are the velocities of the finite element node in the i-direction and j-direction at the current moment, respectively.

[0139] Gaussian point stress The result is obtained through constitutive equation updating, as shown below:

[0140]

[0141] In the formula, This is the constitutive relation function. For the updated deformation gradient, The deformation gradient rate of change It is an internal variable;

[0142] Phase field of Gaussian point inside the unit and phase field gradient As shown below:

[0143] ,

[0144] In the formula, Let K be the phase field at the current moment for a finite element node. The phase field value at the Gaussian point at the current moment, The gradient of the phase field at the Gaussian point at the current moment. Let g be the interpolation function of the Gaussian point g at the node K of the finite element at the current moment.

[0145] Example 6:

[0146] The main structure of this embodiment is the same as any one of embodiments 1-5. Further, in step 5), the internal forces at the finite element nodes... and external forces at finite element nodes As shown below:

[0147] ,

[0148] In the formula, The stress at the Gaussian point at the next moment. Let Jacobian matrix be the value of the Gaussian point at the current time. Gaussian point stamina;

[0149] Maximum tensile strain energy at the current moment at the Gaussian point As shown below:

[0150]

[0151]

[0152] In the formula, The tensile strain energy at the Gaussian point. and The Lamé coefficient is given. , and To adapt to changing circumstances ;

[0153] Phase field driving force of finite element nodes at the current moment Geometric resistance of phase field As shown below:

[0154]

[0155]

[0156] In the formula, It is a degenerate function; , The derivative of the degenerate function. The critical energy release rate of the Gaussian point crack. This represents the historical maximum tensile strain energy at the Gaussian point at the current moment.

[0157] Example 7:

[0158] The main structure of this embodiment is the same as any one of embodiments 1-6. Further, in step 6), the finite element node boundary condition refers to, for fixed boundary conditions, satisfying the external force of the finite element node at the current moment. And the momentum of the nodal points of the finite element .

[0159] Example 8:

[0160] The main structure of this embodiment is the same as any one of embodiments 1-7. Further, in step 7), the quality of the background mesh nodes at the current moment... Momentum of background mesh nodes As shown below:

[0161] , ,

[0162] In the formula, Let the momentum of the finite element node at the current moment be denoted as . Let K be the interpolation function of the finite element node K at the current moment on the background mesh node I, where I is the background mesh node number. This represents the total number of finite element nodes;

[0163] The internal forces at the nodes of the background mesh at the current moment and nodal external forces As shown below:

[0164] ,

[0165] Phase field coefficients of finite element nodes and the phase field driving force at the current moment. Geometric resistance of phase field Mapped onto the background mesh nodes, calculate the phase field coefficients of the background mesh nodes at the current time. Phase field driving force Geometric resistance of phase field As shown below:

[0166] , ,

[0167] In the formula, It is the interpolation function of the finite element node K at the current moment on the background mesh node I.

[0168] Example 9:

[0169] The main structure of this embodiment is the same as any one of embodiments 1-8. Further, in step 8), the momentum increment of the background mesh nodes... and updated background mesh node momentum As shown below:

[0170] ,

[0171] In step 9), the velocity increment of the background mesh nodes and the updated background mesh node speed As shown below:

[0172] ,

[0173] In step 10), the phase field increment of the background mesh nodes. As shown below:

[0174]

[0175] In the formula, For time step.

[0176] Example 10:

[0177] The main structure of this embodiment is the same as any one of embodiments 1 to 9. Further, in step 11), the speed of the updated background mesh nodes is... With speed increment Map back to finite element node and update the finite element node velocity. and location As shown below:

[0178] ,

[0179] In the formula, , These represent the velocity and position of the finite element node at the current moment. This represents the total number of nodes in a single background grid.

[0180] In step 12), the boundary conditions applied again to the finite element nodes are, for fixed boundaries: , ;

[0181] In step 13), update the phase field of the finite element nodes. As shown below:

[0182]

[0183] In the formula, for the updated The scope should be If it appears For cases less than 0, take =0; if it appears For cases greater than 1, take =1.

[0184] Example 11:

[0185] The present invention aims to provide a continuum damage and fracture simulation algorithm that combines the finite element method, the material point method and the phase field fracture method, in order to solve the defects of the traditional material point method in the application of boundary conditions, numerical fracture, Gaussian point migration, crack path tracking during fracture simulation, low computational efficiency and inability to simulate complex crack propagation problems.

[0186] To achieve the above objectives, this invention provides a method for simulating continuum phase field damage and fracture by combining finite element method and material point method, specifically including the following steps:

[0187] Step 101: Construct the Eulerian background mesh and finite element mesh using the material point method. Initialize the physical material parameters of the finite element mesh and the Gaussian point parameters within the elements. The physical material parameters of the finite element mesh include the elastic modulus E and density. Compared to Poisson The Gaussian point parameters inside the element include position. Gaussian point stress Gaussian point strain Gaussian point phase field Critical energy release rate of material cracks Phase field viscous damping Phase field characteristic width Maximum crack driving force in historical phase field It should be noted that bold letters represent tensors.

[0188] Step 102: Calculate the mass of the finite element nodes Phase field coefficient matrix of nodes

[0189] In step 102, the mass of the finite element node is,

[0190] (1)

[0191] In the formula, e is the finite element number, g is the Gaussian point number inside the element, and K is the node number of the finite element. For the density of a continuous material, The total number of Gaussian points within a finite cell. This is the interpolation of the Gaussian point g at the node K of the finite element. The weights of the Gaussian points, Let the Jacobian matrix be a Gaussian point;

[0192] In step 102, the phase field coefficient matrix of the finite element nodes for,

[0193] (2)

[0194] In the formula, It is a Gaussian point phase field viscous damping.

[0195] In step 102, the finite element node mass and the node phase field coefficient matrix remain unchanged throughout the calculation process.

[0196] Step 103: Initialize the finite element nodal parameters, including the finite element nodal internal forces at the current moment. External forces at finite element nodes Finite element node phase field driving force Geometric resistance of phase field at finite element nodes Initialize the background mesh node parameters, including the internal forces of the background mesh nodes at the current moment. External forces at background mesh nodes Background mesh node quality Momentum of background mesh nodes Background grid node phase field lumped coefficient matrix Background mesh node phase field driving force Geometric resistance of phase field with background mesh nodes .

[0197] Step 104: Based on the current velocity of the finite element node, update the stress and strain of the Gaussian point inside the finite element; and use the current phase field of the finite element node to calculate the phase field and phase field gradient of the Gaussian point inside the element.

[0198] In step 104, the strain at the Gaussian point is updated as follows:

[0199] (3)

[0200] In the formula, Indicates the next moment. For the present moment, For time step, For the Gaussian point strain updated in the next moment, For the strain at the Gaussian point at the current moment, Let i be the total number of nodes in a single finite element, and j be the spatial dimensions. and These are the interpolation functions of the Gaussian point g at the current time at the finite element node K in the i and j directions, respectively. The derivative, and Let be the velocities of the finite element node in the i and j directions at the current moment, respectively.

[0201] In step 104, the Gaussian point stress is updated using the constitutive equation. .

[0202] (4)

[0203] In the formula, This is the constitutive relation function. For the updated deformation gradient, The deformation gradient rate of change It is an internal variable.

[0204] In step 104, the Gaussian point phase field and phase field gradient for,

[0205] , (5)

[0206] In the formula, Let K be the phase field at the current moment for a finite element node. The phase field value at the Gaussian point at the current moment, The gradient of the phase field at the Gaussian point at the current moment. Let g be the interpolation function of the Gaussian point g at the node K of the finite element at the current moment.

[0207] Step 105: Calculate the internal and external forces at the finite element nodes using the updated Gaussian point stress and strain, and calculate the geometric resistance and driving force of the phase field at the finite element nodes using the phase field and phase field gradient at the current moment.

[0208] In step 105, the internal forces at the nodal points of the finite element... and nodal external forces ,

[0209] , (6)

[0210] In the formula, The stress at the Gaussian point at the next moment. Let Jacobian matrix be the value of the Gaussian point at the current time. Gaussian point stamina;

[0211] In step 105, the historical maximum tensile strain energy at the Gaussian point at the current moment is calculated. ,

[0212] (7)

[0213] (8)

[0214] In the formula, The tensile strain energy at the Gaussian point. and The Lamé coefficient is given. , and To adapt to changing circumstances .

[0215] In step 105, the phase field driving force of the finite element node at the current moment. Geometric resistance of phase field for,

[0216] (9)

[0217] (10)

[0218] In the formula, It is a degenerate function; , The derivative of the degenerate function. The critical energy release rate of the Gaussian point crack. This represents the historical maximum tensile strain energy at the Gaussian point at the current moment.

[0219] Step 106: Apply boundary conditions: For fixed boundaries, satisfy the external forces at the finite element nodes at the current moment. And the momentum of the nodal points of the finite element .

[0220] Step 107: Map the current finite element node information to the background mesh node and update the background mesh node parameters.

[0221] In step 107, the finite element node mass is... and momentum Mapped onto background mesh nodes, calculate the quality of the background mesh nodes at the current time. and momentum .

[0222] , , (11)

[0223] In the formula, Let the momentum of the finite element node at the current moment be denoted as . Let K be the interpolation function of the finite element node K at the current moment on the background mesh node I, where I is the background mesh node number. This represents the total number of finite element nodes.

[0224] In step 107, the internal forces at the nodes of the finite element are... and external forces Map the values ​​onto the background mesh nodes and calculate the internal forces of the background mesh nodes at the current moment. and external forces .

[0225] , (12)

[0226] In step 107, the phase field coefficients of the finite element nodes are... and the phase field driving force at the current moment. Geometric resistance of phase field Mapped onto the background mesh nodes, calculate the phase field coefficients of the background mesh nodes at the current time. Phase field driving force Geometric resistance of phase field

[0227] , , (13)

[0228] Step 108: Calculate the momentum increment of the background mesh nodes Update the momentum of background mesh nodes. .

[0229] , (14)

[0230] Step 109: Calculate the velocity increment of background mesh nodes And the updated background mesh node speed .

[0231] , (15)

[0232] Step 110: Calculate the phase field increment at the background mesh nodes .

[0233] (16)

[0234] Step 111: Update the speed of the background mesh nodes With speed increment Map back to finite element node and update the finite element node velocity. and location .

[0235] , (17)

[0236] In the formula, This represents the total number of nodes in a single background grid.

[0237] Step 112: Apply boundary conditions again to the updated finite element nodes. For fixed boundaries: , .

[0238] Step 113: Increment the phase field of the background mesh nodes Map back to finite element nodes and update the phase field of the finite element nodes. .

[0239] (18)

[0240] In the formula, for the updated The scope should be If it appears For cases less than 0, take =0; if it appears For cases greater than 1, take =1.

[0241] Step 114: Return to step 103 and continue iteratively solving until the continuous crack propagation fails, the continuous breaks, the internal stress approaches 0, the fracture energy is exhausted, and the energy is completely released.

[0242] Example 12:

[0243] This embodiment discloses a continuum phase field damage and fracture simulation method combining finite element method and material point method, based on its mapping scheme as follows: Figure 1 And the overall calculation process, such as Figure 2 For simulating the fracture process of engineering-scale or laboratory specimens.

[0244] Figure 1 In this code, 1 represents finite element discretization, 2 represents calculating the phase field value, velocity, and position of the Gaussian point, 3 represents calculating the internal forces, phase field geometric resistance, and driving force of the finite element nodes, 4 represents mapping from finite element nodes to background mesh nodes, 5 represents solving for the velocity, velocity increment, and phase field increment of the background mesh nodes, 6 represents mapping from background mesh nodes to finite element nodes, and 7 represents updating the velocity, position, and phase field of the finite element nodes.

[0245] The simulation steps are as follows:

[0246] 1) Based on on-site engineering or laboratory specimens, establish a numerical continuum for calculation using methods such as 3D scanning, for example... Figures 3-4 , Figures 8-9 and Figures 13-14 .

[0247] 2) Utilize indoor testing methods to measure the fundamental physical and mechanical parameters of the continuum, including the elastic modulus E and density. Poisson's ratio and critical energy release rate .

[0248] 3) Construct the Eulerian background mesh and finite element mesh using the material point method, and initialize the physical material parameters and Gaussian point parameters within the elements of the finite element mesh. The physical material parameters of the finite element mesh include the elastic modulus E and density. Compared to Poisson The Gaussian point parameters inside the element include position. Gaussian point stress Gaussian point strain Gaussian point phase field Critical energy release rate of material cracks Phase field viscous damping Phase field characteristic width Maximum crack driving force in historical phase field It should be noted that bold letters represent tensors.

[0249] 4) Calculate the mass of finite element nodes. Phase field coefficient matrix of nodes It should be noted that the nodal particles and nodal phase field coefficient matrices of the finite element remain unchanged throughout the entire calculation process.

[0250] , (1)

[0251] In the formula, e is the finite element number, g is the Gaussian point number inside the element, and K is the node number of the finite element. For the density of a continuous material, The total number of Gaussian points within a finite cell. This is the interpolation of the Gaussian point g at the node K of the finite element. The weights of the Gaussian points, Let be the Jacobian matrix of the Gaussian point. It is a Gaussian point phase field viscous damping.

[0252] 5) Initialize the finite element nodal parameters, including the finite element nodal internal forces at the current moment. External forces at finite element nodes Finite element node phase field driving force Geometric resistance of phase field at finite element nodes Initialize the background mesh node parameters, including the internal forces of the background mesh nodes at the current moment. External forces at background mesh nodes Background mesh node quality Momentum of background mesh nodes Background grid node phase field lumped coefficient matrix Background mesh node phase field driving force Geometric resistance of phase field with background mesh nodes .

[0253] 6) Use the velocity of the finite element node at the current moment to update the stress and strain of the Gaussian point inside the finite element, and use the phase field of the finite element node at the current moment to calculate the phase field and phase field gradient of the Gaussian point inside the element.

[0254] 6-1) First update the strain at the Gaussian point. Then, the Gaussian point stress is updated using the constitutive equation. .

[0255] (2)

[0256] (3)

[0257] In the formula, Indicates the next moment. For the present moment, For time step, For the Gaussian point strain updated in the next moment, For the strain at the Gaussian point at the current moment, Let i be the total number of nodes in a single finite element, and j be the spatial dimensions. and These are the interpolation functions of the Gaussian point g at the current time at the finite element node K in the i and j directions, respectively. The derivative, and Let be the velocities of the finite element node in the i-direction and j-direction at the current moment, respectively. This is the constitutive relation function. For the updated deformation gradient, The deformation gradient rate of change It is an internal variable.

[0258] 6-2) Phase field at Gauss point and phase field gradient for:

[0259] , (4)

[0260] In the formula, Let K be the phase field at the current moment for a finite element node. The phase field value at the Gaussian point at the current moment, The gradient of the phase field at the Gaussian point at the current moment. Let g be the interpolation function of the Gaussian point g at the node K of the finite element at the current moment.

[0261] 7) Using the updated Gaussian point stress and strain, calculate the internal and external forces at the finite element nodes. Using the phase field and phase field gradient at the current moment, calculate the phase field geometric resistance and phase field driving force at the finite element nodes.

[0262] 7-1) Using the updated Gaussian point stress and strain, calculate the nodal forces of the finite element at the current moment. and nodal external forces .

[0263] , (5)

[0264] In the formula, The stress at the Gaussian point at the next moment. Let Jacobian matrix be the value of the Gaussian point at the current time. Let the Gaussian point have the physical strength.

[0265] 7-2) Calculate the historical maximum tensile strain energy at the current moment of the Gaussian point. .

[0266] (6)

[0267] (7)

[0268] In the formula, The tensile strain energy at the Gaussian point. and The Lamé coefficient is given. , and To adapt to changing circumstances .

[0269] 7-3) Calculate the phase driving force of the finite element node at the current moment using the historical maximum tensile strain energy, Gaussian point phase field, and phase field gradient. Geometric resistance of phase field for:

[0270] (8)

[0271] (9)

[0272] In the formula, It is a degenerate function; , The derivative of the degenerate function. The critical energy release rate of the Gaussian point crack. This represents the historical maximum tensile strain energy at the Gaussian point at the current moment.

[0273] 8) Apply boundary conditions: For fixed boundaries, satisfy the external forces at the finite element nodes at the current moment. And the momentum of the nodal points of the finite element .

[0274] 9) Map the current finite element node information to the background mesh node and update the background mesh node parameters.

[0275] 9-1) The mass of finite element nodes and momentum Mapped onto background mesh nodes, calculate the quality of the background mesh nodes at the current time. and momentum .

[0276] , , (10)

[0277] In the formula, Let the momentum of the finite element node at the current moment be denoted as . Let K be the interpolation function of the finite element node K at the current moment on the background mesh node I, where I is the background mesh node number. This represents the total number of finite element nodes.

[0278] 9-2) Internal forces at the nodes of a finite element and external forces Map the values ​​onto the background mesh nodes and calculate the internal forces of the background mesh nodes at the current moment. and external forces .

[0279] , (11)

[0280] 9-3) The phase field coefficients of the finite element nodes and the phase field driving force at the current moment. Geometric resistance of phase field Mapped onto the background mesh nodes, calculate the phase field coefficients of the background mesh nodes at the current time. Phase field driving force Geometric resistance of phase field .

[0281] , , (12)

[0282] 10) Calculate the momentum increment of the background mesh nodes using the internal and external forces at the background mesh nodes. Update the momentum of background mesh nodes. .

[0283] , (13)

[0284] 11) Calculate the velocity increment of background mesh nodes. And the updated background mesh node speed .

[0285] , (14)

[0286] 12) Calculate the phase field increment at the background mesh nodes using the phase field geometric resistance and phase field crack driving force at the background mesh nodes. .

[0287] (15)

[0288] 13) Update the speed of the background mesh nodes. With speed increment Map back to finite element node and update the finite element node velocity. and location .

[0289] , (16)

[0290] In the formula, This represents the total number of nodes in a single background grid.

[0291] 14) Apply boundary conditions again to the updated finite element nodes. For fixed boundaries: , .

[0292] 15) Increase the phase field of the background mesh nodes. Map back to finite element nodes and update the phase field of the finite element nodes. .

[0293] (17)

[0294] In the formula, for the updated The scope should be If it appears For cases less than 0, take =0; if it appears For cases greater than 1, take =1.

[0295] Return to step (5) and continue iteratively solving until the continuum crack propagates, the continuum breaks, the internal stress approaches 0, the fracture energy is exhausted, and the energy is completely released, as shown below. Figures 5-7 , Figures 10-12 and Figure 15-17 .

Claims

1. A method for simulating damage and fracture in a continuum phase field combining finite element method and material point method, characterized in that, Includes the following steps: Step 1) Construct the Eulerian background mesh and finite element mesh using the material point method, and initialize the physical material parameters and Gaussian point parameters within the elements of the finite element mesh; Step 2) Based on the physical material density and phase field viscous damping of the finite element mesh, calculate the nodal mass and phase field coefficient matrix of the finite element nodes using the Gaussian point parameters inside the element. Step 3) Initialize the finite element node parameters and background mesh node parameters; Step 4) Based on the current velocity of the finite element node, update the stress and strain of the Gaussian points inside the finite element. Using the current phase field of the finite element node, calculate the phase field and phase field gradient of the Gaussian points inside the element. Step 5) Based on the updated Gaussian point stress and strain, calculate the internal forces and external forces of the finite element nodes. Using the phase field and phase field gradient of the Gaussian point inside the element at the current moment, calculate the geometric resistance of the phase field of the finite element nodes and the driving force of the phase field of the finite element nodes. Step 6) Apply boundary conditions to the finite element nodes; Step 7) Map the finite element node information to the background mesh node, and calculate the background mesh node mass, background mesh node momentum, background mesh node internal force, background mesh node external force, background mesh node phase field coefficient, background mesh node phase field driving force, and background mesh node phase field geometric resistance at the current moment. Step 8) Based on the internal forces and external forces of the background mesh nodes at the current moment, calculate the momentum increment of the background mesh nodes and update the momentum of the background mesh nodes; Step 9) Based on the internal forces, external forces, and mass of the background mesh nodes at the current moment, calculate the velocity increment of the background mesh nodes and update the velocity of the background mesh nodes; Step 10) Calculate the phase field increment of the background grid nodes based on the phase field coefficients, phase field geometric resistance, and phase field driving force of the background grid nodes at the current moment; Step 11) Update the velocity and position of the finite element nodes based on the updated background mesh node velocities and the background mesh node velocity increments; Step 12) Apply boundary conditions to the finite element nodes again; Step 13) Update the phase field of the finite element nodes based on the phase field increment of the background mesh nodes; Step 14) Return to Step 2), substitute the updated parameters and continue iteratively solving until the crack propagates, the continuum breaks, the internal stress approaches 0, the fracture energy is exhausted, and the energy is completely released.

2. The method for simulating continuum phase field damage and fracture combining finite element method and material point method according to claim 1, characterized in that, Step 1) involves constructing the material point method Eulerian background mesh and the finite element mesh, including: Step 1.1) Establish a numerical model of the continuum to be simulated through 3D scanning; Step 1.2) Measure the basic physical and mechanical parameters of the continuum; The basic physical and mechanical parameters of the continuum include the elastic modulus E and density. Poisson's ratio and critical energy release rate ; Step 1.3) Based on the numerical model of the continuum to be simulated and the basic physical and mechanical parameters of the continuum, construct the Eulerian background mesh and the finite element mesh using the material point method; The physical material parameters of the finite element mesh include elastic modulus E and density. Compared to Poisson ; The Gaussian point parameters within the unit include position. Gaussian point stress Gaussian point strain Gaussian point phase field Critical energy release rate of material cracks Phase field viscous damping Phase field characteristic width Maximum crack driving force in historical phase field .

3. The method for simulating continuum phase field damage and fracture combining finite element method and material point method according to claim 2, characterized in that, In step 2), the finite element node mass As shown below: In the formula, e is the finite element number, g is the Gaussian point number inside the element, and K is the node number of the finite element. For the density of a continuous material, The total number of Gaussian points within a finite cell. This is the interpolation of the Gaussian point g at the node K of the finite element. The weights of the Gaussian points, Let the Jacobian matrix be a Gaussian point; Finite element node phase field coefficient matrix As shown below: In the formula, It is a Gaussian point phase field viscous damping.

4. The method for simulating continuum phase field damage and fracture combining finite element method and material point method according to claim 1, characterized in that, In step 3), the initialization of finite element nodal parameters includes the internal forces of the finite element nodules at the current moment. External forces at finite element nodes Finite element node phase field driving force Geometric resistance of phase field at finite element nodes ; Initialize the background mesh node parameters, including the internal forces of the background mesh nodes at the current moment. External forces at background mesh nodes Background mesh node quality Momentum of background mesh nodes Background grid node phase field lumped coefficient matrix Background mesh node phase field driving force Geometric resistance of phase field with background mesh nodes .

5. The method for simulating continuum phase field damage and fracture combining finite element method and material point method according to claim 3, characterized in that, In step 4), the strain update method for Gaussian points within the finite element is as follows: In the formula, Indicates the next moment, For the current moment, For time step, For the Gaussian point strain updated in the next moment, For the strain at the Gaussian point at the current moment, Let i be the total number of nodes in a single finite element, and j be the spatial dimensions. and These are the interpolation functions of the Gaussian point g at the current time at the finite element node K in the i and j directions, respectively. The derivative, and These are the velocities of the finite element node in the i-direction and j-direction at the current moment, respectively. Gaussian point stress The result is obtained through constitutive equation updating, as shown below: In the formula, This is the constitutive relation function. For the updated deformation gradient, The deformation gradient rate of change It is an internal variable; Phase field of Gaussian point inside the unit and phase field gradient As shown below: ， In the formula, Let K be the phase field at the current moment for a finite element node. The phase field value at the Gaussian point at the current moment, The gradient of the phase field at the Gaussian point at the current moment. Let g be the interpolation function of the Gaussian point g at the node K of the finite element at the current moment.

6. The method for simulating continuum phase field damage and fracture combining finite element method and material point method according to claim 5, characterized in that, In step 5), the internal forces at the nodal points of the finite element and external forces at finite element nodes As shown below: ， In the formula, The stress at the Gaussian point at the next moment. Let Jacobian matrix be the value of the Gaussian point at the current time. Gaussian point stamina; Maximum tensile strain energy at the current moment at the Gaussian point As shown below: In the formula, The tensile strain energy at the Gaussian point. and The Lamé coefficient is given. , and To adapt to changing circumstances ; Phase field driving force of finite element nodes at the current moment Geometric resistance of phase field As shown below: In the formula, It is a degenerate function; , The derivative of the degenerate function. The critical energy release rate of the Gaussian point crack. This represents the historical maximum tensile strain energy at the Gaussian point at the current moment.

7. The method for simulating continuum phase field damage and fracture combining finite element method and material point method according to claim 4, characterized in that, In step 6), the finite element node boundary conditions refer to the conditions that satisfy the external forces at the finite element nodes at the current moment, for fixed boundary conditions. And the momentum of the nodal points of the finite element .

8. The method for simulating continuum phase field damage and fracture combining finite element method and material point method according to claim 6, characterized in that, In step 7), the quality of the background mesh nodes at the current moment. Momentum of background mesh nodes As shown below: ， ， In the formula, Let the momentum of the finite element node at the current moment be denoted as . Let K be the interpolation function of the finite element node K at the current moment on the background mesh node I, where I is the background mesh node number. This represents the total number of finite element nodes; The internal forces at the nodes of the background mesh at the current moment and nodal external forces As shown below: ， Phase field coefficients of finite element nodes and the phase field driving force at the current moment. Geometric resistance of phase field Mapped onto the background mesh nodes, calculate the phase field coefficients of the background mesh nodes at the current time. Phase field driving force Geometric resistance of phase field As shown below: ， ， In the formula, It is the interpolation function of the finite element node K at the current moment on the background mesh node I.

9. The method for simulating continuum phase field damage and fracture combining finite element method and material point method according to claim 8, characterized in that, In step 8), the momentum increment of the background mesh nodes and updated background mesh node momentum As shown below: ， In step 9), the velocity increment of the background mesh nodes and the updated background mesh node speed As shown below: ， In step 10), the phase field increment of the background mesh nodes. As shown below: In the formula, For time step.

10. A method for simulating continuum phase field damage and fracture combining finite element method and material point method according to claim 9, characterized in that, In step 11), the speed of the updated background mesh nodes is... With speed increment Map back to finite element node and update the finite element node velocity. and location As shown below: ， In the formula, , These represent the velocity and position of the finite element node at the current moment. This represents the total number of nodes in a single background grid. In step 12), the boundary conditions applied again to the finite element nodes are, for fixed boundaries: , ; In step 13), update the phase field of the finite element nodes. As shown below: In the formula, Let K be the interpolation function of the finite element node K on the background mesh node I.

Images