Explicit substance point method for analyzing large deformation power of weak compressible substance

Abstract

The invention belongs to the technical field of weak compressible large deformation of soft materials, provides a dynamic analysis explicit substance point method for large deformation of weak compressible materials, and provides a brand new numerical calculation method for large deformation research of the weak compressible materials. According to the method, a weak-form control equation is discretized in a total Lagrangian substance point method framework, and challenges caused by grid distortion when a traditional method is used for solving a large deformation problem are avoided. Besides, the invention provides a volume self-locking processing technology based on an F-bar method, the problem of volume locking of the weakly compressible substance in numerical simulation can be simply, conveniently and effectively solved, and the volume self-locking processing technology can be conveniently embedded into an existing substance point method program. Meanwhile, a construction method of a smooth kernel interpolation function is developed, and the problem that the solving precision of a traditional substance point method in a boundary region is reduced can be effectively solved. According to the method, the problem of large deformation of the weak compressible material can be solved in a stable and efficient manner, and a new technical path is provided for research and application in related fields.

Description

Technical Field

[0001] The present invention belongs to the technical field of large deformation analysis of soft matter, and particularly relates to an explicit material point method for large deformation dynamic analysis of weakly compressible materials. Background Art

[0002] Soft matter widely exists 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. The numerical modeling of these materials faces significant challenges due to their ability to undergo extremely large recoverable elastic deformations. Traditional mesh-based methods (such as the finite element method) solve problems by discretizing the problem domain into a finite number of small elements, where the accuracy of the solution highly depends on the mesh quality. However, under large deformation conditions, severe mesh distortion reduces the computational accuracy and requires the application of computationally expensive remeshing techniques. To overcome these challenges, various meshless methods have been developed, including smoothed particle hydrodynamics, peridynamics, and the material point method. Among these methods, the material point method shows great potential because it effectively combines the advantages of the Lagrangian method and the Eulerian method. In the material point method, the physical domain is discretized by placing material points (particles) within a background mesh. These Lagrangian particles move on a fixed Eulerian mesh and the governing equations are efficiently solved within it.

[0003] As an effective large deformation modeling method, the Material Point Method (MPM) was first proposed by Sulsky et al. and has been widely verified. Standard MPM usually uses linear shape functions to transfer physical properties between particles and background grid nodes. However, when particles cross the background grid (i.e., the element crossing phenomenon), so-called "cross-grid noise" will occur due to discontinuous derivatives at the background grid boundaries, which may even lead to energy oscillation problems or energy non-conservation, and further result in simulation failures. To solve this problem and improve the accuracy of large deformation simulations, a variety of 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, etc. It should be noted that the above interpolation methods are usually constructed based on regular background grids, but the accuracy will be reduced when dealing with the problem that the geometric boundary does not overlap with the background grid. In addition, designing a background grid adapted to a specific problem geometry is somewhat challenging because the movement of Lagrangian particles within the grid may invalidate the carefully designed grid. This invalidation increases the complexity of maintaining accuracy and stability in dynamic scenarios where particle positions change significantly over time. With the development of the Total Lagrangian Material Point Method (TLMPM), this difficult problem is expected to be solved. Different from standard MPM, TLMPM solves the governing equations in the initial configuration, which means that particles no longer cross the grid. And this method has been successfully applied to the simulation of large deformation problems and achieved good results. In addition, TLMPM only needs to construct the shape function once during the calculation process, thus improving the computational efficiency of MPM.

[0004] When dealing with nearly incompressible materials, the solution framework of the Material Point Method (MPM) is prone to volumetric locking, i.e., it exhibits overly rigid behavior, leading to errors in the calculation of strain and stress fields. This problem also exists in the Finite Element Method (FEM) and related methods, because the overly strong incompressibility constraint at the integration points will reduce the accuracy of element stiffness calculation. In the field of finite elements, to alleviate the locking phenomenon, a variety of methods have been proposed to relax the kinematic incompressibility constraint of elements in different ways. Among them, only those methods compatible with finite deformation kinematics and movable integration points can be adapted to MPM. To solve the volumetric locking problem, a series of techniques and strategies have been proposed within the MPM framework, such as the mixed displacement-pressure (u-p) formulation and the mixed displacement-pressure-Jacobi (u-p-J) formulation. However, the above locking mitigation methods usually require significant modifications to the existing MPM framework and are often optimized for specific basis functions or material types. In addition, the mixed formulation requires changing the standard governing equations, which will significantly increase the cost of implementation and use. Other methods also include operator splitting algorithms and B-bar methods, but their existing versions usually require a large number of modifications to the time stepping scheme or the basis functions for interacting with adjacent elements when applied in MPM. In contrast, the F-bar method has received extensive attention due to its simple derivation and has been successfully applied to the updated Lagrangian MPM to alleviate the volumetric locking problem. This volumetric locking treatment strategy based on the F-bar method calculates the assumed deformation using the standard particle-grid transfer scheme in MPM and shows good potential in the dynamic large deformation simulation of weakly incompressible materials.

[0005] Therefore, the present invention will propose a fully Lagrangian explicit material point method for weakly compressible materials to evaluate the dynamic deformation behavior of soft materials under large deformation conditions. Among them, to improve the calculation efficiency and avoid element penetration errors, the fully Lagrangian material point method is used for the discretization and integration of the momentum equation in the physical domain. In addition, 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 grids or arbitrarily arranged nodes. In addition, the F-bar method is further developed within the fully Lagrangian MPM framework to solve the volumetric locking problem. This method not only significantly improves the calculation accuracy of the material point method at the physical domain boundary, but also can efficiently solve the volumetric locking problem with less modification to the existing calculation framework. Summary of the Invention

[0006] For an efficient numerical method for the dynamic analysis of large deformations of soft materials, the present invention innovatively proposes an explicit material point method for the dynamic analysis of large deformations of weakly compressible materials (abbreviation in English: STLMPM). The aims are as follows: First, to improve the computational efficiency of the material point method and solve the numerical oscillation problem that may be caused when material points cross the background grid, the present invention solves the discrete form of the governing equations based on the initial configuration, and simultaneously completes the calculation of the interpolation function and its gradient at the beginning of the calculation; Second, to solve the problem of reduced calculation accuracy of the material point method near the physical domain boundary, a Gaussian kernel function with kernel correction based on the weighted least squares method (WLS) is used to establish the interpolation formula between particles and the background grid; In addition, to alleviate the volume locking problem caused by the nearly incompressible nature of soft materials, the F-bar method is further developed within the framework of TLMPM to relieve the over-stiff displacement field and stress field caused by material incompressibility; Finally, the present invention aims to overcome the limitation that the existing finite element method has a serious decline in numerical accuracy due to mesh distortion when analyzing large deformation problems of soft materials, and at the same time make up for the deficiencies of the existing explicit material point method in terms of low computational efficiency when analyzing weakly compressible problems and the need for a large number of modifications to the constitutive equation or the original calculation framework.

[0007] The technical solution of the present invention: An explicit material point method for the dynamic analysis of large deformations of weakly compressible materials, the specific steps are as follows:

[0008] Step 1, define the Eulerian background grid, and denote the initial coordinates of the Eulerian background grid nodes as X I ; establish a discrete material point model and define the physical material parameters of the discrete material points, and initialize the material point variables including the initial material point position vector X p , the initial material point displacement vector u p , the initial material point velocity vector the material point mass m p , the initial material point volume V p 0 , the initial material point deformation gradient tensor F p , the material point shear modulus μ and Poisson's ratio υ;

[0009] Step 2, establish the mapping relationship between the discrete material points and the Eulerian background grid through the smooth kernel function, and map the physical material parameters on the discrete material points to the Eulerian background grid nodes;

[0010] Adopt the total Lagrangian calculation framework, calculate the interpolation function and its gradient on the initial configuration (undeformed configuration), and denote them as and where Denote the spatial gradient operator on the initial configuration. The subscript Ip represents the values of the interpolation function and its gradient established at the material point p at the grid node I. The specific calculation method is as follows:

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

[0012]

[0013] where, in the present invention, the parameter is selected as k = 1, r is the "unit radial distance", and α is 2.0 - 2.5;

[0014] The unit radial distance r and the nodal weight are calculated as:

[0015]

[0016] where, ‖■‖ represents the 2-norm of ■, R p represents the radius of the support domain of the material point, and the subscripts I and p represent the physical quantities at the background grid node and the material point respectively;

[0017] According to the chain rule of differentiation, the spatial gradient of the nodal weight on the initial configuration is written as:

[0018]

[0019] Reconstruct the interpolation function by the weighted least squares method; the reconstructed interpolation function and its spatial gradient on the initial configuration are respectively denoted as and are expressed as:

[0020]

[0021] where, the correction matrix C(X, X I ) uses the weighted least squares scheme and is written as:

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

[0023]

[0024] where n I is the total number of background grid nodes; then the inverse matrix 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:

[0029]

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

[0031] When performing explicit dynamic numerical simulation of weakly compressible materials, the time step Δt needs to satisfy the CFL condition to ensure the stability of the results, that is, to ensure that within a single-step simulation, the propagation length of the stress wave does not exceed 1 grid, from which the maximum time step of the simulation can be determined; at the same time, the dynamic effects also need to be considered;

[0032] First, the material wave speed c in the CFL condition can be calculated as:

[0033]

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

[0035]

[0036] where C r is a scalar less than 1, serving as a scaling parameter for the time increment; h g is the Eulerian background grid size defined in Step 1; represents the velocity of the material point at the current moment; c is the material sound speed;

[0037] Step 4: Map the information at the material points to the background grid through the interpolation function and its gradient, and calculate the nodal momentum The nodal mass m I , the nodal internal force and calculate the nodal external force through the boundary conditions

[0038] Map the momentum and mass of the material point at the current moment to the background grid, which is specifically written as:

[0039]

[0040] Among them, represents the initial density of the material point, and V p 0 represents the initial volume of the material point, and n p represents the total number of material points, represents the momentum of the grid node, represents the velocity of the grid node at the current moment, and m I represents the mass of the grid node, and its specific calculation formula is written as:

[0041]

[0042] Among them, m p represents the mass of material point p;

[0043] The internal force of the node is calculated through the stress information carried on the material point at the current moment, and is calculated by the first Piola-Kirchhoff stress as:

[0044]

[0045] Among them, represents the first Piola-Kirchhoff stress on the material point at the current moment; According to the applied external force boundary conditions, the external force of the grid node is calculated as:

[0046]

[0047] Among them, represents the force per unit volume, represents the force per unit area, is the thickness of the virtual boundary layer, which is used for the force application boundary conditions in two-dimensional problems;

[0048] Step 5, apply the Dirichlet boundary condition, and set the current moment grid node momentum and node internal force of the grid node at the boundary where the boundary condition is applied to 0;

[0049] Step 6, update the grid node velocity field from the grid node internal force, external force and velocity at the current moment; and update the velocity, position and deformation gradient tensor of the material point through the updated node velocity field, and the specific operations are as follows;

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

[0051]

[0052] Among them, Calculated by ; then the velocity and position of the material point are updated to:

[0053]

[0054] and

[0055]

[0056] Then, according to the basic theory of continuum mechanics, the deformation gradient tensor of the material point is calculated from the velocity of the grid node at the current moment and the current time step:

[0057]

[0058] where represents the deformation gradient tensor of material point p at time t;

[0059] The Jacobian and volume of the material point are respectively updated to:

[0060]

[0061] Step 7, perform incompressible treatment, use the F-bar method combining the incremental total format to average the volume part of the deformation gradient tensor of the material point, and substitute the averaged deformation gradient tensor into the constitutive equation to update the stress of the material point;

[0062] First, map the Jacobian on the material point to the Eulerian background grid node weighted by the initial volume to obtain the Jacobian of the Eulerian background grid node, and the specific form is written as:

[0063]

[0064] where β is a control parameter, is the Jacobian of the average deformation gradient tensor at the current moment, is the Jacobian of the deformation gradient increment of the current time step, and the specific forms are respectively expressed as:

[0065]

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

[0067] Remap the Jacobian on the Eulerian background grid node back to the material point, calculate the average Jacobian on the material point, and the specific form is written as:

[0068]

[0069] Then, is 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 to replace the deformation gradient F in the strain energy density function. Therefore, the strain energy density function is reconstructed as The first Piola-Kirchhoff stress of the material point is updated to:

[0072]

[0073] Then, the Cauchy stress is calculated as:

[0074]

[0075] where the superscript T represents the transpose of the matrix;

[0076] Step 8: Store and output the information of relevant variables, return to Step 3, enter the next time step, and continue until the calculation is completed.

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

[0078] the smooth kernel function constructs an interpolation function that adapts to background grids of any shape and background grid nodes arranged according to any rule; in the incompressibility treatment, it includes the total average format and the incremental average format of the volume part;

[0079] wherein, 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 used to solve the discretized governing equations. When the boundary of the solid domain does not coincide with the background grid, the calculation result accuracy will decrease due to initial geometric modeling errors and inaccurate application of boundary conditions; the Gaussian kernel function is used to assign grid 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, which is expressed as:

[0080]

[0081] where n I represents the total number of grid nodes; correspondingly, the conditions that its spatial gradient should satisfy are written as:

[0082]

[0083] Using weighted least squares for Reconstructed, and calculated through Equations (8) and (9) to obtain and By using the reconstructed interpolation function The nodal values u of the continuous function u(X) I = u(X I ) is approximated as:

[0084]

[0085] When u(X) = B(X - X p ), according to Equations (4)-(6) and Equation (35), we get:

[0086]

[0087] Therefore, if the linear basis B(X - X p ) is used for the reconstructed interpolation function, it satisfies the normalization condition, i.e., Equation (31), and the linear field reproduction condition, i.e., Equation (32); for two-dimensional problems, at least 3 background grid nodes are equipped for each material point as interpolation nodes, and for three-dimensional problems, at least 4 grid nodes are equipped for each material point as interpolation nodes;

[0088] The incompressible treatment formula mentioned above adopts a mixed format of the total average format and the incremental average format: in the F-bar method, the volumetric part of the deformation gradient tensor is averaged. In the calculation framework of the material point method, the volumetric part of the deformation gradient tensor is weighted and mapped onto the background grid, and then mapped back from the background grid to the material point to calculate the average deformation gradient tensor; first, the deformation gradient tensor is decomposed into a deviatoric part and a volumetric part:

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

[0090] where, F d and F v respectively represent the deviatoric part and the volumetric part of the deformation gradient tensor, and are respectively defined as:

[0091] F d = J -1 / dim F F v = J 1 / dim I (38) where, I represents the identity tensor; through the decomposition of the deformation gradient tensor by Equation (38), the obtained deviatoric part and volumetric part 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 volumetric part with the average volumetric part and is expressed as:

[0094]

[0095] where represents the average volumetric part, represents the Jacobian of the average deformation gradient tensor and is calculated as

[0096] For the total average format, the Jacobian of the material point is mapped to the background grid to obtain the nodal Jacobian The specific operation is as follows:

[0097]

[0098] where the subscript tot represents the nodal Jacobian obtained using the total average format;

[0099] For the incremental average format, the relationship between the incremental deformation gradient ΔF and F and their corresponding average forms is written as:

[0100]

[0101] Therefore, and The Jacobians of are expressed as:

[0102]

[0103] Substituting Eqs. (42) and (43) into Eq. (40), we get

[0104]

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

[0106]

[0107] By controlling the parameter β and combining the total format and the incremental format, the Jacobian of the background grid nodes is calculated as:

[0108]

[0109] The average Jacobian of the material point is calculated by Eq. (27); then, the average deformation gradient tensor is calculated by Eq. (28)

[0110] The construction of the interpolation function and the specific implementation process of the weakly compressible treatment technology are as follows:

[0111] The specific implementation process of the construction of 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 coordinates X of the material points p and the initial coordinates X of the grid nodes I ;

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

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

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

[0117] Step 6: Calculate the interpolation function and the gradient of the interpolation function

[0118] The specific implementation process of weakly compressible is as follows:

[0119] Step 1: Calculate the increment of the deformation gradient tensor and its Jacobian

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

[0121] Step 3: Calculate the average deformation gradient tensor of the material points according to equation (28)

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

[0123]

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

[0125]

[0126] where, b0 and t0 are the body force and surface force respectively, and v represents the velocity. By taking the variation of the total energy function, the strong form of the governing equations and boundary conditions can be written as

[0127]

[0128] where ü represents the second - order derivative of u with respect to time t, is the boundary of Ω0, and n0 is the unit outer normal vector of the boundary in the initial configuration.

[0129] Then, the governing equation in weak form on the initial configuration can be written as

[0130]

[0131] where w is the weight function of the displacement field. In the MPM calculation framework, the governing equation in integral form of the weak form is further written as the governing equation in summation form of discrete form

[0132]

[0133] where n p is the total number of material points, and the subscripts p and I represent the variables related to the material points and the background grid nodes respectively. Considering the arbitrariness of, and using the lumped mass matrix, the discrete - form equation for each background grid node in the TLMPM framework can be written as

[0134]

[0135] According to the above - mentioned theoretical derivation, a detailed solution format of the total Lagrangian material point method under the large - deformation framework proposed by the present invention is given. It can realize the explicit material point method for large - deformation dynamic analysis of weakly compressible materials by constructing interpolation functions for arbitrary - structure background grids through Gaussian kernel functions, the F - bar incompressibility correction technique, and the forward Euler time - discretization strategy.

[0136] According to Figure 6 the calculation flow chart of the STLMPM method proposed by the present invention as shown, its specific implementation process is as follows:

[0137] 1) Establish a material - point discrete model 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 formula (1) and formula (2) to calculate the weights of the grid nodes within the support domain for each material point

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

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

[0141] 2.4) Initialize the background grid: Interpolate the function and the gradient

[0142] 3) Time step loop while t < tfinal

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

[0144] 3.2) Calculate the nodal momentum mass m I and internal force and external force

[0145] 3.3) Apply the Dirichlet boundary condition;

[0146] 3.4) Update the nodal momentum and calculate the nodal velocity

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

[0148] 3.6) Calculate the increment of the deformation gradient tensor and its Jacobian

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

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

[0151] 3.9) Calculate the average Jacobian of the material point 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 completed;

[0155] where tfinal represents the total simulation time.

[0156] Advantages of the present invention: (1) An explicit material point method for large deformation dynamic analysis of weakly compressible materials provided by the present invention provides a brand-new numerical calculation method for the study of large deformation dynamics of weakly compressible materials. Since this method is equipped with the explicit material point method, compared with traditional grid-based methods, it can effectively overcome the problem of grid distortion. Therefore, its significant advantage lies in being able to handle strong nonlinear problems such as large deformations well. And this method solves the control equations in the initial configuration, and the interpolation function only needs to be calculated once at the beginning of the simulation, greatly improving the calculation efficiency of the traditional material point method. Moreover, this method can be applied to various constitutive models and extended to the incompressible large deformation dynamic analysis of other materials, and can be extended to complex multi-field coupling analysis by embedding multi-physics field coupling theory, such as solid large deformation dynamic fracture analysis;

[0157] (2) An explicit material point method for large deformation dynamic analysis of weakly compressible materials provided by the present invention. In this method, a construction method of an interpolation function modified by a smooth kernel function and weighted least squares is developed. It assigns weights to the background grid nodes within the support domain of each material point through the smooth kernel function, and then uses the weighted least squares method to correct 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 for improving the calculation accuracy of the material point method at the physical domain boundary. And due to being equipped with the total Lagrangian calculation framework, this interpolation function construction method hardly increases the total simulation time and is easy to be nested into the existing material point program. The proposal of this method not only provides a new perspective for the numerical simulation of large deformations of weakly compressible materials in theory, but also shows important application potential in engineering practice.

[0158] (3) The explicit material point method for large deformation dynamic analysis of weakly compressible materials provided by the present invention uses the F-bar method to average the volume part of the material point deformation gradient tensor, effectively solving the volume locking caused by the incompressible property of the material. Compared with processing methods such as multi-field variational methods, it only requires a set of interpolation functions, and the numerical implementation is more convenient without major modifications to the calculation framework of the traditional material point method. In addition, by replacing the original deformation gradient with the averaged deformation gradient tensor, this method can be directly applied to the constitutive model operation, and thus can be easily extended to any constitutive model and adaptive algorithm, showing extremely high adaptability and flexibility. In addition, this incompressibility processing method has good generality and is not only applicable to the material point method, but can also be naturally compatible with a variety of numerical calculation methods, such as the finite element method, the peridynamics method, etc., providing an efficient, stable and easy-to-implement solution for the large deformation dynamics research of incompressible materials. Description of the Drawings

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

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

[0161] Figure 3 It is a schematic diagram of the update of the material point configuration under the total Lagrangian material 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] □ is the background grid node; ● is the material point;

[0163] Figure 4 It is a flowchart of the operation of the explicit material point method (STLMPM) for large deformation dynamic analysis of weakly compressible materials of the present invention;

[0164] Figure 5 It is a schematic diagram of the structure and boundary conditions of Example 1 for two-dimensional incompressible large deformation dynamic analysis of the present invention and the time history curve of the vertical displacement of point A when 9 grid points are used as interpolation nodes; where (a) is a schematic diagram of the geometric shape and boundary conditions, and (b) is the displacement time history curve of point A;

[0165] Figure 6 It is a convergence schematic diagram of the displacement of point A at t = 7s obtained by eight numerical methods in Example 1 of the present invention with the increase of grid density; where (a) is the displacement simulated by using a smooth kernel function, and (b) is the displacement simulated by using a B-spline function;

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

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

[0168] Figure 9 Convergence schematic diagram 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 the increase of mesh density; where (a) is the tetrahedral mesh division, and (b) is the hexahedral mesh division;

[0169] Figure 10 Schematic diagram of the comparison of the time history displacement curve of point A from 0 to 2 s with the results of Kadapa et al. when using tetrahedral mesh division in Embodiment 2 of the present invention;

[0170] Figure 11 Pressure field contour maps at time 0.075 s, 0.2375 s, 0.425 s, 1.125 s, 1.5 s and 2 s in Embodiment 2 of the present invention; where (a)-(f) are tetrahedral meshes, and (m)-(r) are hexahedral meshes. Detailed implementation mode

[0171] The performance of the present invention will be further described in detail below in conjunction with the drawings and embodiments. The following embodiments are used to illustrate the present invention, but cannot be used to limit the scope of application of the present invention.

[0172] In order to make the purpose, technical solution and specific implementation effect of the present invention more clearly displayed, the accuracy, reliability and excellent performance of the STLMPM method proposed by the present invention will be further described in detail through three specific embodiments in conjunction with the attached Figure 5-11 drawings.

[0173] First, we will implement a two-dimensional large deformation dynamics analysis of soft materials to illustrate the accuracy of the developed incompressible treatment technology. In addition, a three-dimensional large deformation dynamics analysis example is carried out to further illustrate the reliability and excellent performance of the STLMPM method proposed by the present invention in the case of weakly compressible large deformation problems. The reference solutions in the two embodiments are the works of Zhang et al. and Kadapa et al. respectively. All embodiments do not consider the gravity effect.

[0174] (1) Embodiment 1: Dynamic response of Cook's membrane (attached Figures 5 to 7 )

[0175] Example 1 studied the dynamic Cook's membrane problem involving large deformations. According to the 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 shear modulus μ0 = 83.4 Pa, Poisson's ratio υ = 0.499, and density ρ0 = 1 kg / m 3 . The left side of the model was fully fixed, and a tensile force of magnitude t0 = 6.25 Pa was applied on the right side. The geometric parameters and boundary conditions are as shown in Figure 5 (a). To analyze the accuracy and convergence of the incompressibility treatment method, STLMPM and TLMPM were used to model the following 4 different cases respectively:

[0176] Case 1: Use the incompressibility correction algorithm 2, with 9 active nodes configured for each particle (for TLMPM, that is, using B-spline shape functions).

[0177] Case 2: Use the incompressibility correction algorithm 2, with 4 active nodes configured for each particle (for TLMPM, that is, using linear shape functions).

[0178] Case 3: Do not use the incompressibility correction algorithm 2, with 9 active nodes configured for each particle.

[0179] Case 4: Do not use the incompressibility correction algorithm 2, with 4 active nodes configured for each particle.

[0180] It should be noted that in this numerical example, square grids were used in all cases because only the efficiency of the proposed incompressibility correction method was verified. The time increment was set to 1e-5 s. Figure 5 (b) shows the comparison of the time history curve of the vertical displacement of point A obtained by STLMPM in Case 1 with the results of Zhang et al. when N = 44. Additionally, Figure 6 (a) shows the convergence behavior of the displacement of point A in the X2 direction at t = 7 s for 4 different cases using STLMPM under different N (i.e., the number of particles along the left side edge). Figure 6 (b) shows the simulation results of TLMPM. The results indicate that as the grid is gradually refined (i.e., N gradually increases), the displacement of point A in the X2 direction tends to converge, but for Case 3 and Case 4, the convergence rate is slower. The results show that the results of the cases using the incompressibility formula are close to each other, while the cases without using the incompressibility formula have larger differences. In addition, Figure 7 The pressure (p = κ0(J - 1)) distributions of the four cases at t = 7 s are given, where the pressure distribution in Case 4 is not shown due to excessive pressure fluctuations. It can be seen that the pressure distributions of the two cases using the incompressibility correction method are smoother, while the pressure distributions of the other two cases have larger fluctuations. This indicates that the proposed method can well solve the volume locking problem in large deformation dynamic problems.

[0181] (2) Example 2: Torsion of a three-dimensional cylinder (Appended Figures 8 to 11 )

[0182] Example 2 considered the simulation of the torsion behavior of an approximately incompressible cylinder as shown in Figure 8 . The geometric model and boundary conditions of the problem are as shown in Figure 8 (a). A fixed boundary condition is 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 has a position dependence, and its expression is: The parameters are set as shear modulus μ0 = 5.67 MPa, Poisson's ratio υ = 0.499, and density ρ0 = 1100 kg / m 3 . This model is simulated by STLMPM. The tetrahedral background mesh for analysis is as shown in Figure 8 (b), and the hexahedral mesh is a cube. For the tetrahedral mesh and the hexahedral mesh, 5 and 8 activation nodes are assigned to each particle respectively. During the generation of the tetrahedral mesh, a node is added at the center of the hexahedral mesh and it is subdivided into 12 tetrahedra. Figure 9 Shows the displacement of point A in the X3 direction under different mesh sizes. The results of Kadapa et al. are also given in the figure for comparison. The results show that as the mesh size decreases, the displacement-time curves obtained by the two mesh types converge rapidly within the first 0.5 s and are in good agreement with the reference solution. As the mesh is further refined, the long-term response of the displacement-time curve at point A gradually converges. As shown in Figure 10 , the time-history displacement curve of point A within 0 - 2 s when using the tetrahedral mesh is given. In addition, when the mesh size of the tetrahedral mesh h = 0.0625 m and the mesh size of the hexahedral mesh h = 0.05 m, the pressure distributions of the cylinder at 0.075 s, 0.2375 s, 0.425 s, 1.125 s, 1.5 s, and 2 s are as shown in Figure 11 . Among them, Figure 11 (a - f) are the calculation results of the tetrahedral mesh, Figure 11 (m - r) are the calculation results of the hexahedral mesh. The calculation results show that during the deformation of the cylinder, the pressure distributions obtained by the tetrahedral mesh and the hexahedral mesh have good consistency. In addition, the pressure distribution shows excellent smoothness, verifying the effectiveness of the incompressible formula proposed in this study.

[0183] In summary, through the multi-level verification of the above two embodiments, the accuracy and effectiveness of the explicit material point method (STLMPM) for large deformation dynamic analysis of weakly compressible materials proposed by the present invention are comprehensively demonstrated. These verifications not only highlight the significant advantages of the method from the perspectives of computational efficiency and computational scale, but also prove the necessity of its further development. In addition, the results show that the explicit material point method developed by the present invention can be efficiently coupled with other interpolation methods and can meet the solution requirements for the large deformation dynamic response of soft materials. Therefore, the explicit material point method (STLMPM) for large deformation dynamic analysis of weakly compressible materials proposed by the present invention is an innovative numerical calculation tool with both high performance and broad development prospects.

[0184] The embodiments of the present invention are given for purposes of illustration and description, and are not exhaustive or limit the invention to the disclosed form. Many modifications and variations are obvious to those of ordinary skill in the art. The embodiments are chosen and described in order to best explain the principles of the invention and its practical application, and to enable those of ordinary skill in the art to understand the invention and design various embodiments with various modifications suitable for specific purposes.

Claims

1. An explicit material point method for large deformation dynamic analysis of weakly compressible materials, characterized in that: The specific steps are as follows: Step 1: Define the Euler background grid. The initial coordinates of the Euler background grid nodes are marked as 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 Material point mass m p , initial material point volume Initial material point deformation gradient tensor F p , material point shear modulus μ and Poisson's ratio υ; Step 2: Establish a mapping relationship between discrete material points and the Euler background grid through a smooth kernel function, and map the physical material parameters on the discrete material points to the Euler background grid nodes; use the total Lagrangian calculation framework, interpolation function and its gradient to perform calculations on the initial configuration; use the Gaussian kernel function to construct the interpolation function from the material point to the Euler background grid; First, the weights of the material points on the Euler background grid nodes are calculated by the Gaussian kernel function and are recorded as The expression of Gaussian kernel function is: Among them, the parameter selection is k=1, r is the "unit radial distance", and α is 2.0~2.5; then, the unit radial distance r and the node weight for: Among them, ‖■‖ represents the second norm of ■, R p Represents the radius of the material point support domain, and the subscripts I and p represent the physical quantities on the background grid nodes and material points respectively; the weighted least squares method is used for reconstruction to obtain the interpolation function Its spatial gradient in, represents the spatial gradient operator on the initial configuration. The specific reconstruction method is as follows: Among them, X is the spatial coordinate of the initial configuration, C1 is a scalar; C2 is a vector of length dim, dim represents the dimension of the problem; C3 is a dim×dim matrix, which is: C(X,X I )=B(X-X p )M(X p ) -1 B T (X I -X p ) (5) Where B is the kernel basis function vector; the three-dimensional linear kernel basis function expression is B(X) = [1, X1, X2, X3], X1, X2, X3 are the lengths in three dimensions under the initial configuration, n I is the total number of nodes of the Euler background grid; Step 3, calculate the time step of the explicit material point method; When performing numerical simulation of the explicit dynamics of weakly compressible materials, the time step Δt needs to meet the CFL condition to ensure the stability of the results, that is, the length of stress wave propagation should not exceed 1 grid in a single-step simulation, thereby determining the maximum time step of the simulation; at the same time, considering the dynamic effect, the time step that meets the CFL condition is: Where ρ is the initial density of the material point; K is the bulk modulus of the material point; μ is the shear modulus of the material; C r is a scalar less than 1, serving as a scaling parameter for the time increment; h g is the size of the Euler background grid defined in step 1; represents the velocity of the material point at the current moment; c is the material sound speed; Step 4: Map the information on the material point to the Euler background grid through the interpolation function and its gradient, and calculate the variables required to solve the control equation; specifically, the mass m of the material point p and speed Calculating node momentum The mass of the material point m p Calculate the node mass m I , stress at a material point Node internal forces And calculate the node external forces through boundary conditions Step 5, applying Dirichlet boundary conditions, and setting the current moment mesh node momentum and node internal force of the mesh node where the boundary conditions are applied to 0; Step 6: Internal force of background mesh nodes external force and mass m I Calculate the current speed of the node The Euler forward difference scheme is used to calculate the velocity at the next time step. The velocity update formats in the material point method are divided into FLIP format and PIC format. The FLIP format is used for velocity update, and the background grid node velocity is mapped back to the material point to update the material point velocity. And update the material point position through the node velocity and the deformation gradient tensor In step seven, incompressible processing is performed, the volume part of the deformation gradient tensor of the material point is averaged using the invented F-bar method combined with the incremental full-volume format, and the averaged deformation gradient tensor is brought into the constitutive equation to update the material point stress; First, the Jacobian on the material point is mapped to the Euler background grid node according to the initial volume weighting to obtain the Jacobian of the Euler background grid node, which is written as: Among them, n p represents the total number of material points, β is the control parameter, is the Jacobian of the average deformation gradient tensor at the current moment, and its initialization value for each material point is 1. is the Jacobian of the deformation gradient increment at the current time step, and the specific forms are expressed as: in, represents the average deformation gradient tensor of material point p at time t; Remap the Jacobian on the Euler background grid node back to the material point and calculate the average Jacobian on the material point. The specific form is written as: Then, is calculated as: By replacing the deformation gradient tensor F in the material strain energy density function with the average deformation gradient tensor Update PK1 stress of material points; Step 8: Store and output relevant variable information, return to step 3, and enter the next time step until the calculation is completed.

2. The explicit material point method for large deformation dynamic analysis of weakly compressible materials according to claim 1 is characterized in that: The smooth kernel function constructs an interpolation function to adapt to background grids of arbitrary shapes and background grid nodes arranged according to arbitrary rules; The smooth kernel function adopts the Gaussian kernel function: in the material point method, any solid domain is discretized into a series of material points, and a set of background grids are used to solve the discrete form of the control equation. When the solid domain boundary does not coincide with the background grid, the accuracy of the calculation result will decrease due to the initial geometric modeling error and the inaccurate application of the boundary conditions; the Gaussian kernel function is used to assign grid nodes for physical information mapping to the material points and calculate the initial weights And use weighted least squares Reconstruct the interpolation function Its spatial gradient For two-dimensional problems, each material point is equipped with at least three background grid nodes as interpolation nodes. For three-dimensional problems, each material point is equipped with at least four grid nodes as interpolation nodes.

3. The explicit material point method for large deformation dynamic analysis of weakly compressible materials according to claim 2 is characterized in that: The incompressible processing includes the full average format and the incremental average format of the volume part: in the F-bar method, the volume part of the deformation gradient tensor is averaged, and in the material point method calculation framework, the volume part of the deformation gradient tensor is weighted mapped to the background grid, and then mapped from the background grid back to the material point to calculate the average deformation gradient tensor; The full-volume average format maps the Jacobian of the material point to the background grid to obtain the node Jacobian The specific operations are: Wherein, the subscript tot represents the node Jacobian obtained using the full-average format; Incremental average format, the relationship between the deformation gradient increments ΔF and F and their corresponding average forms is written as: therefore, and The Jacobian representation of is: According to equations (17) and (18), we can obtain Then, the node Jacobian in incremental format is defined as: By controlling the parameter β and combining the full format and incremental format, the Jacobian calculation of the background grid node is: The average Jacobian of the material point is calculated by equation (14); then, the average deformation gradient tensor is calculated by equation (15):

4. The explicit material point method for large deformation dynamic analysis of weakly compressible materials according to claim 3 is characterized in that: The specific implementation processes of the interpolation function construction are as follows: The specific implementation process of interpolation function construction is as follows: Step 1: Set the support radius R for each material point p ; Step 2: Set the initial coordinate X of the material point p and the initial coordinates of the grid nodes X I ; Step 3: Set the parameters k and α, and use formula (1) and formula (2) to calculate the weight of the grid nodes in the support domain for each material point. Step 4: According to formulas (6) and (7), calculate the matrix M and its inverse matrix M -1 ; Step 5: Calculate C1, C2 and C3 according to formulas (5) and (7); Step 6: According to formulas (3) and (4), calculate the interpolation function And the interpolation function gradient 5. The explicit material point method for large deformation dynamic analysis of weakly compressible materials according to claim 3, characterized in that: The specific implementation process of the weak compressibility is as follows: Step 1: Calculate the increment of the deformation gradient tensor according to formulas (11) and (12): Its Jacobian Step 2: Calculate the Jacobian matrix of the background grid node according to formulas (16), (20) and (21): Step 3: Calculate the average deformation gradient tensor of the material point according to formula (15):

Images