其它数值方法(完整版)

第六章其它数值方法第一节离散单元法（详细）第二节边界单元法（自学）第三节数值方法的耦合（自学）主要内容第一节离散单元法一、概述二、基本原理四、UDEC（UniversalDiscreteElementCode）应用简介三、运动方程的解法推荐参考书：王泳嘉.离散单元法及其在岩土力学中的应用.东北工学院出版社，1991一、概述有限元和边界元法作为一种数值方法，在研究工程地质问题中是以连续介质力学为基础的。对于不连续面的模拟，有限元法一般采用节理单元，边界元法使用位移间断法，有时采用节理化岩体模型。但是这些方法在模拟岩质边坡的失稳运动、矿山开采的岩层移动等这类块体可能沿弱面作大位移或大角度旋转的运动则显得无能为力。离散元法是专门用来解决不连续介质问题的数值模拟方法。该方法把节理岩体视为由离散的岩块和岩块间的节理面所组成,允许岩块平移、转动和变形,而节理面可被压缩、分离或滑动。因此,岩体被看作一种不连续的离散介质。其内部可存在大位移、旋转和滑动乃至块体的分离,从而可以较真实地模拟节理岩体中的非线性大变形特征。离散元法的一般求解过程为:将求解空间离散为离散元单元阵,并根据实际问题用合理的连接元件将相邻两单元连接起来;单元间相对位移是基本变量,由力与相对位移的关系可得到两单元间法向和切向的作用力;对单元在各个方向上与其它单元间的作用力以及其它物理场对单元作用所引起的外力求合力和合力矩,根据牛顿运动第二定律可以求得单元的加速度;对其进行时间积分,进而得到单元的速度和位移。产生与发展：该方法是20世纪70年代发展起来的，目前已在数值模拟理论与工程应用方面取得了很大进展。该方法的基本特征在于允许各离散块体发生平动和转动，甚至发生分离，弥补了有限元法或边界元法的介质连续和小变形的限制，因而特别适合块状及裂隙介质的大变形及破坏问题的分析。该方法首先由P.A.Cundall首先提出，于20世纪80年代中期由王泳嘉等人引入到我国。目前，离散元法发展迅速，美国有二维和三维的离散元程序UDEC和3DEC。我国有2D-BLOCK和3D-BLOCK。应用领域：边坡、巷道与采场、地下开采、地震、爆炸、核废料储存、散体介质运动、断裂、地下水渗流、热传导等。离散单元法所用的求解方法有静态松弛法和动态松弛法。静态松弛法直接寻找块体失去平衡后达到再平衡时力和位移的关系，通过求解联立平衡方程组，并不断迭代使块体的残余余力和力矩趋于零，从而获得结果。这一方法不考虑粘性阻尼和计算时步，避免了动态松弛法确定这些参数时所遇到的困难。但该方法不足之处是在求解联立平衡方程组时，有时会碰到数值奇异或病态问题。因此，目前最常用的仍然是动态松弛法。动态松弛法是把非线性静力学问题化为动力学问题求解的一种数值方法。该方法的实质是对临界阻尼振动方程进行逐步积分。为了保证求得准确静态解，一般采用质量阻尼和刚度阻尼来吸收系统的功能。当阻尼系数取值小于某一临界值时，系统的振动将以尽可能快的速度消失，同时函数收敛于静态解。这种带有阻尼项的动态平衡方程，利用有限差分法按照时步在计算机上迭代求解就是所谓的动态松弛法。由于被求解方程是时间的线性函数，整个计算过程只需要直接代换，即利用前一次迭代的函数值计算新的函数值。因此，对于非线性问题也能加以考虑，这是动态松弛法的最大优点。二、基本原理（一）块体接触模型块体接触的力学模型示意图与刚度有关的粘性阻尼元件Cn、Ct，在接触点处吸收块体单元相对运动的能量。与质量和速度有关的粘性阻尼元件Cm吸收块体单元绝对运动的动能。当块体在接触点C发生切向滑移时，有库仑-莫尔元件U进行阻尼，并解除切向粘性阻尼元件Ct。当块体间有拉力时，则解除接触点切向力和法向力。（二）块体接触本构关系块体接触的本构关系是指块体间接触的力与位移关系。最简单的是边-角接触的库仑-莫尔关系，但最符合实际的是各种节理模型。这里仅讨论边-角接触的库仑-莫尔关系。块体之间的相互作用，可以假想为通过边-角间的“叠合”而发生，即一个块体的边与另一个块体的角之间的叠合。作用力的大小可以用“叠合量”来表示。设块体间的相互作用力与其相对位移成正比，则在块体接触处，由块体相对法向叠合量△un和剪切量△ut引起的力的法向增量△Fn和切向增量△Ft分别为：时刻增量理论：应变增量与应力增量一一对应关系对于块体间不允许出现拉力，故对于剪切力，其稳定状态有库仑-莫尔定量：式中、为接触处块体的内摩擦角和粘聚力。当Ft趋近于时，滑动即趋于发生。块体间的相互作用还包括阻尼力，在接触点其法向和切向分量Dn、Dt与接触点位移增量成正比，即：式中、为接触点法向和切向阻尼系数,与刚度系数有关。可由阻尼参数β与弹性刚度的乘积得到：本构模型与速度无关的接触模型与速度有关的接触模型角-边接触模型边-边接触模型单状态量摩擦模型双状态量摩擦模型法向力：假定块体之间的法向力Fn正比于它们之间沿法向“叠合”δn的大小，即：Fn=Kn

δn,Kn为法向刚度系数。“叠合”量δn是在数值计算时人为假定的一个量，其值的选取与计算精度要求等因素有关。剪切力：由于块体所受的剪切力与块体运动和加载的历史或路径有关，所以对于剪切力要用增量△Ft来表示。设两块体之间的相对位移为δt，则：Ft=Kt

δt,

Kt为切向刚度系数，δt为两块体之间的相对位移。破坏条件：法向力和切向力所表示的力与位移关系为弹性，但在某些情况下弹性关系是不成立的，需要考虑破坏条件。如当岩块受到张力分离时，作用在岩块表面上的法向力和剪切力随即消失。对于塑性剪切破坏的情况，需要在每次迭代时检查剪切力是否超过C+Fntanφ，其中，C为粘结力，φ为内摩擦角。如果剪切力超过该值，此时剪切力就取为该值。角-边接触模型：该模型认为块体之间不存在拉力，且当切向力Ft达到某一最大值时，就会发生塑性剪切滑移，并由下式确定：该模型的力与位移关系分别如下图所示：（a）法向力与法向位移（b）切向力与切向位移刚度系数的确定：对于如图所示的两个接触块体，其长度和宽度分别为a、b，弹性常数为E、μ。可得其法向刚度系数为：块体接触计算模型切向刚度和法向刚度之间的关系为：块体接触计算模型（三）运动方程-牛顿第二运动定律假设某一块体与周围n个块体接触，则其受到n个力作用，将各个力在x和y方向上分解，则在两个方向上的合力与合力矩分别为：块体集合及作用于个别块体上的力假设块体的质量为m，转动惯量为I，考虑重力，则块体质心的运动方程为：（5-72）三、运动方程的解法采用计算机实施数值计算的过程中通常会涉及到动态松弛法，所谓动态松弛法是把非线性静力学问题转化为动力学问题求解的一种数值方法。该方法的实质是对临界阻尼的振动方程进行逐步积分。对于这种带有阻尼项的动态平衡方程，利用有限差分法按时步在计算机上迭代求解，即所谓的动态松弛法。由于被求解的方程是时间的线性函数，整个计算过程只需要直接代换，因此，对于非线性问题也能加以考虑，这是动态松弛法的最大优点。用动态松弛法时，计算循环是以时步△t向前差分进行的。由于时步很小，每个单元在一个时步内只能以很小的位移与其相邻单元作用，而与其他单元无关，所以在一个时步内只能传递一个单元。计算循环示意图：计算循环示意图力-位移的关系力F运动定律F=ma位移u力的边界条件位移边界条件离散元法采用动态松弛法求解，其基本运动方程为：该方法是采用显式中心差分的动态松弛法进行求解，缺点是计算时步要很小，且需要合理确定阻尼系数等。该算法既可以由已知力求位移，也可以由已知位移求力。其总的算法为循环交错求解，首先根据边界条件已知位移u(t)，则可以求出F(t)，然后可确定,再加上前一轮的，对其积分可得u(t+△t/2),如此往复循环，直到位移和力收敛，达到平衡状态，或者达到某一循环数为止。计算循环示意图：离散元法交错循环求解特性图力位移加速度速度④③②①四、UDEC（UniversalDiscreteElementCode）应用简介Itasca公司的主要产品twodimensionalcontinuum,withjointsthreedimensionalcontinuum,withjointsthreedimensionalDEM*spheres&clumpstwodimensionalDEM*disks&clumpstwodimensionalDEM*polygonalbodiesthreedimensionalDEM*polyhedralbodiesUDEC非连续介质被作为是凸块或凹块的集合体，块体可以是刚体或变形体。不连续面被看作是块体之间的边界条件。沿着不连续面的运动是由切向和法向的线性及非线性力-位移关系控制的。地质体或类似材料可以由很多内置的块体和节理本构模型来描述；还有可供开发的自定义模型。是一种非常适合于用来模拟具有非线性力学行为的非连续材料（包含多个相互切割的非连续面）。尤其是，它具有如下特征：5.有平面-应变、平面-应力及轴对称三种问题模型。6.具有用来描述岩石-结构相互作用的结构单元模型，如锚杆、桩、梁、衬砌，喷射混凝土等。7.能够进行静力和动力问题分析8.沿着节理的暂态和稳态流.9.粘弹性和粘塑性（蠕变）模型.10.热分析能力，流固耦合分析能力OverviewofUDECoperation(1)-工程模拟通常由很长的一系列操作组成-UDEC数据文件可以很容易地利用文本编辑器进行修改，而且多个文件还可以连接在一起这种输入文件的方式提供了一个很好的可以用来保存分析稳定记录的方法.

命令流驱动结构允许通过操作程序的输入与输出来开发其前后处理程序UDEC

是命令流驱动程序例：new

(清除内存)

block0,0140,0

140,-99

0,-99

(创建块体)plotblock