Method and apparatus for accelerated computation of discrete element simulation of particles of arbitrary shape

Abstract

The application discloses a discrete element simulation acceleration calculation method for arbitrary shape particles, comprising the following steps: establishing a triangular template grid for arbitrary shape particles; constructing a two-level hierarchical bounding box tree through the triangular template grid; emitting a light ray in the gradient direction of a discrete potential field function at the vertex of the triangular template grid, and performing collision detection with the hierarchical bounding box tree; and calculating inter-particle contact force and updating particle motion state through an energy conservation discrete potential contact model. By combining a discrete potential field function and a light ray tracing algorithm, a light ray tracing discrete element simulation algorithm for arbitrary shape particles is proposed, which can significantly reduce storage, improve calculation efficiency and robustness compared with a table lookup discrete element algorithm; and a GPU light ray tracing core and a general computing core can be organically combined for discrete element simulation acceleration of arbitrary shape particles, so that the calculation resources of the GPU can be fully utilized, and the simulation calculation efficiency is improved.

Description

Technical Field

[0001] This invention relates to the field of high-performance computing technology for industrial software, and in particular to a method and apparatus for accelerating discrete element simulation of particles of arbitrary shape. Background Technology

[0002] Particle shape plays a crucial role in determining the physical properties and mechanical behavior of granular materials. When fluid flow is present in a granular system, particle shape can significantly influence the movement of nearby fluids, such as bedload sediment transport. In recent years, the role of particle shape has received increasing attention, particularly in the field of numerical modeling. Although the Discrete Element Method (DEM) has been widely used to explore the microscopic roots of the mechanical behavior of granular materials, traditional sphere-based DEMs have limitations. Therefore, various variants of DEMs have been developed to simulate particles of different shapes. However, for large-scale or large-scale simulations, there is a need to introduce computationally efficient models of irregular or even realistically complex particle shapes. This would help reveal the physical and mechanical mechanisms of large granular systems and contribute to the design and optimization of engineering systems in civil, chemical, and mining fields.

[0003] Previous research has largely focused on convex shapes or concave shapes composed of convex elements, such as ellipsoids, hyperellipsoids, extended hyperellipsoids, convex / extended polyhedra, and convex non-uniform spline surfaces. These shapes allow for efficient inter-particle contact detection and computation using convex optimization algorithms, such as those based on the GJK and LM algorithms. However, these algorithms are not directly applicable to particles of arbitrary shapes. Generally, for particles of arbitrary shapes, regardless of their shape representation, exhaustive algorithms can be used. While these algorithms have good robustness, they often sacrifice computational efficiency.

[0004] When particles of arbitrary shape are represented by triangular meshes, there are three common methods for contact detection and resolution. The first method calculates the overlap volume between two meshes by dividing the particle into tetrahedra and resolving the tetrahedral overlaps one by one. This method has good numerical stability when using volume-based contact models, but incurs the highest computational cost. Therefore, the second method avoids calculating the overlap volume by calculating the intersection lines between two meshes. This significantly improves computational efficiency as it degenerates from calculating tetrahedral overlaps to calculating triangular intersection lines. The third method is based on point-to-surface algorithms commonly used in DEM variants, such as LS-DEM and SDF-DEM. In these methods, one particle in the contact pair is represented by a level set (LS) or a signed distance field (SDF) function, while the surface of the other particle is represented by discrete nodes. This method further simplifies the intersection problem between triangles into a vertex-triangle potential problem. The key step is to calculate the potential of a vertex inside a particle. A direct approach is to find a suitable vertex-triangle pair, such as the shortest-distance pair, which is called direct computation. However, the need for explicit comparison of distances between all vertices and triangles presents a significant challenge. To overcome this, LS or SDF functions are used to pre-cachise potential values ​​within a particle, transforming the contact problem into querying the potential value of a vertex within another particle. The pre-cached potential values ​​are typically stored as lattice points for fast lookup using lookup algorithms, allowing the potential value at a vertex to be interpolated based on the values ​​of nearby mesh nodes. However, the allocation and access to pre-cached data can lead to efficiency issues related to computer memory. This approach performs poorly in parallel computation on advanced computing systems such as modern GPUs. Specifically, on GPUs, lookup algorithms suffer from performance degradation due to data communication bottlenecks, and powerful GPU cores are underutilized due to limited computational resources.

[0005] Therefore, a new solution is needed. Summary of the Invention

[0006] The main objective of this invention is to provide a method and apparatus for accelerating discrete element simulation of particles of arbitrary shape.

[0007] To achieve the above objectives, this invention provides a method for accelerating discrete element simulation of particles of arbitrary shape, comprising the following steps:

[0008] Create a triangular template mesh for particles of arbitrary shape;

[0009] A two-level hierarchical bounding box tree is constructed using the triangular template mesh;

[0010] Rays are emitted from the vertices of the triangular template mesh in the gradient direction of the discrete potential function, and collision detection is performed with the hierarchical bounding box tree; and

[0011] The contact force between particles is calculated and the particle motion state is updated using an energy-conserving discrete potential contact model.

[0012] In the discrete element simulation acceleration method for particles of arbitrary shape provided by this invention, in the triangular template mesh, the particle surface x s Given by the following formula:

[0013] x s =sR′x m +x p

[0014] Where, x m It is a surface point on the template mesh; x p R is the centroid of the particle's position; R' is the 3x3 matrix along the particle's direction; s is the isotropic scaling factor.

[0015] In the discrete element simulation acceleration method for particles of arbitrary shape provided by the present invention, the two-level hierarchical bounding box tree includes a bottom boundary volume hierarchical structure tree and a top boundary volume hierarchical structure tree. The bottom boundary volume hierarchical structure tree includes axis-aligned bounding boxes of all triangular faces of a particle surface, and the top boundary volume hierarchical structure tree includes axis-aligned bounding boxes of all particles.

[0016] In the discrete element simulation acceleration calculation method for particles of arbitrary shape provided by this invention, the formula for the discrete potential field function is as follows:

[0017]

[0018]

[0019] Where k is the material stiffness; S i S is the area vector of vertex i; p Given the surface area of ​​the particle; Let θ be the vertex normal and θ be the potential order. It is the point of impact. Surface normal.

[0020] In the discrete element simulation acceleration calculation method for arbitrary shaped particles provided by the present invention, in the step of emitting rays from the vertices of the triangular template mesh in the gradient direction of the discrete potential field function and performing collision detection with the two-level hierarchical bounding box tree, the GPU ray tracing core is invoked, and the intersection of the ray and the axis-aligned bounding box is checked by the box intersection evaluator. If the ray intersects with the axis-aligned bounding box, the intersection point of the ray and the triangle face is calculated by the ray-triangle intersection evaluator and the impact distance dh is returned.

[0021] In the discrete element simulation accelerated calculation method for particles of arbitrary shape provided by this invention, in the step of calculating the contact force between particles using an energy-conserving discrete potential contact model, the discrete potential field function is used... and In the middle, the i-th vertex and its point of impact The potential value is used to calculate the normal contact force, vertex-potential normal contact force. Surface-potential normal contact force and combined normal contact force The formulas are as follows:

[0022]

[0023]

[0024]

[0025]

[0026]

[0027]

[0028] Where k is the normal contact stiffness; dh is the impact distance; γ i It is a local geometric measurement of the normalized surface area projection.

[0029] In the discrete element simulation accelerated calculation method for particles of arbitrary shape provided by this invention, in the step of calculating the contact force between particles using an energy-conserving discrete potential contact model...

[0030] Calculate the tangential contact force using the following formula.

[0031]

[0032] △f i =k t s

[0033]

[0034] Where μ is the coefficient of friction, Δf i k represents the increment of the tangential contact force. t It is the tangential contact stiffness. For the previous tangential contact force, n c It is the contact normal.

[0035] Furthermore, to achieve the above objectives, the present invention also provides a discrete element simulation acceleration computing device for particles of arbitrary shape, comprising:

[0036] The template creation module is used to create triangular template meshes for particles of arbitrary shapes;

[0037] A hierarchical bounding box tree construction module is used to construct a two-level hierarchical bounding box tree using the triangular template mesh;

[0038] A collision detection module is configured to emit rays at the vertices of the triangular template mesh in the gradient direction of a discrete potential function and perform collision detection with the hierarchical bounding box tree; and

[0039] The calculation module is used to calculate the contact force between particles and update the particle motion state using the energy-conserving discrete potential contact model.

[0040] Furthermore, to achieve the above objectives, the present invention also provides a discrete element simulation acceleration computing device for particles of arbitrary shape, comprising:

[0041] Memory, used to store computer programs;

[0042] A processor is configured to implement the steps of the discrete element simulation accelerated computation method for arbitrary-shaped particles as described above when executing the computer program.

[0043] Furthermore, to achieve the above objectives, the present invention also provides a computer-readable storage medium storing a computer program that, when executed by a processor, implements the steps of the discrete element simulation accelerated computation method for particles of arbitrary shape as described above.

[0044] In this invention, a ray tracing discrete element simulation algorithm for arbitrary-shaped particles is proposed by combining discrete potential field functions and ray tracing algorithms. Compared with table lookup-based discrete element algorithms, it can significantly reduce storage while improving computational efficiency and robustness. Furthermore, it can organically combine GPU ray tracing cores with ordinary computing cores to accelerate discrete element simulation of arbitrary-shaped particles, thereby making full use of GPU computing resources and improving simulation efficiency. Attached Figure Description

[0045] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on the provided drawings without creative effort:

[0046] Figure 1 The diagram shows a flowchart of the discrete element simulation acceleration calculation method for particles of arbitrary shape provided in Embodiment 1 of the present invention.

[0047] Figure 2 The diagram shown is a schematic of a two-stage acceleration structure.

[0048] Figure 3 The diagram shows (a) the continuous potential field function, (b) the discrete vertex potential field function, and (c) the discrete surface potential field function.

[0049] Figure 4 The following are several different light emission strategies;

[0050] Figure 5 The figures show (a) contact geometry measurements and (b) discrete vertex potential field functions;

[0051] Figure 6 The diagram shown is a simulation of a nut being turned into a bolt.

[0052] Figure 7 The template mesh used in the particle packing simulation is shown below;

[0053] Figure 8 The packing structure of 500 sand particles with an equivalent radius of 5.0 cm is shown;

[0054] Figure 9 and Figure 10 The stacking structure of 500 bananas (equivalent size 5.2 cm) and 500 pears (equivalent size 8.3 cm) is shown;

[0055] Figure 11 The image shows a mixed stack of 400 bananas and 100 pears using the default stacking simulation configuration;

[0056] Figure 12 and Figure 13 The image shows the accumulation of sand particles using a grid template 127v.

[0057] Figure 14 The figure shows a simulation of the packing of sand particles with extreme particle size ratios;

[0058] Figure 15 The image shows the template mesh for simulating a collapsed column and its 3D-printed sample.

[0059] Figure 16 The diagram shows the initial and final accumulation structures of sand particles.

[0060] Figure 17 The image shown is a snapshot of a sand particle column during the collapse process based on the vertex potential contact model.

[0061] Figure 18 The configuration of the extended hyperellipsoidal particles is shown at different time periods. Detailed Implementation

[0062] To facilitate understanding of the present invention, a more complete description will be given below with reference to the accompanying drawings. Typical embodiments of the invention are shown in the drawings. However, the invention can be implemented in many different forms and is not limited to the embodiments described herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete.

[0063] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains. The terminology used herein in the description of the invention is for the purpose of describing particular embodiments only and is not intended to be limiting of the invention.

[0064] Reference Figure 1 , Figure 1 This is a flowchart illustrating an embodiment of the discrete element simulation acceleration method for arbitrary-shaped particles according to the present invention. In one embodiment, the discrete element simulation acceleration method for arbitrary-shaped particles includes:

[0065] Step S10: Create a triangular template mesh for particles of arbitrary shape.

[0066] Specifically, in one embodiment of the present invention, in the triangular template mesh, the particle surface x s Given by the following formula:

[0067] x s =sR′x m +x p

[0068] Where, x m It is a surface point on the template mesh; x p R' is the centroid of the particle's position; R' is a 3x3 matrix along the particle's direction; s is the scaling factor.

[0069] Specifically, in this invention, the surface of particles, regardless of its complexity, can be effectively represented by a triangular mesh with a certain resolution. For simplicity, we use the number of vertices N. v The number of doughs Nf To quantify the resolution of the particle surface. A triangular face is defined by three vertices a, b, and c, and the outward normal of the face is given by the following formula:

[0070]

[0071] The area vector S of a vertex is defined by weighting the area vectors of its adjacent faces, and the formula is as follows:

[0072]

[0073] Where A i and These are the area and outward normal of the i-th adjacent face of a vertex, respectively; N af It is the number of adjacent faces. Another geometric quantity is the vertex normal n. v Defined as an angle-weighted pseudonormal, used to calculate the signed distance field of a polyhedron:

[0074]

[0075] Where α i It is the angle of the i-th face at the vertex; It is the outward normal of the i-th face; N af It is the number of faces adjacent to a vertex.

[0076] Mesh-based particle shape representation methods can flexibly simulate particles of arbitrary shapes, which can be generated by numerical algorithms, such as reconstructing spherical harmonic functions or X-ray CT scan images. Furthermore, when discrete element methods are coupled with computational fluid dynamics, discrete vertices can serve as Lagrangian points in immersive boundary techniques for particle-fluid interactions.

[0077] Furthermore, the mass m and the center point x p The moment of inertia J is a fundamental inertial property used to solve particle motion problems. These measurements can be calculated by decomposing the particle into tetrahedra, each with a triangular facet on the particle's surface as its base and a reference point (such as the origin) as its vertex. The particle's inertial properties can then be determined by summing the measurements of each tetrahedron, i.e.:

[0078]

[0079]

[0080]

[0081] Where ρ is density; a i b i and c iThese are the position vectors of the three vertices of the i-th triangular face (i.e., the base of the i-th tetrahedron); N f This represents the total number of triangular facets on the particle surface. The directed volume of the i-th tetrahedron is given by the following formula:

[0082]

[0083] Three vertices a i b i and c i The markings are counter-clockwise to ensure the face normals point outwards. The moment of inertia of the i-th tetrahedron is obtained relative to the center of mass, i.e.:

[0084]

[0085]

[0086] h = [a α b α c α ]

[0087] t = [2a β +b β +c β a β +2b β +c β a β +b β +2c β ]

[0088] Where α and β are exponents in [0, 1, 2]; ρ is the density. The principal moment of inertia can be obtained through J. (i) The eigenvalues ​​are obtained by eigendecomposition, where the eigenvectors (e1, e2, e3) and eigenvalues ​​(J1, J2, J3) are the principal axes and the corresponding principal moments of inertia, respectively.

[0089] Furthermore, to facilitate shape transformation and save computational memory, a template mesh is introduced. The mesh's center of mass is used as the local origin, and its principal axes are used as local Cartesian coordinate axes. Therefore, the mesh can be transformed to its local coordinate system in the following way:

[0090]

[0091] in and These are the new position vector and the original position vector of the j-th vertex in the mesh, respectively; M is the rotation matrix determined by the principal axes, i.e.:

[0092] M = |e1, e2, e3|

[0093] The sequence of (e1, e2, e3) in M ​​may be different, where the eigenvectors (e1, e2, e3) and eigenvalues ​​(J1, J2, J3) are the principal axes and the corresponding principal moments of inertia, respectively.

[0094] Furthermore, shape transformations, such as scaling, shearing, rotation, and translation, can be applied to the template mesh. This saves memory by sharing mesh data with particles that may have different positions, orientations, or shapes. Additionally, this method avoids the accumulation of particle shape errors during particle motion integration. Furthermore, the transformation matrices S, R, and T are used to illustrate scaling, rotation, and translation transformations, given by the following formulas:

[0095]

[0096] The overall transformation is given by the following matrix:

[0097]

[0098] Regarding the template mesh, the particle surface x s Given by the following formula:

[0099] x s =sR′x m +x p

[0100] Where x m It is a surface point on the template mesh; x p R is the particle's position (center of mass); R' is a 3x3 matrix in the particle's orientation; s is the scaling factor (i.e., Sx = Sy = Sz = s). With the scaling factor, the basic particle geometry / inertia metrics are given by the following:

[0101] m′=ms 3 J′=Js 5

[0102] The particle orientation is usually represented by a quaternion q(w,x,y,z) to facilitate the implementation of rotational integrals. Therefore, the above formula is rewritten as:

[0103] x s =sqx m +x p .

[0104] Step S20: Construct a two-level hierarchical bounding box tree using the triangular template mesh.

[0105] Specifically, in one embodiment of the invention, the most computationally intensive aspect of ray tracing is collision detection between rays and objects. To address this, a complex data structure is needed to organize spatial objects. This structure is often referred to as an accelerated structure in computer graphics. When dealing with objects with complex shapes, their surfaces are typically represented by a triangular mesh. Therefore, collision detection between rays and objects is simplified to checking for potential intersections between rays and the triangular mesh. The invention utilizes a workflow similar to that of general collision detection in ray tracing. Therefore, contact detection can be easily implemented using state-of-the-art ray tracing algorithms, including accelerated structures and hardware acceleration.

[0106] Specifically, the acceleration structure can be built as a bounding box hierarchy (BVH) tree. This invention introduces two levels (bottom and top) of sub-BVH trees to facilitate the management of the acceleration structure. Figure 2 For the bottom layer, the axis-aligned bounding boxes (AABBs) of all triangular faces on a particle surface are grouped into a BVH tree. For the top layer, the AABBs of all particles are grouped into another BVH tree. Using a mesh template, particles with the same mesh template can share a bottom-layer BVH tree but with different scaling, rotation, and translation. In other words, a bottom-layer BVH can be built relative to the mesh template, and for a particle with a given scaling, rotation, and translation, the corresponding BVH tree can be easily inserted into the top-layer BVH tree; this is called a particle instance (of the template mesh).

[0107] Assuming the particle shape remains constant during the simulation, the particle instance needs to be updated whenever its position or orientation is changed during particle dynamics. A two-stage BVH scheme allows for highly efficient updating of the entire acceleration structure based on previous data, particularly useful for discrete element problems that typically involve only small particle motions within a finite time step. It's worth noting that the acceleration structure's performance can degrade significantly after many updates; in such cases, a new acceleration structure can be built instead of updating based on previous data.

[0108] Step S30: Light rays are emitted from the vertices of the triangular template mesh in the gradient direction of the discrete potential function, and collision detection is performed with the hierarchical bounding box tree.

[0109] Specifically, in one embodiment of the present invention, the formula for the discrete potential field function is as follows:

[0110]

[0111]

[0112] Where k is the material stiffness; S iS is the area vector of vertex i; p Given the surface area of ​​the particle; Let θ be the vertex normal and θ be the potential order. It is the point of impact. Surface normal.

[0113] To introduce an energy-conserving contact model for numerical stability, the potential field function needs to be properly defined. In previous discrete element method (DEM) variants, such as level-basis DEM and symbolic distance-field DEM, the level set and symbolic distance field function can be considered as general potential field functions that are continuously defined in space.

[0114] Specifically, given a particle surface S, a potential field function W(S, x) can be defined. p ) and, make point x p The potential value is minimal relative to the surface. Figure 3 Typically, such continuous potential functions are implicit and require pre-caching on the spatial grid to facilitate lookup table operations, as implemented in LS-DEM and SDF-DEM. In this invention, we propose two discrete potential functions, one for the i-th vertex or face: and like Figure 3 As shown in b and c, discrete potential field functions define potential values ​​only at each individual vertex or face, compared to potential fields based on level set functions or signed distance field functions. This feature significantly reduces computational memory usage because it eliminates the need to distribute pre-cached potential values ​​across the particles.

[0115] For ray-surface intersection pairs, this invention will discretize the potential field function W(x) v (x) is defined as:

[0116]

[0117] Where k is the material stiffness; S i S is the area vector of vertex i; p Given the surface area of ​​the particle; Let θ be the vertex normal, θ be the potential order, and d be the given distance.

[0118]

[0119] in It is the point of impact. The surface normal is equal to the normal n of the surface being hit. f Therefore, the two discrete potential functions are rewritten as:

[0120]

[0121]

[0122] Their derivatives with respect to x are given by the following formula:

[0123]

[0124]

[0125] Furthermore, in this embodiment, the GPU ray tracing core is invoked, and the intersection of the ray and the axis-aligned bounding box is checked using the bounding box intersection evaluator. If the ray intersects the axis-aligned bounding box, the intersection point of the ray and the triangle face is calculated using the ray-triangle intersection evaluator, and the hit distance d is returned. h .

[0126] Ray tracing is a widely used technique in computer graphics and games to enhance the realism of rendered images and videos. In actual rendering, millions of rays are projected into a model scene, and when these rays collide with objects, their paths change through reflection or refraction, similar to the real world. Using a similar approach, ray tracing has been used for contact detection in computer graphics. In this invention, we propose a ray tracing-based algorithm for contact detection in discrete-element simulations of particles of arbitrary shapes.

[0127] like Figure 4 As shown, a line of length l i The ray will be projected onto a vertex i along the opposite direction of the normal, and the length of the ray can be determined as follows:

[0128]

[0129] in It is the point of impact on the particle surface, that is, the intersection of the ray and the triangular facet of the impact (see...). Figure 4 a); It is the position vector of vertex i. However, it is not necessary for the ray to hit its own particle surface. Therefore, the ray length is reduced by the weight w (see...). Figure 4 b). Furthermore, the ray length can be the same for all vertices of a particle (see...). Figure 4 c). To simplify calculations, the length of a single ray can be defined as the average of three half-lengths (lx, ly, lz), i.e.:

[0130]

[0131] The half-length can be defined relative to the axis-aligned bounding box (AABB) or the main axis of the particle.

[0132] Given two particles, contact occurs when a ray from one particle strikes a triangular facet of the other. To facilitate contact detection between rays and triangular faces, a bounding box hierarchy (BVH) is introduced to organize the axis-aligned bounding box (AABB) for each triangular facet. Therefore, ray-triangular facet contact detection first checks the intersection of the ray with the AABB; if the ray intersects the AABB, then a precise intersection calculation is performed on the ray-triangular facets.

[0133] The computational framework proposed in this invention is well-suited for parallel computing on traditional multi-core or multi-node CPU systems. Furthermore, it introduces a new paradigm of efficient parallel computing by making extensive use of modern GPU hardware. Unlike traditional GPGPU computing that only uses shader cores (such as NVIDIA's CUDA cores), this invention employs state-of-the-art hardware, namely ray tracing cores, for calculating contact detection and resolution between particles of arbitrary shapes.

[0134] This invention utilizes GPU acceleration, including general acceleration using shaders (such as CUDA acceleration on NVIDIA GPUs) and ray tracing (RT) acceleration using RT cores. Therefore, the method of this invention can be simply referred to as the Ray Tracing Discrete Element Method (RTDEM). In RTDEM, numerous rays emanate from the vertices of a particle and are then queried against the accelerated structure. While this detection task is typically computationally time-consuming, it can generally be accomplished through parallel computation on the shaders of a multi-core CPU or GPU. In addition to GPU shaders for general computation (such as CUDA cores on NVIDIA GPUs), RTDEM also significantly benefits from the ray tracing core, which was first introduced into NVIDIA's Turing architecture in 2019 to accelerate ray tracing tasks. The ray tracing core has two separate physical units: a ray-bounding box intersection estimator and a ray-triangle intersection estimator. The bounding box-ray intersection estimator checks for intersections between rays and bounding boxes, which can be seen as a coarse-screening stage for contact detection. The ray-triangle intersection evaluator is a fine-grained solution stage for contact detection. It calculates the intersection of the ray and the triangle face and returns the hit distance d. h .

[0135] Step S40: Calculate the contact force between particles and update the particle motion state using the energy-conserving discrete potential contact model.

[0136] Specifically, in one embodiment of the present invention, the contact force includes normal contact force, tangential contact force, and motion integration.

[0137] a. Normal contact force

[0138] Given a scalar potential field W(x), the normal contact force f nIt can be obtained in the following ways:

[0139]

[0140] This leads to an energy-conserving contact model. In previous studies, node-to-surface algorithms were used to calculate contact forces. Specifically, nodes (or vertices) on a particle surface are considered intruders.

[0141] like Figure 5 As shown, it is necessary to calculate the potential value of the intruder relative to the surface of another particle. The traditional method is to do this by looking up a pre-cached (pre-calculated and stored) table, which requires pre-caching the potential value and its gradient. This invention eliminates the need for the traditional pre-cached table, which is a major innovation of this invention.

[0142] With our proposed discrete potential function and The i-th vertex and its point of impact The potential value is used to calculate the contact force. Therefore, at the vertex and hitting point The normal contact force is given by the following formula:

[0143]

[0144]

[0145]

[0146] Where k is the normal contact stiffness; dh is the impact distance; γ i It is a local geometric measurement of the normalized area projection and can be regarded as the weight of the entire normal contact force. Therefore, the normal contact force can be defined using the following three equations:

[0147]

[0148]

[0149]

[0150] The third equation is for simplification purposes. and The contact forces are a linear combination. For clarity, we represent these three normal contact forces as the vertex potential normal contact forces. normal contact force of surface potential and combined normal contact force The contact force acts directly on particle B, while its reaction force acts on particle A. The direction of the normal contact force is considered to be the contact normal n. cAnd the contact point It can be defined as:

[0151]

[0152] Even when θ = 2, the overall normal contact force often does not exhibit a linear relationship with the impact distance. For γ i and In the extreme case where all vertices remain consistent, such as in a plane, d h With the same embedding depth as traditional models, the contact model is equivalent to a volume-based linear contact model. The volume-based linear contact model also holds true for other symmetrical shapes such as spheres.

[0153] b. Tangential contact force

[0154] Regarding the tangential contact force, we combine a linear contact model with a Coulomb sliding model. The calculation of the tangential contact force follows the traditional implementation method. At the contact point x... c The relative displacement s is defined as:

[0155] s = (v AB -(v AB ·n c )n c )Δt

[0156] Where, n c It is the contact normal, equal to -n v ;v AB Δt is the relative velocity at the moment of contact; Δt is the time step. To maintain consistency, the particle with ray projection is denoted as particle A, and the particle with surface impact is denoted as particle B. Therefore, the relative velocity v AB Given by the following formula:

[0157]

[0158] Where v, ω and x p These represent the tangential velocity, rotational velocity, and particle position, respectively; subscripts A and B denote the corresponding particles A and B, respectively. Therefore, the tangential contact force Δf i The increment is given by the linear displacement-force law, that is:

[0159] △f i =k t s

[0160] Where kt is the tangential contact stiffness. Unlike the calculation of the normal contact force, the previous tangential contact force needs to be rotated to the current tangential plane. Therefore, two orthogonal axes n, orthogonal to the contact normal, are introduced. oa and torsion axis n ta Rotation of a convenient force is defined as:

[0161] n oa =n′ c ×n c

[0162]

[0163] Where n c` It is the previous contact normal. Therefore, the previous tangential contact force... for:

[0164]

[0165]

[0166] In addition, the previous tangential contact force It can be projected onto the current contact plane, that is:

[0167]

[0168] Therefore, tangential contact force The result obtained is:

[0169]

[0170] Its size is affected by the Coulomb condition, namely:

[0171]

[0172] Where μ is the coefficient of friction. Additionally, the tangential contact force... It is not necessarily parallel to the relative tangential velocity when in contact in three-dimensional space.

[0173] Therefore, the tangential contact force is calculated using the following formula.

[0174]

[0175] △f i =k t s

[0176]

[0177] Where μ is the coefficient of friction, Δf i k represents the increment of the tangential contact force. t It is the tangential contact stiffness. For the previous tangential contact force, n c It is the contact normal.

[0178] c. Motor integration

[0179] The motion of particles follows the Newton-Euler equations, which are listed here for ease of description:

[0180]

[0181]

[0182] Where v and ω are the translational and rotational velocities of the particle, respectively; m and J are the mass and moment of inertia of the particle, respectively; F (b) It is physical strength; F (d) and T (d) These are damping force and torque, respectively; F (c) and T (c) These are the contact force and torque on the particles, respectively, given by the following formulas:

[0183]

[0184]

[0185] in, and These are the number of vertices within the contact overlap area of ​​particles A and B, respectively; N P x is the number of neighbors of particle A; A It is the position vector of particle A.

[0186] The present invention also provides a discrete element simulation acceleration computing device for particles of arbitrary shape, comprising:

[0187] The template creation module is used to create triangular template meshes for particles of arbitrary shapes;

[0188] A hierarchical bounding box tree construction module is used to construct a two-level hierarchical bounding box tree using the triangular template mesh;

[0189] A collision detection module is configured to emit rays at the vertices of the triangular template mesh in the gradient direction of a discrete potential function and perform collision detection with the hierarchical bounding box tree; and

[0190] The calculation module is used to calculate the contact force between particles and update the particle motion state using the energy-conserving discrete potential contact model.

[0191] Those skilled in the art will understand that the above are embodiments of the discrete element simulation acceleration computing device for arbitrary-shaped particles provided by the present invention. This system and device belong to the same inventive concept as the above-described discrete element simulation acceleration computing method for arbitrary-shaped particles. For details not described in detail in the embodiments of the discrete element simulation acceleration computing device for arbitrary-shaped particles, please refer to the embodiments of the above-described discrete element simulation acceleration computing method for arbitrary-shaped particles.

[0192] This invention also provides an accelerated computing device for discrete element simulation of particles of arbitrary shape, which may include:

[0193] Memory, used to store computer programs;

[0194] When a processor executes a computer program stored in the aforementioned memory, it can perform the following steps:

[0195] A triangular template mesh for particles of arbitrary shape is established; a two-level hierarchical bounding box tree is constructed using the triangular template mesh; light rays are emitted from the vertices of the triangular template mesh in the gradient direction of the discrete potential field function and collision detection is performed with the hierarchical bounding box tree; and the contact force between particles is calculated and the particle motion state is updated using an energy-conserving discrete potential contact model.

[0196] This invention also provides a computer-readable storage medium storing a computer program, which, when executed by a processor, can perform the following steps;

[0197] A triangular template mesh for particles of arbitrary shape is established; a two-level hierarchical bounding box tree is constructed using the triangular template mesh; light rays are emitted from the vertices of the triangular template mesh in the gradient direction of the discrete potential field function and collision detection is performed with the hierarchical bounding box tree; and the contact force between particles is calculated and the particle motion state is updated using an energy-conserving discrete potential contact model.

[0198] The computer-readable storage medium may include various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0199] Numerous specific details are set forth in the specification provided herein. However, it will be understood that embodiments of the invention may be practiced without these specific details. In some instances, well-known methods, structures, and techniques have not been shown in detail so as not to obscure the understanding of this specification.

[0200] The following are several examples of simulations using the algorithm provided in this invention.

[0201] 1. Nut to Bolt Simulation

[0202] This test involves the frictionless rotation of a nut onto a bolt under gravity. The rotation of the nut involves complex contact detection within the nut-bolt system, which has been selected as a benchmark for multibody dynamic simulations in previous studies. The energy conservation equation for the nut is as follows:

[0203]

[0204] The constraint equation is dh / dt = v z and 2π / ω z =H / v z Analyze the rotational speed ω z Given by the following formula

[0205]

[0206] Where H is the vertical displacement of the nut around the bolt by a rotation of 2π. The simulation is as follows... Figure 6 As shown, the bolt is fixed, and the nut rotates freely under its own weight. The nut's mass m and principal inertia I... z They were 38.95 kg and 0.4 kg·m, respectively. 2 The vertical span between two adjacent bolt tracks H is 15.07 mm, resulting in an analytical rotational acceleration of 2.137 rad / s². 2 Three potential-based contact force models are employed: vertex potential energy, surface potential energy, and combined potential energy model, along with two θ( Used in conjunction with 2). In the simulation, all contacts are frictionless. To ensure numerical stability, a time step of 1×10⁻⁶ is used. -4 At time s, the normal contact stiffness kn is set to 1×10 6 N / m 1 / 2 and 1×10 7 N / m, With θ = 2, a small local damping coefficient of 0.025 is used. Aside from numerical stability, the value of kn does not affect the results because it is simplified during the calculation of the analytical solution.

[0207] like Figure 6 As shown, the surface potential contact model outperforms the other two models, almost matching the analytical solution results. The combined potential contact model performs better than the vertex potential model in most cases. Furthermore, the results for only the nut generating rays are almost identical to those for both the nut and bolt generating rays. The average simulation speeds are 15813 steps / second and 9724 steps / second for only the nut generating rays and for both the nut and bolt generating rays, respectively. This indicates that having only one contact particle generating rays can significantly improve computational performance, and the results are comparable to those for both contact particles generating rays. However, for… When using the vertex potential or combined potential contact model, the simulation results for generating the two types of rays show some differences. Furthermore, choosing which particle to use for ray generation in a multi-particle system remains an issue that requires further investigation. The remaining simulations in this embodiment will use the projected rays from both contacting particles. In summary, regardless of the ray projection, the model order parameter θ and the surface-potential contact model are the most stable.

[0208] 2. Particle packing simulation

[0209] This section will examine the robustness and performance of RTDEM in simulating particle systems with different shapes. We used sand grain template meshes (127v and 502v) to study the effect of mesh resolution. Two other template meshes, namely banana and pear, are also used. Figure 7 As shown.

[0210] 2.1 Stacking Method

[0211] In the particle packing simulation, unless otherwise specified, the following simulation strategy is used for particle packing. The initial packing configuration is a 5×5×20 grid with a spacing of 0.11 meters. Then, the particles are allowed to freely pack under gravity into a square box with a base of 1.42×1.42m. For particle-particle contact, the normal contact stiffness kn is set to 1×10⁻⁶. 6 N / m 1 / 2 The tangential contact stiffness ks is set to 1×10 5 For the contact between the particle and the box, the normal contact stiffness kn is set to 1×10 N / m. 8 N / m 1 / 2 The tangential contact stiffness ks is set to 1×10 5 The frictional force between particles and between particles and the box is set to 0.3 N / m. The material density of the particles is 2650 kg / m³. 3 A local damping coefficient of 0.3 is used to promote energy dissipation. The accumulation simulation runs for 30,000 steps with a fixed time step of 1 × 10⁻⁶. -4 Second.

[0212] 2.2 Simulation of Sand Particle Accumulation at Different Mesh Resolutions

[0213] Figure 8 The stacked structure of 500 sand particles with an equivalent radius of 5.0 cm is shown. The total number of vertices exceeds 63,500 and 251,000 respectively in the stacks using 127v and 502v meshes. Although there are modest differences in the final structures between the two mesh resolutions, the stacking simulation is numerically stable even with a mesh of only 127 vertices.

[0214] 2.3 Simulation of the stacking of bananas and pears

[0215] Figure 9 and Figure 10 The stacking structures of 500 bananas (equivalent size 5.2 cm) and 500 pears (equivalent size 8.3 cm) are shown respectively. The banana and pear stacks have a total of over 126,000 and 375,000 vertices respectively. The particle stacking simulation is also numerically stable for these two distorted shapes.

[0216] Figure 11 The image shows a mixed stack of 400 bananas and 100 pears using the default stacking simulation configuration.

[0217] 2.4 Simulation of the packing of sand particles with different gradations

[0218] To further test the robustness of RTDEM in simulating polydisperse particle systems with different particle sizes, two additional examples of sand particle packing using a 127v mesh template are presented: one is a particle packing with a uniform particle size distribution between 2.5 and 5.0 cm (referred to as GSD packing), and the other is a binary packing with sizes of 1.0 and 5.0 cm (referred to as binary mixed packing). For the binary mixed packing, the normal contact stiffness kn is set to 1 × 10⁻⁶. 5 N / m 1 / 2 and 1×10 4 N / m 1 / 2 This corresponds to the contact between the particle box and the contact between the fine particles. In addition, there are 100 coarse particles and 9375 fine particles. The rest of the simulation settings are the same as the previous example. Figure 12 and 13 The stacking configurations for GSD stacking and binary hybrid stacking are shown at different time points. It can be observed that the simulated values ​​are stable during the stacking process.

[0219] 2.5 Summary of computational efficiency

[0220] Table 1 shows the time taken for particle packing simulation using RTDEM (NVIDIA RTX 2080Ti GPU).

[0221]

[0222] Computational efficiency is another important metric for evaluating RTDEM performance. Table 1 summarizes the time consumption of all simulations, where the simulation scale can be roughly quantified by the total number of vertices for intuitive comparison. Particle shape and size distribution are also major factors affecting computational efficiency. From the table, we can conclude that in most simulations, the surface potential contact model is slightly less computationally efficient than the vertex potential contact model. This is because the surface potential contact model requires additional computation to access the surface normals. Considering overall performance, including accuracy, robustness, and efficiency, we recommend using the surface potential contact model.

[0223] 2.6 Simulation of Packing of Sand Particles with Extreme Particle Size Ratios

[0224] As described above, the proposed computational framework can efficiently simulate the dynamics of monodisperse and polydisperse particle systems with arbitrary shapes. An extreme case is provided here to demonstrate the performance of the proposed framework. In the simulation, 200,000 fine particles are deposited on top of coarse particles within a box. The equivalent sizes of the coarse and fine particles are 1 meter and 1 centimeter, respectively, implying an extreme size ratio of 100:1. All particles have the same shape as the template mesh sand 127v. The surface potential contact model has normal contact stiffnesses of 1×10⁻⁶ for coarse and fine particles, respectively. 8 N / m 1 / 2 and 1×10 4 N / m 1 / 2 Other simulation parameters were the same as in previous simulations. The simulation ran for 30,000 steps with a fixed time step of 1×10⁻⁶. -4 A snapshot of the particle structure in seconds, such as Figure 14 On a gaming GPU (NVIDIA RTX 2080 Ti GPU with 11GB of RAM), the entire process took 29.4 minutes, with GPU memory usage of 2323MB. This demonstrates that the proposed framework provides an efficient solution for modeling indeterminate-shaped particles, even those with large aspect ratios.

[0225] 3. Collapse Column Simulation

[0226] To verify the simulation results in RTDEM, an experimental test of column collapse was conducted. Two sets of particle columns were prepared using monodisperse particles: one set consisted of 1540 sand particles introduced in previous sections, and the other set consisted of extended hyperellipsoidal particles. The extended hyperellipsoidal particles have smooth surfaces. The particles were 3D printed using a Formlabs Form 3L printer with a resolution of 50 micrometers. Figure 15 The template mesh and its 3D printed sample are shown.

[0227] 3D-printed particles were poured into an acrylic container with a base measuring 4 cm × 3.5 cm. After removing one sidewall of the container, the particle column collapsed to near equilibrium. Both the initial and final structures are drawn on [the provided text]. Figure 16 middle.

[0228] It can be seen that, due to their smoother surfaces, the extended hyperellipsoidal particles travel further than the sand particles. The slopes of the sand particles and the extended hyperellipsoidal particles are 27.7° and 21.9°, respectively. The final profile was used as the baseline for the numerical simulation. Additionally, a moderate deviation is inevitable between the replicated profiles.

[0229] In the proposed RTDEM numerical model, model parameters such as contact stiffness, friction coefficient, and damping coefficient may affect the final profile after the collapse of the particle column. To avoid any possible interference, a comprehensive parameter study is not presented here. Furthermore, detailed calibration of these model parameters is required to quantitatively reproduce the experimental results, but this is beyond the scope of this embodiment. For simplicity, the friction coefficient μ is used as a variable for sensitivity analysis, while the values ​​of other model parameters are chosen as follows: mass density 1200 kg / m³. 3 For particle-particle contact, the normal and tangential contact stiffnesses are 5 × 10⁻⁶. 2 N / m 1 / 2 and 1×10 4 N / m. For particle-box contact, the normal and tangential contact stiffness are 5 × 10⁻⁶ N / m. 3 N / m 1 / 2 and 1×10 5 N / m, damping coefficient is 0.1.

[0230] Figure 17 The vertex-potential contact model shows a snapshot of the sand particle column during collapse. Using two friction coefficients, μ = 0.3 and μ = 0.2, the resulting slopes in the final state are 29.4° and 26.3°, respectively. The simulated final profile with μ = 0.2 closely approximates the experimental profile.

[0231] Figure 18 The configuration of the extended hyperellipsoidal particles at different time points is shown. With friction coefficients μ = 0.2 and μ = 0.15, the final slopes are 25.2° and 21.0°, respectively. Clearly, by slightly increasing the friction coefficient, the simulated configuration can be made to match the experimental results. This also demonstrates that RTDEM can reproduce experimental results without requiring extensive adjustments to model parameters.

Claims

1. A method for accelerating discrete element simulation of particles of arbitrary shape, characterized in that, Includes the following steps: Establish Triangular template mesh for particles of arbitrary shape; A two-level hierarchical bounding box tree is constructed using the triangular template mesh; Rays are emitted from the vertices of the triangular template mesh in the gradient direction of the discrete potential function, and collision detection is performed with the hierarchical bounding box tree; and The contact force between particles is calculated and the particle motion state is updated using an energy-conserving discrete potential contact model. The formula for the discrete potential field function is as follows: Where k is the material stiffness; For the i-th vertex The area vector; Given the surface area of ​​the particle; Let θ be the vertex normal and θ be the potential order. It is the point where light hits the particle surface. Surface normal.

2. The method for accelerating discrete element simulation of particles of arbitrary shape as described in claim 1, characterized in that, In the triangular template mesh, the particle surface Given by the following formula: in, These are surface points on the template mesh; It is the position of the particle's centroid; It is a 3x3 matrix along the particle direction; s is the scaling factor.

3. The method for accelerating discrete element simulation of particles of arbitrary shape as described in claim 2, characterized in that, The two-level hierarchical bounding box tree includes a bottom bounding box hierarchical structure tree and a top bounding box hierarchical structure tree. The bottom bounding box hierarchical structure tree includes axis-aligned bounding boxes of all triangular faces on a particle surface, and the top bounding box hierarchical structure tree includes axis-aligned bounding boxes of all particles.

4. The method for accelerating discrete element simulation of particles of arbitrary shape as described in claim 3, characterized in that, In the step of emitting rays from the vertices of the triangular template mesh in the gradient direction of the discrete potential function and performing collision detection with the two-level hierarchical bounding box tree, the GPU ray tracing core is invoked. The bounding box intersection evaluator checks the intersection of the ray with the axis-aligned bounding box. If the ray intersects with the axis-aligned bounding box, the ray-triangle intersection evaluator calculates the intersection point of the ray and the triangle face and returns the collision distance. .

5. The method for accelerating discrete element simulation of particles of arbitrary shape as described in claim 4, characterized in that, In the step of calculating the interparticle contact force using the discrete potential contact model based on energy conservation, the discrete potential field function is used... and In the middle, the i-th vertex and its point of impact The potential value is used to calculate the normal contact force; the vertex potential normal contact force. Surface potential normal contact force and combined normal contact force The formulas are as follows: Where k is the normal contact stiffness; It is the distance between the vertex of the emitted ray and the point of impact; It is a local geometric measurement of normalized area projection.

6. The method for accelerating discrete element simulation of particles of arbitrary shape as described in claim 5, characterized in that, In the step of calculating the interparticle contact force using the discrete potential contact model based on energy conservation, the tangential contact force is calculated using the following formula. : Where μ is the coefficient of friction. For the increment of tangential contact force, It is the tangential contact stiffness. For the previous tangential contact force, It is the contact normal.

7. A discrete element simulation acceleration computing device for particles of arbitrary shape, characterized in that, include: The shape template creation module is used to create triangular template meshes for particles of arbitrary shapes; A hierarchical bounding box tree construction module is used to construct a two-level hierarchical bounding box tree using the triangular template mesh; A collision detection module is used to emit light rays at the vertices of the triangular template mesh in the gradient direction of the discrete potential function and perform collision detection with the hierarchical bounding box tree; as well as The calculation module is used to calculate the contact force between particles and update the particle motion state using the energy-conserving discrete potential contact model. The formula for the discrete potential field function is as follows: Where k is the material stiffness; For the i-th vertex The area vector; Given the surface area of ​​the particle; Let θ be the vertex normal and θ be the potential order. It is the point where light hits the particle surface. Surface normal.

8. A discrete element simulation acceleration computing device for particles of arbitrary shape, characterized in that, include: Memory, used to store computer programs; A processor, configured to implement the steps of the discrete element simulation accelerated computation method for arbitrary shaped particles as described in any one of claims 1 to 6 when executing the computer program.

9. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores a computer program that, when executed by a processor, implements the steps of the discrete element simulation accelerated computation method for particles of arbitrary shape as described in any one of claims 1 to 6.

Images