CN109408973A - One kind is based on apart from potential function two-dimensional variable shape convex polygon block discrete element method - Google Patents
One kind is based on apart from potential function two-dimensional variable shape convex polygon block discrete element method Download PDFInfo
- Publication number
- CN109408973A CN109408973A CN201811276841.XA CN201811276841A CN109408973A CN 109408973 A CN109408973 A CN 109408973A CN 201811276841 A CN201811276841 A CN 201811276841A CN 109408973 A CN109408973 A CN 109408973A
- Authority
- CN
- China
- Prior art keywords
- finite elements
- discrete
- time
- apart
- displacement
- 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
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
The invention discloses one kind based on apart from potential function two-dimensional variable shape convex polygon block discrete element method, comprising: establishes deformable discrete unit system, determines the time step of system；The contact force of each osculating element and object element is calculated in current time step, contact force is converted to the equivalent node force vector of load；The dynamic governing equations that system incremental form is established by the equivalent node force vector of load obtain the displacement of finite elements；According to the displacement of finite elements, each finite elements coordinate information is updated；It computes repeatedly until all time steps have been calculated.This method realizes the deformable of discrete element system, realize the contact detection and contact force computational problem of different size, form unit, reduce the quantity of actual division unit, improve computational efficiency, so that discrete element analysis reacts the stress and strain situation of bulk inner more accurately, it can be used for simulating more engineering problems.
Description
Technical field
The present invention relates to one kind based on apart from potential function two-dimensional variable shape convex polygon block discrete element method, belongs to variable
Shape discrete element technical field.
Background technique
Discrete element elements method is the method for numerical simulation for being specifically used to solve the problems, such as discontinuous media, and this method can be accurate
Capture block system separation, sliding rupture, the Discontinuous Deformations characteristics such as rotation of toppling.And deformable discrete element can be pressed
Contracting, separation or sliding.Discrete element can more accurately simulate the large deformation feature in rock mass, therefore discrete element method is ground in theory
Study carefully and obtains tremendous development with numerous areas such as engineer applications.
Finite discrete elements method is proposed by Britain professor A.MUNJIZA at present, it is equal by the way that research object is divided into size
Even tetrahedron block unit, and establish potential function based on unit centroid and define with the contact force between this computing unit.
Professor A.MUNJIZA proposes the deformable discrete element based on potential function method, in conjunction with discrete element method and limited list
First method solves deformable discrete element problem.Munjiza utilizes explicit scheme solving finite element, and it is non-thread to avoid solving finite element
The iterative process of property equation group.But there are still some problems: the uniform tetrahedron element of application size, one side model and reality
Situation is not inconsistent, on the other hand in practical application, the unit size and simplest unit form of homogenization can greatly increase
The quantity of Rigid Body Element is divided, computational efficiency is reduced.Solve these problems based on the discrete element method apart from potential function, still
There is no the deformable of consideration discrete element, therefore not fully meet engineering reality.
Summary of the invention
The technical problem to be solved by the present invention is in the prior art based on apart from any convex polygon block of potential function from
It is non-deformable to dissipate elements method, provides one kind based on apart from potential function two-dimensional variable shape convex polygon block discrete element method, changes
The deficiency of kind discrete element method makes numerical simulation be more in line with reality.
The present invention uses following technical scheme to solve above-mentioned technical problem:
The present invention provides one kind to be based on apart from potential function two-dimensional variable shape convex polygon block discrete element method, including with
Lower step:
Step 1: establishing deformable discrete unit system, and the system includes multiple discrete units and by discrete unit
The finite elements formed after subdivision grid；
Step 2: the time step Δ t of system is determined；
Step 3: contact detection is carried out to one layer of the grid cell in discrete unit periphery；
Step 4: step 3 is calculated to the resultant force shape of the contact force acted on osculating element and object element
Function is converted to the equivalent node force vector of grid cell current time load
Step 5: by the equivalent node force vector for the load that step 4 is calculated, the power of system incremental form is established
Governing equation and solve obtain the displacement of subsequent time t+ Δ t finite elements；
Step 6: according to the displacement of step 5 finite elements, updating the information of the coordinate on each grid cell vertex, completes
The calculating at current time；
Step 6: it repeats step 2 to step 6 and calculates, until all time steps have been calculated.
Further, in step 2, the time step Δ t must meet:
Δ t=min (Δ tD,Δts),
Δts≤ L/C,
Wherein, Δ tDFor the calculating time step of discrete unit；ξ be discrete unit damping ratio andM is
Discrete unit block quality, c are damped coefficient, and k is stiffness coefficient, Δ tsFor the time step of finite elements, L is all limited
The minimum side length of unit, the value of C are 10000.
Further, it in step 4, is used when calculating the contact force of current time each osculating element and object element
Contact force algorithm based on the discrete element method apart from potential function.
Further, in step 5, the motion control equation of system incremental form is established, and solve and obtain subsequent time
The displacement of t+ Δ t finite elements, specifically includes:
Step 5.1: establishing reference frame, be configured as reference configuration before selection deformation；
Step 5.2: it is discrete using broad sense Newmark method progress time-domain, predict current time mechanical quantity；
Step 5.3: establishing the dynamic governing equations of system incremental formSolution is worked as
Acceleration increment of the preceding moment t to subsequent time t+ Δ tWherein, M is the mass matrix of finite elements, and D is finite elements
Damping matrix, K is the stiffness matrix of finite elements,It is the acceleration increment of finite elements,It is the speed of finite elements
Increment, Δ u are the displacement increment of finite elements, the equivalent node force vector R of grid cell current time load；
Step 5.4: it is discrete by broad sense Newmark method progress time-domain again, calculate each finite elements subsequent time t+
It is every that subsequent time t+ Δ t is calculated since grid cell is consistent with the node coordinate of finite elements in the displacement of Δ t
The displacement of a grid cell；
Further, in step 6, the coordinate formula on the vertex of each grid cell is updated are as follows:
X (t+ Δ t)=x (t)+(r (t+ Δ t))x
Y (t+ Δ t)=y (t)+(r (t+ Δ t))y
Wherein, ((t+ Δ t) is on coordinate x (t), the y (t) on the vertex of current time step grid cell are to x by t+ Δ t), y
Coordinate (the r (t+ Δ t)) on the vertex of one time step Rigid Body Elementx、(r(t+Δt))yRespectively the displacement of grid cell is in x, the side y
To component.
Present invention technical effect achieved: this method realizes the deformable of discrete element system, so that discrete element analysis
The more accurate stress and strain situation that must react bulk inner, can be used for simulating more engineering problems；It realizes not
The contact detection and contact force computational problem of same size, form unit, reduce the quantity of actual division unit, improve calculating
Efficiency.
Detailed description of the invention
Two discrete units of Fig. 1 the method for the present invention specific embodiment are further split into the schematic diagram of finite elements；
Contact of Fig. 2 the method for the present invention specific embodiment osculating element with object element is overlapped schematic diagram；
Fig. 3~Fig. 6 embodiment rock side slope difference calculates the process schematic that the landslide at moment is destroyed；
State when wherein Fig. 3 is the Rock Slope Stability of specific embodiment；
Fig. 4~Fig. 6 is the motion process that side slope generates landslide along sliding surface in specific embodiment under the influence of by external condition.
Specific embodiment
Technical solution of the present invention is described in further detail with reference to the accompanying drawing:
The present invention provides one kind to be based on apart from potential function two-dimensional variable shape convex polygon block discrete element method, including with
Lower step:
Step 1 establishes deformable discrete unit system, as shown in Figure 1, this system includes multiple discrete units, and
The finite elements that will be formed after discrete unit subdivision grid；Each mesh definition after discrete unit subdivision grid is grid list
Member, the node coordinate of this grid cell and the node coordinate of finite elements are consistent, wherein the parameter of discrete unit include from
The node coordinate of throwaway member, quality, damping ratio, rigidity, the parameter of finite elements include the node coordinate of finite elements, moment of mass
Battle array, damping matrix, stiffness matrix；
Step 2, step-length of fixing time, calculating time step Δ t must meet:
Δ t=min (Δ tD,Δts)
Δts≤L/C
Wherein, Δ tDFor the calculating time step of discrete element；ξ is the damping ratio of system,M is bulk single
First quality, c are damped coefficient, and k is stiffness coefficient, Δ tsFor the time step of finite element, L is the minimum edge of all finite elements
Long, the value of C is 10000.
Step 3, one layer of the grid cell in discrete unit periphery carry out contact and detect to determine object element and osculating element
Method is that one layer of the discrete unit periphery grid cell for generating contact is defined as object element, outside the discrete unit contacted
It encloses one layer of grid cell and is then defined as osculating element.As shown in Figure 2 according to based on apart from any convex polygon block of potential function
The definition of discrete element method calculates current time step and acts on the normal direction contact force of object element and osculating element, tangential contact
Power；
The resultant force shape of the contact force acted on osculating element and object element is calculated in step 3 by step 4
Function is converted to the equivalent node force vector of grid cell current time load
Wherein,It is the equivalent node force vector that grid cell currently carves load load,WithIt is current time respectively
Grid cell physical strength and face power load vector, N is the shape function of grid cell node, V0It is the volume of grid cell, A0t
It is the surface area of current t moment grid cell, A0It is the surface area of grid cell；
Step 5, the power control of the equivalent node force vector solving system incremental form for the load being calculated by step 4
Equation processed obtains the displacement of finite elements；
Step 5.1: establishing reference frame, be configured as reference configuration before selection deformation；
Step 5.2: it is discrete using broad sense Newmark method progress time-domain, predict current time mechanical quantity；
Step 5.3: establishing the dynamic governing equations of system incremental formSolution is worked as
Acceleration increment of the preceding moment t to subsequent time t+ Δ tWherein, M is the mass matrix of finite elements, and D is finite elements
Damping matrix, K is the stiffness matrix of finite elements,It is the acceleration increment of finite elements,It is the speed of finite elements
Increment, Δ u are the displacement increments of finite elements；
Step 5.4: it is discrete by broad sense Newmark method progress time-domain again, calculate each finite elements subsequent time t+
It is every that subsequent time t+ Δ t is calculated since grid cell is consistent with the node coordinate of finite elements in the displacement of Δ t
The displacement of a grid cell；
In other specific embodiments, the dynamic governing equations for the system incremental form that can be proposed through the invention are asked
Solve the speed or acceleration of for example each grid cell t+ time Δt of other variate-values of finite elements.
Step 6 updates the geometry such as the coordinate on each finite elements vertex letter according to the displacement of finite elements in step 5
Breath completes the calculating of current time step；
Update the coordinate formula on the vertex of each grid cell are as follows:
X (t+ Δ t)=x (t)+(r (t+ Δ t))x
Y (t+ Δ t)=y (t)+(r (t+ Δ t))y
Wherein, ((t+ Δ t) is on coordinate x (t), the y (t) on the vertex of current time step grid cell are to x by t+ Δ t), y
Coordinate (the r (t+ Δ t)) on the vertex of one time step Rigid Body Elementx、 (r(t+Δt))yRespectively the displacement of grid cell is in x, y
The component in direction.
Step 7 repeats step 2 to six calculating future time steps, until all time steps have been calculated.
Embodiment:
Certain rock side slope under the influence of by external conditions such as earthquake, rainfalls, may lead since there are weak intercalated layers inside it
Send out the geological disaster on landslide.Using method provided by the invention, Discrete element model is established for rock side slope, as shown in figure 3, will
Discrete unit slip mass and slope body divide finite elements 130 altogether, so that grid cell is also 130.
State when Fig. 3 is Rock Slope Stability.
Fig. 4 to Fig. 6 is side slope generates landslide along sliding surface motion process under the influence of by external condition.When slip mass is discrete
It when unit is started sliding by gravity, is just in contact with slope body, between slip mass discrete unit and slip mass discrete unit, comes down
Contact force is generated due to contacting between body discrete unit and slope body discrete unit, then contact force is sweared by the equivalent node of load
Amount is transformed into finite elements up, so that the displacement of finite elements changes, has updated the node coordinate of grid cell, thus
The node coordinate of discrete unit is updated, and carries out contact detection again, and new contact occurs, and generates new contact force, has
The displacement of limit unit changes again, is recycled with this until movement terminates.Based on provided by the present invention apart from potential function two
Deformable block discrete-time epidemic model rock side slope process of landslides is tieed up, rock side slope can be explicitly described to be influenced by unfavorable load
Under, it is whether safe under load action to can be very good analysis rock side slope for the process destroyed along sliding surface, if it is broken to generate landslide
Form, volume, the scale of accumulation body etc. that bad process and slip mass are formed can intuitively be shown very much.
The above is only a preferred embodiment of the present invention, it is noted that for the ordinary skill people of the art
For member, without departing from the technical principles of the invention, several improvement and deformations can also be made, these improvement and deformations
Also it should be regarded as protection scope of the present invention.
Claims (7)
1. one kind is based on apart from potential function two-dimensional variable shape convex polygon block discrete element method, which is characterized in that including following
Step:
Step 1: establishing deformable discrete unit system, the deformable discrete unit system include several discrete units with
And the finite elements that will be formed after discrete unit subdivision grid；Step 2: the time step Δ t of system is determined；
Step 3: contact detection is carried out to the grid cell of all one layer of discrete unit peripheries:
Step 4: step 3 is calculated to the resultant force shape function of the contact force acted on osculating element and object element
It is converted to the equivalent node force vector of grid cell current time load
Step 5: by the equivalent node force vector for the load that step 4 is calculated, the dynamic Control of system incremental form is established
Equation and solve obtain the displacement of subsequent time t+ Δ t finite elements；
Step 6: according to the displacement of step 5 finite elements, updating the information of the coordinate on each grid cell vertex, completes current
The calculating at moment；
Step 7: it repeats step 2 to step 6 and calculates, until all time steps have been calculated.
2. one kind according to claim 1 is based on apart from potential function two-dimensional variable shape convex polygon block discrete element method,
It is characterized in that, the time step Δ t must meet:
Δ t=min (Δ tD,Δts),
Δts≤ L/C,
Wherein, Δ tDFor the calculating time step of discrete unit；ξ be discrete unit damping ratio andM is discrete list
First block quality, c are damped coefficient, and k is stiffness coefficient, Δ tsFor the time step of finite elements, L is all finite elements
Minimum side length, the value of C are 10000.
3. one kind according to claim 1 is based on apart from potential function two-dimensional variable shape convex polygon block discrete element method,
It is characterized in that in step 3, using based on apart from gesture when calculating the contact force of current time each osculating element and object element
The contact force algorithm of the discrete element method of function.
4. one kind according to claim 1 is based on apart from potential function two-dimensional variable shape convex polygon block discrete element method,
It is characterized in that, the method for contact detection includes:
By discrete unit periphery one layer generate contact grid cell be defined as object element, contacted with the object element from
The grid cell of one layer of throwaway member periphery is then defined as osculating element, and calculates current time each osculating element and target list
The contact force of member.
5. one kind according to claim 1 is based on apart from potential function two-dimensional variable shape convex polygon block discrete element method,
It is characterized in that, establishing the motion control equation of system incremental form, and solve and obtain subsequent time t+ in the step 5
The displacement of Δ t finite elements, specifically includes:
Step 5.1: establishing reference frame, be configured as reference configuration before selection deformation；
Step 5.2: it is discrete using broad sense Newmark method progress time-domain, predict current time mechanical quantity；
Step 5.3: establishing the dynamic governing equations of system incremental formWhen solution obtains current
Carve the acceleration increment Delta ü of t to subsequent time t+ Δ t, wherein M is the mass matrix of finite elements, and D is the resistance of finite elements
Buddhist nun's matrix, K are the stiffness matrix of finite elements, and Δ ü is the acceleration increment of finite elements,It is the speed increasing of finite elements
Amount, Δ u is the displacement increment of finite elements；
Step 5.4: it is discrete by broad sense Newmark method progress time-domain again, calculate each finite elements subsequent time t+ Δ t's
Displacement, since grid cell is consistent with the node coordinate of finite elements, is calculated each grid of subsequent time t+ Δ t
The displacement of unit.
6. one kind according to claim 1 is based on apart from potential function two-dimensional variable shape convex polygon block discrete element method,
It is characterized in that, updating the coordinate formula on the vertex of each finite elements in the step 5 are as follows:
X (t+ Δ t)=x (t)+(r (t+ Δ t))x
Y (t+ Δ t)=y (t)+(r (t+ Δ t))y
Wherein, ((t+ Δ t) is that coordinate x (t), the y (t) on the vertex of current time step grid cell were the upper time to x by t+ Δ t), y
Walk the coordinate (r (t+ Δ t)) on the vertex of Rigid Body Elementx、(r(t+Δt))yRespectively the displacement of grid cell is in x, point in the direction y
Amount.
7. one kind according to claim 1 is based on apart from potential function two-dimensional variable shape convex polygon block discrete element method,
It is characterized in that, the parameter of discrete unit includes the node coordinate of discrete unit, quality, damping ratio, rigidity；The ginseng of finite elements
Number includes node coordinate, mass matrix, damping matrix, the stiffness matrix of finite elements.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811276841.XA CN109408973B (en) | 2018-10-30 | 2018-10-30 | Distance potential function based two-dimensional deformable convex polygon block discrete unit method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811276841.XA CN109408973B (en) | 2018-10-30 | 2018-10-30 | Distance potential function based two-dimensional deformable convex polygon block discrete unit method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109408973A true CN109408973A (en) | 2019-03-01 |
| CN109408973B CN109408973B (en) | 2021-06-08 |
Family
ID=65470307
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201811276841.XA Active CN109408973B (en) | 2018-10-30 | 2018-10-30 | Distance potential function based two-dimensional deformable convex polygon block discrete unit method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109408973B (en) |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6099573A (en) * | 1998-04-17 | 2000-08-08 | Sandia Corporation | Method and apparatus for modeling interactions |
| CN104809301A (en) * | 2015-05-07 | 2015-07-29 | 北京理工大学 | Time domain substructure method for reflecting nonlinear rigid-soft hybrid junction characteristics |
| CN105912852A (en) * | 2016-04-08 | 2016-08-31 | 河海大学 | Arbitrary convex polygon block discrete unit method based on distance potential function |
| CN106529146A (en) * | 2016-11-03 | 2017-03-22 | 河海大学 | Three-dimensional random convex polygon block discrete element method based on distance potential function |
-
2018
- 2018-10-30 CN CN201811276841.XA patent/CN109408973B/en active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6099573A (en) * | 1998-04-17 | 2000-08-08 | Sandia Corporation | Method and apparatus for modeling interactions |
| CN104809301A (en) * | 2015-05-07 | 2015-07-29 | 北京理工大学 | Time domain substructure method for reflecting nonlinear rigid-soft hybrid junction characteristics |
| CN105912852A (en) * | 2016-04-08 | 2016-08-31 | 河海大学 | Arbitrary convex polygon block discrete unit method based on distance potential function |
| CN106529146A (en) * | 2016-11-03 | 2017-03-22 | 河海大学 | Three-dimensional random convex polygon block discrete element method based on distance potential function |
Non-Patent Citations (3)
| Title |
|---|
| LANHAO ZHAO ET AL: "A novel discrete element method based on the distance potential for arbitrary 2D convex elements", 《INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING》 * |
| 刘勋楠 等: "距离势离散单元法", 《岩石力学与工程学报》 * |
| 曹源: "三维离散元法软件的改进研究", 《中国优秀硕士学位论文全文数据库电子期刊 基础科学辑》 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109408973B (en) | 2021-06-08 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| JP5911077B2 (en) | Coupling calculation apparatus, coupling calculation method, and coupling calculation program for air, water and soil skeleton | |
| CN105676903B (en) | A kind of vibration optimal control system design method based on Multidisciplinary systems optimization | |
| CN104281730B (en) | A kind of finite element method of the plate and shell structure dynamic response of large rotational deformation | |
| CN106529146A (en) | Three-dimensional random convex polygon block discrete element method based on distance potential function | |
| Simeonov et al. | Nonlinear analysis of structural frame systems by the state‐space approach | |
| CN110717216A (en) | Method for forecasting rolling response of helicopter with flexible air bag under irregular wave | |
| CN105404758A (en) | Numerical simulation method of solid continuum deformation based on finite element method | |
| CN104298857A (en) | Mechanism reliability calculating method based on multi-factor coupling | |
| CN112949065A (en) | Double-scale method, device, storage medium and equipment for simulating mechanical behavior of layered rock mass | |
| CN113377014A (en) | Robust stabilization control method and system for mechanical arm system | |
| Barbosa et al. | Discrete finite element method | |
| Nguyen-Thoi et al. | Development of the cell-based smoothed discrete shear gap plate element (CS-FEM-DSG3) using three-node triangles | |
| CN109408973A (en) | One kind is based on apart from potential function two-dimensional variable shape convex polygon block discrete element method | |
| CN109408977A (en) | One kind is based on apart from the deformable three-dimensional convex polyhedron block discrete element method of potential function | |
| Riley et al. | Implementation issues and testing of a hybrid sliding isolation system | |
| Doolin | Unified displacement boundary constraint formulation for discontinuous deformation analysis (DDA) | |
| CN110826153B (en) | Water acting force simulation and realization method applied to helicopter water stability calculation | |
| Hung et al. | Predicting dynamic responses of frame structures subjected to stochastic wind loads using temporal surrogate model | |
| Betts et al. | Direct transcription solution of inequality constrained optimal control problems | |
| CN107092730A (en) | Suitable for the three-dimensional infinite element Artificial Boundaries method for building up of Explicit Analysis | |
| Qi et al. | Identifying critical loads of frame structures with equilibrium equations in rate form | |
| CN109977498A (en) | A method of FGM ladder beam dynamic response is calculated based on HOC | |
| Bibikov et al. | Development of mathematical models for calculating statics and dynamics of membrane sensitive elements of microelectromechanical control systems | |
| CN109284537B (en) | Deformable two-dimensional arbitrary rounding convex polygon discrete unit method | |
| Park et al. | Numerical Prediction of Soil Displacement Based on Imu Sensor Measurement |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |