A fast particle packing method based on triangulated mesh and breadth-first search algorithm

Abstract

A fast particle stacking method for discrete element model based on triangular mesh and breadth-first search algorithm is proposed. The method comprises the following steps: S1: Constructing an initial triangular mesh. The target geometric region is divided using a triangular mesh partitioning algorithm, and mesh node information and triangular mesh element information are saved to complete the geometric model initialization; S2: Mesh heterogeneity processing; S3: Generating mesh traversal paths. Based on the perturbed triangular mesh topology information, the generation path of particle circles inside the mesh is obtained; S4: Particle generation inside the mesh. The triangular mesh elements are divided into boundary meshes and non-boundary meshes. Based on the constraints of circle tangency to circle and circle tangency to mesh edges, particle circles are generated inside each triangular mesh; S5: Mesh vertex particle filling. The boundary vertices and internal vertices of the triangular mesh are distinguished. The maximum-size filling particle circle is calculated and generated, constrained by the particle circles not overlapping and being tangent to surrounding generated particle circles; S6: Post-processing of the algorithm.

Description

Technical Field

[0001] This invention belongs to the field of particle system generation, such as smooth particle hydrodynamics (SPH) and discrete element method (DEM), and specifically focuses on a two-dimensional dense particle packing algorithm based on geometric construction. Background Technology

[0002] The Discrete Element Method (DEM), as a numerical method for studying the macroscopic and microscopic mechanical behavior of soil and rock masses and granular materials, has demonstrated irreplaceable value in key geotechnical engineering problems. However, limitations of traditional methods in terms of efficiency, contact quality, boundary adaptability, and heterogeneity hinder the widespread application of the DEM in engineering.

[0003] Existing particle packing algorithms are mainly divided into two categories: dynamic methods and geometric construction methods. Dynamic methods generate samples by simulating the real physical processes of particles under external forces such as gravity and compression. Their advantage lies in achieving extremely high density, but their computational efficiency is extremely low. Geometric construction methods aim to quickly place particles through pure geometric calculations, avoiding the time-consuming dynamic relaxation process, but they suffer from many problems such as low contact quality and coordination number, difficulty in handling complex boundaries, and lack of heterogeneity. Therefore, the development of numerical simulation in geotechnical engineering urgently requires a method that can efficiently, with high quality, and with high fidelity generate complex discrete element samples.

[0004] Against this backdrop, the particle fast packing method based on triangular mesh and breadth-first search algorithm specifically addresses the problems related to low coordination number and low density, and heterogeneous and irregular modeling in existing particle packing algorithms for geotechnical engineering. Specifically, this is reflected in the following aspects:

[0005] (1) Through precise mathematical solutions, the close contact between adjacent particles is fundamentally guaranteed, the coordination number and packing density of particles are increased, and the generated sample is closer to the pre-compacted state.

[0006] (2) The introduction of triangular mesh technology solves the modeling problem of irregular boundaries in geotechnical engineering. At the same time, by introducing random perturbation, the regularity of the initial mesh can be broken, thereby introducing controllable randomness in particle size and spatial distribution, and better simulating the heterogeneous characteristics of natural materials. Summary of the Invention

[0007] Objective: The core contribution of this invention lies in proposing a geometric stacking algorithm based on triangular meshes and breadth-first search. The method of this invention ensures strict contact between adjacent particles by dynamically establishing and solving the circle-edge and circle-circle tangency equations for each triangular element, significantly improving the coordination number. Supporting key technologies such as vertex filling and random mesh perturbation further guarantee the algorithm's high adaptability to arbitrarily complex geometric domains, the heterogeneity of the generated results, and high efficiency.

[0008] The technical solution of this invention is: a fast stacking method for discrete element particles based on triangular mesh and breadth-first search algorithm; including the following steps:

[0009] 1) Constructing the initial triangular mesh: Construct the initial triangular mesh using a triangular meshing algorithm.

[0010] 2) Mesh heterogeneity: Add noise to the triangular mesh to make it heterogeneous.

[0011] 3) Obtaining the path for generating circles inside the triangular mesh: Path generation is achieved using the breadth-first search algorithm.

[0012] 4) Generation of circles inside the mesh: Fill the triangular mesh with particles sequentially according to the generation path.

[0013] 5) Vertex filling: Divide the mesh vertices into internal vertices and boundary vertices, and fill the internal vertices with particles.

[0014] 6) Post-processing of the algorithm: Save the data and plot it.

[0015] Specifically, the following steps are included:

[0016] Step 10: Create an initial mesh as needed, and save the mesh nodes and triangular mesh information for use in the concretization of the geometric model.

[0017] Step 11: Based on the divided grid size, randomly perturb the coordinates of non-boundary grid nodes to change the node coordinates and reset the model.

[0018] Step 12: Obtain the internal circle generation path using a breadth-first search algorithm. Based on the mesh information of the perturbed model, select the geometric center mesh as the initial mesh. Starting from the initial mesh, use a breadth-first search algorithm to obtain the internal circle generation path.

[0019] Step 13: Following the path obtained by the breadth-first search algorithm in the previous step, generate circles sequentially within the grid. Divide adjacent grids into two categories: boundary grids and internal grids, and calculate and generate the internal circles for each grid, while saving the information of the generated internal circles.

[0020] Step 14: Use the function to screen out the boundary vertices and internal vertices of the triangular mesh. For internal vertices, calculate the maximum number of particles that can be filled at the vertex according to the principle of maximum and non-overlapping, and generate a circle at the vertex to simulate vertex filling, thereby reducing porosity. At the same time, save the information of the generated vertex filling circle.

[0021] Step 15: Use the information of the inscribed circle and the vertex filling circle to draw the final particle generation diagram and save the particle information. Then, use the saved information to draw the gradation distribution curve, analyze the particle size distribution characteristics, and draw the particle size distribution diagram.

[0022] This invention deeply integrates computational geometry and high-performance scientific computing, proposing a rapid particle packing method for discrete element method (DEM) simulations based on triangular meshes and breadth-first search algorithms. Utilizing techniques such as triangular mesh generation, mesh heterogeneity, breadth-first path traversal, and hierarchical particle filling, it balances computational efficiency with particle packing quality. Its applications cover multiple fields including geotechnical engineering, geological hazard simulation, and mining engineering, serving as a core processing tool before discrete element method (DEM) numerical simulations. This method is adaptable to various geometric models, including regular regions, irregular boundaries, and locally refined modeling. It can efficiently prepare particle packing models with low porosity, high contact accuracy, and reasonable coordination numbers, meeting the needs of applications ranging from small-scale sample test simulations to large-scale engineering numerical simulations. Its core service is for the DEM preprocessing stage in geotechnical and materials science fields where rapid particle packing is urgently needed.

[0023] The beneficial effects of this invention are as follows: Particle packing based on triangular meshes and breadth-first search algorithms significantly improves particle packing quality, such as coordination number and contact accuracy, while ensuring computational speed. The complete algorithm system has strong practicality and is applicable to complex geometries, providing an efficient particle generation tool for particle numerical simulation. It is a powerful tool for efficiently generating high-quality samples in large-scale discrete element simulations, and is particularly suitable for engineering applications that require handling complex geometric boundaries and local refinement. This invention integrates computational geometry and ergonomic search algorithms, balancing particle generation efficiency and packing quality, effectively improving the coordination number and contact accuracy of the particle system, reducing porosity, and adapting to complex geometric boundary models. It is suitable for the rapid preparation of particle samples for large-scale discrete element simulations in the fields of geotechnical engineering and materials science. Attached Figure Description

[0024] Figure 1 Flowchart of numerical simulation of particle packing based on triangular mesh and breadth-first search algorithm;

[0025] Figure 2A , Figure 2B These are comparison images before and after mesh heterogeneity;

[0026] Figure 3A , 3B These are comparison charts of mesh generation quality; namely, quality statistics charts of the original mesh and the heterogeneous network.

[0027] Figure 4 Schematic diagram of boundary grid and non-boundary grid;

[0028] Figure 5 The five figures in the Chinese imagery (AE) illustrate different categories of internal circle generation. One type of boundary mesh occurs when two adjacent meshes already contain internal circles; in this case, a new circle must be tangent to the internal circles of both neighboring meshes. (See [link to image]). Figure 5a. Boundary Mesh: If only one neighboring mesh has already generated an inner circle, then the new circle must be tangent to the inner circle of one of the neighboring meshes when it is generated. See section a. Figure 5 b.

[0029] Non-boundary meshes: If three neighboring meshes have already generated internal circles, then the new circle must be tangent to the internal circles of the three neighboring meshes when it is generated. See [link to documentation]. Figure 5 c.

[0030] Non-boundary meshes: If two neighboring meshes have already generated internal circles, then the new circle generated must be tangent to the internal circles of the two neighboring meshes. See [link to documentation]. Figure 5 d.

[0031] Non-boundary meshes: If only one neighboring mesh has generated an internal circle, then the new circle must be tangent to the internal circle of a neighboring mesh and tangent to the other two edges of the mesh it belongs to. See [link to documentation]. Figure 5 e.

[0032] Figure 6 Selecting the vertex-filled circle. Detailed Implementation

[0033] The present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. These embodiments are implemented based on the technical solution of the present invention and provide a detailed implementation process. However, the scope of protection of the present invention is not limited to the following embodiments.

[0034] This embodiment discloses a rapid particle packing method for discrete element numerical simulations based on triangular meshes and a breadth-first search algorithm. It is mainly applied to the preprocessing stage of discrete element numerical simulations in geotechnical engineering and materials science, enabling the rapid generation of high-quality particle packing samples, adaptable to complex geometric boundary models. The overall process is described below. Figure 1 This is a flowchart of a numerical simulation of particle packing based on a triangular mesh and a breadth-first search algorithm.

[0035] This embodiment selects a conventional rectangular geometric region as the simulation object, with the goal of generating a spherical particle packing system with high density, reasonable coordination number, and controllable particle size distribution. The complete execution steps are as follows:

[0036] Step 10: Use the triangulation algorithm to initially divide the triangular mesh and save the mesh and node information.

[0037] Construct the initial triangular mesh to complete the model initialization;

[0038] The rectangular target region is uniformly meshed using the general Delaunay triangulation algorithm. The mesh size is set according to the simulation accuracy requirements. After the meshing is completed, the two-dimensional coordinate information and node number of all mesh nodes are uniformly stored, as well as the three node numbers and element numbers contained in each triangular mesh element, to establish a complete geometric model data file, providing basic geometric data for subsequent calculations.

[0039] Step 11: Mesh Heterogeneity Processing and Mesh Quality Assessment

[0040] The node random perturbation iterates through all grid nodes, filtering out non-boundary grid nodes. Based on the grid size, the coordinates of the non-boundary grid nodes are randomly perturbed using the average distance from the node to all remaining nodes containing that node, thereby changing the node coordinates. The node random perturbation pi is given by the following formula:

[0041] In the formula, The noise intensity coefficient, The length of the local feature (calculated as the average distance to adjacent nodes). The vector is a random direction vector, with its components uniformly distributed within the interval [−0.5, 0.5]. The perturbation displacement is superimposed onto the original node coordinates, updating the positions of all non-boundary nodes, thus completing the heterogeneous transformation of the triangular mesh and obtaining a new heterogeneous geometric model. Visualization and mesh quality evaluation are performed by plotting the original triangular mesh before perturbation and the heterogeneous triangular mesh after perturbation (as shown in Figure 2). Using the degree to which the triangular mesh elements approximate equilateral triangles as an evaluation index, the mesh quality coefficient of each triangular element is quantitatively calculated, generating a mesh partitioning quality comparison diagram (as shown in Figure 3). The heterogeneous mesh breaks the uniformity of the regular mesh, more closely resembling the disordered characteristics of natural particle packing, while the overall mesh quality meets the computational requirements.

[0042] Step 12: Obtaining the path for generating the inner circle

[0043] Utilizing the node information of the perturbed model, an edge is represented by two node numbers. The node numbers contained within a triangular mesh are used to determine whether an edge belongs to a particular triangular mesh, thus obtaining the triangular mesh number to which each edge belongs. Furthermore, by sharing edges, the neighboring mesh information of each triangular mesh can be obtained. This invention uses a breadth-first search algorithm to obtain the generation path of the circle within the mesh. To improve computational efficiency, the geometric center mesh is selected as the initial mesh. Based on the coordinates of the three nodes of the mesh, the inscribed circle information of the triangular mesh is calculated. Using this triangular mesh as the base point, its neighboring meshes are searched. If a neighboring mesh does not exist, no operation is performed; instead, the process jumps directly to the next remaining neighboring mesh. If a neighboring mesh exists, its number is sequentially appended to the base point triangular mesh number to form a numbering sequence. After the neighboring meshes of the base point triangular mesh have been searched, the same operation is performed on the remaining triangular meshes in the sequence according to the newly formed numbering sequence, ultimately generating a sequence containing all triangular mesh numbers, i.e., the generation path of the circle within the mesh.

[0044] Step 13: Following the triangular mesh numbering order obtained in the previous step, use the obtained neighbor mesh information to determine whether each triangular mesh has neighbor meshes, and then divide the triangular meshes into boundary meshes and non-boundary meshes. Figure 4 Since the inner circle has already been generated in the base point triangular mesh, and the generation of the inner circle is carried out according to the generation path obtained in step 3, the inner circle of each triangular mesh, except for the base point triangular mesh, must satisfy the condition that it is tangent to the circle in the previous triangular mesh when it is generated. Then, based on whether there is an inscribed circle in the remaining neighboring meshes, the boundary mesh is further divided into two cases: two neighboring meshes have inscribed circles and only one neighboring mesh has an inscribed circle. The non-boundary mesh is further divided into three cases: three neighboring meshes have inscribed circles, only two neighboring meshes have inscribed circles, and only one neighboring mesh has an inscribed circle. Specifically, the following cases are identified:

[0045] (1) Boundary Mesh: If two adjacent meshes have already generated internal circles, the new circle must be tangent to the internal circles (thick gray solid lines) of the two neighboring meshes and to the remaining edge (thick gray solid line) of the mesh it belongs to. See details below. Figure 5 a.

[0046] (2) Boundary Mesh: If only one neighboring mesh has generated an inner circle, then when generating a new circle, it must satisfy the following conditions: it must be tangent to the inner circle (thick gray solid line circle) of a neighboring mesh and tangent to the other two edges (thick gray solid lines) of the mesh it belongs to. See the details below. Figure 5 b.

[0047] (3) Non-boundary mesh: If the inner circles of the three neighboring meshes have already been generated, the new circle must be tangent to the inner circles (thick gray solid line circles) of the three neighboring meshes when it is generated. See the details below. Figure 5 c.

[0048] (4) Non-boundary meshes: If two neighboring meshes have already generated internal circles, the new circle must be tangent to the internal circles (thick gray solid lines) of the two neighboring meshes and to the remaining edge (thick gray solid line) of the mesh it belongs to. See details below. Figure 5 d.

[0049] (5) Non-boundary mesh: If only one neighboring mesh has generated an inner circle, the new circle must be tangent to the inner circle (thick gray solid line circle) of a neighboring mesh and to the other two edges (thick gray solid lines) of the mesh it belongs to. See details. Figure 5 e.

[0050] For different situations, the internal circles of the triangular mesh are generated based on constraints such as circles being tangent to each other and circles being tangent to edges. Figure 5 The gray dashed border is round in the middle of the image.

[0051] Step 14: Represent edges as node pairs, with smaller node numbers first and larger node numbers second. Extract the three edges of all numbered triangular meshes, count the frequency of each edge, and identify the boundary edges that appear only once. Obtain the boundary edge node numbers. Further combining these with all node numbers yields the internal node numbers, which are the vertex numbers of the internal triangular meshes. For each internal node, find all triangular meshes containing that node, obtain the inscribed circle information from these meshes, and filter out circles with radii greater than zero. Then, iterate through all combinations of three circles, generating node-filling circles with the constraint of being tangent to all three circles and not overlapping with any previously generated circles. Compare the radius of each generated circle with the radius of the previously generated circle, and select the circle with the largest radius as the final node-filling circle. Specific selection is as follows... Figure 6 In this context, thick gray solid circles represent circles that have already been generated in neighboring meshes, thin black dashed circles represent circles that may be generated but do not meet the requirements, and circles with black borders and gray fills represent circles that meet the requirements.

[0052] Step 15: Use the information of the inscribed circle and the vertex filling circle to draw the final particle generation diagram, and save the particle information. Then use the saved information to draw the gradation distribution curve, analyze the particle size distribution characteristics, and draw the particle size distribution diagram.

[0053] Implementation Results Description

[0054] This embodiment, after generating particle stacking using the above method, has the following advantages compared to traditional random placement and regular grid filling methods:

[0055] Combining breadth-first search with ordered traversal, the algorithm is computationally efficient and can quickly generate particles for large-scale grid regions, making it suitable for large-scale discrete element models.

[0056] The combination of mesh heterogeneity and partition constraints generates particles with strong particle packing disorder, significantly improving coordination number and particle contact accuracy, and better reflecting the real packing state of natural soil and powder materials.

[0057] By using a dual approach of mesh interior filling and vertex filling, the system porosity is effectively reduced, resulting in higher packing density.

[0058] The algorithm has a unified logic and can be flexibly adapted to complex geometric models such as polygons and irregular boundaries, making it highly versatile.

[0059] The specific embodiments described in this specification are only used to explain the technical solutions of the present invention, and are not intended to limit it. The present invention is not limited to the above-described embodiments, and its scope of protection covers equivalent transformations, parameter substitutions, step fine-tuning, structural improvements, and new implementation methods that can be obtained by those skilled in the art based on the core technical ideas of the present invention without creative effort.

Claims

1. A fast stacking method for discrete element particles based on triangular mesh and breadth-first search algorithm, characterized in that, Includes the following steps: S1: Construct the initial triangular mesh. Use a triangular meshing algorithm to divide the target geometric region, save the mesh node information and triangular mesh cell information, and complete the geometric model initialization; S2: Mesh heterogeneity processing. Randomly perturb the coordinates of the non-boundary mesh nodes of the initial triangular mesh, reset the geometric model, and realize triangular mesh heterogeneity; S3: Generate mesh traversal path. Based on the perturbed triangular mesh topology information, select the triangular mesh at the geometric center as the initial mesh, and use a breadth-first search algorithm to traverse all triangular mesh cells to obtain the generation path of the internal particle circles of the mesh; S4: Mesh Internal Particle Generation. Following the generation path, each triangular mesh cell is traversed sequentially. Combining the particle generation state of adjacent mesh cells, the triangular mesh cells are divided into boundary meshes and non-boundary meshes. Based on the constraints of circles being tangent to each other and circles being tangent to mesh edges, particle circles are generated inside each triangular mesh, and the internal particle circle information is saved. S5: Mesh Vertex Particle Filling. The boundary vertices and internal vertices of the triangular mesh are distinguished. For internal vertices, the maximum size filling particle circle is calculated and generated based on the constraints that the particle circles do not overlap and are tangent to the surrounding generated particle circles. The vertex filling particle circle information is saved. S6: Post-processing of the algorithm: Based on the information of all internal particle circles and vertex-filled particle circles, draw the particle stacking effect diagram, save all particle data, and draw particle size distribution curves and particle size distribution maps based on the particle data to complete the particle size distribution feature analysis.

2. The method for rapid stacking of discrete element particles according to claim 1, characterized in that, The random perturbation of node coordinates in step S2 is calculated as follows: Δpi = α·di·ri, where Δpi is the displacement vector of non-boundary grid node pi; α is the noise intensity coefficient; di is the local feature length of node pi, which is the average distance from the node to all adjacent nodes; and ri is a random direction vector, with each component of the vector uniformly distributed in the interval [−0.5, 0.5].

3. The method for rapid stacking of discrete element particles according to claim 1, characterized in that, After step S2 is completed, the following steps are also included: mesh quality assessment and visualization: visualize the triangular mesh before and after the disturbance, evaluate the mesh partitioning quality based on the degree to which the triangular mesh approaches an equilateral triangle, and output a mesh quality comparison chart.

4. The method for rapid stacking of discrete element particles according to claim 1, characterized in that, The specific process of obtaining triangular mesh topology information and generating traversal path in step S3 is as follows: 3.1 Represent mesh edges by a combination of two node numbers. Based on the node numbers contained in the triangular mesh unit, determine the triangular mesh unit to which each mesh edge belongs, and then determine the neighboring mesh units of each triangular mesh unit; 3.2 Select the geometric center mesh as the initial mesh for breadth-first search, and add the initial mesh number to the traversal sequence; 3.3 Search all neighboring mesh units of the current mesh unit in sequence, and append the effective neighboring mesh unit numbers to the traversal sequence in order; 3.4 Repeatedly perform neighbor search and number appending operations on the remaining mesh units in the traversal sequence until all triangular mesh units are included in the traversal sequence, and finally obtain the complete particle circle generation path.

5. The rapid stacking method for discrete element particles according to claim 1, characterized in that, In step S4, the triangular mesh is divided into boundary meshes and non-boundary meshes. Based on the particle generation status of adjacent meshes, five scenarios are used to generate internal particle circles: 4.1 Scenario 1: Boundary mesh, with two adjacent meshes that have already generated particle circles; the new particle circle must be tangent to the particle circles of both adjacent meshes and to one remaining edge of the current mesh; 4.2 Scenario 2: Boundary mesh, with only one adjacent mesh that has already generated particle circles; the new particle circle must be tangent to the particle circle of that adjacent mesh and to the two remaining edges of the current mesh; 4.3 Scenario 3: Non-boundary mesh, with particle circles generated in all three adjacent meshes; the new particle circle must be tangent to the particle circles of all three adjacent meshes; 4.4 Scenario 4: Non-boundary mesh, with two adjacent meshes that have already generated particle circles; the new particle circle must be tangent to the particle circles of both adjacent meshes and to one remaining edge of the current mesh; 4.5 Condition 5: Non-boundary mesh, and there is only one adjacent mesh with a generated particle circle; the new particle circle must be tangent to the particle circle of the adjacent mesh and tangent to the remaining two mesh edges of the current mesh; based on the constraints of the corresponding condition, the particle circles inside each triangular mesh are solved.

6. The method for rapid stacking of discrete element particles according to claim 1, characterized in that, The method for distinguishing between boundary vertices and internal vertices in step S5 is as follows: represent the grid edges in the form of node pairs with the smaller node number first and the larger node number second, and count the frequency of occurrence of each grid edge; the grid edge that appears only once is the boundary edge, and the boundary vertex is obtained from the node corresponding to the boundary edge, and the remaining nodes are the internal vertices.

7. The method for rapid stacking of discrete element particles according to claim 1, characterized in that, The specific steps for generating filled particle circles at the internal vertices in step S5 are as follows: 5.1 For a single internal vertex, extract all triangular mesh cells containing that vertex and filter out the generated particle circles with a radius greater than 0 in the corresponding mesh. 5.2 Traverse all three-circle combinations of the filtered particle circles, and solve for the candidate filling particle circles with the constraint that the new particle circle is tangent to the three particle circles in the combination and does not overlap with the particle circles already generated in the field; 7.3 Compare the radii of all candidate filling particle circles, and select the particle circle corresponding to the maximum radius as the final filling particle circle of the vertex.

Images