岩石力学数值计算方法研究

高宏 庞毅力

摘要：对几种常用岩石力学数值计算方法的优越性和局限性进行探讨，分析了岩石力学数值计算方法存在的瓶颈和当前解决这些问题的途径。

关键词：数值方法，岩石力学

中图分类号：O3文献标识码： A 文章编号：

1 引言

当前岩石力学数值计算方法得到迅猛发展，出现了有限差分、有限元、边界元、离散元、块体元、无限元、流形元及其混合应用等各种数值模拟技术，使复杂岩石力学工程问题的设计发生了根本性的变化[1]。不同数值计算方法的结合,更能发挥各种数值方法优势互补的作用，如有限元—边界元的混合、有限元—离散元的混合、有限元—无限元和有限元—块体元的混合采用等。然而，由于岩体具有非连续、非均质、各向异性、天然初始地应力影响、地下水影响及复杂边界条件处理等诸多复杂性使得当前岩石力学数值计算仍然是一个值得探讨的问题。

2 常用岩石力学数值计算方法应用分析

2. 1 有限元法

1966 年,布理克[2] (W. Blake)最先应用有限元法解决地下工程岩石力学问题。目前，在岩石力学数值计算方面，有限元法主要用来求解线弹性、弹塑性、粘弹塑性、粘塑性等问题，是地下工程岩体应力-应变分析最常用的方法。其优点是可以部分地考虑地下结构岩体的非均质和不连续性，可以给出岩体的应力、变形大小和分布，并可近似地依据应力、应变规律去分析地下结构的变形破坏机制。

2. 2 边界元法

边界元法在20 世纪70年代中期得到迅速发展，在处理半无限域或无限域问题方面非常方便，适用于解决岩石力学问题尤其是岩石力学中地下开采的有关问题。该法只在求解区域的边界上进行离散(剖分单元)，这样就把考虑问题的维数降低了一维，这也是边界元法的优点。例如,在线弹性区域或无限域、半无限域采用边界元法，在非线性的区域采用有限元法，充分发挥各自的优势，使计算效率、计算精度得到提高和改进，这对工程实际应用是很有意义的[2 ,4]。王泳嘉[4]等人讨论了边界元的应力不连续法和直接边界积分法，并用应力不连续法给出了位移不连续时的解。

2. 3 离散单元法

离散单元法(Distinct Element Method)是20世纪70年代后发展起来的一种用于节理岩体应力分析的数值方法，特别适用于节理岩体及其与锚杆(索)的应力分析。该方法以结构面切割而成的离散体为基本单元，其几何形状取决于岩体结构中不连续面的空间位置及其产状，应用牛顿运动定律描述各块体的运动过程，块体可以发生有限移动与转动，体现了变形和应力的不连续性。沈宝堂等学者通过两个模型实验的结果，并与离散元法数值模拟的结果相比较，验证了离散元法用于边坡稳定性分析的可行性[5]。郭爱民[6]等人利用离散单元法研究矿山边坡的破坏机理，笪盍等人[7]利用离散单元法对盘石镍矿边坡进行稳定性分析和计算，计算结果与现场岩体位移进行比较表明两种结果吻合较好。

2. 4 关键块理论

关键块理论KB T ( Key Block Theory)是在1985年首先由Goodman 教授和石根华博士提出并用于工程稳定性分析的另一种数值计算方法。块体理论认为，在开挖面上所揭露的块体，可以分为可能产生向开挖面运动的块体和不可能向开挖面运动的块体。不可能向开挖面运动的块体即为稳定块体。石根华[8]等人提出搜寻关键块体一般方法，并介绍了传统关键块体理论的近期发展、应用和局限。

2. 5DDA法

不连续变形分析方法DDA (Discontinuous Deformation Analysis)是由石根华博士首创，基于岩体介质非连续性发展起来的，以模拟复杂加载条件下离散块体系统不连续大变形的力学行为为目的的平行于有限元法的一种数值方法，与有限元法不同之处是可以计算不连续面的位错、滑移、开裂和旋转等大位移的静力和动力问题。孙亚东[9]等人利用DDA法分析岩质边坡的倾倒破坏，并与Goodman 和Bray提出的基于极限平衡原理的分析方法进行比较。邬爱清、丁秀丽等学者应用DDA法研究了复杂地质条件下地下厂房围岩的变形与破坏特征[10]。

2. 6 FLAC方法

Cundall根据有限差分法原理,提出了FLAC( FastLagrangian Analysis of Continuum)分析方法。该方法采用了混合离散方法、动态松弛方法和显示差分方法，不形成刚度矩阵。它的求解方法虽同离散元法的显式按时步迭代求解，但是结点的位移连续，本质上仍属于求解连续介质范畴的方法。

2. 7 数值流形法

数值流形法是利用现代数学—“流形”的有限覆盖技术建立起来的一种新的数值方法，将有限元、不连续变形分析(DDA)和解析方法统一到一种计算方法中，它吸收了有限元、DDA和解析法等的优点，通过采用分片光滑的覆盖函数，对连续和非连续问题建立了统一的计算格式，是一种十分适合于岩土工程分析的数值方法[12] 。它最早由石根华博士在1991年提出并率先应用在块体与节理岩体的变形模拟中，在国内外学者的共同努力下，二维数值流形方法已经拥有一套比较完善的理论，而且应用方面的探讨也日渐增多[13]。郑榕明、张勇慧等也曾经对DDA法作了大量的研究工作[14 ]，在此基础上发展了二维流行元程序，并应用在地下洞室开挖的模拟中[15]，刘红岩等人[16]利用数值流形方法对一层状岩石边坡的倾倒破坏过程进行了模拟。

3岩石力学数值方法应用中的几个瓶颈与展望

3. 1 计算参数的取值问题

由于岩体性态与环境的复杂性,准确确定这些参数并非易事,因而数值分析手段至今仍不能最终为工程设计和工程决策提供可靠依据,幸运的是至今仍火热的反演分析方法有望在这方面为原始参数的取值提供一种新的途径[17]。中国近年来在反分析方面进行了大量研究工作，已由简单的线弹性反演问题发展到非线性、粘弹塑性反分析，从单一的毛洞围岩到考虑支护结构体系的反演,从有限元位移反演到边界元位移反演，从确定性反演到非确定性随机反演等等。而且反分析的目的已不仅仅是得到模型参数，更重要的是应用这些参数进行相应的时间序列值分析以及从参数估计发展到模型识别进而建立新的模型，以便对工程效果做出更合理的评价和有依据的预测。联邦德国对计算参数的取值问题也十分重视，他们认为计算输入的参数必须源于客观实际,地质勘探、岩石力学和数值计算必须紧密地结合起来，甚至从事计算的人需要自己动手到现场去取得第一手资料，而不只是单纯依靠委托单位提供的参数，这一点同主观臆断假定参数或依靠委托单位提供参数而不深入实际研究参数可靠性等脱离实际的作法形成鲜明的对比[18]，是我们应该值得学习和重视的。

3. 2 本构关系的选取问题

从事数值模拟计算的人都知道岩体性质比迄今为止人们所熟知的任何工程材料都要复杂得多，它既非连续介质,又非松散介质，既非弹性体，又非塑性、粘性体，从而导致计算中采用的本构关系很难准确把握。尽管用数值模拟对岩体结构进行力学分析的方法得到了广泛的应用，并且取得了许多成果，但是不敢断言这种方法在将来是否会对这样一类问题的研究有新的突破，至少目前还不可能将这一类问题的研究提高到一个全新的高度[19]。基于这样的原因，人们也在力求寻找其它的弥补途径，有学者改变单一的确定性正向分析方法，适应岩体的非确定性特征，将数值模拟方法与反分析方法、随机方法、模糊方法、人工智能等结合起来，这或许是数值模拟方法的一个正确发展方向[19]。

3. 3 单一数值方法的局限性问题

为了吸取各种计算方法的长处以弥补其不足,近年来涌现出大量的各种数值方法的耦合计算,这种思路进一步反映了岩体工程的计算特点，清华大学的周维恒先生在1993年就断言这种思路对岩石力学数值计算的发展是十分有益的[20] 。离散元、块体元和有限元、边界元、无穷元之间的结合又可提出若于种耦合方法，以发挥离散元和块体元在模拟不连续岩体方面的长处，并利用有限元、边界元、无穷元在模拟连续介质方面成熟的理论和计算技术，使应力分析、破坏、垮落及运动整个过程的数值模拟有可能实现。数值方法同解析方法或半解析方法的结合则是另一条可行的途径，这种结合的特点是用相应的数学推导得到更精确的(也更复杂了)单元模式，再通过离散化求解，解题效率及精度提高，不足之处是应用的局限性也随之增加。何翔[19]等人在Biot固结方程的基础上,引入介质渗透系数张量随应力张量的变化函数,建立能反映介质中应力场与渗流场非线性耦合作用的微分方程组,并在此基础上进行渗流场与应力场耦合问题的有限元求解,采用了精细积分方法进行时间离散。

参考文献:

[1] 谢和平,刘夕才,王金安,等. 关于21 世纪岩石力学发展战略的思考[J].岩土工程学报,1996 , 18 (4) : 98 - 102.

[2] 王泳嘉. 边界元法及其在岩石力学中的应用[J].东北工学院学报, 1984 (1) : 1 - 12.

[3] 刘红岩,秦四清. 层状岩石边坡倾倒破坏过程的数值流形方法模拟[J].水文地质工程地质, 2006 (5) :22 - 24.

[4] 王泳嘉, 冯夏庭. 关于计算岩石力学发展的几点思考[J].岩土工程学报. 1996 ,18 (4) :103 - 104.

[5] 刁心宏, 冯夏庭, 杨成祥,等. 岩石工程中数值模拟的关键问题及其发展[J].金属矿山, 1999 (6) :5 - 7.

[6] 周维坦. 岩体力学数值计算方法的现状与展望[J] . 岩石力学与工程学报, 1993 , 12 (1) :84 - 88.

注：文章内所有公式及图表请用PDF形式查看。