udec数值方法学习课程

1、 离散元法是专门用来解决不连续介质问题的数值模拟方法。离散元法是专门用来解决不连续介质问题的数值模拟方法。该方法把节理岩体视为由离散的岩块和岩块间的节理面所组成该方法把节理岩体视为由离散的岩块和岩块间的节理面所组成,允许岩块平移、转动和变形允许岩块平移、转动和变形,而节理面可被压缩、分离或滑动。而节理面可被压缩、分离或滑动。因此因此,岩体被看作一种不连续的离散介质。其内部可存在大位岩体被看作一种不连续的离散介质。其内部可存在大位移、旋转和滑动乃至块体的分离移、旋转和滑动乃至块体的分离,从而可以较真实地模拟节理从而可以较真实地模拟节理岩体中的非线性大变形特征。岩体中的非线性大变形特征。一、概述第

2、1页/共42页第一页，编辑于星期日：一点 五十九分。 离散元法的一般求解过程为离散元法的一般求解过程为:将求解空间离散为离散元单元将求解空间离散为离散元单元阵阵,并根据实际问题用合理的连接元件将相邻两单元连接起来并根据实际问题用合理的连接元件将相邻两单元连接起来;单元间相对位移是基本变量单元间相对位移是基本变量,由力与相对位移的关系可得到两单元由力与相对位移的关系可得到两单元间法向和切向的作用力间法向和切向的作用力;对单元在各个方向上与其它单元间的作对单元在各个方向上与其它单元间的作用力以及其它物理场对单元作用所引起的外力求合力和合力矩用力以及其它物理场对单元作用所引起的外力求合力和合力矩,根

3、据牛顿运动第二定律可以求得单元的加速度根据牛顿运动第二定律可以求得单元的加速度;对其进行时间积分对其进行时间积分,进而得到单元的速度和位移。进而得到单元的速度和位移。 离散单元法所用的求解方法有静态松弛法和动态松弛法。离散单元法所用的求解方法有静态松弛法和动态松弛法。第2页/共42页第二页，编辑于星期日：一点 五十九分。第3页/共42页第三页，编辑于星期日：一点 五十九分。产生与发展：产生与发展：该方法是该方法是20世纪世纪70年代发展起来的，目前已在数值模拟年代发展起来的，目前已在数值模拟理论与工程应用方面取得了很大进展。该方法的基本特征在于允许各离理论与工程应用方面取得了很大进展。该方法的

4、基本特征在于允许各离散块体发生平动和转动，甚至发生分离，散块体发生平动和转动，甚至发生分离，弥补了有限元法或边界元法弥补了有限元法或边界元法的介质连续和小变形的限制，的介质连续和小变形的限制，因而特别适合块状及裂隙介质的大因而特别适合块状及裂隙介质的大变形及破坏问题的分析。变形及破坏问题的分析。 该方法首先由首先提出，于该方法首先由首先提出，于20世纪世纪80年代中期由王泳嘉等人引入年代中期由王泳嘉等人引入到我国。目前，离散元法发展迅速，美国有二维和三维的离散元程序到我国。目前，离散元法发展迅速，美国有二维和三维的离散元程序UDEC和和3DEC。我国有。我国有2D-BLOCK和和3D-BLOC

5、K。应用领域应用领域：边坡、巷道与采场、地下开采、地震、爆炸、核废：边坡、巷道与采场、地下开采、地震、爆炸、核废料储存、散体介质运动、断裂、地下水渗流、热传导等。料储存、散体介质运动、断裂、地下水渗流、热传导等。第4页/共42页第四页，编辑于星期日：一点 五十九分。二、基本原理（一）块体接触模型第5页/共42页第五页，编辑于星期日：一点 五十九分。（二）块体接触本构关系 块体接触的本构关系是指块体间接触的力与位移关系。最简单块体接触的本构关系是指块体间接触的力与位移关系。最简单的是边的是边-角接触的库仑角接触的库仑-莫尔关系，但最符合实际的是各种节理模型。这莫尔关系，但最符合实际的是各种节理模

6、型。这里仅讨论边里仅讨论边-角接触的库仑角接触的库仑-莫尔关系。莫尔关系。FnFnun(a)ut(b)FtFt时刻t的位置时刻t+t的位置第6页/共42页第六页，编辑于星期日：一点 五十九分。 块体之间的相互作用，可以假想为通过边块体之间的相互作用，可以假想为通过边-角间的角间的“叠合叠合”而而发生，即一个块体的边与另一个块体的角之间的叠合。作用力的大小发生，即一个块体的边与另一个块体的角之间的叠合。作用力的大小可以用可以用“叠合量叠合量”来表示。来表示。FnFnun(a)ut(b)FtFt时刻t的位置时刻t+t的位置第7页/共42页第七页，编辑于星期日：一点 五十九分。 设块体间的相互作用力

7、与其相对位移成正比，则在块体接触设块体间的相互作用力与其相对位移成正比，则在块体接触处，由块体相对法向叠合量处，由块体相对法向叠合量un和剪切量和剪切量ut引起的力的法向增引起的力的法向增量量Fn和切向增量和切向增量Ft分别为：分别为：tttnnnKFKF)()()()()()(tFtFttFtFtFttFtttnnntt时刻时刻增量理论：应变增量理论：应变增量与应力增量增量与应力增量一一对应关系一一对应关系第8页/共42页第八页，编辑于星期日：一点 五十九分。0nFjjnttcFFFtanmaxjjc对于块体间不允许出现拉力，故对于块体间不允许出现拉力，故 对于剪切力，其稳定状态有库仑对于剪

8、切力，其稳定状态有库仑-莫尔定量：莫尔定量： 式中式中 、为接触处块体的内摩擦角和粘聚力。为接触处块体的内摩擦角和粘聚力。 当当Ft趋近于趋近于 时，滑动即趋于发生。时，滑动即趋于发生。 maxtF第9页/共42页第九页，编辑于星期日：一点 五十九分。 块体间的相互作用还包括阻尼力，在接触点其法向和切向分量块体间的相互作用还包括阻尼力，在接触点其法向和切向分量Dn、Dt与接触点位移增量成正比，即：与接触点位移增量成正比，即： tttnnnCDCDnCtC 式中式中 、为接触点法向和切向阻尼系数，与刚度系为接触点法向和切向阻尼系数，与刚度系数数 有关。可由阻尼参数有关。可由阻尼参数与弹性刚度的乘

9、与弹性刚度的乘积得到：积得到：ttnnKCKC第10页/共42页第十页，编辑于星期日：一点 五十九分。本构模型本构模型与速度无关的与速度无关的接触模型接触模型与速度有关的接与速度有关的接触模型触模型角角- -边接触模型边接触模型边边- -边接触模型边接触模型单状态量摩擦模型单状态量摩擦模型双状态量摩擦模型双状态量摩擦模型第11页/共42页第十一页，编辑于星期日：一点 五十九分。法向力：法向力：假定块体之间的法向力假定块体之间的法向力Fn正比于它们之间沿法向正比于它们之间沿法向“叠叠合合”n的大小，即：的大小，即：Fn=Kn n, Kn为法向刚度系数。为法向刚度系数。 “叠合叠合”量量n是在数值

10、计算时人为假定的一个量，其值的选取是在数值计算时人为假定的一个量，其值的选取与计算精度要求等因素有关。与计算精度要求等因素有关。第12页/共42页第十二页，编辑于星期日：一点 五十九分。剪切力：剪切力：由于块体所受的剪切力与块体运动和加载的历史或路由于块体所受的剪切力与块体运动和加载的历史或路径有关，所以对于剪切力要用增量径有关，所以对于剪切力要用增量F Ft t来表示。设两块体之间来表示。设两块体之间的相对位移为的相对位移为t t，则：，则：Ft=Ktt, , Kt为切向刚度系数，为切向刚度系数， t为两块为两块体之间的相对位移。体之间的相对位移。第13页/共42页第十三页，编辑于星期日：一

11、点 五十九分。破坏条件：破坏条件：法向力和切向力所表示的力与位移关系为弹性，法向力和切向力所表示的力与位移关系为弹性，但在某些情况下弹性关系是不成立的，需要考虑破坏条件。但在某些情况下弹性关系是不成立的，需要考虑破坏条件。如当岩块受到张力分离时，作用在岩块表面上的法向力和如当岩块受到张力分离时，作用在岩块表面上的法向力和剪切力随即消失。对于塑性剪切破坏的情况，需要在每次剪切力随即消失。对于塑性剪切破坏的情况，需要在每次迭代时检查剪切力是否超过迭代时检查剪切力是否超过c+Fntan，其中，其中，c为粘结力，为粘结力，为内摩擦角。如果剪切力超过该值，此时剪切力就取为内摩擦角。如果剪切力超过该值，此

12、时剪切力就取为该值。为该值。第14页/共42页第十四页，编辑于星期日：一点 五十九分。角角- -边接触模型：边接触模型：该模型认为块体之间不存在拉力，且当切向力该模型认为块体之间不存在拉力，且当切向力Ft达到某一最大值时，就会发生塑性剪切滑移，并由下式确定：达到某一最大值时，就会发生塑性剪切滑移，并由下式确定： maxtantntFFF 该模型的力与位移关系分别如下图所示：该模型的力与位移关系分别如下图所示：（a a）法向力与法向位移）法向力与法向位移（b b）切向力与切向位移）切向力与切向位移第15页/共42页第十五页，编辑于星期日：一点 五十九分。刚度系数的确定：刚度系数的确定：对于如图所

13、示的两个接触块体，其对于如图所示的两个接触块体，其长度和宽度分别为长度和宽度分别为a、b，弹性常数为，弹性常数为E、。可得其法。可得其法向刚度系数为：向刚度系数为： bEakabESKEnnnn2S12个单位，则设块体厚度为块体接触计算模型块体接触计算模型第16页/共42页第十六页，编辑于星期日：一点 五十九分。切向刚度和法向刚度之间的关系为：切向刚度和法向刚度之间的关系为：为剪切模量所以，又Gbakkktnt2G)2(1EG12块体接触计算模型块体接触计算模型第17页/共42页第十七页，编辑于星期日：一点 五十九分。（三）运动方程- -牛顿第二运动定律 假设某一块体与周围假设某一块体与周围n

14、 n个块体接触，则其受到个块体接触，则其受到n n个力作用，将各个力在个力作用，将各个力在x x和和y y方向上分解，则在两个方向上分解，则在两个方向上的合力与合力矩分别为：方向上的合力与合力矩分别为：nixixff1niyiyff1niixiiyiyyfxxfM100块体集合及作用于个别块体上的力块体集合及作用于个别块体上的力第18页/共42页第十八页，编辑于星期日：一点 五十九分。 假设块体的质量为假设块体的质量为m，转动惯量为，转动惯量为I，考虑重力，则，考虑重力，则块体质心的运动方程为：块体质心的运动方程为：mfuxx mfuyy IM niiniiyiyyyniixixxxMcImg

15、DFucumDFucum1.1.1.)()(（5-72）第19页/共42页第十九页，编辑于星期日：一点 五十九分。三、运动方程的解法采用计算机实施数值计算的过程中通常会涉及到采用计算机实施数值计算的过程中通常会涉及到动态松弛法，所谓动态松弛法是把非线性静力学问动态松弛法，所谓动态松弛法是把非线性静力学问题转化为动力学问题求解的一种数值方法。该方法题转化为动力学问题求解的一种数值方法。该方法的实质是对临界阻尼的振动方程进行逐步积分。对的实质是对临界阻尼的振动方程进行逐步积分。对于这种带有阻尼项的动态平衡方程，利用有限差分于这种带有阻尼项的动态平衡方程，利用有限差分法按时步在计算机上迭代求解，法按

16、时步在计算机上迭代求解，即所谓的动态松弛即所谓的动态松弛法法。由于被求解的方程是时间的线性函数，整个计由于被求解的方程是时间的线性函数，整个计算过程只需要直接代换，因此，对于非线性问题也算过程只需要直接代换，因此，对于非线性问题也能加以考虑，能加以考虑，这是动态松弛法的最大优点这是动态松弛法的最大优点。第20页/共42页第二十页，编辑于星期日：一点 五十九分。 用动态松弛法时，计算循环是以时步用动态松弛法时，计算循环是以时步t t向前差分进行向前差分进行的。由于时步很小，每个单元在一个时步内只能以很小的。由于时步很小，每个单元在一个时步内只能以很小的位移与其相邻单元作用，而与其他单元无关，所以

17、在的位移与其相邻单元作用，而与其他单元无关，所以在一个时步内只能传递一个单元。一个时步内只能传递一个单元。第21页/共42页第二十一页，编辑于星期日：一点 五十九分。计算循环示意图：计算循环示意图：计算循环示意图计算循环示意图力力- -位移的关系：位移的关系：K力力F F运动定律运动定律F=maF=ma位移位移u u力的边界条件力的边界条件位移边界条件位移边界条件第22页/共42页第二十二页，编辑于星期日：一点 五十九分。离散元法采用动态松弛法求解，其基本运动方程为：离散元法采用动态松弛法求解，其基本运动方程为：该方法是采用显式中心差分的动态松弛法进行求解，缺点该方法是采用显式中心差分的动态松

18、弛法进行求解，缺点是计算时步要很小，且需要合理确定阻尼系数等。是计算时步要很小，且需要合理确定阻尼系数等。该算法既可以由已知力求位移，也可以由已知位移求力。该算法既可以由已知力求位移，也可以由已知位移求力。其总的算法为循环交错求解，首先根据边界条件已知位移其总的算法为循环交错求解，首先根据边界条件已知位移u(t)u(t)，则可以求出则可以求出F(t)F(t)，然后可确定，然后可确定 , ,再加上前一轮再加上前一轮的的 ，对其积分可得，对其积分可得u(t+u(t+t/2),t/2),如此往复循环，如此往复循环，直到位移和力收敛，达到平衡状态，或者达到某一循环数为止。直到位移和力收敛，达到平衡状态

19、，或者达到某一循环数为止。 tftkutuctum tu 2ttu第23页/共42页第二十三页，编辑于星期日：一点 五十九分。计算循环示意图：计算循环示意图：离散元法交错循环求解特性图离散元法交错循环求解特性图力力位移位移加速度加速度速度速度第24页/共42页第二十四页，编辑于星期日：一点 五十九分。UDEC1.1. 非连续介质被作为是凸块或凹块的集合体，块体可以是刚体非连续介质被作为是凸块或凹块的集合体，块体可以是刚体或变形体。或变形体。2.2. 不连续面被看作是块体之间的边界条件。不连续面被看作是块体之间的边界条件。3.3. 沿着不连续面的运动是由切向和法向的线性及非线性力沿着不连续面的运

20、动是由切向和法向的线性及非线性力- -位移关系位移关系控制的。控制的。4.4. 地质体或类似材料可以由很多内置的块体和节理本构模型来描述地质体或类似材料可以由很多内置的块体和节理本构模型来描述；还有可供开发的自定义模型。；还有可供开发的自定义模型。是一种非常适合于用来模拟具有非线性力学行为的是一种非常适合于用来模拟具有非线性力学行为的非连续材料（包含多个相互切割的非连续面）。尤非连续材料（包含多个相互切割的非连续面）。尤其是，它具有如下特征：其是，它具有如下特征：四、UDEC（Universal Discrete Element Code）应用简介第25页/共42页第二十五页，编辑于星期日：一

21、点 五十九分。5. 5. 有平面有平面- -应变、平面应变、平面- -应力及轴对称三种问题模型。应力及轴对称三种问题模型。6. 6. 具有用来描述岩石具有用来描述岩石- -结构相互作用的结构单元模型，如锚杆、结构相互作用的结构单元模型，如锚杆、桩、梁、衬砌，喷射混凝土等。桩、梁、衬砌，喷射混凝土等。7. 7. 能够进行静力和动力问题分析能够进行静力和动力问题分析8. 8. 沿着节理的暂态和稳态流沿着节理的暂态和稳态流. .9. 9. 粘弹性和粘塑性（蠕变）模型粘弹性和粘塑性（蠕变）模型. .10. 10. 热分析能力，流固耦合分析能力热分析能力，流固耦合分析能力第26页/共42页第二十六页，编

22、辑于星期日：一点 五十九分。Overview of UDEC operation (1)- 工程模拟通常由很长的一系列操作组成工程模拟通常由很长的一系列操作组成- - UDECUDEC数据文件可以很容易地利用文本编辑器进数据文件可以很容易地利用文本编辑器进 行修改，而且多个文件还可以连接在一起行修改，而且多个文件还可以连接在一起- 这种输入文件的方式提供了一个很好的可以用来这种输入文件的方式提供了一个很好的可以用来 保存分析稳定记录的方法保存分析稳定记录的方法. .- 命令流驱动结构允许通过操作程序的输入与输出命令流驱动结构允许通过操作程序的输入与输出 来开发其前后处理程序来开发其前后处理程序UDEC 是命令流驱动程序第27页/共42页第二十七页，编辑于星