Efficient different-diameter particle search algorithm for discrete element numerical calculation
By optimizing background grid construction and particle numbering management, efficient and accurate contact relationship detection in the reducer particle system is achieved, and the problem of inefficient calculation efficiency in traditional methods is solved, and it is suitable for large-scale particulate matter mechanical analysis.
Patent Information
- Application Number
- CN202510341070.1
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2025-03-21
- Publication Date
- 2025-07-08
- Estimated Expiration
- Not applicable · inactive patent
AI Technical Summary
In discrete element numerical calculation, when traditional methods deal with reducer particles, the number of particles in the background grid increases, resulting in low contact search efficiency, high calculation cost, and difficult to adapt to the computing needs of large-scale systems.
An efficient neighborhood search algorithm based on background mesh optimization is adopted to establish a cube background mesh cell, set the mesh size based on the minimum particle size, and quickly retrieve candidate particles in the particle influence domain through the mesh index, combining particle number management and geometric contact detection to ensure efficient and accurate contact recognition.
It significantly reduces invalid contact judgment and improves calculation efficiency by about 60%-80%, which is suitable for particulate matter mechanical analysis in large-scale systems.
Smart Images
Figure CN120277974A_ABST
Abstract
Description
Technical Field
[0001] The present invention relates to the technical field of particulate computational mechanics and numerical calculation, and particularly relates to a method for searching for particles with different diameters for discrete element numerical calculation. Background Technique
[0002] The discrete element method (DEM) is a widely used numerical technique for simulating the behavior of particulate materials, powders, and other particulate systems. DEM simulations involve modeling individual particles and their interactions, which are governed by contact mechanics, force, and motion equations. These simulations are of great significance in industries such as mining, pharmaceuticals, civil engineering, and additive manufacturing. Therefore, understanding particulate dynamics is crucial for optimizing processes and designs. However, in DEM simulations, how to efficiently identify adjacent particles and calculate their interactions is a key challenge. Traditional search methods usually determine the background grid size based on the maximum particle size and check whether particles are in contact within adjacent grids. When dealing with particles of different diameters, this method will cause a significant increase in the number of particles in the background grid, resulting in a substantial reduction in the contact search efficiency, making the calculation expensive and not suitable for the calculation of large-scale systems. Summary of the Invention
[0003] Aiming at the technical bottleneck of the sharp increase in the computational amount of contact detection and the significant decline in the algorithm efficiency in the discrete element numerical simulation of large-scale multi-diameter particle systems, the present invention innovatively proposes an efficient neighborhood search algorithm based on background grid optimization.
[0004] An efficient search algorithm for particles with different diameters for discrete element numerical calculation proposed by the present invention specifically includes the following key innovative steps:
[0005] S1. Spatial discretization preprocessing: Construct an initial particle system with different diameters, set the size of the cubic background grid unit based on the minimum particle diameter in the particle system with different diameters, and establish a background space grid system covering the entire computational domain; uniquely number all particles in the particle system with different diameters, and establish a mapping relationship between the spatial coordinate information of each particle and the background space grid system;
[0006] S2. Determine the particle influence domain: According to the geometric center coordinates of each particle, determine the adjacent grid area corresponding to its spatial influence domain; quickly retrieve all candidate particles that may come into contact within this influence domain through grid indexing, and establish a fast search mechanism based on grid topology;
[0007] S3. Generate candidate contact pairs: Traverse the influence domain grids of each particle, store the particle numbers located in the same or adjacent grids into the candidate pair list, and automatically exclude duplicate pairings by comparing the size relationship of particle numbers to ensure that each contact pair is generated only once;
[0008] S4. Precise Contact Judgment: Perform geometric contact detection on the particle pairs in the candidate pair list, calculate the distance between the centers of the two spheres, and compare it with the sum of the radii of the two spheres to determine whether the particle pairs are in contact, and output the final contact list.
[0009] Preferably, step S1 includes:
[0010] S11. Continuously number the particles in the non-uniform particle system from 1 to n;
[0011] S12. Determine the minimum particle diameter D in the non-uniform particle system min , and set the side length L of the cubic background grid unit to: L = 1.1D min ;
[0012] S13. According to the particle center position and the particle radius R i , determine the background grid region boundary as follows:
[0013]
[0014] where i ∈ (1,..., n); X min , Y min , Z min respectively represent the minimum boundary coordinates of the region, and X max , Y max , Z max respectively represent the maximum boundary coordinates of the region;
[0015] S14. Number the background grid according to the background grid region boundary to form a three-dimensional grid index system of (x i , y i , z i ):
[0016] Number sequentially from 1 to x in the x direction i , where,
[0017] Number sequentially from 1 to y in the y direction i , where,
[0018] Number sequentially from 1 to z in the z direction i , where,
[0019] According to the set background grid, combined with the center position of each particle, determine the position relationship between the particle and the background space grid system.
[0020] Preferably, step S2 includes:
[0021] S21. According to the particle center position Particle radius R i and the size of the background grid to determine the control domain range of the particle. For particle P i The grid range covered by the control domain Ω is calculated as follows:
[0022]
[0023] In the formula is the floor function, where "·" represents the formula within, are the minimum coordinates of the control domain in the x, y, and z directions respectively, are the maximum coordinates of the control domain in the x, y, and z directions respectively;
[0024] S22. Use the spatial coverage judgment criterion to label the particles for each background grid cell covered by the particle control domain;
[0025] S23. Based on the background grid system (indexed as x i , y i , z i ) constructed based on the three-dimensional Cartesian coordinate system, through the grid traversal algorithm, sequentially count the set of particles marked in each grid cell, and construct a data structure containing the grid spatial coordinates and the associated particle list, that is, the grid-particle spatial association data.
[0026] Preferably, step S3 includes:
[0027] Based on the grid-particle spatial association data generated in step S23, use the following progressive screening strategy to construct a list of non-repetitive contact candidate pairs:
[0028] S31. Perform a spatial traversal on the contact detection control domain for each target particle P i :
[0029] Obtain all the background grid cells associated with the control domain where particle P i is located, traverse each associated unit grid, collect all the particle sets stored in the unit, pair all the particles stored in the unit with the target particle P i and record the pairing information in the preliminary contact candidate list. This list is a (n×m) dynamic array, where n is the total number of system particles and m is the total number of particles paired with particle P i ;
[0030] S32. Establish a dynamic marking anti-duplication mechanism:
[0031] Establish a particle status marking array. First, initialize an array Mark[n] containing n elements, where n is the total number of system particles, and the initial value of each element is 0;
[0032] When the particle P j forms a pairing relationship with P i the array element Mark[j]=1 is updated synchronously;
[0033] During pairing screening, a new contact pair is allowed to be established only when Mark[j]==0. After each pairing operation, the corresponding flag bit is updated immediately to ensure real-time status synchronization;
[0034] S33. State reset mechanism:
[0035] After all neighborhood detections of the current target particle P i are completed, the Mark array is cleared to ensure entering the processing flow of the next particle;
[0036] Finally, the pairing relationship is established for the particle numbers marked in each background grid by the binary method, that is, only the relationships that meet the condition of particle number i<j are paired. After traversing all the grids, a candidate pair list is formed.
[0037] Preferably, step S4 includes:
[0038] For the candidate pair list established in step S3, perform an accurate contact determination. The specific steps are as follows:
[0039] S41. Traverse the paired particles of P i Compare the distance between the particle centers and the sum of the radii of the two particles R i +R j If the particles are marked as in contact, otherwise the particles not in contact are removed;
[0040] S42. Traverse the paired arrays of each particle in turn, repeat step S41, obtain the contact particle pairs of each particle, and output the final contact list
[0041] Beneficial effects:
[0042] This method is mainly used to evaluate and analyze the contact relationships of particulate matter in complex mechanical environments. Especially in the discrete element simulations involving particles of different sizes (heterogeneous particles), it can efficiently and accurately detect the contact relationships between particles. The present invention can be widely applied to scenarios such as mechanical analysis of particulate matter, particulate flow simulation, and multiphase flow numerical calculation in the fields of civil engineering, geotechnical engineering, chemical engineering, materials science, etc.
[0043] The innovation of this method lies in constructing a background grid structure based on the minimum particle size of the system, and realizing the rapid positioning of contacting particle pairs by establishing a grid space indexing mechanism. Compared with the traditional global search algorithm, it can reduce about 60%-80% of the ineffective contact judgments, significantly improving the operation efficiency while ensuring the calculation accuracy. Description of the Drawings
[0044] Figure 1 is the flow chart of the search and inspection working principle for the non-uniform particle system;
[0045] Figure 2 is the two-dimensional schematic diagram of the particle control domain division,
[0046] In the figure, the black circles are particles, the red dotted boxes are the particle control domains set to prevent missed detections, and the shaded parts are the contact detection control domains associated with the red dotted boxes;
[0047] Figure 3 is the background grid control domain occupied by particle P i (particle number i),
[0048] In the figure, j, k, e, f, k, g are the candidate contacting particles associated within the background grid. Detailed Implementation Manner
[0049] The technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings in the embodiments of the present invention.
[0050] The present invention discloses an efficient search algorithm for non-uniform particles used in discrete element numerical calculations, as Figure 1 shown, which shows the main operation logic of this algorithm:
[0051] This algorithm specifically includes the following steps:
[0052] Step S1: Perform spatial discretization preprocessing, specifically:
[0053] S11: Construct an initial non-uniform particle system, assign unique numbers to all particles in the system, the total number of particles is n, and the particles are continuously numbered from 1 to n;
[0054] S12: Set the size of the cubic background grid unit based on the minimum particle diameter in the non-uniform particle system, and establish a background space grid system covering the entire calculation domain; determine that the minimum particle diameter in the particle system is D min , and set the side length L of the cubic background grid unit to be slightly larger than the minimum particle diameter: L = 1.1D min . As Figure 2 shows the two-dimensional schematic diagram of the particle control domain division, ensuring that the grid size is slightly larger than the minimum particle diameter, so as to provide redundant space for subsequent contact search;
[0055] S13. To ensure that the background grid covers all particles and to establish a particle control domain subsequently, the background grid area needs to be appropriately enlarged, that is, the background grid area is extended, so as to avoid missing detection caused by particles moving beyond the computational domain. Extract the center coordinates of all particles and the particle radius R i , calculate the range of the extended grid area; according to the particle center position and the particle radius R i , extend the grid area outward by the length dimension of two grids (2L).
[0056]
[0057] where i ∈ (1,..., n); X min , Y min , Z min respectively represent the minimum boundary coordinates of the area, and X max , Y max , Z max respectively represent the maximum boundary coordinates of the area.
[0058] S14. Based on the extended grid boundary, systematically divide and number the background grid, and number the background grid according to the background grid area boundary: uniformly divide the cubic unit along the x, y, and z directions with side length L Figure 2 is a two-dimensional schematic diagram, and the background grid is numbered using two-dimensional Cartesian coordinates (x i , y i ), where x i and y i are the index values of the grid in the x and y directions respectively. In practical applications, the z i direction numbering needs to be added according to the same rule in three-dimensional space to form a three-dimensional grid index system of (x i , y i , z i ).
[0059] Number sequentially from 1 to x n in the x direction, where
[0060] Number sequentially from 1 to y n in the y direction, where
[0061] Number sequentially from 1 to z n in the z direction, where
[0062] According to the set background grid, combined with the center position of each particle Determine the positional relationship between the particles and the background grid, that is, establish a mapping relationship between the spatial coordinate information of each particle and the background spatial grid system.
[0063] Since the side length of the background grid is 1.1 times the minimum particle diameter, particles often occupy multiple background grids; to better locate the positional relationship between the grid and the particles, the method of the present invention uses the method of particle control domain for processing, see step S2 for details:
[0064] S2. According to Figure 2 The two-dimensional schematic diagram of particle control domain division shown, to avoid missing detection of boundary particles, it is necessary to expand on the basis of the particle control domain. Specifically:
[0065] S21. According to the geometric center coordinates of each particle, determine the adjacent grid area corresponding to its spatial influence domain; according to the particle center position particle radius R i and the size of the background grid to determine the range of the particle control domain. To reduce the number of particle contact detections in the discrete element calculation, the particle control domain should be appropriately increased. For the convenience of the method of the present invention, a cube-shaped control domain is used so that it can fit the background grid, and the grid control domain is expanded outward by half a grid size (L / 2) in each of the x, y, and z directions. Then, for particle P i The grid range covered by the control domain Ω is calculated as follows:
[0066]
[0067] In the formula is the floor function, where "·" represents the formula inside, and the control domain expansion amount L is used to avoid missing detection of boundary particles. are the minimum coordinates of the particle control domain in the x, y, and z directions respectively. are the maximum coordinates of the particle control domain in the x, y, and z directions respectively. The particle calculation domain obtained from the above calculation corresponds to Figure 2 the red dashed box in, and for the convenience of retrieving particles through the grid subsequently, this area needs to be further topologically mapped to the covered background grid, which is the contact detection control domain of the particle, such as Figure 2 the shaded area around each particle in.
[0068] S22. Use the spatial coverage judgment criterion to label the particle numbers for each background grid cell covered by the particle contact detection control domain, that is, mark the particle numbers corresponding to each grid within the corresponding particle control domain (the grid covered by the shaded area). The specific implementation is as Figure 2 shown: For particle P i (number i) and particle P jThe grid cells covered by the contact detection control field of (serial number j) are respectively marked with their corresponding particle numbers i and j.
[0069] S23. A background grid system (indexed as x i , y i , z i ) is constructed based on a three-dimensional Cartesian coordinate system. Through a grid traversal algorithm, the particle sets marked in each grid cell are sequentially counted, and a data structure containing the grid space coordinates and the associated particle list is constructed, that is, the grid-particle space association data. For example, Figure 2 for the particles P i and P j in [reference], the detection control fields overlap at the grid cells (x = 3, y = 5) and (x = 4, y = 5). Therefore, the associated particle lists of these two grid cells will respectively record the index numbers i and j.
[0070] S3. Pre-screening of neighboring particles and generation of candidate contact pairs: Traverse the influence domain grids of each particle, store the particle numbers located in the same or adjacent grids into the candidate pair list, and automatically exclude duplicate pairings by comparing the size relationship of particle numbers through the binary search method to ensure that each contact pair is generated only once;
[0071] Based on the grid-particle space association data generated in step S23, the following progressive screening strategy is adopted to construct a non-repetitive contact candidate pair list:
[0072] S31. Perform a spatial traversal on the contact detection control field of each target particle P i to pre-screen the particle set that may come into contact with it:
[0073] Obtain all the background grid cells associated with the control field where the particle P i is located, traverse each associated unit grid, collect all the particle sets stored in the unit, pair all the particles stored in the unit with the target particle P i , and enter the pairing information into the preliminary contact candidate list. This list is a dynamic array of (n × m), where n is the total number of particles in the system and m is the total number of particles paired with the particle P i . For example, as shown in Figure 3 for the particle P i (serial number i), the detection control field contains a particle set (j, k, e, g, f) that may come into contact with it, and these numbers will establish a temporary association with Pi in the traversal order;
[0074] S32. Establish a dynamic marking anti-duplication mechanism:
[0075] To achieve the uniqueness control of the pairing process, a particle status marking array is established. First, initialize an array Mark[n] containing n elements, where n is the total number of particles in the system, and the initial value of each element is 0;
[0076] When the particle P j and P i establish a pairing relationship, synchronously update the array element Mark[j]=1;
[0077] During pairing screening, a new contact pair is allowed to be established only when Mark[j]==0. After each pairing operation, the corresponding flag bit is immediately updated to ensure real-time status synchronization;
[0078] S33. Status reset mechanism:
[0079] After completing all neighborhood detections of the current target particle P i perform the operation of clearing the Mark array to ensure entering the processing flow of the next particle.
[0080] Finally, use the binary method to establish a pairing relationship for the particle numbers marked in each background grid, that is, only pair the relationships that meet the condition of particle number i<j. After traversing all grids, form a list of pairing relationships established with particle numbers, that is, the candidate pair list.
[0081] S4. Accurate contact determination: Perform geometric contact detection on the particle pairs in the candidate pair list, calculate the distance between the two sphere centers and compare it with the sum of the radii of the two spheres to determine whether the particle pairs are in contact, and output the final contact list.
[0082] For the candidate pair list established in step S3, perform accurate contact determination. The specific steps are as follows:
[0083] S41. Traverse the paired particles of P i particles, compare the distance between the particle centers and the sum of the radii R i +R j between them. If then the particle is marked as in contact, otherwise the particle is not in contact and is excluded;
[0084] S42. Traverse the paired arrays of each particle in turn, repeat step S41, obtain the contact particle pairs of each particle, and output the final contact list.
Claims
1. An efficient search algorithm for particles with different diameters in discrete element numerical calculations, characterized in that It includes the following steps: S1. Spatial discretization preprocessing: Construct an initial non-uniform particle system, set the size of the cubic background grid cell based on the minimum particle size in the non-uniform particle system, and establish a background space grid system covering the entire computational domain; uniquely number all the particles in the non-uniform particle system, and establish a mapping relationship between the spatial coordinate information of each particle and the background space grid system; S2. Determine the particle influence domain: According to the geometric center coordinates of each particle, determine the adjacent grid area corresponding to its spatial influence domain; quickly retrieve all candidate particles that may come into contact within this influence domain through grid indexing, and establish a fast search mechanism based on grid topology; S3. Generate candidate contact pairs: Traverse the influence domain grid of each particle, store the particle numbers located in the same or adjacent grids into the candidate pair list, and automatically exclude duplicate pairings by comparing the size relationship of the particle numbers to ensure that each contact pair is generated only once; S4. Accurate contact determination: Perform geometric contact detection on the particle pairs in the candidate pair list, calculate the distance between the two sphere centers and compare it with the sum of the two sphere radii to determine whether the particle pair is in contact, and output the final contact list.
2. A high-efficiency search algorithm for non-uniform particles in discrete element numerical calculation according to claim 1, characterized in that: Step S1 includes: S11. Continuously number the particles in the non-uniform particle system from 1 to n; S12. Determine the minimum particle diameter D in the non-uniform particle system min , and set the side length L of the cubic background grid cell as: L = 1.1D min ; S13. Determine the background grid region boundary according to the particle center position and the particle radius R i as follows: where \(i\in(1,\ldots,n)\); \(X\) min , \(Y\) min , \(Z\) min respectively represent the minimum boundary coordinates of the region, \(X\) max , \(Y\) max , \(Z\) max respectively represent the maximum boundary coordinates of the region; S14. Number the background grid according to the background grid area boundary to form a three-dimensional grid index system of (x i , y i , z i ): Number them sequentially from 1 to x in the x direction i , where Number sequentially from 1 to y in the y direction i , where Number them sequentially from 1 to z in the z direction i , where According to the set background grid, combined with the center position of each particle, determine the position relationship between the particle and the background space grid system.
3. A high-efficiency search algorithm for non-uniform particles in discrete element numerical calculation according to claim 2, characterized in that: Step S2 includes: S21. Determine the control domain range of the particle according to the particle center position Particle radius R i and the size of the background grid, and the control domain range of the particle P i The grid range covered by the control domain Ω is calculated as follows: x direction: y direction: z direction: In the formula is the floor function, where "·" represents the formula within, are respectively the minimum coordinates of the control domain in the x, y, and z directions, are respectively the maximum coordinates of the control domain in the x, y, and z directions; S22. Use the spatial coverage judgment criterion to label the particle numbers for each background grid cell covered by the particle control domain; S23. Background grid system (indexed as x i , y i , z i ) ， Based on the three-dimensional Cartesian coordinate system, by means of the grid traversal algorithm, successively count the set of marked particles in each grid cell, and construct a data structure containing the grid space coordinates and the associated particle list, that is, the grid-particle space association data.
4. A high-efficiency search algorithm for non-uniform particles in discrete element numerical calculation according to claim 3, characterized in that: Step S3 includes: Based on the grid-particle spatial association data generated in step S23, use the following progressive screening strategy to construct a non-repetitive contact candidate pair list: S31. Perform a spatial traversal on each target particle P i for the contact detection control region: Obtain particle P i All background grid cells associated with the control domain where it is located. Traverse each associated cell grid, collect all the particle sets stored in the cell, and pair all the particles stored in the cell with the target particle P i Perform pairing, and enter the pairing information into the preliminary contact candidate list. This list is a (n×m) dynamic array, where n is the total number of system particles, and m is the total number of particles paired with particle P i The total number of paired particles; S32. Establish a dynamic marking anti-duplication mechanism: Establish a particle state marking array. First, initialize an array Mark[n] containing n elements, where n is the total number of particles in the system, and the initial value of each element is 0; When the particle P j and P i establish a pairing relationship, synchronously update the array element Mark[j]=1; When pairing and screening, only when Mark[j] == 0 is a new contact pair allowed to be established. Immediately update the corresponding marked bit after each pairing operation to ensure real-time state synchronization; S33. State reset mechanism: Complete the detection of all the neighborhoods of the current target particle P i After that, perform the operation of clearing the Mark array to ensure entering the processing flow of the next particle; Finally, use the binary method to establish a pairing relationship for the particle numbers marked in each background grid, that is, only pair the relationships that meet the condition of particle number i < j. After traversing all the grids, form a candidate pair list.
5. A high-efficiency search algorithm for non-uniform particles in discrete element numerical calculation according to claim 4, characterized in that: Step S4 includes: For the candidate pair list established in step S3, perform accurate contact determination. The specific steps are as follows: S41. Traverse P i The paired particles of the particles, and compare the center distance of the particles with the sum R of the radii of the two particles i +R j to determine the size relationship therebetween. If the particles are marked as in contact; otherwise, the particles that are not in contact are removed. S42. Traverse the pairing array of each particle in turn, repeat step S41, obtain the contact particle pairs of each particle, and output the final contact list.