SPH particle pair rapid search method suitable for explosion scene

Abstract

The invention discloses an SPH particle pair rapid search method suitable for an explosion scene, and relates to the field of smoothed particle fluid dynamics numerical simulation. The particle information of the fluid particles within the potential interaction range around the fluid particles is stored to construct the neighbor list, so that the stored neighbor list can be reused to search interaction particle pairs, global particle pair search does not need to be performed every time, the number of times of performing global search can be greatly reduced, and the search efficiency is improved. The time consumed by the particles for searching is reduced, so that the calculation efficiency of the SPH method is improved, and the method is particularly suitable for scenes with substance splashing such as explosion and penetration.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of smoothed particle hydrodynamics numerical simulation, and in particular to a method for quickly searching for particle pairs in an explosion scenario. BACKGROUND

[0002] The smoothed particle hydrodynamics (SPH) method is currently the most in-depth and mature meshless method, and has been widely used in many fields. Due to the flexibility of meshless, the SPH method performs excellently in simulations of scenarios such as multi-phase interface capture and solid extreme deformation.

[0003] However, due to the lack of fixed grid connection, the SPH method needs to search for interacting particle pairs in each time step iteration, which increases the time cost and makes the SPH method less efficient than the finite element method. This high time cost greatly restricts the promotion and use of the SPH method. Therefore, finding a fast particle pair search method has become an important issue in the SPH field. The simplest and most intuitive particle pair search method is to traverse all particles, which is called direct search method. However, its quadratic time complexity makes it difficult to use in large-scale practical calculations. The more commonly used particle pair search methods are linked list search method and tree search method.

[0004] The linked list search method generates a background grid at each time step, and the grid length is the particle interaction radius. In this way, the particles in the current grid will only interact with the particles in the adjacent grid, which will greatly reduce the number of judgments in the particle pair search process and save calculation time. However, this method has some significant problems when simulating explosion scenarios. First, explosions are usually accompanied by the spatter of matter, which means that the volume of space occupied by particles will inevitably increase as the simulation process progresses. Therefore, the number of background grids will inevitably increase, resulting in an increase in memory usage and a decrease in computational efficiency. Another problem is that for explosion scenarios, the rapid expansion of the explosion gas often causes particles to be unable to obtain enough neighboring particles. Therefore, adaptive smoothing length and adaptive particle splitting techniques are often used to ensure simulation accuracy. When the particle smoothing lengths are not consistent, the background grid must take the largest smoothing length of all particles to ensure the correctness of the results. Therefore, the number of background grids decreases, which ultimately leads to a decrease in algorithm efficiency.

[0005] In contrast, the tree search method can partially solve the above problems. At each time step, the tree search method constructs a tree structure based on the particle positions and then searches for particle pairs based on this tree structure. Since this method does not depend on the background mesh, it avoids the problems associated with the linked list search method when simulating explosion problems. However, this method is slower and not well-suited for large-scale computations. Summary of the Invention

[0006] To address the aforementioned problems and technical requirements, this application proposes a fast SPH particle pair search method suitable for explosion scenarios. The technical solution of this application is as follows:

[0007] A fast SPH particle pair search method suitable for explosion scenarios is proposed. This method involves dividing the computational domain of the explosion scenario into multiple sub-regions, creating a corresponding process for each sub-region, and having multiple processes execute in parallel. Each process executes the following methods:

[0008] Obtain the physical information of each fluid particle in the corresponding sub-region and its outer region at the current time step. The physical information of the fluid particles includes three-dimensional coordinates, fluid physical quantities, and particle radius of action.

[0009] Based on the physical information of each fluid particle in the corresponding sub-region and its outer region at the first time step, a particle space layout octree structure is constructed. Based on the particle space layout octree structure, a nearest neighbor list of each fluid particle in the corresponding sub-region is constructed. The nearest neighbor list of each fluid particle in the corresponding sub-region is used as the search range to search for SPH particle pairs that interact within the current time step, and the nearest neighbor list of each fluid particle in the corresponding sub-region is stored. The nearest neighbor list of each fluid particle includes the global ID of all fluid particles within the potential interaction range of the fluid particle.

[0010] Starting from the second time step, the physical information of each fluid particle in the corresponding sub-region at the current time step is used to check whether the nearest neighbor list update condition is met. If it is determined that the nearest neighbor list update condition is not met, the nearest neighbor list that has been determined and stored in the historical time steps of each fluid particle in the corresponding sub-region is used as the search range to search for interacting SPH particle pairs in the current time step. If it is determined that the nearest neighbor list update condition is met, the particle space layout octree structure is reconstructed based on the physical information of each fluid particle in the corresponding sub-region and its extended region at the current time step, and the nearest neighbor list of each fluid particle in the corresponding sub-region is updated. The updated nearest neighbor list of each fluid particle in the corresponding sub-region is used as the search range to search for interacting SPH particle pairs in the current time step.

[0011] The beneficial technical effects of this application are:

[0012] The application discloses a SPH particle pair fast searching method suitable for explosion scenes.

[0013] The method uses a tree searching method as a preprocessing to construct the neighbor list, and updates the neighbor list in time when the neighbor list updating condition is met, so as to guarantee the effectiveness and reliability of searching, can effectively reduce the searching time consumption, and has better stability compared with a linked list searching method. The method is easy to parallelize, can be conveniently transplanted when a parallelization based on region decomposition is used, and the application of the hash array in the searching process can assist in improving the searching efficiency. The compressed transmission of the neighbor list can save the time consumption of inter-process information exchange, not only improves the information exchange efficiency, but also reduces the communication cost. BRIEF DESCRIPTION OF DRAWINGS

[0014] Figure 1 FIG. 1 is a method flowchart of the SPH particle pair fast searching method in an embodiment of the application.

[0015] Figure 2 FIG. 3 is a two-dimensional projection schematic diagram of four sub-regions obtained by dividing the calculation domain in an example.

[0016] Figure 3 FIG. 5 is a schematic diagram of constructing the octree structure of the particle space layout in an example.

[0017] Figure 4 FIG. 7 is a method flowchart of constructing the neighbor list of the fluid particle in an embodiment of the application.

[0018] Figure 5 FIG. 9 is a two-dimensional projection schematic diagram of the search domain of a single fluid particle when the linked list searching method and the searching method of the application are used in an example.

[0019] Figure 6 FIG. 11 is a calculation configuration schematic diagram in an underwater explosion example.

[0020] Figure 7 Figure 6 ​The comparison chart of the time consumption of the particle pair search step in the underwater explosion example of the traditional tree search method, the linked list search method and the search method of the present application with the change of the time step.

[0021] Figure 8 is Figure 6 The comparison chart of the time consumption of the particle pair search step in the underwater explosion example of the traditional tree search method, the linked list search method and the search method of the present application with the change of the time step. DETAILED DESCRIPTION

[0022] The specific embodiments of the present application will be further described in combination with the accompanying drawings.

[0023] Considering that the motion of the fluid particles in the SPH method in space is continuous, the core idea of the SPH particle pair fast search method suitable for explosion scenarios disclosed in the present application is to save the information of other fluid particles around each fluid particle, and the subsequent can directly reuse these already saved historical information, thereby being conducive to improving the search efficiency to adapt to large-scale calculation.

[0024] In order to reuse the information of other fluid particles around each fluid particle, a neighbor list is saved for each fluid particle, and the global id of all fluid particles located in the potential interaction range of the fluid particle is included in the neighbor list of each fluid particle, and the global id of each fluid particle is used to uniquely identify a fluid particle. Considering that the fluid particles will move in space, therefore the potential interaction range of each fluid particle is determined on the basis of the particle interaction radius of the fluid particle and also considers the motion process of the fluid particle, and the potential interaction range is slightly larger than the particle interaction radius of the fluid particle.

[0025] In this way, if no fluid particle outside the neighbor list enters the interaction range of a fluid particle within a certain time, the interacting particle pair can be directly found in the neighbor list without the need for global search, so the time consumption of the particle pair search step can be greatly reduced, thereby improving the search efficiency. As can be seen from the above introduction, the prerequisite for directly reusing the neighbor list of the fluid particle for particle pair search is that the fluid particles that may interact are all in the neighbor list, but as the fluid particles continuously move in space, other fluid particles outside the neighbor list will inevitably enter the interaction range of the fluid particles, at which time the information recorded in the already stored neighbor list is no longer accurate and becomes invalid, and the invalid neighbor list can no longer be used for particle pair search, so it needs to be dynamically updated.

[0026] Since the SPH method itself is slow in calculation, parallelization is a natural choice when facing large-scale problems, so the whole application adopts a parallelization method based on domain decomposition, which can maximize the preservation of the original algorithm format and is widely used. In the SPH domain decomposition method, the calculation domain in the explosion scene is divided into multiple sub-regions, so that the fluid particles in the calculation domain are attributed to different sub-regions, and a corresponding process is created for each sub-region. When dividing the calculation domain, the calculation domain is divided into multiple sub-regions based on the number of fluid particles, so that the number of fluid particles contained in each sub-region is as similar as possible, thereby optimizing the parallel effect. In the explosion scene, as the time step advances, the volume occupied by the fluid particles will differ significantly, so the sub-region division can be re-performed according to the actual distribution of the fluid particles at different time steps.

[0027] Multiple processes perform the SPH particle pair fast search method in parallel based on the above design idea, and the method performed by each process includes the following steps, please refer to the flowchart shown in Figure 1 .

[0028] In step 110, the physical information of each fluid particle in the sub-region and its extended region corresponding to the process at the current time step is obtained, and the physical information of each fluid particle includes the three-dimensional coordinates of the fluid particle, the fluid physical quantity and the particle action radius. The fluid physical quantity of the fluid particle commonly includes particle velocity, acceleration, density, specific internal energy, smoothing length, etc.

[0029] The extended region of each sub-region includes the region between the boundary of the sub-region and the other adjacent sub-regions, and the distance between the boundary of the current sub-region and the other adjacent sub-regions is less than or equal to r max , r max is the maximum value of the particle action radius of all fluid particles in the entire calculation domain at the current time step.

[0030] Each process maintains the physical information of each fluid particle in its corresponding sub-region, but the process also needs to use the physical information of the fluid particles in the extended region of its corresponding sub-region when performing the SPH particle pair fast search method, and the physical information of the fluid particles in the extended region is maintained by other processes, so there is an information exchange process between different processes. Each process obtains the physical information of the fluid particles in the extended region of its corresponding sub-region from other processes through the information exchange process. In actual implementation, each process will sequentially traverse the physical information of each fluid particle in its corresponding sub-region, and screen out each fluid particle whose distance between the boundary of the adjacent sub-region is less than or equal to r max . The screened fluid particles are the fluid particles in the extended region of the adjacent sub-region, and then the physical information of the screened fluid particles is sent to the process corresponding to the adjacent sub-region to realize information exchange.

[0031] The calculation domain in the explosion scenario to which the present application is directed is a three-dimensional calculation domain, and thus each sub-region obtained by division is also a three-dimensional region. In order to facilitate illustration, a two-dimensional projection of fluid particles contained in each sub-region obtained by division of the calculation domain in one instance at one time step is shown in FIG. 1. In this instance, the calculation domain is divided into four sub-regions, which are respectively denoted as sub-region 1 to sub-region 4. A corresponding process is constructed for each of the four sub-regions, which are respectively denoted as process 1 to process 4. Taking process 2 as an example, the sub-region corresponding to process 2 and its expanded region are shown by the shaded area in FIG. 2. Figure 2 Figure 2

[0032] Process 1 sequentially traverses the physical information of each fluid particle in the sub-region 1 corresponding to process 1, and screens out fluid particle p7 whose distance from the boundary with the adjacent sub-region 2 is less than or equal to r max , and sends the physical information of fluid particle p7 to process 2 corresponding to the adjacent sub-region 2. Similarly, process 4 screens out fluid particles p8 and p9 in the sub-region 4 corresponding to process 4 whose distance from the boundary with the adjacent sub-region 2 is less than or equal to r max , and sends the physical information of fluid particles p8 and p9 to process 2 corresponding to the adjacent sub-region 2. Sub-region 2 also sends the physical information of fluid particles in the expanded region to process 1 and process 4, which will not be described here.

[0033] Thus, process 2 maintains the physical information of fluid particles p1, p2, p3, p4, p5, and p6 in the sub-region 2 corresponding to process 2, and obtains the physical information of fluid particle p7 from process 1 through the information exchange process, and obtains the physical information of fluid particles p8 and p9 from process 4 through the information exchange process, and finally obtains the physical information of all nine fluid particles p1 to p9 in the sub-region corresponding to process 2 and its expanded region at the current time step.

[0034] At step 120, in the first time step, the near-neighbor list of each fluid particle has not been stored, and thus the particle space layout octree structure is constructed based on the physical information of each fluid particle in the sub-region corresponding to the process and its expanded region in the first time step, and the near-neighbor list of each fluid particle in the sub-region corresponding to the process is constructed based on the particle space layout octree structure. After the near-neighbor list of each fluid particle in the sub-region corresponding to the process is constructed, the search range is set as the near-neighbor list of each fluid particle in the corresponding sub-region to search for the SPH particle pairs that exist in interaction at the current time step.

[0035] In addition, the process also stores the near-neighbor list of each fluid particle in the corresponding sub-region, and also records the three-dimensional coordinates X old ​​That is, in the present application, each process maintains particle information of fluid particles in the corresponding sub-region, the particle information including physical information of the fluid particles, a neighbor list, and three-dimensional coordinates X old Since the length of the neighbor list of each fluid particle is difficult to determine, a variable-length storage container such as a Vector is used to store the particle information, and the particle information of fluid particles in the corresponding sub-region and its extended region at the current time step and the respective local id are stored in the storage container, so as to minimize the memory usage.

[0036] For the explosion scene calculation domain of a three-dimensional space, the search tree is an octree, and the particle space layout octree structure includes a plurality of tree nodes, each tree node corresponding to a three-dimensional node control region, the node control region having a center coordinate and a control region radius. Each tree node corresponds to a unique fluid particle, and the three-dimensional coordinates of the fluid particle corresponding to each tree node are located in the node control region of the tree node. When the tree node has child nodes, the three-dimensional coordinates of the fluid particles corresponding to all child nodes of the tree node are also located in the node control region of the tree node.

[0037] In an example, the sub-region and its extended region corresponding to the process form a cubic region, and the distribution of four fluid particles m1-m4 in the sub-region and its extended region is shown in FIG. 1. Figure 3 The method for constructing the particle space layout octree structure includes the following steps:

[0038] First, a tree node n0 is initialized as a root node, the node control region corresponding to the tree node n0 being the cubic region formed by the sub-region and its extended region corresponding to the process, the center coordinate of the node control region corresponding to the tree node n0 being the center coordinate of the cubic region, and the control region radius of the node control region being half of the length of the cubic region. A fluid particle m1 is taken, and it is determined that the three-dimensional coordinates of the fluid particle m1 are located in the node control region corresponding to the tree node n0, and the tree node n0 does not have a corresponding fluid particle at this time. Therefore, the fluid particle m1 is added to the current tree structure as the fluid particle corresponding to the tree node n0.

[0039] Continue to take a fluid particle m2, from the root node of the current tree structure, in turn, first determine the three-dimensional coordinates of the fluid particle m2 in the node control area corresponding to the tree node n0, but at this time the tree node n0 has a corresponding fluid particle m1, so the node control area of the tree node n0 is divided into 8 quadrants, the 8 quadrants are the node control areas corresponding to the 8 child nodes of the tree node n0, and the control area radius of the node control areas corresponding to the 8 child nodes is half of the edge length of the node control area of the tree node n0. Assuming that the three-dimensional coordinates of the fluid particle m2 are determined to be located in the node control area corresponding to the tree node n11, the fluid particle m2 is added to the current tree structure as the fluid particle corresponding to the tree node n11.

[0040] Continue to take a fluid particle m3, from the root node of the current tree structure, in turn, first determine the three-dimensional coordinates of the fluid particle m3 in the node control area corresponding to the tree node n0, but at this time the tree node n0 has a corresponding fluid particle m1. Continue to traverse the node control areas corresponding to the 8 child nodes of the tree node n0. Assuming that the three-dimensional coordinates of the fluid particle m3 are determined to be located in the node control area corresponding to the tree node n12, and at this time the tree node n12 does not have a corresponding fluid particle, the fluid particle m3 is added to the current tree structure as the fluid particle corresponding to the tree node n12.

[0041] Continue to take a fluid particle m4, and traverse according to the above process, and determine that the three-dimensional coordinates of the fluid particle m4 are located in the node control area corresponding to the tree node n12, but at this time the tree node n12 has a corresponding fluid particle m3, so the node control area of the tree node n12 is further divided into 8 quadrants, the 8 quadrants are the node control areas corresponding to the 8 child nodes of the tree node n12, and the control area radius of the node control areas corresponding to the 8 child nodes is half of the edge length of the node control area of the tree node n12. Assuming that the three-dimensional coordinates of the fluid particle m4 are determined to be located in the node control area corresponding to the tree node n21, the fluid particle m4 is added to the current tree structure as the fluid particle corresponding to the tree node n21.

[0042] In the actual application process, in the case of larger data size, it can continue to be constructed in the same way. The particle space layout octree structure constructed in this way contains all the fluid particles in the corresponding sub-region and the expanded region of the process, each tree node corresponds to a node control area, and each tree node has and only has one corresponding fluid particle, which means that the fluid particle and all the fluid particles of its child nodes are placed in the current node control area. For example, in the example of Figure 3 , the fluid particle m3 corresponding to the tree node n12 and the fluid particle m4 corresponding to the child node n21 of the tree node n12 are placed in the node control area of the tree node n12.

[0043] After the particle space layout octree structure is constructed, the three-dimensional coordinates of each fluid particle in the process corresponding sub-region and its outer expansion region become easy to grasp, and the neighbor list of each fluid particle in the process corresponding sub-region can be constructed based on the particle space layout octree structure.

[0044] Since the fluid particles are constantly moving, the motion process of the fluid particles also needs to be considered to construct the neighbor list of each fluid particle. In an embodiment, according to the maximum value v of the particle velocity of all fluid particles in the entire calculation domain at the current time step max The buffer distance b is calculated according to the following formula:

[0045] b = 2av max dt

[0046] Where a is an input parameter used to determine the update frequency of the neighbor list. dt is the time step.

[0047] The buffer distance b calculated in this way represents the maximum value of the relative distance reduction of two fluid particles within a time step, and therefore, based on the particle space layout octree structure combined with the buffer distance b, the neighbor list of each fluid particle in the process corresponding sub-region can be constructed, which can determine all fluid particles in the potential interaction range of each fluid particle considering the motion process of the fluid particles.

[0048] The construction of the neighbor list of any fluid particle i includes the following steps, please refer to the flowchart shown in Figure 4 .

[0049] Step 121, starting from the root node, sequentially traverse the tree nodes in the particle space layout octree structure according to the predetermined traversal order. The predetermined traversal order can use the existing tree structure traversal method, and generally uses the depth-first traversal method.

[0050] Step 122, for any tree node node currently traversed, according to the physical information of the fluid particle i combined with the buffer distance b, detect whether the scope of the fluid particle i intersects with the node control region of the tree node node.

[0051] Specifically, when the following conditions are met, it is determined that the scope of the fluid particle i intersects with the node control region of the tree node node, otherwise it is determined that the scope of the fluid particle i does not intersect with the node control region of the tree node node.

[0052]

[0053] Where (x i , y i , z i) is the three-dimensional coordinate of fluid particle i at the current time step, (x node ,y node ,z node ) is the center coordinate of the node control region corresponding to the tree node node, d node is the region radius of the node control region corresponding to the tree node node, r i is the particle action radius of fluid particle i.

[0054] Step 123, when it is determined that the action scope of fluid particle i intersects with the node control region of tree node node, it is continued to detect whether the fluid particle j corresponding to the tree node node is in the potential interaction range of fluid particle i. When it is determined that the fluid particle j corresponding to the tree node node is located in the potential interaction range of fluid particle i, the global id of the fluid particle j is added to the neighbor list of fluid particle i, and it is continued to traverse the tree nodes in the particle space layout octree structure in the predetermined traversal order one by one. Otherwise, when it is determined that the fluid particle j corresponding to the tree node node is not in the potential interaction range of fluid particle i, it is directly continued to traverse the tree nodes in the particle space layout octree structure in the predetermined traversal order one by one.

[0055] For any fluid particle j corresponding to the tree node node, when the distance D now between the three-dimensional coordinate X now of fluid particle j at the current time step and the three-dimensional coordinate X ij of fluid particle i at the current time step satisfies the following formula, it is determined that the fluid particle j is in the potential interaction range of fluid particle i, otherwise it is determined that the fluid particle j is not in the potential interaction range of fluid particle i:

[0056] D ij < r ij + b

[0057] In the above formula, r ij is the interaction radius between fluid particle i and fluid particle j. In the explosion scene to which the present application is directed, considering that adaptive time step or adaptive particle splitting method will be used in the simulation, it is expected that each fluid particle has its own particle action radius in the calculation process, so in order to ensure the symmetry of the interaction between fluid particles, the interaction radius r ij between fluid particle i and fluid particle j is calculated according to the following formula:

[0058] r ij = (r i + r j ) / 2

[0059] Wherein, r i is the particle action radius of fluid particle i, rj is the particle interaction radius of fluid particle j. This formula means that there are different interaction radii r ij between fluid particle i and different fluid particle j

[0060] As can be seen from the formula, the greater the value of parameter a, the lower the update frequency of the neighbor list, but the greater the value of the buffer distance b, which will result in a larger size of the neighbor list of each fluid particle, resulting in greater memory usage. Therefore, parameter a needs to be adjusted by the user according to the actual problem and demand.

[0061] Since the final demand result is the interacting particle pair, if both fluid particles save the other in the neighbor list, it will cause redundancy, so in another embodiment, when it is determined that the scope of fluid particle i intersects with the node control area of tree node node, first compare the particle interaction radius r j of fluid particle j corresponding to tree node node with the particle interaction radius r i of fluid particle i:

[0062] When r i > r j , perform the step of detecting whether the fluid particle j corresponding to the tree node node is in the potential interaction range of the fluid particle i, and when it is determined that the fluid particle j is in the potential interaction range of the fluid particle i, add the global id of the fluid particle j to the neighbor list of the fluid particle i.

[0063] When r i ≤ r j , even if the fluid particle j and the fluid particle i are in the potential interaction range of each other, the global id of the fluid particle j is not recorded in the neighbor list of the fluid particle i, at this time, the fluid particle j is no longer detected whether it is in the potential interaction range of the fluid particle i, but directly continues to traverse the tree nodes in the particle space layout octree structure in turn according to the predetermined traversal order. In this case, in the process of constructing the neighbor list of the fluid particle j, the global id of the fluid particle i will be recorded in the neighbor list of the fluid particle j. In this way, only the global id of the other fluid particle is saved in the neighbor list of the fluid particle with a larger particle interaction radius, which can effectively reduce the memory usage.

[0064] Step 124, when it is determined that the scope of fluid particle i does not intersect with the node control area of tree node node, it can be directly determined that the fluid particles corresponding to the tree node node and all its child nodes are not possible in the potential interaction range of the fluid particle i, without the need to continue to compare. At this time, the subsequent traversal of the tree node node and all its child nodes is skipped, and the tree nodes in the particle space layout octree structure are traversed in turn according to the predetermined traversal order.

[0065] For example, in the example of Figure 3 In the example of

[0066] At step 130, the process has stored the neighbor list of each fluid particle in the corresponding sub-region since the second time step.

[0067] As the fluid particles move, the fluid particles contained in each sub-region will change at different time steps, and fluid particles can move from other sub-regions to the sub-region corresponding to the process, and similarly, fluid particles in the sub-region corresponding to the process can also move to other sub-regions. Therefore, the fluid particles in the sub-region corresponding to the process at the current time step can have been in the sub-region corresponding to the process at the previous time step, or can have moved to the sub-region corresponding to the process from other sub-regions at the current time step. For those fluid particles that have been in the sub-region corresponding to the process at the previous time step, the process itself maintains the particle information of these fluid particles. For those fluid particles that have just moved to the sub-region corresponding to the process from other sub-regions, the process needs to obtain the particle information of these fluid particles from other processes through an information exchange process. In the present application, since the neighbor list of the fluid particle needs to be continuously saved at different time steps, the neighbor list also needs to be exchanged between different processes along with the fluid particles, that is, when a fluid particle m moves from the sub-region Q corresponding to process Q to the sub-region P corresponding to process P, the process Q needs to send the particle information of the fluid particle m maintained by the process Q to the process P, and the particle information of the fluid particle m also includes the neighbor list of the fluid particle m that has been determined and stored at the historical time steps and the three-dimensional coordinates of the fluid particle m when the neighbor list is constructed. For example, in the example of Figure 2In the example of Fig. 1, when fluid particle p10 moves from sub-region 4 to sub-region 2, process 4 needs to send the particle information of fluid particle p10 to process 2.

[0068] When process Q sends the neighbor list containing the particle information of fluid particle m to process P, because the number of fluid particles contained in the neighbor list of different fluid particles is inconsistent, in the transmission of the neighbor list of multiple fluid particles, the traditional method often establishes a two-dimensional array based on the length of the neighbor list containing the largest number of fluid particles, and then transmits the two-dimensional array to process P. However, this will result in a considerable amount of data exchange, and the two-dimensional array will carry a lot of meaningless information, affecting the exchange efficiency. Therefore, in order to reduce communication consumption, compression operation is considered. For multiple fluid particles m moving into sub-region P, process Q will first concatenate the neighbor list of each fluid particle m determined and stored in the history time step to convert it into a one-dimensional list array. In addition, a number array is generated to record the number of global ids contained in the neighbor list of each fluid particle m screened out. Then process Q sends the one-dimensional list array and the number array to process P, and process P decompresses the list array according to the number of global ids contained in the neighbor list of each fluid particle m in the number array to obtain the neighbor list of each fluid particle m moving into sub-region P from sub-region Q, which has been determined and stored in the history time step. The list array obtained by conversion reduces the transmission of redundant information compared with the traditional two-dimensional array, thereby saving the time-consuming and communication consumption of information exchange between different processes.

[0069] In the process of information exchange between processes, each process has obtained the particle information of each fluid particle in the corresponding sub-region, including the neighbor list. As analyzed above, the precondition for being able to search for the particle pair that interacts with fluid particle i using the neighbor list of any fluid particle i is that the fluid particles that may interact with fluid particle i are all in the neighbor list of fluid particle i. However, if the movement distance of fluid particle i is greater than half of the buffer distance b, and there is a fluid particle outside the neighbor list that is moving opposite to fluid particle i, then this fluid particle outside the neighbor list may also interact with fluid particle i in the current time step. In order to avoid this situation, first, it is necessary to detect whether the neighbor list update condition is met according to the physical information of each fluid particle in the corresponding sub-region in the current time step. Similarly, detection needs to be made based on the movement process of the fluid particle, so the buffer distance b needs to be combined to detect whether the neighbor list update condition is met. Specifically:

[0070] Each process detects the three-dimensional coordinates X now (i) the three-dimensional coordinates Xold (i) between the moving distance AX(i) of the fluid particle i and the moving distance AX(i) of the fluid particle i+1. When the process determines that there exists a moving distance AX(i) in the corresponding sub-region of the fluid particle i, it determines that there exists a target fluid particle in the corresponding sub-region of the process. And when all fluid particles in the corresponding sub-region of the process satisfy , the above risk does not occur, and it is determined that there is no target fluid particle in the corresponding sub-region of the process.

[0071] Each process respectively detects whether there is a target fluid particle in the corresponding sub-region. When there is a sub-region containing a target fluid particle in the entire calculation domain, the above risk of interaction with the relative motion of the fluid particles outside the neighbor list may occur, and the stored neighbor list cannot be directly used for particle pair search. At this time, it is determined that the neighbor list update condition is met. When there is no sub-region containing a target fluid particle in the entire calculation domain, the stored neighbor list can be directly used for particle pair search, and it is determined that the neighbor list update condition is not met.

[0072] Step 140, when it is determined that the neighbor list update condition is not met, the neighbor list of each fluid particle in the corresponding sub-region that has been determined and stored in the historical time step is directly used as the search range to search for the SPH particle pair that exists in the current time step.

[0073] Step 150, when it is determined that the neighbor list update condition is met, each process reconstructs the particle space layout octree structure based on the physical information of each fluid particle in the corresponding sub-region and its extended region in the current time step, and updates the neighbor list of each fluid particle in the corresponding sub-region. The specific method is similar to step 120, which will not be repeated here. At the same time, the three-dimensional coordinates X old (i) of each fluid particle in the current time step in which the neighbor list is updated will also be updated. Then the updated neighbor list of each fluid particle in the corresponding sub-region is used as the search range to search for the SPH particle pair that exists in the current time step.

[0074] ​The method of searching for SPH particle pairs that interact in the current time step in steps 120, 140, and 150 above is the same, which includes: for any fluid particle p in the corresponding sub-region of the process, sequentially traversing each global id in the neighbor list of the fluid particle p, and detecting whether the fluid particle p and the fluid particle q corresponding to the global id form an SPH particle pair according to the physical information of the fluid particle q and the physical information of the fluid particle p. In this process, the physical information of the fluid particle p and the fluid particle q needs to be obtained, and as described above, the process stores the particle information and the local id of each fluid particle in the corresponding sub-region and the extended region thereof in the storage container, and the global id of the fluid particle is recorded in the neighbor list. In order to correspond the local id of the fluid particle in the storage container with the global id, the process constructs a hash table to record the mapping relationship between the global id of each fluid particle in the storage container and the local id thereof in the storage container. Then the process sequentially traverses each global id in the neighbor list of each fluid particle p in the storage container P, so that when any global id in the neighbor list of the fluid particle p is traversed, the local id corresponding to the global id in the hash table is directly detected, and the physical information of the corresponding fluid particle q in the storage container can be read according to the local id. When it is determined that the traversed global id is not included in the hash table, it can be determined that the fluid particle corresponding to the global id will not form an SPH particle pair with the fluid particle p, and then the traversal of other global ids in the neighbor list of the fluid particle p is continued. In this way, it is not necessary to perform a loop search on the storage container for each global id in each neighbor list, which is beneficial to improve the efficiency.

[0075] The improvement of the calculation efficiency of the method compared with the traditional tree search method is obvious, because the number of search times of the method is greatly reduced compared with the traditional tree search method. Therefore, only one example is used to illustrate the speed difference between the method of the present application and the most commonly used linked list method in the field of SPH at present. In the example, it is assumed that the fluid particles are uniformly distributed, and the size of the search domain searched by a single fluid particle can reflect the speed of the search. For the red particle in Figure 5 , when the linked list search method is used, the side length of the background grid is the particle interaction radius r of the red particle, and then all fluid particles in the grid where the red particle is located and the surrounding adjacent grids need to be searched when the linked list search method is used. The search area in the two-dimensional plane projection is a rectangular area as shown by the light blue square in Figure 5 . For three-dimensional, the search domain volume is:

[0076] V1=27r 3

[0077] However, when using the method of this application, it is necessary to search for all fluid particles within a spherical range of r+b around the red particle. The search area for the two-dimensional plane projection is as follows: Figure 5 The circular region r+b in the middle has a search domain volume of: (in 3D)

[0078]

[0079] The buffer distance b is typically taken as around 0.3r, then V2 = 9.2r. 3 .

[0080] In comparison, the search domain volume V1 = 27r for the linked list search method in this example. 3 And the search domain volume V2 = 9.2r of the method in this application. 3 It can be seen that the search domain volume using the method of this application is only one-third of the search domain volume of the linked list search method. Therefore, the search speed of the method of this application is higher than that of the linked list search method.

[0081] Furthermore, the effectiveness of the proposed method is verified through an underwater explosion case test. The calculation uses the classic SPH method, and its update format is as follows:

[0082]

[0083] In the formula, ρ is the particle density, m is the particle mass, v is the particle velocity, e is the particle specific energy, and v ij Let be the velocity difference between fluid particles i and j, t be time, p be particle pressure, Π be the artificial viscosity term, N be the number of particles, and W be... ij Here is the kernel function. Since it's an explosion problem, an adaptive smooth length update is used to ensure simulation accuracy. The update formula is as follows:

[0084]

[0085] h is the smooth radius of the particle, which differs from the particle interaction radius mentioned earlier by a constant k, r i =k*h i .

[0086] The schematic diagram of the configuration profile used in the calculation is as follows: Figure 6 As shown, the outermost layer consists of dummy wall particles, while the inner water area is a cube with sides of 0.45m x 0.45m x 0.45m, described using the Mie-Gruneisen equations of state. A sphere with a central radius of 0.04m is designated as a "medicine pack," described using the Jones-Wilkins-Lee equations. All particles are arranged in a simple cubic format within their respective geometric regions, with a particle spacing of 0.01m and a smoothing radius of 0.012m. The computational configuration will contain approximately 800,000 particles, and the simulation will use a fixed time step of 10.-8 s, the algorithm parameter a is 200.

[0087] The total time consumption of SPH program in the first 8000 time steps using the traditional tree search method, the linked list search method and the present method is compared as shown in Figure 7 The time consumption of particle pair search step in the first 8000 time steps of SPH program using the present method is compared as shown in Figure 8 It can be seen from Figure 7 and Figure 8 that the calculation time of the linked list search method increases exponentially with the progress of the calculation, while the calculation time of the tree search method and the present method remains stable with the progress of the simulation, and it can be noted that the calculation time using the present method is much lower than using the tree search method, the average time consumption of the tree search method is 7.2s, while the present method is only 2.2s. It can be seen that the present method greatly improves the efficiency of SPH calculation.

[0088] The above is only the preferred embodiment of the present application, and the present application is not limited to the above embodiments. It can be understood that other improvements and changes directly derived or thought by those skilled in the art without departing from the spirit and concept of the present application should be considered to be included in the protection scope of the present application.

Claims

1. A fast search method for SPH particle pairs suitable for explosion scenarios, characterized in that, The SPH particle pair fast search method includes dividing the computational domain in the explosion scenario into multiple sub-regions, creating a corresponding process for each sub-region, and having multiple processes execute in parallel. The methods executed by each process include: Obtain the physical information of each fluid particle in the corresponding sub-region and its outer region at the current time step. The physical information of the fluid particles includes three-dimensional coordinates, fluid physical quantities, and particle radius of action. Based on the physical information of each fluid particle in the corresponding sub-region and its outer region at the first time step, a particle space layout octree structure is constructed. Based on the particle space layout octree structure, a nearest neighbor list of each fluid particle in the corresponding sub-region is constructed. The nearest neighbor list of each fluid particle in the corresponding sub-region is used as the search range to search for SPH particle pairs that interact within the current time step, and the nearest neighbor list of each fluid particle in the corresponding sub-region is stored. The nearest neighbor list of each fluid particle includes the global ID of all fluid particles located within the potential interaction range of the fluid particle. Starting from the second time step, the physical information of each fluid particle in the corresponding sub-region at the current time step is used to check whether the nearest neighbor list update condition is met. If it is determined that the nearest neighbor list update condition is not met, the nearest neighbor list that has been determined and stored in the historical time steps of each fluid particle in the corresponding sub-region is used as the search range to search for interacting SPH particle pairs in the current time step. If it is determined that the nearest neighbor list update condition is met, the particle space layout octree structure is reconstructed based on the physical information of each fluid particle in the corresponding sub-region and its extended region at the current time step, and the nearest neighbor list of each fluid particle in the corresponding sub-region is updated. The updated nearest neighbor list of each fluid particle in the corresponding sub-region is used as the search range to search for interacting SPH particle pairs in the current time step.

2. The SPH particle pair fast search method according to claim 1, characterized in that, The SPH particle pair fast search method also includes: Based on the maximum particle velocity v of all fluid particles in the entire computational domain at the current time step. max The buffer distance b = 2av is obtained. max dt, where a is a parameter and dt is the time step; Based on the octree structure of the particle space layout and the buffer distance b, a nearest neighbor list of each fluid particle in the corresponding sub-region is constructed. Furthermore, based on the physical information of each fluid particle in the corresponding sub-region at the current time step and the buffer distance b, it is detected whether the nearest neighbor list update condition is met.

3. The SPH particle pair fast search method according to claim 2, characterized in that, Based on the physical information of each fluid particle in the corresponding sub-region at the current time step, combined with the buffer distance b, it is determined whether the nearest neighbor list update conditions are met, including: Each process detects the three-dimensional coordinates X of any fluid particle i within the corresponding sub-region at the current time step. now (i) The three-dimensional coordinates X of the fluid particle i at the time step of the last update of the nearest neighbor list. old The distance ΔX(i) between (i); when there is fluid particle i, the distance ΔX(i) moves. When the process is in progress, it is determined that a target fluid particle exists within the corresponding sub-region; when all fluid particles satisfy the condition... At that time, it is determined that there are no target fluid particles in the sub-region corresponding to the process; If a subregion in the computational domain contains the target fluid particle, the nearest neighbor list update condition is determined to be met; otherwise, the nearest neighbor list update condition is determined not to be met.

4. The fast search method for SPH particle pairs according to claim 2, characterized in that, Based on the particle space layout octree structure and combined with the buffer distance b, a nearest neighbor list for any fluid particle i within the corresponding sub-region is constructed, including: Starting from the root node, the tree nodes in the particle space layout octree structure are traversed sequentially according to a predetermined traversal order. Each tree node in the particle space layout octree structure corresponds to a three-dimensional node control region, and each tree node corresponds to a unique fluid particle. The three-dimensional coordinates of the fluid particles corresponding to each tree node and all its child nodes are within the node control region of the current tree node. For any currently traversed tree node, based on the physical information of fluid particle i and the buffer distance b, detect whether the scope of fluid particle i intersects with the node control area of ​​tree node. When it is determined that the scope of fluid particle i intersects with the node control region of tree node node, it is checked whether the fluid particle j corresponding to tree node node is within the potential interaction range of fluid particle i. When it is determined that fluid particle j is within the potential interaction range of fluid particle i, the global ID of fluid particle j is added to the nearest neighbor list of fluid particle i. Continue to execute the step of traversing the tree nodes in the particle space layout octree structure in a predetermined traversal order. When it is determined that the domain of fluid particle i does not intersect with the node control region of tree node, skip the subsequent traversal of tree node and all its child nodes, and continue to execute the step of traversing the tree nodes in the particle space layout octree structure in the predetermined traversal order.

5. The fast search method for SPH particle pairs according to claim 4, characterized in that, Detecting whether the fluid particle j corresponding to the tree node is within the potential interaction range of the fluid particle i includes: When fluid particle j is at the current time step's three-dimensional coordinate X now (j) and the three-dimensional coordinates X of fluid particle i at the current time step now The distance D between (i) ij <r ij When +b, it is determined that fluid particle j is within the potential interaction range of fluid particle i; otherwise, it is determined that fluid particle j is not within the potential interaction range of fluid particle i; where r ij It is the interaction radius between fluid particle i and fluid particle j.

6. The fast search method for SPH particle pairs according to claim 5, characterized in that, Each fluid particle has its own particle interaction radius, and the interaction radius r between fluid particle i and fluid particle j is... ij =(r i +r j ) / 2; where r i r is the particle radius of action of fluid particle i. j It is the particle radius of action of fluid particle j.

7. The fast search method for SPH particle pairs according to claim 4, characterized in that, The nearest neighbor list for any fluid particle i within a corresponding sub-region, constructed based on the particle space layout octree structure and the buffer distance b, also includes: When the domain of fluid particle i intersects with the node control region of tree node and r i >r j At that time, the step of detecting whether the fluid particle j corresponding to the tree node is within the potential interaction range of the fluid particle i is executed; When the domain of fluid particle i intersects with the node control region of tree node and r i ≤r j At that time, continue to execute the step of traversing the tree nodes in the particle space layout octree structure in the predetermined traversal order; Where, r i r is the particle radius of action of fluid particle i. j It is the particle radius of action of fluid particle j.

8. The fast search method for SPH particle pairs according to claim 4, characterized in that, Based on the physical information of fluid particle i and the buffer distance b, it is determined whether the domain of fluid particle i intersects with the node control region of the tree node, including: When the condition is satisfied If the domain of fluid particle i intersects with the node control region of tree node, then determine that the domain of fluid particle i does not intersect with the node control region of tree node. Among them, (x i ,y i ,z i (x) represents the three-dimensional coordinates of fluid particle i at the current time step. node ,y node ,z node ) represents the center coordinates of the control region corresponding to the tree node, d node r is the radius of the control region corresponding to the tree node. i Let be the particle radius of fluid particle i.

9. The fast search method for SPH particle pairs according to claim 1, characterized in that, The process stores the particle information and local ID of each fluid particle in the corresponding sub-region and its surrounding region at the current time step in a storage container. The particle information includes physical information, a nearest neighbor list, and the three-dimensional coordinates (X) of the fluid particle when constructing the nearest neighbor list. old The outer region of each sub-region includes other adjacent sub-regions surrounding the current sub-region whose distance from the boundary of the current sub-region is less than or equal to r. max The region, r max It is the maximum value of the particle interaction radius of all fluid particles in the entire computational domain at the current time step; The search is conducted using the nearest neighbor list of each fluid particle to find interacting SPH particle pairs within the current time step, including: Construct a hash table to record the global ID of each fluid particle in the storage container and its mapping relationship with its local ID in the storage container; For any fluid particle p within the corresponding sub-region, iterate through each global ID in the nearest neighbor list of fluid particle p, detect the local ID corresponding to the global ID in the hash table, and read the physical information of the corresponding fluid particle q in the storage container based on the local ID. Then, iterate through the physical information of fluid particle p and fluid particle q to determine whether fluid particle p and fluid particle q constitute an SPH particle pair. If the global ID is not included in the hash table, continue to iterate through other global IDs in the nearest neighbor list of fluid particle p.

10. The fast search method for SPH particle pairs according to claim 1, characterized in that, The SPH particle pair fast search method also includes: When fluid particle m moves from sub-region Q corresponding to process Q to sub-region P corresponding to process P, when process Q sends the nearest neighbor lists of each fluid particle m that have been determined and stored in the historical time step to process P, it first concatenates the nearest neighbor lists of each fluid particle m to convert them into a one-dimensional list array, and generates a count array to record the number of fluid particles with global IDs in the nearest neighbor lists of each fluid particle m. Then, it sends the list array and the count array to process P. Process P determines the number of fluid particles with global IDs in the nearest neighbor lists of each fluid particle m based on the count array, and splits and decompresses the list array according to the number of global IDs to obtain the nearest neighbor lists of each fluid particle m that have been determined and stored in the historical time step.

Images