Explicit material point method for dynamic analysis of large deformation of weakly compressible materials

Abstract

The present application belongs to the technical field of soft material weakly compressible large deformation, and proposes a weakly compressible material large deformation dynamic analysis explicit material point method, which provides a new numerical calculation method for the study of weakly compressible material large deformation. The method discretizes the weak form control equation in the total Lagrangian material point method framework, avoiding the challenge of solving the large deformation problem due to grid distortion in the traditional method. In addition, the present application proposes a volume self-locking processing technology based on the F-bar method, which can simply and effectively solve the volume locking problem of weakly compressible material in numerical simulation, and can be easily embedded into the existing material point method program. At the same time, the construction method of the smoothing kernel interpolation function is developed, which can effectively improve the problem of the decline of the solving accuracy of the traditional material point method in the boundary area. The method can solve the large deformation problem of weakly compressible material in a stable and efficient way, and provides a new technical path for the research and application of related fields.

Description

Technical Field

[0001] This invention belongs to the field of large deformation analysis technology for soft matter, and particularly relates to an explicit material point method for dynamic analysis of large deformation of weakly compressible matter. Background Technology

[0002] Soft matter is widely found in nature, such as biological tissues, and in advanced engineering applications, such as elastomers and polymer gels used in stretchable electronics, tire manufacturing, food packaging, cosmetics, biomedicine, drug delivery systems, and soft actuators. Numerical modeling of these materials presents significant challenges due to their ability to undergo large, recoverable elastic deformations. Traditional mesh-based methods (such as the finite element method) solve the problem by discretizing the problem domain into a finite number of small elements, where the accuracy of the solution is highly dependent on the mesh quality. However, under large deformation conditions, severe mesh distortion degrades computational accuracy and necessitates the application of computationally expensive re-meshing techniques. To overcome these challenges, various meshless methods have been developed, including smoothed particle hydrodynamics, ambient dynamics, and the matter point method. Among these methods, the matter point method shows great potential because it effectively combines the advantages of Lagrangian and Eulerian methods. In the matter point method, the physical domain is discretized by matter points (particles) arranged within a background mesh. These Lagrangian particles move on a fixed Eulerian mesh, where the governing equations are solved efficiently.

[0003] As an effective large deformation modeling method, the Material Point Method (MPM) was first proposed by Sulsky et al. and has been widely validated. Standard MPMs typically use linear shape functions to transfer physical properties between particles and background mesh nodes. However, when particles cross the background mesh (i.e., element crossing), so-called "cross-mesh noise" arises due to the discontinuity of derivatives at the background mesh boundaries. This can even lead to energy oscillations or energy non-conservation, resulting in simulation failure. To address this issue and improve the accuracy of large deformation simulations, various interpolation methods have been developed within the MPM computational framework, such as the Generalized Interpolation Material Point Method (GIMP), the Convective Particle Domain Interpolation Method (CPDI), and the B-spline Material Point Method. It is important to note that these interpolation methods are usually based on regular background mesh constructions, which can reduce accuracy when dealing with problems where geometric boundaries do not overlap with the background mesh. Furthermore, designing a suitable background mesh for a specific problem geometry is challenging because the motion of Lagrangian particles within the mesh can cause even a carefully designed mesh to fail. This failure increases the complexity of maintaining accuracy and stability in dynamic scenes where particle positions change significantly over time. With the development of the Total Lagrangian Matter Point Method (TLMPM), this problem is expected to be solved. Unlike the standard MPM, TLMPM solves the governing equations on the initial configuration, meaning that particles no longer traverse the mesh. Furthermore, this method has been successfully applied to the simulation of large deformation problems with good results. In addition, TLMPM only requires the construction of first-order shape functions during the computation, thus improving the computational efficiency of the MPM.

[0004] When dealing with approximately incompressible materials, the solution framework of the Material Point Method (MPM) is susceptible to volume locking, exhibiting excessive stiffness that leads to errors in strain and stress field calculations. This problem also exists in the Finite Element Method (FEM) and related methods, because overly strong incompressible constraints at integration points reduce the accuracy of element stiffness calculations. In the FEM domain, various methods have been proposed to alleviate locking by relaxing the incompressible constraints of element kinematics in different ways. Only those methods compatible with finite deformation kinematics and movable integration points are suitable for MPM. To address volume locking, a series of techniques and strategies have been proposed within the MPM framework, such as the hybrid displacement-pressure (up) formula and the hybrid displacement-pressure-Jacobi (upJ) formula. However, these locking mitigation methods typically require significant modifications to the existing MPM framework and are often optimized for specific basis functions or material types. Furthermore, hybrid formulas require changes to the standard governing equations, which significantly increases the cost of implementation and use. Other methods include operator splitting algorithms and the B-bar method, but their existing versions typically require significant modifications to the time-stepping scheme or basis functions interacting with adjacent elements when applied to MPMs. In contrast, the F-bar method has gained widespread attention due to its simple derivation and has been successfully applied to updating Lagrangian MPMs to alleviate volume locking problems. This volume locking handling strategy based on the F-bar method uses the standard particle-mesh transfer scheme in MPMs to calculate assumed deformations and has shown good potential in simulating dynamic large deformations of weakly incompressible materials.

[0005] Therefore, this invention proposes a fully Lagrangian explicit matter point method for weakly compressible materials to evaluate the dynamic deformation behavior of soft matter under large deformation conditions. To improve computational efficiency and avoid element crossing errors, the fully Lagrangian matter point method is used for discretization and integration of the momentum equation in the physical domain. Furthermore, an interpolation function is constructed using a Gaussian kernel function based on weighted least squares (WLS) kernel correction, enabling the method to handle irregular background meshes or arbitrarily arranged nodes. In addition, the F-bar method is further developed within the fully Lagrangian MPM framework to address the volume locking problem. This method not only significantly improves the computational accuracy of the matter point method at the physical domain boundary but also efficiently solves the volume locking problem, while requiring minimal modification to the existing computational framework. Summary of the Invention

[0006] To address the need for efficient numerical methods for dynamic analysis of large deformation in soft materials, this invention innovatively proposes an explicit material point method (STLMPM) for the dynamic analysis of large deformation in weakly compressible materials. The objectives are twofold: First, to improve the computational efficiency of the material point method and address the numerical oscillation problem that may be caused when material points cross the background grid, this invention solves the discrete governing equations based on the initial configuration, and simultaneously calculates the interpolation function and its gradient in a single step at the beginning of the computation. Second, to solve the problem of reduced computational accuracy of the material point method near the physical domain boundary, a weighted least squares (WLS) method is employed. A Gaussian kernel function with kernel correction is used to establish the interpolation formula between particles and the background mesh. In addition, to alleviate the volume locking problem caused by the near-incompressibility of soft materials, the F-bar method is further developed within the framework of TLMPM to alleviate the overly stiff displacement and stress fields caused by the incompressibility of materials. Finally, this invention aims to overcome the limitations of existing finite element methods in analyzing large deformation problems of soft materials due to mesh distortion, while also making up for the shortcomings of existing explicit material point methods in analyzing weakly compressible problems, such as low computational efficiency and the need for extensive modifications to constitutive equations or the original computational framework.

[0007] The technical solution of this invention: an explicit material point method for dynamic analysis of large deformation of weakly compressible materials, the specific steps of which are as follows:

[0008] Step 1: Define the Euler background mesh. The initial coordinates of the Euler background mesh nodes are X. I Establish a discrete material point model and define the physical material parameters of the discrete material points. Initialize the material point variables, including the initial material point position vector X. p Initial material point displacement vector u p Initial material point velocity vector The mass of a substance point m p Initial particle volume V p 0 Initial material point deformation gradient tensor F p The point shear modulus μ and Poisson's ratio υ of the substance;

[0009] Step 2: Establish the mapping relationship between discrete material points and the Eulerian background mesh through a smoothing kernel function, and map the physical material parameters on the discrete material points to the nodes of the Eulerian background mesh.

[0010] Using a total Lagrangian computational framework, the interpolation function and its gradient are calculated on the initial configuration (the configuration without deformation), and are denoted as follows: and in This represents the spatial gradient operator on the initial configuration. The subscript Ip indicates the interpolation function established at material point p and its gradient value at grid node I. The specific calculation method is as follows:

[0011] Using the Gaussian kernel function as the kernel function for constructing the interpolation function, its expression is:

[0012]

[0013] In this invention, the parameters are selected as k = 1, r is the "unit radial distance", and α is 2.0 to 2.5.

[0014] Unit radial distance r and node weight Calculated as:

[0015]

[0016] Where ||■|| represents the L2 norm of the symbol, R p The radius of the support domain of the material point is represented by I, and the subscripts I and p represent the physical quantities on the background mesh node and the material point, respectively.

[0017] According to the chain rule, the spatial gradient of the node weights on the initial configuration is written as:

[0018]

[0019] The interpolation function is reconstructed using the weighted least squares method; the reconstructed interpolation function and its spatial gradient on the initial configuration are denoted as follows: and It is represented as:

[0020]

[0021] Wherein, the correction matrix C(X,X) I Using the weighted least squares approach, it can be written as:

[0022] C(X,X I ) = B(XX p )M(X p ) -1 B T (X I -X p (5) Where X is the spatial coordinate in the initial configuration, and B is the kernel basis function vector; the three-dimensional linear kernel basis function expression is B(X) = [1, X1, X2, X3], where X1, X2, X3 are the lengths in the three dimensions under the initial configuration; matrix M is calculated as:

[0023]

[0024] Where, n I This represents the total number of background grid nodes; then the inverse of matrix M is denoted as:

[0025]

[0026] Where C1 is a scalar; C2 is a vector of length dim, where dim represents the dimension of the problem; C3 is a dim × dim matrix; the reconstructed interpolation function is calculated as:

[0027]

[0028] Its corresponding spatial gradient is calculated as follows:

[0029]

[0030] Step 3: Calculate the time step for the explicit material point method;

[0031] When performing numerical simulations of the explicit dynamics of weakly compressible materials, the time step Δt needs to satisfy the CFL condition to ensure the stability of the results. That is, the length of stress wave propagation should not exceed one grid in a single simulation step, which can determine the maximum time step of the simulation. At the same time, dynamic effects also need to be considered.

[0032] First, the matter wave velocity c in the CFL condition can be calculated as follows:

[0033]

[0034] Where ρ is the initial density of the substance point; K is the bulk modulus of the substance point, calculated as... μ is the shear modulus of the material. Then, the time increment Δt of the current time step is calculated as:

[0035]

[0036] Among them, C r h is a scalar less than 1, used as a scaling parameter for the time increment; g It is the size of the Euler background mesh defined in step one; The velocity of the material point is represented by c; c is the speed of sound of the material.

[0037] Step four: Map the information at the material points onto the background mesh using the interpolation function and its gradient, and calculate the nodal momentum. Node mass m I Node internal forces Calculate nodal external forces using boundary conditions.

[0038] Mapping the current moment's momentum and mass of the matter point onto the background mesh is specifically written as:

[0039]

[0040] in, V represents the initial density of a point mass. p 0 n represents the initial volume of a point mass. p This represents the total number of matter points. This represents the momentum of a grid node. The velocity m represents the current velocity of a grid node. I The quality of a mesh node is represented by the following formula:

[0041]

[0042] Where, m p This represents the mass of the substance point p;

[0043] The nodal internal forces are calculated using the stress information carried by the material point at the current moment, and are calculated using the first Pikachuv stress as follows:

[0044]

[0045] in, This represents the first Pikachuff stress at the material point at the current moment; based on the applied external force boundary conditions, the external forces at the mesh nodes are calculated as follows:

[0046]

[0047] in, This represents the force per unit volume. This represents the force per unit area. It is the virtual boundary layer thickness, used for applying force boundary conditions in two-dimensional problems;

[0048] Step 5: Apply Dirichlet boundary conditions and set the current moment of the mesh node momentum and nodal internal forces at the point where the boundary conditions are applied to 0.

[0049] Step 6: Update the mesh node velocity field based on the internal and external forces of the mesh nodes and the velocity at the current moment; and update the velocity, position, and deformation gradient tensor of the material points through the updated node velocity field. The specific operations are as follows;

[0050] The velocity of the grid node at the next moment Calculated as:

[0051]

[0052] in, pass The calculation is obtained; then the velocity and position of the material point are updated as follows:

[0053]

[0054] and

[0055]

[0056] Then, based on the fundamental theory of continuum mechanics, the deformation gradient tensor of the material point is calculated using the current velocity of the grid node and the current time step:

[0057]

[0058] in, Let represent the deformation gradient tensor of the material point p at time t;

[0059] The Jacobian and volume of the matter point are updated as follows:

[0060]

[0061] Step 7: Perform incompressible processing, use the invented F-bar method combining incremental full format to average the volume part of the deformation gradient tensor of the material point, and bring the averaged deformation gradient tensor into the constitutive equation to update the stress of the material point.

[0062] First, the Jacobian at the material point is weighted according to the initial volume and mapped onto the Euler background mesh nodes to obtain the Jacobian of the Euler background mesh nodes, which is written in the following form:

[0063]

[0064] Where β is the control parameter, It is the Jacobian of the average deformation gradient tensor at the current moment. It is the Jacobian of the deformation gradient increment at the current time step, specifically expressed as follows:

[0065]

[0066] in, This represents the average deformation gradient tensor of the material point p at time t;

[0067] Remap the Jacobian on the Euler background mesh nodes back to the material points, and calculate the average Jacobian on the material points, in the following form:

[0068]

[0069] Then, Calculated as:

[0070]

[0071] The strain energy density function of the material is φ, which is a function of the deformation gradient tensor and is denoted as φ = φ(F). In the F-bar method, the average deformation gradient is used. Replacing the deformation gradient F in the strain energy density function, the strain energy density function is thus reconstructed as follows: The first Pikachuff stress at the material point is updated to:

[0072]

[0073] Then, Cauchy stress Calculated as:

[0074]

[0075] Where the superscript T denotes the transpose of the matrix;

[0076] Step 8: Store and output the relevant variable information, return to Step 3, proceed to the next time step, until the calculation is complete.

[0077] 2. The explicit material point method for dynamic analysis of large deformation of weakly compressible materials according to claim 1, characterized in that,

[0078] The smooth kernel function constructs an interpolation function that adapts to background meshes of arbitrary shapes and background mesh nodes arranged according to arbitrary rules; the incompressible processing includes full-volume averaging and incremental averaging formats for the volume portion.

[0079] The smoothing kernel function used is a Gaussian kernel function. In the material point method, an arbitrary solid domain is discretized into a series of material points, and a background mesh is used to solve the discrete form of the governing equations. When the boundary of the solid domain does not coincide with the background mesh, the accuracy of the calculation results will decrease due to initial geometric modeling errors and inaccurate application of boundary conditions. The Gaussian kernel function is used to assign mesh nodes for physical information mapping to the material points and calculate the initial weights. It does not satisfy the normalization condition and the linear field reproduction condition, as expressed as:

[0080]

[0081] Where, n I This represents the total number of grid nodes; correspondingly, the condition that its spatial gradient should satisfy is written as:

[0082]

[0083] Using weighted least squares pairs Reconstruct the data and calculate the results using equations (8) and (9). and By using the reconstructed interpolation function The node values ​​u of the continuous function u(X) I =u(X) I Approximately:

[0084]

[0085] When u(X) = B(XX) p According to equations (4)-(6) and (35), we get:

[0086]

[0087] Therefore, if a linear basis B(XX) is used p The reconstructed interpolation function satisfies the normalization condition (31) and the linear field reproduction condition (32). For two-dimensional problems, each material point is equipped with at least 3 background grid nodes as interpolation nodes, and for three-dimensional problems, each material point is equipped with at least 4 grid nodes as interpolation nodes.

[0088] The aforementioned incompressible processing formula employs a hybrid format combining full-volume averaging and incremental averaging: In the F-bar method, the volume component of the deformation gradient tensor is averaged. Within the material point method computation framework, the volume component of the deformation gradient tensor is weighted and mapped onto the background mesh, and then mapped back from the background mesh to the material points to calculate the average deformation gradient tensor. First, the deformation gradient tensor is decomposed into a bias component and a volume component:

[0089] F = F d ·F v (37)

[0090] Among them, F d and F v Let represent the biased and volumetric components of the deformation gradient tensor, respectively, and define them as follows:

[0091] F d =J -1 / dim FF v =J 1 / dim I (38) Where I represents the unit tensor; by decomposing the deformation gradient tensor using equation (38), the partial and volume components obtained satisfy the following conditions:

[0092] det(F d )=1 det(F v )=det(F)=J (39)

[0093] Then, the average deformation gradient tensor is obtained by replacing the volume part with the average volume part, and is expressed as:

[0094]

[0095] in, This represents the average volume portion. Let Jacobian represent the average deformation gradient tensor and be calculated as...

[0096] The full average format maps the Jacobian of the material points to the background mesh to obtain the node Jacobian. The specific steps are as follows:

[0097]

[0098] Wherein, the subscript tot indicates the node Jacobian obtained using the full average format;

[0099] The incremental averaging scheme, the relationship between the deformation gradient increments ΔF and F, and their corresponding averaging forms, can be written as follows:

[0100]

[0101] therefore, and Jacobi is represented as:

[0102]

[0103] Substituting equations (42) and (43) into equation (40), we get

[0104]

[0105] Then, the node Jacobian in the incremental format is defined as:

[0106]

[0107] By controlling the parameter β and combining the full and incremental formats, the Jacobian of the background mesh nodes is calculated as follows:

[0108]

[0109] The average Jacobian of the material points is calculated using equation (27); then, the average deformation gradient tensor is calculated using equation (28).

[0110] The specific implementation processes for constructing the interpolation function and the weak compressibility processing technique are as follows:

[0111] The specific implementation process for constructing the interpolation function is as follows:

[0112] Step 1: Set the support domain radius R for each material point. p ;

[0113] Step 2: Set the initial X coordinates of the material point p and the initial X coordinates of the grid nodes I ;

[0114] Step 3: Set parameters k and α, and use formulas (1) and (2) to calculate the weight of the grid nodes in the support domain for each material point.

[0115] Step 4: Calculate matrix M according to formula (6) and calculate its inverse matrix M. -1 ;

[0116] Step 5: Calculate C1, C2 and C3 according to formula (7);

[0117] Step 6: Calculate the interpolation function according to formulas (8) and (9) respectively. and the gradient of the interpolation function

[0118] The specific implementation process for weak compressibility is as follows:

[0119] Step 1: Calculate the increment of the deformation gradient tensor according to formulas (24) and (25). and Jacobi

[0120] Step 2: Calculate the Jacobian matrix of the background mesh nodes according to formulas (41), (45) and (46).

[0121] Step 3: Calculate the average deformation gradient tensor of the material point according to formula (28).

[0122] Based on the finite strain assumption, such as Figure 1 As shown, the midpoint X of the reference configuration Ω will be mapped to the current configuration Ω at time t through the deformation mapping function x = x(X,t). t At the midpoint x, the deformation gradient tensor is

[0123]

[0124] By introducing the strain energy density function φ and considering the work done by kinetic energy and external forces, the total energy function can be written as:

[0125]

[0126] Where b0 and t0 represent the volume force and surface force, respectively, and v represents the velocity. By performing variational analysis on the total energy function, the strong-form governing equations and boundary conditions can be written as follows:

[0127]

[0128] Where ü represents the second derivative of u with respect to time t. It is the boundary of Ω0. And n0 is the unit outward normal vector of the boundary in the initial configuration.

[0129] Then, the weak form governing equations on the initial configuration can be written as follows:

[0130]

[0131] Here, w is the weighting function of the displacement field. Within the MPM computational framework, the weakly integral form of the governing equations is further rewritten as a discrete summation form of the governing equations.

[0132]

[0133] Where, n p This represents the total number of matter points, where the subscripts p and I denote variables related to matter points and background grid nodes, respectively. Considering the arbitrariness of and using a lumped mass matrix, the discrete form equation for each background grid node in the TLMPM framework can be written as:

[0134]

[0135] Based on the above theoretical derivation, a detailed solution format for the total Lagrangian material point method under the large deformation framework proposed in this invention is given. It can realize the explicit material point method for large deformation dynamic analysis of weakly compressible materials to be developed in this invention by constructing an interpolation function for arbitrary structure background mesh using Gaussian kernel function, F-bar incompressibility correction technique and forward Euler time discretization strategy.

[0136] according to Figure 6 The flowchart of the STLMPM method proposed in this invention is shown below, and its specific implementation process is as follows:

[0137] 1) Establish a discrete model of material points and define material parameters (μ, υ, X) p ,u p ,v p,0 ,m p V p 0 ,F p ,P p ,β);

[0138] 2.1) Set parameters k and α, and use formulas (1) and (2) to calculate the weight of the grid nodes in the support domain for each material point.

[0139] 2.2) Calculate matrix M and its inverse matrix M according to formula (6). -1 ;

[0140] 2.3) Calculate C1, C2 and C3 according to formula (7);

[0141] 2.4) Initialize the background mesh: Interpolate the function according to formulas (8) and (9) respectively. and gradient

[0142] 3) Time step loop whilet <tfinal

[0143] 3.1) Calculate the time increment Δt of the current time step according to formula (11);

[0144] 3.2) Calculate the nodal momentum according to formulas (12)-(15) mass m I Internal strength and external forces

[0145] 3.3) Apply Dirichlet boundary conditions;

[0146] 3.4) Update the nodal momentum according to formula (16). And calculate the node velocity

[0147] 3.5) Update the substance point variables according to formulas (17) and (18): and

[0148] 3.6) Calculate the increment of the deformation gradient tensor according to formulas (24) and (25). and Jacobi

[0149] 3.7) Update the deformation gradient tensor of the material point and volume

[0150] 3.8) Calculate the nodal Jacobian according to formulas (41), (45) and (46).

[0151] 3.9) Calculate the average Jacobian of the material points according to formula (27).

[0152] 3.10) Calculate the average deformation gradient tensor of the material point according to formula (28).

[0153] 3.11) Update the material point stress according to formula (29)

[0154] 4) Return to step 3) until the calculation is complete;

[0155] Where tfinal represents the total simulation time.

[0156] The beneficial effects of the present invention are as follows: (1) The present invention provides an explicit material point method for dynamic analysis of large deformation of weakly compressible materials, which provides a new numerical calculation method for the study of large deformation dynamics of weakly compressible materials. Since the method is equipped with the explicit material point method, it can effectively overcome the mesh distortion problem compared with the traditional mesh-based method. Therefore, its significant advantage is that it can handle large deformation and strong nonlinear problems well. Moreover, the method solves the control equation on the initial configuration, and the interpolation function only needs to be calculated once at the beginning of the simulation, which greatly improves the calculation efficiency of the traditional material point method. Furthermore, the method can be applied to the dynamic analysis of large deformation of incompressible materials by extending various constitutive models, and can be extended to complex multi-field coupling analysis by embedding multi-physics coupling theory, such as dynamic fracture analysis of large deformation of solids.

[0157] (2) The present invention provides an explicit material point method for dynamic analysis of large deformation of weakly compressible materials. In this method, an interpolation function construction method based on a smooth kernel function and weighted least squares correction is developed. The smooth kernel function assigns weights to the background grid nodes in the support domain of each material point, and the weighted least squares method corrects the weights to satisfy the normalization condition and the linear field reproduction condition. This interpolation function is applicable to any irregular background grid nodes, providing an efficient and feasible method to improve the computational accuracy of the material point method at the physical domain boundary. Furthermore, due to the inclusion of a total Lagrangian computational framework, this interpolation function construction method hardly increases the total simulation time and is easy to embed into existing material point programs. The proposed method not only provides a new perspective for the numerical simulation of large deformation of weakly compressible materials in theory, but also shows important application potential in engineering practice.

[0158] (3) The present invention provides an explicit material point method for dynamic analysis of large deformation of weakly compressible materials. This method uses the F-bar method to average the volume portion of the deformation gradient tensor at the material point, effectively solving the volume locking problem caused by the incompressibility of the material. Compared to multi-field variational methods, it only requires one set of interpolation functions, making numerical implementation simpler and requiring no major modifications to the computational framework of the traditional material point method. Furthermore, by replacing the original deformation gradient with the averaged deformation gradient tensor, this method is directly applied to constitutive model calculations, thus easily extending to any constitutive model and adaptive algorithm, demonstrating extremely high adaptability and flexibility. In addition, this incompressibility handling method has good versatility, not only applicable to the material point method but also naturally compatible with various numerical calculation methods, such as the finite element method and near-field dynamics methods, providing an efficient, stable, and easily implemented solution for the study of large deformation dynamics of incompressible materials. Attached Figure Description

[0159] Figure 1 This is a schematic diagram of the continuum deformation of the present invention;

[0160] Figure 2 This is a schematic diagram of a discrete continuum in the material point calculation framework of the present invention;

[0161] Figure 3 This is a schematic diagram of matter point configuration update under the total Lagrange matter point calculation framework of the present invention; where (a) is the initial configuration, (b) is the current configuration, and (c) is the current configuration on the initial configuration;

[0162] □ represents background mesh nodes; ● represents material points;

[0163] Figure 4 This is a flowchart illustrating the operation of the Explicit Material Point Method (STLMPM) for dynamic analysis of large deformation of weakly compressible materials according to the present invention.

[0164] Figure 5 This is a schematic diagram of the structure and boundary conditions of Example 1 of the two-dimensional incompressible large deformation dynamic analysis of the present invention, and the displacement time history curve of point A in the vertical direction when using 9 grid points as interpolation nodes; wherein, (a) is a schematic diagram of the geometry and boundary conditions, and (b) is the displacement time history curve of point A;

[0165] Figure 6 This is a convergence diagram of the displacement of point A at t=7s obtained by eight numerical methods in Embodiment 1 of the present invention as the mesh density increases; wherein, (a) is the displacement obtained by simulation using a smooth kernel function, and (b) is the displacement obtained by simulation using a B-spline function;

[0166] Figure 7The pressure contour maps obtained by eight numerical methods in Embodiment 1 of the present invention are shown, wherein (a)-(d) are smooth kernel function interpolation methods, and (e)-(h) are displacement contour maps obtained by B-spline interpolation function methods.

[0167] Figure 8 This is a schematic diagram of the second embodiment of the three-dimensional incompressible large deformation dynamic analysis of the present invention, where (a) is the structural boundary condition and (b) is a schematic diagram of the tetrahedral background mesh division.

[0168] Figure 9 This is a schematic diagram showing the convergence of the time-history displacement curves of point A in two mesh divisions of Embodiment 2 of the present invention from 0 to 0.5 s with increasing mesh density; wherein, (a) is a tetrahedral mesh division and (b) is a hexahedral mesh division;

[0169] Figure 10 This is a schematic diagram comparing the time-history displacement curve of point A from 0 to 2 s with the results of Kadapa et al. when using tetrahedral meshing in Embodiment 2 of the present invention.

[0170] Figure 11 The pressure field cloud diagrams for Embodiment 2 of the present invention at times of 0.075s, 0.2375s, 0.425s, 1.125s, 1.5s and 2s are shown; where (a)-(f) are tetrahedral meshes and (m)-(r) are hexahedral meshes. Detailed Implementation

[0171] 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.

[0172] To make the objectives, technical solutions, and specific implementation effects of this invention clearer, three specific embodiments are described below in conjunction with the appendix. Figure 5-11 The accuracy, reliability, and superior performance of the STLMPM method proposed in this invention will be further described in detail.

[0173] First, we will conduct a two-dimensional dynamic analysis of large deformation in soft materials to demonstrate the accuracy of the developed incompressible treatment technique. Furthermore, we will further illustrate the reliability and superior performance of the proposed STLMPM method for weakly compressible large deformation problems by conducting a three-dimensional dynamic analysis of large deformation in soft materials. The reference solutions in the two examples are from the work of Zhang et al. and Kadapa et al., respectively. Gravity effects are not considered in any of the examples.

[0174] (1) Example 1: Kinetic response of Cook membrane (with appendix) Figures 5-7 )

[0175] Example 1 investigated the dynamic Cook membrane problem involving large deformation. Based on simulations by Castanar et al., an approximately incompressible Neo-Hookean material model was selected. The parameters for the large deformation dynamic problem were set as follows: shear modulus μ0 = 83.4 Pa, Poisson's ratio υ = 0.499, and density ρ0 = 1 kg / m³. 3 The left side of the model is completely fixed, while a tensile force of t0 = 6.25 Pa is applied to the right side. Geometric parameters and boundary conditions are as follows: Figure 5 As shown in (a). To analyze the accuracy and convergence of the incompressible processing method, STLMPM and TLMPM were used to model the following four different cases:

[0176] Case 1: Using incompressible correction algorithm 2, each particle is configured with 9 activation nodes (for TLMPM, i.e., using B-spline function).

[0177] Case 2: Using the incompressible correction algorithm 2, each particle is configured with 4 activation nodes (for TLMPM, i.e., using linear shape functions).

[0178] Case 3: Without using the incompressible correction algorithm 2, each particle is configured with 9 activation nodes.

[0179] Case 4: Without using the incompressible correction algorithm 2, each particle is configured with 4 activation nodes.

[0180] It should be noted that square grids are used in all cases in this numerical example, as the efficiency of the proposed incompressible correction method is only being verified. The time increment is set to 1e-5 s. Figure 5 (b) shows a comparison between the time history curve of the vertical displacement of point A obtained using STLMPM in case 1 when N=44 and the results of Zhang et al. Additionally, Figure 6 (a) shows the convergence behavior of the displacement of point A in the X2 direction at t = 7s using STLMPM for four different cases with different N (i.e., the number of particles along the left edge). Figure 6 (b) shows the simulation results of TLMPM. The results indicate that the displacement in the X2 direction of point A tends to converge as the mesh is gradually refined (i.e., N gradually increases), but the convergence rate is slower for cases 3 and 4. The results show that the results using the incompressible formula are similar, while the results without the incompressible formula differ significantly. Furthermore, Figure 7 The pressure (p = κ0(J-1)) distributions for four cases at t = 7s are presented, with case 4 omitted due to excessive pressure fluctuations. It can be seen that the two cases using the incompressible correction method exhibit smoother pressure distributions, while the other two cases show larger pressure fluctuations. This demonstrates that the proposed method can effectively solve the volumetric locking problem in large deformation dynamics.

[0181] (2) Example 2: Three-dimensional cylinder torsion (with appendix) Figures 8-11 )

[0182] Example 2 considers the following: Figure 8 Simulation of the torsional behavior of the approximately incompressible cylinder shown. The geometric model and boundary conditions of the problem are as follows: Figure 8 As shown in (a), fixed boundary conditions are applied to the bottom surface of the model, and the origin of the coordinate system is set at the center of the bottom surface. The initial velocity of the cylinder is position-dependent, and its expression is: The parameters were set as follows: shear modulus μ0 = 5.67 MPa, Poisson's ratio υ = 0.499, and density ρ0 = 1100 kg / m³. 3 The model is simulated using STLMPM. The tetrahedral background mesh used for analysis is as follows: Figure 8 As shown in (b), the hexahedral mesh is represented by a cube. For the tetrahedral and hexahedral meshes, each particle is assigned 5 and 8 active nodes, respectively. During the generation of the tetrahedral mesh, a node is added at the center of the hexahedral mesh, subdividing it into 12 tetrahedra. Figure 9 The displacement in the X3 direction at point A is shown under different mesh sizes, and the results of Kadapa et al. are also presented for comparison. The results show that as the mesh size decreases, the displacement-time curves obtained from both mesh types converge rapidly within the first 0.5 seconds and agree very well with the reference solution. With further mesh refinement, the long-term response of the displacement-time curve at point A gradually converges, as shown... Figure 10 As shown, the time-history displacement curves of point A within the range of 0–2 s are presented when using a tetrahedral mesh. Furthermore, when the mesh size of the tetrahedral mesh is h = 0.0625 m and the mesh size of the hexahedral mesh is h = 0.05 m, the pressure distribution of the cylinder at times 0.075 s, 0.2375 s, 0.425 s, 1.125 s, 1.5 s, and 2 s are shown below. Figure 11 As shown. Among them, Figure 11 (af) represents the calculation result of the tetrahedral mesh. Figure 11 (mr) represents the calculation result for the hexahedral mesh. The calculation results show that the pressure distribution obtained from both tetrahedral and hexahedral meshes exhibits good consistency during the deformation of the cylinder. Furthermore, the pressure distribution demonstrates excellent smoothness, validating the effectiveness of the incompressible formula proposed in this study.

[0183] In summary, the multi-level verifications of the two embodiments above comprehensively demonstrate the accuracy and effectiveness of the Explicit Material Point Method (STLMPM) for dynamic analysis of large deformation in weakly compressible materials proposed in this invention. These verifications not only highlight the significant advantages of this method in terms of computational efficiency and scale, but also prove the necessity for its further development. Furthermore, the results show that the explicit material point method developed in this invention can be efficiently coupled with other interpolation methods and can efficiently meet the solution requirements for the dynamic response of large deformation in soft materials. Therefore, the Explicit Material Point Method (STLMPM) for dynamic analysis of large deformation in weakly compressible materials proposed in this invention is an innovative numerical computation tool that combines high performance with broad development prospects.

[0184] 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. An explicit material point method for dynamic analysis of large deformation of weakly compressible material, characterized by, The specific steps are as follows: Step one, define Euler background mesh, the initial coordinates of Euler background mesh nodes are marked as ; establish discrete material point model and define physical material parameters of discrete material points, initialize material point variables including material point initial coordinates , initial material point displacement vector , initial material point velocity vector , material point mass , initial material point volume , initial material point deformation gradient tensor , material point shear modulus and Poisson's ratio ; Step two, the mapping relationship between the discrete material points and the Euler background grid is established by the smooth kernel function, the physical material parameters on the discrete material points are mapped to the Euler background grid nodes, the total Lagrangian calculation framework, the interpolation function and its gradient are calculated on the initial configuration, and the Gaussian kernel function is used to construct the interpolation function of the material points to the Euler background grid; Firstly, the weight of material point on the Euler background grid node is calculated by the Gaussian kernel function, and is respectively recorded as The expression of the Gaussian kernel function is ； wherein is the "unit radial distance", is 2.0 to 2.5; then, the unit radial distance and the node weight is: ； in, express The 2-norm, Represents the radius of the supporting domain of a matter point, subscript and The background mesh nodes and material points are represented respectively; the interpolation function is obtained by reconstructing the data using weighted least squares. and its spatial gradient ,in, The spatial gradient operator represents the initial configuration, and the specific reconstruction method is as follows: ； wherein, is a spatial coordinate in the initial configuration, is a scalar; is a vector of length dim, dim denoting the dimension of the problem; is a dim x dim matrix, which is: ； wherein, is the nuclear basis function vector; the linear kernel basis function expression in three dimensions is , is the length in three dimensions under the initial configuration, is the total number of Euler background grid nodes; Step three, calculate the explicit material point method time step; When performing numerical simulation of explicit dynamics of weakly compressible materials, the time step The CFL condition needs to be satisfied to ensure the stability of the results, i.e. to make the length of stress wave propagation not exceed 1 grid within a single step simulation, thereby determining the maximum time step of the simulation; meanwhile, considering the dynamic effect, the time step in the CFL condition is: ； wherein, is the initial density of the material point; is the bulk modulus of the material point; is the shear modulus of the material point; is a scalar less than 1 as a scaling parameter for the time increment; is the Euler background mesh size defined in step one; denotes the velocity of the material point at the current time; is the material sound speed; Step four, map the information on the material points to the Eulerian background grid by the interpolation function and its gradient, calculate the variables needed to solve the governing equations; specifically, calculate the node momentum and velocity by the mass and velocity of the material points , calculate the node internal force by the stress of the material points , and calculate the node external force by the boundary conditions ; Step five, apply Dirichlet boundary conditions, set the current time grid node momentum and node internal force of the grid nodes to which the boundary conditions are applied to 0; Step six, based on the internal forces of the background mesh nodes. ,external force and quality Calculate the current velocity of the node The velocity at the next time step is calculated using the Euler forward difference scheme. In the material point method, the velocity update format is divided into FLIP format and PIC format. FLIP format is used for velocity updates, mapping the background mesh node velocities back to the material points to update the material point velocities. And update the position of the matter point through the node velocity. and deformation gradient tensor ; Step seven, perform incompressible processing, average the volume part of the deformation gradient tensor of the material points in combination with the F-bar method of the incremental total quantity format, and bring the averaged deformation gradient tensor into the constitutive equation to update the stress of the material points; First, the Jacobian on the material points is mapped to the Euler background grid nodes according to the initial volume to obtain the Euler background grid node Jacobian, which is specifically written as: ； in, Representing a point of matter The initial volume, Represents the total number of matter points. For control parameters, It is the Jacobian of the average deformation gradient tensor at the current moment, and its initial value for each material point is 1. It is the Jacobian of the deformation gradient increment at the current time step, specifically expressed as follows: ； wherein denotes the average deformation gradient tensor of the material point at the instant The Jacobian on the Euler background grid nodes is mapped back to the material points to calculate the average Jacobian on the material points, which is specifically written as: ； Then, is calculated as: ； By replacing the deformation gradient tensor at time t by the average deformation gradient tensor at time t, the PK1 stress of the material point is updated as ​ Step eight, store and output relevant variable information, return to step three, enter the next time step, and continue until the calculation is completed.

2. The weakly compressible material large deformation dynamic analysis explicit material point method according to claim 1, wherein The smooth kernel function constructs the interpolation function, which is suitable for background grids of any shape and background grid nodes arranged according to any rule; The smooth kernel function adopts a Gaussian kernel function: in the material point method, an arbitrary solid domain is discretized into a series of material points, and a set of background grids are matched to solve the discrete form of the control equation. When the solid domain boundary does not coincide with the background grid, the calculation result accuracy will be reduced due to the initial geometric modeling error and inaccurate boundary condition application; the Gaussian kernel function is used to distribute the grid nodes for physical information mapping to the material points and calculate the initial weight And using weighted least squares to Reconstruction to get interpolation function And its spatial gradient ; For a two-dimensional problem, at least 3 background grid nodes are provided for each material point as interpolation nodes, and for a three-dimensional problem, at least 4 grid nodes are provided for each material point as interpolation nodes.

3. The weakly compressible material large deformation dynamic analysis explicit material point method according to claim 2, wherein In the incompressible processing, the total quantity average format and the incremental average format contain the volume part: in the F-bar method, the volume part of the deformation gradient tensor is averaged, the volume part of the deformation gradient tensor is weighted and mapped to the background grid under the material point method calculation framework, and then the average deformation gradient tensor is calculated by mapping back from the background grid to the material points; Full-mass averaging format, the Jacobian of the material point is mapped to the background grid to get the nodal Jacobian The specific operation is as follows: ； Wherein, the subscript tot represents the node Jacobian obtained by using the total quantity average format; Incremental average format, deformation gradient increment and Instantaneous deformation gradient tensor The relationship between them and their corresponding average forms is written as: ； Thus, and The Jacobian of the representation of the system is given by: ； According to formula (17) and (18), the following is obtained ； Then, the incremental format node Jacobian is defined as: ； By controlling parameters In conjunction with the full and delta formats, the Jacobian computation for background mesh nodes is: ； The average Jacobian of the material points is computed by equation (14); then, the average deformation gradient tensor is computed by equation (15) .

4. The weakly compressible material large deformation dynamic analysis explicit material point method of claim 3, wherein, The specific implementation process of the interpolation function construction is as follows: The specific implementation process of the interpolation function construction is as follows: Step 1, Set the support domain radius for each material point ; Step 2, Set initial coordinates of material points and initial coordinates of Euler background grid nodes ; Step 3, Set parameters and The weight of each grid node within the support domain of each material point is calculated using equation (1) and equation (2) ; Step 4, Calculate the matrix according to formula (6) and (7) and calculate its inverse matrix ; Step 5, calculate according to formulas (5) and (7) , and ; Step 6, Calculate the interpolation function and the interpolation function gradient according to equations (3) and (4), respectively.

5. The weakly compressible material large deformation dynamic analysis explicit material point method of claim 3, wherein, The specific implementation process of the weakly compressible is as follows: Step 1, compute the increment of deformation gradient tensor according to equations (11) and (12) and its Jacobian ; Step 2, compute the background grid node Jacobian matrix according to equations (16), (20) and (21) ; Step 3. Calculate the average deformation gradient tensor of the material point according to formula (15) .

Images