JP5467262B2 - Particle simulation apparatus and particle simulation method - Google Patents

Description

  The present invention relates to a particle simulation apparatus and a particle simulation method.

  DEM (Discrete Element Method), which originated from the civil engineering research field, was later applied to powder engineering, and nowadays it is very diverse that it is difficult to handle not only as a particulate material but also as a continuum. As a numerical analysis method for this problem, the demand is increasing. On the other hand, many techniques other than DEM such as MD (Molecular Dynamics), LGA (Lattice Gas Automaton), and SPH (Smoothed Particle Hydrodynamics) have been developed for numerically tracking discrete particle motion. However, among the many numerical particle analysis methods, the reason why DEM is widely applied is to solve for Newton's equation of motion by calculating the interaction force between particles while taking into account the geometric information of each particle such as the size and shape of the particle. The calculation can be performed very quickly. Of course, DDA (Discontinuous Deformation Analysis) other than DEM is a method that can take into account the size and shape. However, since the calculation load is originally high, it is difficult to handle a large number of particles practically.

  By the way, the actual granular materials and powders we see like sand and wheat flour (hereinafter referred to as “powder particles” together) are not subject to special treatment such as particle size adjustment by sieving. In general, particles of various sizes are present in a mixed state (hereinafter referred to as a “particle size distribution” state). Actually, this particle size distribution is a major factor governing the mechanical properties such as homogeneity, deformation, and flow of the granular material. Homogenization and mixing efficiency, transport and filling control are always a problem. Therefore, it is indispensable to correctly understand the influence of the particle size distribution on the behavior of the entire particle group in order to handle powder particles.

Since DEM is a technique that can take into account the size of the particles as described above, it is natural to want to investigate the micromechanics brought about by the particle size distribution with DEM. However, when the particle diameter becomes 1 / n, the volume of one particle becomes (1 / n) ³ , and n ³ particles are required to fill the volume of one particle of the original size. (Strictly speaking, since voids are generated, the order is less than n ³ but the order does not change). Therefore, in reality, a very large particle size distribution cannot be handled due to the amount of memory of the computer and the calculation time. For example, in natural soil, the particle size is distributed over about five digits from gravel of several centimeters to clay particles of several microns. Then, when all sizes occupy the same volume, even if the number of particles having the maximum particle size is one, the total number of particles is on the order of (10 ⁵ ) ³ .

In addition, when dealing with the particle size distribution by DEM, in addition to the above-mentioned problem of the number of particles, an unbalance of the number of contact points due to the difference in particle size also becomes a problem. For example, when the ratio of the maximum particle size to the minimum particle size is r: 1, the case where the surface of the maximum particle is filled with the minimum particle and the case where the surface of the minimum particle is filled with the maximum particle are considered. the number of particles are respectively ^{4 (r + 1) 2/} 1 2 and ^{4 (r + 1) 2 /} r 2. That is, the imbalance of the number of contact points is on the order of r ² : 1 at the maximum. The unbalance of the number of contact points of each particle leads to an unbalance between the number of contact force calculations for each particle and the total calculation amount of contact, which is a big barrier to parallel calculation. In other words, unlike the case of uniform particle size, in DEM, when the particle size ratio increases, the number of particles is expected to increase simply from the geometric information of the particle size, and the number of contact points is unbalanced. As a result, it can be seen that the calculation load is considerably high. In addition to this, even in the processing inside the actual DEM program, if there is a particle size distribution, the search for the contact candidate to be contacted and the processing for storing the contact candidate pair are performed with a uniform particle size. As compared with the case of, the calculation load is significantly increased.

  Several DEM programming techniques have been proposed for speeding up processing related to extraction and storage of contact candidate pairs for the problem of this particle size distribution (for example, Patent Document 1).

  On the other hand, recently, speeding up of numerical simulation using a GPU (Graphics Processing Unit), which has a much larger ratio of calculation performance to price than a CPU, has been actively performed. Along with this, for particle system simulation, a faster algorithm for GPUs has been proposed that can perform massively parallel computation equivalent to vector computation and is cheaper than vectorization processing using an expensive supercomputer (for example, patents). Reference 2).

JP 2007-286514 A JP 2009-69930 A

  However, the conventional DEM acceleration method such as Patent Document 1 is a method of first optimizing the cell size and cell search range after storing the particles in the cell, and further increases the large-scale acceleration by parallelization or vectorization. Not conscious. As a specific example, if the process of storing particles in a cell with a particle number loop is performed by SIMD (Single Instruction Multiple Data) parallel processing such as vectorization, two or more particles are stored in the same cell. This is very inconvenient because it causes memory write contention. This memory write contention problem occurs when multiple particles are contained in one cell even in the case of a single particle size component, but when there is a particle size distribution, there is a large distribution in the number of particles contained in one cell itself. Therefore, when the program speed is increased by parallelization or vectorization, it becomes an even more troublesome problem.

  Further, in the GPU acceleration algorithm such as Patent Document 2, in the detection process of the contact candidate particle pair, a method of searching for a cell after storing particles in the cell is used in the same manner as a CPU-based program. However, as described above, even when a plurality of particles are stored in the same cell in the GPU, memory write contention still occurs, and several measures have been proposed in order to avoid this problem. One is a method of repeating particle storage and confirmation operations a plurality of times, and the other is a method using an exclusive processing function (atomic function) attached to the compiler. In addition, in order to prevent memory write competition in the interparticle contact force calculation routine, a device that avoids the problem by calculating the contact force between the particles ij twice for each of the particles i and j. Has been made.

  However, as described above, if there is a particle size distribution, the number of particles stored by the cell tends to be greatly biased, and the number of contact particles for large particles is much larger than that for small particles. It is easy to damage greatly. For example, when a GPU calculation is performed on a binary particle group using a three-dimensional DEM using these methods, the calculation time is 3 as compared with a single component when the particle size ratio of the two components is only 1: 4. It has been reported that it will be nearly double.

  And when performing a DEM simulation of a granular material in consideration of the particle size distribution, the number of neighboring particles that should consider the contact force increases rapidly compared to the case of the same number of particles with a uniform particle size. There has been a problem that the amount of memory used that must be secured to realize the DEM calculation increases.

  The GPU calculation cannot be expanded in hardware like the CPU calculation, and can only be processed within a limited memory capacity mounted on the graphic card. Of course, a method of relying on the CPU memory for the shortage of the GPU memory is also conceivable. However, since the transfer speed between the GPU memory and the CPU memory is very slow, the high calculation efficiency of the GPU is not utilized at all and it is not recommended. As described above, when the GPU is used for the DEM calculation, since the memory capacity of the current GPU is about 4 GB at the maximum, the problem of the memory usage is a limitation in applying the GPU to the particle simulation. . This problem is not limited to the GPU, but also hardware (for example, a vector processor) that performs the SIMD (Single Instruction Multiple Data) parallel processing of the shared memory type similar to the GPU (hereinafter referred to as “shared memory type parallel computer”). This is also a problem.

  The present invention has been made to solve the above-described problems, and suppresses the memory usage when executing a DEM (Discrete Element Method) calculation of a granular material having a particle size distribution using a shared memory parallel computer. An object of the present invention is to provide a particle simulation apparatus and a particle simulation method that can be performed.

In order to solve the above problems, the particle simulation apparatus according to the present invention calculates the position / velocity based on the contact force with other particles for a particle group having a particle size distribution in the work space, and simulates the behavior of the particles. A particle simulation apparatus for inputting particle information including position information and velocity information for each particle of a particle group, particle information holding means for holding particle information, and each of the particles as position information. A granting means for granting specific information to be specified in an order according to the selected particle, in order of ascending order of the specific information, from the particles, and a particle from the target particle that may be in contact with the target particle a contact candidate selecting means for selecting a contact candidate pair with small other particles diameters for each contact candidate pairs selected by contacting a candidate selecting means, the contact The contact force between the two particles of the complementary pair is calculated based on the position information and velocity information of these particles, and the calculated contact force information is determined according to the order of the contact candidate pairs corresponding to the contact force. A reference for referring to contact force calculation means stored in a table and contact force of each particle with a particle larger than the particle diameter of other particles located in the vicinity of the particle from the contact force table. Contact force reference table creating means for creating a contact force reference table including information, and for each of the particles, the number of contact candidates indicating the number of particles smaller than the particle diameter of other particles located in the vicinity of the particle The cumulative contact candidate number calculating means for calculating the cumulative contact candidate number, which is the cumulative value of each, and the contact force between each particle and a particle having a particle diameter larger than the particle is referred to in the contact force reference table. Extract from the contact force table, extract the contact force with a particle having a particle diameter smaller than the particle by specifying the storage position of the contact force table using the accumulated number of contact candidates, and use these extracted contact forces A contact force sum calculating means for calculating the total contact force for each particle; a calculating means for calculating new position information and velocity information of the particles based on the contact force for each particle calculated by the contact force sum calculating means; Particle information updating means for updating the particle information held in the particle information holding means using the new position information and velocity information calculated by the calculating means.

  Similarly, in order to solve the above problem, the particle simulation method according to the present invention calculates the position / velocity based on the contact force with other particles for a particle group having a particle size distribution in the work space, and A particle simulation method in a particle simulation apparatus for simulating behavior, wherein particle information including position information and velocity information for each particle in a particle group is input, and the particle information is held in a particle information holding unit; For each of the particles, an assigning step for assigning specific information to be specified in an order corresponding to the position information, and sequentially selecting the target particle from the particles in ascending order of the specific information, and the possibility of being in contact with the target particle A contact candidate selection step for selecting a contact candidate pair with another particle having a particle diameter smaller than that of the target particle, and a contact candidate selection step For each of the contact candidate pairs selected in the above, the contact force between the two particles of the contact candidate pair is calculated based on the position information and velocity information of these particles, and the calculated contact force information is calculated for the contact force pair. A contact force calculation step for storing in the contact force table according to the order of the contact candidate pairs corresponding to the force, and for each of the particles, among the other particles located in the vicinity of the particle, a particle larger than the particle diameter of the particle A contact force reference table creating step for creating a contact force reference table including reference information for referring to the contact force from the contact force table, and for each of the particles, among the other particles located in the vicinity of the particle An integrated contact candidate number calculating step of calculating an integrated contact candidate number that is an integrated value of the number of contact candidates indicating the number of particles smaller than the particle diameter, and for each particle, The contact force with a particle having a larger particle diameter is extracted from the contact force table using the reference information in the contact force reference table, and the contact force with a particle having a particle diameter smaller than that particle is contacted using the cumulative number of contact candidates. Identify and extract the storage position of the force table, and use these extracted contact forces to calculate the sum of contact forces for each particle, and contact for each particle calculated in the contact force summation step Based on the force, a calculation step for calculating new position information and velocity information of the particles, and the particle information held in the particle information holding means are updated using the new position information and velocity information calculated in the calculation step. And a particle information updating step.

  According to such a particle simulation apparatus and particle simulation method, for each of the particles in the work space, the particle diameter of the particle among the other particles that are located in the vicinity of the particle and are likely to come into contact with the particle For smaller particles, only the accumulated contact candidate number indicating the number of particles is retained, and the storage position of the contact force table is specified using this accumulated contact candidate number to contact the particle with other particles. Extract power. Thus, the contact force reference table does not need to include information on particles smaller than the particle diameter of the particles among other particles located in the vicinity of the particles, and only information on particles larger than the particle diameter of the particles. If you hold it. When a particle group having a particle size distribution is targeted for simulation, the number of particles that can be contacted at the same time increases when considering contact to particles with a particle size smaller than the self, so it should be reserved in advance for the contact force reference table You have to increase the memory area. On the other hand, according to the above configuration of the present invention, it is only necessary to consider contact with particles having a particle size larger than that of the self, and the number of particles considering contact can be reduced. The memory area to be allocated can be reduced, and as a result, the memory usage can be suppressed. By suppressing the memory usage in this way, it is possible to execute a particle simulation using a shared memory parallel computer such as a GPU.

Further, in the particle simulation apparatus, the work space is divided into a plurality of cells, each cell is assigned a cell number, and the particle information is specified based on the position information, and the cell in which the particle is arranged. The cell number is included, and the assigning means identifies each of the particles in the order according to the cell number, and further identifies the specific information identified in the order according to the particle diameter between the particles arranged in the same cell. to grant.

Similarly, in the particle simulation method, the work space is divided into a plurality of cells, each cell is assigned a cell number, and the particle information is a cell in which the particle is specified, which is specified based on position information. Cell number, and the assigning step specified for each of the particles in the order according to the cell number, and further in the order according to the particle diameter between the particles arranged in the same cell to grant specific information.

  With this configuration, each particle is specified in the order corresponding to the cell number, and the specific information specified in the order corresponding to the particle diameter is given between the particles arranged in the same cell. Sorting is performed in the order of numbers, and further in the order of particle diameter in the same cell. For this reason, when creating the contact force reference table, it is possible to easily recognize particles having a small particle diameter based on the specific information, and to determine whether or not particles having a small particle diameter can be included in the contact force reference table. There is no need to do it. As a result, calculation efficiency can be improved in the creation of the contact force reference table.

Further, in the particle simulation apparatus, the contact candidate selection unit selects the target particle from the particle group in ascending order of the specific information, and may be in contact with the target particle. The particle diameter is smaller than the target particle. Create a contact candidate list that includes contact candidate pair information consisting of pairs with other particles and contact candidate numbers assigned in ascending order to this set for all particles in the work space. For each piece of contact candidate pair information included in the candidate list, a contact force between two particles related to the contact candidate pair information is calculated based on the position information and velocity information of the particle, and the calculated contact force and contact candidate are calculated. The contact force information associated with the number is stored in the contact force table, and the contact force reference table creating means is more likely to be in contact with the particle for each particle. As reference information for referring to the contact force with other particles having a large size from the contact force table, the contact force that holds the other particle and the contact candidate number related to the particle in association with the other particle and the particle. A reference table is created, and the cumulative contact candidate number calculation means calculates the number of contact candidates indicating the number of other particles having a particle diameter smaller than the particle, which may be in contact with the particle, for each of the particles. The accumulated contact candidate number s_jgi [i] (i indicates the specific information of the particle), which is a value accumulated in the order of specific information, is calculated, and the contact force sum calculation means refers to the contact force reference table for each particle. The contact force is extracted from the contact force table using information, and the contact force with contact candidate numbers s_jgi [i-1] to s_jgi [i] -1 is calculated using the accumulated contact candidate number s_jgi [i]. Extract from the table and use these extracted contact forces It summation contact force for each particle.
In the particle simulation method, in the contact candidate selection step, the target particle is selected from the particle group in ascending order of the specific information, and the particle diameter may be smaller than the target particle that may be in contact with the target particle. Create a contact candidate list that includes contact candidate pair information consisting of pairs with other particles and contact candidate numbers assigned in ascending order to this set for all particles in the work space. For each piece of contact candidate pair information included in the candidate list, a contact force between two particles related to the contact candidate pair information is calculated based on the position information and velocity information of the particle, and the calculated contact force and contact candidate are calculated. The contact force information associated with the number is stored in the contact force table, and each of the particles is in contact with the particle in the contact force reference table creation step. As reference information for referring to the contact force with other particles having a particle diameter larger than that of the particles from the contact force table, the contact candidate numbers for the other particles and the particles are used as the reference information. And a contact force reference table to be held in association with the particle, and in the cumulative contact candidate number calculation step, each of the particles may be in contact with the particle, and may have another particle size smaller than the particle. An integrated contact candidate number s_jgi [i] (i indicates the specific information of the particle), which is a value obtained by integrating the number of contact candidates indicating the number of particles in the order of the specific information, For each of these, the contact force is extracted from the contact force table using the reference information of the contact force reference table, and the contact candidate numbers s_jgi [i-1] to s_jgi [ i] -1 contact force Extracting from the table, and using these extracted contact forces, the total contact force for each particle is calculated.

  According to the particle simulation apparatus and the particle simulation method according to the present invention, when the DEM (discrete element method) calculation of the granular material having the particle size distribution is executed by the shared memory parallel computer, the memory usage is suppressed. Can do.

It is a functional block diagram of the particle simulation device concerning one embodiment of the present invention. It is the schematic which shows the Voigt model applied to a particle | grain model in this embodiment. It is a figure which shows the work area | region in this embodiment. It is a figure which shows the flow of a particle number change. It is a figure which shows the process which changes the particle number in a work area, and the process which memorize | stores the maximum value and minimum value of a particle number in each cell. It is a figure which shows an example of a structure of a near particle table. It is a figure which shows an example of a structure of a contact candidate list | wrist. It is a figure which shows the process which investigates whether it was a pair one step before about the pair particle | grains of the present contact candidate list | wrist. It is a figure which shows an example of a structure of a contact force reference table. It is a figure which shows the correspondence of a contact force reference table, the number of accumulation contact candidates, and a contact force box. It is a flowchart which shows the process performed in the particle | grain simulation apparatus of this embodiment. It is a figure which shows the memory usage-amount of this embodiment at the time of changing the ratio of the largest particle diameter of a particle | grain model, and the minimum particle diameter, and a comparative example.

  Hereinafter, embodiments of the present invention will be described in detail with reference to the drawings. In the description of the drawings, the same or equivalent elements are denoted by the same reference numerals, and redundant description is omitted.

  FIG. 1 is a functional block diagram of a particle simulation apparatus 10 according to an embodiment of the present invention. As shown in FIG. 1, the particle simulation apparatus 10 includes a particle information holding unit (input unit, particle information holding unit) 11, a particle information acquisition unit 12, a contact candidate list update determination unit 13, and a particle number changing unit (granting unit). 14, a proximity particle table creation unit 15, a contact candidate list creation unit (contact candidate particle selection unit, integrated contact candidate number calculation unit) 16, a contact force reference table creation unit (contact force reference table creation unit) 17, a contact determination unit ( A contact force calculating means) 18, a contact force calculating section (contact force sum calculating means) 19, and a particle information updating section (calculating means, particle information updating means) 20 are provided.

  The particle simulation apparatus 10 physically includes a GPU (Graphics Processing Unit), a RAM (Random Access Memory) and a ROM (Read Only Memory) that are main storage devices, an auxiliary storage device such as a hard disk device, an input key that is an input device, and the like. And a communication module as a data transmission / reception device such as a network card. Each function of the particle simulation apparatus 10 to be described below with reference to FIG. 1 includes an input device and an output device under the control of the GPU by reading predetermined computer software on hardware such as a GPU and a RAM. This is realized by operating the communication module and reading and writing data in the RAM and the auxiliary storage device.

  Here, the particle model applied in this embodiment will be described first. The particle simulation apparatus 10 of the present embodiment is intended for DEM simulation of a granular material group having a particle size distribution. That is, a particle model in which particles having various particle sizes are mixed is applied.

ここでｘとｖは相対変位および相対速度ベクトル、Ｋは線形バネ７２における弾性係数、ηはダッシュポット７３の粘性減衰係数であり、それぞれの添字ｎ，ｔにより法線方向成分及び接線方向成分を表す。μ_tは滑り摩擦係数、ｍｉｎは絶対値が小さい方の値を表す関数である。粘性減衰係数ηは、反発係数eと次式の関連がある。

ここでｍ_ｉ，ｍ_ｊは、それぞれ粒子ｉ，ｊの質量を表す。添字ｑは、法線方向成分の場合はｑ＝ｎとなり、接線方向成分の場合にはｑ＝ｔとなる。 In such a particle model, the normal direction and tangential direction components F _n and F _t of the inter-particle contact force generated between the particles are expressed by the following equation using the Voigt model 71 shown in FIG. The

Here, x and v are relative displacement and relative velocity vectors, K is an elastic coefficient in the linear spring 72, η is a viscous damping coefficient of the dashpot 73, and the normal direction component and the tangential direction component are expressed by the subscripts n and t, respectively. Represent. mu _t is sliding friction coefficient, min is a function representing the value of the direction the absolute value is smaller. The viscous damping coefficient η is related to the restitution coefficient e by the following equation.

Here, m _i and m _j represent the masses of the particles i and j, respectively. The subscript q is q = n for the normal direction component and q = t for the tangential direction component.

ここで、ｕ_ｐとω_ｐは粒子の速度および角速度、ｍ_ｐとＩ_ｐは粒子の質量および慣性モーメント、Ｒは粒子中心から接触点へ向かう位置ベクトルであり、添字ｐは、粒子ｉの場合はｐ＝iとなり、粒子ｊの場合はｐ＝ｊとなる。また、（５）式の右辺第二項は転がり摩擦力を表し、μ_ｒは転がり摩擦係数、ｂは接触面幅である。（１）、（２）式のＶｏｉｇｔモデル７１は、後述する接触判定部１８で粒子間接触力を算出するのに用いられ、また（４）、（５）式の運動方程式は、粒子情報更新部２０で粒子の座標を算出するのに用いられる。 Therefore, the equation of motion of translation and rotation around the center of gravity handled by the DEM is expressed by the following equation.

Here, u _p and ω _p are the velocity and angular velocity of the particle, m _p and I _p are the mass and moment of inertia of the particle, R is a position vector from the particle center to the contact point, and the subscript p is for the particle i P = i, and for particle j, p = j. Also, (5) represents the frictional force right side second term rolling type, mu _r is the rolling friction coefficient, b is a contact surface width. The Voigt model 71 of the equations (1) and (2) is used to calculate the interparticle contact force by the contact determination unit 18 described later, and the equations of motion of the equations (4) and (5) The unit 20 is used to calculate the coordinates of the particles.

In addition, as shown in FIG. 3, the work area 51, which is an area in which particles move in the present embodiment, is divided into cubic cells whose side size is csize. The cell size csize is given by the following equation using the parameter φ (0 ≦ φ ≦ 1).
csize = 2.0 {(1.0−φ) r _max + φr _min )} (1.0 + α) (6)
Here, r _max is the maximum particle size of the particle model arranged in the work area, r _min is the minimum particle size, and α is a parameter for adjusting the update frequency, for example α = 0.2. Each cell in the work area 51 is assigned a cell number cell_id. The cell number is one-dimensional as shown in FIG. 3, and is represented by the following equation, for example.
cell_id = i.x + id.x x i.y + id.x x id.y xiz (7)
im = {0, 1, 2,…, id.m-1} (m = x, y, z)
Here, cell_id is a cell number, im is a cell position in the m direction, and id.m is the number of cells in the m direction. That is, the particle simulation apparatus 10 automatically calculates and sets the cell size csize of the work area 51 when the maximum particle diameter r _max , the minimum particle diameter r _min , and the parameters φ and α are set as simulation conditions. In addition, cell numbers can be assigned to further defined cells.

  The particle information holding unit 11 holds particle information for each of the plurality of particles in the work space 51. The particle information includes coordinates (position information), translation speed (speed information), rotation speed (speed information), particle radius, mass, particle number, and cell number. The coordinates are the three-dimensional coordinates (pos.x, pos.y, pos.z) in the three-dimensional work space 51 shown in FIG. The particle number is a number for identifying each particle. The cell number is the cell number of the cell to which the particle belongs. When the simulation is started, the particle information holding unit 11 sets the particle radius and mass of each particle according to the particle size distribution applied to the simulation, and randomly sets the coordinates, velocity, and particle number to each particle. assign. The cell number is assigned with the number of the cell to which the particle belongs, which is specified based on the coordinates randomly assigned to each particle. Each piece of information generated in this way is input and held as particle information in the initial state. Further, during the execution of the simulation, when each information of the particle information is updated by the particle information updating unit 20 described later, the particle information holding unit 11 updates the held particle information by the particle information updating unit 20. Replace with a new one.

  The particle information acquisition unit 12 acquires particle information from the particle information holding unit 11. The particle information acquisition unit 12 acquires the initial particle information set in advance from the particle information holding unit 11 at the start of the particle simulation, and is updated each time the particle information holding unit 11 is updated during the simulation operation. Accurate particle information is acquired from the particle information holding unit 11. The particle information acquisition unit 12 transmits the acquired particle information to the contact candidate list update determination unit 13.

  The contact candidate list update determination unit 13 determines whether or not to update a contact candidate list 53 described later. For example, the update determination is performed when the particle radius is r and there is a particle having an accumulated movement distance of the particle larger than rα. Here, α is a parameter for adjusting the update frequency. For example, α = 0.2. Further, the accumulated movement distance of the particles is, for example, every time the particle information acquisition unit 12 acquires new particle information from the particle information holding unit 11, based on the difference of the coordinate data with the particle information of one step before this one step. The movement distance of the particles at is calculated, and this is calculated for each particle. The contact candidate list update determination unit 13 stores the particle information of the previous step, the parameter α, and the cumulative movement distance of the particles. When new particle information is received from the particle information acquisition unit 12, the cumulative movement distance Is updated, and it is determined whether or not there is a particle having the updated accumulated movement distance larger than rα. Then, when there are particles whose accumulated movement distance is equal to or greater than a predetermined value, it is determined that the contact candidate list is to be updated, and new particle information is transmitted to the particle number changing unit 14. The integrated value of the particle movement distance is initialized every time the contact candidate list 53 is updated. If there is no particle whose updated accumulated movement distance is larger than rα, it is determined that the contact candidate list is not updated, and particle information is not transmitted to the particle number changing unit 14 and the contact determining unit 18 is newly updated. The correct particle information.

  The particle number changing unit 14 specifies each of the particles in the work area 51 in the order corresponding to the cell number (that is, specified in the order corresponding to the position information), and further between the particles arranged in the same cell. Thus, the particle information (specific information) specified in the order corresponding to the particle diameter is given. In other words, the particle number changing unit 14 sorts the particle numbers by changing the particle number of each particle based on the cell to which each particle in the work area 51 belongs and its particle diameter.

  When the particle number changing unit 14 receives the particle information from the contact candidate list update determining unit 13, as shown in FIG. 4, the cell number of the cell to which the particle belongs and the particle number of the particle belong to the particle number of each particle. The diameters are associated with each other to form one information set, and these information sets are stored in order of particle number (FIG. 4A). Next, these information sets are rearranged so that the cell numbers are in ascending order (FIG. 4B), and the particle numbers are reassigned in the order of the cell numbers (FIG. 4C). Subsequently, the particle diameters of the particle groups having the same cell number are compared, and the particle numbers are rearranged in descending order of the particle diameters in these particle groups (FIG. 4D), and the particle numbers are renumbered in the order of the particle diameters (FIG. 4). (E)).

  By the process of changing the particle number, for example, the particle number of the particles arranged in the work space 51 as shown in FIG. 5A is aligned with the cell number and the particle diameter as shown in FIG. 5B. Changed to particle number. The particle number changing unit 14 transmits the particle information in which the change of the particle number is reflected to the neighboring particle table creating unit 15.

  The neighboring particle table creating unit 15 creates a neighboring particle table 52 including information on neighboring particles located in the vicinity of each particle in the work space 51. The neighboring particle table creating unit 15 is a particle having a particle diameter larger than itself, among particles existing in a cell to which the target particle belongs and a cell adjacent to this cell or a neighboring cell further outside the adjacent cell, or a particle When the diameters are the same, particles having a particle number larger than that of the initial arrangement and having an interparticle distance of a certain value or less are referred to as “neighboring particles”. Details will be described below.

  The nearby particle table creation unit 15 reserves a memory area for creating the nearby particle table 52 in advance. The memory area to be secured is the same size for each neighboring particle table of each particle in consideration of the number of neighboring particles that may come into contact with each particle.

  Upon receiving the particle information from the particle number changing unit 14, the neighboring particle table creating unit 15 stores the minimum value Pnum_in_cellmin and the maximum value Pnum_in_cellmax among the particles having the same cell number in the cell. (FIG. 5C). In FIG. 5C, the maximum particle number is stored in the upper left of each cell, and the minimum particle number is stored in the lower right. That is, for each cell cell_id, only an array for storing information on the minimum value Pnum_in_cellmin [cell_id] and the maximum value Pnum_in_cellmax [cell_id] of the number of particles included in this cell is retained. As a result, it can be seen that the particles from the Pnum_in_cellmin [cell_id] number to the Pnum_in_cellmax [cell_id] number belong to the cell of each cell number cell_id in the descending order of the particle diameter.

The number nei_cell_id of the cell adjacent to the cell with the number cell_id is expressed by the following equation.
nei_cell_id [n] = cell_id + con [n] (8)
Here, con is a constant determined by id.x and id.y. Since the particle numbers are sorted by cell number, it can be seen that particles of numbers from Pnum_in_cellmin [nei_cell_id] to Pnum_in_cellmax [nei_cell_id] exist continuously in each adjacent cell.

Next, the neighboring particle table creation unit 15 prepares a box (array) nei_pn [i] [box_id] for storing neighboring particles as shown in FIG. 6 for each particle, and searches around the particle i ( Among the particles included in the definition (to be described later), (a) the interparticle distance with the particle i is smaller than dis_lim,
(B) larger than the diameter of the particle i,
(C) When the diameter is the same as the diameter of the particle i, the particles satisfying the condition that the particle number (i = p_id) at the initial placement is larger than the particle i are packed from the box with the smaller number (box_id) and stored. Here, dis_lim = (r _i + r _j ) (1.0 + α), and r _i and r _j are the particle radii of the neighboring particle pair i and j. α is a parameter for adjusting the number of contact candidate updates. For example, α = 0.2. r _i and r _j are included in the particle information received from the particle number changing unit 14, α is a value previously held by the neighboring particle table creating unit 15, and the neighboring particle table creating unit 15 stores these data. To calculate dis_lim.

  Particles satisfying the above conditions (a) to (c) are referred to as “neighboring particles” for the particle i. A set of a series of boxes shown in FIG. 6 prepared for all the particles i in the work space 51 is referred to as a neighborhood particle table 52 in this embodiment. The neighboring particle table 52 is a two-dimensional array nei_pn [i] [box_id] regarding the particle i and the number box_id of neighboring particles of the particle i.

  By this method, only the particle j having a particle diameter equal to or larger than the particle i is stored in the neighboring particle table 52 for the target particle i, and thus the increase in the memory usage necessary for the neighboring particle table 52 can be suppressed. it can.

Here, in order to determine the neighboring particles of the particle i, the neighboring particle table creating unit 15 needs to search for a particle at a distance of (r _max + r _i ) (1.0 + α) from the particle i. In addition to the 26 cells adjacent to the cell to which the particle i belongs, the search range is further expanded to the outer cell according to the radius r _max of the maximum particle of the particle model.

The range of neighboring cells searched by the neighboring particle table creation unit 15 includes the cell size csize calculated by the above equation (6), the maximum particle radius r _max indicating the particle diameter of the largest particle in the work area, and the focused particle radius. It is calculated _| required by following Formula from ri.
nei_cell = 2 · ceil {2.0 (r _max + r _i ) (1.0 + α) / csize} −1 (9)
ceil is a function for rounding up the fractional part, and is held in the neighboring particle table creation unit 15 in advance. nei_cell represents the number of cells searched in the + -direction for each of the x, y, and z directions from the cell to which the particle i belongs.

As a result, there is a concern that the number of cells to be searched increases rapidly as the cell size decreases, and the calculation load increases. However, for the cells outside the adjacent 26 cells, it is possible to improve the efficiency of searching for neighboring particles by using the arrangement order of the particles sorted by the particle diameter. In other words, particles from Pnum_in_cellmin [cell id] number to Pnum_in_cellmax [cell id] number belong to each cell in descending order of particle diameter, and when checking whether the particles in each cell are neighboring particles It will be examined in order of particle size. At this time, since the particle diameter that does not satisfy the above condition (a) is obtained from the positional relationship between the cell to which the particle i belongs and its neighboring cells, it is not necessary to investigate all the interparticle distances up to Pnum_in_cellmax [cell id] number. The search can be stopped when the particle number has a particle diameter that does not satisfy the condition (a). Here, the particle diameter rj that does not satisfy the condition (a) for the particle i satisfies the following expression.
2 ceil {2.0 (r _j + r _i ) (1.0 + α) / csize} −1 <dis cell (10)
dis cell represents the number of cells from the cell to which the particle i belongs to the cell to which the particle j belongs, and represents the smallest number of cells in the x, y, and z directions.

  Furthermore, the neighborhood particle table creation unit 15 calculates the number of neighborhood particles n_jli [i] for the particle i and the number of particles n_jgi [i] satisfying the above condition (a) but not satisfying (b) and (c), respectively. Store in a one-dimensional array.

In the example of FIG. 6, for example, as a result of searching the cell for the 4th particle (i = 4), the 5th, 6th, 7th, and 8th particles satisfy the condition (a). If the particles of Nos. 8 and 5 and 6 do not satisfy the conditions (b) and (c),
nei_pn [4] [0] = 7
nei_pn [4] [1] = 8
n_jgi [4] = 2
n_jli [4] = 2
It becomes.

  As described above, when receiving the particle information from the particle number changing unit 14, the neighboring particle table creating unit 15 demarcates a cell of the search range nei_cell for each particle of interest i, and other particles included in the search range. Among them, the particle numbers of the particles satisfying the above conditions (a) to (c) are stored in the array nei_pn [i] [box_id] of the neighboring particle table 52. Further, the number of neighboring particles is stored in the array n_jli [i], and the number of particles satisfying the above condition (a) but not satisfying (b) and (c), that is, the particle diameter of the particles existing in the search range is larger than the particle of interest. The number of small particles that are not stored in the neighboring particle table 52 is stored in the array n_jgi [i]. Then, information on these nei_pn [i] [box_id], n_jli [i], and n_jgi [i] is transmitted to the contact candidate list creating unit 16 together with the particle information.

  In the present embodiment, the neighboring particle table 52 holds only particles having a diameter equal to or larger than the particle size of the target particle as neighboring particles. In a multi-component system composed of a group of particles having a plurality of particle sizes, there is a possibility that the number of small particles that are candidates for contact with large particles becomes enormous. That is, the number of particles that can be simultaneously contacted with particles having a particle diameter smaller than that of the target particle is extremely large. For this reason, if the neighboring particles include particles having a particle diameter smaller than the target particle, the number of arrays prepared for the neighboring particle table in advance must be increased in consideration of the maximum value. In that case, the amount of memory necessary for the arrangement of the neighboring particle table 52 becomes enormous, and there is a risk that the global memory of the GPU is crushed. In this embodiment, since only particles larger than the self are stored in the nearby particle table 52, the memory usage of the array prepared for the nearby particle table 52 can be suppressed.

  The contact candidate list creation unit 16 sequentially selects the target particles from all the particles in the work space 51 in ascending order of the particle numbers (specific information), and from the target particles that may be in contact with the target particles. A contact candidate pair with another particle having a small particle size is selected. In other words, the contact candidate list creation unit 16 creates a contact candidate list 53 including particle pairs that may be in contact with all the particles in the work space 51. Specifically, the contact candidate list 53 is a set of one-dimensional arrays list_i [Lnum] and list_j [Lnum] that store the particle numbers of two particles i and j that may be in contact with each other. Here, Lnum is a contact candidate number assigned to each contact particle pair.

When the contact candidate list creation unit 16 receives information (neighbor particle table 52nei_pn [i] [box_id], number of neighboring particles n_jli [i], n_jgi [i], and particle information) from the neighborhood particle table creation unit 15, first, As shown in FIG. 7, the prefix sum s_jgi of n_jgi is obtained from the following equation.

  This prefix sum s_jgi [i] is present in the number of particles not stored in the nearby particle table 52 of the particle i, that is, in a predetermined search range around the particle i (neighboring cells defined by the above equation (9)). The number n_jgi [i] (number of contact candidates) of the particles smaller than the particle diameter of the particles i is a one-dimensional array representing the integrated value up to the particle number i. " In other words, the cumulative contact candidate number s_jgi [i] is a value obtained by integrating the contact candidate number n_jgi [i] in the order of the particle number i.

Next, the contact candidate list creation unit 16 sequentially selects the target particle i in ascending order of the particle number, and for each particle i, selects a cell in the search range around the particle i defined by the above equation (9). Searching again, among the particles included in the search range around the particle i, a) the distance between the particles i is smaller than dis_lim, b) smaller than the diameter of the particles i,
c) When the diameter is the same as the diameter of the particle i, the particle j satisfying the condition that the particle number at the initial placement is smaller than the particle i is extracted. Then, contact candidate numbers Lnum are given to these pairs of particles i and j in ascending order, and the particle numbers of the respective particles are stored in the one-dimensional arrays list_i and list_j to create a contact candidate list 53.

  Satisfying these conditions a) to c) is equivalent to satisfying the condition (a) for extracting the neighboring particles stored in the above-mentioned neighboring particle table and not satisfying (b) and (c). That is, a particle pair that is not stored in the neighboring particle table 52 (for example, a pair with i = 5 and 6 particles in the case of i = 4 of interest particles) is a contact candidate pair.

つまり、接触候補リスト５３は、作業空間５１の全粒子について、接触している可能性がある粒子ペアlist_i，list_jに番号Lnumを付加した接触候補ペア情報を含むものである。 For example, in the example of FIG. 7, the contents of the contact candidate list 53 are as shown in Table 1 below.

That is, the contact candidate list 53 includes contact candidate pair information obtained by adding the number Lnum to the particle pairs list_i and list_j that may be in contact with all particles in the work space 51.

  Since the amount of spring expansion and contraction required when obtaining the contact force in the tangential direction is given as an integrated value, if the contact is made in the previous step and contact is made in the current step, the tangential spring is updated when the contact candidate pair list is updated. It is necessary to take over the amount of expansion and contraction. Then, next, the contact candidate pair Blist_i, Blist_j one step before is searched, and it is checked whether the current contact candidate pair list_i, list_j was a pair one step before. If it is a pair, the normal vector in the particle model shown in FIG. 2 and the force (spring expansion / contraction amount) in the tangential direction are stored in a new candidate list number array.

Here, the notation of each parameter and array one step before and now is as follows.

  When there is a contact candidate pair with particle numbers i = list_i [Lnum] and j = list_j [Lnum], the particle numbers of the paired particles one step before are Bi = Bpn [i] and Bj = Bpn [j]. Contact candidates Blist_j [BLnum] in the range from Bs_jgi [Bi−1] to Bs_jgi [Bi] −1 are checked to see if they match Bj. If they match, the list vector BLnum is used to substitute the normal vector one step before and the tangential contact force BSpring [BLnum] into the current list number array Spring [Lnum].

  With reference to FIG. 8, this process will be described with reference to FIG. 8 for checking whether a pair particle of i = 4 and j = 6 in the current contact candidate list number Lnum = 2 is a pair one step before.

  First, the particle number one step before i = 4 and j = 6 is obtained as Bi = 2 = Bpn [4] and Bj = 3 = Bpn [6] using the Bpn [] function. Here, the relationship between Bi and Bj always satisfies the conditions a) to c) of the contact candidate list 53, like the relationship between i and j. This is because the conditions a) to c) use the particle diameter that does not change by sorting and the particle number at the initial placement. Therefore, the contact candidate particle numbers for Bi are stored in Blist_j [BLnum] where the list number BLnum is in the range of Bs_jgi [Bi−1] to Bs_jgi [Bi] −1 (1 to 1 in this example). Therefore, if Blist_j [BLnum] matches Bj within this BLnum range, Blist_j [1] = 3 and Bj = 3 when BLnum = 1.

  Therefore, it can be seen that the candidate pair with the current list number Lnum = 2 was a pair one step before, and the list number is BLnum = 1. By substituting the normal vector of BLnum = 1 and the contact force value in the tangential direction into the array of Lnum = 2 (Spring [2] = BSpring [1]), it is possible to take over the previous contact state. .

  In this manner, the contact candidate list creation unit 16 receives various information reports (neighbor particle table 52nei_pn [i] [box_id], neighboring particle numbers n_jli [i], n_jgi [i], and particle information from the neighboring particle table creation unit 15). ) For each particle of interest i, a number Lnum is assigned to a pair with the particle j that satisfies the conditions a) to c) among the other particles included in the search range, and the particle pair i , J are stored in arrays list_i [Lnum] and list_j [Lnum] of the contact candidate list 53, respectively. Further, the cumulative contact candidate number s_jgi [i], which is an integrated value of the number n_jgi [i] of particles that are present in the peripheral cell of the target particle i and are smaller than the particle diameter of the particle i, is calculated. Further, from the current contact candidate pairs list_i [Lnum], list_j [Lnum], a pair that was a contact candidate one step before is extracted, and the contact force of these pairs one step before is taken over from BSpring [BLnum], Store in Spring [Lnum]. Then, the derived contact candidate list 53list_i [Lnum], list_j [Lnum] and the neighboring particle table 52nei_pn [i] [box_id] are transmitted to the contact force reference table creation unit 17, and the contact candidate list 53list_i [Lnum], list_j [Lnum], the accumulated contact candidate number s_jgi [i], and information on the contact force Spring [Lnum] of the contact candidate pair inherited from one step before are transmitted to the contact determination unit 18 together with the particle information.

  The contact force reference table creation unit 17 creates a function Ref [i] [box_id] for calculating the corresponding candidate list number Lnum from box_id. In the present embodiment, the two-dimensional array Ref [i] [box_id] representing this function is referred to as a contact force reference table 54. As shown in FIG. 9, the contact force reference table 54 refers to the contact force Fij [Lnum] received by the i particles in the contact candidate list 53 from the j particles for each arbitrary particle i (i = 4 in FIG. 9). (Lnum) (reference information). In particular, in the present embodiment, the contact force reference table 54 includes a reference destination (Lnum) of contact force with a neighboring particle stored in the neighboring particle table 52 relating to the target particle i, that is, a particle having a larger particle diameter than the target particle i. Only stored. Here, the reference destination (Lnum) of the contact force refers to two types of one-dimensional arrays Fijt [Lnum] and Fijr [Lnum] that hold a translational component and a rotational component of the interparticle contact force described later, and will be described later (13). This is a contact candidate pair list number (contact candidate number) Lnum used to acquire the corresponding contact force from the contact force box composed of the contact force Fij [Lnum] composed of β [Lnum] used in the equation.

  The contact force reference table creating unit 17 reserves a memory area for creating the contact force reference table 54 in advance. The memory area to be secured is the same size for each contact force reference table 54 of each particle in consideration of the number of neighboring particles that may come into contact with each particle.

  Specifically, when the contact force reference table creation unit 17 receives the contact candidate list 53list_i [Lnum], list_j [Lnum] and the neighboring particle table 52nei_pn [i] [box_id] from the contact candidate list creation unit 16, the contact candidate For the particle i (= list_i [Lnum]) and the particle j (= list_j [Lnum]) which are the Lnum-th particle pair in the list 53, the neighborhood where the particle i is stored by searching the neighborhood particle table 52 of the particle j The particle table number k (= box_id) is acquired. As described above, because the particle j of the contact candidate pair has a smaller particle diameter than the particle i, the particle i is stored in the neighboring particle table 52 of the particle j as a neighboring particle of the particle j. Then, the contact candidate list number Lnum of the particles i and j is stored in the k-th array Ref [j] [k] of the contact force reference table 54 for the particle j.

In addition, the creation process of the contact candidate list 53 by the contact candidate list creation unit 16 and the creation process of the contact force reference table 54 by the contact force reference table creation unit 17 can be performed by the following algorithm. “Contact candidate pair condition == true” in the following algorithm means that the conditions a) to c) of the contact candidate list 53 are satisfied.

  In the present embodiment, the contact force reference table 54 holds only particles having a diameter equal to or larger than the particle size of the target particle. For this reason, the memory usage of the arrangement prepared for the contact force reference table 54 can be suppressed for the same reason as the above-mentioned neighboring particle table 52.

  The contact force reference table creation unit 17 transmits the derived contact force reference table 54Ref [i] [box_id] to the contact force calculation unit 19.

  The contact determination unit 18 refers to the contact candidate list 53 to determine contact between two particles for all paired particles, and determines contact forces Fijr and Fijt between the particles when determined to be in contact. calculate.

Specifically, the contact determination unit 18 includes the contact candidate list 53list_i [Lnum], list_j [Lnum], the accumulated contact candidate number s_jgi [i] from the contact candidate list creation unit 16, and the contact candidate pair inherited from the previous step. When the contact force Spring [Lnum] and the particle information are received or the particle information is received from the contact candidate list update unit 13, the contact candidate list 53 is scanned to determine the contact of the particles i and j of the contact candidate pair. Do. In the contact determination, for example, the interparticle distance between the particle i and the particle j is calculated, and it is determined that the two particles are in contact when the interparticle distance is equal to or less than the sum of the particle diameters of both particles, r _i + r _j. To do.

In the contact determination, when it is determined that both particles are in contact, the contact forces Fn and Ft are calculated using the above-described equations (1) and (2) of the Voigt model. In the above formula (2), the current tangential contact force by the spring is calculated as K _t x _t = spring [] + ΔFt. As described above, spring [] is a tangential contact force by a spring one step before, and ΔFt is a tangential contact force by a spring increased (or decreased) from one step before to the present time. ΔFt is calculated by multiplying the relative speed between particles by the time Δt of one step and the elastic coefficient Kt of the spring of the tangential component. The elastic coefficient Kn, Kt of the spring 72 of the normal direction component and the tangential direction component and the repulsion coefficient e between the particles are given values set before the simulation, and the viscous damping coefficients μn, μt of the dashpot 73 are set. Is derived from Equation (3) using given Kn, Kt, e and masses mi, mj included in the particle information. Relative displacements xn, xt and relative velocities vn, vt between particles are calculated based on particle coordinates (position information) and translation / rotation speed (velocity information) included in the particle information. Fn and Ft are calculated by performing the arithmetic processing shown in the equations (1) and (2) using the variables and parameters acquired in this way. Then, using the calculated Fn and Ft, the translational component Fn + Ft and the rotational component Ri × Ft of the contact force are calculated, and these components are arrayed with the contact candidate pair list number (contact candidate number) Lnum as an element. Fijt [Lnum] and Fijr [Lnum] (arrays included in the contact force box described above) are stored.

  Further, in the contact determination, when it is determined that the contact is not made, the contact force is 0, so the arrays Fijt [Lnum] and Fijr [Lnum] having the contact candidate pair list number Lnum as an element (the contact force box described above) ) Stores 0.

The contact force calculation process by the contact determination unit 18 can be expressed by the following algorithm.

  In this way, the contact determination unit 18 creates / updates the contact force box so as to be associated with the contact candidate pair list number Lnum every time the contact candidate list 53 is created / updated. Even if the contact candidate list is not updated, if the particle information is updated, the contact force box is created / updated using the existing contact candidate list.

  The contact determination and contact force are calculated by threading the contact candidate pair list number. The global memory of the GPU can obtain the best access efficiency when the access pattern is in the order of threads. However, since the number of the pair particle is random with respect to the list number, the global memory in which the particle information is stored is accessed at random, and the memory access speed is greatly reduced. On the other hand, since the texture memory has a cache function, a decrease in memory access speed due to random access can be suppressed. Therefore, the global memory area storing the particle information is set so that it can be referred to from the texture memory. Furthermore, since the particles are numbered so as to be in the order of the cell numbers to which they belong, the particle numbers of the contact candidate pairs are close to each other. In addition, the contact candidate pair list is arranged in the list in ascending order of the particle number i. ing. Therefore, by threading with the contact candidate pair list number, the cache hit rate is increased, and the routine for determining contact between particles and calculating the contact force is accelerated.

  The contact determination unit 18 transmits the created / updated contact force box information to the contact force calculation unit 19 together with the accumulated contact candidate number s_jgi [i] and the particle information.

The contact force calculator 19 calculates the translation component Ft _all [i] and the rotation component Fr _all [i] of the total contact force acting on each particle i using the contact force reference table 54 and the accumulated contact candidate number s_jgi. When the contact force reference table 54 is received from the contact force reference table creation unit 17 and the accumulated contact candidate number is received from the contact determination unit 18, first, for each particle i of interest, the contact force reference table 54 and the accumulated contact candidate number s_jgi are obtained. The contact force with the particle in contact with the target particle is extracted from the contact force box. Specifically, for a particle having a particle size larger than the target particle i, based on the box number (Lnum) (reference information) held in the contact force reference table 54, the contact force Fij [Lnum] corresponding to this box number is stored. Is extracted from the contact force box. For a particle having a particle diameter smaller than the particle of interest i, the storage position of the contact force table is specified and the contact force is extracted using the accumulated contact candidate number s_jgi []. More specifically, contact forces Fij [s_jgi [i-1]] to Fij [s_jgi [i] −1 corresponding to the box numbers (contact candidate numbers) of s_jgi [i-1] to s_jgi [i] −1. ] Is extracted from the contact force box. Here, the contact force Fij [] stored in the contact force box is specifically the translation component Fijt [Lnum], the rotation component Fijr [Lnum], and β [Lnum] used in the equation (13) described later. It holds three one-dimensional arrays as a set.

  The process of extracting the contact force from the contact force box will be specifically described with reference to FIG. FIG. 10 shows an example of the target particle i = 4, as in FIGS. The contact force reference table 54 includes reference information of contact force with particles (i = 7, 8) having a particle diameter larger than the particle of interest (i = 4). The box number (= contact candidate list number) Lnum of the contact force box for storing the contact force between the target particle of i = 4 and the particle of i = 7 is 4, and the contact force with the particle of i = 8 is stored. The box number Lnum of the contact force box is 7. The contact force calculation unit 19 uses the box number included in the contact force reference table 54 to contact the contact force Fij [4] with the particle with i = 7 and the particle with i = 8 from the contact force box. Contact force Fij [7] is extracted.

  Further, as shown in FIG. 10, the cumulative contact candidate number s_jgi [i] of the target particle i = 4 is s_jgi [4] = 3, and the particle one before the particle number of the target particle (i in FIG. 10) = 3) The cumulative contact candidate number s_jgi [i-1] is s_jgi [3] = 3. Therefore, for the particles having a particle size smaller than the target particle i = 4 (i = 5, 6 in this example), s_jgi [i-1] to s_jgi [i] −1, that is, the first and second box numbers The contact forces Fij [1] and Fij [2] are extracted from the contact force box.

Next, using the contact force Fij [] (= Fijt [], Fijr [], β []) extracted from the contact force box for each particle i of interest, the translational component Ft _{all of the} total contact force acting on the particle i [i] and the rotation component Fr _all [i] are calculated.

The total contact force (translational component) Ft _all [i] acting on the i-th particle is obtained by adding the contact forces of contact candidate list numbers s_jgi [i-1] to s_jgi [i] -1 as they are. Further, the contact forces referred to from the tables of reference table numbers 0 to n_jli−1 using the contact force reference table 54 (Ref) are added with the opposite sign. And it is expressed as:

Similarly, the total contact force (rotational component) Fr _all [i] acting on the i-th particle is also added. The contact forces with contact candidate list numbers s_jgi [i-1] to s_jgi [i] -1 are referred to as they are, and the contact force reference table 54 (Ref) is used to obtain the reference table numbers 0 to n_jli-1 from the table. The contact force to be referred to is expressed by the following equation by adding β after multiplying it with the opposite sign. Here, β is a value obtained by dividing the distance from the center of the j particle to the contact point by the distance from the center of the i particle to the contact point.

  Here, the contact force of s_jgi [i-1] to s_jgi [i] -1 is the contact force that the particle i receives from the particle j, whereas the contact force reference table 54 of 0 to n_jli-1 The contact force referred to is the contact force that the particle i gives to the particle j. Therefore, with respect to the contact force reference table 54 of No. 0 to n_jli−1, the translation component of the contact force is added after the signs are reversed, and the rotation component of the contact force is multiplied by β [box id]. Add together. In this way, the array Fjt [], Fijr [] in which the contact force is directly applied using s_jgi [i] without using the contact force reference table 54Ref [i] [box id] is applied to the force that the particle i receives from the particle j. ], It is possible to add the contact forces received from the small particles that are a huge number of contact candidate pairs to the large particles without increasing the memory usage of the contact force reference table 54. It was.

In the calculation using the computer, specifically, the particle number i is threaded and a loop is performed with box_id to obtain Ft _all and Fr _all . Here, memory access to Fij is random access, so texture memory is used. As described above, since the texture memory has a high cache hit rate, it can be processed at a higher speed than when the global memory is directly accessed. Regarding the contact force reference table 54, continuous access is possible by declaring the array so that the B element size of the array Ref [A] [B] is a multiple of 16.

The contact force calculation unit 19 sends the translation component Ft _all [i] and the rotation component Fr _all [i] of the total contact force acting on each particle i calculated as described above to the particle information update unit 20 together with the particle information. Send.

When the particle information update unit 20 receives the translation component Ft _all [i] and the rotation component Fr _all [i] of the total contact force acting on each particle i from the contact force calculation unit 19 and the particle information, Based on the contact force received from the position, speed, and contact force calculation unit 19, the coordinates (position information) of the particles, the translation speed and the rotation speed (speed information) are calculated, and the calculated coordinates, translation speed and rotation speed are calculated. To update the particle information held in the particle information holding unit 11. Specifically, the translational velocity and particle coordinates at a new time are obtained by solving the translational equation of motion (the above equation (4)) using the leap-flog method. Furthermore, the equation of motion of the rotation (the above equation (5)) is similarly solved using the leap-flog method, and the rotation speed obtained here is used as a predicted value. A rolling frictional force is calculated with respect to the predicted rotational speed, and a rotational speed at a new time is obtained by solving a rotational equation of motion in consideration of the rotational force due to the contact force and the rolling frictional force.

  The particle information update unit 20 updates the coordinate, translation speed, and rotation speed of the particle information using the calculated position coordinates, translation speed, and rotation speed of each particle, and transmits the updated particle information to the particle information holding unit 11. Then, the particle information in the particle information holding unit 11 is updated to the latest one.

  Next, with reference to FIG. 11, the process performed in the particle simulation apparatus 10 will be described, and the particle simulation method according to the present embodiment will be described. FIG. 11 is a flowchart showing processing executed in the particle simulation apparatus 10 of the present embodiment. The particle information inputted and held in the particle information holding unit 11 is acquired from the particle information holding unit 11 by the particle information acquiring unit 12 (S101: input step). The contact candidate list update determination unit 13 determines whether or not the contact candidate list needs to be updated (S102). For example, the update determination is performed when the particle radius is r and there is a particle having an accumulated movement distance of the particle larger than rα. Here, α is a parameter for adjusting the update frequency. For example, α = 0.2. When it is determined that updating is necessary, the process proceeds to step S103, and when it is determined that updating is not necessary, the process proceeds to step S108.

  Next, the particle number is changed by the particle number changing unit 14 (S103: giving step). The particle number changing unit 14 arranges the particles in the order of the cell numbers to which the particles belong, reassigns the particle numbers in that order, reorders the particles having the same cell numbers in the order of the particle diameters, and reassigns the particle numbers in that order. Subsequently, the neighboring particle table creating unit 15 creates a neighboring particle table 52 including information on neighboring particles located in the vicinity of each particle in the work space 51 (S104). The conditions for storing as neighboring particles in the neighboring particle table 52 are the above-mentioned conditions (a) to (c).

  Next, by the contact candidate list creation unit 16, the cumulative contact candidate number s_jgi that is the prefix sum of the number of particles n_jgi [i] that satisfies the condition (a) of the neighboring particle table 52 but does not satisfy (b) and (c). [i] is calculated (S105: step of calculating the total number of contact candidates), and then a contact candidate list 53 including information on particle pairs that may be in contact with all particles in the work space 51 is created ( S106: Contact candidate selection step). The conditions stored in the contact candidate list 53 are the above-mentioned conditions a) to c). The contact force reference table creation unit 17 creates a contact force reference table 54 (S107: contact force reference table creation step).

  Next, the contact determination unit 18 performs contact determination between two particles with reference to the contact candidate list 53 (S108), and when it is determined that they are in contact, the contact force Fijr between the particles, Fijt is calculated using equations (1) and (2) and stored in an array of contact force boxes (S109: contact force calculation step). When it is determined that the particles are not in contact, the contact forces Fijr and Fijt are set to 0 and stored in the contact force box array.

  The contact force calculation unit 19 refers to the contact force reference table 54 and the accumulated contact candidate number s_jgi [i], and selects other particles that come into contact with the target particle, and makes contact from the corresponding arrays Fijr and Fijt from the contact force box. Power is acquired. Then, the contact force is summed using equations (12) and (13), and the contact forces Fr and Ft acting on each particle are calculated (S110: contact force summation step). Based on the contact force calculated in step S109 and the current particle information, the particle information update unit 20 calculates the coordinates (position information) of each particle, and the translational speed and the rotational speed (speed information). The particle information held in the particle information holding unit 11 is updated using the calculated coordinates, translation speed, and rotation speed (S111: particle information update step).

  Then, it is determined whether the particle simulation has been completed (S112). In the end determination, for example, it is determined that the simulation has ended when a predetermined time has elapsed. If it is determined that the particle simulation has not ended, the process returns to step S101, and the next step is entered. If it is determined that the particle simulation has ended, the process ends.

  Next, with reference to FIG. 12, the memory usage by the particle simulation apparatus 10 of the present embodiment will be described. FIG. 12 shows the memory usage of the present embodiment (solid line) and the comparative example (broken line) when the ratio between the maximum particle size and the minimum particle size of the particle model is changed. The number of particles in the particle model was 1,000,000, and a two-component system in which particles having the maximum particle diameter (maximum particles) and particles having the minimum particle diameter (minimum particles) were mixed. The comparative example is different from the present embodiment in that information on particles smaller than the particle diameter of the target particle i is also stored in the nearby particle table 52 and the contact force reference table 54. The horizontal axis of FIG. 12 represents “maximum particle size / minimum particle size”, and the vertical axis represents “memory material amount (MB)”.

In the case of the comparative example, the memory usage required for the neighborhood particle table 52nei_pn [i] [box_id] and the contact force reference table 54Ref [i] [box_id] is 4 bytes for one data, the number of particles × the maximum contact candidate particle It becomes a number x 4 bytes. Here, the maximum number of contact candidate particles is the maximum number of neighboring particles stored in the neighboring particle table 52 and the contact force reference table 54, and it is only necessary to consider a state in which the surface of the largest particle is filled with the smallest particle. ,
Maximum number of contact candidate particles = 4π (r _max + r _min ) ² / πr _min ² (14)
It becomes. Here, r _max is the particle size of the largest particle, and r _min is the particle size of the smallest particle. That is, in the case of the comparative example, as shown in FIG. 12, the maximum number of contact candidate particles increases as the maximum particle size increases with respect to the minimum particle size, that is, as the maximum particle size / minimum particle size increases. Memory usage increases. For example, under the condition of a particle size ratio of 1:10 (maximum particle size / minimum particle size = 10), the number of contact candidate particles per particle increases to a maximum of about 500, so that the required memory amount is as small as 2 GB. Become.

On the other hand, in the present embodiment, particles having a particle size smaller than the target particle are not considered as neighboring particles and are not stored in the neighboring particle table 52 and the contact force reference table 54. ) Formula does not hold, and the maximum number of contact candidate particles in this embodiment does not depend on the maximum and minimum particle diameters. In the present embodiment, the maximum number of contact candidate particles is only the parameter α included in dis_lim = (r _i + r _j ) (1.0 + α) that is the threshold of the interparticle distance in the above-described neighboring particle determination condition (a). It depends on the particle size distribution and is constant without depending on the particle size distribution. In the case of α = 0.2, it is sufficient that the maximum number of contact candidate particles is 16 under any particle size distribution condition. Therefore, here, the memory usage was calculated with the maximum contact candidate particle number set to 16. As a result, as shown in FIG. 12, the amount of memory used in this embodiment was a constant value of the number of particles × the maximum number of contact candidate particles × 4 bytes = 1000000 × 16 × 4 = 64 MB.

  According to the particle simulation apparatus 10 and the particle simulation method according to the present embodiment, for each of the particles in the work space 51, other particles that are located in the vicinity of the particle i and are likely to come into contact with the particle i. Of the particles smaller than the particle diameter of the particle i, only the cumulative contact candidate number s_jgi [i] indicating the number of particles is retained. And the contact determination part 18 specifies the storing position of a contact force table using this integrated contact candidate number s_jgi [i], and extracts the contact force of the said particle | grain and another particle. Thereby, the contact force reference table 54 does not need to include information on particles smaller than the particle diameter of the particles i among other particles located in the vicinity of the particle i, and relates to particles larger than the particle diameter of the particles i. It only needs to retain information. When a particle group having a particle size distribution is a simulation target, the number of particles that can be contacted at the same time increases when considering contact to particles having a particle size smaller than that of the particle group. The memory area to be increased must be increased. On the other hand, with the above configuration, it is sufficient to consider only contact with particles having a particle diameter larger than that of the self, the number of particles considering contact can be reduced, and memory allocated in advance for creating the contact force reference table 54 The area can be reduced, and as a result, the memory usage can be suppressed. By suppressing the memory usage in this way, it is possible to execute a particle simulation using a shared memory parallel computer such as a GPU.

  Further, the particle number changing unit 14 specifies all the particles in the work space 51 in the order corresponding to the cell numbers, and further specifies the particles specified in the order corresponding to the particle diameter between the particles arranged in the same cell. A number is assigned and each particle is sorted in order of cell number, and further in the order of particle size within the same cell. For this reason, when creating the contact force reference table, particles having a small particle diameter can be easily recognized based on the particle number, and whether or not particles having a small particle diameter can be included in the contact force reference table is determined. There is no need to do it. As a result, calculation efficiency can be improved in the creation of the contact force reference table.

  While the present invention has been described with reference to the preferred embodiment, the present invention is not limited to the above embodiment. For example, in the above-described embodiment, the particle number changing unit 14 sorts (1) the particles in the order of the cell numbers to which the particles belong, and (2) the particles having the same cell numbers in the order of the particle diameters. (1) You may comprise so that only the process which arranges a particle in order of the cell number to which it belongs may be performed.

  In the above-described embodiment, the GPU is cited as an application target of the present invention, but hardware (for example, a vector processor) that performs SIMD (Single Instruction Multiple Data) parallel processing of the shared memory type similar to the GPU, that is, a shared memory It can also be applied to type parallel computers.

  DESCRIPTION OF SYMBOLS 10 ... Particle simulation apparatus, 11 ... Particle information holding part (input means, particle information holding means), 12 ... Particle information acquisition part, 13 ... Contact candidate list update determination part, 14 ... Particle number change part (granting means), 15 ... Near particle table creation unit, 16 ... Contact candidate list creation unit (contact candidate selection unit, integrated contact candidate number calculation unit), 17 ... Contact force reference table creation unit (contact force reference table creation unit), 18 ... Contact determination unit (Contact force calculation means), 19 ... contact force calculation section (contact force total calculation means), 20 ... particle information update section (calculation means, particle information update means), 51 ... work area, 52 ... neighborhood particle table, 53 ... Contact candidate list, 54 ... contact force reference table.

Claims (2)

A particle simulation device that calculates the position / velocity based on the contact force with other particles for a particle group having a particle size distribution in the work space, and simulates the behavior of the particles,
Input means for inputting particle information including position information and velocity information for each particle of the particle group;
Particle information holding means for holding the particle information;
For each of the particles, an imparting unit that imparts specific information that is identified in an order according to the position information;
The target particles are sequentially selected from the particles in the ascending order of the specific information, and a contact candidate pair with another particle having a smaller particle diameter than the target particle that may be in contact with the target particle is selected. Contact candidate selection means;
For each contact candidate pair selected by the contact candidate selection means, the contact force between two particles of the contact candidate pair is calculated based on the position information and the velocity information of these particles, and the calculated contact Contact force calculation means for storing force information in a contact force table according to the order of contact candidate pairs corresponding to the contact force;
For each of the particles, a contact force reference table including reference information for referring to the contact force with a particle larger than the particle diameter of the other particles located in the vicinity of the particle from the contact force table is created. Means for creating a contact force reference table,
For each of the particles, an integrated contact candidate number calculation that calculates an integrated contact candidate number that is an integrated value of the number of contact candidates indicating the number of particles smaller than the particle diameter of other particles located in the vicinity of the particle. Means,
For each of the particles, the contact force with a particle having a larger particle diameter than the particle is extracted from the contact force table using the reference information of the contact force reference table, and the contact force with a particle having a smaller particle diameter than the particle A contact force total calculating means for specifying and extracting the storage position of the contact force table using the accumulated contact candidate number, and calculating the contact force for each particle using the extracted contact force;
Calculation means for calculating new position information and velocity information of the particles based on the contact force for each particle calculated by the contact force total calculation means;
Particle information updating means for updating particle information held in the particle information holding means using the new position information and velocity information calculated by the calculating means;
Equipped with a,
The work space is divided into a plurality of cells, each cell is assigned a cell number,
The particle information includes a cell number of a cell in which the particle is specified, which is specified based on the position information,
The assigning means, for each of the particles, is specified in the order according to the cell number, and between the particles arranged in the same cell, the specific information specified in the order according to the particle diameter is given,
The contact candidate selection unit selects the target particle from the particle group in ascending order of the specific information, and may be in contact with the target particle, and other particles having a particle diameter smaller than the target particle A contact candidate list including contact candidate pair information consisting of a set of and contact candidate numbers assigned in ascending order to this set for all particles in the work space,
The contact force calculation means, for each of the contact candidate pair information included in the contact candidate list, the contact force between the two particles related to the contact candidate pair information in the position information and the velocity information of the particle Calculated based on the contact force information that associates the calculated contact force with the contact candidate number in the contact force table,
The contact force reference table creating means refers to the contact force table for each of the particles, which may be in contact with the particle, with other particles having a particle diameter larger than the particle. As the reference information, a contact force reference table that holds the other particle and the contact candidate number related to the particle in association with the other particle and the particle is created,
The cumulative contact candidate number calculation means calculates the number of contact candidates indicating the number of other particles having a particle diameter smaller than the particle in contact with the particle in the order of the specific information. Calculate the accumulated contact candidate number s_jgi [i] (i indicates the specific information of the particle) that is the accumulated value;
The contact force sum calculating means extracts contact force from the contact force table using the reference information of the contact force reference table for each of the particles, and uses the accumulated contact candidate number s_jgi [i]. The contact candidate numbers s_jgi [i-1] to s_jgi [i] -1 are extracted from the contact force table, and the contact force for each particle is calculated using the extracted contact forces.
A particle simulation apparatus characterized by that. A particle simulation method in a particle simulation apparatus for calculating a position / velocity based on contact force with other particles for a particle group having a particle size distribution in a work space, and simulating particle behavior,
An input step of inputting particle information including position information and velocity information for each particle of the particle group, and holding the particle information in a particle information holding unit;
For each of the particles, an imparting step for imparting specific information that is identified in an order according to the position information;
The target particles are sequentially selected from the particles in the ascending order of the specific information, and a contact candidate pair with another particle having a smaller particle diameter than the target particle that may be in contact with the target particle is selected. A contact candidate selection step;
For each contact candidate pair selected in the contact candidate selection step, a contact force between two particles of the contact candidate pair is calculated based on the position information and the velocity information of these particles, and the calculated contact A contact force calculation step for storing force information in a contact force table according to the order of contact candidate pairs corresponding to the contact force; and
For each of the particles, a contact force reference table including reference information for referring to the contact force with a particle larger than the particle diameter of the other particles located in the vicinity of the particle from the contact force table is created. A contact force reference table creation step,
For each of the particles, an integrated contact candidate number calculation that calculates an integrated contact candidate number that is an integrated value of the number of contact candidates indicating the number of particles smaller than the particle diameter of other particles located in the vicinity of the particle. Steps,
For each of the particles, the contact force with a particle having a larger particle diameter than the particle is extracted from the contact force table using the reference information of the contact force reference table, and the contact force with a particle having a smaller particle diameter than the particle A contact force summation calculating step for identifying and extracting the storage position of the contact force table using the accumulated contact candidate number, and summing up the contact force for each particle using the extracted contact force;
A calculation step of calculating new position information and velocity information of the particles based on the contact force for each particle calculated in the contact force total calculation step;
A particle information update step for updating the particle information held in the particle information holding means using the new position information and velocity information calculated in the calculation step;
Only including,
The work space is divided into a plurality of cells, each cell is assigned a cell number,
The particle information includes a cell number of a cell in which the particle is specified, which is specified based on the position information,
In the assigning step, for each of the particles, identify in an order according to the cell number, and further provide specific information identified in an order according to the particle diameter between the particles arranged in the same cell,
In the contact candidate selection step, the target particle is selected from the particle group in ascending order of the specific information, and may be in contact with the target particle, and other particles having a particle diameter smaller than the target particle A contact candidate list including contact candidate pair information consisting of a set of and contact candidate numbers assigned in ascending order to this set for all particles in the work space,
In the contact force calculation step, for each of the contact candidate pair information included in the contact candidate list, the contact force between two particles related to the contact candidate pair information is used as the position information and the velocity information of the particle. Calculated based on the contact force information that associates the calculated contact force with the contact candidate number in the contact force table,
In the contact force reference table creation step, for referring to the contact force with each of the particles, which may be in contact with the particle, with other particles having a particle diameter larger than the particle, from the contact force table. As the reference information, a contact force reference table that holds the other particle and the contact candidate number related to the particle in association with the other particle and the particle is created,
In the cumulative contact candidate number calculating step, for each of the particles, the number of contact candidates indicating the number of other particles having a particle diameter smaller than the particle that may be in contact with the particle in the order of the specific information. Calculate the accumulated contact candidate number s_jgi [i] (i indicates the specific information of the particle) that is the accumulated value;
In the contact force sum calculation step, for each of the particles, the contact force is extracted from the contact force table using the reference information of the contact force reference table, and the accumulated contact candidate number s_jgi [i] is used. The contact candidate numbers s_jgi [i-1] to s_jgi [i] -1 are extracted from the contact force table, and the contact force for each particle is calculated using the extracted contact forces.
A particle simulation method characterized by the above.

Images