1、离散单元法,曹 平,岩体是一种地质材料，它经受长期的地质构造作用，在一定的地质环境中形成一定的结构，显现出宽广和多变的材料响应范围。岩体与一殷工程材黔相比，其最大特点是一般具有结构上的不连续比多为层面和节理面等弱面所切割。实际k，岩体结构的连续性和不连续性那是相对的。如图l1所示，当取样范图较小时，可能碰到完整岩石，逐渐扩大取样范围，则是节理岩体。,一般从宏观意义上说，岩体可以视为连续介质，从而可以用弹性力学或塑性力学的方法来进行分析和计算。但在某些情况下，岩体却不能视做连续介质，如地下节理岩体中的巷道(见图12)，这时，就不宜用处理连续介质的力学方法来进行计算。于尼离散单元法作为一种处理节理
2、岩体的数值方法就应运而生。,离散元技术的威力可从下图模拟理想轴向压缩试验中看到。左边的承台推动试件撞击另一个固定的承台当试件中的应力超过强度极限时，试件出现裂纹。进一次压缩引起大变形，试件断裂，伴随着高速度出现一些飞离的碎片、应当注意断裂过程不用用户干预，断裂方向、新单元的产生及新的相互作用力都足由程序自动计算的。如用常规的有限元或有限差分法进行这种比较直观的离散元分析则是非常困难的。,工程中所见到的岩体其形态常呈非连续结构，所形成的岩石块体运动和受力情况多是几何或材料非线性问题，离散单元法充分考虑到岩体结构的不连续性，适用于解决节理岩石力学问题。离散单元法除了用于边坡、采场和巷道的稳定性研究
3、以及颗粒介质微观结构的分析外，已扩展到用于研究地震、爆炸等动力过程和地下水渗流、热传导等物理过程。,离散单元法的基本思想离散单元法也将区域划分成单元(见图13(a)。但是，单元因受节理等不连续面控制，在以后的运动过程中，单元结点可以分离，即一个单元与其邻近单元可以接触，也可以分开(见图13(b)。 单元之间相互作用的力可以根据力和位移的关系求出，而个别单元的运动则完全根据该单元所受的不平衡力和不平衡力矩的大小按牛顿运动定律确定。,离散单元法是一种显式求解的数值方法。该方法与在时域中进行的共他显式计算相似，例如与解抛物线型偏微分方程的显式差分格式相似。 “显式”是针对一个物理系统进行数值计算时所
4、用的代数方程式的性质而言。在用显式法计算时，所有方程式一侧的量都是已知的，而另一测的量只要用简单的代入法就可求根。这与隐式法不同，隐式法必须求解联立方程组。,2 刚性块体模型离散学元法是将所研究的区域划分成一个个分立的多边形块体单元，单元之间可以看成是角角接触、角边接触或边边接触，而且随着单元的平移和转孔允许调整各个单元之间的接触关系。最终，块体单元可能达到平衡状态，也可能一直运动下去。离题单元法的理论基础是结合不同本构关系的牛顿第二定律，因而可以采用动态松弛法或静态松弛法进行求解。离散单元法的原理比较简单，解决非连续介质大变形问题时非常实用。,离散单元法的单元，从性质上分可以是刚性的，也可以
5、是非刚性的；从几何形状上分，可以是任意多边形，也可以是圆形。下面将要介绍的是刚性块体模型，块体可以是任意多边形。刚性假设对于应力水平比较低的问题(如边坡稳定和放矿模拟等)是合理的。,22 离散单元法的基本方程在解决连续介质力学问题时，除了边界条件外，还有3个方程必须满足，即平衡方程、变形协调方程和本构方程。变形协调方程保证介质的变形连续；本构方程即物理方程，它表征介质应力和应变间的物理关系。对于离散单元法而言，由于介质一开始就假定为离散块体的集合(图21(a)，故块与块之间没有变形协调的约束，但平衡方程需要满足。,运动方程牛顿第二运动定律根据岩块的几何形状及其与邻近岩块的关系，可以利用上面讲述
6、的原理，计算出作用在其一特定岩块上的一组力，由这一组力不难计算出它们的合力和合力矩，并可以根据牛顿第二运动定律确定块体质心的加速度和角加速度，进而可以确定在时步dt内的速度和角速度以及位移和转动量。,对于块体沿y方向的运动及其转动有类似的算式。,23 离散单元法的计算机实施离散单元法的计算原理虽然很简单，但在计算机上实施起来却非常复杂，涉及到很多问题。本节将讨论离散单元法中的4个主要问题，即动态松弛法、力和位移的计算循环、分格检索以及数据结构，而将时步、阻尼以及前后处理等问题留在以后的章节中讨论。,231 动态松弛法离散单元法中所用的求解方法有“静态松弛法”和“动态松弛法”两种。松弛法作为解联
7、立方程组的一种方法在力学中有着重要的应用。其中，动态松弛法是把非线性静力学问题化为动力学问题求解的一种数值方法，该方法的实质是对临界阻尼振动方程进行逐步积分。,为了保证求得准静解，一般采用质量阻尼和刚度阻尼来吸收系统的动能，当阻尼系数取值稍小于某一临界值时，系统的振动将以尽可能快的速度消失，同时函数收敛于静态值。这种带有阻尼项的动态平衡方程，利用有限差分法按时步在计算机上迭代求解就是所谓的动态松弛法。,由于被求解方程是时间的线性函数，整个计算过程只需要直接代换，即利用前一迭代的函数值计算新的函数值，因此对于非线性问题也能加以考虑，这是动态松弛法的最大优点。其具体求解方法可以通过下面的简革例子来
8、说明。离散单元法的基本运动方程为：,从以上介绍不难看出，离散单元法利用中心差分法进行动态松弛求解，是一种显式解法。它不需要解大型矩阵，计算比较简单，也节省计算时间，并且允许单元发生很大的平移和转动，因此克服了以往有限单元法和边界单元法的小变形假设，可以用来求解一些非线性问题。,232 力和使移的计算循环在用动态松弛法时,计算循环(见图2-3)是以时步向前差分进行。由于时步选取得非常小，每个单元在一个时步内只能以很小的位移与其相邻接的单元作用，而与较远的单元无关系，所以力在一个时步内只能传递到一个单元(见图21)。用于求每个单元运动大小的基本方程则是牛顿第二运动定律。,由于本章中所讲的单元是刚性
9、的，因此单元之间的相对位移增量等完全由单元的几何尺寸、质心平移和单元绕其质心转角的大小来决定。如图25所示为在一个循环中所历经的计算过程。,离散元方法一变形体模型 41变形体动力学在对变形体单元进行数学推导之前，让我们先来熟悉要研究的问题。考察受有常力F作用的单个变形体单元(如图41),给定单元的几何、密度和弹性模量，需要确定它的运动特性。我们知道用牛顿定律可以确定质心的运动和绕质心的转动。但不能确定单元相对于质心的变形。,在变形体单元中，单元质心是在运动的，固定在单元质心处的一组坐标既有平动又有转动。为了得到相对于这样一组坐标系的运动变形方程有必要将通常惯性坐标系中的动力平衡方程转换成用单元
10、主轴表示的形式，单元主轴既作平动又作转动，如图42。,建立于这种坐标系中的单元具有零相对线动量和角动量，其变形可精确地用单元固有弹性模态的展开式来表示。模态法完全可用于大转动和大应变问题。对后一种情况,我们必须严格选择适当的应力应变度量方法，并注意变形前后的旋转路径。为避免复杂，我们在这里仅限于小应变问题。,动力平衡方程向非惯性坐标系的转换以下将讨论在惯性系中作任意大转动的弹性体所须遵循的方程。变形体单元的运动可分解为平均刚体运动和相对运动，刚体运动定义了一个非惯性动力参考系，相对运动指相对于该非惯性参考系的变形。,选择满足零相对线动量和角动员的Tisserand条件的动力参考系。对相对位移的
11、均方值而不是对相对动能取极小值，将得到线性形式的零角动量方程反之将得到非线性方程。这就使得相对位移能精确地用物体固有弹性模态的展开式来表示。运用哈米顿原理及广义坐标的概念，可以推导出所有的刚体运动方程和变形方程，在每一时步中采用小应变的假设，则对简化这些方程。这些方程不仅适用于小变形，同时也能正确地处理大变形的问题。,无论所选择的模态是特征模态还是非特征模态，应当注意的是在任何时候模态都必须以相对于绕单元转动的局部坐标系的形式给定,若模态方程写成用整体坐标表示，则当单元转动时。质量矩阵不再是对角矩阵。455 模态方程的旋转假定所得到的解藕后的单元模态方程是相对于固定的整体坐标系的。具体地说就是
12、让长方形单元的边与坐标铀对齐。,46 本构关系461 单元本构关系正如在以前几节中所看到的那样，离散元可完全用与有限元相同的方法构造。材料模型对两种方法也是相同的，因此任何适用于有限元的本构模型也适用于离散元。例如、离散元可描述具有任意形式的屈服面或塑性势的弹粘塑性性质。,除普通塑性破坏外,离散元还能反映脆性破坏,当单元中的应力状态超过了用户所定义的极限时，脆性断裂发生。单元分成两部分。新单元的生成是自动进行的，出于不需要形成总刚度矩阵，因此在重新组织数据方面不会发生问题(单元也可很容易地删除)；,采用不连续法模拟煤层的结果如图49和固410所示。图49表示煤层中初应力场的主应力分布，图410表示最终的府力分布，图4.ll表示假设煤层为连续介质时分析的结果。在不连续法分析中，由于滑动，与连续法的分析比较，通道顶部只能承受较小的水平应力。为了能足够地承受水平推力，随着深度的增加，应力也在增大。,
