Dynamic impact / contact toughness metal fracture analysis explicit cohesive force phase field material point method

Abstract

The present application belongs to the technical field of ductile metal fracture analysis, and proposes an explicit cohesive phase field material point method, which provides a new numerical calculation method for the study of metal dynamic fracture failure. In the method, for the large deformation problem and high-speed impact problem in the ductile metal layer cracking and penetration failure, the cohesive phase field material point method and the elastic ball compliant contact algorithm are innovatively proposed. The cohesive phase field model can well describe the damage softening and fracture behavior of ductile metal. The elastic ball compliant contact algorithm effectively suppresses the velocity over-shock of the impact contact surface and the difficulty in calculating the outer normal vector of the contact surface. Then, the control equation of the ductile metal dynamic fracture is derived by using the Lagrange equation, and the control equation of the coupling displacement and cohesive phase field model is discretized in the material point method framework. Finally, the phase field-displacement field staggered solution strategy is adopted to stably and efficiently solve the strong nonlinear large-scale ductile fracture failure problem such as contact and large deformation.

Description

Technical Field

[0001] This invention belongs to the field of metal ductile fracture analysis technology, and relates to an explicit cohesive phase field material point method for dynamic impact / contact ductile fracture analysis. Background Technology

[0002] Tough metals possess excellent plastic deformation capabilities, effectively handling complex stress and strain conditions, and are therefore widely used in military, automotive, construction, and aerospace fields. However, under extreme loads, these metals can also fail, such as through high-speed bullet penetration and delamination caused by dynamic tension under impact loading. This failure process involves multi-scale analysis across spatial and temporal dimensions, representing a complex interdisciplinary challenge at the intersection of physics, mechanics, and materials science. While many researchers have conducted extensive experiments, the high cost and limitations of measurement techniques have led to the adoption of numerical simulation as a more efficient and economical method to study complex failure mechanisms, particularly for problems involving high-speed impact and extreme deformation. Accurate numerical simulation of high-speed impact and extreme deformation processes depends on material constitutive models, contact algorithms, damage models, and numerical methods. Because these factors are influenced by complex nonlinear interactions between cracks, material behavior, geometry, and contact conditions, accurate and reliable numerical simulation of high-speed impact and extreme deformation processes remains a significant challenge.

[0003] Crack modeling methods can be broadly categorized into two types: discrete crack models and diffuse crack models. Discrete crack models represent cracks as distinct entities, such as the removed element method, virtual crack closure technique (VCCT), extended finite element method (XFEM), and crack-driven configuration force method. However, these methods typically require additional degrees of freedom (as in XFEM) and complex algorithms to track crack branching, merging, and propagation. In contrast, diffuse crack models, such as peri-field dynamics (PD), phase-field fracture models, and crack particle method (CPM), represent cracks using continuous field variables, thus eliminating the need for explicit crack tracking algorithms. Among these, the phase-field fracture model has shown great potential in solving complex fracture problems. This model has been successfully applied to various fracture problems, such as brittle fracture, quasi-static and dynamic fracture, fatigue analysis, and fracture in porous media. Although phase-field fracture analysis methods for brittle fracture are very mature, extending them to ductile fracture, especially ductile fracture under finite deformation conditions, remains an active research area. For example, Hofacker, Miehe, and Ulmer et al., in "Aphase field model for ductile to brittle failure mode transition," extended phase-field modeling to ductile fracture by using elastic and plastic energies as crack driving forces. Building on this, Wu, in "Aunified phase-field theory for the mechanics of damage and quasi-brittle failure," established a unified phase-field theory applicable to both brittle fracture and quasi-brittle damage, and subsequently proposed the cohesive phase-field model (PF-CZM), which has shown significant advantages in predicting complex crack propagation in ductile metals. However, most of the aforementioned work relies on the finite element method (FEM) to solve the governing equations. While the finite element method is one of the most commonly used numerical methods in mechanical simulations, its application in large deformation problems is hindered by mesh deformation issues, which greatly increases computational complexity. Furthermore, this complexity increases further when contact problems are incorporated into the finite element framework.

[0004] The Material Point Method (MPM) is a particle-based computational framework that combines Eulerian and Lagrangian descriptions, and has become a powerful tool for solving mesh deformation problems while maintaining computational efficiency in large deformation analysis. Originally proposed by Sulsky et al. in "Aparticle method for history-dependent materials," this method has achieved great success in various fields of finite deformation mechanics. For example, Hu et al. proposed an explicit phase-field coupled material point method (ePF-CCPDI) in "Coupling explicit phase-field MPM for two-dimensional hydromechanical fracture in poro-elastoplastic media" to simulate two-dimensional elastoplastic fracture in porous media. Zhang et al. proposed an explicit phase-field total Lagrangian material point method (TLMPM) in "Explicit phase-field total Lagrangian material point method for the dynamic fracture of hyperelastic materials," successfully simulating the large deformation dynamic fracture behavior of hyperelastic solids. Sun et al. proposed a phase-field material point method in "Explicit phase-field material point method for thermally induced fractures" to simulate the large deformation fracture of circular thin ceramic specimens under thermal shock. Despite these advances, contact algorithms within the material point method framework remain a significant challenge, particularly in impact simulations. Contact algorithms are crucial for accurately capturing the interactions between objects undergoing large deformations. However, existing material point contact algorithms still have limitations that affect the accuracy and stability of numerical simulations: (1) Mesh-based contact algorithms face numerical difficulties in determining the normal vectors of the contact surfaces, especially for complex geometries undergoing severe deformation. (2) While particle-based contact algorithms do not require the calculation of normal vectors, the limitations of non-penetrating constraints can lead to over-oscillation of contact forces, thereby disrupting stress wave propagation. (3) Frictional contact is highly sensitive to the discrete size and manner of material points at geometric boundaries, making convergence difficult in frictional problems. These issues highlight the urgent need to develop novel material point contact algorithms.

[0005] Given the limitations of existing contact algorithms, a deeper understanding of contact behavior is needed. Collision behavior is typical of nonlinear and discontinuous nonsmooth systems, characterized by sudden changes in system velocity caused by collisions. Current methods describing the contact behavior of colliding bodies can be broadly categorized into two types: nonsmooth dynamics methods and regularized methods. Nonsmooth dynamics methods assume that the contacting bodies remain unchanged throughout the collision. In contrast, regularized methods allow deformation within the contact region and define the contact force as a function of the deformation. When dealing with contact problems considering friction, nonsmooth dynamics methods may lead to multiple solutions or no solutions, as well as energy non-conservation. In contrast, regularized methods, also known as compliant contact force models (CCFM), avoid many problems encountered in rigid body methods. This method is widely used in multibody dynamics software, including ADAMS, RecurDyn, and EDEM, due to its excellent performance. Therefore, combining compliant contact force models with particle-based contact methods can effectively address the shortcomings of existing material point contact algorithms.

[0006] This invention proposes a novel coupled explicit cohesive phase-field material point method. To predict crack propagation and material failure behavior of ductile metals under dynamic loads, this paper establishes a crack model using the cohesive phase-field method and derives the coupled control equations from an energy perspective using the Lagrangian method. Then, within the framework of the material point method, the coupled equations of the displacement field and the cohesive phase field are solved to address the fracture damage behavior of ductile metals. Simultaneously, a novel material point contact algorithm based on a compliant contact force model, the Ball-Ball Compliant Contact Algorithm (PCCA), is proposed. Compared with traditional contact algorithms, this algorithm exhibits higher stability when solving friction and high-speed impact problems. Finally, the algorithm employs an alternating scheme and an explicit solution strategy to achieve robust and efficient solutions to the coupled field control equations. Summary of the Invention

[0007] To address the efficient numerical analysis of dynamic impact / contact ductile fracture in metals, this invention innovatively proposes an explicit cohesive phase-field material point method (PF-CZM-MPM-PCCA) based on a novel pinball compliant contact algorithm for dynamic impact / contact ductile fracture analysis. The objectives are twofold: First, to overcome the limitations of traditional damage models in accurately capturing crack propagation paths, and to address the inefficiencies of implicit phase-field analysis in describing material softening behavior during ductile fracture and large-scale crack propagation, this invention derives an explicit cohesive phase-field fracture model suitable for ductile fracture analysis based on the Lagrange equation. Second, to overcome the difficulties in solving the out-of-contact normal vector of the contact surface and the inaccurate application of friction force when solving friction problems using traditional material point contact algorithms. This invention proposes a ball-and-ball compliant contact algorithm. Furthermore, it employs convective particle domain interpolation (CPDI) to eliminate numerical noise caused by material points crossing the mesh in the traditional material point method, and utilizes the ball-and-ball compliant contact algorithm to analyze highly nonlinear impact / contact fracture problems. Finally, this invention aims to address the shortcomings of existing cohesive phase-field finite element methods in analyzing large deformation ductile fracture problems due to severe decreases in numerical accuracy caused by mesh distortion, and the low computational efficiency of existing implicit cohesive phase-field material point methods in analyzing dynamic fracture and large-scale crack propagation problems.

[0008] The technical solution of this invention is an explicit cohesive phase field material point method for dynamic impact / contact toughness fracture analysis based on a novel ball-ball compliant contact algorithm. The specific steps are as follows:

[0009] Step 1: Define the Eulerian background mesh, establish a discrete material point model, define the physical material parameters of the discrete material points, and initialize the material point variables. The physical material parameters include elastic parameters, plastic parameters, and phase field parameters. Elastic parameters include Grüneisen equation of state parameters and shear modulus G; plastic parameters include yield strength σ. Y Hardening modulus E p Phase field parameters include the phase field model regularization parameter l0 and the critical fracture energy release rate G. c Phase field viscosity parameter η; initialization of material point variables including initial material point position vector x p Initial material point displacement vector u p Initial material point velocity vector Initial matter point phase field d p Mass m of a substance point p Initial particle volume V pInitial material point deformation gradient tensor F p and the particle domain vector r of matter points 1,0 r 2,0 and r 3,0 ;

[0010] Step 2: Initialize the Euler background mesh and calculate the interpolation function φ using the CPDI interpolation technique. Ip and the gradient of the interpolation function Establish a mapping relationship between discrete material points and Eulerian background mesh, and map the physical material properties on discrete material points to Eulerian background mesh nodes;

[0011] In the CPDI interpolation technique, each discrete material point is defined to have a parallelepiped particle domain, and the deformation of the particle domain is updated to r. 1,n+1 =F p,n+1 ·r 1,0 r 2,n+1 =F p,n+1 ·r 2,0 and r 3,n+1 =F p,n+1 ·r 3,0 , where (r 1,0 ,r 2,0 ,r 3,0 ) and (r 1,n+1 ,r 2,n+1 ,r 3,n+1 F represents the particle field vector of the matter point at the initial time and the current time, respectively. p,n+1 The deformation gradient tensor represents the current time step; the generalized interpolation function φ represents the geometry at the current time step. Ip Its gradient analytical expression is as follows:

[0012]

[0013] in, Let be the interpolation basis function for the particle domain corner point i of discrete material point p with respect to the Euler background grid. This represents the gradient operator under the current configuration, J is the particle domain vector Jacobi matrix, the subscript I indicates the background mesh node, and N is the coefficient matrix;

[0014] Step 3: Identify the material point using the ball-ball compliant contact algorithm and calculate the contact force. And map it to the background grid node;

[0015]

[0016] in, F represents the total number of particles in contact with particle p. pk and These represent the normal contact force and tangential force exerted on particle p by particle k, respectively.

[0017] Step four: Calculate the internal forces and resultant external forces of the background mesh nodes. Based on the explicit time integration method, the governing equations for the displacement field and phase field are calculated separately using an alternating solution strategy. This is achieved by considering the phase field d at the material points. p The coupled explicit cohesive phase-field model is used to calculate the nodal accelerations of the displacement field. I and node speed The rate of change of the phase field at the nodal points is calculated using the historical strain field function H of the displacement field. The time interval [0, t] from 0 to t is discretized into several time increments Δt. These time increments satisfy the Courant-Friedrichs-Lewy (CFL) condition to ensure the stability of the results. The specific operation method for solving the discrete governing equations of the phase field and displacement field using the staggered solution strategy is as follows:

[0018] For the displacement field, the displacement vector of the material point at time t is... The velocity vector of a matter point Background node velocity vector and background mesh node momentum Update the displacement vector of the material point at time t+Δt The velocity vector of a matter point Node velocity vector The specific steps are as follows:

[0019]

[0020]

[0021] in, For node-centric quality, and Let N represent the nodal external force vector and the nodal internal force vector of the displacement field, respectively. I This represents the total number of background mesh nodes. For nodal momentum, velocity is updated in FLIP mode when κ = 1 and in PIC mode when κ = 0. In this invention, κ = 0.99.

[0022]

[0023] in, N represents the current discrete matter point density. p Represents the total number of discrete material points. Let G represent the Cauchy stress tensor, G represent the gravitational acceleration vector, and h represent the boundary thickness.

[0024] For the phase field, the phase field of the matter point at time t is... Phase field damping on background grid nodes Historical variable H and phase field driving force Update the matter point at time t+Δt.

[0025]

[0026] Where η represents the phase field viscosity coefficient, It is a degenerate function. Let a1, a2, and a3 be the crack geometry function, and a1, a2, and a3 be the cohesive phase field parameters. Furthermore, to ensure that the phase field value d ≤ 1, when... season

[0027] Step four: After calculating the discrete governing equations for the displacement field and phase field, the relevant information on the background mesh is mapped back to the material points to update the relevant variables d of the material points. p , u p F p and H;

[0028] according to The historical strain field function is updated as follows:

[0029]

[0030] in, Let V represent the hydrostatic pressure, V represent the volume, and H(x) be the Heaviside function, which means that when x>0, H(x)=1, and otherwise H(x)=0.

[0031] Step 5: Store and output the relevant variable information, return to Step 2, proceed to the next time step, until the calculation is complete.

[0032] The ductile metal is described by the Grüneisen equation of state and elastoplastic constitutive model for its mechanical behavior under high-speed impact. Equivalent stress is used... It is divided into two parts, one of which is affected by phase field damage. and the part unaffected by phase field damage

[0033]

[0034] Equivalent pressure Calculated using Grüneisen's equations of state, in the following form:

[0035]

[0036] Where μ = ρ / ρ0-1 is the compressibility coefficient, ρ is the current density, ρ0 is the initial density, ζ is the internal energy per unit reference volume, C, γ0, α and S1 are material constants, and H(x) is the Heaviside function, that is, when x>0, H(x)=1, and when x≤0, H(x)=0.

[0037] Equivalent stress The calculation formula is as follows:

[0038]

[0039] Where G is the shear modulus, ∈ is the strain tensor, and e dev For the partial strain tensor, This refers to the plastic component of the deviatoric strain tensor. For equivalent deviatoric stress The second variable, It is the plastic flow factor. It represents the rate of change of a certain value.

[0040] In addition, the von Mises yield function is used for plastic evolution.

[0041]

[0042] The cohesive phase field and displacement field coupled model uses the Lagrange formula to derive the control equations for dynamic fracture of ductile metals.

[0043] L = T + W v -(Φ-W) (26)

[0044] Where, T, Φ, W v W and W represent kinetic energy, potential energy, phase field dissipation energy, and external force potential energy, respectively.

[0045] Potential energy can be divided into fracture energy Φ d Elastic energy Φ el Plasticity Φ pl ,

[0046]

[0047] in, It is a plastic internal variable.

[0048] The definition and decomposition expression of elastic energy are as follows:

[0049]

[0050] in, Let σ be the equivalent deviatoric stress, σ be the equivalent stress, and ∈ be the strain tensor. Let be the deviatoric strain tensor. The elastic strain energy is decomposed into tensile and compressive components, and the expression for the elastic energy after damage is:

[0051]

[0052] Considering elastic energy as the sole driving force for the phase field, the variational form of Φ is derived using the divergence theorem as follows:

[0053]

[0054] in, Let represent the generalized variable, n represent the outward unit normal vector of the boundary, u represent the displacement field, and σ represent the Cauchy stress. Since plastic strain is obtained from the yield function, the terms related to plasticity in the above equation are ignored in the derivation.

[0055] Furthermore, the expressions for kinetic energy, external potential energy, and phase field dissipation energy are as follows:

[0056]

[0057] in, G represents the velocity field, and G represents gravitational acceleration. Represents surface force.

[0058] Based on Lagrange functionals, the governing equations are derived from the Lagrange equations.

[0059]

[0060] Considering the condition that the crack cannot heal, a historical variable H is introduced instead. The governing equations are as follows:

[0061]

[0062] Where Ω is the computational domain of the object. Object boundary.

[0063] The specific implementation process of the described compliant contact model is as follows:

[0064] Search for potential contact particles based on background grid nodes and calculate particle radii.

[0065]

[0066] Where V is the current particle volume and a is the expansion coefficient of the particle radius.

[0067] The contact criterion e is used to determine whether two particles are in contact.

[0068] e = ||X rs ||-R r -Rs <0 (44)

[0069] Among them, X rs R represents the vector pointing from the center of particle r to the center of particle s. r and R s Let represent the sphere radii at the current time t for particles r and s. When the contact criterion is met, contact is considered to have occurred between the two particles, and the contact force is calculated. Assume the bulk modulus K of particle r is... r Compared to the bulk modulus K of particle s s If the particle size is small, then particle r is selected as the primary contactor and particle s as the secondary contactor.

[0070] First, based on Newton's third law, we can obtain...

[0071]

[0072] Among them, S r and S s V represents the contact area between particles r and s, respectively. r and V r ′ Let V be the volume of particle r before and after the collision. s and V s ′ Let R′ be the volume of particle s before and after the collision. r and R′ s Let Δx represent the radius of the sphere at time t+Δt for particles r and s. r and Δx s This represents the change in the radius of the sphere at time t+Δt for particles r and s, i.e., |e| = Δx r +Δx s .

[0073] Substituting formulas (46)-(49) into formula (45) yields

[0074]

[0075] Expand the above equation

[0076]

[0077] Since e is a very small quantity, only retain and After obtaining the first-order terms, we get Δx. r The expression is

[0078]

[0079] Assuming the contact area between the two particles is equal, that is The above equation can then be rewritten as:

[0080]

[0081] For nonlinear materials, since their nature is unknown, we choose to solve for them. make The following calculation formula can be obtained.

[0082]

[0083]

[0084] Among them, EOS r and EOS s Representing the equations of state for particles r and s respectively, ρ r and ρ r ′ represents the initial and current density of particle r, and μ is the relative density increment. Therefore, the normal contact force is...

[0085]

[0086] In order to satisfy the non-penetration condition, the contact area is taken as a conservative value, that is...

[0087] For linear elastic materials, the bulk modulus K r and K s It is known.

[0088]

[0089] To further suppress the velocity oscillation problem during the contact process, a smoothing function is introduced. Normal contact force F rs The calculation formula can be written as

[0090]

[0091] Based on Coulomb's law of friction, frictional force The calculation formula is

[0092]

[0093] Where v represents the coefficient of friction. Represents relative tangential velocity, v s and v r These represent the velocities of particle r and particle s, respectively.

[0094] Based on the above theoretical derivation, the detailed solution formats of the explicit cohesive phase field material point method and the ball-ball compliant contact algorithm proposed in this invention are given. The explicit cohesive phase field material point method for dynamic impact / contact toughness fracture analysis proposed in this invention can be realized by coupling the explicit cohesive phase field model, the convective particle domain interpolation technique and the staggered solution strategy.

[0095] according to Figure 1 The flowchart of the calculation process of the PF-CZM-MPM-PCCA method proposed in this invention is shown below, and its specific implementation process is as follows;

[0096] 1) Establish a discrete model of material points and define material parameters (elastic-plastic: ρ0,C,S1,γ0,a,G,σ). Y Cohesive phase field: G c a2, a3, l0, d p H p Other: x p ,u p ,F p ,…);

[0097] 2) Time step loop: while t <t end ,do

[0098] 2.1) Initialize the background mesh and calculate the interpolation function φ Ip and its gradient

[0099] 2.2) Map the information carried by the material points onto the background mesh;

[0100] 2.3) Perform contact detection and calculate contact force between objects in different background meshes.

[0101] 2.4) Map the contact forces to the background mesh nodes and calculate the internal forces of the nodes. and combined external forces

[0102] 2.5) Update the mesh node velocities and apply Dirichlet boundary conditions;

[0103] 2.6) Update the velocity of matter points Displacement Deformation gradient wait;

[0104] 2.7) Mapping the phase field of the matter points to the background mesh nodes yields d. I ;

[0105] 2.8) Update the material point history variable H p ;

[0106] 2.9) Calculate the nodal phase field increments and update the phase field values ​​of the material points.

[0107] 3) Return to step 2) until the calculation is complete;

[0108] The beneficial effects of the present invention are as follows: (1) The present invention provides an explicit cohesive phase field material point method for dynamic impact / contact tough fracture analysis, which provides a new numerical calculation method for the study of tough fracture failure of metallic materials. Since the method is based on the explicit material point method, it can effectively overcome the mesh distortion problem compared with the traditional mesh method. Therefore, its significant advantage is that it can handle large deformation and contact and other strong nonlinear fracture failure problems well. In addition, with the help of the cohesive phase field model, it can automatically handle complex crack bifurcation, intersection and free expansion in three-dimensional space. Furthermore, the method can be extended to the fracture failure analysis of other materials by replacing the plastic constitutive model, and can be extended to complex multi-field coupled fracture failure analysis, such as thermo-coupled tough fracture analysis, by embedding multi-physics coupling theory.

[0109] (2) This invention provides an explicit cohesive phase field material point method for dynamic impact / contact ductile fracture analysis. In this method, an explicit cohesive phase field fracture model is developed, and the coupling control equation is derived from the energy perspective based on the Lagrange equation, which can effectively predict the ductile fracture behavior of metals.

[0110] (3) This invention provides a ball-and-ball compliant contact algorithm, which combines the advantages of the ball-and-ball algorithm and the compliant force contact algorithm. It avoids the difficulty of calculating the external normal vector of the contact surface in the contact algorithm based on the background mesh under complex deformation, effectively suppresses the problem of excessive oscillation of the contact surface velocity during the contact process, reduces the sensitivity to the discrete size and discrete method of the contact surface, and has good convergence for friction problems.

[0111] (4) The present invention provides an explicit cohesive phase field material point method for dynamic impact / contact toughness fracture analysis. It adopts a phase field-displacement field staggered solution strategy combined with an explicit time integration scheme. Compared with the implicit staggered iterative solution scheme, it greatly reduces the difficulty of numerical implementation and improves the computational efficiency, providing a feasible solution for parallel computational research on large-scale complex fracture failure problems. Attached Figure Description

[0112] Figure 1 This is a schematic diagram of the PF-CZM-MPM-PCCA calculation loop proposed in this invention;

[0113] Figure 2 This is a schematic diagram of the collision of two hyperelastic rings in Embodiment 1 of the present invention;

[0114] Figure 3The equivalent stress cloud diagrams of different contact algorithms at different times in Embodiment 1 of the present invention are shown. (a)-(c), (d)-(f) and (g)-(i) represent the equivalent stress cloud diagrams of the traditional mesh contact algorithm, the traditional particle contact algorithm and the ball-ball compliant contact algorithm at t=0ms, t=0.8ms and t=2ms, respectively.

[0115] Figure 4 The kinetic energy time history curves for different contact algorithms in Embodiment 1 of the present invention are shown below.

[0116] Figure 5 The strain energy time history curves for different contact algorithms in Embodiment 1 of the present invention are shown.

[0117] Figure 6 The total energy time history curves for different contact algorithms in Embodiment 1 of the present invention are shown below.

[0118] Figure 7 This is a schematic diagram of the geometry and boundary conditions of Embodiment 2 of the present invention;

[0119] Figure 8 The velocity contour maps of different times and grid scales in Embodiment 2 of the present invention are shown in (a)-(c) and (d)-(f), which represent the velocity contour maps at t=0μs, t=0.2μs and t=0.6μs respectively when h=5.5μm and h=4.1μm.

[0120] Figure 9 The following are phase field contour maps and equivalent plastic strain contour maps at different grid scales in Embodiment 2 of the present invention. (a)-(c) represent phase field contour maps at h=8.2μm, h=5.5μm and h=4.1μm at t=0.6μs, respectively. (d)-(f) represent equivalent plastic strain contour maps at h=8.2μm, h=5.5μm and h=4.1μm at t=0.6μs, respectively.

[0121] Figure 10 The free surface velocity time history curves under different mesh sizes are shown in Embodiment 3 of the present invention. Detailed Implementation

[0122] The performance of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. The following embodiments are used to illustrate the present invention, but should not be used to limit the scope of application of the present invention.

[0123] To make the objectives, technical solutions, and specific implementation effects of this invention clearer, two specific embodiments are described below in conjunction with the accompanying drawings. Figures 2-10This paper provides a further detailed explanation of the accuracy, reliability, and superiority of the PF-CZM-MPM-PCCA method proposed in this invention. The first embodiment verifies the accuracy and effectiveness of the ball-and-rock compliant contact model proposed in this invention. The second embodiment simulates the delamination problem caused by the collision of two three-dimensional cylinders to demonstrate the ability of the algorithm developed in this invention to handle dynamic tensile fracture of ductile metals.

[0124] (1) Example 1: Collision of two hyperelastic rings

[0125] Example 1 considers the collision of two hyperelastic rings. The reliability of PCCA in terms of energy conservation and collision deformation is verified by comparing the numerical solution obtained with that obtained from a traditional contact algorithm. The geometry and boundary conditions of Example 1 are as follows: Figure 2 As shown. The ring is made of compressible Neo-Hookean material with the following parameters: Young's modulus E = 73 MPa, Poisson's ratio v = 0.4, and density ρ = 10¹⁰ × 10⁻⁶. -12 kg / mm 3 outer radius r out =40mm, inner radius r in =30mm, l=200mm, particle radius expansion coefficient a=2. The traditional mesh contact algorithm, the traditional particle contact algorithm, and PCCA were used for solving the problem. In all simulations, the mesh size was set to 1mm. At the Gaussian intersection of the mesh, four particles were inserted into each mesh to form a loop. The time step Δt was set to 0.2μs, and the total computation time was 4ms.

[0126] The von Mises contour plots obtained from the numerical simulation at each time point are as follows: Figure 3 As shown, at 0.8ms and 2ms, due to the presence of spurious contact, a noticeable gap exists between the two rings in the results obtained using the traditional mesh contact algorithm. However, there is no significant difference between the results obtained using PCCA and the traditional particle contact algorithm. Furthermore, the stress fields obtained by all three contact algorithms are very smooth, which is due to the use of CPDI interpolation, which greatly reduces numerical noise between meshes.

[0127] In addition, from Figure 4-6 It can be seen that PCCA is superior to the other two methods in maintaining energy conservation, and the results obtained in this paper agree well with those calculated by Vaucorbeil et al. using TLMPM. In summary, the above results demonstrate that the contact algorithm proposed in this invention is universal and also illustrate that the assumptions made in the derivation of contact forces are reasonable.

[0128] (2) Example 2: Delamination caused by collision of two cylinders

[0129] The purpose of this embodiment is to verify the accuracy of the cohesive phase field material point method in solving the problem of metal delamination. This embodiment consists of three parts: the first part is the simulation of delamination under abrupt pressure load; the second part is the simulation of delamination under two-dimensional flyer impact; and the third part is the simulation of delamination under three-dimensional cylinder impact. All simulations use the material parameters of oxygen-free high-conductivity copper, and the specific parameters are shown in Table 1.

[0130] Table 1. Parameters of Oxygen-Free High Conductivity Copper Materials

[0131]

[0132] This embodiment considers the problem of delamination caused by the collision of two cylinders, aiming to test the reliability of the method developed in this invention in solving three-dimensional delamination problems, and further demonstrate the ability of this invention to handle large-scale problems. Based on the experiments of Peng et al., the geometric dimensions and boundary conditions are set as follows: Figure 7 As shown. The cylinders have lengths of 0.42 mm and 0.82 mm, and radii of 0.328 mm. The collision of the two cylinders occurs within a sleeve tightly fitted to the cylinders. The left cylinder collides with the stationary right cylinder at a speed of 246 m / s. The background mesh sizes are h = 8.2 μm, h = 5.5 μm, and h = 4.1 μm. The initial positions of the matter points are located at the eight Gaussian points of the 3D mesh. The number of matter points corresponding to the three mesh sizes are 381728, 1283912, and 3039520, respectively. The expansion coefficient of the particle radius is set to a = 1. The time step is set to 0.2 ns, and the total computation time is t = 0.8 μs. The phase field viscosity coefficient is η = 1 × e -7 N·s / mm 2 Critical energy release rate G c =0.8 N / mm, material strength f t =1.3 GPa. The phase field characteristic width l0 is maintained at three times the grid width, which are 0.0246 mm, 0.0165 mm and 0.0123 mm respectively.

[0133] Figure 8 The diagrams show velocity contours at different times for two scenarios: h = 5.5 μm and h = 4.1 μm. As can be seen from the diagrams, the velocity propagates axially after the two cylinders collide. Upon reaching the free surface, the velocity is reflected, forming a tensile wave that creates a crack at the midpoint. Figure 9(a) The phase field of the final state of the cylinder at three mesh scales is presented. As can be seen from the figure, the intermediate layer crack region shrinks continuously with decreasing mesh scale and l0. Simultaneously, the mottled damage on the cylinder increases with decreasing mesh scale. This is partly because, with the increase in the number of material points, the overlap of reflected waves after delamination is more accurately characterized, capturing more damage. On the other hand, this is due to the oscillation caused by the weakening of PCCA compliance when the particle radius decreases. Furthermore, the equivalent plastic strain cloud after delamination is as follows... Figure 9 As shown in (b), the equivalent plasticity of the fracture region increases proportionally as the mesh size decreases. This phenomenon can be attributed to the increased tensile force applied to the particles in the fracture region as the fracture region shrinks.

[0134] at last, Figure 10 The velocity-time curves of the freeform surface at three grid scales are presented. As can be seen from the figures, the velocity-time curves of the freeform surface at each grid scale exhibit good characteristics of small-block fracturing and good convergence.

[0135] In summary, the two embodiments described above have verified the accuracy and effectiveness of the PF-CZM-MPM-PCCA proposed in this invention from different perspectives and levels, and have demonstrated the significant advantages of this method in terms of accuracy and computational scale, showcasing the broad application prospects of this invention.

[0136] The embodiments of the present invention are given for illustrative and descriptive purposes only, and are not intended to be exhaustive or to limit the invention to the forms disclosed. Many modifications and variations will be apparent to those skilled in the art. The embodiments were chosen and described to better illustrate the principles and practical application of the invention, and to enable those skilled in the art to understand the invention and design various embodiments with various modifications suitable for a particular purpose.

Claims

1. A dynamic impact / contact toughness metal fracture analysis explicit cohesive force phase field material point method, characterized in that, The steps are as follows: Step one, defining Euler background mesh, establishing discrete material point model and defining discrete material point variables, initializing material point variables; discrete material point variables include physical material parameters, elastic parameters, plastic parameters and phase field parameters, wherein, physical material parameters include elastic parameters, plastic parameters and phase field parameters; elastic parameters include Grüneisen equation of state parameters and shear modulus G; plastic parameters include yield strength σ Y and hardening modulus E p ; phase field parameters include phase field model regularization parameter l0, critical fracture energy release rate G c and phase field viscosity parameter η; initializing material point variables include initial material point position vector x p,0 , initial material point displacement vector u p,0 , initial material point velocity vector initial material point phase field d p , material point mass m p , initial material point volume V p , initial material point deformation gradient tensor F p and material point particle domain vectors r1, r2 and r3; Step two, initialize the Euler background grid, calculate the generalized interpolation function φ by CPDI interpolation technology Ip and interpolation function gradient The mapping relationship between the discrete material point variable and the Euler background grid is established, and the physical material parameters on the discrete material point variable are mapped to the Euler background grid nodes. In the CPDI interpolation technique, each discrete material point variable is defined with a parallelepiped material point particle domain, whose deformation update is given by 1,n+1 = F p,n+1 · r 1,0 , r 2,n+1 = F p,n+1 · r 2,0 , and r 3,n+1 = F p,n+1 · r 3,0 , where (r 1,0 , r 2,0 , r 3,0 ) and (r 1,n+1 , r 2,n+1 , r 3,n+1 ) represent the material point particle domain vectors at the initial and current time instants, respectively, and F p,n+1 denotes the deformation gradient tensor at the current time instant; the analytical expressions of the generalized interpolation function φ Ip and its gradient at the current time geometry are given by wherein is the interpolation basis function of the particle domain angular point i of the discrete material point variable p with respect to the Euler background grid, is the gradient operator at the current configuration, J is the material point particle domain vector Jacobi matrix, the subscript I denotes the background grid node, N is the coefficient matrix, S T is the transpose of the interpolation function matrix; Step three, distinguish discrete material point variables by billiard soft contact algorithm, calculate contact force and map it to the Euler background grid nodes; wherein, represents the total number of particles in contact with the discrete material point variable, i.e. particle p, F pk and respectively represent the normal and tangential contact forces on particle p by particle k; Step four, calculate the internal force and the resultant external force of Euler background grid nodes; based on the explicit time integration method, the control equations of displacement field and phase field are calculated respectively through staggered solution strategy; wherein the node acceleration of displacement field is calculated by considering the coupling of discrete material point phase field d p and the explicit cohesive force phase field model I and the node velocity of the node phase field change rate is calculated by the history strain field function H of the displacement field Disperse the time interval [0, t] from 0 to t into several time increments Δt, and the time increment satisfies the CFL condition to ensure the stability of the results; the specific operation mode of solving the discrete control equations of phase field and displacement field respectively through staggered solution strategy is as follows: For the displacement field, the material point displacement vector at time t is updated by the material point velocity vector at time t material point velocity vector background node velocity vector and background mesh node momentum the material point displacement vector at time t+At is updated material point velocity vector node velocity vector The specific operation is as follows: in, For node-centric quality, and Let N represent the nodal external force vector and the nodal internal force vector of the displacement field, respectively. I This represents the total number of background mesh nodes. For nodal momentum; when κ=1, velocity is updated using the FLIP method, and when κ=0, velocity is updated using the PIC method; in, N represents the current discrete matter point density. p Represents the total number of discrete material points. Let G denote the Cauchy stress tensor and G denote the gravitational acceleration vector. For surface force, h represents the boundary thickness. Let p be the volume of particle p. For the Laplace operator; For the phase field, the phase field of the matter point at time t is... Phase field damping on background grid nodes Historical variables and phase field driving force Update the matter point at time t+Δt; Where η represents the phase field viscosity coefficient, It is a degenerate function. Let G be the crack geometry function, a1, a2, and a3 be the cohesive phase field parameters, and G be the crack geometry function. c Let l0 be the critical energy release rate, and l0 be the characteristic length of the crack-dispersion region; furthermore, to ensure that the phase field value d ≤ 1, when... season Step four: After calculating the discrete governing equations for the displacement field and phase field, map the relevant information on the Eulerian background grid back to the discrete material point variables to update the discrete material point variable d. p , u p F p and H p ; according to The historical strain field function is updated as follows: in, Let V be the hydrostatic pressure, V represent the volume, and H be the Heaviside function, i.e., when... hour, Conversely, Step 5: Store and output the relevant variable information, return to Step 2, proceed to the next time step, until the calculation is complete.

2. The explicit cohesive phase field material point method for dynamic impact / contact toughness metal fracture analysis according to claim 1, characterized in that, Step four, the cohesive phase-field material point method, combines the cohesive phase-field model and the material point method, as detailed below: First, based on the cohesive phase-field model, the governing equations for dynamic fracture of ductile metals are derived using the Lagrange formula. L = T + W v - (Φ - W) (18) where T, Φ, W v and W represent kinetic energy, potential energy, phase field dissipation energy and external force potential energy, respectively. The potential energy is divided into the breaking energy Φ d , the elastic energy Φ el and the plastic energy Φ pl : in, It is a plastic internal variable; Elastic strain energy Φ el The definition and its decomposition expression are: in, Let σ be the equivalent deviatoric stress, σ be the equivalent stress, and ∈ be the strain tensor. Let be the deviatoric strain tensor; the elastic strain energy is decomposed into tensile and compressive components, and the expression for the elastic energy after damage is: Where ω(d) is the phase field degradation function; Considering the elastic energy as the sole driving force for the phase field, the variational form of Φ is derived using the divergence theorem as follows: in, Represents a generalized variable, n represents the outward unit normal vector of the boundary, and u represents the displacement field; Furthermore, the expressions for kinetic energy, external potential energy, and phase field dissipation energy are as follows: in, Represents the velocity field, and ρ represents the current density of the material. Represents surface force, The rate of change of the phase field; Based on the Lagrange functional, the governing equations are derived from the Lagrange equations: Considering the condition that the crack cannot heal, a historical variable H is introduced instead. The governing equations are as follows: On the domain Ω (31) At the border (32) On the domain Ω (33) At the border (34) Where Ω is the computational domain of the object. Object boundary.

3. According to the explicit cohesive phase field material point method for dynamic impact / contact toughness metal fracture analysis as described in claim 1, the specific implementation process of the bouncy contact model in step three is as follows: The contact criterion e is used to determine whether two particles are in contact: e = ||X rs ‖-R r -R s <0(35) in, X rs denotes the vector from the sphere center of particle r to the sphere center of particle s, R r and R s denotes the sphere radius of particle r and particle s at the current time t; when the contact criterion is satisfied, it is considered that contact occurs between the two particles, and the contact force is calculated; it is assumed that the bulk modulus K r of particle r is smaller than the bulk modulus K s of particle s, particle r is selected as the main contact body, and particle s is selected as the secondary contact body; First, based on Newton's third law, we get: where S r and S s are the contact areas of the particle r and the particle s, respectively, V r and V r ′ are the volumes of the particle r before and after the collision, V s and V s ′ are the volumes of the particle s before and after the collision, K r and K s are the bulk modulus of the particle r and the particle s, respectively, R ′ r and R s ′ denote the spherical radii of the particle r and the particle s at the time t+Δt, Δx r and Δx s denote the spherical radius change amounts of the particle r and the particle s at the time t+Δt, i.e., |e| = Δx r + Δx s ; Substituting formulas (37)-(40) into formula (36), we get: Expanding the above equation: Since e is a very small quantity, only retain and After obtaining the first-order terms, we get Δx r The expression is: Assuming the contact area between the two particles is equal, that is... The above formula can then be rewritten as: For nonlinear materials, since their nature is unknown, we choose to solve for them. make The following calculation formula is obtained: where EOS r and EOS s represent the equation of state of particle r and particle s, respectively, and p r and p r ′ are the initial and current densities of particle r, respectively. The normal contact force is then given by: f rs = S r EOS 1 (μ r ) (47) In order to satisfy the non-penetration condition, the contact area is taken as a conservative value, that is... For linear elastic materials, the bulk modulus K r and K s are known; To further suppress the velocity oscillation problem during the contact process, a smoothing function is introduced. Normal contact force F rs The calculation formula can be written as Based on Coulomb's law of friction, frictional force The calculation formula is Where ν represents the coefficient of friction, Represents relative tangential velocity, v s and v r These represent the velocities of particle r and particle s, respectively.

Images