摘要
     本文提出了一种二维变形体离散元与时域边界元的耦合模型，这一模型可以将非连续体的模拟与无限域的模拟统一在一个模型中，可用于在地震波动输入条件下，考虑辐射阻尼的岩体边坡或地下结构等的动力稳定和变形分析，拓宽了离散元动力分析的领域。算例分析表明本耦合分析模型具有较高的精度。
    关键词：离散元、边界元、耦合

    离散元法是一种模拟离散介质的计算方法，自Cundall在70年代提出以来，在岩石力学、土力学、结构分析等领域的数值模拟中得到广泛应用，是一种新兴的非连续体分析方法[1,2,3]。在动力分析中，许多分析均表明辐射阻尼对分析结果影响很大，应该充分重视 [4]，动力边界元方法由于在其基本解中包含了无限远处的辐射条件，在处理辐射阻尼的影响时十分方便，是分析地下结构动力响应的一种有力的工具[5]。动力边界元法可以分为频域边界元法[5]与时域边界元法[6,7]，前者主要适用于线性问题，是发展较早、相对成熟的动力边界元方法。时域边界元方法直接在时域内求解，能够解决非线性问题，根据求解的问题具有非线性的特点，本文选择动力时域边界元方法与可变形体离散元耦合，提出了一个二维时域动力边界元-可变形体离散元耦合模型，充分发挥离散元与边界元的优点，将非连续体的模拟与无限介质辐射阻尼的模拟统一到一个模型中，为地下结构和岩质边坡的抗震稳定分析提供了全新的手段。

1. 离散元-边界元耦合模型
1.1 二维可变形体离散元原理
    可变形体离散元法将模拟的计算区域看作是若干可变形块体的组合，见图1。这些块体可以任意平移、旋转，块体之间的相互作用力，用法向和切向弹簧表示，称为接触，接触力的大小由块体的相对位置决定。每个块体又划分为三角形差分网格以模拟变形，每个差分三角形的顶点称为节点，其应变假定为常数，可以由节点的位置确定。因为一旦节点位置确定，便可以得到差分三角形乃至整个块体的变形和应力，从而得到所有的响应历程，因此，所有的计算将围绕节点进行，边界条件也可以通过给定边界节点位移或节点力来实现。

    离散元计算采用显式步进的方法，首先将计算的过程分为若干等长的时步，静力分析作为动力分析的特例，可以采用临界阻尼以增加收敛速度，此时静力分析的时步不具有真实的物理意义，可看作迭代步。在每一时步内，对所有块体的所有节点分别进行循环计算，每一节点循环计算的主要步骤为：

 (1) 根据上一循环的结果或边界条件，确定本时步初节点的位置和速度等运动量。
 (2) 根据本节点与相邻节点的位置可确定差分网格的应变，根据本构关系求出应力，再积分得到本节点所受弹性力。若为块体边界上的节点，可根据相邻块体的位置确定接触力，从而得到本节点在本时步所受合力F'。 以上具体公式和算法可参见文献[1,2,3]。
 (3) 由牛顿第二定律或动力平衡方程
    (1) 式中， 是接触力Fc 、弹性力Fe与刚度阻尼力Fd=b(Fc + Fe)之和； 代表质量阻尼力；a,b 是Rayleigh阻尼系数；g代表重力加速度，最后一项表示重力的影响。由式(1)通过中心差分，可得到本时步(第i时步)的速度
    (2) 进而可以得到本时步末本节点的位置。返回第(1)步，继续下一轮循环计算，最终求得所有时刻所有节点的解。

    由于任意块体间在运动过程中都有发生各种接触的可能，在计算中需要对所有的块体之间进行接触判断并计算接触力，这是一项十分耗时的工作，也是离散元法的关键技术之一，为减少计算工作量，针对岩体中构造面的特点，已有一些较好的算法，如Cundall提出的域算法，充分利用生成离散块体时的信息，将接触检索局限于初始顶点构成的域中，特别适合于由构造面切割形成的离散块体系统，由于篇幅所限，在此不能详细介绍，请参阅有关文献[1,2]。

1.2 二维时域边界元基本原理
    二维全平面时域动力边界元方程[6,7,8 ]  (3)式中， 为二维全平面时域动力基本解，其具体的表达式可参见文献[6]和[8]，u，p分别表示位移与面力，uI 表示从无限远处入射的位移场，S，Q分别表示源点及场点。同时对时间、空间进行离散，可以得到时域边界元方程可写作：[6,7 ,8]  (4)其中，{u}L, {p}L分别是L时刻的位移与面力分量，  (5)其中，[H]l,[G]l,(l=1,L)共2L个矩阵均为系数矩阵，他们可由各边界单元的子矩阵集成而得，具体表达式可参见文献[6]和[8]。

    由于{u}l, {p}l (l=1, 2, …, L-1)是1，2，…，L -1时刻的量，在求解L 时刻时 为已知量，而在{u}L, {p}L中有N(总的自由度数)个值已知，另N个值未知，因此可解此方程组，可求得L时刻的N个未知量，依此类推，即可求得全部时程的历程反应。

1.3耦合模型的实现
    为了实现时域边界元和离散元的耦合，必须解决以下问题：(1)离散元模拟的离散体与边界元分析的连续体实现耦合；(2)求解方式的统一是实现耦合的关键；(3)离散元计算时步与边界元计算时步的协调及离散元节点与边界元单元节点的对应可以进一步提高耦合计算的效率。以下，分别介绍解决的方法。

    (1) 首先，我们将计算区域划分为两个区域，见图2，即主要考虑无限域辐射阻尼和地震输入的连续体区域，称为边界元域，用时域边界元(BEM)模拟，其模拟的重点是无限域，地表水平自由边界在离散一段距离后截断，经过收敛计算证明可以很好地模拟远域的影响[8]。另一区域，称为离散元域，用离散元(DEM)模拟，模拟的重点是离散体。在远域与近域之间设立一个界面块(Interface Block)，它是一个完整的可变形离散块体，可进一步划分为三角差分网格，本身是连续的，差分节点同时也可以是边界元节点，边界元与离散元的耦合完全在过渡块上进行，这样，就避免了一个边界元节点可能与两个块体相连而造成的非连续问题。

    (2) 在求解方式上，将整个边界元区域看成是一个离散元法的块体，离散元法计算的核心是根据位置、位移、速度等运动量求作用力，为此，我们改写边界元方程(4)，可以得到  (6)式中，[R]是将面力{p}L转换成作用力的转换矩阵。这样，只要知道过渡块中与边界元相接的节点的位移，就可以根据(6)式求出边界元区域作用在这些节点的作用力，进而用(2)式计算节点速度、位移，从而完成一次循环，实现了两种模型的耦合。

    (3) 为了保证计算效率和计算稳定性，时域边界元方法的计算时步不能太小，通常为离散元方法计算时步的几十到上百倍。为解决这一难题，我们采用了异步计算的办法。设离散元法的时步为Dt，则取时域边界元法的时步Dt=K1*K2*Dt，即将每个边界元时步划分为K1个小时步，每个小时步由K2个离散元时步组成，并认为一个小时步内边界元域的作用力不变，既每K2个离散元时步才调用一次(6)式更新作用力。在一个边界元时步内，式(6) 的各系数[G] L，[H] L，{B} L皆保持不变，{u}L将不断按下式更新：
    k =1,2,…,K1 (7)
式中， 是第k个小时步中第一个离散元时步初的节点位移。
当一个边界元时步结束后，(6)式中的面力向量也相应按下式更新（6）进而更新{B}，再进行下一个边界元时步的计算。

    为保证计算效率，边界元的节点通常要大大少于耦合边界上离散元的节点，我们采用一个边界单元对应若干离散元节点的办法，一个边界元内相互作用面力按线性分布，再根据力和力矩的平衡，分别计算分配给每个离散元节点的相互作用力 。这样，通过空间和时间的异步，大大提高了计算效率和改善计算稳定性，从而实现了时域边界元方法和离散元方法的耦合。

2.模型验证
为验证耦合模型的精度，计算了多个算例，因篇幅所限，仅列出以下两个算例。其他算例可见文献[9]。
2.1 岩柱算例
     如图3所示的一个岩柱，分为离散元块体模拟的4个岩块和一个边界元区域模拟的顶部岩块，最下一个离散岩块固定，在边界元岩块最上的边界作用单位阶跃荷载1.0MN/m2，岩块的弹模E = 30 GPa，密度r = 2000kg/m3，阻尼比x = 0.1，Kn= Ks = 1011 N/m， f =1.18， c = 0。由于在计算条件下，所有构造面保持完整接触，能保证岩柱的连续性，所以可同时采用时域边界元法对整个岩柱按连续体进行分析。

    图4为岩柱上各点的位移响应及其与时域边界元法的计算结果对比，在整个计算时间内，两者吻合很好。这一算例显示了当离散块体紧密结合成一个连续体时，耦合模型能够与连续体模型分析结果吻合，并可清楚地看到波动在岩柱中的传播过程，波动能够准确地在边界元区域与离散元区域之间和离散块体之间传播，从一个侧面证实了耦合模型及软件的正确性。

2.2 半平面的离散块受入射波作用
    为进一步验证耦合模型在处理波入射动力问题的计算精度,首先计算了如图5所示的半平面上的离散块体在竖直向上的SV波入射情况下的响应，离散元域除界面块外仅有一个块体，并且假定离散块体之间的摩擦系数与抗拉强度足够大以保证他们之间不会滑动，边界元域共采用了26个单元，单元长度均为10米，其中6个单元与界面块体接触，在SV波入射的位移历程见图6。

    图7示出了A点的水平位移响应与完全采用时域边界元方法计算结果进行的比较。两种方法计算结果的对比再次说明耦合模型具有较高的精度，同时说明了耦合模型能够模拟波动从无限远入射并正确输入到离散元域。其他证明耦合模型能够模拟离散块体的开合等大变形行为的算例，因为篇幅所限和耦合模型与普通离散元模型在这些方面并无本质区别，所以不在列出，可以参考文献[9]。

    致 谢
    本文的研究工作是在国家"九五"攻关课题的资助下完成的，得到了国家电力公司成都勘测设计研究院肖百云和王仁坤两位总工的支持和帮助，清华大学水利水电工程系张楚汉教授和徐艳杰博士也给予了大量的帮助，在此一并表示感谢。

参考文献
[1] Cundall PA, Hart RD. Development of Generalized 2-D and 3-D Distinct Element Programs for Modelling Jointed Rock[R]. ITASCA Consulting Group，Misc. Paper SL-85-1，1985.
[2] 鲁军, 张楚汉, 王光纶, 金峰. 岩体动静力稳定分析的三维离散单元法[J]. 清华大学学报，1996, (10).
[3] Zhang Chuhan, O. A. Pekau, Jin Feng and Wang Guanglun. Application of distinct element method in dynamic analysis of high rock slopes and blocky structures.[J] Soil Dyn. Earthq. Eng. 1997,(12).
[4] 张楚汉. 结构－地基动力相互作用问题[C]. 结构与介质相互作用理论及其应用，南京，河海大学出版社, 1993.
[5] Niwa Y, Kobayashi S and Azuma N. Application of Integral Equation Method to Some Geomechanical Problems[C]. Numerical Method in Geomechanics, 120-131，ASCE，New York，1976.
[6] Mansur WJ and Brebbia, CA. Topics in Boundary Element Research[C]，Chap.4，Springer-Verlag World Publishing Company，87-123，1985.
[7] 金峰，张楚汉，王光纶. 有阻尼的时域边界元方法[J]. 力学学报, 1997,29(15).
[8] 任允涛，各向同性与各向异性介质波动问题边界元法及其工程应用[D],清华大学博士论文,1996.
[9] 贾伟伟，离散元-边界元动力耦合模型研究及其工程应用[D],清华大学硕士论文，1999.

Coupling model of Distinct Element-Boundary Element
Jin Feng1  Wang Guanglun1   Jia Weiwei1
(1. Department of Hydraulic Engineering, Tsinghua University, Beijing, 100084)
Abstract
 A coupling model of 2-dimensional deformable block distinct element method and time domain boundary element method is presented. This model can simulate the static and dynamic responses of discontinuous rock and the effects of infinite domain, simultaneously. It can be used in static and dynamic stability and/or deformation analysis of rocky slopes and/or underground structures, especially when wave propagation input of earthquake and radiation damping must be considered. The application field of distinct element method is developed. The analysis result of benchmark problems shows that the precision of presented method is high.