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