离散单元法及其工程应用进展研究

张晨辰，程 静

（浙江同济科技职业学院，浙江 杭州 311231）

离散单元法及其工程应用进展研究

张晨辰，程 静

（浙江同济科技职业学院，浙江 杭州 311231）

离散单元法是20世纪70年代发展起来的一种解决非连续介质问题的数值计算方法。分析了离散单元法的基本原理及其在国内外的研究进展，探讨了离散单元法在岩土工程、结构工程以及动力分析等方面的应用。

离散单元法；基本原理；研究进展；工程应用

1 离散单元法的原理及发展

离散单元法 （Discrete Element Method，又称Distinct Element Method，简称DEM）是1971年Cundall为解决岩石系统的大规模运动提出来的一种解决非连续介质问题的数值计算方法［1～2］。DEM把岩体视为被节理切割而成的若干个块体的组合体，认为岩体的变形主要依赖于软弱结构面 （如裂隙、节理及断层等），并假定岩块为刚性，以刚性元及其周界的几何、运动和本构方程为基础，采用动态松弛迭代格式，建立了求解节理岩块非连续介质大变形的差分方程，以此来模拟岩体的断裂和破坏过程。它特别适合研究在准静力或动力条件下的节理系统或块体集合的力学问题［3］。它在边坡、基础、巷道稳定以及放矿力学研究等各方面得到了应用。

1.1 离散单元法的基本原理

离散单元法的基本原理与有限单元法有些类似，也是将整个研究区域划分成许多微小的单元，再进行求解。它与有限元法的区别是：单元因受不连续面控制，在以后的运动过程中，单元结点可以分离。即一个单元与其邻近单元可以接触，也可以分开。单元之间相互作用的力可以根据力和位移的关系求出，个别单元的运动则完全根据该单元所受的不平衡力和不平衡力矩的大小按牛顿运动定律确定［2～3］。简而言之，离散单元法的基本思想是，结合不同本构关系的牛顿第二定律。因而，可以采用动态松弛法或静态松弛法进行求解。离散单元法的原理虽然比较简单，但在解决非连续介质大变形问题时却非常实用。

离散单元法是一种按显式求解的不连续数值计算方法。该方法与在时域中进行的其他显式计算相似。由于显式法不需要形成矩阵，所以，它可以在考虑大的位移和非线性的同时，不花费额外的计算时间。假定在每一迭代时步内，每个块体单元仅对其相邻的块体单元产生力的影响，但是为了使显式算法稳定，时步就需要取的足够小。其基本方程有两组，一组是运动方程（牛顿运动方程），另一组是物理方程（力和位移关系）。在进行计算时，先利用中心差分法将运动方程进行离散，然后在每个时步△t内进行一次迭代。根据前一次迭代所得到的单元位置，由物理方程求出接触力，作为下一次迭代的出发点，再用来求出单元的新位置。如此反复进行迭代，直至最后达到平衡状态为止。最后的解一般表现为各个单元（或整个结构）的时间-位移（变形）关系。

离散单元法将岩体看成是由断层、节理、裂隙等结构面切割而成的、一个个刚性的、或者可以变形的块体。块体与块体之间通过角、面或者边进行接触。块体可以平移、转动或者变形；节理面可以被压缩、分离、滑动。所有块体镶嵌排列。在某一时刻，当给定块体一个外力或者边界位移约束时，各个块体就会产生力和力矩的作用。由牛顿第二定律可以得到各个块体的加速度，然后对时间进行积分，就可以依次求出块体的速度、位移，最后得到块体的变形量。在位移矢量的方向，块体会发生调整。这样，又会产生力和力矩的作用。如此循环，直到所有块体达到一种平衡状态或者处于某种运动状态之下。因此，离散单元法比较适合于模拟节理系统或者离散颗粒组合体在准静态或者动态下的变形过程。

1.2 离散单元法的研究进展

最初的离散单元法是一种二维的角－边（面）接触模型，用来分析准静力或者动力条件下岩石边坡的运动。由于其允许单元之间改变原有的接触关系，可以很有效地描述岩体等不连续介质以及颗粒散体，甚至对接触面有很大的变形、位移脱落的情况都很有效［3］。1978年，Cundall和Strack合作开发了二维圆形块体的BALL［4］程序，主要是为了研究颗粒介质的力学行为，同时，也证明了离散单元法是研究散体颗粒本构关系的有效工具。BALL程序随后在颗粒介质本构方程的研究中起到了不可替代的作用。1980年，Cundall开始研究块体受力后的变形以及根据破坏准则允许断裂的离散单元法［5］。这显然是对将块体视为刚体的离散单元法的一个进步，Cundall称这种方法为UDEC［6］。UDEC现已广泛用于岩土力学和采矿工程，被公认为对节理岩体进行数值模拟的一种行之有效的方法。1983年，Lemos开发了离散单元法与边界单元法耦合的半平面程序［7］，并用于计算节理和断裂介质中的应力分布问题。同年，Dowding［7］等在离散单元法二维接触模型的基础上，提出了边－边（面－面）接触模型。这一模型的特点是：在单元接触关系不变的前提下，允许单元之间分离，不必检索单元的接触关系。这样，不仅提高了计算效率，也避免了角—边接触模型的“锁定状态”。同时，他们还开发出基于刚性块体理论的二维离散元程序，即NURBM－2D，用于分析岩洞、隧道等地下岩体工程的稳定性问题以及地震荷载下的运动情况。1984年，Lorig开发了包括前处理和后处理的离散单元法与边界单元法耦合程序［8］。1985年，他又修改了原先的程序。修改后的程序文本称为HYDEBE，功能更强，包括一个前处理程序CREATE（类似于有限元程序中的自动划分网格），一个与边界单元法耦合程序BOUND和一个离散单元法程序BLOCK。1988年，Gilbert提出了面－面接触模型的二维离散元基本方程和算法［7］，同时开发了FORTRAN文本的二维程序NURBM－3D。日本学者Kawai等提出了固体力学模型中新的离散模型。该模型中的刚性单元用弹簧单元相连接，可用来进行一维和二维的任意边界条件下的裂隙附近弹塑性应力场分析。

众所周知，当将研究对象由二维转到三维时，会使数据结构变得更为复杂，这不仅要求计算机具有较大的容量，而且对计算结果的图形显示也提出了更高的要求。所以，三维离散单元法程序发展要滞后一些。三维离散单元法程序3DEC是Cundall与ITASCA咨询集团于1986年开发出来的，其基本原理同UDEC一样，但在数据结构方面作了较大的改进［6］。随着计算机技术的飞速发展，三维离散元程序已经比较多了，并且在实际工程中的应用也取得了较好的成果。

离散单元法在我国的研究和应用大概始于20世纪80年代中期。1986年，在第一届全国岩石力学数值计算及模型试验讨论会上，王泳嘉和剑万禧首次向我国岩石力学与工程界介绍了离散单元法的基本原理，并举了几个应用例子［9］。自此，与离散单元法 （主要以块体单元为主）相关的研究论文相继出现。

在离散单元法的应用中，边坡稳定分析的研究是最多的，有对山体滑坡这种典型的非连续大位移变形的应用，也有用其评价大坝滑坡体的稳定性，还有运用离散元数值模拟方法研究黄土坡滑坡与深部地质环境的关系以及评价滑坡体稳定性等等［10］。

离散单元法原先是为研究节理岩体的边坡稳定和巷道稳定而开发的，我国学者充分利用离散单元法中个别块体可以脱离母体而冒落的特点，将其应用于放矿的数值模拟和自然崩落法崩落机制及底部结构的稳定性研究中，并取得了很好的成绩。同时，应用离散单元法对不连续介质的研究也比较多，用其求解在港口和海岸工程中的堆石结构的非连续介质问题（如重力式码头后方抛石体的土压力问题、在波浪荷载作用下堆石防波堤、护岸的稳定性问题），都取得了令人满意的结果。

姚建国将离散单元法用于岩层移动的研究，认为这是一种比较符合矿山工程实际的数值方法［11］。刘建武根据Stewar和Brown提出的静态松弛法，在我国首先开发出了静态松弛离散元法程序，并结合二滩水电站模型试验资料进行计算，结果比较满意。1998年，王泳嘉和北京软脑软件有限公司合作开发了名为2D－BLOCK的离散元软件［12］，填补了我国的离散元相关软件的空白。随后，北京大学、清华大学、中国科技大学、东北大学、河海大学等著名大学和中国科学院力学研究所、中国科学院武汉岩土力学研究所、中国铁道研究院等著名科学研究部门均有人从事离散单元法的研究和应用工作，成果颇丰［13］。目前，离散单元法已在我国的采矿工程、岩土工程以及水利水电工程等得到广泛应用。

2 离散单元法在工程实施中的应用

鉴于离散单元法在解决非连续介质问题方面的独特优势，一经提出，很快在边坡、基础、巷道稳定，放矿力学研究等方面得到了应用。与此同时，许多学者注意到了离散单元法的独特长处，将其引入到结构工程领域的分析研究工作中，并取得了良好的效果。下面就离散单元法的在岩土工程、结构工程以及动力分析中的应用进行探讨。

2.1 边坡工程

边坡的失稳破坏运动是一个动态过程，存在岩土体的滑动、平移、转动和岩体的断裂及松散等复杂过程，具有宏观上的不连续性和单个岩块体运动的随机性［13～15］。而离散单元法考虑岩土介质的各向异性、非连续性，允许单元介质发生大的滑移和岩体运动过程中块体的松散甚至空洞的存在，故比较适合做边坡稳定性分析。很多学者在这一方面做了大量的工作。沈宝堂等［14］通过模型试验连续拍摄了边坡失稳瞬间的照片，并将其与离散单元法的计算结果对照，结果十分吻合。这表明，离散单元法的数值模拟可以代替昂贵的相似材料模型试验。文献［15］将裂隙岩体渗流模型与离散元力学模型进行耦合，根据裂隙岩体主干裂隙与网络状裂隙的特点，建立了离散介质与连续介质耦合渗流模型，并将该模型用于岩质高边坡进行稳定性分析，取得了比较理想的成果。文献［16］用离散元法分析了陡峭岩体边坡在多种假想条件下的稳定性，并对加固系统的性能进行评估，将模拟结果与现场实测结果进行比较，两者基本吻合。有的学者将离散元与强度折减结合起来，确定岩质边坡的滑动面和安全系数，论证了这是一种可靠、有效的方法。还有学者用离散单元法研究了边坡加锚节理岩体的力学特性，所得结论与采用物理试验结果相符。这验证了用数值模拟方法进行岩石力学试验的可行性和可靠性。

2.2 地下工程

对于地下工程而言，岩体开挖以后，洞室周围一定范围内的岩块在节理、裂隙等结构面的影响下都会松散、滑落。对于这种状况，基于连续介质数值方法是没有办法模拟的，而离散单元法可以对这种现象进行比较好地模拟。有学者运用离散单元法对地下洞室群进行了稳定性分析，得出了可能产生的破坏块体，可以预见潜在的破坏发生规模以及具体的分布特征等等。文献［17］将 3DGIS与3DEC相结合，对复杂的地下工程施工、稳定与受力状态进行了可视化、透明化处理，为设计咨询提供了方便。文献［18］将UDEC用于喜马拉雅山下的一个大型洞室变形机制研究，探讨了不同输入参数，如节理间距、节理本构模型等对洞室变形行为和位移的影响。文献［19］将离散单元法应用于地下铁矿开采时的围岩稳定研究中，研究了不同开采方法对边坡和围岩的影响，以及围岩的变形规律，得出了有意义的结论。

2.3 结构工程

在20世纪80年代，日本学者Meguro K.和Hakuno M.在Cundall的离散元模型的基础上，提出了一种扩展离散单元模型，用于模拟钢筋混凝土等连续体结构的破坏和倒塌［17］。当时，他们用这种模型来模拟墨西哥地震中的房屋倒塌情况，结果与实际情况很相似；他们曾演示因某层结构出现问题而使房屋压缩掉一层的破坏形式，这种破坏形式在1995年1月17日的日本阪神地震中得到了验证。在国外，还有用离散单元法研究钢筋混凝土结构在冲击荷载作用下和在高应变率情况下的破坏（局部破坏）形式。比如，加拿大的F.V.Donze等人就用此模型来分析混凝土材料在高应变率情况下强度提高的原因。

在国内，把离散单元法应用在结构工程中的人不是很多（至少没有用在岩土工程中的人多），所研究的范围主要是在混凝土结构的动态破坏过程中，比如说在地震、爆炸等力作用下的破坏形式。刘凯欣等利用离散单元法对冲击荷载下混凝土平板变形、破坏的全过程进行了数值模拟和动画显示，并分析了不同强制位移速率下的不同破坏形式，验证了应变率效应对混凝土破坏机理的影响是很大的［19］；同济大学的秦东等还把此模型应用于桥梁结构中，分析桥墩等结构在地震作用下的倒塌过程［20］。

离散单元法虽然在岩土工程中的应用已经比较广泛，而且也有相应的系统软件，但其在结构工程中的应用并不多，可以说是刚处于起步阶段。所以，离散单元法在结构工程中的应用还需要进一步研究。

2.4 动力分析

岩体所受的力是复杂的，除了静力以外，还有地震和施工过程中的爆破所带来的动力。由于动力具有突然性、不稳定性，其对岩体的稳定性影响要大大超过静力所带来的影响。对于节理岩体，这种影响尤为明显。因此，对节理岩体的动力稳定性进行研究，给工程加固和管理一个科学的依据，是当前节理岩体研究中的一个热点问题。文献［21］研究了具有有限剪切强度的单一岩石节理对垂直入射的一维剪切波传播的影响，用UDEC模拟求得了透射、反射系数，所得结果与解析结果完全一致。文献［22］模拟了不同药量的节理岩体中集中药包的爆破过程，将数值模拟与试验数据进行了比较，吻合较好。因此，离散单元法可以比较好地模拟节理岩体动力响应问题。

3 结语

经过几十年的发展，离散单元法的相关理论有了很大的发展，算法也日趋成熟。因为它在运动、受力、变形3个要素都是假设，所以其理论还不完善。同时，离散单元法在算法以及前后处理上还需要一定的改进，以便被更广泛应用。

［1］Cundall P A.The Measurement and Analysis on Rock Slopes［D］.London：University of London，Imperial College of Science and Technology，1971.

［2］Cundall P A.A Computer Model for Simulating Progressive Large Scale Movement in Block Rock System［J］.Symposium ISRM，1971（2）：129-136.

［3］王泳嘉，邢纪波.离散单元法及其在岩土力学中的应用［M］.沈阳：东北工学院出版社，1991：110－125.

［4］方韬.离散单元法的研究及其在结构工程中的应用［D］.杭州：浙江大学，2004.

［5］Cundall P A.A Generalized Distinct Element Program for Modeling Jointed Rock［J］.U.S.Army，1980（3）：80.

［6］刘连峰，王泳嘉.三维节理岩体计算模型的建立［J］.岩石力学与工程学报，1997，16（2）：36－42.

［7］柏署，张磊，杨燕.离散单元法程序研制和应用［J］.西部探矿工程，2010（2）：8－12.

［8］Lorig L J.A Hybrid Computational Model for Excavation and Support［M］.Design in M innesota，1984.

［9］王泳嘉.离散单元法：一种适用于节理岩石力学分析的数值方法［C］／／中国岩石力学与工程学会.第一届全国岩石力学数值计算及模型试验讨论会论文集，中国岩石力学与工程学会1986：32－37.

［10］沈宝堂，王泳嘉.边坡破坏机制的离散单元法研究［J］.东北工学院学报，1989，10（4）：349－354.

［11］姚建国.开采影响下岩层移动变形的块体单元数值模拟［C］／／徐文焕.第一届全国计算岩土力学研讨会论文集（一），成都：西南交通大学出版社，1987：232－237.

［12］王泳嘉，宋文洲，赵艳娟.离散单元法软件系统2DBlock的现代化特点 ［J］.岩石力学与工程学报，2000（19）：1057－1060.

［13］周先齐，徐卫亚，钮新强，等.离散单元法研究进展及应用综述［J］.岩土力学，2007（S1）：400－416.

［14］王洪涛.裂隙网络渗流与离散元耦合分析充水岩质高边坡的稳定性［J］.水文地质工程地质，2000（2）：30－33.

［15］刘礼领，殷坤龙.离散单元法在水库库岸滑坡稳定性分析中的应用 ［J］.水文地质工程地质，2003（4）：65－68.

［16］王涛，陈晓玲，杨建.基于3DGIS和3DEC的地下洞室围岩稳定行研究 ［J］.岩石力学与工程学报，2005，24（19）：3476－3481.

［17］李承.基于离散单元法的钢筋混凝土框架结构爆破拆除计算机仿真分析［D］.上海：同济大学，2000.

［18］刘辉，陈文胜，冯夏庭，等.大冶铁矿露天转地下开采的离散元数值模拟研究 ［J］.岩土力学，2004，25（9）：1413－1417.

［19］刘凯欣，高凌天，郑文刚.混凝土动态破坏过程的数值模拟［J］.工程力学，2000（增刊）：470－474.

［20］秦东，范立础.钢筋混凝土结构倒塌全过程数值模拟［J］.同济大学学报：自然科学版，2001，29（6）：80－83.

［21］程东幸，潘炜，刘大安，等.锚固节理岩体等效力学参数二维离散元模拟 ［J］.岩土力学，2006，27（12）：2128－2132.

［22］郭易圆，李世海.离散元法在节理岩体爆破振动分析中的应用［J］.岩石力学与工程学报：自然科学版，2002，21（2）：2408－2412.

[责任编辑 杨明庆]

Research on Discrete Element M ethod and Its Engineering Application

ZHANG Chen-chen,CHENG Jing
(Zhejiang Tongji Vocational College of Science and Technology,Hangzhou 311231,Zhejiang,China)

Discrete element method is a numerical computation method for solving the non-continuum problem which developed in the 1970s.It analyzes the basic principles of discrete element method and its research progress at home and abroad as well,it also discusses its application in geotechnical engineering, structural engineering and dynamic analysis,etc.

Discrete element method;basic principles;research progress;engineering application

TU457

A

1008-486X（2015）02-0030-04

2015-01-28

张晨辰（1988-），女，浙江诸暨人，助教，硕士研究生，主要从事工程施工专业的教学与研究工作，研究方向：高坝结构及复杂地基。