数值模拟实验在岩土工程中的应用与展望

彭岩岩 刘宇航 王天佐 杜 伟 郑质彬

(1.绍兴文理学院 土木工程学院，浙江 绍兴312000；2.绍兴文理学院 岩土工程与地质灾害实验中心，浙江 绍兴312000)

岩土是一般材料，也是一种地质结构体，它具有非连续、非均质、非线性的特性及复杂的加卸载条件和边界条件，这使得岩土工程问题通常无法简单求解.因此数值模拟方法成了解决某些岩土工程问题的有效工具之一.从20世纪50年代开始，人们就利用数值模拟方法对岩土工程问题进行大量研究，经过60多年的发展，针对岩土工程问题的数值模拟方法逐渐成熟，成为岩土工程学科一个重要的研究方向，使复杂岩土工程问题的设计发生了根本性的变化.岩土工程数值模拟不仅取代了传统线弹性力学实验，而且也在岩石工程非线性实验中显示出极大的优势[1].

岩土工程的数值模拟实验是对岩土工程活动和自然环境变化过程中岩体及工程结构的力学行为进行数值模拟的一种手段.在进行岩土工程数值模拟的时候，需要对岩体进行分类：一是基于连续岩体力学的数值模拟，二是基于不连续岩体力学的数值模拟.有限元法、边界元法、有限差分法把岩体看作连续介质进行模拟，而离散元法、不连续变形分析法、数值流形法把岩体看作不连续介质进行模拟，本文会对上述方法和应用进行总结介绍.

1 数值模拟实验在岩土工程中的优势

岩体本身具有非均质、非连续、非线性及复杂的加卸载条件和边界条件的特点，加之岩体所处环境也比较复杂，岩体工程开挖前就受地应力、地下水、周围温度等耦合作用，所以很难建立完善的地质力学模型[2].

1.1 传统研究方法的局限性

由于岩体的不连续性与非均质性，实验室岩石试样具有明显的“尺寸效应”，并且通过实验得到的结论往往与实际工程相差甚大[3].

通过现场原位实验的方式取得数据，也非常艰难.岩体力学原位实验一般耗资巨大，且受地形、地质、施工条件限制，得到的实验结果不具代表性，难以推广到其他工程[2].

1.2 数值模拟在岩土工程中的优越性

数值模拟在岩土工程中的适用范围非常广，并且节约资金，它不仅能模拟岩体复杂结构特性，还能研究岩土工程活动对周围环境的影响，并对工程灾害进行预报.通过对现场原位实验的实测与反分析，可以获得节理岩体的等效力学模型，逐步取代成本较高的原位实验，加快工程进度，且可应用于各种地形、地质与施工条件，推广实验结果的应用范围与使用条件[3].因此用数值模拟的方法来解决岩土工程问题是行之有效的.

2 岩土工程中常用的数值模拟方法

2.1 有限元法

2.1.1 有限元法原理

有限元法是利用变分的原理去求解数学物理问题的一种数值模拟方法.有限元法最早由布理克(W.Blake)在1966年引入岩土工程领域，用来解决岩土工程问题.有限元法基于最小总势能原理通过解方程组的方法来求解，是目前岩土工程领域中应用最广泛的数值模拟方法.

有限元法是用多个彼此相联系的单元体所组成的近似等价物理模型来代替实际的结构或者连续物体，通过结构及连续体力学的基本原理及单元的物理特性建立起表征力和位移的关系去建立方程组，解方程求其基本未知物理量，并由此求得各个单元的应力、应变及其他辅助值[4].

2.1.2 有限元法的应用

有限元法由线性发展到非线性和大变形问题的应用(二维发展到三维)，目前还可考虑流变、温度与应力场耦合，损伤、渗流、断裂以及波动和动力效应[5].

刘庭金等[6]利用有限元分析矿山、地铁等地下工程由于洞室开挖引起的围岩卸载过程中，洞室孔壁围岩附近发生的损伤演化和应力场调整全过程进行分析.郑颖人等[7]对有限元强度折减法的计算精度和影响因素进行了详细分析，包括屈服准则、流动法则、有限元模型本身以及计算参数对安全系数计算精度的影响，并给出了提高计算精度的具体措施.应用于岩质边坡的稳定分析，得到了岩质边坡的滑动面和安全系数，开创了求节理岩质边坡滑动面与稳定安全系数的先例.

2.2 边界元法

2.2.1 边界元法原理

边界元法同有限差分法和有限元法一样，都是一种用来解决边值问题的数值分析方法，它可以用来解决弹性力学、塑性力学，以及热传导、地下水力学等方面的问题，发展历史悠久， 但直到20世纪60年代后期在计算机技术得到发展时， 边界元法才成为实际可行的一种数值模拟方法[8].由于岩土工程问题的复杂性，边界元法在1976年才引入到岩土工程中.

边界元法是通过求解边界积分的方法来求解边值问题，在边界元上划分单元，求边界积分方程的解，进而求出区域内任意点的场变量，所以边界元法也称边界积分方程法.边界元法又分以互等功原理为基础的直接法，和以叠加原理建立起来的间接法.

2.2.2 边界元法的应用

边界元法和有限元法比起来，可以用降维的方法来简化计算(三维问题二维化，二维问题一维化)，不但计算起来方便，而且计算精度高，但是面对非连续、非线性介质问题边界元法则比较难适应.目前边界元法主要在地下工程开挖、土体结构相互作用及地下水流动过程的一般应力和变形分析有着应用.虽然边界元法的适用范围有限，但是和其他数值方法的联合使用能充分发挥其优越性，为解决岩土工程问题开辟了新的途径.例如，在线弹性区域或无限域、半无限域可以采用边界元法，在非线性的区域采用有限元法，发挥两种算法各自的优势，使计算效率及精度得到提高，对工程实际应用有很大的帮助[9].

马天寿等[10]用边界元法对页岩地层井眼坍塌问题进行了分析并得出弹性模量各向异性、水平地应力差异和钻井液密度等对井壁应力分布影响较大，而泊松比各向异性的影响较小的结论，如图1.

图1 边界元方法解出的井周应力分布图[10]

2.3 有限差分法

2.3.1 有限差分法原理

有限差分方法是以最小势能原理，通过解方程的方式进行求解.这种方法是一种最古老的求解方程组的数值方法，在计算机出现以前一般的手摇计算器也可求解.20世纪80年代末由美国ITASCA公司开发的FLAC程序广泛采用差分方法进行求解，并且在岩土工程数值计算中得到了广泛应用[5].

2.3.2 FLAC3D的应用

鉴于有限差分法单独在岩土工程中的应用并不多，一般都是基于有限差分法的FLAC程序进行岩土工程计算，所以这里讲的是FLAC程序的应用.FLAC程序可用来模拟地质材料的大变形、失稳、动力、流变、支护、建造及开挖等问题，同时还可以模拟渗流场和温度场对岩土工程的影响.

李为腾等[11]解决了 FLAC3D中 CABLE 单元无法实现锚杆(索)破断失效的问题，并采用Fish语言编程，将修正模型嵌入到 FLAC3D主程序中，实现锚杆破断失效的“单元化”，如图2.

图2 FLAC模拟修正影响系数随地应力变化曲线[11]

2.4 离散元法

2.4.1 离散元法原理

离散元法是Cundall在1971年所提出来的，后经Voegele等人的发展， 成为一种新的数值模拟方法[12].离散元法是以牛顿运动定律的显示求解的数值方法，离散元法也要将区域划分为单元，但是单元因受节理、劈理等不连续面的控制，在以后的运动过程中，单元节点可以分离，即一个单元与其邻近单元可以接触，也可以分开.单元之间相互作用的力可以根据力和位移的关系求出，而个别单元的运动则完全可以根据单元所受的不平衡力和不平衡力矩的大小按牛顿运动定律确定[9].

2.4.2 离散元法的应用

离散单元法可分为动态松弛法和静态松弛法两种.目前常用的大多是动态松弛法.动态松弛法用解决动力学问题的方法去求解非线性静力学问题，用显式中心差分法近似对运动方程进行积分计算，并假设块体在运动时将动能转化成热能耗散掉，把人工黏性阻尼引入计算，使系统达到平衡状态、运动趋于稳定[13].

常晓林等[14]用离散元的方法模拟了岩石介质从小变形到大变形再到破坏的全过程.李晓柱等[15]用离散元的方法分析了堆石坝现场碾压实验，验证离散元数值模拟方法应用于堆石坝碾压特性研究的可行性，更直观地从细观角度解释堆石体碾压过程中宏观参数(如干密度等)的变化规律，为大坝I区堆石料选取科学合理的碾压施工参数，为堆石体碾压特性研究提供新的途径.王贵君[16]针对国内某高速公路隧道工程，应用离散单元法对节理裂隙岩体中不同埋深无支护暗挖隧洞的稳定性及其机理进行了数值模拟.

2.5 不连续变形分析法(DDA法)

2.5.1 不连续变形分析法原理

不连续变形分析方法一般简称DDA法(Discontinuous Deformation Analysis)，由石根华首创，基于岩体介质非连续性发展起来的，以模拟复杂加载条件下离散块体系统不连续大变形的力学行为为目的的一种数值方法[1]，该方法以最小势能原理为基础，通过解方程的方式求解问题，适用于发生大变形，岩体发生非连续破坏的情况.

石根华在1988至1989年期间开发了二维DDA程序，用来进行前处理、正分析、反分析和后处理的计算，程序包括DDACUT，DDA FORWARD，DDA BACKWARD和DDAGRAPH，4个部分[1].

2.5.2 不连续变形分析法的应用

DDA法在岩土工程当中主要应用于两个方面：一是在爆破工程研究中的应用，二是在地下岩体开挖工程中的应用.

张丽娟等[17]采用DDA 方法对台阶爆破抛掷的全过程进行了模拟分析， 给出了人工边界、爆破荷载的处理算法.邬爱清等[18]利用DDA的方法分析了块体稳定性验证及其在岩质边坡性分析中的应用，如图3.吴建宏等[19]将DDA法应用在岩石边坡失稳数值仿真当中.

图3 DDA计算后块体系统变形与破坏[18]

2.6 数值流形法

2.6.1 数值流形法原理

在1995年，石根华提出数值流形法(NMM)，数值流形法是利用现代数学中“流形”的有限覆盖技术建立起来的一种新的数值计算方法，将有限元、不连续变形分析(DDA)和解析方法统一到一种计算方法中，它吸收了有限元、DDA 和解析法各自的优点，通过分片光滑的覆盖函数，对连续和非连续问题统一了计算格式，是一种十分适合于岩土工程分析的数值方法[13].

2.6.2 数值流形法的应用

数值流形法可同时处理非连续问题与连续问题.但由于网格的连接与单元划分的限制，数值流形法在开裂计算方面仍存在一定的困难.目前研究主要集中在连续与非连续问题的求解和裂纹扩展的模拟[20].

李树忱等[21]充分利用数值流形方法中两套独立网格的特点，采用围线积分法来计算应力强度因子和最大周向应力确定裂纹的扩展角，模拟裂纹的扩展过程.钱莹[22]等利用三角形元素来构成数学网格，利用流变法来分析爆破振动对边坡动态稳定性的影响，如图4.

图4 流形元模型网格[22]

3 数值模拟实验存在的问题及发展

尽管数值模拟方法取得了很多研究成果，但从当前数值模拟在岩土工程中的应用来看，主要存以下几个问题：

(1)参数选取的问题.无论哪种数值模拟方法都必须准确选定岩石或岩体的物理力学参数.但由于岩体本身与所处环境的复杂性，确定这些参数并非易事[1]，因此得到的数据和工程实际也有出入.

(2)计算机储存量不足和计算速度受限的问题[23].由于越来越多大型工程的兴建，当前岩土工程的数值模拟迫切需要发展并行计算方法，用多核或多联计算机将原有的程序做重构，以提高计算速度.

(3)岩体进入破坏模式后力学机制发生转变，数值模拟中的相关破坏准则难以准确描述破坏后的岩体力学行为[24].

岩土工程中的数值模拟方法有很多，每种模拟方法都有其优点，但也各局限性.要想用单一的方法解决岩土工程问题是不现实的，为了更好地发挥每一种方法的长处，多种数值模拟方法综合应用成了近年来岩土工程数值模拟实验的新趋势.如在围岩影响区采用非连续分析法，在围岩影响区之外的原岩应力区用连续介质分析法.但是要解决两种不同算法在围岩影响区和原岩应力区交界面的位移和应力的协调问题.遗憾的是，目前专注于上述问题的学者有限，尚未有突破性进展.

4 总 结

本文介绍了数值模拟实验在岩土工程中的研究进展，列举了先进的数值模拟实验方法，并分别介绍了相关原理及在岩土工程中的应用，最后指出了目前数值模拟实验存在的问题及其发展方向.总而言之，随着岩土工程规模日益扩大，对岩土工程建设的科研设计水平和建设精度的要求越来越高，用数值模拟去解决岩土工程问题是一种有效手段.数值模拟实验对掌握岩土工程围岩应力变化规律和变形破坏规律有重要意义，能够对实际工程的支护设计提供理论依据.