CN104239277A - Rapid discrete-element numerical-value calculating method applicable to GPU (Graphics Processing Unit) scalar-matrix operation - Google Patents
Rapid discrete-element numerical-value calculating method applicable to GPU (Graphics Processing Unit) scalar-matrix operation Download PDFInfo
- Publication number
- CN104239277A CN104239277A CN201410405053.1A CN201410405053A CN104239277A CN 104239277 A CN104239277 A CN 104239277A CN 201410405053 A CN201410405053 A CN 201410405053A CN 104239277 A CN104239277 A CN 104239277A
- Authority
- CN
- China
- Prior art keywords
- matrix
- particle
- adjacent particles
- particles
- discrete
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 239000011159 matrix material Substances 0.000 title claims abstract description 286
- 238000000034 method Methods 0.000 title claims abstract description 35
- 238000012545 processing Methods 0.000 title abstract description 4
- 239000002245 particle Substances 0.000 claims abstract description 322
- 238000004364 calculation method Methods 0.000 claims abstract description 18
- 230000033001 locomotion Effects 0.000 claims abstract description 16
- 238000004088 simulation Methods 0.000 claims abstract description 11
- 238000012216 screening Methods 0.000 claims abstract description 4
- 238000011068 loading method Methods 0.000 claims description 22
- 230000002153 concerted effect Effects 0.000 claims description 9
- 238000006073 displacement reaction Methods 0.000 claims description 9
- 230000000694 effects Effects 0.000 claims description 8
- 238000011156 evaluation Methods 0.000 claims description 8
- 238000013178 mathematical model Methods 0.000 claims description 6
- 238000004422 calculation algorithm Methods 0.000 claims description 4
- 238000007796 conventional method Methods 0.000 claims description 3
- 239000008187 granular material Substances 0.000 claims description 3
- 238000001914 filtration Methods 0.000 abstract 1
- 230000001131 transforming effect Effects 0.000 abstract 1
- 238000011439 discrete element method Methods 0.000 description 8
- 238000010586 diagram Methods 0.000 description 3
- 239000000463 material Substances 0.000 description 2
- 230000001133 acceleration Effects 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000015556 catabolic process Effects 0.000 description 1
- 230000001351 cycling effect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000018109 developmental process Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 238000005272 metallurgy Methods 0.000 description 1
- 238000005065 mining Methods 0.000 description 1
- 239000011435 rock Substances 0.000 description 1
- 230000011218 segmentation Effects 0.000 description 1
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
The invention relates to a rapid discrete-element numerical-value calculating method which is applicable to GPU (Graphics Processing Unit) scalar-matrix operation. The rapid discrete-element numerical-value calculating method comprises the following steps of (1) establishing an adjacent particle matrix and a particle discrete-element accumulating model; storing adjacent particle numbers which are possibly in contact with particles into the corresponding lines of the adjacent particle matrix Pn from NO.1 to NO.m, and filling line-length differences with m+1 dummy particle numbers; (2) realizing the scalar-matrix iterative calculation of particle stress; transforming the coordinates and the attributes of adjacent particles into a m*n matrix mode corresponding to the adjacent particle matrix on the basis of the adjacent particle matrix; obtaining a particle primary-stress matrix Fn0 (the size of the matrix is m*n) through matrix calculation in discrete-element iterative operation; (3) filtering stress calculating results by using a contact-relation matrix to finish iterative calculation; calculating a contact-relation boolean matrix Bc according to factors, such as stress, screening out a practical stress unit in the Fn0 by utilizing the Bc to obtain a practical particle-stress matrix Fn, calculating a resultant force and finishing particle motion simulation.
Description
Technical field
The present invention relates to granular discrete-element simulation, be particularly useful for the quick calculation method of GPU Scalar matrix computing.
Background technology
Distinct element method is a kind of Main Numerical analogy method solving discontinuous media problem, can the large deformation of dynamic similation rock mass, rupture failure, achieve in fields such as ground, mining and metallurgy, environment, pharmacy and apply comparatively widely.For a long time, due to the calculated amount that distinct element method is huge, its count particles, usually at several thousand to several ten thousand, significantly limit development and the application of discrete element.In recent years, GPU (graphics processing unit) is by means of its powerful dense matrix data computing power, and the application in general-purpose computations field is more and more extensive, as imaging of medical, oil-gas exploration and scientific algorithm etc.In floating-point operation, matrix computations etc., and even GPU can provide decades of times hundreds of times in the performance of CPU, can significantly improve numerical evaluation speed.But because discrete element simulation relates to large deformation and rifting, discrete element particle annexation and the reason such as connection quantity is fixing, its computation process lacks Scalar matrix computing feature, is realizing the quick Scalar matrix computing of GPU has difficulties.Therefore, need to adopt specific method to realize the computing of discrete element GPU Scalar matrix, and improve the counting yield of distinct element method.
Summary of the invention
In order to the calculated amount overcoming distinct element method existence is large, the problem that counting yield is low.The present invention seeks to, propose a kind of fast discrete unit numerical computation method being applicable to the computing of GPU Scalar matrix.By setting up specific adjacent particles matrix, realize the GPU computing of Scalar matrix.The method all calculates and all realizes by matrix operation, effectively can improve the counting yield of distinct element method, realize the simulation of large-scale granular discrete-element, and ultra-large granular discrete-element parallel computation can be further used for, realize the Fast simulation of complicated geotechnical engineering problems.
The present invention is the technical scheme solving the quick computational problem proposition of discrete element, a kind of fast discrete unit numerical computation method being applicable to the computing of GPU Scalar matrix.Its step comprises:
(1) adjacent particles matrix and granular discrete-element Mathematical Model of heaped-up is set up; Particle is numbered by 1 to m, the adjacent particles that may contact with each particle (centrophyten) is numbered in the corresponding line being stored in adjacent particles matrix Pn, namely matrix Pn i-th line item be numbered the centrophyten of i adjacent particles numbering (each particle can be centrophyten, it is relative to adjacent particles, and namely these adjacent particles are all near Current central particle); The matrix different rows difference in length empty particle numbering of m+1 is filled;
Each particle can be centrophyten, and it is relative to adjacent particles; Namely these adjacent particles are all near current " centrophyten "; The i.e. adjacent particles numbering of Pn i-th line item numbering i particle.
(2) Scalar matrix iterative computation numerical density is realized; Based on adjacent particles matrix, adjacent particles coordinate and attribute are changed into the m*n matrix form corresponding with adjacent particles matrix.In discrete element interative computation, obtain the preliminary stressed matrix F n of particle by matrix computations
0(matrix size m*n).
(3) use contact relation matrix to filter Force Calculation result, complete iterative computation.Calculate contact relation Boolean matrix Bc according to the factor such as stressed, utilize Bc to filter out Fn
0in actual loading unit, obtain particle actual loading matrix, calculate make a concerted effort and complete movement of particles simulate.
Step 10 builds initial granular discrete-element Mathematical Model of heaped-up, particle by 1 open numbering to m; According to particle numbering, the geometric parameter of particle and mechanics parameter are stored in each array;
Step 11 by conventional method obtain each particle adjacent particles, namely may contact particle, and it is stored in line by line in adjacent particles matrix Pn; The adjacent particles numbering of each particle is stored in adjacent particles matrix corresponding line, and matrix size is m*n, i.e. the capable n row of m;
Particle habit is changed into the m*n form corresponding with Pn matrix by step 12, and carries out Scalar matrix discrete element numerical evaluation (specifically seeing step 20-26), obtains the preliminary stressed matrix F n of particle
0(m*n); Particle numbering one_to_one corresponding in preliminary stressed matrix and adjacent particles matrix, have recorded the primary action power between particle and its adjacent particles;
Step 13 uses contact relation Boolean matrix Bc to record connection and the contact condition of particle; Contact relation matrix and adjacent particles matrix one_to_one corresponding, whether contact relation matrix record adjacent particles contacts with centrophyten, namely whether there is the effect of power.Count particles contact relation Boolean matrix Bc (specifically seeing step 30-36), because in adjacent particles matrix, only partial particulate and centrophyten have contact relation, by preliminary stressed matrix F n
0be multiplied by Bc item by item and obtain particle actual loading matrix F n (m*n matrix).
Step 14 calculates with joint efforts according to discrete element particle and motion calculation carries out movement of particles simulation.Ask for suffered by each particle according to actual loading matrix and make a concerted effort, and by the motion of classical mechanics method count particles.(specifically seeing step 30-36)
Step 15, by continuous iteration, realizes the dynamic similation of granular discrete-element.
Step 16 obtains numerical simulation net result.
Step 20-26 describes the preliminary stressed matrix procedures of step 11-12 Scalar matrix iterative computation particle in detail:
Step 20 inputs the adjacent particles array of each particle, and the attribute array of all particles, comprises radius (R), coordinate (P [X, Y, Z]), positive stiffness factor (Kn) and quality (M), array size is m*1, and namely m capable 1 arranges.
Step 21 builds adjacent particles matrix Pn.The adjacent particles array of particle is stored in adjacent particles matrix by order line by line by number.Particle matrix Pn is m*n matrix, namely has m capable and n row, and the i-th line number value representative is numbered the adjacent particles numbering of i particle.Because the adjacent particles number of each particle is different, the width n of Pn matrix is determined by the adjacent particles array row grown most.There is dummy cell in the row that in Pn, length is shorter, fills with empty particle numbering m+1;
Step 22 increases the empty particle habit value of m+1 at attribute array end, and calculates attribute matrix corresponding to adjacent particles matrix (m*n).Owing to there is m+1 element in Pn, need in attribute array finally for the empty particle of numbering m+1 increases a property value, to ensure matrix algorithms true(-)running; Pn is updated in attribute array and obtains adjacent particles attribute matrix.
Step 23 is according to particle coordinates matrix and adjacent particles matrix, and obtain the x coordinate difference value matrix dPnX of centrophyten and adjacent particles, dPnY, dPnZ, its size is m*n.
Step 24 calculates centrophyten and adjacent particles spacing matrix dPn:
DPn=sqrt (dPnX.^2+dPnY.^2+dPnZ.^2); (by gauged distance formula)
Step 25 obtains relative displacement Sn between particle further:
Sn＝dPn-(R*ones(1,n)+R’(Pn))；
Wherein R*ones (1, n) be particle radius array (m*1); R ' (Pn) is the radius array increasing empty particle; The right Section 2 R*ones (radius sum of particle and adjacent particles centered by 1, n)+R ' (Pn).One binomial difference is the distance of particle surface, namely obtains relative displacement matrix (size m*n) between particle.
Preliminary positive force (Fn between step 26 count particles
0), i.e. the product item by item of positive relative displacement matrix (Sn) between stiffness factor matrix (Kn) and particle:
Fn
0＝Kn.*Sn；
In formula, Kn and Sn is m*n matrix, and in the two matrix, element is multiplied item by item, obtains the preliminary positive force matrix F n that size is m*n
0.
Step 30-36 describes step 13-14 contact relation Boolean matrix in detail and calculates actual loading and discrete element calculation process:
Step 30 inputs the preliminary positive force matrix F n of particle
0boolean matrix Bb (being m*n matrix) is connected with particle.Bb and adjacent particles matrix one_to_one corresponding, have recorded between particle and adjacent particles whether have connecting forces, namely whether allow pulling force effect.
Step 31 is according to Fn
0upgrade particle connection matrix Bb, when pulling force exceedes particular value (stretching resistance matrix value) when between particle, between particle, connect fracture.Realized by matrix Boolean calculation:
Bb=(Bb & Fn
0<Fn
max), the instruction of Matlab matrix, lower with;
Wherein Fn
maxfor stretching resistance matrix (m*n) between particle.
Step 32 is according to Fn
0acting force (pressure is negative) Boolean matrix Bp is pressed between count particles:
Bp＝(Fn
0<0)；
Step 33 count particles Contact relation Boolean matrix Bc, namely indicates whether have actual force between particle.Between particle, actual force comprises pressure and pulling force, and therefore, Bc obtains by getting union to particle connection matrix Bb and pressure acting force matrix Bp:
Bc＝(Bb|Bp)；
Step 34 is according to Fn
0actual loading Fn between count particles.In adjacent particles matrix, also not all adjacent particles all has contact relation with present granule, obtains actual forward power by the screening of contact relation Boolean matrix:
Fn＝Fn
0.*Bc；
In contact relation matrix, have contact to be 1, contactless is 0.
Fn is decomposed into FnX by step 35, FnY, FnZ vector form.Coordinate difference value matrix dPnX between the particle obtained in integrating step 23, dPnY, dPnZ (i.e. positive force direction), can be decomposed into coordinate vector form D nX by Fn, FnY, FnZ.Particle is suffered makes a concerted effort for adjacent particles is to its acting force sum, namely get on actual loading row matrix direction and:
FX=sum (FnX, 2); The instruction of Matlab matrix, lower same;
FY＝sum(FnY,2)；
FZ＝sum(FnZ,2)；
Numerical density is F=[FX, FY, FZ] (m*3 matrix).
Step 36 carrys out the component motion of count particles in xyz direction according to classical mechanics.
Distinct element method is a kind of Computational Mechanics method, and finite element method is similar.Material breakdown is become a series of and has specific discrete element particle by it, and by Newtonian mechanics come between count particles by force and motion, thus the deformation and failure of simulation material.
Step 20-26 is the refinement of step 11-12, and step 30-36 is the refinement of step 13-14.Fig. 1 is overall procedure.
The invention has the beneficial effects as follows: realize the calculating of all discrete elements and all complete with Scalar matrix form of calculation, to realize quick GPU computing, significantly improve the efficiency that discrete element calculates.This method avoids in conventional discrete unit computing method need constantly to carry out cycling to calculate each particle contact, by force and motion etc., but by achieving distinct element method Scalar matrix rapid computations based on adjacent particles matrix and contact relation Boolean matrix.By applying these methods, making discrete element calculate the speed that can reach general business software tens times, single personal computer can realize the dynamic similation of 1,000,000 particles.And based on the discrete element numerical evaluation of Scalar matrix, computed segmentation can be carried out easily, be conducive to the parallel computation of further ultra-large discrete element.
The present invention, by obtaining the actual loading matrix of particle, has calculated particle dynamic similation finally by conventional discrete unit method.In numerical procedure, all properties data are all recorded in corresponding matrix based on adjacent particles matrix, and adopt Scalar matrix computing.Therefore, the method is applicable to the matrix operation of GPU (graphics processing unit) high speed, significantly improves the speed of discrete element numerical method, can realize the dynamic similation of 1,000,000 particles on single personal computer.
Accompanying drawing explanation
Fig. 1 fast discrete unit's numerical evaluation general structure and principle of work process flow diagram;
The computing of Fig. 2 Scalar matrix obtains the preliminary stressed matrix procedures figure of particle;
Fig. 3 contact relation matrix computations actual loading and distinct element method calculation flow chart.
Embodiment
Fig. 1 is general structure and the principle of work of this method enforcement.Be characterized in that realizing Scalar matrix discrete element in conjunction with adjacent particles matrix Pn and contact relation Boolean matrix Bc calculates.
Step 10 builds initial granular discrete-element Mathematical Model of heaped-up, particle by 1 open numbering to m.According to numbering, the geometric parameter of particle and mechanics parameter are stored in each array.
Step 11 obtains the adjacent particles (namely may contact particle) of each particle by conventional method, and it is stored in line by line in adjacent particles matrix Pn.Accompanying drawing 1a gives adjacent particles matrix, and the adjacent particles numbering of each particle is stored in matrix corresponding line, and matrix size is m*n, i.e. the capable n row of m.
Particle habit is changed into the m*n form corresponding with Pn matrix by step 12, and carries out Scalar matrix discrete element numerical evaluation (specifically seeing step 20-26), obtains the preliminary stressed matrix F n of particle
0(m*n).Particle numbering one_to_one corresponding in preliminary stressed matrix and adjacent particles matrix, have recorded the primary action power between particle and its adjacent particles.
Step 13 count particles contact relation Boolean matrix Bc, and obtain particle actual loading matrix F n.Due in adjacent particles matrix only partial particulate and centrophyten have contact relation, contact relation Boolean matrix Bc must be used to record connection and the contact condition of particle.As accompanying drawing 1b, contact relation matrix and adjacent particles matrix one_to_one corresponding, and record adjacent particles and whether contact with centrophyten, namely whether there is the effect of power.Further, by preliminary stressed matrix F n
0be multiplied by Bc item by item and obtain particle actual loading matrix F n (m*n matrix).
Step 14 is made a concerted effort and motion calculation according to discrete element particle.Ask for suffered by each particle according to actual loading matrix and make a concerted effort, and by the motion of classical mechanics method count particles.
Step 15, by continuous iteration, realizes the dynamic similation of granular discrete-element.
Step 16 obtains numerical simulation net result.
Fig. 2 is the preliminary stressed matrix procedures figure of Scalar matrix iterative computation particle.Be characterized in building corresponding particle habit matrix according to adjacent particles matrix, and obtain particle and adjacent particles intermolecular forces matrix by Scalar matrix computing.Step 20-26 is the refinement explanation of step 11-12.
Step 20 inputs the adjacent particles array of each particle, and the attribute array of all particles, such as radius (R), coordinate (P [X, Y, Z]), positive stiffness factor (Kn) and quality (M) etc., array size is m*1, and namely m capable 1 arranges.
Step 21 builds adjacent particles matrix Pn.The adjacent particles array of particle is stored in adjacent particles matrix by order line by line by number.Particle matrix Pn is m*n matrix, namely has m capable and n row, and the i-th line number value representative is numbered the adjacent particles numbering of i particle.Because the adjacent particles number of each particle is different, the width n of Pn matrix is determined by the adjacent particles array row grown most.There is dummy cell in the row that in Pn, length is shorter, fills with empty particle numbering m+1.As accompanying drawing 2a, particle 2 (the second row) has maximum adjacent particles, and its last element is 5.Other row relevant position is then filled with numbering m+1.Because particle numbering is by 1 to m, in further matrix operation, numbering m+1 can be distinguished.
Step 22 increases the empty particle habit value of m+1 at attribute array end, and calculates attribute matrix corresponding to adjacent particles matrix (m*n).Owing to there is m+1 element in Pn, need in attribute array finally for the empty particle of numbering m+1 increases a property value, to ensure matrix algorithms true(-)running.Further, Pn is updated in attribute array obtains attribute matrix.Such as particle x coordinate array is PX ' (m+1 element), then corresponding adjacent particles x coordinates matrix is PX ' (Pn) (instruction of Matlab matrix).
Step 23 is according to particle coordinates matrix and adjacent particles matrix, and obtain the x coordinate difference value matrix dPnX of centrophyten and adjacent particles, dPnY, dPnZ, its size is m*n.Specifically obtain (for dPnX) by following formula:
DPnX=PX ' (Pn)-PX*ones (1, n); (Matlab matrix instruction, lower same)
On the right of formula, Section 1 obtains adjacent particles coordinate in the x direction; Section 2 obtains centrophyten coordinate in the x direction.PX is m*1 matrix, ones (1, n) obtain being worth be entirely 1 1*n matrix, the two is multiplied and obtains m*n matrix, the x coordinate figure of its every behavior centrophyten; Namely one binomial difference obtains adjacent particles and centrophyten matrix of differences in the x direction.Similarly, adjacent particles and centrophyten matrix of differences dPnY in the y and z directions and dPnZ can be obtained.
Step 24 calculates centrophyten and adjacent particles spacing matrix dPn:
dPn＝sqrt(dPnX.^2+dPnY.^2+dPnZ.^2)；
Step 25 obtains relative displacement Sn between particle further:
Sn＝dPn-(R*ones(1,n)+R’(Pn))；
Wherein R is particle radius array (m*1); R ' is for increasing the radius array of empty particle; The radius sum of particle and adjacent particles centered by the Section 2 of the right.One binomial difference is the distance of particle surface, namely obtains relative displacement matrix (size m*n) between particle.
Preliminary positive force (Fn between step 26 count particles
0), i.e. the product item by item of positive relative displacement matrix (Sn) between stiffness factor matrix (Kn) and particle:
Fn
0＝Kn.*Sn；
In formula, Kn and Sn is m*n matrix, and in the two matrix, element is multiplied item by item, obtains the preliminary positive force matrix F n that size is m*n
0.Accompanying drawing 2b is preliminary positive force matrix schematic diagram, itself and adjacent particles matrix one_to_one corresponding, and have recorded the positive acting power between particle and adjacent particles.Preliminary positive force matrix is intergranular possibility acting force, because the unit in adjacent particles matrix may not contact with centrophyten (namely not having pressure and tension force effect), need to adopt contact relation Boolean matrix to filter out actual acting force further.
Fig. 3 is contact relation matrix computations actual loading and discrete element calculation flow chart.Be characterized in setting up contact relation Boolean matrix to filter preliminary positive force matrix, obtain intergranular actual force, and complete discrete element numerical evaluation.Step 30-36 is the refinement explanation of step 13-14.
Step 30 inputs the preliminary positive force matrix F n of particle
0boolean matrix Bb (being m*n matrix) is connected with particle.Bb and adjacent particles matrix one_to_one corresponding, have recorded between particle and adjacent particles whether have connecting forces, namely whether allow pulling force effect.
Step 31 is according to Fn
0upgrade particle connection matrix Bb, when pulling force exceedes particular value when between particle, between particle, connect fracture.Realize by matrix Boolean calculation:
Bb=(Bb & Fn
0<Fn
max) (Matlab matrix instruction, lower same)
Wherein Fn
maxfor stretching resistance matrix (m*n) between particle.
Step 32 is according to Fn
0acting force (pressure is negative) Boolean matrix Bp is pressed between count particles:
Bp＝(Fn
0<0)；
Step 33 count particles Contact relation Boolean matrix Bc, namely indicates whether have actual force between particle.Between particle, actual force comprises pressure and pulling force, and therefore, Bc obtains by getting union to particle connection matrix Bb and pressure acting force matrix Bp:
Bc＝(Bb|Bp)；
Accompanying drawing 3a is Bc matrix schematic diagram, and its size is m*n, with adjacent particles matrix element one_to_one corresponding, and have recorded the effect that whether there is contact and power between particle.
Step 34 is according to Fn
0actual loading Fn between count particles.In adjacent particles matrix, also not all adjacent particles all has contact relation with present granule, obtains actual forward power by the screening of contact relation Boolean matrix:
Fn＝Fn
0.*Bc；
In contact relation matrix, have contact to be 1, contactless is 0.Preliminary positive force Fn in accompanying drawing 2b
0be multiplied item by item with the particle contacts Boolean matrix in accompanying drawing 3a, then obtain actual positive force matrix F n between the particle in accompanying drawing 3b.
Fn is decomposed into FnX by step 35, FnY, FnZ vector form.Coordinate difference value matrix dPnX between the particle obtained in integrating step 23, dPnY, dPnZ (i.e. positive force direction), can be decomposed into coordinate vector form D nX by Fn, FnY, FnZ.Particle is suffered makes a concerted effort for adjacent particles is to its acting force sum, namely get on actual loading row matrix direction and:
FX=sum (FnX, 2); (Matlab matrix instruction, lower same)
FY＝sum(FnY,2)；
FZ＝sum(FnZ,2)；
Numerical density is F=[FX, FY, FZ] (m*3 matrix).
Step 36 carrys out the component motion of count particles in xyz direction according to classical mechanics.As component is calculated as in the x direction:
AX＝FX./M；
VX＝VX+AX*dT；
PX＝PX+VX*dT；
In formula, M is granular mass array, and AX is x directional acceleration, and VX is x direction speed, and PX is particle x direction coordinate, and dT is iteration time step-length.
Example above gives the method realizing the calculating of discrete element positive force and the computing of interative computation Scalar matrix.Similarly, between particle, tangential force and torque etc. also can be calculated by similar approach, thus realize the Scalar matrix computing of granular discrete-element.
Claims (4)
1. be applicable to the fast discrete unit numerical computation method of GPU Scalar matrix computing, it is characterized in that step comprises:
(1) adjacent particles matrix and granular discrete-element Mathematical Model of heaped-up is set up; By particle by 1 to m numbering, may number in the corresponding line being stored in adjacent particles matrix Pn with the adjacent particles of particle contacts, the line length difference empty particle numbering of m+1 is filled;
(2) Scalar matrix iterative computation numerical density is realized; Based on adjacent particles matrix, adjacent particles coordinate and attribute are changed into the m*n matrix form corresponding with adjacent particles matrix; In discrete element interative computation, obtain the preliminary stressed matrix F n of particle by matrix computations
0, matrix size m*n;
(3) use contact relation matrix to filter Force Calculation result, complete iterative computation; Calculate contact relation Boolean matrix Bc according to the factor such as stressed, utilize Bc to filter out Fn
0in actual loading unit, obtain particle actual loading matrix F n, calculate make a concerted effort and complete movement of particles simulate.
2. the fast discrete unit numerical computation method being applicable to the computing of GPU Scalar matrix according to claim 1, is characterized in that (1) sets up the step of adjacent particles matrix and granular discrete-element Mathematical Model of heaped-up:
Step 10 builds initial granular discrete-element Mathematical Model of heaped-up, particle by 1 open numbering to m; According to particle numbering, the geometric parameter of particle and mechanics parameter are stored in each array;
Step 11 by conventional method obtain each particle adjacent particles, namely may contact particle, and it is stored in line by line in adjacent particles matrix Pn; The adjacent particles numbering of each particle is stored in adjacent particles matrix corresponding line, and matrix size is m*n, i.e. the capable n row of m;
Particle habit is changed into the m*n form corresponding with Pn matrix by step 12, and carries out Scalar matrix discrete element numerical evaluation, obtains the preliminary stressed matrix F n of particle
0, m*n rank; Particle numbering one_to_one corresponding in preliminary stressed matrix and adjacent particles matrix, have recorded the primary action power between particle and its adjacent particles;
Step 13 uses contact relation Boolean matrix Bc to record connection and the contact condition of particle; Contact relation matrix and adjacent particles matrix one_to_one corresponding, whether contact relation matrix record adjacent particles contacts with centrophyten, namely whether there is the effect of power; Count particles contact relation Boolean matrix Bc, because in adjacent particles matrix, only partial particulate and centrophyten have contact relation, by preliminary stressed matrix F n
0be multiplied by Bc item by item and obtain particle actual loading matrix F n, m*n matrix;
Step 14 calculates with joint efforts according to discrete element particle and motion calculation carries out movement of particles simulation; Ask for suffered by each particle according to actual loading matrix and make a concerted effort, and by the motion of classical mechanics method count particles; Specifically see step 30-36;
Step 15, by continuous iteration, realizes the dynamic similation of granular discrete-element.
Step 16 obtains numerical simulation net result.
3. the fast discrete unit numerical computation method being applicable to the computing of GPU Scalar matrix according to claim 1, is characterized in that the preliminary stressed matrix procedures of Scalar matrix iterative computation particle; Build corresponding particle habit matrix according to adjacent particles matrix, and obtain particle and adjacent particles intermolecular forces matrix by Scalar matrix computing;
Step 20 inputs the adjacent particles array of each particle, and the attribute array of all particles, comprises radius (R), coordinate (P [X, Y, Z]), positive stiffness factor (Kn) and quality (M), array size is m*1, and namely m capable 1 arranges;
Step 21 builds adjacent particles matrix Pn; The adjacent particles array of particle is stored in adjacent particles matrix by order line by line by number; Particle matrix Pn is m*n matrix, namely has m capable and n row, and the i-th line number value representative is numbered the adjacent particles numbering of i particle; Because the adjacent particles number of each particle is different, the width n of Pn matrix is determined by the adjacent particles array row grown most; There is dummy cell in the row that in Pn, length is shorter, fills with empty particle numbering m+1; Step 22 increases the empty particle habit value of m+1 at attribute array end, and calculates attribute matrix corresponding to adjacent particles matrix (m*n).Owing to there is m+1 element in Pn, need in attribute array finally for the empty particle of numbering m+1 increases a property value, to ensure matrix algorithms true(-)running; Pn is updated in attribute array and obtains adjacent particles attribute matrix; Step 23 is according to particle coordinates matrix and adjacent particles matrix, and obtain the x coordinate difference value matrix dPnX of centrophyten and adjacent particles, dPnY, dPnZ, its size is m*n; Step 24 calculates centrophyten and adjacent particles spacing matrix dPn:
dPn＝sqrt(dPnX
2+dPnY
2+dPnZ
2)；
Step 25 obtains relative displacement Sn between particle further:
Sn＝dPn-(R*ones(1,n)+R’(Pn))；
Wherein R*ones (1, n) be particle radius array (m*1); R ' (Pn) is the radius array increasing empty particle; The right Section 2 R*ones (radius sum of particle and adjacent particles centered by 1, n)+R ' (Pn); One binomial difference is the distance of particle surface, namely obtains relative displacement matrix (size m*n) between particle;
Preliminary positive force (Fn between step 26 count particles
0), i.e. the product item by item of positive relative displacement matrix (Sn) between stiffness factor matrix (Kn) and particle:
Fn
0＝Kn.*Sn；
In formula, Kn and Sn is m*n matrix, and in the two matrix, element is multiplied item by item, obtains the preliminary positive force matrix F n that size is m*n
0.
4. the fast discrete unit numerical computation method being applicable to the computing of GPU Scalar matrix according to claim 1, it is characterized in that adopting contact relation Boolean matrix to filter out actual acting force, contact relation Boolean matrix calculates actual loading and discrete element calculation process;
Set up contact relation Boolean matrix to filter preliminary positive force matrix, obtain intergranular actual force, and complete discrete element numerical evaluation:
Step 30 inputs the preliminary positive force matrix F n of particle
0boolean matrix Bb (being m*n matrix) is connected with particle; Bb and adjacent particles matrix one_to_one corresponding, have recorded between particle and adjacent particles whether have connecting forces, namely whether allow pulling force effect;
Step 31 is according to Fn
0upgrade particle connection matrix Bb, when pulling force exceedes particular value (stretching resistance matrix value) when between particle, between particle, connect fracture; Realized by matrix Boolean calculation:
Bb=(Bb & Fn
0<Fn
max), the instruction of Matlab matrix, lower with;
Wherein Fn
maxfor stretching resistance matrix (m*n) between particle;
Step 32 is according to Fn
0acting force (pressure is negative) Boolean matrix Bp is pressed between count particles:
Bp＝(Fn
0<0)；
Step 33 count particles Contact relation Boolean matrix Bc, namely indicates whether have actual force between particle; Between particle, actual force comprises pressure and pulling force, and therefore, Bc obtains by getting union to particle connection matrix Bb and pressure acting force matrix Bp:
Bc=(Bb|Bp); Step 34 is according to Fn
0actual loading Fn between count particles; In adjacent particles matrix, also not all adjacent particles all has contact relation with present granule, obtains actual forward power by the screening of contact relation Boolean matrix:
Fn＝Fn
0.*Bc；
In contact relation matrix, have contact to be 1, contactless is 0;
Fn is decomposed into FnX by step 35, FnY, FnZ vector form; Coordinate difference value matrix dPnX between the particle obtained in integrating step 23, dPnY, dPnZ (i.e. positive force direction), be decomposed into coordinate vector form D nX, FnY, FnZ by Fn; Particle is suffered makes a concerted effort for adjacent particles is to its acting force sum, namely get on actual loading row matrix direction and:
FX=sum (FnX, 2); The instruction of Matlab matrix, lower same;
FY＝sum(FnY,2)；
FZ＝sum(FnZ,2)；
Numerical density is F=[FX, FY, FZ] (m*3 matrix);
Step 36 carrys out the component motion of count particles in xyz direction according to classical mechanics.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201410405053.1A CN104239277B (en) | 2014-08-15 | 2014-08-15 | Rapid discrete-element numerical-value calculating method applicable to GPU (Graphics Processing Unit) scalar-matrix operation |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201410405053.1A CN104239277B (en) | 2014-08-15 | 2014-08-15 | Rapid discrete-element numerical-value calculating method applicable to GPU (Graphics Processing Unit) scalar-matrix operation |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN104239277A true CN104239277A (en) | 2014-12-24 |
| CN104239277B CN104239277B (en) | 2017-04-12 |
Family
ID=52227375
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201410405053.1A Active CN104239277B (en) | 2014-08-15 | 2014-08-15 | Rapid discrete-element numerical-value calculating method applicable to GPU (Graphics Processing Unit) scalar-matrix operation |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN104239277B (en) |
Cited By (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105866855A (en) * | 2015-01-22 | 2016-08-17 | 中国石油天然气股份有限公司 | Analysis method of geological tectonic evolution and deformation process |
| CN106201931A (en) * | 2016-08-10 | 2016-12-07 | 长沙中部芯空微电子研究所有限公司 | A kind of hypervelocity matrix operations coprocessor system |
| CN107330227A (en) * | 2017-07-31 | 2017-11-07 | 南京大学 | Consider the discrete Meta Model and method for numerical simulation of the Rock And Soil triaxial test of film effect |
| CN109977505A (en) * | 2019-03-13 | 2019-07-05 | 南京大学 | The discrete element pore system method for fast searching calculated based on GPU matrix |
| CN110929430A (en) * | 2019-12-31 | 2020-03-27 | 浙江交通职业技术学院 | Modulus calculation method for randomly filling silicon carbide particles in diesel engine particle catcher |
| WO2021088099A1 (en) * | 2019-11-04 | 2021-05-14 | 南京大学 | Total matrix computation-based discrete element neighbor search and solution finding method and system |
| CN113656656A (en) * | 2021-07-21 | 2021-11-16 | 南京南力科技有限公司 | Efficient neighbor retrieval method and system for wide-gradation discrete element particle system |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20130297270A1 (en) * | 2010-09-15 | 2013-11-07 | Commonwealth Scientific And Industrial Research Organisation | Discrete element method |
| CN103984829A (en) * | 2014-05-22 | 2014-08-13 | 李笑宇 | Discrete element method based method for improving particle discrete contact detection efficiency |
-
2014
- 2014-08-15 CN CN201410405053.1A patent/CN104239277B/en active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20130297270A1 (en) * | 2010-09-15 | 2013-11-07 | Commonwealth Scientific And Industrial Research Organisation | Discrete element method |
| CN103984829A (en) * | 2014-05-22 | 2014-08-13 | 李笑宇 | Discrete element method based method for improving particle discrete contact detection efficiency |
Non-Patent Citations (4)
| Title |
|---|
| LIN ZHANG ETC.: "A fast scalable implementation of the two-dimensional triangular Discrete Element Method on a GPU platform", 《ADVANCES IN ENGINEERING SOFTWARE》 * |
| 刘凯欣等: "离散元法研究的评述", 《力学进展》 * |
| 孙 鹏等: "基于离散元理论的振动筛分数值模拟程序开发", 《中国制造业信息化》 * |
| 常新正: "基于GPU的颗粒离散元计算方法研究", 《中国优秀硕士学位论文全文数据库 信息科技辑》 * |
Cited By (12)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105866855A (en) * | 2015-01-22 | 2016-08-17 | 中国石油天然气股份有限公司 | Analysis method of geological tectonic evolution and deformation process |
| CN105866855B (en) * | 2015-01-22 | 2018-04-06 | 中国石油天然气股份有限公司 | A kind of analysis method of geologic-tectonic evolution and deformation process |
| CN106201931A (en) * | 2016-08-10 | 2016-12-07 | 长沙中部芯空微电子研究所有限公司 | A kind of hypervelocity matrix operations coprocessor system |
| CN106201931B (en) * | 2016-08-10 | 2019-01-18 | 长沙中部翼天智能装备科技有限公司 | A kind of hypervelocity matrix operation coprocessor system |
| CN107330227A (en) * | 2017-07-31 | 2017-11-07 | 南京大学 | Consider the discrete Meta Model and method for numerical simulation of the Rock And Soil triaxial test of film effect |
| CN109977505A (en) * | 2019-03-13 | 2019-07-05 | 南京大学 | The discrete element pore system method for fast searching calculated based on GPU matrix |
| CN109977505B (en) * | 2019-03-13 | 2024-02-09 | 南京大学 | Discrete element pore system quick searching method based on GPU matrix calculation |
| WO2021088099A1 (en) * | 2019-11-04 | 2021-05-14 | 南京大学 | Total matrix computation-based discrete element neighbor search and solution finding method and system |
| CN110929430A (en) * | 2019-12-31 | 2020-03-27 | 浙江交通职业技术学院 | Modulus calculation method for randomly filling silicon carbide particles in diesel engine particle catcher |
| CN110929430B (en) * | 2019-12-31 | 2023-04-11 | 浙江交通职业技术学院 | Modulus calculation method for randomly filling silicon carbide particles in diesel engine particle catcher |
| CN113656656A (en) * | 2021-07-21 | 2021-11-16 | 南京南力科技有限公司 | Efficient neighbor retrieval method and system for wide-gradation discrete element particle system |
| CN113656656B (en) * | 2021-07-21 | 2024-01-26 | 南京南力科技有限公司 | Efficient neighbor retrieval method and system for wide-grading discrete element particle system |
Also Published As
| Publication number | Publication date |
|---|---|
| CN104239277B (en) | 2017-04-12 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN104239277A (en) | Rapid discrete-element numerical-value calculating method applicable to GPU (Graphics Processing Unit) scalar-matrix operation | |
| Zhang et al. | GPU-accelerated smoothed particle finite element method for large deformation analysis in geomechanics | |
| Recuero et al. | A high-fidelity approach for vehicle mobility simulation: Nonlinear finite element tires operating on granular material | |
| JP6009075B2 (en) | Particle flow simulation system and method | |
| Xiao et al. | A novel CNN-based Poisson solver for fluid simulation | |
| Huang et al. | Implementation of multi-GPU based lattice Boltzmann method for flow through porous media | |
| Ren et al. | Dual-support smoothed particle hydrodynamics in solid: variational principle and implicit formulation | |
| US10354099B2 (en) | Particle simulation device, particle simulation method, and particle simulation program | |
| CN109033537A (en) | The calculation method and system of rock-fill concrete casting process numerical simulation | |
| Koric et al. | Evaluation of massively parallel linear sparse solvers on unstructured finite element meshes | |
| Cui et al. | Metal forming analysis using the edge-based smoothed finite element method | |
| Gharti et al. | Simulation of multistage excavation based on a 3D spectral-element method | |
| CN103631981A (en) | Designing a modeled volume represented by dexels | |
| Lu et al. | Time-discontinuous material point method for transient problems | |
| Liu et al. | A new discrete element-embedded finite element method for transient deformation, movement and heat transfer in packed bed | |
| Beyabanaki et al. | Non-rigid disk-based DDA with a new contact model | |
| Bisson et al. | Multiscale hemodynamics using GPU clusters | |
| CN113011072B (en) | Discrete element complex model identification method based on MIDAS-PFC3D | |
| Lai et al. | Machine‐learning‐enabled discrete element method: Contact detection and resolution of irregular‐shaped particles | |
| Barreto et al. | A guide to modeling the geotechnical behavior of soils using the discrete element method | |
| Hsieh et al. | ESFM: An essential software framework for meshfree methods | |
| Wang et al. | Lattice Boltzmann method with immersed spring boundaries for flow around deformable porous media | |
| Wen et al. | Accelerating Polyhedral Discrete Element Method with CUDA | |
| Va et al. | Parallel Cloth Simulation Using OpenGL Shading Language. | |
| Pantalé et al. | Development of an object-oriented finite element program: application to metal-forming and impact simulations |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |