CN118965832A - A parallel optimal transport meshless impact collision process simulation method and its application - Google Patents
A parallel optimal transport meshless impact collision process simulation method and its application Download PDFInfo
- Publication number
- CN118965832A CN118965832A CN202411442797.0A CN202411442797A CN118965832A CN 118965832 A CN118965832 A CN 118965832A CN 202411442797 A CN202411442797 A CN 202411442797A CN 118965832 A CN118965832 A CN 118965832A
- Authority
- CN
- China
- Prior art keywords
- node
- nodes
- material point
- time step
- formula
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种并行最优输运无网格的冲击碰撞过程模拟方法，属于计算力学与数值仿真技术领域，其包括以下步骤：建立几何模型；对几何模型进行空间离散，将时间域离散为若干时间步；建立物质点的拓扑连接关系图并分解，确边界物质节点和共享节点；物质点子集分别分配给各并行处理器，并建立共享数据交换表；各处理器内力计算和接触作用力，得到局部质量矩阵、局部节点力及局部节点加速度；共享数据交换表进行数据交换和组装，得到全局计算结果；进入下一个时间步，进行迭代计算或输出计算结果。本发明结合了图拓扑结构分析和分布式多进程并行化，能够有效的利用超级计算机对OTM模拟进行加速，有效节省计算时间。
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 a local mass matrix, local node forces and local node accelerations; the shared data exchange table performs data exchange and assembly, and obtains a global calculation result; enters the next time step, performs iterative calculation or outputs the calculation result. The present invention combines graph topology structure analysis and distributed multi-process parallelization, and can effectively use supercomputers to accelerate OTM simulations, effectively saving calculation time.
Description
Technical Field
The invention belongs to the technical field of computational mechanics and numerical simulation, and particularly relates to a parallel optimal transport gridless impact collision process simulation method.
Background
The impact collision problem is a hot spot problem in the research of the fields of aerospace and the like, relates to complex characteristics such as extremely large deformation, dynamic crack initiation and expansion, geometric distortion, material fission and the like, in particular to the high-speed and ultra-high-speed impact collision problem with the speed of more than 1 km/s, even the complex processes such as material phase change and the like can be involved, and the numerical simulation is difficult to carry out by using a traditional numerical method based on grids. Traditional numerical methods based on lagrangian and euler grids, such as finite element methods, finite volume methods, finite difference methods, etc., have limited their application in practical engineering problems involving high-speed impact collisions due to their inherent grid characteristics. When solving the high-speed and ultra-high-speed impact collision problem, the Lagrangian-based grid solving method is easy to encounter the problems of grid distortion, grid entanglement and the like caused by grid large deformation, and needs to be solved by a method that errors are easy to be introduced and calculation load is greatly increased through grid reconstruction and the like; the euler mesh-based solution has difficulty in dealing with material dynamic response and crack propagation problems related to loading history and in determining free surfaces, deformable boundaries and interfaces for moving materials. Therefore, in order to solve the above-mentioned difficulties of the conventional grid-based numerical method, a new generation numerical calculation method, a grid-free method, has been widely and rapidly developed. The goal of the grid-less approach is to eliminate the dependency on the grid completely (at least in part) by discretizing the entire problem domain in a node-wise fashion. The main advantages are as follows: ① The calculation process does not involve grids, has no grid distortion, grid reconstruction and other difficulties, and is very suitable for processing the problem of large deformation; ② The implementation of H-type and P-type adaptation becomes relatively simple; ③ The shape function of the continuity of any order can be easily constructed; ④ Convergence results tend to be much better than grid-based approaches; ⑤ The problem of volume self-locking can be alleviated by adjusting the expansion parameters of the shape function.
The optimal transport gridless method (Optimal Transportation Meshfree method, OTM method) effectively utilizes the spatial dispersion of material points, the local maximum entropy interpolation function and the optimal transport theoretical time dispersion, avoids the instability of calculated tensile stress, meets the 'weak' Cronecker attribute, and can accurately apply displacement boundary conditions. In addition, the OTM method can dynamically adjust the connection relation between the object points and the nodes, can automatically detect the dynamic contact surface without other complex interface tracking technologies, and has great advantages in solving the high-speed/ultra-high-speed collision problem. However, the OTM method requires a large computational load in solving the impact collision problem, which results in low computational efficiency and slow simulation speed. How to improve the calculation efficiency of the OTM method in solving the impact collision process is one of the problems that the OTM method needs to solve in a major way.
Disclosure of Invention
The invention aims to provide a parallel optimal transport gridless impact collision process simulation method, which aims to solve the problem of insufficient calculation efficiency in the impact collision simulation problem of the conventional OTM method.
In order to achieve the above object, the technical scheme of the present invention is as follows:
The invention relates to a parallel optimal transport gridless impact collision process simulation method, which comprises the following steps:
S1, setting a problem domain for solving an impact collision problem, and further establishing a geometric model;
s2, performing space dispersion on the geometric model, dispersing the geometric model into a group of node sets and a group of substance node sets by utilizing a finite element grid, initializing the node sets and the substance node sets, and dispersing a time domain into a plurality of time steps;
S3, establishing a topological connection relation diagram of object points, decomposing the topological connection relation diagram into substance point subsets equal to the number of the parallel processors, and determining boundary substance nodes and shared nodes in each substance point subset;
S4, respectively distributing the substance point subsets to each parallel processor, and establishing a shared data exchange table;
S5, each processor calculates internal force based on the distributed substance point subsets, and calculates contact acting force of nodes on two objects in the impact collision problem through a contact algorithm to obtain a local mass matrix, local node force and local node acceleration;
s6, carrying out data exchange and assembly according to the shared data exchange table, and calculating a global quality matrix, global node force and global acceleration;
S7, updating the material node and the shared data exchange table, entering the next time step, judging whether the time step is the last time step set, and if not, returning to S5 to enter the next iterative calculation; if yes, outputting the calculation result of the step S6, completing the impact collision process simulation, and obtaining the final kinematic information and the material physical information.
Preferably, in the step S2, the initial kinematic information is stored in the node, and the initial physical information of the material is stored in the object point set; the kinematic information comprises displacement, speed, acceleration and temperature, and the physical information of the material comprises deformation, stress and internal parameters of the material;
The specific step of S5 comprises the following steps:
s5.1, initializing nodes and object points, and calculating a shape function and a shape function derivative of the initial nodes;
the shape function adopts a local maximum entropy shape function, and the calculation formula of the shape function is as follows:
，
Wherein N a represents a shape function of a node with index number a, x represents a position of any one node in the node set, x a represents a position of the node with index number a, beta represents a Paret optimization parameter, And Z is an introduced variable related to node x, satisfying the following equation:
，
，
Wherein λ represents a variable related to a neighborhood search radius;
The form of the derivative of the shape function #, N a (x), is calculated as follows:
，
where J is an introduced variable related to node x, following the form:
，
wherein, The method is characterized by comprising the steps of calculating a bias derivative of a variable lambda, wherein r represents a gradient vector of an objective function, N represents the number of nodes in the neighborhood of the node, and the expression of the gradient vector is as follows:
；
s5.2, calculating deformation gradient on object particles, wherein the calculation formula is as follows:
，
In the formula, F p,k is the deformation gradient of a substance point with an index number p in the kth time step, k is the kth time step, and x p,k is the position of the substance point with the index number p in the kth time step;
S5.3, updating the node positions and calculating the deformation gradient of the updated substance points, wherein the updating formula of the node positions is as follows:
，
In the formula, x a,k+1 is the position of a node at the (k+1) th time step, v a,k is the motion speed of the node with the index number a at the (k) th time step, Δt is the time step, and a a,k is the acceleration of the node with the index number a at the (k) th time step;
The update formula of the deformation gradient of the object point is as follows:
，
wherein F p,k+1 is the deformation gradient of the material point in the (k+1) th time step, and phi k→k+1 is the transmission mapping from the k step to the k+1 step;
S5.4, decomposing the material point deformation gradient F p,k+1 into a pure rotation tensor R and a pure stretching tensor V, and calculating strain according to the pure stretching tensor V;
the formula for decomposing the material point deformation gradient F p,k+1 into a pure rotation tensor R and a pure stretching tensor V is:
，
The calculation formula for calculating the material strain according to the pure tensile tensor V is:
，
wherein, Representing material strain;
s5.5, calculating stress of material points based on the strain of the material, wherein a calculation formula of yield stress of the material points is as follows:
，
Where σ y is the yield stress of the object point, σ 0 is the initial quasi-static yield stress of the object point, AndFor an effective plastic strain and strain rate,AndFor reference effective plastic strain and strain rate, T is absolute temperature, n is the strengthening index, m is the strain rate index, T 0 is the reference temperature, T m is the melting temperature, and l is the heat softening index; the plastic strain and the plastic strain rate are calculated to obtain the material strainObtaining by iterative loop;
S5.6, judging whether the material points are damaged and generate cracks, if the material points are not damaged, interpolating the stress of the material points to each node in the neighborhood of the material points to obtain the local node force of each node caused by material deformation; the nodal force caused by deformation of the material is calculated by the following formula:
，
wherein, Representing node force caused by material deformation, M representing the number of object points in the neighborhood of the node, w representing the volume of the object point, and P p representing the first Piola-Kirchhoff stress of the object point;
s5.7, calculating the tensile force applied on the boundary and the contact acting force generated during collision, and calculating to obtain the local node force of each processor.
Preferably, the S5.1 initializing nodes and substance points includes initializing mass, speed and node force of each node; the calculation formula of the acceleration of the S5.2 node is as follows:
，
In the formula, f a is the node force of the node with index a.
Preferably, when the S5.6 determines whether the material point is damaged and a crack is generated, calculating an equivalent energy release rate of the material point by adopting a EigenFracture model, and when the equivalent energy release rate of the material point is not lower than a critical energy release rate, determining that the material point is damaged and a crack is generated;
the calculation formula of the equivalent energy release rate of the substance points is as follows:
，
In the formula, G (x p) is the equivalent energy release rate of an object point with index number p, ch p is the size of the object point, V ϵ is the volume of the object point, W is the elastic energy stored by the object point, x q is the point size of the object point, In order to be a field of averaging the number of the devices,Let u be the displacement gradient and V q be the energy parameter.
Preferably, the calculation formula of the contact acting force generated during the S5.7 collision is as follows:
，
Where f c is the contact force, K p is the contact stiffness, For the contact depth of the contact pair, F 1 and F 2 are equivalent contact forces calculated when two nodes are in contact;
the calculation formula of the contact depth of the contact pair is as follows:
，
In the formula, R i and R j respectively represent the contact diameters of two nodes, and R i j represents the distance between the two nodes;
the equivalent contact force of the two nodes is calculated by the following formulas:
，
，
in the formula, m' is the average mass of two nodes, c is the speed of sound in two objects, and Δt is the time step.
Preferably, in S5.7, the local node force caused by the deformation of the material, the tensile force on the boundary, and the contact force occurring at the time of collision of each node are added to obtain the local node force, that is:
，
wherein, The local node force for the kth time step calculated for processor number I, t is the applied pull force on the boundary and f c is the contact force.
Preferably, the step S6 of exchanging and assembling data according to the shared data exchange table, and calculating the global quality matrix, the global node force and the global acceleration comprises the following specific steps:
S6.1, each processor transmits local node force, local mass matrix and local acceleration of the shared node contained in the data exchange table to other processors, and receives calculation results of the shared node in the other processors;
s6.2, each processor performs data synchronization and data assembly on own data and received data of other processors to obtain a global quality matrix, global node force, global acceleration and global speed.
Preferably, the calculation formula of the S6.2 global quality matrix is:
，
In the formula, m a,k is the global quality matrix of the kth time step, The local quality matrix calculated for the processor numbered I at the kth time step, Q being the shared node, C k being the shared data exchange table,The local quality matrix obtained by exchanging the nodes is shared by other processors in the kth time step;
The calculation formula of the global node force is as follows:
，
in the formula, f a,k is the global node force of the kth time step, The node force calculated for the processor numbered I at the kth time step,Sharing node force obtained by node exchange for other processors in the kth time step;
the calculation formula of the global acceleration is as follows:
，
the calculation formula of the global speed is as follows:
。
preferably, the specific step of establishing the shared data exchange table in S4 is:
s4.1, calculating the neighborhood of each object point on the processor;
s4.2, determining boundary material points on each processor according to the obtained topological structure diagram;
S4.3, determining nodes in the neighborhood of each boundary material point according to the boundary material points of each processor to obtain shared nodes;
And S4.4, storing the shared node into a shared data exchange table to obtain the shared data exchange table.
The invention also relates to application of the grid-free impact collision process simulation method for parallel optimal transportation, which is applied to simulation of the impact collision process of two objects.
Compared with the prior art, the technical scheme provided by the invention has the following beneficial effects:
1. The invention relates to a parallel optimal transport gridless impact collision process simulation method, which is characterized in that a topological connection relation diagram of object points is established on the basis of the existing OTM method, the topological connection relation diagram is decomposed into substance point subsets equal to processors according to the number of parallel processors, boundary substance nodes and shared nodes in each substance point subset are determined, then the substance point subsets are respectively distributed to each parallel processor, a shared data exchange table is established, each processor obtains a local mass matrix, local node force and local node acceleration based on the distributed substance point subsets, and finally data exchange and assembly are carried out according to the shared data exchange table, and global mass matrix, global node force and global acceleration are obtained. The method combines graph topological structure analysis and distributed multi-process parallelization, can effectively utilize the supercomputer to accelerate OTM simulation, and effectively saves calculation time.
2. The parallel optimal transport gridless impact process simulation method is a gridless method for simulating impact problems at different speeds, can overcome the defects of grid entanglement or grid reconstruction and the like caused by grids in the traditional numerical simulation method, can accurately describe the problems of oversized deformation, fracture crushing, material phase change and the like of materials in the impact process, especially in the high-speed and ultra-high-speed impact process, has higher theoretical precision, can be infinitely close to an analytic solution as a result, is suitable for impact simulation under different conditions such as low-speed impact, medium-speed impact, high-speed impact and ultra-high-speed impact, and has greatly improved calculation efficiency and high calculation stability.
Drawings
FIG. 1 is a flow chart of a parallel optimal transport gridless impact collision process simulation method in accordance with the present invention;
FIG. 2 is a schematic illustration of spatial discretization of a geometric model by combining object points with nodes;
FIG. 3 is a schematic diagram of establishing topological links according to material point sets and dividing a plurality of material point sets according to topological links;
FIG. 4 is a schematic diagram of a point-to-point contact algorithm according to the present invention;
FIG. 5 is a diagram showing an impact process of a simulated steel ball on a first aluminum plate in an application example;
FIG. 6 is a graph showing the initial form of the steel ball in the application example, the form of the steel ball after striking the first aluminum plate, and the form after striking;
FIG. 7 is a schematic diagram of the material response of the target plate during the impact of the steel ball on the first aluminum plate in the application example;
FIG. 8 is a graph comparing the size of a hole formed in a first aluminum plate after a steel ball impacts the first aluminum plate with experimental results;
FIG. 9 is a schematic view showing the contact range between the broken cloud formed after the steel ball impacts the first aluminum plate and the second aluminum plate in the application example;
Fig. 10 shows the effect of the simulated copper flyer impact stiffness plane to approximate linear acceleration.
Detailed Description
The invention will be further understood by reference to the following examples which are given to illustrate the invention but are not intended to limit the scope of the invention.
Examples: referring to fig. 1, the method for simulating the parallel optimal transport gridless impact collision process comprises the following steps:
S1, setting a problem domain for solving an impact collision problem, and further establishing a geometric model;
S2, initializing nodes and substance points: and performing space dispersion on the geometric model, dispersing the geometric model into a group of node sets and a group of substance node sets by utilizing a finite element grid, initializing the node sets and the substance node sets, and dispersing the time domain into a plurality of time steps. The generation mode and initial position of the object points and the nodes can be determined by a user according to different algorithms, and the kinematic information of the neighborhood is calculated, wherein the kinematic information comprises displacement, speed, acceleration, temperature and the like and is stored on the nodes; physical information of the material, such as deformation, stress, material internal parameters, etc., will be stored on the object point.
In the specific implementation process, for the two-dimensional problem, the geometric model is divided into first-order triangle unit grids, for the three-dimensional problem, the geometric model is divided into first-order tetrahedron unit grids, and object points and nodes in the OTM method are initialized by using the grids. As shown in fig. 2, the open dots represent node x a,k (a represents the index of the node in the discrete domain, represents the number of nodes, k represents the number of time steps, time steps and number of steps are controlled by the user), and the solid triangle dots represent object dot x p,k (p represents the index of the object dot in the discrete domain, represents the number of object dots, and k represents the number of time steps).
S3, establishing a topological connection relation diagram of object points, decomposing the topological connection relation diagram into substance point subsets equal to the number of the parallel processors, and determining boundary substance nodes and shared nodes in each substance point subset; as shown in fig. 3, substance points corresponding to adjacent cells in the original finite element mesh are connected, a topological connection relationship between the substance points is established, a graph is formed, the nodes of the graph are the substance points, and edges of the graph are the topological connection relationship of the substance points; meanwhile, according to the number P of parallel processors, the topological structure diagram is decomposed into P object point subsets, so that the broken edge number of the topological structure diagram is minimum, and meanwhile, the number of the object points in each object point subset is guaranteed to be similar; nodes in all substance point neighbors in the substance point subset also belong to the substance point set, so that the OTM space discrete grouping is completed; defining material points with topological connection relation with material points on other material point subsets as boundary material points, wherein nodes in all boundary material point neighbors are shared nodes, and processors to which the broken edges are connected when the topological connection relation diagram is grouped are processors which need to communicate with each other.
S4, respectively distributing the substance point subsets to each parallel processor, and establishing a shared data exchange table, wherein the specific steps of establishing the shared data exchange table are as follows:
s4.1, calculating the neighborhood of each object point on the processor;
s4.2, determining boundary material points on each processor according to the obtained topological structure diagram;
S4.3, determining nodes in the neighborhood of each boundary material point according to the boundary material points of each processor to obtain shared nodes;
s4.4, storing the shared node into a shared data exchange table to obtain the shared data exchange table;
The processor information and the shared node information to be communicated are stored in a shared data exchange table stored on each processor.
S5, each processor calculates internal force based on the distributed substance point subsets, and calculates contact acting force of nodes on two objects in the impact collision problem through a contact algorithm to obtain a local mass matrix, local node force and local node acceleration, wherein the method specifically comprises the following steps:
S5.1, initializing nodes and substance points, wherein the initialization time step k is 0, and initializing the quality, speed and node force of each node; calculating shape function of initial material point Sum function derivative；
The shape function employed in the present invention is a local maximum entropy shape function. Given a set of discrete nodes x i of a continuous medium, the shape function of node a is，
Wherein N a represents a shape function of a node with index number a, x represents a position of any one node in the node set, x a represents a position of the node with index number a, beta represents a Paret optimization parameter,And Z is an introduced variable related to node x, satisfying the following equation:
，
，
wherein λ represents a variable related to a neighborhood search radius; the numerical solution is obtained by newton-laprison iteration. Is a Pareto optimization parameter used to measure the maximum entropy and the influence of other nodes in the neighborhood. In practice, to take into account the influence of the node distribution, a size parameter γ is used to control the shape function of the LME, γ satisfying the following equation:，
wherein h is a feature scale of the node set, and can be generally set as the shortest distance of the calculation neighborhood;
The form of the derivative of the shape function is calculated as follows:
，
wherein, The form of (2) is as follows:
，
wherein, The method is characterized by comprising the steps of calculating a bias derivative of a variable lambda, wherein r represents a gradient vector of an objective function, N represents the number of nodes in the neighborhood of the node, and the expression of the gradient vector is as follows:
；
S5.2, after the shape function is obtained, calculating the deformation gradient on the object points, wherein the deformation gradient on the object points is calculated by using the following formula:
，
In the formula, F p,k is the deformation gradient of a substance point with an index number p in the kth time step, k is the kth time step, and x p,k is the position of the substance point with the index number p in the kth time step;
s5.3, updating the node positions and calculating the updated substance point deformation gradient;
The update formula of the node position is:
，
In the formula, x a,k+1 is the position of a node at the (k+1) th time step, v a,k is the motion speed of the node with the index number a at the (k) th time step, Δt is the time step, and a a,k is the acceleration of the node with the index number a at the (k) th time step;
the calculation formula of the acceleration of the node is as follows:
，
In the formula, f a is the node force of the node with the index number a;
The update formula of the deformation gradient of the object point is as follows:
，
wherein F p,k+1 is the deformation gradient of the material point in the (k+1) th time step, and phi k→k+1 is the transmission mapping from the k step to the k+1 step;
S5.4, decomposing the material point deformation gradient F p,k+1 into a pure rotation tensor R and a pure stretching tensor V, and calculating strain according to the pure stretching tensor V;
Wherein decomposing the point deformation gradient F p,k+1 into a pure rotation tensor R and a pure stretching tensor V follows the following formula:
。
The material strain is calculated according to the pure tensile tensor V, and the material strain under large deformation can be obtained by the following formula:
，
wherein, Representing material strain;
S5.5, in the embodiment, stress of the material points is calculated by adopting a J2 Power Law constitutive model, stress of the material points is calculated based on the strain of the material points, and a calculation formula of yield stress of the material points is as follows:
，
Where σ y is the yield stress of the object point, σ 0 is the initial quasi-static yield stress of the object point, AndFor an effective plastic strain and strain rate,AndFor reference effective plastic strain and strain rate, T is absolute temperature, n is the strengthening index, m is the strain rate index, T 0 is the reference temperature, T m is the melting temperature, and l is the heat softening index; the plastic strain and the plastic strain rate are calculated to obtain the material strainObtaining by iterative loop;
S5.6, judging whether the material points are damaged and cracks are generated, in the embodiment, calculating the equivalent energy release rate of the material points by adopting a EigenFracture model, and judging that the material points are damaged and cracks are generated when the equivalent energy release rate of the material points is not lower than the critical energy release rate; the calculation formula of the equivalent energy release rate of the substance points is as follows:
，
In the formula, G (x p) is the equivalent energy release rate of an object point with index number p, ch p is the size of the object point, V ϵ is the volume of the object point, W is the elastic energy stored by the object point, x q is the point size of the object point, In order to be a field of averaging the number of the devices,Let u be the displacement gradient and V q be the energy parameter.
When the calculated equivalent energy release rate of the material point is greater than or equal to the critical energy release rate G c (T, Z) which can be borne by the material, namely G (X p)≥Gc (T, Z)), the material represented by the material point is damaged and generates a free surface, namely the material point is in failure of the material point, the stress is zero, otherwise, if the material represented by the material point is not damaged, the stress of the material point is interpolated to each node in the neighborhood of the material point, and the local node force of each node caused by the material deformation is obtained, wherein the node force caused by the material deformation is calculated by the following formula:
，
wherein, Representing node force caused by material deformation, M representing the number of object points in the neighborhood of the node, w representing the volume of the object point, and P p representing the first Piola-Kirchhoff stress of the object point;
S5.7, calculating the tensile force on the boundary and the contact acting force generated during collision, and adding the local node force of each node caused by material deformation, the tensile force on the boundary and the contact acting force generated during collision to obtain the local node force;
，
wherein, The kth time step calculated for processor number I, t is the applied tension on the boundary and f c is the contact force.
The calculation of the contact acting force adopts a point-to-point collision contact model, as shown in fig. 4, the nodes on each object are defined as a contact node set, the contact range of the points is set, contact pairs are arranged between any two nodes from different contact node sets, and the contact pairs are failure modes in the initial state. When the distance between two nodes from different contact node sets is smaller 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 acting force generated by the contact pair is calculated through the contact depth and is directly applied to the two nodes contained in the contact pair.
The calculation formula of the contact acting force is as follows:
，
Where f c is the contact force, K p is the contact stiffness, For the contact depth of the contact pair, F 1 and F 2 are equivalent contact forces when two nodes are in contact, and the formula indicates that both F 1 and F 2 are selected larger;
the calculation formula of the contact depth of the contact pair is as follows:
，
In the formula, R i and R j respectively represent the contact diameters of two nodes, and R i j represents the distance between the two nodes;
the equivalent contact force of the two nodes is calculated by the following formulas:
，
，
in the formula, m' is the average mass of two nodes, c is the speed of sound in two objects, and Δt is the time step.
S6, carrying out data exchange and assembly according to the shared data exchange table, wherein the data exchange table specifically comprises the following steps:
S6.1, each processor uses the local node force of the shared node contained in the data exchange table Local mass matrixLocal accelerationSending the result to other processors and receiving the calculation result of the sharing node in the other processors;
s6.2, each processor performs data synchronization and data assembly on own data and received data of other processors to obtain a global quality matrix, global node force, global acceleration and global speed.
The calculation formula of the global quality matrix is as follows:
，
In the formula, m a is a global quality matrix, The local quality matrix calculated for the processor numbered I, Q is the shared node, C k is the shared data exchange table,Exchanging the obtained local quality matrix for the sharing nodes of the rest processors;
The calculation formula of the global node force is as follows:
，
in the formula, f a,k is the global node force of the kth time step, The node force calculated for the processor numbered I at the kth time step,Sharing node force obtained by node exchange for other processors in the kth time step;
the calculation formula of the global acceleration is as follows:
，
the calculation formula of the global speed is as follows:
。
s7, updating the material node and the shared data exchange table, including updating the material point coordinates Updating object point neighborhoodUpdating object shape functionsUpdating the derivative of the object shape functionAnd updating the shared node set by using the neighborhood updated by each boundary material point, adding a new node into the shared data exchange table C k, deleting nodes which are not in the neighborhood of the boundary material point, and resetting the shared node data exchange table C k. Then, entering the next time step, judging whether the time step is the last time step set, if not, returning to S5 to enter the next iterative calculation; if yes, outputting the calculation result of the step S6, completing the impact collision process simulation, and obtaining final kinematic information and material physical information to obtain physical information and kinetic data of the material points and nodes in the time period of t 0→tn, wherein the physical information and the kinetic data comprise deformation, stress, density, displacement, speed, acceleration and temperature, and completing the dynamic response analysis of the material.
Application example: the parallel optimal transport gridless impact collision process simulation method in the embodiment is applied to simulation of two-object impact collision process simulation. The application example takes the ultra-high-speed impact of the steel balls on the two layers of aluminum plates as an example to illustrate the parallel optimal transportation gridless method for simulating the impact collision process.
The application example simulates the condition that a steel ball with the diameter of 3 mm strikes two aluminum plates which are 100 mm apart at the speed of 6.09 km/s. The thickness of the aluminum plate is 5mm a. The case contains 1,368,048 object points, 832,682 degrees of freedom, and 64 CPUs are used for calculation on a distributed high-performance calculation server. During the impact process, the steel balls invade into the first layer of aluminum plate, as shown in (a) - (c) in fig. 5, the temperature gradually rises, flanging bulges are gradually formed on the surface of the aluminum plate along with the increase of plastic deformation, and the steel balls gradually change from spherical to oblate spherical along with the deeper invasion depth, as shown in (a) (b) (c) in fig. 6; as the penetration depth increases, the steel balls start to fail with increased stress, the back of the target cracks, forming larger petal-shaped hole edges, as shown in fig. 5 (d) - (h), and eventually the steel balls leave the surface of the first layer of aluminum plate, as shown in fig. 5 (i).
When the steel ball and the aluminum plate come into contact, the bottom of the steel ball is stressed and softened after plastic deformation, the steel ball gradually becomes elliptical, the deformation of the surface, which is in contact with the aluminum plate, is blocked along with the invasion of the aluminum plate, the middle part of the steel ball gradually concaves inwards to form a flattened sphere shape, and a certain number of failure material particles are formed, as shown by the white semitransparent material particles in fig. 6. As the depth of penetration increases until the first plate is penetrated, some of the material on the projectile adheres to the plate and some of the material forms a reverse cloud of fragments.
After the steel balls invade the first laminate, the temperature of the deformed part is increased along with the increase of plastic deformation, the softening effect of the material is obvious, and flanging bulges can be formed on the incident surface of the projectile. As the invasion increases, the steel ball is flattened, and the tensile stress is increased due to the stress wave reflected by the lower surface, and the flat plate forms spalling near the lower surface, so that the failure range of the bottom surface part is enlarged, and the penetration aperture is larger than the flanging aperture of the upper surface. Fig. 7 shows the material response of the first aluminum plate in the simulation, and fig. 8 shows the hole shape after the projectile penetrated the first aluminum plate. In fig. 8, (a) is a hole in the first laminate surface, having a diameter of about 12 mm a and a cuff height of about 4.4: 4.4 mm a; (b) a first lamina bottom surface penetration diameter of about 17.6 mm; (d) the perforation morphology in the experiment.
Fig. 9 shows the contact range of the fragment cloud generated after the steel ball hits the first target plate and the second target plate. At 31.48 mus, the range of action of the projectile and the first target plate fragment cloud on the second target plate was seen to be approximately 54 mm. Assuming that after the projectile penetrates through the first target plate, the fragment clouds do not interfere with each other and the movement track is kept to be a straight line, the action range is enlarged by 2.5 times, namely when the distance between the two target plates is 100 mm, the action range is about 135 mm; on the basis, the distance between the two target plates is adjusted to 100 mm for verification, and under the working condition, the action range of the projectile and the first target plate fragment cloud on the second target plate is approximately 136 mm.
Effect example: to simulate the copper flyer strike the rigid plane at a speed of 270 m/s. The embodiment is used for testing the acceleration performance of the parallel optimal transport gridless method for simulating the impact collision process. The present embodiment includes 1,116,024 object points and 683,772 degrees of freedom, and calculates on a distributed high-performance calculation server by using 1,2,4,8, 16, 32, and 64 CPUs, so as to obtain an effect similar to linear acceleration, and the parallel efficiency is very high, as shown in fig. 10.
The present invention has been described in detail with reference to the embodiments, but the description is only the preferred embodiments of the present invention and should not be construed as limiting the scope of the invention. All equivalent changes and modifications within the scope of the present invention should be considered as falling within the scope of the present invention.
Claims (10)
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202411442797.0A CN118965832B (en) | 2024-10-16 | 2024-10-16 | A parallel optimal transport meshless impact collision process simulation method and its application |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202411442797.0A CN118965832B (en) | 2024-10-16 | 2024-10-16 | A parallel optimal transport meshless impact collision process simulation method and its application |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN118965832A true CN118965832A (en) | 2024-11-15 |
| CN118965832B CN118965832B (en) | 2025-04-01 |
Family
ID=93407873
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202411442797.0A Active CN118965832B (en) | 2024-10-16 | 2024-10-16 | A parallel optimal transport meshless impact collision process simulation method and its application |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN118965832B (en) |
Citations (15)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20060149516A1 (en) * | 2004-12-03 | 2006-07-06 | Andrew Bond | Physics simulation apparatus and method |
| CN101657712A (en) * | 2007-04-12 | 2010-02-24 | 新日本制铁株式会社 | Breakage prediction method, calculation processing device, program, and recording medium |
| JP2014170556A (en) * | 2005-08-18 | 2014-09-18 | D.E. Shaw Research LLC | Parallel calculation architecture for calculating particle interaction |
| CN106446432A (en) * | 2016-09-30 | 2017-02-22 | 云翼超算（北京）软件科技有限公司 | Optimal transportation meshless method for solving large deformation of material |
| US20180285502A1 (en) * | 2014-05-27 | 2018-10-04 | Livermore Software Technology Corporation | Dynamically-Positioned Search Domain Used In Numerical Simulation Of An Impact Event Between Two Objects |
| CN110457785A (en) * | 2019-07-25 | 2019-11-15 | 江西理工大学 | A material information mapping method based on material point method for large deformation response of structure |
| CN112116128A (en) * | 2020-08-20 | 2020-12-22 | 河海大学 | Simulation prediction method for structural spalling and multiple spalling under impact load action |
| CN112862942A (en) * | 2021-02-08 | 2021-05-28 | 腾讯科技（深圳）有限公司 | Physical special effect simulation method and device, electronic equipment and storage medium |
| US20220067246A1 (en) * | 2020-09-01 | 2022-03-03 | Fujitsu Limited | Simulation method and system |
| EP4009216A1 (en) * | 2020-12-01 | 2022-06-08 | Fraunhofer-Gesellschaft zur Förderung der angewandten Forschung e.V. | Method for modelling a collision of two vehicles |
| CN114912263A (en) * | 2022-04-30 | 2022-08-16 | 云翼超算(苏州)软件科技有限公司 | Parallel optimal transport grid-free processing method and system based on graph topological structure |
| CN114970226A (en) * | 2021-02-26 | 2022-08-30 | 湖南大学 | A simulation calculation method for material fracture damage based on peridynamics |
| CN115803744A (en) * | 2020-05-14 | 2023-03-14 | 牛津大学科技创新有限公司 | Computational analysis of physical systems |
| CN117034689A (en) * | 2023-08-02 | 2023-11-10 | 大连理工大学 | Soil liquefaction large deformation analysis method based on gridless RBF mapping technology |
| CN118673731A (en) * | 2024-08-20 | 2024-09-20 | 云翼(嘉兴)软件科技有限公司 | Impact large deformation process simulation method based on gridless interpolation technology and application |
-
2024
- 2024-10-16 CN CN202411442797.0A patent/CN118965832B/en active Active
Patent Citations (15)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20060149516A1 (en) * | 2004-12-03 | 2006-07-06 | Andrew Bond | Physics simulation apparatus and method |
| JP2014170556A (en) * | 2005-08-18 | 2014-09-18 | D.E. Shaw Research LLC | Parallel calculation architecture for calculating particle interaction |
| CN101657712A (en) * | 2007-04-12 | 2010-02-24 | 新日本制铁株式会社 | Breakage prediction method, calculation processing device, program, and recording medium |
| US20180285502A1 (en) * | 2014-05-27 | 2018-10-04 | Livermore Software Technology Corporation | Dynamically-Positioned Search Domain Used In Numerical Simulation Of An Impact Event Between Two Objects |
| CN106446432A (en) * | 2016-09-30 | 2017-02-22 | 云翼超算（北京）软件科技有限公司 | Optimal transportation meshless method for solving large deformation of material |
| CN110457785A (en) * | 2019-07-25 | 2019-11-15 | 江西理工大学 | A material information mapping method based on material point method for large deformation response of structure |
| CN115803744A (en) * | 2020-05-14 | 2023-03-14 | 牛津大学科技创新有限公司 | Computational analysis of physical systems |
| CN112116128A (en) * | 2020-08-20 | 2020-12-22 | 河海大学 | Simulation prediction method for structural spalling and multiple spalling under impact load action |
| US20220067246A1 (en) * | 2020-09-01 | 2022-03-03 | Fujitsu Limited | Simulation method and system |
| EP4009216A1 (en) * | 2020-12-01 | 2022-06-08 | Fraunhofer-Gesellschaft zur Förderung der angewandten Forschung e.V. | Method for modelling a collision of two vehicles |
| CN112862942A (en) * | 2021-02-08 | 2021-05-28 | 腾讯科技（深圳）有限公司 | Physical special effect simulation method and device, electronic equipment and storage medium |
| CN114970226A (en) * | 2021-02-26 | 2022-08-30 | 湖南大学 | A simulation calculation method for material fracture damage based on peridynamics |
| CN114912263A (en) * | 2022-04-30 | 2022-08-16 | 云翼超算(苏州)软件科技有限公司 | Parallel optimal transport grid-free processing method and system based on graph topological structure |
| CN117034689A (en) * | 2023-08-02 | 2023-11-10 | 大连理工大学 | Soil liquefaction large deformation analysis method based on gridless RBF mapping technology |
| CN118673731A (en) * | 2024-08-20 | 2024-09-20 | 云翼(嘉兴)软件科技有限公司 | Impact large deformation process simulation method based on gridless interpolation technology and application |
Non-Patent Citations (2)
| Title |
|---|
| BAZEIA, D.等: "Scattering of compactlike structures", EUROPEAN PHYSICAL JOURNAL C, vol. 79, no. 12, 11 December 2019 (2019-12-11), pages 1 - 11 * |
| 王长生等: "基于MATLAB矢量化物质点法的车身防撞梁碰撞分析", 计算机辅助工程, vol. 33, no. 01, 31 March 2024 (2024-03-31), pages 28 - 36 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN118965832B (en) | 2025-04-01 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Badcock et al. | Elements of computational fluid dynamics on block structured grids using implicit solvers | |
| CN105843073B (en) | A kind of wing structure aeroelastic stability analysis method not knowing depression of order based on aerodynamic force | |
| CN105302974B (en) | A kind of real-time cutting simulation method of flexible article analyzed based on finite element and time-varying modal | |
| CN110457790A (en) | Peridynamic Discontinuous Galerkin Finite Element Method for Structural Deformation Analysis | |
| Sawant et al. | Application of adaptively refined unstructured grids in DSMC to shock wave simulations | |
| Gaitonde et al. | A three-dimensional moving mesh method for the calculation of unsteady transonic flows | |
| CN117725859A (en) | Fatigue crack growth prediction method and system for aircraft fatigue damage key part | |
| CN118965832B (en) | A parallel optimal transport meshless impact collision process simulation method and its application | |
| CN116306246A (en) | Large wallboard riveting deformation prediction and optimization method | |
| CN118673731A (en) | Impact large deformation process simulation method based on gridless interpolation technology and application | |
| CN103310049A (en) | Optimization method of two-dimension solid-phase phonon crystal mixed mode band gap | |
| CN111475940B (en) | A Flexible Baseline Dynamic Prediction Method Based on Fiber Bragg Grating Sensors and Wing Modes | |
| CN117688807B (en) | Structural deformation and damage analysis method based on indirect force bond finite element method | |
| CN118917192A (en) | Flexible spacecraft layout optimization method based on multi-target mixed particle swarm | |
| CN109598059B (en) | Thermal protection system optimization design method and design system based on proxy model | |
| CN114912263A (en) | Parallel optimal transport grid-free processing method and system based on graph topological structure | |
| Banerjee et al. | Vibration analysis and prediction natural frequencies of cracked timoshenko beam by two optimization techniques-Cascade ANN and ANFIS | |
| Zhang et al. | Reliability-based design optimization for flexible mechanism with particle swarm optimization and advanced extremum response surface method | |
| CN114186442B (en) | Lattice material parameter optimization method based on neural network model and numerical simulation | |
| Li et al. | A parallel material-point method with application to 3D solid mechanics | |
| CN118052167B (en) | Method for constructing flow field model with multidimensional correlation response | |
| CN115828711B (en) | Method and system for predicting ice breaking residual speed of underwater structure | |
| CN118734555B (en) | Improved triangle mesh numerical manifold method for simulating progressive damage of jointed rock slope | |
| CN113505430B (en) | Large-scale space flexible film dynamic parameter calculation method | |
| CN118133641B (en) | Deformation and damage analysis method of solid structure based on quasi-bond finite element method |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant |