1、3DEC数值模拟方法,2013 武汉依泰斯卡 (武汉）咨询有限公司,提 纲3DEC(5.0)软件特点3DEC(5.0)理论基础3DEC(5.0)新功能3DEC (5.0)算例分析,提 纲3DEC(5.0)软件特点3DEC(5.0)理论基础3DEC(5.0)新功能3DEC (5.0)算例分析,What is 3DEC?,A Distinct Element code that allows the simulation of the interaction of blockse.g. jointed rock mass, masonry wall,What is 3DEC?,Discontinu
2、ities are regarded as distinct boundary interactions between blocksContacts are automatically updated as calculation progresses large displacements are allowedNon-linear force displacement behavior may be assigned to contacts (joints),What is 3DEC?,3DEC blocks are rigid or deformable,Deformable bloc
3、ks are automatically meshed with tetrahedra to simulate a non-linear continuum,What is 3DEC?,Solves the full dynamic equations of motion even for quasi-static problemsAllows for large displacements, rotations and dynamic analyses,What is 3DEC?,3DEC is best suited to modeling discontinuous materials
4、(containing many intersecting discontinuities in 3D),that exhibit nonlinear behavior,矿层分布模型,地质构造分布,开挖后对地表建筑物的影响,屈服矿柱分布情况,矿山开采稳定性分析研究,当前开采条件下矿柱几乎没有屈服断层影响不大,矿藏分布三维模型,拉底步骤,矿岩三维节理模型,拉底后矿岩出露模型,矿岩崩落数值模拟,对崩落速度和规模进行预测,崩落过程模拟,Nuclear Waste Disposal,Shotcrete,Concrete,Canister cover,Localization in claystone
5、at 6 months,Claystone,Reproduced with permission of ANDRA,Hydraulic Fracture,Discrete Fracture Network,Fractured reservoir model,Hydraulic Fracture,Fluid pressure in joints during injectionTop portion of reservoir hidden for ease of viewing,Microseismicity,Fractured rock volume loaded axially with 1
6、0 MPa confining stress,100 m,Fundamental Rock Mechanics,Simulated unconfined compression test on crystalline rock sample,Displacements,Grain structure,Exaggerated deformation,Slope Stability,Automatic Factor of Safety calculation for 3D slopesUses Shear Strength Reduction method,Abutment stability a
7、ssessment at Hongrich arch dam1,1Koliji et al, 2001. Int. J. Hydropower & Dams, 18, 56-61,Toppling,Normal faulting in the upper part of the slope,S-shaped deformation of rock lamellae,Flexural toppling,Normal faulting,S-shaped deformed rock lamellae,提 纲3DEC(5.0)软件特点3DEC(5.0)理论基础3DEC(5.0)新功能3DEC (5.0
8、)算例分析,离散元方法,离散元的理论基础是牛顿第二定律。离散元法是将所研究的区域划分成一个个分立的多边形块体单元，单元与单元之间具有一定的初始接触状态，随着单元的平移和转动，调整各个单元之间的接触关系。最终，块体单元可能达到平衡状态,也可能一直运动下去。因此，离散元法适用于研究节理系统或块体集合在准静力或动力条件下的运动问题。离散元法以各种块体为研究对象，块体与块体之间没有变形协调的约束，只需要满足力的平衡条件。构成单元的节点可以分离。单元之间的相互作用力可以根据力和位移的关系求出，而单个单元的运动则完全根据该单元所受的不平衡力和不平衡力矩的大小按牛顿运动定律确定。,Types of Disc
9、rete Element Methods for Discontinuum Analysis,Distinct ElementUse explicit time-marching scheme to solve equations of motion directly. Bodies may be rigid or deformable, contacts are deformableModal MethodsSimilar to distinct element method in the case of rigid blocks. For deformable bodies, modal
10、superposition is used so non-linearity is difficult to implementDiscontinuous Deformation AnalysisAssumes contacts are rigid bodies and bodies may be rigid or deformable. No-penetration is achieved by iterationMomentum Exchange MethodsAssume both the contacts and bodies to be rigid. Momentum is exch
11、anged between two contacting bodies during collision. Can represent friction sliding.,DEM Definitions,The name “Discrete Element Method” (DEM) should be applied to a method only if it*:allows finite displacements and rotations of discrete bodies; including complete detachmentrecognizes new interacti
12、ons (contact) automatically as the calculation progressesA discrete element code will embody an efficient algorithm for detecting and classifying contacts. It will maintain a data structure and memory allocation scheme that can handle many hundreds or thousands of discontinuities or contacts.,The na
13、me “Distinct Element Method” is used for a DEM that uses an explicit dynamic solution to Newtons laws of motion.,*Cundall,P.A., and R.D. Hart.”Numerical Modeling of Discontinua,” Engineering Computations, 9(2), 101-113 (1992),Finite element codes for modeling “discontinua” are often modified continu
14、um programs, which cannot handle general interaction geometry (e.g. many intersecting joints). Their efficiency may degenerate drastically when connections are broken repeatedly.,Distinct Element Method,Three aspects .GeometryContact mechanicsSolid body mechanics,Distinct Element Method,GeometrySpec
15、ification of shapes in 2 and 3 dimensionsInteraction of pairs of contacting blocks or particlesIdentification of contact character between 2 blocks,Distinct Element Method,Specification of shapes in 2 and 3 dimensions2 & 3 dimensions PFC disks & spheres2 dimensions UDECarbitrary polygons convex and
16、concave with rounded corners3 dimensions 3DECarbitrary polyhedra concave bodies are constructed of several convex bodies attached together,Distinct Element Method,Interaction of pairs of contacting blocks or particlesIf we attempt to identify neighboring blocks by an exhaustive scan (i.e. each block
17、 tested against each other), then the search-time is proportional to N2 where N is the number of blocks.Two methods to reduce search time:Cell-mapping, used in UDEC & 3DECCirculating data structure, that mimics the topology of the system as used in UDEC,Cell-mapping,The solution-space is covered by
18、rectangular cells. Each block deposits a marker in all the cells that it overlaps this process takes N proportional time. Each block can find all of its neighbors by looking in just those cells that it overlaps - this process is also N proportional, if the cell size is of a similar order to the bloc
19、k size,Circulating data structure,Linked lists that describe block boundaries also trace out the void spaces automatically. A block needs only to search its local void spaces to find its neighbors. This scheme breaks down if many blocks become disconnected, but it is well suited to model fluid flow
20、in joints,closed domain,can restrict search to the closed domain any new contacts must be there,Distinct Element Method,Identification of contact character between 2 blocksWe need to know:type of contact (e.g. corner-to-corner, corner-to-edge, etc.)direction of normal to sliding directiongap between
21、 blocks, or contact overlap,Block 2,Joint,un2,us2,Block 1,us1,un1,Initial Positionof block 2,x,y,block centroid,Contact Between two Rigid Blocks,A contact is created at each corner interacting with a corner or edge of an opposing block.,Rounded Corners in UDEC,Corner rounding scheme with constant le
22、ngth d,r,d,d=r,d r,r,d,d=distance to the cornerr=radius of the rounded corner,r,d,d=r,d,r,d r,Corner rounding scheme with constant radius r, showingthat small angles in the corner leads to large distances d,Definition of contact normal,Rounded corner-to-edge contact,Rounded corner-to-corner contact,
23、L1,L2,L3,1,2,3,Element,Nodes,1,2,3,Corner-Edge contacts,L1, L2, L3,Lengths associatedto the contacts,Contacts and Domains between Two Deformable Blocks,D1,D2,D1,D2,Domains,BLOCK1,BLOCK2,“Common plane” logic,Algorithm executed in parallel with mechanical calculation: “maximize the gap between the com
24、mon plane (c-p) and the closest vertex”,c-p,1,2,“Common plane” logic,If we assume a c-p can be found, blocks are convexthen.test for contacts is easytest corners of both blocks for contact with c-p (16 dot products)if test for both blocks are positive then blocks are touchingdetermining gap is easys
25、um of min. distances from each blocks corner to c-pcontact normal is .the normal of the c-p is always defined,Common plane logic,What about assumptions?Use the algorithm “maximize the smallest gap between each block (of a contact-pair) and the plane”Construct concave blocks from several convex block
26、s (the “join” does not exist , as far as physics is concerned),shift, then.,rotation,Common plane logic,Modes of contact can be found easily with contact-plane logic“count the number of corners that touch the common plane”mode of contact is determined by counts: e.g. 1- 4 = “corner-face”note: ”touch
27、” implies gap 0,Block A1 corner2 edge 3 face,Block B1 corner2 edge 3 face,Distinct Element Method,Contact mechanicsAll contacts are assumed to be “soft”- i.e. contact forces are directly related to the deformations or “overlaps” at contacts. For point-contact, Hertz/Mindlin contact laws can be used.
28、 For edge-to-edge contact, such as rock joints, various constitutive laws can be used e.g., elastic/Coulomb slip. UDEC & 3DEC also have the ”continuously yielding” model, which employs continuous functions for the force displacement relations. There is also internal damage accumulation.,44/160,Disti
29、nct Element Method,Solid body mechanicsThere are two formulationsRigid body translation and rotationDeformable body mechanics,45/160,Rigid body translation and rotation,If most of the movement in a system takes place in a discontinuous way (e.g. sliding, opening, relative rotation, interlocking), th
30、en the bodies may be assumed to be rigid.,Distinct Element Method,46/160,Deformable body mechanics,If there is appreciable deformation of the intact material, compared to discontinuous motion, then the bodies must be taken as deformable.,Distinct Element Method,47/160,Deformable body mechanics,In or
31、der to model deformable blocks, there are automatic generators in UDEC to divide a body into triangles, and in 3DEC to divide bodies into tetrahedra. The finite-difference formulation for these internal elements is identical to that for FLAC. Several linear and non-linear constitutive models can be
32、used in the internal elements.,Distinct Element Method,48/160,The Explicit Dynamic Solution Schemeapplied with the Distinct Element Method,49/160,UDEC solves the full dynamic equations of motion even for quasi-static problems. This has advantages for problems that involve physical instability, such
33、as collapse. To model the “static” response of a system, a relaxation scheme is used in which damping absorbs kinetic energy. This approach can model collapse problems in a more realistic and efficient manner than other schemes, e.g., matrix-solution methods.00,Basis of the Solution Scheme,movie,50/
34、160,Overview of DEM & explicit, dynamic solution scheme,The formulation is very simple. For example, for a ball impacting a wall,mass,One time step,unknowns,knowns,(all contacts, in general),(all particles, in general),Full dynamic equations (integration of Newtons 2nd law),Explicit solution scheme,
35、Three consequences of this formulation are as follows ,(central difference 2nd order accurate),51/160,1. Treating each body as discrete (DEM) allows discontinuous material (such as a rock mass) to be modeled easily.,2. Full dynamic equations of motion allow the evolution of unstable systems to be si
36、mulated realistically.,3. Explicit solution scheme makes the task of handling nonlinearity trivial. Examples of nonlinearities are: (a) contact making e.g.,INPUT,OUTPUT,The explicit scheme uses a time step so small that information cannot propagate between neighbors in one step.,Thus, each element i
37、s isolated during one step, enabling,m,t,k,D ,S,52/160,Computation Cycle in UDEC,All the contacts,CONSTITUTIVE,At the centroid,ALL THE BLOCKS,MOVEMENT,zone,node,At the element,At the node,MOVEMENT,ALL THE BLOCKS,RIGIDBLOCKS,DEFORMABLEBLOCKS,Go to,53/160,TIMESTEP,The solution scheme used in DEM is co
38、nditionally stable. A limiting timestep must satisfy both the calculation for internal block deformation and inter-block relative displacement. For stability of the block deformation computation:where mi is the mass associated with block node i; and ki is the measure of stiffness of the elements sur
39、rounding the node. For calculations of the inter-block relative displacement, the limiting timestep iswhere Mmin is the mass of the smallest block in the system; and Kmin is the maximum contact stiffness. The controlling timestep is,54/160,MECHANICAL DAMPING,Mechanical damping is used in UDEC for st
40、atic (non-inertial) solutions and for dynamic solutions. For static analysis, the approach is similar to “dynamic relaxation.” The equations of motion are damped to reach a force equilibrium state as quickly as possible under the applied initial and boundary conditions.,55/160,DYNAMIC RELAXATION,In
41、dynamic relaxation blocks (and gridpoints) are moved according to Newtons law of motion. The equations of motion are damped to reach a force equilibrium state as quickly as possible Damping is velocity-proportional i.e., the magnitude of the damping force is proportional to the velocity of the block
42、s. Velocity-proportional damping introduces body forces that can affect the solution. UDEC provides two different types of damping to minimize this problem for static analysis:adaptive global damping (DAMP auto)local damping (DAMP local),56/160,ADAPTIVE GLOBAL DAMPING,Adaptive global damping is a se
43、rvo-mechanism used to adjust the damping constant automatically such that the power absorbed by damping is proportional to the rate of change of kinetic energy in the system.Adjustment of the viscosity constant is made by a numerical servo-mechanism that seeks to keep the ratio, R, equal to a given
44、ratio (e.g., 0.5).,Damping forces are introduced to the equations of motion:where Fi is the unbalanced force,57/160,LOCAL DAMPING,The damping force, Fd is:,In UDEC the unbalanced force ratio (ratio of unbalanced force, Fi , to the applied force magnitude, Fm) is monitored to determine the static sta
45、te. By default, when Fi / Fm 10-5 in UDEC then the model is considered to be in an equilibrium state.,The damping force at a gridpoint is proportional to the magnitude of the unbalanced force with the sign set to ensure that vibrational modes are damped.,58/160,BASIS OF STATIC SOLUTIONS IN UDEC,The
46、unbalanced force ratio (ratio of unbalanced force, Fi , to the applied force magnitude, Fm) is monitored to determine the static state.Both the adaptive global damping and local damping dynamic solution methods provide a static solution (with the effect of inertial forces minimized) provided the unb
47、alanced force ratio reaches a small value.This is comparable to the “level of residual error” or “convergence criterion” defined for matrix solution methods used in many finite element programs. In UDEC, the level of error is quantified by the unbalanced force ratio. In both UDEC and FE solutions, t
48、he static solution process terminates when the error is below a desired value.,59/160,Comparison of adaptive global damping and local damping,A simple example: ten blocks slide on a rough base*,Velocity proportional to out-of-balance force with proportionality constantDynamic relaxation with standard velocity-proportional dampingAdaptive global dampingLocal damping,*,60/160,
