A parallel optimal transport meshless impact collision process simulation method and its application

Abstract

The present invention discloses a parallel optimal transport gridless impact collision process simulation method, which belongs to the field of computational mechanics and numerical simulation technology, and includes the following steps: establishing a geometric model; discretizing the geometric model in space, discretizing the time domain into a number of time steps; establishing a topological connection relationship diagram of material points and decomposing them, confirming boundary material nodes and shared nodes; material point subsets are respectively assigned to each parallel processor, and a shared data exchange table is established; each processor calculates internal forces and contact forces, and obtains local mass matrix, local node force and local node acceleration; the shared data exchange table exchanges and assembles data, and obtains global calculation results; enters the next time step, performs iterative calculation or outputs calculation results. The present invention combines graph topology analysis and distributed multi-process parallelization, and can effectively use supercomputers to accelerate OTM simulations, effectively saving calculation time.

Description

Technical Field

[0001] The invention belongs to the technical field of computational mechanics and numerical simulation, and in particular relates to a parallel optimal transport gridless impact collision process simulation method. Background Art

[0002] Impact collision is a hot topic in the fields of aerospace, involving complex features such as extreme deformation, dynamic crack initiation and propagation, geometric distortion, and material fission. In particular, high-speed and ultra-high-speed impact collision problems with speeds above 1 km / s may even involve complex processes such as material phase change, which are difficult to simulate using traditional numerical methods based on meshes. Traditional numerical methods based on Lagrangian meshes and Euler meshes, such as finite element method, finite volume method, and finite difference method, are limited in their application in practical engineering problems involving high-speed impact collisions due to their inherent mesh characteristics. When solving high-speed and ultra-high-speed impact collision problems, the Lagrangian mesh-based solution method is prone to grid distortion and grid entanglement caused by large grid deformation, which requires methods such as grid reconstruction that easily introduce errors and greatly increase the computational load to solve; the Euler mesh-based solution method is difficult to handle the dynamic response of materials related to the loading history and crack propagation problems, and it is difficult to determine the free surface, deformable boundary, and interface of moving materials. Therefore, in order to solve the above difficulties of traditional numerical methods based on meshes, a new generation of numerical calculation methods, meshless methods, have been widely and rapidly developed. The goal of the meshless method is to completely (at least partially) eliminate the dependence on the mesh by discretizing the entire problem domain in a node-based manner. Its main advantages are: ① The calculation process does not involve the mesh, there is no mesh distortion, mesh reconstruction and other difficulties, which is very suitable for dealing with large deformation problems; ② H-type and P-type adaptive implementation becomes relatively simple; ③ It is easy to construct shape functions of arbitrary order continuity; ④ The convergence results are often much better than those based on the mesh; ⑤ The volume self-locking problem can be alleviated by adjusting the shape function expansion parameters.

[0003] The Optimal Transportation Meshfree method (OTM method) effectively utilizes the spatial discretization of material points, the local maximum entropy interpolation function and the time discretization of the optimal transport theory, avoids the instability of the calculated tensile stress, satisfies the "weak" Kronecker property, and can accurately impose displacement boundary conditions. In addition, the OTM method can dynamically adjust the connection relationship between material points and nodes, and can automatically detect dynamic contact surfaces without other complex interface tracking technologies, which has great advantages in solving high-speed / ultra-high-speed collision problems. However, the computational load required by the OTM method to solve impact collision problems is large, which leads to low computational efficiency and slow simulation speed. How to improve the computational efficiency of the OTM method in solving impact collision processes is one of the key issues that the OTM method needs to solve. Summary of the invention

[0004] The purpose of the present invention is to provide a parallel optimal transport gridless impact collision process simulation method to solve the problem of insufficient computational efficiency of the existing OTM method in simulating impact collision problems.

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

[0006] The present invention relates to a parallel optimal transport gridless impact collision process simulation method, which comprises the following steps:

[0007] S1. Setting a problem domain for solving the impact and collision problem, and then establishing a geometric model;

[0008] S2. Discretize the geometric model spatially, discretize the geometric model into a set of node sets and a set of material point sets using a finite element grid, initialize the node set and the material point set, and discretize the time domain into a number of time steps;

[0009] S3. Establish a topological connection relationship diagram of material points, decompose the topological connection relationship diagram into material point subsets equal to the number of processors according to the number of parallel processors, and determine the boundary material nodes and shared nodes in each material point subset;

[0010] S4. Allocate the material point subsets to each parallel processor and establish a shared data exchange table;

[0011] S5. Each processor calculates internal forces based on the assigned material point subsets, calculates the contact forces of the nodes on the two objects in the impact collision problem through the contact algorithm, and obtains the local mass matrix, local node forces and local node accelerations;

[0012] S6. perform data exchange and assembly according to the shared data exchange table, and calculate the global mass matrix, global node force and global acceleration;

[0013] S7. Update the material nodes and shared data exchange table, enter the next time step, and determine whether this time step is the last time step set. If not, return to S5 to enter the next iterative calculation; if so, output the calculation result of S6, complete the impact collision process simulation, and obtain the final kinematic information and material physics information.

[0014] Preferably, in S2, the initial kinematic information is stored in the node, and the initial material physical information is stored in the material point set; the kinematic information includes displacement, velocity, acceleration and temperature, and the material physical information includes deformation, stress and material internal parameters;

[0015] The specific steps of S5 include:

[0016] S5.1. Initialize nodes and material points, and calculate shape functions and shape function derivatives of initial nodes;

[0017] The shape function adopts the local maximum entropy shape function, and the calculation formula of the shape function is:

[0018] ,

[0019] in, N a Indicates that the index number is a The shape function of the node, x Represents the position of any node in the node set. x a Indicates that the index number is a The location of the node, β represents the Pareto optimization parameter, and Z For the introduced node x The relevant variables must satisfy the following formula:

[0020] ,

[0021] ,

[0022] in, λ Represents variables related to neighborhood search radius;

[0023] Shape function derivative ▽N a ( x ) is calculated as follows:

[0024] ,

[0025] in, J The variables introduced with respect to node x follow the following form:

[0026] ,

[0027] in, represents the partial derivative of the variable λ, r represents the gradient vector of the objective function, N Represents the number of nodes in the neighborhood of the node, and the expression of the gradient vector is:

[0028] ;

[0029] S5.2. Calculate the deformation gradient at the material point using the formula:

[0030] ,

[0031] In the formula, F p,k The index number is p The material point is k The deformation gradient of the time step, k For the k time steps, x p,k The index number is p The material point is k The position of the time step;

[0032] S5.3. Update the node position and calculate the updated material point deformation gradient. The update formula of the node position is:

[0033] ,

[0034] In the formula, x a,k+1 For the k +The position of the node at 1 time step, v a,k The index number is a The node in k The speed of the movement in time steps, △t is the time step, a a,k The index number is a The node in k The acceleration of the time step;

[0035] The update formula of the material point deformation gradient is:

[0036] ,

[0037] in, F p,k+1 For the k+1 deformation gradient of the material point at time step, φ k→k+1 For k steps to k +1-step transfer mapping;

[0038] S5.4. Deformation gradient of material point F p,k+1 Decomposed into a pure rotation tensor R and a pure stretch tensor V , according to the pure stretch tensor V Calculate strain;

[0039] Deformation gradient of material point F p,k+1 The formula for decomposing into a pure rotation tensor R and a pure stretch tensor V is:

[0040] ,

[0041] According to the pure tensor V The calculation formula for calculating material strain is:

[0042] ,

[0043] in, represents material strain;

[0044] S5.5. Based on the above material strain, the stress of the material point is calculated. The calculation formula of the yield stress of the material point is:

[0045] ,

[0046] in, σ y is the yield stress of the material point, σ 0 is the initial quasi-static yield stress of the material point, and is the effective plastic strain and strain rate, and For reference effective plastic strain and strain rate, T is the absolute temperature, n For the reinforcement index, m is the strain rate exponent, T 0 is the reference temperature, T m is the melting temperature, l is the thermal softening index; the plastic strain and plastic strain rate are calculated by the material strain Iterate and loop to obtain;

[0047] S5.6. Determine whether the material point is damaged and cracks are generated. If it is not damaged, interpolate the stress of the material point to each node in its neighborhood to obtain the local node force caused by material deformation at each node; the node force caused by material deformation is calculated by the following formula:

[0048] ,

[0049] in, represents the nodal forces due to material deformation, M represents the number of material points in the node neighborhood, w represents the volume of a material point, P p Represents the first Piola-Kirchhoff stress of a material point;

[0050] S5.7. Calculate the tension applied on the boundary and the contact force occurring during collision, and calculate the local node force of each processor.

[0051] Preferably, the S5.1 initialization of nodes and material points includes initializing the mass, velocity and node force of each node; the calculation formula of the acceleration of the node in S5.2 is:

[0052] ,

[0053] In the formula, f a The index number is a The nodal forces at the nodes.

[0054] Preferably, when determining whether a material point is damaged and cracks are generated in S5.6, an EigenFracture model is used to calculate the equivalent energy release rate of the material point. When the equivalent energy release rate of the material point is not lower than the critical energy release rate, it is determined that the material point is damaged and cracks are generated;

[0055] The calculation formula of the material point equivalent energy release rate is:

[0056] ,

[0057] In the formula, G ( x p ) is the index number p The equivalent energy release rate of the material point is ch p is the material point size, V ϵ is the material point volume, W is the elastic energy stored in the material point, x q is the material point position, For the average area, is the intrinsic strain, ▽ u is the displacement gradient, V q is the energy parameter.

[0058] Preferably, the calculation formula of the contact force generated during the collision in S5.7 is:

[0059] ,

[0060] In the formula, f c is the contact force, K p is the contact stiffness, is the contact depth of the contact pair, F 1 and F 2 is the equivalent contact force calculated when two nodes are in contact;

[0061] The calculation formula of the contact depth of the contact pair is:

[0062] ,

[0063] In the formula, R i and R j Represent the contact diameters of the two nodes, r i j Represents the distance between two nodes;

[0064] The equivalent contact forces of the two nodes are calculated using the following formulas:

[0065] ,

[0066] ,

[0067] In the formula, m' is the average mass of the two nodes, c is the speed of sound in the two objects, △t is the time step.

[0068] Preferably, in S5.7, the local node force of each node caused by material deformation, the tension on the boundary and the contact force occurring during the collision are added together to obtain the local node force, that is:

[0069] ,

[0070] in, The first kThe local nodal forces at each time step, t is the tensile force applied on the boundary, f c is the contact force.

[0071] Preferably, the specific steps of performing data exchange and assembly according to the shared data exchange table in S6 and calculating the global mass matrix, global node force and global acceleration are:

[0072] S6.1. Each processor sends the local node force, local mass matrix, and local acceleration of the shared node contained in the data exchange table to other processors, and receives the calculation results of the shared node in other processors;

[0073] S6.2. Each processor synchronizes and assembles its own data and the data received from other processors to obtain the global mass matrix, global node force, global acceleration and global velocity.

[0074] Preferably, the calculation formula of the S6.2 global mass matrix is:

[0075] ,

[0076] In the formula, m a,k For the k The global mass matrix for time steps, For the k The local mass matrix calculated by processor number I at time step, Q For shared nodes, C k To share data exchange tables, No. k The local mass matrix obtained by exchanging nodes is shared by other processors at the time step;

[0077] The calculation formula of the global nodal force is:

[0078] ,

[0079] In the formula, f a,k For the k The global nodal forces for each time step, For the k The node force calculated by processor number I at the time step, For the k The remaining processors share the node forces obtained by node exchange during the time step;

[0080] The calculation formula of the global acceleration is:

[0081] ,

[0082] The calculation formula of the global velocity is:

[0083] .

[0084] Preferably, the specific steps of S4 establishing a shared data exchange table are:

[0085] S4.1. Calculate the neighborhood of each material point on the processor;

[0086] S4.2. Determine the boundary material points on each processor according to the obtained topological structure diagram;

[0087] S4.3. According to the boundary material points of each processor, determine the nodes in the neighborhood of each boundary material point to obtain the shared nodes;

[0088] S4.4. Save the shared nodes to the shared data exchange table to obtain the shared data exchange table.

[0089] The present invention also relates to an application of the above parallel optimal transport gridless impact collision process simulation method, which is applied to the simulation of the impact collision process between two objects.

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

[0091] 1. The parallel optimal transport gridless impact collision process simulation method involved in the present invention establishes a topological connection relationship diagram of material points on the basis of the existing OTM method, decomposes the topological connection relationship diagram into material point subsets equal to the number of processors according to the number of parallel processors, determines the boundary material nodes and shared nodes in each material point subset, and then allocates the material point subsets to each parallel processor and establishes a shared data exchange table. Each processor obtains a local mass matrix, local node force and local node acceleration based on the allocated material point subset, and finally exchanges and assembles data according to the shared data exchange table to obtain a global mass matrix, global node force and global acceleration. This method combines graph topology analysis and distributed multi-process parallelization, can effectively use supercomputers to accelerate OTM simulations, and effectively saves computing time.

[0092] 2. The parallel optimal transport gridless impact collision process simulation method involved in the present invention is a gridless method for simulating impact collision problems at different speeds. It can overcome the drawbacks of grid entanglement or grid reconstruction caused by grids in traditional numerical simulation methods, and can accurately describe the problems such as large deformation, fracture, and material phase change experienced by materials in the impact collision process, especially in high-speed and ultra-high-speed impact collision processes. Theoretically, it has higher accuracy and its results can be infinitely close to analytical solutions. The parallel optimal transport gridless method for simulating impact collision processes is suitable for impact collision simulation under different conditions such as low-speed impact, medium-speed impact, high-speed impact and ultra-high-speed impact, and its calculation efficiency is greatly improved and its calculation stability is high. BRIEF DESCRIPTION OF THE DRAWINGS

[0093] FIG1 is a flow chart of a parallel optimal transport gridless impact collision process simulation method according to the present invention;

[0094] Figure 2 A schematic diagram of spatial discretization of a geometric model by combining material points with nodes;

[0095] Figure 3 A schematic diagram for establishing topological connections based on material point sets and dividing several material point sets based on topological connection relationships;

[0096] Figure 4 is a schematic diagram of a point-to-point contact algorithm in the present invention;

[0097] Figure 5 This is a diagram simulating the impact process of a steel ball hitting the first layer of aluminum plate in the application example;

[0098] Figure 6 The comparison diagrams are of the initial shape of the steel ball, the shape during impact with the first layer of aluminum plate, and the shape after impact in the application example;

[0099] Figure 7 This is a schematic diagram of the material response of the target plate when the steel ball hits the first aluminum plate in the application example;

[0100] Figure 8 This is a comparison chart between the hole size formed on the first aluminum plate after the steel ball hits the first aluminum plate and the experimental results;

[0101] Fig. 9 This is a schematic diagram of the contact range between the debris cloud formed after the steel ball hits the first aluminum plate and the second aluminum plate;

[0102] Fig.10 This example simulates the effect of copper flying pieces hitting a rigid plane to obtain an approximately linear acceleration effect. DETAILED DESCRIPTION

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

[0104] Example: Refer to the attached Figure 1 As shown, the present invention relates to a parallel optimal transport gridless impact collision process simulation method comprising the following steps:

[0105] S1. Setting a problem domain for solving the impact and collision problem, and then establishing a geometric model;

[0106] S2. Initialize nodes and material points: Discretize the geometry model in space, use the finite element mesh to discretize the geometry model into a set of node sets and a set of material point sets, initialize the node sets and material point sets, and discretize the time domain into several time steps. The generation method and initial position of material points and nodes can be determined by the user according to different algorithms, and the kinematic information of the calculation neighborhood, including displacement, velocity, acceleration and temperature, will be stored on the nodes; the physical information of the material, such as deformation, stress, internal parameters of the material, will be stored on the material points.

[0107] In the specific implementation process, for two-dimensional problems, the geometric model is divided into a first-order triangular unit grid, and for three-dimensional problems, the geometric model is divided into a first-order tetrahedral unit grid, which is used to initialize the material points and nodes in the OTM method. Figure 2 As shown, the hollow points represent nodes x a,k ( a is the node index number in the discrete domain, representing the node number. k represents the time step, the time step length and number of steps are controlled by the user), the solid triangle points represent the material points x p,k ( p is the material point index number in the discrete domain, representing the material point number. k represents the number of time steps).

[0108] S3. Establish a topological connection relationship diagram of material points, decompose the topological connection relationship diagram into material point subsets equal to the number of processors according to the number of parallel processors, and determine the boundary material nodes and shared nodes in each material point subset; Figure 3As shown, the material points corresponding to adjacent units in the original finite element grid are connected, and a topological connection relationship between the material points is established to form a graph, in which the nodes are material points and the edges are the topological connection relationships of the material points. At the same time, according to the number of parallel processors P, the topological structure graph is decomposed into P material point subsets, so that the number of disconnected edges in the topological structure graph is minimized, while ensuring that the number of material points in each material point subset is similar. All nodes in the neighborhood of the material point subset also belong to the material point set, thereby completing the discrete grouping of the OTM space. The material points that have a topological connection relationship with the material points on other material point subsets are defined as boundary material points, and all nodes in the neighborhood of the boundary material points are shared nodes. When the topological connection relationship graph is grouped, the processors to which the material point sets connected by the disconnected edges belong are the processors that need to communicate with each other.

[0109] S4. Allocate the material point subsets to each parallel processor and establish a shared data exchange table. The specific steps of establishing the shared data exchange table are:

[0110] S4.1. Calculate the neighborhood of each material point on the processor;

[0111] S4.2. Determine the boundary material points on each processor according to the obtained topological structure diagram;

[0112] S4.3. According to the boundary material points of each processor, determine the nodes in the neighborhood of each boundary material point to obtain the shared nodes;

[0113] S4.4. Save the shared node to the shared data exchange table to obtain the shared data exchange table;

[0114] The processor information and shared node information that need to communicate are stored in a shared data exchange table, and the shared data exchange table is stored on each processor.

[0115] S5. Each processor calculates the internal force based on the assigned material point subset, calculates the contact force of the nodes on the two objects in the impact collision problem through the contact algorithm, and obtains the local mass matrix, local node force and local node acceleration, which are specifically:

[0116] S5.1. Initialize nodes and material points, including initializing the time step k To 0, initialize the mass, velocity and nodal force of each node; calculate the shape function of the initial material point and shape function derivatives ;

[0117] The shape function used in the present invention is the local maximum entropy shape function. Given a discrete node set of a continuous medium x i ,node a The shape function is ,

[0118] in, N a Indicates that the index number is a The shape function of the node, x Represents the position of any node in the node set. x a Indicates that the index number is a The location of the node, β represents the Pareto optimization parameter, and Z The variables introduced with respect to node x must satisfy the following formula:

[0119] ,

[0120] ,

[0121] in, λ Represents the variable associated with the neighborhood search radius; the numerical solution is obtained through Newton-Raphson iteration. is a Pareto optimization parameter that measures the maximum entropy and the influence of other nodes in the neighborhood. In fact, in order to consider the influence of node distribution, a size parameter is used γ To control the shape function of LME, γ Satisfy the following formula: ,

[0122] in, h It is the characteristic scale of the node set, which can generally be set to the shortest distance of the calculated neighborhood;

[0123] The form of the shape function derivative after calculation is as follows:

[0124] ,

[0125] in, The form is as follows:

[0126] ,

[0127] in, Represents a variable λ Find the partial derivative, r represents the gradient vector of the objective function, N Represents the number of nodes in the neighborhood of the node, and the expression of the gradient vector is:

[0128] ;

[0129] S5.2. After the shape function is obtained, the deformation gradient at the material point is calculated. The deformation gradient at the material point is calculated using the following formula:

[0130] ,

[0131] In the formula, F p,k The index number is p The material point is k The deformation gradient of the time step, k For the k time steps, x p,k The index number is p The material point is k The position of the time step;

[0132] S5.3. Update the node position and calculate the updated material point deformation gradient;

[0133] The update formula of the node position is:

[0134] ,

[0135] In the formula, x a,k+1 For the k +The position of the node at 1 time step, v a,k The index number is a The node in k The speed of the movement in time steps, △t is the time step, a a,k The index number is a The node in k The acceleration of the time step;

[0136] The calculation formula of the node acceleration is:

[0137] ,

[0138] In the formula, f a is the nodal force of the node with index number a;

[0139] The update formula of the material point deformation gradient is:

[0140] ,

[0141] in, F p,k+1 For the k +1 deformation gradient of the material point at time step, φk→k+1 for k Step to k +1-step transfer mapping;

[0142] S5.4. Deformation gradient of material point F p , k+1 Decomposed into a pure rotation tensor R and a pure stretch tensor V , according to the pure stretch tensor V Calculate strain;

[0143] Among them, the material point deformation gradient F p , k+1 Decomposed into a pure rotation tensor R and a pure stretch tensor V Follow this formula:

[0144] .

[0145] According to the pure tensor V Calculate the material strain. The material strain under large deformation can be obtained by the following formula:

[0146] ,

[0147] in, represents material strain;

[0148] S5.5. In this embodiment, the stress of the material point is calculated using the J2 Power Law constitutive model. The stress of the material point is calculated based on the above material strain. The calculation formula for the yield stress of the material point is:

[0149] ,

[0150] in, σ y is the yield stress of the material point, σ 0 is the initial quasi-static yield stress of the material point, and is the effective plastic strain and strain rate, and For reference effective plastic strain and strain rate, T is the absolute temperature, n For the reinforcement index, m is the strain rate exponent, T 0 is the reference temperature, T m is the melting temperature, lis the thermal softening index; the plastic strain and plastic strain rate are calculated by the material strain Iterate and loop to obtain;

[0151] S5.6. Determine whether the material point is destroyed and cracks are generated. This embodiment uses the EigenFracture model to calculate the equivalent energy release rate of the material point. When the equivalent energy release rate of the material point is not lower than the critical energy release rate, it is determined that the material point is destroyed and cracks are generated; wherein, the calculation formula of the equivalent energy release rate of the material point is:

[0152] ,

[0153] In the formula, G ( x p ) is the index number p The equivalent energy release rate of the material point is ch p is the material point size, V ϵ is the material point volume, W is the elastic energy stored in the material point, x q is the material point position, For the average area, is the intrinsic strain, ▽ u is the displacement gradient, V q is the energy parameter.

[0154] When the calculated material point equivalent energy release rate is greater than or equal to the critical energy release rate that the material can withstand G c ( T,Z ), that is G ( X p )≥ G c ( T,Z ), it indicates that the material represented by the material point is damaged and a free surface is generated, which is manifested as the initiation or expansion of cracks, the material point fails, and its stress is zero; on the contrary, if the material represented by the material point is not damaged, the stress of the material point is interpolated to each node in its neighborhood to obtain the local nodal force caused by material deformation at each node. The nodal force caused by material deformation is calculated by the following formula:

[0155] ,

[0156] in, represents the nodal forces due to material deformation, M represents the number of material points in the node neighborhood, wrepresents the volume of a material point, P p Represents the first Piola-Kirchhoff stress of a material point;

[0157] S5.7. Calculate the tensile force on the boundary and the contact force occurring during collision. The local nodal force is obtained by adding the local nodal force caused by material deformation, the tensile force on the boundary and the contact force occurring during collision at each node.

[0158] ,

[0159] in, The first k time steps, t is the tensile force applied on the boundary, f c is the contact force.

[0160] The contact force is calculated using the point-to-point collision contact model, such as Figure 4 As shown in the figure, each node on the object is defined as a contact node set, the contact range of the point is set, and a contact pair is set between any two nodes from different contact node sets. The contact pair is in failure mode in the initial state. When the distance between two nodes from different contact node sets is less than the contact range, the contact pair between the two nodes is enabled, the contact depth of the contact pair is calculated, and the contact force generated by the contact pair is calculated by the contact depth and directly applied to the two nodes contained in the contact pair.

[0161] The calculation formula of the contact force is:

[0162] ,

[0163] In the formula, f c is the contact force, K p is the contact stiffness, is the contact depth of the contact pair, F 1 and F 2 is the equivalent contact force when two nodes are in contact. This formula indicates F 1 and F 2 The larger of the two is selected;

[0164] The calculation formula of the contact depth of the contact pair is:

[0165] ,

[0166] In the formula, R i and R j Represent the contact diameters of the two nodes, r i j Represents the distance between two nodes;

[0167] The equivalent contact forces of the two nodes are calculated using the following formulas:

[0168] ,

[0169] ,

[0170] In the formula, m' is the average mass of the two nodes, c is the speed of sound in the two objects, △t is the time step.

[0171] S6. Data exchange and assembly are performed according to the shared data exchange table, specifically:

[0172] S6.1. Each processor exchanges the local node forces of the shared nodes contained in the data exchange table. , local mass matrix , local acceleration Send to other processors, and receive the calculation results of the shared node on other processors;

[0173] S6.2. Each processor synchronizes and assembles its own data and the data received from other processors to obtain the global mass matrix, global node force, global acceleration and global velocity.

[0174] Wherein, the calculation formula of the global mass matrix is:

[0175] ,

[0176] In the formula, m a is the global mass matrix, is the local mass matrix calculated by processor number I, Q For shared nodes, C k To share data exchange tables, The local quality matrix obtained by exchanging nodes is shared by the remaining processors;

[0177] The calculation formula of the global nodal force is:

[0178] ,

[0179] In the formula, f a,k For the k The global nodal forces for each time step, is the node force calculated by processor numbered I at the kth time step, The node force obtained by exchanging the shared nodes of the remaining processors at the kth time step;

[0180] The calculation formula of the global acceleration is:

[0181] ,

[0182] The calculation formula of the global velocity is:

[0183] .

[0184] S7. Update material nodes and shared data exchange tables, including updating material point coordinates , update the material point neighborhood , Update material point shape function , Update the material point shape function derivative , use the updated neighborhood of each boundary material point to update the shared node set and add the new node to the shared data exchange table C k And delete the nodes that are not in the neighborhood of the boundary material point, and then reset the shared node data exchange table C k Then, enter the next time step to determine whether this time step is the last time step set. If not, return to S5 to enter the next iterative calculation; if so, output the calculation result of S6 to complete the impact collision process simulation, and obtain the final kinematic information and material physics information, and obtain the material points and nodes from t 0 →t n The physical information and dynamic data within a time period, including deformation, stress, density, displacement, velocity, acceleration, and temperature, complete the dynamic response analysis of the material.

[0185] Application example: The parallel optimal transport gridless impact collision process simulation method in the embodiment is applied to the simulation of the impact collision process between two objects. This application example takes the ultra-high-speed impact of a steel ball on two layers of aluminum plates as an example to illustrate the parallel optimal transport gridless method of the present invention for simulating the impact collision process.

[0186] This application example simulates the condition where a steel ball with a diameter of 3 mm hits two aluminum plates 100 mm apart at a speed of 6.09 km / s. The thickness of the aluminum plate is 5 mm. This case contains 1,368,048 material points, 832,682 degrees of freedom, and is calculated using 64 CPUs on a distributed high-performance computing server. During the collision, the steel ball penetrates the first aluminum plate, such as Figure 5 In (a)-(c), as the plastic deformation increases, the temperature gradually rises, and a flange bulge gradually forms on the surface of the aluminum plate. As the penetration depth increases, the steel ball gradually changes from a spherical shape to an oblate spherical shape, as shown in Figure 6 As shown in (a), (b), and (c), as the penetration depth increases, the steel ball begins to fail due to the increased force, and the back of the target plate cracks, forming a large petal-shaped hole edge, such as Figure 5 In (d)-(h), the steel ball finally leaves the surface of the first aluminum plate, as shown in Figure 5 As shown in (i) in .

[0187] When the steel ball and the aluminum plate begin to contact, the bottom of the steel ball is subjected to force, and it softens after plastic deformation and gradually becomes an elliptical shape. As it penetrates the aluminum plate, the deformation of the surface in contact with the aluminum plate is hindered, and the middle of the steel ball gradually concave to form an oblate spherical shape, forming a certain number of failure material points, such as Figure 6 As shown by the white translucent material dots in the middle, as the penetration depth increases until the first plate is penetrated, part of the material on the projectile adheres to the flat plate, and a small part of the material forms a reverse debris cloud.

[0188] After the steel ball penetrates the first layer of the plate, as the plastic deformation increases, the temperature of the deformed part rises, the material softens significantly, and a flange bulge is formed on the projectile incident surface. As the penetration increases, the steel ball becomes flatter, and the stress wave reflected back from the lower surface causes the tensile stress to increase, and the plate forms a layer crack near the lower surface. Therefore, the failure range of the bottom part is expanded, and the penetration aperture is larger than the flange aperture of the upper surface. Figure 7 is the material response of the first aluminum plate in the simulation, Figure 8 This is the shape of the hole after the projectile penetrates the first layer of aluminum plate. Figure 8 In the figure, (a) is the hole on the surface of the first layer of board, with a diameter of about 12 mm and a flange height of about 4.4 mm; (b) is the penetration diameter of the bottom surface of the first layer of board, which is about 17.6 mm; (d) is the perforation morphology in the experiment.

[0189] Fig. 9is the contact range between the debris cloud generated by the steel ball hitting the first target plate and the second target plate. When the time is 31.48 μs, it can be seen that the range of action of the projectile and the first target plate debris cloud on the second target plate is about 54 mm. Assuming that after the projectile penetrates the first target plate, there is no mutual interference between the debris clouds and the motion trajectory remains a straight line, then after the range of action is expanded by 2.5 times, that is, when the distance between the two target plates is 100 mm, the range of action is about 135 mm; on this basis, the distance between the two target plates is adjusted to 100 mm for verification. Under this condition, the range of action of the projectile and the first target plate debris cloud on the second target plate is about 136 mm.

[0190] Effect example: This is a simulation of a copper flying piece hitting a rigid plane at a speed of 270m / s. This example is to test the acceleration performance of the parallel optimal transport gridless method proposed by the present invention for simulating the impact collision process. This implementation case includes 1,116,024 material points and 683,772 degrees of freedom. The calculation is performed on a distributed high-performance computing server, using 1, 2, 4, 8, 16, 32, and 64 CPUs respectively, and an effect similar to linear acceleration is obtained. The parallel efficiency is very high, such as Fig.10 shown.

[0191] 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 parallel optimal transport gridless impact collision process simulation method, characterized by: It includes the following steps: S1. Setting a problem domain for solving the impact and collision problem, and then establishing a geometric model; S2. Discretize the geometric model in space, discretize the geometric model into a set of node sets and a set of material point sets using a finite element grid, initialize the node sets and the material point sets, store initial kinematic information in the nodes, and store initial material physical information in the material point sets; the kinematic information includes displacement, velocity, acceleration and temperature, and the material physical information includes deformation, stress and internal parameters of the material; discretize the time domain into a number of time steps; S3. Establish a topological connection relationship diagram of material points, decompose the topological connection relationship diagram into material point subsets equal to the number of processors according to the number of parallel processors, and determine the boundary material nodes and shared nodes in each material point subset; S4. Assign the material point subsets to each parallel processor and establish a shared data exchange table. The specific steps are as follows: S4.

1. Calculate the neighborhood of each material point on the processor; S4.

2. Determine the boundary material points on each processor according to the obtained topological structure diagram; S4.

3. According to the boundary material points of each processor, determine the nodes in the neighborhood of each boundary material point to obtain the shared nodes; S4.

4. Save the shared node to the shared data exchange table to obtain the shared data exchange table; S5. Each processor calculates the internal force based on the assigned material point subset, calculates the contact force of the nodes on the two objects in the impact collision problem through the contact algorithm, and obtains the local mass matrix, local node force and local node acceleration. The specific steps include: S5.

1. Initialize nodes and material points, and calculate shape functions and shape function derivatives of initial nodes; The shape function adopts the local maximum entropy shape function, and the calculation formula of the shape function is: ， in, N a Indicates that the index number is a The shape function of the node, x Represents the position of any node in the node set. x a Indicates that the index number is a The location of the node, β represents the Pareto optimization parameter, and Z For the introduced node x The relevant variables must satisfy the following formula: ， ， in, λ Represents variables related to neighborhood search radius; Shape function derivative ▽N a ( x ) is calculated as follows: ， in, J The variables introduced with respect to node x follow the following form: ， in, represents the partial derivative of the variable λ, r represents the gradient vector of the objective function, N Represents the number of nodes in the neighborhood of the node, and the expression of the gradient vector is: ； S5.

2. Calculate the deformation gradient at the material point using the formula: ， In the formula, F p,k The index number is p The material point is k The deformation gradient of the time step, k For the k time steps, x p,k The index number is p The material point is k The position of the time step; S5.

3. Update the node position and calculate the updated material point deformation gradient. The update formula of the node position is: ， In the formula, x a,k+1 For the k +The position of the node at 1 time step, v a,k The index number is a The node in k The speed of the movement in time steps, △t is the time step, a a,k The index number is a The node in k The acceleration of the time step; The updating formula of the material point deformation gradient is: ， in, F p,k+1 For the k +1 deformation gradient of the material point at time step, φ k→k+1 For k steps to k +1-step transfer mapping; S5.

4. Deformation gradient of material point F p,k+1 Decomposed into a pure rotation tensor R and a pure stretch tensor V , according to the pure stretch tensor V Calculate strain; Deformation gradient of material point F p,k+1 The formula for decomposing into a pure rotation tensor R and a pure stretch tensor V is: ， According to the pure tensor V The calculation formula for calculating material strain is: ， in, represents material strain; S5.

5. Based on the above material strain, the stress of the material point is calculated. The calculation formula of the yield stress of the material point is: ， in, σ y is the yield stress of the material point, σ 0 is the initial quasi-static yield stress of the material point, and is the effective plastic strain and strain rate, and For reference effective plastic strain and strain rate, T is the absolute temperature, n For the reinforcement index, m is the strain rate exponent, T 0 is the reference temperature, T m is the melting temperature, l is the thermal softening index; the plastic strain and plastic strain rate are calculated by the material strain Iterate and loop to obtain; S5.

6. Determine whether the material point is damaged and cracks are generated. If it is not damaged, interpolate the stress of the material point to each node in its neighborhood to obtain the local node force caused by material deformation at each node; the node force caused by material deformation is calculated by the following formula: ， in, represents the nodal forces due to material deformation, M represents the number of material points in the node neighborhood, w represents the volume of a material point, P p Represents the first Piola-Kirchhoff stress of a material point; S5.

7. Calculate the tension applied on the boundary and the contact force occurring during the collision, and calculate the local node force of each processor; S6. perform data exchange and assembly according to the shared data exchange table, and calculate the global mass matrix, global node force and global acceleration; S7. Update the material nodes and shared data exchange table, enter the next time step, and determine whether this time step is the last time step set. If not, return to S5 to enter the next iterative calculation; if so, output the calculation result of S6, complete the impact collision process simulation, and obtain the final kinematic information and material physics information.

2. The parallel optimal transport gridless impact collision process simulation method according to claim 1, characterized in that: The S5.1 node and material point initialization includes initializing the mass, velocity and node force of each node; the calculation formula of the node acceleration in S5.2 is: ， In the formula, f a The index number is a Nodal forces at the nodes.

3. The parallel optimal transport gridless impact collision process simulation method according to claim 1, characterized in that: When determining whether a material point is damaged and cracks are generated in S5.6, the EigenFracture model is used to calculate the equivalent energy release rate of the material point. When the equivalent energy release rate of the material point is not lower than the critical energy release rate, it is determined that the material point is damaged and cracks are generated; The calculation formula of the material point equivalent energy release rate is: ， In the formula, G ( x p ) is the index number p The equivalent energy release rate of the material point is ch p is the material point size, V ϵ is the material point volume, W is the elastic energy stored in the material point, x q is the material point position, For the average area, is the intrinsic strain, ▽ u is the displacement gradient, V q is the energy parameter.

4. The parallel optimal transport gridless impact collision process simulation method according to claim 1, characterized in that: The calculation formula for the contact force generated during the S5.7 collision is: ， In the formula, f c is the contact force, K p is the contact stiffness, is the contact depth of the contact pair, F 1 and F 2 is the equivalent contact force calculated when two nodes are in contact; The calculation formula of the contact depth of the contact pair is: ， In the formula, R i and R j Represent the contact diameters of the two nodes, r i j Represents the distance between two nodes; The equivalent contact forces of the two nodes are calculated using the following formulas: ， ， In the formula, m’ is the average mass of the two nodes, c is the speed of sound in the two objects, △t is the time step.

5. The parallel optimal transport gridless impact collision process simulation method according to claim 1, characterized in that: In S5.7, the local node force of each node caused by material deformation, the tension on the boundary and the contact force occurring during collision are added together to obtain the local node force, that is: ， in, The first k The local nodal forces at each time step, t is the tensile force applied on the boundary, f c is the contact force.

6. The parallel optimal transport gridless impact collision process simulation method according to claim 1, characterized in that: The specific steps of S6 performing data exchange and assembly according to the shared data exchange table and calculating the global mass matrix, global node force and global acceleration are as follows: S6.

1. Each processor sends the local node force, local mass matrix, and local acceleration of the shared node contained in the data exchange table to other processors, and receives the calculation results of the shared node in other processors; S6.

2. Each processor synchronizes and assembles its own data and the data received from other processors to obtain the global mass matrix, global node force, global acceleration and global velocity.

7. The parallel optimal transport gridless impact collision process simulation method according to claim 6, characterized in that: The calculation formula of the S6.2 global mass matrix is: ， In the formula, m a,k For the k The global mass matrix for time steps, For the k The local mass matrix calculated by processor number I at time step, Q For shared nodes, C k To share data exchange tables, No. k The local mass matrix obtained by exchanging nodes is shared by other processors at the time step; The calculation formula of the global nodal force is: ， In the formula, f a,k For the k The global nodal forces for the time step, For the k The node force calculated by processor number I at time step, For the k The remaining processors share the node forces obtained by node exchange during the time step; The calculation formula of the global acceleration is: ， The calculation formula of the global velocity is: 。 8. Application of the parallel optimal transport gridless impact collision process simulation method according to claim 1, which is applied to the simulation of the impact collision process of two objects.

Images