Simulation method and application of large deformation process in impact based on meshless interpolation technology

Abstract

The present invention provides a method and application for simulating a large deformation process of an impact based on a meshless interpolation technique, which belongs to the field of computational mechanics and numerical simulation technology, including model parameter setting, discretization of the problem domain and the time domain, iterative updating of the variables to be determined in the problem domain using the regularized Newton method, solving the local deformation gradient at the next moment, solving the node acceleration, updating the node coordinates, updating the material point coordinates and physical information, completing the simulation of the impact process, and the like. The present invention uses the regularized Newton method for the objective function of the local maximum entropy, and by constructing a simple shape near the calculation point, provides a method for quickly, accurately, and efficiently determining the initial point, thereby ensuring the solution efficiency, and using the Nelder‑Mead algorithm to solve the extreme case iteration point, thereby ensuring the stability of the solution, and is used for the complete simulation of the large deformation process of materials under ultra-high-speed impact, and has the characteristics of high calculation accuracy and strong stability.

Description

Technical Field

[0001] The present invention belongs to the technical field of computational mechanics and numerical simulation, and in particular relates to a method and application for simulating a large impact deformation process based on a meshless interpolation technique. Background Art

[0002] In the past few decades, traditional mesh methods have developed rapidly and have been widely used in the study of scientific and engineering problems. However, due to the inherent mesh properties, these methods often encounter difficulties in solving large impact deformation problems. For example, Lagrangian meshes (such as finite element method) are plagued by mesh distortion and cannot solve problems such as large deformation and crack propagation. Euler meshes (such as finite difference method and finite volume method) are difficult to track deformable boundaries and interfaces of moving materials.

[0003] The meshless method does not require mesh discretization, which to a certain extent solves some of the computational difficulties caused by mesh dependence. It is beneficial for solving large deformation processes caused by impact and is widely used in explosion problems.

[0004] At present, the commonly used meshless methods include Smoothed Particle Hydrodynamics (SPH), Element-Free Garlerkin (EFG), Reproducing Kernel Particle Method (RKPM) and Material Point Method (MPM). A large number of researchers have used the above methods and their improved methods to simulate the impact large deformation problem. Although they have achieved rich research results, they still cannot solve the problems of accurately setting displacement boundary conditions, low calculation efficiency, unstable tensile stress and effective calculation of dynamic contact with friction. At the same time, they lack rigorous convergence and error theory analysis. The Optimal Transportation Meshfree (OTM) is an explicit incremental update Lagrangian meshless method proposed for dynamic impact, penetration, explosion and other problems. The OTM method has the following advantages when dealing with impact large deformation problems:

[0005] 1. This method uses the Local Maximum Entropy (LME) interpolation function, which overcomes the essential defect that the interpolation function in the general meshless method (such as the kernel function approximation in SPH and the original moving least squares interpolation in EFG) does not satisfy the Kronecker-delta property, and solves the problem that the traditional meshless method is difficult to accurately impose displacement boundary conditions;

[0006] 2. This method uses a combination of nodes and material points for discretization, and uses material points as integration points to effectively avoid the instability of the calculation results under tensile loads;

[0007] 3. This method adopts a crack propagation algorithm based on the energy criterion and uses the material energy release rate of the material point as the failure criterion. It can simulate the damage of the material well and overcome the shortcomings of the traditional crack propagation numerical method, such as the crack propagation path grid correlation, non-convergence, and inability to clearly characterize the real deformation and failure physical mechanism of the material.

[0008] Compared with traditional grid-based methods, a major difference of meshless methods lies in the processing of interpolation functions. The interpolation function in traditional grid methods is bound to the unit and represented in the form of a specific function, which can be directly obtained by function assignment; while the interpolation function of the meshless method needs to be recalculated by fitting or optimization after each update of the neighborhood particles, and the calculation cost is large. Therefore, in the calculation process of the meshless method, efficient and stable interpolation technology is particularly important. The local maximum entropy interpolation function used by the OTM method needs to be obtained by solving a convex optimization problem, but the use of general algorithms to solve it is inefficient and unstable when simulating large deformation problems. For example, when the material deformation is too large, the number of neighborhood particles decreases, and it is often impossible to obtain a converged shape function value at the material point within a single time step, resulting in non-physical calculation results, which further causes the crash of the calculation program. Therefore, this patent is committed to improving the solution process of the local maximum entropy interpolation function, solving the bottleneck problem of the interpolation function in the OTM method, and providing an efficient and stable simulation strategy to deal with the problem of large impact deformation.

[0009] Local Maximum Entropy Interpolation Function N (x) (in convex optimization problems, the probability p (x) representation) can be obtained by solving a convex optimization problem:

[0010] ,

[0011] Where x is the coordinate of the material point to be interpolated, which should be located in the point set X The convex set conv X Inside, x a ( a =1,…, N ) is the node coordinate, p(x)=( p 1(x) , p 2(x) ,… p N (x) T , is the interpolation function to be found, β is the Pareto optimization parameter, which can control the locality of the shape function and is a given constant.u h (x) can be used with node values u a = u a (x a ) interpolation yields:

[0012] ,

[0013] Using convex optimization theory for rigorous derivation, the dual problem of convex optimization problem is equivalent to:

[0014] ,

[0015] Representative l Minimize a function for a variable The value of This is the objective function of the problem to be solved; This is the unconstrained convex optimization problem that needs to be solved in the end. At this time, the interpolation function p(x) can be expressed as:

[0016] ,

[0017] in, , is the partition function in the sense of statistical mechanics, is the norm operator, and the 2-norm is generally used here. It can be proved that the objective function F ( l )'s Hessian matrix J( l ) is positive definite, so F ( l ) is R d A strictly convex function on , whose minimum exists and is unique.

[0018] Since the objective function is too complex, it is not possible to directly use the extreme value theorem to find an analytical solution. Generally, for this unconstrained convex optimization problem, the Newton method can be used to iteratively solve it. F ( l ). Newton's method is a descent method based on the second-order Taylor expansion to solve the extreme point of the function. This method can find a descent direction that tends to the extreme value at each calculation time. Following this direction, a point closer to the extreme value than the current calculation point can always be obtained. Taking the kth calculation time step as an example, different time steps can be obtained by selecting different search algorithms. t k , combined with the descending direction The calculation point of the step can be calculated , where r and J are the objective functions F ( l )’s gradient and Hessian matrix.

[0019] However, the local convergence of Newton's method inevitably brings two problems:

[0020] First, when the objective function is non-convex, Newton's method often fails, which is reflected in the fact that the determinant of the Hessian matrix J of some points in the space is 0, and it is impossible to find the inverse, so it is impossible to find the Newton direction. Therefore, in convex optimization theory, it is often tested whether the Hessian matrix of the objective function is positive definite;

[0021] Second, the convergence of Newton's method depends on the choice of initial values. In areas far away from the extreme points, Newton's method cannot reach the second-order convergence speed. Inappropriate initial values ​​will cause oscillation or even divergence of the solution. Therefore, in applications, the Lipschitz continuity condition is often used to test the quadratic state of the local points of the objective function.

[0022] The formula to be solved will often encounter the following problems:

[0023] A specific problem of solving the local maximum entropy shape function: Figure 1 As shown, the two-dimensional point set X Take {(0,0),(1,0), (1,1) , (0,1)}, the calculation point x is (0.8,0.2), and the parameter is 1.6, we can draw Two-dimensional surface (attached Figure 2 ), it is not difficult to see that the image is indeed R 2 Furthermore, the surface is a smooth arc near the extreme point, showing obvious quadraticity, while the area far from the extreme point is divided into four obvious planes, that is, when || l When ||→∞, F ( l ) approaches a linear function, then || l ||=0, that is, the Hessian matrix is ​​singular and the determinant is 0 (see Appendix Figure 3 ), there is no inverse matrix, and the Newton direction cannot be directly obtained. More specifically, when the initial value is selected in these regions, or the iterative value falls into these regions, the traditional Newton method will lose its effect.

[0024] To solve the above problems, people have proposed different methods to correct the Hessian matrix, among which the more famous ones are the damped Newton method and the quasi-Newton method. The damped Newton method introduces a damping parameter, also known as a penalty factor, to correct the original Hessian matrix J to J+ m I, ensures the existence of the descending direction. m When it is large, this method actually degenerates Newton's method into the steepest descent method with only first-order convergence rate, and the same constant is used in the whole domain. m, which greatly affects the computational efficiency. The basic idea of ​​the quasi-Newton method is to "approximate the Hessian matrix (or inverse)", making it faster than the steepest descent method, simpler to calculate than the Newton method, and with good overall convergence. Different quasi-Newton methods can be derived from the quasi-Newton conditions, among which the DFP algorithm and the BFGS algorithm are more widely used, but the computational efficiency is still problematic. Summary of the invention

[0025] The purpose of the present invention is to provide a method for simulating large deformation processes of impact based on meshless interpolation technology, so as to solve the problems of failure, low calculation efficiency, etc. when iteratively solving the minimum parameter value of F(λ) by Newton's method in the process of solving the impact problem of local maximum entropy shape function.

[0026] In order to achieve the above object, the technical solution of the present invention is as follows:

[0027] The present invention relates to a method for simulating a large impact deformation process based on a meshless interpolation technique, comprising the following steps:

[0028] S1. Set the impact large deformation problem domain, set the global parameters and model parameters including material model;

[0029] S2. Discretize the problem domain into a set of material points and a set of nodes, initialize the neighborhood node set for each material point, and initialize the physical variables for the nodes. The physical variables include the mass, linear momentum and node force of the node set.

[0030] S3. Discrete time domain, with time increments [ t n , t n+1 ] as a unit cycle to perform S4-S8;

[0031] S4. For each material point, construct the objective function on the neighborhood node set F ( l ), and the regularized Newton method is used to solve the variables in the problem domain. l Perform iterative solution and save the final variables l * , calculate the local maximum entropy interpolation function and its derivative;

[0032] S5. Update the node concentrated mass based on the local maximum entropy interpolation function, update the node coordinates at the next moment, calculate the local deformation gradient increment based on the node coordinates of the next node and the derivative of the maximum entropy interpolation function, and calculate the local deformation gradient at the next moment based on the local deformation gradient increment;

[0033] S6. Substitute the local deformation gradient increment and the local deformation gradient obtained in S5 into the material model to obtain the stress of the material point;

[0034] S7. Integrate the stress and boundary conditions of the material points obtained in S6, update the node force and linear momentum, and determine whether the calculation is to the last time step. If not, go to S8. If yes, go to S9.

[0035] S8. Update the material point coordinates and the neighborhood node set, and return to S4;

[0036] S9. Complete the impact process simulation and output the results.

[0037] Preferably, the regularized Newton method is used in S4 to calculate the variables in the problem domain. l The specific steps for iterative solution are:

[0038] S4.1. Based on the node set and initial values ​​of the physical variables to be determined in the initialization phase l 0, set the accuracy threshold d , let the physical variable to be determined be λ=λ 0;

[0039] S4.2. Calculate the partition function value at each node Z a , the calculation formula is:

[0040] (1),

[0041] In the formula, β is the Pareto optimization parameter, x represents the material point to be calculated, and x a Represents a node set {x a}, ||·|| represents norm operation;

[0042] The partition function values ​​of each node are summed up to get ;

[0043] S4.3. Calculation of shape functions N a , the calculation formula is:

[0044] (2) ;

[0045] S4.4. Objective function F ( l ) is regularized and the corrected objective function is , whose expression is:

[0046] (3),

[0047] in, l is the regularization correction point, that is, the calculation point at the current moment, v is the independent variable of the regularized correction function, ▽F ( l ) represents the gradient vector of the objective function;

[0048] S4.5. Calculate the gradient vector of the objective function ▽ F ( l ), the calculation formula is:

[0049] (4),

[0050] In the formula, N represents the number of nodes, a Indicates the node index number;

[0051] S4.6. Perform norm operation on the gradient vector to obtain accuracy and compare it with the set accuracy threshold d Compare and if the calculated progress is not greater than the accuracy threshold, output the variable at this time as the final variable l * , proceed to S4.9, otherwise proceed to S4.7- S4.8;

[0052] S4.7. Calculate the Hessian matrix H of the objective function ( l ), the regularized Newton descent direction is calculated based on the Hessian matrix and the gradient vector,

[0053] The calculation formula of the Hessian matrix is:

[0054] (5),

[0055] In the formula, r(x, l ) represents the gradient vector, I represents the identity matrix;

[0056] The calculation formula for the regularized Newton descent direction is:

[0057] (6),

[0058] In the formula, rg( l ) represents the direction of regularized Newton descent;

[0059] S4.8. Iteration of the physical variable to be determined , and order , return to step S4.2,

[0060] The iteration formula is:

[0061] (7),

[0062] In the formula, t ( l ) represents the search step size in the direction of regularized Newton descent;

[0063] S4.9. The final variable l * Substitute it into the shape function formula of S4.3 to calculate the local maximum entropy interpolation function. The derivative of the local maximum entropy interpolation function is calculated according to the local maximum entropy interpolation function. The calculation formula of the derivative of the local maximum entropy interpolation function is:

[0064] (8),

[0065] in, (9),

[0066] (10),

[0067] In the formula, r * Represents the gradient vector of the original objective function, J * Represents the Hessian matrix of the original objective function.

[0068] Preferably, the iterative process in S4.8 regularizes the search step length in the Newton descent direction t ( l ) is initially set to 1. After iterating the physical variable to be determined, S4.8 calculates the objective function values ​​before and after the iteration. If the objective function value after the iteration is less than the objective function value before the iteration, it directly returns to step S4.2. If the objective function value after the iteration is not less than the objective function value before the iteration, the optimal step length in the regularized Newton descent direction is calculated by the precise line search algorithm. t ( l ), and iterate again according to the iteration formula, and return to step S4.2 after the re-iteration.

[0069] Preferably, the initial value of the physical variable to be determined in S4.1 is l 0 is obtained by:

[0070] S4.1.1. Determine whether the problem to be solved in the problem domain is a three-dimensional problem or a two-dimensional problem. If it is a three-dimensional problem, select the four nodes closest to the material point x in the node set to form a tetrahedron, and calculate the centroid coordinates of the material point x with respect to the four vertices of the tetrahedron. If it is a two-dimensional problem, select the three points closest to the material point x in the node set to form a triangle, and calculate the centroid coordinates of the material point x with respect to the three vertices of the triangle.

[0071] S4.1.2. Let the shape function of the material point be equal to the components of the corresponding barycentric coordinates, and substitute them into the calculation formula of the shape function to obtain a set of variables l The linear equations are solved to obtain the initial values ​​of the variables. l 0.

[0072] Preferably, the S4.1.2 solves to obtain the initial value of the variable l 0, the partition function values ​​of each node are summed up to get Z ( l 0), determine whether the value is less than the minimum threshold or greater than the maximum threshold. If so, use the Nelder-Mead algorithm to correct the centroid coordinates obtained in S4.1.1. The correction method is as follows:

[0073] S4.1.2.1. l Select one in the space with l 0 associated simple shape, set tolerance e and the maximum number of comparisons NMax, set the initial number of operations Nfunk to 0;

[0074] S4.1.2.2. Sort the simple shape vertices according to the objective function value, mark the worst vertex as W, the second worst vertex as S, and the best vertex as B, and calculate the objective function value of the variable at the worst vertex W F (W) and the objective function value of the variable at the best vertex B F (B) if F (W)- F (B)> e Or if Nfunk>NMax, then output the corrected center of gravity coordinates, otherwise go to S4.1.2.3 to start correction;

[0075] S4.1.2.3. Discard the worst vertex W, find the center M of the face or line formed by the remaining vertices, find the symmetric point R of the worst vertex W about the center M, and increase the number of operations Nfunk by 2;

[0076] S4.1.2.4. Comparison of objective function values ​​of variables at symmetric points F (R) and the objective function value of the variable at the worst vertex W F (W), if F (R)< F (W), then proceed to S4.1.2.5, otherwise proceed to S4.1.2.7;

[0077] S4.1.2.5. Comparison of objective function values ​​of variables at symmetric points F (R) and the objective function value of the variable at the best vertex B F (B) size, if F (R)≤ F (B), then determine point E on the extension line of the center of gravity M and the symmetry point R. The calculation method of point E is E=M+2(RM). Further compare the objective function value of the variable at point E. F (E) Objective function with variables at symmetric points F (R) size, ifF (E)< F (R), replace the worst vertex W with point E, return to S4.1.2.2, if F (E)≥ F (R), replace the worst vertex W with the symmetric point R, and return to S4.1.2.2; if F (R)> F (B), then proceed to S4.1.2.6;

[0078] S4.1.2.6. Objective function for comparing variables at symmetric points F (B) The objective function of the variable at the secondary difference vertex S F (S), if F (R)≤ F (S), replace the worst vertex W with the symmetric point R, reduce the number of operations Nfunk by 1, and return to S4.1.2.2. F (R)> F (S), go to S4.1.2.7;

[0079] S4.1.2.7. Replace the worst vertex W with the symmetric point R and take the midpoint C of the line connecting the centroid M and the symmetric point R;

[0080] S4.1.2.8. Further comparison of the objective function values ​​of the variables at point C F (C) The objective function value of the variable at the worst vertex W F (W), if F (C)> F (W), replace the worst vertex W with point C, return to S4.1.2.2, if F (C)≤ F (W), shrink the simple shape to the best vertex B, output the shrunk vertex, increase the number of operations Nfunk by 3, and return to S4.1.2.2.

[0081] Preferably, the Pareto optimization parameter in S4.2 is calculated by the following formula:

[0082] (11),

[0083] In the formula, h is the local characteristic length of the convex set formed by the node set, taking the minimum node spacing, c is a dimensionless parameter;

[0084] After the neighborhood node set in S8 is updated, the Pareto optimization parameters are recalculated.

[0085] Preferably, the specific steps of S5 include:

[0086] S5.1. Update the node concentrated mass, the formula is:

[0087] (12),

[0088] In the formula,, N a represents the local maximum entropy interpolation function, m represents the concentrated mass at the node, a Indicates the node index number, n express n time, p Indicates the material point index number, N H (x a,n ) represents the node x a The neighborhood material point set of ;

[0089] S5.2. Update the node coordinates of the next node, expressed as:

[0090] (13),

[0091] In the formula, l a,n Indicates that the index number is a The node is n The linear momentum at time, f a,n Indicates that the index number is a The node is n Nodal force at the moment;

[0092] S5.3. Calculate the local deformation gradient increment based on the node coordinates of the next node and the derivative of the maximum entropy interpolation function. The calculation formula is:

[0093] (14),

[0094] In the formula, F p,n→n+1 Indicates from n The local deformation gradient increment from moment to the next moment, N H (x p,n ) represents each material point x p The set of neighboring nodes;

[0095] S5.4. Calculate the local deformation gradient of the material point at the next moment. The calculation formula is:

[0096] (15),

[0097] In the formula, F p,n express n The local deformation gradient at time, F p,n+1 Represents the local deformation gradient at the next moment.

[0098] Preferably, the updated node force in S7 includes node internal force and node external force, and the calculation formula of the node force is:

[0099] (16),

[0100] In the formula, f int a,n+1 Indicates that the index number is a The node internal force at the next moment, f ext a,n+1 Indicates that the index number is a The node external force at the next moment, f a,n+1 Indicates that the index number is a The node force of the node at the next moment;

[0101] The calculation formula for updating the linear momentum in S7 is:

[0102] (17),

[0103] In the formula, l a,n Indicates that the index number is a The node is n The linear momentum at time, l a,n+1 Indicates that the index number is a The linear momentum of the node at the next moment.

[0104] Preferably, the formula for updating the material point coordinates in S8 is:

[0105] (18),

[0106] In the formula, x a,n+1 Indicates the node coordinates at the next moment, x p,n+1 represents the coordinates of the material point at the next moment, N H (x p,n ) represents the material point x p The set of neighboring nodes;

[0107] Update the neighboring node set at the next moment N H (x p,n+1 ), which is obtained by the following method: given a search radius, determine the nodes within the search radius of each material point as the neighborhood nodes of the material point, and input the neighborhood nodes as the neighborhood node set. Traverse the neighborhood node set of the material point N H (x p,n+1 ), construct the neighborhood material point set of the node N H (x a,n+1 ).

[0108] The invention also relates to the application of a method for simulating a large impact deformation process based on a gridless interpolation technique, which is applied to simulating a large impact deformation process of a bullet impacting a target plate.

[0109] Compared with the prior art, the technical solution provided by the present invention has the following beneficial effects:

[0110] 1. The method for simulating large deformation process of impact based on meshless interpolation technology of the present invention adopts regularized Newton method to calculate the physical variables to be determined in the problem domain. l Iterative updates are performed to repair the area where the determinant of the Hessian matrix of the original objective function has a value of zero. After the repair, the regularized Newton direction can be constructed at any point in the space for solution without causing the Newton method to lose its effect.

[0111] 2. The method for simulating large deformation process of impact based on meshless interpolation technology of the present invention adopts regularized Newton method to solve the physical variables in the problem domain. l When performing iterative updates, for three-dimensional problems, select the four nodes closest to the material point x in the node set to form a tetrahedron, and calculate the centroid coordinates of the material point x about the four vertices of the tetrahedron. For two-dimensional problems, select the three points closest to the material point x in the node set to form a triangle, and calculate the centroid coordinates of the material point x about the three vertices of the triangle. Then let the shape function of the material point be equal to the components of the corresponding central coordinates, and substitute them into the calculation formula of the shape function respectively. After sorting, a set of variables is obtained. l The linear equations are solved to obtain the initial values ​​of the variables. l 0. The initial value obtained by this method is close to the optimal solution and can converge quickly after several iterations, thus increasing the iteration efficiency.

[0112] 3. After the initial value is obtained by the impact large deformation process simulation method based on the meshless interpolation technology involved in the present invention, for certain extreme cases, the Nelder-Mead method is used to shrink the iteration point to the extreme point, ensuring that the calculation will not collapse when the iteration value falls into the extreme case, that is, ensuring the stability of the entire algorithm. BRIEF DESCRIPTION OF THE DRAWINGS

[0113] Figure 1 Schematic diagram of point set and calculation point distribution;

[0114] Figure 2 It is a two-dimensional surface diagram of the objective function (including an enlarged diagram of the extreme point area);

[0115] Figure 3 is the determinant graph of the Hessian matrix of the objective function;

[0116] Figure 4It is the determinant diagram of the Hessian matrix under the regularized Newton method (including the enlarged diagram of the extreme point area);

[0117] Figure 5 This is a schematic diagram of point selection in three dimensions using the Nelder-Mead algorithm used in the present invention;

[0118] Figure 6 The algorithm flow chart of the Nelder-Mead method adopted in the present invention;

[0119] Figure 7 The overall process of solving the local maximum entropy shape function of the present invention;

[0120] Figure 8 Schematic diagram of neighboring nodes of a single material point in the particle search step in an embodiment of the present invention. DETAILED DESCRIPTION

[0121] In order to further understand the content of the present invention, the present invention is described in detail in conjunction with examples. The following examples are used to illustrate the present invention but are not used to limit the scope of the present invention.

[0122] Embodiment 1: The present invention relates to a method for simulating a large impact deformation process based on a meshless interpolation technique, which comprises the following steps:

[0123] S1. Setting the impact large deformation problem domain , set the global parameters of the algorithm and the model parameters including the material model, which also include the crack propagation algorithm and the contact algorithm;

[0124] S2. Discretize the problem domain into a set of material points {x p}( p represents the material point index number) and a set of node sets {x a}( a Represents the node index number), the physical information of the material (including deformation, stress, internal parameters, etc.) is stored on the material point, the dynamic information (including displacement, velocity, acceleration and temperature, etc.) is stored on the node, and the neighborhood node set is initialized for each material point N H (x p,0 ), initialize physical variables for nodes (node ​​concentrated mass m a,0 , linear momentum l a,0 and the nodal force f a,0 ), where 0 represents the initial time t 0;

[0125] S3. Discrete time domain, with time increments [ t n , t n+1] as a unit cycle to perform S4-S8;

[0126] S4. For each material point x p In the neighborhood node set N H (x p,n ) to construct the objective function F ( l ), and the regularized Newton method is used to solve the variables in the problem domain. l Perform iterative solution and save the final variables l * , and calculate the local maximum entropy interpolation function and its derivatives; the specific steps are as shown in the attached Figure 7 As shown, the following steps are included:

[0127] S4.1. Determine the initial values ​​of the node set and the physical variables to be determined l 0, the node set is obtained in the following way: given a search radius, determine the nodes of each material point within the search radius as the neighborhood nodes of the material point, and input the neighborhood nodes as a node set, such as Figure 8 As shown; the initial value of the physical variable to be determined l 0 is obtained by:

[0128] S4.1.1. Determine whether the problem to be solved in the problem domain is a three-dimensional problem or a two-dimensional problem. If it is a three-dimensional problem, select the four nodes closest to the material point x in the node set to form a tetrahedron, and calculate the centroid coordinates of the material point x with respect to the four vertices of the tetrahedron. If it is a two-dimensional problem, select the three points closest to the material point x in the node set to form a triangle, and calculate the centroid coordinates of the material point x with respect to the three vertices of the triangle. Then, in the node set {x a |x a ∈R d , a =1,…, n} Select the d+1 points closest to the calculated point x to form R d Simple shapes in space, from which we can find the centroid coordinates of the calculation point x about these vertices, which are also the coordinates of each vertex x i Weight w i ;

[0129] Let the shape functions of these points p a = w i ( a For vertex x i The number in the node set), at this time, the shape functions of the other points are all 0, and the partition function value on each node Z a The calculation formula is:

[0130] (1),

[0131] Shape Function N a The calculation formula is:

[0132] (2),

[0133] Expanding formula (2) yields:

[0134] (19),

[0135] Divide the last d formulas by the first formula ( N 1 is not 0, which can be ensured by calculating that point x does not fall on the edge corresponding to vertex x1) and sorting it out, we can get the system of equations:

[0136] (20),

[0137] This linear system of equations contains d equations and d unknowns, and has a unique solution. l 0.

[0138] Table 1 is attached Figure 1 In the specific embodiment described, different parameters β Next, the initial value l Comparison table of guessed and exact values ​​of 0,

[0139] Table 1 Comparison of initial guess values ​​and exact values ​​under different parameters

[0140]

[0141] From Table 1, we can see that the initial value l Although there is a difference between the guessed value and the exact value of 0, the difference is not large, indicating that the initial value l 0 is very close to the extreme point;

[0142] S4.1.2. Let the shape function of the material point be equal to the components of the corresponding barycentric coordinates, and substitute them into the calculation formula of the shape function to obtain a set of variables l The linear equations are solved to obtain the initial values ​​of the variables. l 0, determine whether the value is less than the minimum threshold or greater than the maximum threshold. If so, use the Nelder-Mead algorithm to correct the initial value obtained in S4.1.2, refer to Figure 5 and Figure 6 As shown, the correction method in three-dimensional case is as follows:

[0143] S4.1.2.1. lTake a simple shape in the space (such as {(0,0,0), ( l 0x ,0,0), (0, l 0y ,0), (0,0, l 0z )}), set the tolerance e The maximum number of comparisons NMax and the initial number of operations Nfunk are set to 0;

[0144] S4.1.2.2. Sort the simple shape vertices according to the objective function value, mark the worst vertex as W (Worst), the second worst vertex as S (Second to Worst), the best vertex as B (Best), and the other vertex as N, and calculate F (W)- F (B) If F (W)- F (B)> e Or Nfunk>NMax, then output the coordinates of the center of gravity of the corrected simple shape, which is: (W+S+N+B) / 4, otherwise go to S4.1.2.3 to start correction;

[0145] S4.1.2.3. Discard the worst vertex W, find the centroid M of the triangle formed by the remaining vertices, find the symmetric point R of the worst vertex W about the centroid M, and increase the number of operations Nfunk by 2;

[0146] S4.1.2.4. Comparison of objective function values ​​of variables at symmetric points F (R) and the objective function value of the variable at the worst vertex W F (W), if F (R)< F (W), and proceed to S4.1.2.5, otherwise, proceed to S4.1.2.7;

[0147] S4.1.2.5. Comparison of objective function values ​​of variables at symmetric points F (R) and the objective function value of the variable at the best vertex B F (B) is the size of F (R)≤ F (B), then determine point E on the extension line of the center of gravity M and the symmetry point R. The calculation method of point E is E=M+2(RM). Further compare the objective function value of the variable at point E. F (E) Objective function with variables at symmetric points F (R) size, if F (E)< F (R), replace the worst vertex W with point E, return to S4.1.2.2, if F (E)≥F (R), replace the worst vertex W with the symmetric point R, and return to S4.1.2.2; if F (R)> F (B), then proceed to S4.1.2.6;

[0148] S4.1.2.6. Objective function for comparing variables at symmetric points F (B) The objective function of the variable at the secondary difference vertex S F (S), if F (R)≤ F (S), replace the worst vertex W with the symmetric point R, reduce the number of operations Nfunk by 1, and return to S4.1.2.2. F (R)> F (S), proceed to S4.1.2.7;

[0149] S4.1.2.7. Replace the worst vertex W with the symmetric point R and take the midpoint C of the line connecting the centroid M and the symmetric point R;

[0150] S4.1.2.8. Further comparison of the objective function values ​​of the variables at point C F (C) The objective function value of the variable at the worst vertex W F (W), if F (C)> F (W), replace the worst vertex W with point C, return to S4.1.2.2, if F (C)≤ F (W), shrink the tetrahedron to vertex B, output the shrunken vertex, increase the number of operations Nfunk by 3, and return to S4.1.2.2.

[0151] It should be noted that the present invention does not require the Nelder-Mead algorithm to find the global optimal point, which will take a lot of time. It only needs to pull the initial value back to the neighborhood of the extreme point, so the maximum number of comparisons will not be set to a particularly large number, and it is generally effective when set to 50; in addition, the Nelder-Mead algorithm is rarely called when actually solving the LME problem, but its existence ensures that the calculation will not crash when the iterative value falls into an extreme case, that is, it ensures the stability of the entire algorithm; even, the Nelder-Mead algorithm is not only used to determine whether the initial value is appropriate, but can also be applied to any iteration point.

[0152] After completing the calculation and correction of the above initial values, set the accuracy threshold d , let the physical variable to be determined be λ=λ 0;

[0153] S4.2. Calculate the partition function value at each node Z a , the calculation formula is:

[0154] (1),

[0155] In the formula, β is the Pareto optimization parameter, x represents the material point to be calculated, and x a Represents a node set {x a}, ||·|| represents the norm operation, and the 2-norm is generally used;

[0156] The Pareto optimization parameters β Calculated by the following formula:

[0157] (11),

[0158] In the formula, h is the local characteristic length of the convex set formed by the node set, taking the minimum node spacing, c It is a dimensionless parameter and is generally set to 1.8;

[0159] The partition function value at each node Z a Then, the partition function values ​​of each node are summed up to get Z ,

[0160] ,

[0161] S4.3. Calculation of shape functions N a , the calculation formula is:

[0162] (2);

[0163] S4.4. Objective function F ( l ) is regularized and the corrected objective function is , whose expression is:

[0164] (3),

[0165] in, l is the regularization correction point, that is, the calculation point at the current moment, v is the independent variable of the regularization correction function, v ∈R d , the function satisfies ,▽ F ( l ) represents the gradient vector of the objective function;

[0166] It can be seen from the modified objective function that In Rd has strong convexity, so There must be an extreme value with the following three properties:

[0167] 1. ;

[0168] 2. ;

[0169] 3. ;

[0170] Among them, J and H are the Hessian matrices of the original objective function and the correction function respectively. F ( l )| l ∈R d Using the traditional Newton method, we can get v = l Newton's descent direction at 0 - H( l 0) -1 ▽ F ( l 0), take this direction as the original problem arg min{ F ( l )| l ∈R d}exist l = l The regularized Newton direction rg at 0 ( l 0), it can be proved that the regularized direction rg( l ) has a smaller condition number, and in the area near the extreme point, the algorithm still maintains the second-order convergence rate.

[0171] Depend on The gradient vector r and Hessian matrix J of the original objective function can be obtained:

[0172] ,

[0173] ,

[0174] Where x is the coordinate of the material point whose interpolation function is to be calculated, x a are the node coordinates (a=1,…,N);

[0175] Furthermore, the regularized Newton direction rg ( l ) is expressed as:

[0176] .

[0177] S4.5. Calculate the gradient vector of the objective function ▽ F ( l ), the calculation formula is:

[0178] (4),

[0179] In the formula, N represents the number of nodes, a Indicates the node index number;

[0180] S4.6. Perform norm operation on the gradient vector to obtain accuracy and compare it with the set accuracy threshold d Compare and if the calculated progress is not greater than the accuracy threshold, output the variable at this time as the final variable l * , otherwise proceed to S4.7- S4.8;

[0181] S4.7. Calculate the Hessian matrix H ( l ), the regularized Newton descent direction is calculated based on the Hessian matrix and the gradient vector,

[0182] The calculation formula of the Hessian matrix is:

[0183] (5),

[0184] In the formula, H( l ) represents the Hessian matrix, r represents the gradient vector, that is, r(x, l )=▽ F ( l ), I represents the identity matrix;

[0185] The calculation formula for the regularized Newton descent direction is:

[0186] (6),

[0187] In the formula, rg ( l ) represents the direction of regularized Newton descent;

[0188] Compare with Figure 3 and Figure 4 It can be seen that the regularized Newton method repairs the area where the determinant of the Hessian matrix of the original objective function is zero. After the repair, the regularized Newton direction can be constructed at any point in the space for solution;

[0189] S4.8. Iteration of the physical variable to be determined , and order , return to step S4.2,

[0190] The iteration formula is:

[0191] (7),

[0192] In the formula, t ( l ) represents the optimal step size in the direction of regularized Newton descent;

[0193] The optimal step size in the direction of regularized Newton descent during each iteration t ( l ) is initially set to 1. After the physical variable to be determined is iterated, S4.8 calculates the objective function values ​​before and after the iteration. If the objective function value after the iteration is less than the objective function value before the iteration, that is, , then directly return to step S4.2. Otherwise, if the objective function value after iteration is not less than the objective function value before iteration, the optimal step length in the regularized Newton descent direction is calculated by the exact line search algorithm. t ( l ), and iterate again according to the iteration formula, and after re-iteration, return to step S4.2:

[0194] After the above iterations, the final variable is obtained l * ;

[0195] S4.9. The final variable l * Substituting into the shape function formula of S4.3, the local maximum entropy interpolation function is calculated, that is,

[0196] ,

[0197] in, , is the partition function in the sense of statistical mechanics, ||·|| is the norm operator, and the 2-norm is generally used;

[0198] The derivative of the local maximum entropy interpolation function is calculated according to the formula for calculating the derivative of the local maximum entropy interpolation function:

[0199] (8),

[0200] in, (9),

[0201] (10),

[0202] In the formula, r * Represents the gradient vector of the original objective function, J * Represents the Hessian matrix of the original objective function.

[0203] S5. Update the node concentrated mass based on the local maximum entropy interpolation function, update the node coordinates at the next moment, calculate the local deformation gradient increment based on the node coordinates of the next node and the derivative of the maximum entropy interpolation function, and calculate the local deformation gradient at the next moment based on the local deformation gradient increment. The specific steps include:

[0204] S5.1. Update the node concentrated mass, the formula is:

[0205] (12),

[0206] In the formula, N a represents the local maximum entropy interpolation function, m represents the concentrated mass at the node, a Indicates the node index number, n express n time, p Indicates the material point index number, N H (x a,n ) represents node x a The neighborhood material point set of ;

[0207] S5.2. Update the node coordinates of the next node, expressed as:

[0208] (13),

[0209] In the formula, l a,n Indicates that the index number is a The node is n The linear momentum at time, f a,n Indicates that the index number is a The node is n Nodal force at the moment;

[0210] S5.3. Calculate the local deformation gradient increment based on the node coordinates of the next node and the derivative of the maximum entropy interpolation function. The calculation formula is:

[0211] (14),

[0212] In the formula, F p,n→n+1 Indicates from n The local deformation gradient increment from moment to the next moment, N H (x p,n ) represents each material point x p The set of neighboring nodes;

[0213] S5.4. Calculate the local deformation gradient of the material point at the next moment. The calculation formula is:

[0214] (15),

[0215] In the formula, F p,n express n The local deformation gradient at time, F p,n+1Represents the local deformation gradient at the next moment.

[0216] S6. The local deformation gradient increment F obtained in S5 is p,n→n+1 and the local deformation gradient F p,n+1 Substitute into the material model (specific material model specific analysis), and obtain the stress P of the material point n+1 ;

[0217] S7. Combine the P obtained in S6 n+1 and boundary conditions, calculate the node internal stress f caused by the incremental motion at time n+1 int a,n+1 and the node external force f ext a,n+1 , update the node force f a,n+1 and linear momentum l a,n+1 ,Right now:

[0218] (16),

[0219] In the formula, f int a,n+1 Indicates that the index number is a The node internal force at the next moment, f ext a,n+1 Indicates that the index number is a The node external force at the next moment, f a,n+1 Indicates that the index number is a The node force of the node at the next moment;

[0220] (17),

[0221] In the formula, l a,n Indicates that the index number is a The node is n The linear momentum at time, l a,n+1 Indicates that the index number is a The linear momentum of the node at the next moment;

[0222] After the update, determine whether the calculation has reached the last time step. If not, go to S8. If so, go to S9.

[0223] S8. Update material point coordinates, neighborhood node sets and other physical information:

[0224] The formula for updating material point coordinates is:

[0225] (18),

[0226] In the formula, x a,n+1 Indicates the node coordinates at the next moment, x p,n+1represents the coordinates of the material point at the next moment, N H (x p,n ) represents the material point x p The set of neighboring nodes;

[0227] Update the next set of neighbor nodes N H (x p,n+1 ), which is obtained by: given a search radius, determining the nodes of each material point within the search radius as the neighborhood nodes of the material point, and inputting the neighborhood nodes as the neighborhood node set; traversing the neighborhood node set of the material point N H (x p,n+1 ), construct the neighborhood material point set of the node N H (x a,n+1 ).

[0228] 9. Complete the impact process simulation and output the results.

[0229] Embodiment 2: This embodiment takes the bullet hitting the target plate as an example to explain the content of the present invention in detail.

[0230] S1. Model setting: Set the impact large deformation problem domain: The bullet hits the target plate vertically with an initial speed of 970m / s. Use the parametric modeling method to complete the modeling of the bullet (including the shell and the inner core) and the target plate (including the impact area and the external area). Use Femap software for meshing. Because the axial impact simulation model of the bullet meets the symmetry condition, in order to save calculation time, take one-eighth of the complete model for calculation. Set global parameters and model parameters including material model. Model parameters include material model, crack extension algorithm and contact algorithm. The specific steps are as follows:

[0231] Material model setting: The bullet case and core select the parameters of copper and lead in the J2 viscoelastic-plastic model respectively, and use the Mie-Gruneisen state equation; the layer plate selects the high-strength steel parameters in the J2 viscoelastic-plastic model, and also uses the Mie-Gruneisen state equation; the J2 viscoelastic-plastic model has effects such as strain hardening, high-temperature softening and strain rate correlation. The Mie-Gruneisen state equation links the pressure and volume of a solid at a given temperature, and is used to determine the pressure in an impact-compressed solid. The combination of the two can effectively simulate the material behavior under high-pressure conditions.

[0232] The material failure model uses the EigenErosion crack propagation algorithm, which realizes the capture of crack tip propagation based on energy and variational principles, and can effectively avoid the grid or discrete correlation of the crack propagation path.

[0233] The contact algorithm between the bullet and the target plate adopts the Seizing multi-body dynamic contact algorithm, which automatically captures the contact by dynamically updating the neighborhood, and can well simulate the erosion, penetration and other problems in the process of projectile penetration.

[0234] S2. Discretize the problem domain into a set of material points {x p} and a set of nodes {x a}, in which the physical information of the material (including deformation, stress, internal parameters, etc.) is stored at the material points, and the dynamic information (including displacement, velocity, acceleration and temperature, etc.) is stored at the nodes. The neighborhood node set is initialized for each material point, and the physical variables are initialized for the nodes. The physical variables include node concentrated mass, linear momentum and node force.

[0235] S3. Time domain discretization: The simulation time is discretized into different time intervals according to the global parameters { t 0, t 1, t 2, … , t n , … , T}, for each time increment [ t n , t n+1 ] loop S4-S8;

[0236] S4. For each material point, construct the objective function on the neighborhood node set F ( l ), and the regularized Newton method is used to solve the variables in the problem domain. l Perform iterative solution and save the final variables l * , calculate the local maximum entropy interpolation function and its derivative, specifically including the following steps:

[0237] S4.1. Determine the initial values ​​of the node set and the physical variables to be determined l 0, the node set is obtained by the following method: The specific method of neighboring particle search is: given the search radius, determine the neighborhood nodes of each material point, as shown in the attached Figure 8 As shown, the big dots represent nodes, the small dots represent material points, and the dotted lines represent the search range (here the material point is the center of the circle and the second layer particle distance is the radius search). The nodes that fall within the search range (including the boundary) are the neighborhood nodes, and the initial values ​​of the physical variables to be calculated are obtained. l 0:

[0238] Initialize calculation parameters: Input the N neighboring nodes of material point x as a node set X , calculated by the node set X The local characteristic length of the convex seth (generally the minimum node spacing), the form of the local maximum entropy interpolation function is determined by the dimensionless parameter c This parameter is generally set to 1.8, according to the formula Calculate the Pareto optimal parameters β , used to solve the objective function of the local maximum entropy problem, each time the node set is updated (that is, the local feature length is updated h After) need to be updated β value;

[0239] The initial value of the physical variable to be determined is obtained l The specific steps of 0 are:

[0240] S4.1.1. Determine that the problem domain is a three-dimensional problem. X Select the 4 points closest to the material point x to form a tetrahedron, and then find the centroid coordinates of the material point x about these four vertices ( w 1, w 2, w 3, w 4);

[0241] S4.1.2. Let the shape functions of these points p i Equal to the corresponding barycentric coordinate components w i , substitute into the expression of shape function N a =Z a / Z and sort them out to get a set of physical variables to be determined l The linear equations can be solved to get the initial value l 0.

[0242] Nelder-Mead correction: First determine the situation at the iteration point. If the function value Z ( l 0) Very large (for example, greater than 10 32 ) or very small (e.g. less than 10 -32 ), the position of the iteration point needs to be corrected by the Nelder-Mead algorithm. The correction algorithm is shown in the attached Figure 5 and Figure 6 For convenience, the four initial points are usually taken as {(0,0,0), ( l x ,0,0), (0, l y ,0), (0,0, l z )},in l x , ly and l z Iteration points l 0, after a finite number of iterative calculations, the output is l 0 is used as the initial point of the regularized Newton method. The iterative process is the same as steps S4.1.2.1-S4.1.2.8 of Example 1. This example will not be described in detail. This point is not the global optimal point, but is generally located in the neighborhood of the extreme point.

[0243] S4.2. Node set determined according to the initialization phase X , Pareto optimal parameters β and initial value l 0 Calculate the partition function value at each node , sum to get Z (λ0);

[0244] S4.3. Calculation of shape functions N a =Z a / Z;

[0245] S4.4. Objective function F ( l ) for regularization correction;

[0246] S4.5. Calculate the gradient vector ▽ of the objective function according to formula (4) F ( l );

[0247] S4.6. Perform norm operation on the gradient vector to obtain accuracy and compare it with the set accuracy threshold. If the calculated progress is not greater than the accuracy threshold, output the variable at this time as the final variable l * , otherwise proceed to S4.7- S4.8;

[0248] S4.7. According to formula (5), the modified Hessian matrix H is obtained ( l ), according to formula (6), the regularized Newton descent direction rg ( l );

[0249] S4.8. Use formula (7) to get the temporary iteration point. In this process, if , then the search step length in the regularized Newton descent direction is t ( l ) takes 1, otherwise, if , then an exact line search algorithm is used to find the search step size in the direction of regularized Newton descent t ( l ),Right now , numerically, it can be implemented using the binary search method;

[0250] S4.9. The final variable l * Substituting into the shape function formula of S4.3, the local maximum entropy interpolation function is calculated, that is,

[0251] ,

[0252] in, , is the partition function in the sense of statistical mechanics, ||·|| is the norm operator, and the 2-norm is generally used;

[0253] The derivative of the local maximum entropy interpolation function is calculated according to the formula for calculating the derivative of the local maximum entropy interpolation function:

[0254] (8),

[0255] in, (9),

[0256] (10),

[0257] In the formula, r * Represents the gradient vector of the original objective function, J * Represents the Hessian matrix of the original objective function.

[0258] S5. Update the node concentrated mass based on the local maximum entropy interpolation function, update the node coordinates at the next moment, calculate the local deformation gradient increment based on the node coordinates of the next node and the derivative of the maximum entropy interpolation function, and calculate the local deformation gradient at the next moment based on the local deformation gradient increment. The specific steps include:

[0259] S5.1. Update the node concentrated mass, the formula is:

[0260] (12),

[0261] In the formula, m represents the concentrated mass at the node, a Indicates the node index number, n express n time, p Indicates the material point index number, N H (x a,n ) represents the node x a The neighborhood material point set of ;

[0262] S5.2. Update the node coordinates of the next node, expressed as:

[0263] (13),

[0264] In the formula, l a,n Indicates that the index number is a The node is n The linear momentum at time, f a,n Indicates that the index number is a The node is n Nodal force at the moment;

[0265] S5.3. Calculate the local deformation gradient increment based on the node coordinates of the next node and the derivative of the maximum entropy interpolation function. The calculation formula is:

[0266] (14),

[0267] In the formula, F p,n→n+1 Indicates from n The local deformation gradient increment from moment to the next moment, N H (x p,n ) represents each material point x p The set of neighboring nodes;

[0268] S5.4. Calculate the local deformation gradient of the material point at the next moment. The calculation formula is:

[0269] (15),

[0270] In the formula, F p,n express n The local deformation gradient at time, F p,n+1 Represents the local deformation gradient at the next moment;

[0271] S6. Perform Hencky processing on the local deformation gradient obtained in S5 , and perform elastic-plastic decomposition The offspring is introduced into the J2 viscoelastic-plastic model and the stress of the material point is solved iteratively. The equivalent yield stress of the J2 viscoelastic-plastic model with power-law hardening and strain rate dependence can be expressed as: ,in is the plastic strain, is the plastic strain rate, is the initial quasi-static yield stress at the reference temperature, and are the reference plastic strain and reference plastic strain rate, n is the hardening index, and m is the rate sensitivity index, and the values ​​are related to the specific material;

[0272] S7. Combining the material point stress and boundary conditions obtained in S6, calculate the node internal stress f caused by the incremental motion at time n+1 int a,n+1 and the node external force fext a,n+1 , update the nodal force and linear momentum according to formula (16) and formula (17);

[0273] S8. Update material point coordinates , Neighborhood Node Set and other physical information, updating the neighborhood material point set of the node;

[0274] S9. Calculate until the last time step to complete the simulation of the bullet hitting the target plate and output the results.

[0275] The present invention is described in detail above in conjunction with the embodiments, but the contents described are only preferred embodiments of the present invention and cannot be considered to limit the scope of implementation of the present invention. All equivalent changes and improvements made according to the scope of application of the present invention should still fall within the scope of the patent coverage of the present invention.

Claims

1. A method for simulating large deformation during impact based on meshless interpolation technology, characterized in that: The following steps are involved: S1. Set the impact large deformation problem domain, set the global parameters and model parameters including material model; S2. Discretize the problem domain into a set of material points and a set of nodes, initialize the neighborhood node set for each material point, and initialize the physical variables for the nodes. The physical variables include the mass, linear momentum and node force of the node set. S3. Discrete time domain, with time increments [ t n , t n+1 ] as a unit cycle to perform S4-S8; S4. For each material point, construct the objective function on the neighborhood node set F ( λ ), and the regularized Newton method is used to solve the variables in the problem domain. λ Perform iterative solution and save the final variables λ * , calculate the local maximum entropy interpolation function and its derivative; The variables in the problem domain are regularized using the Newton method. λ The specific steps for iterative solution are: S4.

1. Based on the node set and initial values ​​of the physical variables to be determined in the initialization phase λ 0, set the accuracy threshold δ , let the physical variable to be determined be λ=λ 0; the initial value of the physical variable to be determined λ 0 is obtained by: S4.1.

1. Determine whether the problem to be solved in the problem domain is a three-dimensional problem or a two-dimensional problem. If it is a three-dimensional problem, select the four nodes closest to the material point x in the node set to form a tetrahedron, and calculate the centroid coordinates of the material point x with respect to the four vertices of the tetrahedron. If it is a two-dimensional problem, select the three points closest to the material point x in the node set to form a triangle, and calculate the centroid coordinates of the material point x with respect to the three vertices of the triangle. S4.1.

2. Let the shape function of the material point be equal to the components of the corresponding barycentric coordinates, and substitute them into the calculation formula of the shape function to obtain a set of variables λ The linear equations are solved to obtain the initial values ​​of the variables. λ 0； S4.

2. Calculate the partition function value at each node Z a , the calculation formula is: （1）， In the formula, β is the Pareto optimization parameter, x represents the material point to be calculated, and x a Represents a node set {x a }, ||·|| represents norm operation; The partition function values ​​of each node are summed up to get ; S4.

3. Calculation of shape functions N a , the calculation formula is: （2） ； S4.

4. Objective function F ( λ ) is regularized and the corrected objective function is , whose expression is: （3）， in, λ is the regularization correction point, that is, the calculation point at the current moment, v is the independent variable of the regularized correction function, ▽ F ( λ ) represents the gradient vector of the objective function; S4.

5. Calculate the gradient vector of the objective function ▽ F ( λ ), the calculation formula is: （4）， In the formula, r(x, λ ) represents the gradient vector, N represents the number of nodes, a Indicates the node index number; S4.

6. Perform norm operation on the gradient vector to obtain accuracy and compare it with the set accuracy threshold δ Compare and if the calculated progress is not greater than the accuracy threshold, output the variable at this time as the final variable λ * , proceed to S4.9, otherwise proceed to S4.7-S4.8; S4.

7. Calculate the Hessian matrix H of the objective function ( λ ), the regularized Newton descent direction is calculated based on the Hessian matrix and the gradient vector, The calculation formula of the Hessian matrix is: （5）， In the formula, r(x, λ ) represents the gradient vector, I represents the identity matrix; The calculation formula for the regularized Newton descent direction is: （6）， In the formula, rg( λ ) represents the direction of regularized Newton descent; S4.

8. Iteration of the physical variable to be determined , and order , return to step S4.2, The iteration formula is: （7）， In the formula, t ( λ ) represents the search step size in the direction of regularized Newton descent; S4.

9. The final variable λ * Substitute it into the shape function formula of S4.3 to calculate the local maximum entropy interpolation function. The derivative of the local maximum entropy interpolation function is calculated according to the local maximum entropy interpolation function. The calculation formula of the derivative of the local maximum entropy interpolation function is: （8）， in, (9), （10）， In the formula, r * Represents the gradient vector of the original objective function, J * Represents the Hessian matrix of the original objective function; S5. Update the node concentrated mass based on the local maximum entropy interpolation function, update the node coordinates at the next moment, calculate the local deformation gradient increment based on the node coordinates of the next node and the derivative of the maximum entropy interpolation function, and calculate the local deformation gradient at the next moment based on the local deformation gradient increment; S6. Substitute the local deformation gradient increment and the local deformation gradient obtained in S5 into the material model to obtain the stress of the material point; S7. Integrate the stress and boundary conditions of the material points obtained in S6, update the node force and linear momentum, and determine whether the calculation is to the last time step. If not, go to S8. If yes, go to S9. S8. Update the material point coordinates and the neighborhood node set, update the neighborhood material point set of the node and return to S4; S9. Complete the impact process simulation and output the results.

2. The method for simulating large impact deformation process based on meshless interpolation technology according to claim 1 is characterized in that: The iterative process in S4.8 regularizes the search step size in the Newton descent direction t ( λ ) is initially set to 1. After iterating the physical variable to be determined, S4.8 calculates the objective function values ​​before and after the iteration. If the objective function value after the iteration is less than the objective function value before the iteration, it directly returns to step S4.

2. If the objective function value after the iteration is not less than the objective function value before the iteration, the optimal step length in the regularized Newton descent direction is calculated by the precise line search algorithm. t ( λ ), and iterate again according to the iteration formula, and return to step S4.2 after the re-iteration.

3. The method for simulating large impact deformation process based on meshless interpolation technology according to claim 1 is characterized in that: S4.1.2 solves the initial values ​​of the variables λ 0, the partition function values ​​of each node are summed up to get Z ( λ 0), determine whether the value is less than the minimum threshold or greater than the maximum threshold. If so, use the Nelder-Mead algorithm to correct the centroid coordinates obtained in S4.1.

1. The correction method is as follows: S4.1.2.

1. λ Select one in the space with λ 0 associated simple shape, set tolerance ε and the maximum number of comparisons NMax, set the initial number of operations Nfunk to 0; S4.1.2.

2. Sort the simple shape vertices according to the objective function value, mark the worst vertex as W, the second worst vertex as S, and the best vertex as B, and calculate the objective function value of the variable at the worst vertex W F (W) and the objective function value of the variable at the best vertex B F (B) if F (W)- F (B)> ε Or if Nfunk>NMax, then output the corrected center of gravity coordinates, otherwise go to S4.1.2.3 to start correction; S4.1.2.

3. Discard the worst vertex W, find the center M of the face or line formed by the remaining vertices, find the symmetric point R of the worst vertex W about the center M, and increase the number of operations Nfunk by 2; S4.1.2.

4. Comparison of objective function values ​​of variables at symmetric points F (R) and the objective function value of the variable at the worst vertex W F (W), if F (R)< F (W), then proceed to S4.1.2.5, otherwise proceed to S4.1.2.7; S4.1.2.

5. Comparison of objective function values ​​of variables at symmetric points F (R) and the objective function value of the variable at the best vertex B F (B) size, if F (R)≤ F (B), then determine point E on the extension line of the center of gravity M and the symmetry point R. The calculation method of point E is E=M+2(RM). Further compare the objective function value of the variable at point E. F (E) Objective function with variables at symmetric points F (R) size, if F (E)< F (R), replace the worst vertex W with point E, return to S4.1.2.2, if F (E)≥ F (R), replace the worst vertex W with the symmetric point R, and return to S4.1.2.2; if F (R)> F (B), then proceed to S4.1.2.6; S4.1.2.

6. Objective function for comparing variables at symmetric points F (B) The objective function of the variable at the secondary difference vertex S F (S), if F (R)≤ F (S), replace the worst vertex W with the symmetric point R, reduce the number of operations Nfunk by 1, and return to S4.1.2.

2. F (R)> F (S), go to S4.1.2.7; S4.1.2.

7. Replace the worst vertex W with the symmetric point R and take the midpoint C of the line connecting the centroid M and the symmetric point R; S4.1.2.

8. Further comparison of the objective function values ​​of the variables at point C F (C) The objective function value of the variable at the worst vertex W F (W), if F (C)> F (W), replace the worst vertex W with point C, return to S4.1.2.2, if F (C)≤ F (W), shrink the simple shape to the best vertex B, output the shrunk vertex, increase the number of operations Nfunk by 3, and return to S4.1.2.

2.

4. The method for simulating large impact deformation process based on meshless interpolation technology according to claim 1 is characterized in that: The Pareto optimization parameters in S4.2 are calculated by the following formula: （11）， In the formula, h is the local characteristic length of the convex set formed by the node set, taking the minimum node spacing, γ is a dimensionless parameter; After the neighborhood node set in S8 is updated, the Pareto optimization parameters are recalculated.

5. The method for simulating large impact deformation process based on meshless interpolation technology according to claim 1 is characterized in that: The specific steps of S5 include: S5.

1. Update the node concentrated mass, the formula is: （12）， In the formula, N a represents the local maximum entropy interpolation function, m represents the concentrated mass at the node, a Indicates the node index number, n express n time, p Indicates the material point index number, N H (x a,n ) represents node x a The neighborhood material point set of ; S5.

2. Update the node coordinates of the next node, expressed as: （13）， In the formula, l a,n Indicates that the index number is a The node is n The linear momentum at time, f a,n Indicates that the index number is a The node is n Nodal force at the moment; S5.

3. Calculate the local deformation gradient increment based on the node coordinates of the next node and the derivative of the maximum entropy interpolation function. The calculation formula is: （14）， In the formula, F p,n→n+1 Indicates from n The local deformation gradient increment from moment to the next moment, N H (x p,n ) represents each material point x p The set of neighboring nodes; S5.

4. Calculate the local deformation gradient of the material point at the next moment. The calculation formula is: （15）， In the formula, F p,n express n The local deformation gradient at time, F p,n+1 Represents the local deformation gradient at the next moment.

6. The method for simulating large impact deformation process based on meshless interpolation technology according to claim 5 is characterized in that: The updated node force in S7 includes node internal force and node external force, and the calculation formula of the node force is: （16）， In the formula, f int a,n+1 Indicates that the index number is a The node internal force at the next moment, f ext a,n+1 Indicates that the index number is a The node external force at the next moment, f a,n+1 Indicates that the index number is a The node force of the node at the next moment; The calculation formula for updating the linear momentum in S7 is: （17）， In the formula, l a,n Indicates that the index number is a The node is n The linear momentum at time, l a,n+1 Indicates that the index number is a The linear momentum of the node at the next moment.

7. The method for simulating large impact deformation process based on meshless interpolation technology according to claim 5 is characterized in that: The formula for updating the material point coordinates in S8 is: （18）， In the formula, x a,n+1 Indicates the node coordinates at the next moment, x p,n+1 represents the coordinates of the material point at the next moment, N H (x p,n ) represents the material point x p The set of neighboring nodes; Update the next set of neighbor nodes N H (x p,n+1 ), which is obtained by: given a search radius, determining the nodes of each material point within the search radius as the neighborhood nodes of the material point, and inputting the neighborhood nodes as the neighborhood node set; traversing the neighborhood node set of the material point N H (x p,n+1 ), construct the neighborhood material point set of the node N H (x a,n+1 ).

8. An application of the method for simulating a large impact deformation process based on a meshless interpolation technique according to claim 1, wherein the method is applied to simulating a large impact deformation process when a bullet hits a target plate.

Images