CN104239277B - 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
- CN104239277B CN104239277B CN201410405053.1A CN201410405053A CN104239277B CN 104239277 B CN104239277 B CN 104239277B CN 201410405053 A CN201410405053 A CN 201410405053A CN 104239277 B CN104239277 B CN 104239277B
- Authority
- CN
- China
- Prior art keywords
- matrix
- granule
- adjacent particles
- discrete
- particle
- 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.)
- Active
Links
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, is particularly suited for the quick calculation method of GPU Scalar matrix computings.
Background technology
Distinct element method is a kind of Main Numerical analogy method for solving the problems, such as discontinuous media, can be with dynamic analog rock mass
Large deformation, rupture failure, in fields such as ground, mining and metallurgy, environment, pharmacy relatively broad application is achieved.For a long time, due to
The huge amount of calculation of distinct element method, it calculates granule generally thousand of to tens of thousands of, significantly limit the development of discrete element with
Using.In recent years, GPU (graphics processing unit) relies on its powerful dense matrix data computing capability, in general-purpose computations neck
The application in domain is more and more extensive, such as imaging of medical, oil-gas exploration and scientific algorithm.At aspects such as floating-point operation, matrix calculus,
GPU can provide the performance of decades of times or even up to a hundred times of CPU, can significantly improve numerical computations speed.But, due to discrete
Unit's simulation is related to large deformation and abstriction, discrete element granule annexation and unfixed etc. reason of connection quantity, and it was calculated
Journey lacks Scalar matrix computing feature, has difficulties on the quick Scalar matrix computings of GPU are realized.Accordingly, it would be desirable to using specific side
Method improves the computational efficiency of distinct element method to realize discrete element GPU Scalar matrix computings.
The content of the invention
For computationally intensive, the low problem of computational efficiency for overcoming distinct element method to exist.The present invention seeks to, propose a kind of
Suitable for the fast discrete unit numerical computation method of GPU Scalar matrix computings.By setting up specific adjacent particles matrix, realize pure
The GPU computings of matrix.The method is all calculated and can realized by matrix operationss, can effectively improve the calculating effect of distinct element method
Rate, realizes large-scale granular discrete-element simulation, and can be further used for ultra-large granular discrete-element parallel computation, realizes
The Fast simulation of complicated geotechnical engineering problems.
The present invention is to solve the technical scheme that the quick computational problem of discrete element is proposed, a kind of suitable for GPU Scalar matrix fortune
The fast discrete unit numerical computation method of calculation.Its step includes：
(1) adjacent particles matrix and granular discrete-element Mathematical Model of heaped-up are set up；By granule by 1 to m numbering, it would be possible to each
The adjacent particles numbering of granule (centrophyten) contact is stored in the corresponding line of adjacent particles matrix Pn, i.e. the rows of matrix Pn i-th
For the adjacent particles numbering of the centrophyten of i, (each granule can be centrophyten to record number, and it is relative to neighbouring
Grain, i.e., these adjacent particles are all near Current central granule)；The m+1 void granule numbering fillings of matrix different rows difference in length；
Each granule can be centrophyten, and it is relative to adjacent particles；I.e. these adjacent particles are all current
Near " centrophyten "；I.e. the adjacent particles of the i-th rows of Pn record number i granule are numbered.
(2) realize that Scalar matrix iterates to calculate numerical density；Based on adjacent particles matrix, adjacent particles coordinate and attribute are turned
Chemical conversion m*n matrix forms corresponding with adjacent particles matrix.In discrete element interative computation, at the beginning of obtaining granule by matrix calculus
Step stress matrix F n0(matrix size m*n).
(3) Force Calculation result is filtered using contact relation matrix, completes iterative calculation.According to factors such as stress
Contact relation Boolean matrix Bc is calculated, using Bc Fn is filtered out0In actual loading unit, obtain granule actual loading matrix,
Calculate and make a concerted effort and complete granule motion simulation.
Step 10 builds initial granular discrete-element Mathematical Model of heaped-up, and granule is by 1 open numbering to m；Numbered according to granule, will
The geometric parameter and mechanics parameter of granule is stored in each array；
Step 11 is obtained the adjacent particles of each granule by conventional method, may contact granule, and it is stored line by line
In adjacent particles matrix Pn；The adjacent particles numbering of each granule is stored in adjacent particles matrix corresponding line, matrix size
For m*n, i.e. m rows n row；
Particle habit is changed into m*n forms corresponding with Pn matrixes by step 12, and carries out Scalar matrix discrete element numerical computations
(being specifically shown in step 20-26), obtains preliminary stress matrix F n of granule0(m*n)；In preliminary stress matrix and adjacent particles matrix
Granule numbering is corresponded, and have recorded the adjacent intergranular primary action power of granule；
Step 13 records the connection of granule and contact condition using contact relation Boolean matrix Bc；Contact relation matrix with
Adjacent particles matrix is corresponded, and whether contact relation matrix record adjacent particles contact with centrophyten, i.e., with the presence or absence of power
Effect.Calculate granule contact relation Boolean matrix Bc (being specifically shown in step 30-36), due in adjacent particles matrix only part
Grain has contact relation with centrophyten, by preliminary stress matrix F n0Bc is multiplied by item by item obtains granule actual loading matrix F n (m*n
Matrix).
Step 14 is calculated with joint efforts according to discrete element granule and motion calculation carries out granule motion simulation.According to actual loading square
Battle array is asked for making a concerted effort suffered by each granule, and calculates granule motion by classical mechanics method.(being specifically shown in step 30-36)
Step 15 realizes the dynamic analog of granular discrete-element by continuous iteration.
Step 16 obtains numerical simulation final result.
Step 20-26 describes the preliminary stress matrix procedures of step 11-12 Scalar matrix iterative calculation granule in detail：
Step 20 is input into the adjacent particles array of each granule, and the attribute array of all granules, including radius (R),
Coordinate (P [X, Y, Z]), positive stiffness factor (Kn) and quality (M), array size is m*1, i.e. m rows 1 and arranges.
Step 21 builds adjacent particles matrix Pn.By number order is line by line stored in the adjacent particles array of granule neighbouring
In granule matrix.Granule matrix Pn is m*n matrixes, that is, have m rows and n row, and the i-th line number value represents neighbouring numbered as i granules
Grain numbering.Because the adjacent particles number of each granule is different, the width n of Pn matrixes is true by the adjacent particles array row most grown
It is fixed.In Pn there is dummy cell in the shorter row of length, be filled with empty granule numbering m+1；
Step 22 increases m+1 void particle habit values at attribute array end, and calculates the corresponding attribute square of adjacent particles matrix
Battle array (m*n).Due to there is m+1 elements in Pn, need finally to increase a property value for numbering m+1 void granule in attribute array,
To ensure that matrix algorithm correctly runs；Pn is updated in attribute array and obtains adjacent particles attribute matrix.
Step 23 obtains centrophyten poor with the x coordinate of adjacent particles according to granule coordinates matrix and adjacent particles matrix
Value matrix dPnX, dPnY, dPnZ, its size is m*n.
Step 24 is calculated distance matrix dPn between centrophyten and adjacent particles：
DPn=sqrt (dPnX.^2+dPnY.^2+dPnZ.^2)；(by gauged distance formula)
Step 25 further obtains relative displacement Sn between granule：
Sn=dPn- (R*ones (1, n)+R ' are (Pn))；
Wherein R*ones (1, n) it is particle radius array (m*1)；R ' is (Pn) the radius array that increases empty granule；The right
Section 2 R*ones (1, n)+R ' (Pn) centered on granule and adjacent particles radius sum.One binomial difference is particle surface
Distance, that is, obtain relative displacement matrix (size m*n) between granule.
Step 26 calculates intergranular preliminary positive force (Fn0), i.e., positive stiffness factor matrix (Kn) is relative and between granule
The product item by item of transposed matrix (Sn)：
Fn0=Kn.*Sn；
In formula, Kn and Sn is m*n matrixes, and element is multiplied item by item in the two matrix, obtains size and just makees for the preliminary of m*n
Firmly matrix F n0。
Step 30-36 describes step 13-14 contact relation Boolean matrix in detail and calculates actual loading and discrete element calculating stream
Journey：
Step 30 is input into preliminary positive force matrix F n of granule0With granule connection Boolean matrix Bb (being m*n matrixes).
Whether Bb is corresponded with adjacent particles matrix, be have recorded and whether have between granule and adjacent particles connecting forces, i.e., allow pulling force to make
With.
Step 31 is according to Fn0Granule connection matrix Bb is updated, when pulling force exceedes particular value (stretching resistance matrix value) between granule
When, connection fracture between granule.Realized by matrix Boolean calculation：
Bb=(Bb＆Fn0<Fnmax), Matlab matrixes instruction, similarly hereinafter；
Wherein FnmaxFor stretching resistance matrix (m*n) between granule.
Step 32 is according to Fn0Calculate and press between granule active force (pressure is negative) Boolean matrix Bp：
Bp=(Fn0<0)；
Step 33 calculates interparticle contact relation Boolean matrix Bc, that is, to indicate and whether have actual force between granule.Granule
Between actual force include pressure and pulling force, therefore, Bc to granule connection matrix Bb and pressure active force matrix Bp by taking union
Obtain：
Bc=(Bb | Bp)；
Step 34 is according to Fn0Calculate actual loading Fn between granule.In adjacent particles matrix and not all adjacent particles all with
Present granule has contact relation, and by the screening of contact relation Boolean matrix actual forward power is obtained：
Fn=Fn0.*Bc；
In contact relation matrix, there is contact to be 1, contactless is 0.
Fn is decomposed into FnX, FnY, FnZ vector form by step 35.With reference to coordinate difference between the granule obtained in step 23
Matrix dPnX, dPnY, dPnZ (i.e. positive force direction), can be decomposed into coordinate vector form D nX, FnY, FnZ by Fn.Granule institute
By being with joint efforts adjacent particles to its active force sum, i.e., take on actual loading matrix line direction and：
FX=sum (FnX, 2)；Matlab matrixes are instructed, similarly hereinafter；
FY=sum (FnY, 2)；
FZ=sum (FnZ, 2)；
Numerical density is F=[FX, FY, FZ] (m*3 matrixes).
Step 36 calculates component motion of the granule in xyz directions according to classical mechanics.
Distinct element method is a kind of Computational Mechanicses method, similar with FInite Element.It resolves into material a series of with spy
Fixed discrete element granule, and intergranular stress and motion are calculated by Newtonian mechanics, so as to the deformation of simulation material and broken
It is bad.
Step 20-26 is the refinement of step 11-12, and step 30-36 is the refinement of step 13-14.Fig. 1 is ensemble stream
Journey.
The invention has the beneficial effects as follows：Realize that all discrete elements are calculated to complete with Scalar matrix form of calculation, it is fast to realize
Fast GPU computings, significantly improve the efficiency of discrete element calculating.This method avoids needing constantly in conventional discrete unit computational methods
It is circulated operation to calculate contact, stress and motion of each granule etc., but by based on adjacent particles matrix and contact
Relation Boolean matrix realizes distinct element method Scalar matrix rapid computations.Can be reached by using these methods calculating, discrete element
The speed of general tens times of business software, is capable of achieving the dynamic analog of million granules on single personal computer.And based on Scalar matrix
Discrete element numerical computations, can easily carry out calculating segmentation, be conducive to the parallel computation of further ultra-large discrete element.
Actual loading matrix by obtaining granule of the invention, calculates finally by conventional discrete unit method and completes granule dynamic
Simulation.In numerical procedure, all properties data are recorded in corresponding matrix based on adjacent particles matrix, and are adopted
Use Scalar matrix computing.Therefore, the method is applied to GPU (graphics processing unit) high speed matrix operationss, significantly improves discrete
The speed of first numerical method, can realize the dynamic analog of million granules on single personal computer.
Description of the drawings
Fig. 1 fast discretes unit's numerical computations population structure and operation principle flow chart；
The computing of Fig. 2 Scalar matrix obtains the preliminary stress matrix procedures figure of granule；
Fig. 3 contact relation matrix calculus actual loadings and distinct element method calculation flow chart.
Specific embodiment
Fig. 1 is the population structure and operation principle that this method is implemented.It is characterized in combining adjacent particles matrix Pn and contact
Relation Boolean matrix Bc realizes that Scalar matrix discrete element is calculated.
Step 10 builds initial granular discrete-element Mathematical Model of heaped-up, and granule is by 1 open numbering to m.According to numbering, by granule
Geometric parameter and mechanics parameter be stored in each array.
Step 11 obtains the adjacent particles (may contact granule) of each granule by conventional method, and it is stored line by line
In adjacent particles matrix Pn.Accompanying drawing 1a gives adjacent particles matrix, and the adjacent particles numbering of each granule is stored in matrix
In corresponding line, matrix size is m*n, i.e. m rows n row.
Particle habit is changed into m*n forms corresponding with Pn matrixes by step 12, and carries out Scalar matrix discrete element numerical computations
(being specifically shown in step 20-26), obtains preliminary stress matrix F n of granule0(m*n).In preliminary stress matrix and adjacent particles matrix
Granule numbering is corresponded, and have recorded the adjacent intergranular primary action power of granule.
Step 13 calculates granule contact relation Boolean matrix Bc, and obtains granule actual loading matrix F n.Due to neighbouring
Only partial particulate has contact relation with centrophyten in grain matrix, and the company of granule must be recorded using contact relation Boolean matrix Bc
Connect and contact condition.Such as accompanying drawing 1b, contact relation matrix is corresponded with adjacent particles matrix, and record adjacent particles whether with
Centrophyten is contacted, i.e., with the presence or absence of the effect of power.Further, by preliminary stress matrix F n0Bc is multiplied by item by item obtains granule reality
Border stress matrix F n (m*n matrixes).
Step 14 is made a concerted effort and motion calculation according to discrete element granule.Close according to suffered by actual loading matrix asks for each granule
Power, and granule motion is calculated by classical mechanics method.
Step 15 realizes the dynamic analog of granular discrete-element by continuous iteration.
Step 16 obtains numerical simulation final result.
Fig. 2 is the preliminary stress matrix procedures figure of Scalar matrix iterative calculation granule.It is characterized according to adjacent particles matrix structure
Corresponding particle habit matrix is built, and granule and adjacent particles intermolecular forces matrix are obtained by Scalar matrix computing.Step 20-26
For the refinement explanation of step 11-12.
Step 20 is input into the adjacent particles array of each granule, and the attribute array of all granules, such as radius (R),
Coordinate (P [X, Y, Z]), positive stiffness factor (Kn) and quality (M) etc., array size is m*1, i.e. m rows 1 and arranges.
Step 21 builds adjacent particles matrix Pn.By number order is line by line stored in the adjacent particles array of granule neighbouring
In granule matrix.Granule matrix Pn is m*n matrixes, that is, have m rows and n row, and the i-th line number value represents neighbouring numbered as i granules
Grain numbering.Because the adjacent particles number of each granule is different, the width n of Pn matrixes is true by the adjacent particles array row most grown
It is fixed.In Pn there is dummy cell in the shorter row of length, be filled with empty granule numbering m+1.Such as accompanying drawing 2a, granule 2 (the second row) is gathered around
There are most adjacent particles, its last element is 5.And other row relevant positions are then filled with numbering m+1.Because granule is compiled
Number by 1 to m, in further matrix operationss, numbering m+1 can be distinguished.
Step 22 increases m+1 void particle habit values at attribute array end, and calculates the corresponding attribute square of adjacent particles matrix
Battle array (m*n).Due to there is m+1 elements in Pn, need finally to increase a property value for numbering m+1 void granule in attribute array,
To ensure that matrix algorithm correctly runs.Further, Pn is updated in attribute array and obtains attribute matrix.Such as granule x coordinate
Array is PX ' (m+1 element), then corresponding adjacent particles x coordinate matrix be PX ' (Pn) (instruction of Matlab matrixes).
Step 23 obtains centrophyten poor with the x coordinate of adjacent particles according to granule coordinates matrix and adjacent particles matrix
Value matrix dPnX, dPnY, dPnZ, its size is m*n.Specifically can be obtained by following formula (by taking dPnX as an example)：
DPnX=PX ' (Pn)-PX*ones (1, n)；(Matlab matrixes are instructed, similarly hereinafter)
Formula the right Section 1 obtains adjacent particles coordinate in the x direction；Section 2 obtains centrophyten in the x direction
Coordinate.PX is m*1 matrixes, and ones (1, n) obtain being worth the 1*n matrixes for being all 1, the two multiplications obtains m*n matrixes, its often behavior
The x coordinate value of centrophyten；One binomial difference obtains adjacent particles and centrophyten matrix of differences in the x direction.Equally
Ground, can obtain adjacent particles and centrophyten matrix of differences dPnY in the y and z directions and dPnZ.
Step 24 is calculated distance matrix dPn between centrophyten and adjacent particles：
DPn=sqrt (dPnX.^2+dPnY.^2+dPnZ.^2)；
Step 25 further obtains relative displacement Sn between granule：
Sn=dPn- (R*ones (1, n)+R ' are (Pn))；
Wherein R is particle radius array (m*1)；R ' is the radius array for increasing empty granule；The right Section 2 centered on
Grain and the radius sum of adjacent particles.One binomial difference is the distance of particle surface, that is, obtain relative displacement matrix between granule (big
Little m*n).
Step 26 calculates intergranular preliminary positive force (Fn0), i.e., positive stiffness factor matrix (Kn) is relative and between granule
The product item by item of transposed matrix (Sn)：
Fn0=Kn.*Sn；
In formula, Kn and Sn is m*n matrixes, and element is multiplied item by item in the two matrix, obtains size and just makees for the preliminary of m*n
Firmly matrix F n0.Accompanying drawing 2b is preliminary positive force matrix schematic diagram, and it is corresponded with adjacent particles matrix, and be have recorded
Positive acting power between granule and adjacent particles.Preliminary positive force matrix is intergranular possible active force, due to neighbouring
Unit in granule matrix may not contact (i.e. without pressure and tension force effect) with centrophyten, need further using contact
Relation Boolean matrix is filtering out actual active force.
Fig. 3 is contact relation matrix calculus actual loading and discrete element calculation flow chart.It is characterized in setting up contact relation
Boolean matrix obtains intergranular actual force filtering preliminary positive force matrix, and completes discrete element numerical computations.Step
Rapid 30-36 is the refinement explanation of step 13-14.
Step 30 is input into preliminary positive force matrix F n of granule0With granule connection Boolean matrix Bb (being m*n matrixes).
Whether Bb is corresponded with adjacent particles matrix, be have recorded and whether have between granule and adjacent particles connecting forces, i.e., allow pulling force to make
With.
Step 31 is according to Fn0Granule connection matrix Bb is updated, when pulling force exceedes particular value between granule, is connected between granule disconnected
Split.Can be realized by matrix Boolean calculation：
Bb=(Bb＆Fn0<Fnmax) (Matlab matrixes are instructed, similarly hereinafter)
Wherein FnmaxFor stretching resistance matrix (m*n) between granule.
Step 32 is according to Fn0Calculate and press between granule active force (pressure is negative) Boolean matrix Bp：
Bp=(Fn0<0)；
Step 33 calculates interparticle contact relation Boolean matrix Bc, that is, to indicate and whether have actual force between granule.Granule
Between actual force include pressure and pulling force, therefore, Bc to granule connection matrix Bb and pressure active force matrix Bp by taking union
Obtain：
Bc=(Bb | Bp)；
Accompanying drawing 3a is Bc matrix schematic diagrams, and its size is m*n, is corresponded with adjacent particles matrix element, and be have recorded
With the presence or absence of contact and the effect of power between granule.
Step 34 is according to Fn0Calculate actual loading Fn between granule.In adjacent particles matrix and not all adjacent particles all with
Present granule has contact relation, and by the screening of contact relation Boolean matrix actual forward power is obtained：
Fn=Fn0.*Bc；
In contact relation matrix, there is contact to be 1, contactless is 0.Preliminary positive force Fn in accompanying drawing 2b0With accompanying drawing
Granule contact Boolean matrix in 3a is multiplied item by item, then obtain actual positive force matrix F n between the granule in accompanying drawing 3b.
Fn is decomposed into FnX, FnY, FnZ vector form by step 35.With reference to coordinate difference between the granule obtained in step 23
Matrix dPnX, dPnY, dPnZ (i.e. positive force direction), can be decomposed into coordinate vector form D nX, FnY, FnZ by Fn.Granule institute
By being with joint efforts adjacent particles to its active force sum, i.e., take on actual loading matrix line direction and：
FX=sum (FnX, 2)；(Matlab matrixes are instructed, similarly hereinafter)
FY=sum (FnY, 2)；
FZ=sum (FnZ, 2)；
Numerical density is F=[FX, FY, FZ] (m*3 matrixes).
Step 36 calculates component motion of the granule in xyz directions according to classical mechanics.As in the x direction component is calculated
For：
AX=FX./M；
VX=VX+AX*dT；
PX=PX+VX*dT；
M is granular mass array in formula, and AX is x directional accelerations, and VX is x directions speed, and PX is granule x directions coordinate,
DT is iteration time step-length.
Example above gives the method for realizing that discrete element positive force is calculated with interative computation Scalar matrix computing.Similarly,
Tangential force and torque etc. can also be calculated by similar approach between granule, so as to realize the Scalar matrix computing of granular discrete-element.
Claims (2)
1., suitable for the fast discrete unit numerical computation method of GPU Scalar matrix computings, it is characterized in that step includes：
1）Set up adjacent particles matrix and granular discrete-element Mathematical Model of heaped-up；By granule by 1 to m numbering, it would be possible to contact with granule
Adjacent particles numbering be stored in the corresponding line of adjacent particles matrix Pn, row difference in length with m+1 void granule numbering filling；
2）Realize that Scalar matrix iterates to calculate numerical density；Based on adjacent particles matrix, adjacent particles coordinate and attribute are changed into
M*n matrix forms corresponding with adjacent particles matrix；In discrete element interative computation, obtain granule by matrix calculus and tentatively receive
Torque battle array Fn0、Matrix size m*n；
3）Force Calculation result is filtered using contact relation matrix, completes iterative calculation；Calculated according to factors such as stress
Contact relation Boolean matrix Bc, using Bc Fn is filtered out0In actual loading unit, obtain granule actual loading matrix F n, count
Calculate and make a concerted effort and complete granule motion simulation；
Step 1）The step of setting up adjacent particles matrix and granular discrete-element Mathematical Model of heaped-up：
Step 10 builds initial granular discrete-element Mathematical Model of heaped-up, and granule is by 1 open numbering to m；Numbered according to granule, by granule
Geometric parameter and mechanics parameter be stored in each array；
Step 11 is obtained the adjacent particles of each granule by conventional method, may contact granule, and it is stored in line by line neighbour
In nearly granule matrix Pn；The adjacent particles numbering of each granule is stored in adjacent particles matrix corresponding line, and matrix size is m*
N, i.e. m rows n are arranged；
Particle habit is changed into m*n forms corresponding with Pn matrixes by step 12, and carries out Scalar matrix discrete element numerical computations, is obtained
Preliminary stress matrix F n of granule0, m*n ranks；Preliminary stress matrix is corresponded with granule numbering in adjacent particles matrix, record
Granule adjacent intergranular primary action power；
Step 13 records the connection of granule and contact condition using contact relation Boolean matrix Bc；Contact relation matrix with it is neighbouring
Granule matrix is corresponded, and whether contact relation matrix record adjacent particles contact with centrophyten, i.e., with the presence or absence of the work of power
With；Granule contact relation Boolean matrix Bc is calculated, because only partial particulate contacts pass with centrophyten in adjacent particles matrix
System, by preliminary stress matrix F n0Bc is multiplied by item by item obtains granule actual loading matrix F n, m*n matrix；
Step 14 is calculated with joint efforts according to discrete element granule and motion calculation carries out granule motion simulation；According to actual loading Matrix Calculating
Take each granule suffered with joint efforts, and granule motion is calculated by classical mechanics method；It is specifically shown in step（30）-（36）；
Step（30）Preliminary positive force matrix F n of input granule0Connect Boolean matrix Bb, be m*n matrixes with granule；Bb with
Adjacent particles matrix is corresponded, and be have recorded and whether have between granule and adjacent particles connecting forces, i.e., whether allow pulling force effect；
Step（31）According to Fn0Granule connection matrix Bb is updated, when pulling force is stretching resistance matrix value more than particular value between granule,
Connection fracture between granule；Realized by matrix Boolean calculation：
Bb=(Bb&Fn0<Fnmax), Matlab matrixes instruction, similarly hereinafter；
Wherein FnmaxFor stretching resistance matrix m*n between granule；
Step（32）According to Fn0It is negative Boolean matrix Bp to calculate and press between granule active force, pressure：
Bp=(Fn0< 0)；
Step（33）Interparticle contact relation Boolean matrix Bc is calculated, that is, to be indicated and whether have actual force between granule；Between granule
Actual force includes pressure and pulling force, therefore, Bc to granule connection matrix Bb and pressure active force matrix Bp by taking union and obtaining
Arrive：
Bc=(Bb|Bp)；
Step（34）According to Fn0Calculate actual loading Fn between granule；In adjacent particles matrix and not all adjacent particles all with currently
Granule has contact relation, and by the screening of contact relation Boolean matrix actual forward power is obtained：
Fn=Fn0*Bc；
In contact relation matrix, there is contact to be 1, contactless is 0；
Step（35）Fn is decomposed into into FnX, FnY, FnZ vector form；With reference to coordinate difference between the granule obtained in step 23
Matrix dPnX, dPnY, dPnZ are positive force direction, and Fn is decomposed into into coordinate vector form D nX, FnY, FnZ；Granule institute
By being with joint efforts adjacent particles to its active force sum, i.e., take on actual loading matrix line direction and：
FX=sum(FnX, 2);Matlab matrixes are instructed, similarly hereinafter；
FY=sum(FnY, 2);
FZ=sum(FnZ, 2);
Numerical density is F=[FX, FY, FZ], m*3 matrixes；
Step（36）Component motion of the granule in xyz directions is calculated according to classical mechanics；
Step 15 realizes the dynamic analog of granular discrete-element by continuous iteration；
Step 16 obtains numerical simulation final result.
2. the fast discrete unit numerical computation method suitable for GPU Scalar matrix computings according to claim 1, is characterized in that
The preliminary stress matrix procedures of Scalar matrix iterative calculation granule；Corresponding particle habit matrix is built according to adjacent particles matrix, and
Granule and adjacent particles intermolecular forces matrix are obtained by Scalar matrix computing；
Step 20 is input into the adjacent particles array of each granule, and the attribute array of all granules, including radius R, coordinate P
[X, Y, Z], positive stiffness factor Kn and mass M, array size is m*1, i.e. m rows 1 and arranges；
Step 21 builds adjacent particles matrix Pn；By number the adjacent particles array of granule is stored in line by line adjacent particles by order
In matrix；Granule matrix Pn is m*n matrixes, that is, have m rows and n row, and the i-th line number value represents the adjacent particles volume numbered as i granules
Number；Because the adjacent particles number of each granule is different, the width n of Pn matrixes is determined by the adjacent particles array row most grown；Pn
Be present dummy cell in the shorter row of middle length, filled with empty granule numbering m+1；Step 22 increases empty of m+1 at attribute array end
Grain property value, and calculate the corresponding attribute matrix m*n of adjacent particles matrix；Due to there is m+1 elements in Pn, need in attribute number
Group finally increases a property value for numbering m+1 void granule, to ensure that matrix algorithm correctly runs；Pn is updated to into attribute array
In obtain adjacent particles attribute matrix；Step 23 obtains centrophyten with neighbour according to granule coordinates matrix and adjacent particles matrix
X coordinate matrix of differences dPnX of nearly granule, dPnY, dPnZ, its size is m*n；Step 24 be calculated centrophyten with
Distance matrix dPn between adjacent particles：
dPn=sqrt(dPnX2+ dPnY2+dPnZ2);
Step 25 further obtains relative displacement Sn between granule：
Sn=dPn-(R*ones(1,n)+R’(Pn));
Wherein R* ones (1, n) it is particle radius array m*1；R ' is (Pn) the radius array that increases empty granule；The right second
R*ones (1, n)+R ' (Pn) centered on granule and adjacent particles radius sum；One binomial difference for particle surface away from
From it is m*n to obtain relative displacement matrix, size between granule；
Step 26 calculates intergranular preliminary positive force Fn0, i.e., relative displacement matrix Sn between positive stiffness factor matrix K n and granule
Product item by item：
Fn0=Kn*Sn;
In formula, Kn and Sn is m*n matrixes, and element is multiplied item by item in the two matrix, obtains the preliminary positive force that size is m*n
Matrix F n0。
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 CN104239277A (en) | 2014-12-24 |
| CN104239277B true 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) |
Families Citing this family (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105866855B (en) * | 2015-01-22 | 2018-04-06 | 中国石油天然气股份有限公司 | A kind of analysis method of geologic-tectonic evolution and deformation process |
| 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 |
| CN109977505B (en) * | 2019-03-13 | 2024-02-09 | 南京大学 | Discrete element pore system quick searching method based on GPU matrix calculation |
| CN110781448B (en) * | 2019-11-04 | 2024-03-12 | 南京大学 | Discrete element neighbor searching and solving method and system based on complete matrix calculation |
| CN110929430B (en) * | 2019-12-31 | 2023-04-11 | 浙江交通职业技术学院 | Modulus calculation method for randomly filling silicon carbide particles in diesel engine particle catcher |
| CN113656656B (en) * | 2021-07-21 | 2024-01-26 | 南京南力科技有限公司 | Efficient neighbor retrieval method and system for wide-grading discrete element particle system |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103984829A (en) * | 2014-05-22 | 2014-08-13 | 李笑宇 | Discrete element method based method for improving particle discrete contact detection efficiency |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103229177B (en) * | 2010-09-15 | 2016-08-24 | 联邦科学与工业研究组织 | discrete element method |
-
2014
- 2014-08-15 CN CN201410405053.1A patent/CN104239277B/en active Active
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| 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 |
|---|
| A fast scalable implementation of the two-dimensional triangular Discrete Element Method on a GPU platform;Lin Zhang etc.;《Advances in Engineering Software》;20121120;第70–80页 * |
| 基于GPU的颗粒离散元计算方法研究;常新正;《中国优秀硕士学位论文全文数据库 信息科技辑》;20130915(第9期);摘要、第3章3.1，3.2，第4章4.2，4.3，4.4 * |
| 基于离散元理论的振动筛分数值模拟程序开发;孙 鹏等;《中国制造业信息化》;20120331;第41卷(第5期);第71-73页 * |
| 离散元法研究的评述;刘凯欣等;《力学进展》;20031125;第33卷(第4期);第483-490页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN104239277A (en) | 2014-12-24 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN104239277B (en) | Rapid discrete-element numerical-value calculating method applicable to GPU (Graphics Processing Unit) scalar-matrix operation | |
| Recuero et al. | A high-fidelity approach for vehicle mobility simulation: Nonlinear finite element tires operating on granular material | |
| Yang et al. | Intra: 3d intracranial aneurysm dataset for deep learning | |
| Raji et al. | Model for the deformation in agricultural and food particulate materials under bulk compressive loading using discrete element method. I: Theory, model development and validation | |
| WO2018180263A1 (en) | Information processing device, information processing method, and computer-readable storage medium | |
| JP2015530636A (en) | Particle flow simulation system and method | |
| WO2008021652A9 (en) | Computer simulation of physical processes | |
| JP5644872B2 (en) | Simulation apparatus, simulation method, and program | |
| US20140288906A1 (en) | Analysis device, analysis method, analysis program, and recording medium | |
| Singh et al. | Accurate and efficient solution of bivariate population balance equations using unstructured grids | |
| Chen et al. | Hrnet: Hamiltonian rescaling network for image downscaling | |
| Sierakowski | GPU-centric resolved-particle disperse two-phase flow simulation using the Physalis method | |
| Crous et al. | The potential of graphical processing units to solve hydraulic network equations | |
| CN109753682A (en) | A kind of finite element matrix analogy method based on the end GPU | |
| Recuero et al. | ANCF continuum-based soil plasticity for wheeled vehicle off-road mobility | |
| Va et al. | Parallel Cloth Simulation Using OpenGL Shading Language. | |
| Carrillo et al. | Design of a large scale discrete element soil model for high performance computing systems | |
| Serban et al. | An integrated framework for high-performance, high-fidelity simulation of ground vehicle-tyre-terrain interaction | |
| Cetinaslan | Position‐based simulation of elastic models on the GPU with energy aware gauss‐seidel algorithm | |
| CN110909473B (en) | Dynamic fluid-solid interaction simulation method based on SP H and shape matching mixed model | |
| De Cola et al. | Concurrent adaptive mass-conserving comminution of granular materials using rigid elements | |
| Zhang et al. | Correlation mechanism of friction behavior and topological properties of the contact network during powder compaction | |
| Luo | Dealing with extremely large deformation by nearest-nodes FEM with algorithm for updating element connectivity | |
| Alart et al. | Parallel computational strategies for multicontact problems: Applications to cellular and granular media | |
| Hou et al. | Simulation of heat conduction in fluids on GPU with particle method. |
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 |