建立适用于离散元方法的岩体构造面流变模型，在传统的离散元基本框架的基础上，引入积分型流变本构方程和粘性等效接触力增量并嵌入到离散元的计算流程中，成功模拟了离散介质结构面的流变性质。通过对一个简单算例计算值和理论解的对比，验证了本模型有良好的精度，具有应用于岩质高边坡和离散块体系统流变分析的重要前景。
 
　　关键词  离散单元法，流变，接触

　　离散单元法（DEM）是一种适用于模拟离散介质的数值方法。自Cundall[1]于70年代提出以来，离散元方法在岩石介质数值分析中得到了迅速的发展[2,3]，它与有限元、DDA等方法一起，成为岩体变形与稳定分析的有力工具。由于离散元方法既能模拟接触面的大变形，又能模拟块体内部的连续变形，它特别适合于模拟被节理、断层等构造面切割而成的不连续岩体介质。大量工程实践表明，岩体系统的变形和稳定问题与时间因素密切相关，在长期荷载持续作用下,其流变性态即使在不太高的应力水平下也十分明显，甚至成为最终变形破坏的主导因素。因此，在进行岩质块体结构的变形和稳定分析时，考虑介质的流变属性是十分必要的。当岩块本身坚硬完整时，其时效变形主要体现在节理、裂隙等构造面的流变行为上。

　　离散单元法由于计算时步很小，计算中所施加的阻尼主要是加快收敛速度，不能直接用于考虑较长时间效应的流变计算，作为采用离散元模型分析岩体流变效应的初步成果，本文建立了分析构造面流变的可变形体离散元数值模型。

1．离散单元法基本原理简介

　　可变形体离散单元法的突出特点是它不仅能考虑离散块体本身的变形，又能同时模拟接触面的张开和滑移等力学行为，由于采用显式时步迭代，对应于应力增量，本构方程采用增量形式，离散元计算方法对材料非线形本构关系和大变形问题具有内在的适应性。

　　如图1(a) 所示，岩体被裂隙、断层等构造面切割成块体系统。构造面采用法向和切向弹簧来模拟，构造面的接触力Fc可以根据构造面的法向和切向变形，通过接触弹簧的力-变形关系来确定。由于在不平衡力作用下各个块体不断运动，这些构造面上的接触力的大小、位置、方向均在不断变化，需要在计算过程中不断搜索、判断，这是离散元法的重要技术，正在不断地改进中[1]，具体内容可以参见文献[1,3]。

　　可变形体离散元法中，每个块体本身为可变形体，它的变形通过将块体内部划分成为三角形常应变差分单元模拟，如图1(b) 所示，网格节点N为差分单元B1，B2…Bn的公共节点。三角形差分单元把惯性质量平均分配给其三个节点，即把虚线 所围起的多边形的质量分配给N节点。因为每个三角形差分单元的应力、应变均为常量，可以由三角形各节点的位移和块体本身的本构关系来确定，再通过积分可以得到节点弹性力。

　　离散元计算采用显式时步迭代的方法，静力问题亦采用类似动力分析的松弛迭代求解，阻尼力Fd作为一种吸收能量加速系统收敛的手段被引入方程。离散元系统每个网格节点都应满足Newton第二定律，即 m =Fe + Fc + Fd + Fb+ Fg                                          (1)

　　其中m为网格节点所分配的质量，  是该节点的加速度。Fe，Fc，Fd，Fb，Fg分别为作用在该节点上的弹性力，接触力，阻尼力，给定的外力和重力。

　　离散元法用显式时步迭代进行平衡计算，对于每一计算时步，将对所有节点进行循环，计算步骤（如图2虚框所示）如下：

　　(a) 在迭代的初始时刻，设定边界条件和初始条件，由于边界条件或外加力等的变化，系统处于不平衡状态，开始进行平衡迭代；
　　(b)根据上一个时步的结果或初始条件、边界条件，可得到本时步初节点的位移、速度、坐标等状态量，如果是该节点的运动是给定的边界条件，则跳过(c)和(d)，直接得到新的状态量。否则，
　　(c) 根据节点的状态，利用各种力-位移关系分别确定作用在该节点上的各种作用力。
　　(d) 应用动力平衡方程式(1)，求出本时步该节点的加速度，再通过积分得到时步末的速度和位移，并更新坐标，得到新的状态量。 
　　(e)对每一时步内的每一节点，重复以上(b)~(d)，由于体系内部的阻尼力将不断耗散系统能量，迭代将收敛到静力问题的解。在静力分析中，为尽快耗散能量，可以采用自适应阻尼，对于离散单元法更详细的论述，可参见文献 [1,2,3]。 

　　2. 岩块构造面流变数值模型

　　作为研究离散元流变模型的第一步，针对许多情况下，由于岩块本身的流变变形相对于构造面流变变形小，我们首先建立构造面流变变形的分析模型。构造面变形可以分为瞬时变形和流变粘性时效变形（下称粘性变形）两部分。我们采用最简单的Kelvin模型（一对并联的弹簧和阻尼器）来模拟粘性变形部分，其他复杂的流变模型也可以用本文的方法类似地引入离散元方法中。
若从t＝0时刻开始，构造面上作用有接触力Fc 并保持常数，则Kelvin模型计算的构造面粘性变形δve随时间的变化过程为[4]
                                            
　　其中，E和h分别为Kelvin模型中弹簧的弹性模量和阻尼器的粘性系数。在实际的岩质边坡和大型地下洞室的开挖卸荷过程中，接触力的变化过程将十分复杂，(2)式不能适用。为了分析复杂条件下的流变行为，我们假定构造面上的接触弹簧的法/切向粘性变形δve可以采用积分型流变本构方程[4]：
                            
　　其中，Fc(t)分别代表接触力向量的法/切向分量；J(t)为该接触面的法/切向蠕变柔量，J(0)代表瞬时的法/切向柔量；
对应于Kelvin模型，令 ，有：
                                         (4)
并且有
  ；                              (5)
所以，Kelvin 模型的遗传型积分公式(3)可化为： 
                     (6)
将整个计算时段划分为若干等长的时步，并假设在一个流变计算时步Δt内，接触力和参数E， 保持不变，并注意到ti+1 = ti+Δt，可以得到递推公式
                       (7)
  为方便离散元计算，也可以写成增量形式的递推公式
                     (8)
递推初值
                   
其中， ，为ti 时刻接触面法/切向接触力相对于ti-1时刻的增量； ，为ti  时刻接触面法/切向粘性变形相对于ti-1时刻的增量。

　　3. 流变模型在离散元中的实现方法

　　离散元方法是通过显式时步步进迭代的方法，计算系统的静力响应，为保证计算精度和稳定性，在时步迭代中，采用的迭代时步步长Δt一般为10－4s~10-5s，而且在静力分析中，时步仅代表迭代的步数，并无实际的物理意义。另外，流变是在系统长期变形中才体现出来的属性，我们常常需要计算几天、几个月、甚至几年的粘弹性响应。因此，不可能通过简单的改变结构面的本构关系，直接通过时步步进方法进行流变计算。本文提出一种将流变本构关系在离散单元中实现的方法。首先将整个流变计算时段划分为若干等长的流变时步，假定在每个流变时步内，接触力和参数保持常数。定义接触面法/切向粘性等效接触力增量ΔFc (ti +) 为：
                                            
　　其中，Δδve(ti+1)由式(8)确定，K为接触面法/切向刚度。当ti-1～ti 流变时步的离散元迭代收敛以后，在ti～ti+1 流变时步开始时（ti+时刻），人为地把所有构造面的法/切向接触力由Fc(t)增加ΔFc (ti +)到Fc (ti+ )，即：     Fc (ti +)= Fc (ti)+ΔFc (ti +)                                      (10)

　　原来系统的平衡状态被人为打破，重新开始离散元迭代，直至达到迭代收敛，就可以得到该流变计算时步末（ti+1）的构造面的接触力和粘性变形。从而得到可变形体离散元系统中所有块体节点由构造面粘弹性引起的位移，整个系统的坐标得到更新，然后继续进行下一流变时步的计算。
用该数值模型计算岩体系统变形的步骤说明如下：

　　(a)在t=0时刻，系统由于开挖卸荷等原因，产生瞬时外荷载变化，系统失去平衡，通过离散元系统的平衡迭代(图2中虚框所示)，迭代收敛时得到系统的卸荷瞬时变形响应，节点坐标更新并达到受力平衡状态，接触弹簧产生瞬时变形δe并且卸荷引起了接触力变化。
　　(b)ti~ti+1 (i=0,1,2....) 流变时步应用公式(8) (9) (10)，计算出该时步初系统的接触面法/切向粘性等效接触力增量ΔFc ( ti+ )，人为改变法/切向接触力为Fc( ti+)，系统平衡状态因此被打破，节点受到不平衡力的作用，再次进行系统平衡迭代，直至达到迭代收敛，即得到该流变计算时步末的接触力和各节点和接触面的粘性变形δve，块体系统坐标等状态量得到更新。
　　(c) 重复过程(b)，依次进行下个流变时步计算，得到各个时刻系统的流变变形值。当岩体系统再次经历开挖卸荷时，重复(a)～(c)，计算出该次卸荷瞬时变形以及随后的粘性变形。由于计算模型采用的是遗传积分形式的流变本构方程，粘性位移采用递推形式表达，因此，每一次卸荷引起的粘性变形都将持续下去，在某一时步计算出的粘性变形为该时步以前的各次卸荷引起的粘性变形的叠加，反映了应力应变历史对粘性变形的影响。

　　需要注意的是，计算模型中出现了两个时间尺度的步长，一个是离散元系统进行平衡迭代所需要的步长，一般为10-4s~10-5s，在静力分析中没有实际的物理意义；另一个Δt为流变计算的时步长度，一般以若干天或月为一个步长，具有明确的物理意义。

　　4. 模型的验证

　　为了验证模型的正确性，设计了以下简单岩块结构分级卸荷的算例。岩柱被构造面切割成4个10m×10m的可变形岩块并安放在一个固定块体上，岩块内部进一步划分成三角形差分单元，如图3所示。设岩块密度r =2750kg/m3，杨氏模量E = 4×1010 Pa，泊松比n=0.22。构造面线法向弹性刚度kn=6.9×108Pa。Kelvin模型的流变力学参数（接触面法向）为：线刚度E=6.9×108Pa, 线粘性系数h=6.9×109Pa·d。分别在第0天，第10天，第20天依次卸去上面第1、2、3块，研究剩余块体的粘性变形。在计算中，认为卸荷瞬间完成。取流变计算时步Δt为1d，离散元平衡迭代时步Δt＝1.318× s，每次平衡迭代收敛后，即进行下一流变时步计算。岩块内部A，B，C点的竖向总位移随时间变化的计算值如图4所示，应用流变学理论，可以求出这个简单问题的理论解答，也显示在图4中供比较。

　　瞬时变形求得后，开始粘弹性变形（时效变形）的计算。首先，根据式(9)计算出粘性等效接触力增量，并据此人为改变系统接触面的接触力，此时，图6中I-I接触面的应力从b跃变到b'，由于接触力的改变，这时系统节点将受到不平衡力作用开始产生运动，图5中B点经过系统平衡迭代，位移从b点收敛到c点，b、c两点的位移之差即为该流变时步粘性变形的增加量，对应于图4中，B点位移为从第10d的b点增加到第11d的c点，得到了该流变时步末的位移解。与这一过程相对应，图6中接触面的应力由b'点通过平衡迭代收敛到c。值得指出的是，b和c两点的应力值是相同的。即，通过平衡迭代，结构面发生粘性变形，该变形引起的接触应力增量刚好抵消了接触面应力在计算时步初的人为改变量，接触面的应力得到了恢复，系统重新恢复到新的平衡状态。又重新进行下一流变时步c、d的计算。

　　本算例很好地模拟计算了由于瞬时卸荷引起的岩块和接触面的瞬时弹性变形以及由于接触面的粘性变形，块体位移随时间的变化过程。计算值和理论值符合的很好，验证了本文模型的正确性和精度。对模型进一步的完善和改进，并用于实际工程应用，将是今后需要继续深入研究的课题。

　　参考文献
[1] Cundall PA. UDEC-- A generalized distinct element program for modelling jointed rock[R], Final technical report to European Research Office, U.S. Army, 1980.
[2] Zhang Chuhan, Pekau OA, Jin Feng, Wang Guanglun. Application of distinct element method in dynamic analysis of high rock slopes and blocky structures[J], Soil Dyn. & Earthq. Eng. , 1997(16).
[3] 王刚，离散元流变模型及其在三峡工程中的应用[D]，清华大学，2000.
[4] 杨挺青. 粘弹性力学[M]， 武汉: 华中理工大学出版社，1990.

A Contact Rheological Model Embedded in Distinct Element Method

Wang Gang, Jin Feng & Xu Yanjie
(Department of Hydraulic Engineering, Tsinghua University, Beijing 100084,China)

Abstract  A numerical model embedded in Distinct Element Method(DEM) is presented for simulation of contact rheological behavior in rocks. Based on the basic frame of DEM, the rheological constitutive equation and equivalent incremental contact force are introduced and embedded in the current calculation procedure of DEM. Then the responses of discrete blocks due to rheological deformation of contact interfaces can be simulated. Comparison between the model calculation results and theoretical solutions of a simple benchmark problem shows that the presented model has excellent accuracy, and possesses significant potential for analyzing rheological responses of high rock slopes and blocky structures.

Keywords:  Distinct Element Method (DEM), rheology, contact