基于能量耗散的脆性岩体张开位移数值模拟

郑文华，朱维申，刘德军

（1. 中国建筑科学研究院 建筑安全与环境国家重点实验室，北京 100013；2. 山东大学 岩土与结构工程研究中心，济南 250061；3. 同济大学 地下建筑与工程系，上海 200092）

1 引 言

从2009年4月20日的第2届水力发电技术国际会议获悉，截至2008年底，我国水电装机容量达到1.72亿千瓦，水电能源开发利用率已达到27%。因此，在我国越来越多的大型水利水电工程进入施工阶段。

大量的现场实测数据表明，不少地下厂房的开挖过程中都会产生劈裂裂缝[1－3]及其扩展现象（如图1所示），且裂缝张开和扩展所引起的变形量在实测位移值中占有较大比重。在国外，文献[4]中对16个日本大型洞室进行的研究表明，洞周一倍宽度的范围内，在岩浆岩和沉积岩中往往都存在显著的空腔性应变，并且，该应变在总应变中占了相当大的比例，如表1所示。在国内，彭琦等[5]对地下厂房围岩变形特征进行分析，以施工期变形监测资料为基础，结合地质和施工资料，对地下厂房的围岩变形特征及其机制进行分析。得出围岩的变形主要由控制性结构面的“张开位移”构成，占全部变形的84%～92%。

图1 现场观测的劈裂裂缝现象Fig.1 Crack phenomena of field observation

表1 日本大型地下洞室张开位移的监测结果Table 1 Monitoring results of opening displacement of large caverns in Japan

研究高地应力脆性岩体在开挖过程中裂隙张开和扩展所产生的张开位移，有助于更加精确地分析研究开挖引起的洞周围岩变形，对分析地下洞室群的稳定性有着十分重要的意义。

国内外有关地下洞室围岩变形的研究比较多，其方法有连续体力学方法[6]、非连续体力学方法[7－10]以及现场测试方法等。非连续体力学方法，虽然可以考虑非连续的变形，但有关张开位移的研究尚处于起步阶段。现场测试方法可以测得包含张开位移在内的总位移，例如钻孔电视[11]、多点位移计[12]等测量手段。但现场测试方法仅限于场地条件便于安装测量设备的情况下，且该类方法测量范围有限。本文从能量耗散的角度出发，提出考虑能量耗散作用的变弹性模量弹脆性模型，并通过二次开发把该模型程序化。将该模型应用到锦屏一级水电站的实际开挖模拟中，比较了位移现场实测结果与数值模拟结果，取得了较好的效果，其结果对类似工程提供了一定的参考价值。

2 裂纹扩展的能量变化

2.1 Griffith理论

Griffith曾指出，裂纹扩展时会释放弹性能，同时，在新形成的裂纹表面又会产生新的表面能增加。当释放的弹性能等于增加的表面能时，裂纹处于稳定状态，停止扩展；而当释放的能量大于增加的表面能时，裂纹处于不稳定状态，会发生进一步的扩展现象。

同时，Griffith[13]给出了释放的弹性能计算公式：

式中：c为裂纹长度的一半；E为弹性模量；σ为作用在裂纹长轴投影面上的平均应力，如图2所示。

图2 模型受力示意图Fig.2 The stress diagram of model

2.2 变弹模弹脆性模型

一般而言，脆性岩体最明显的特征就是当岩体受到的应力达到其峰值强度时，随着变形的进一步增加，其应力值会迅速降低到一个比较低的水平。具体来说就是应力-应变关系曲线有一个很陡的应力跌落现象[14－15]。

对于弹脆性模型来说，考虑能量耗散作用下卸荷状态的弹性模量变化时，认为若单元处于弹性阶段则按弹性卸载，若单元超过最大强度值则发生脆性破坏，此时，岩体强度降低。根据以往学者进行的脆性岩体室内试验的结果[16－18]，可将弹脆性模型的加、卸载过程简化，如图3所示。

图3 考虑能量耗散作用的弹脆性应力应变关系曲线Fig.3 Stress-strain curve of elastic-brittle model considering energy dissipation

因此，考虑能量耗散作用时，有

基于上述理论可以得出，考虑能量耗散作用的弹脆性模型的关键点在于加、卸载阶段的弹性模量不一样，卸载阶段的弹性模量变为等效弹模而应力超过岩体最大强度值时，岩体发生脆性破坏，岩体强度降低，此时可近似认为，按极限强度时的弹模进行卸载。

3 FLAC3D计算模型的二次开发

3.1 二次开发平台

FLAC3D[19]是快速拉格朗日分析方法，它采用显式或隐式有限差分，可以模拟岩土或其他材料的力学行为。FLAC3D中的所有本构模型都是以动态链接库文件（dll文件）的形式提供给用户的。在计算过程中主程序会自动调用用户所指定的本构模型的dll文件。该软件允许用户在C++的环境下将自定义的本构模型编译成动态链接库 dll文件，由主程序调用执行，实现本构模型的二次开发。

考虑能量耗散后的模型本构方程是在加载阶段本构方程不变，与FLAC3D中自带的本构模型一致，只是卸载阶段有所变化。卸载阶段的本构方程如下：

3.2 二次开发程序计算流程

FLAC3D的二次开发程序主要分为以下几步进行计算求解：①FLAC3D自带模型的模拟试算，在该计算过程中利用 Fish语言记录各单元弹性能变化，即SAOB和SAOC；②利用Lajtai应力比率法判断围岩破损区；③对于破损区的单元，利用式（2）、（3）求出系数t；④将参数 t代入FLAC3D二次开发程序中，进行考虑能量耗散作用的下的弹脆性模型开挖计算。

图4 程序计算流程图Fig.4 The process of program

4 数值模拟计算

4.1 工程背景

锦屏一级水电站[20]位于四川省凉山彝族自治州木里藏族自治县和盐源县境内，是雅砻江干流下游龙头梯级电站，总装机容量360×104kW，年均发电166.2亿千瓦时。电站主要由拦河坝、右岸泄洪洞、右岸引水发电系统及开关站等组成，拦河坝为混凝土双曲拱坝，最大坝高为305 m，为世界第一高坝。电站地下厂房洞室群规模巨大，主要由引水洞、主厂房、母线洞、主变室、尾水调压室和尾水洞等组成。厂房轴线方向为NW65°，厂内安装6台600 MW机组，厂房全长为276.99 m，吊车梁以下开挖跨度为25.60 m，以上开挖跨度为28.90 m，开挖高度为 69.30 m；主变室长为 201.60 m，宽为19.30 m，总高为32.70 m。尾水调压室采用“三机一室一洞”布置型式，设置两个圆型调压室，直径分别为34.00、38.00 m（下室），高为80.50、81.50 m。

4.2 计算模型及相关参数

以锦屏水电站为工程背景，运用ANSYS有限元分析软件进行建模和单元剖分。模型尺寸大小为：x水平方向为540 m长，厂房区域分布在中间位置；y垂直方向为主厂房底板以下200 m至地表；z厚度方向为3个机组厚度。计算模型如图5所示。地应力数值根据已有的实测地应力点数值进行选取。岩体力学参数见表 2。裂纹的相关参数为：裂纹与最大主应力的夹角为 45°，摩擦系数为 0.50，裂纹半长为0.002 5 m，裂纹密度为0.024。

图5 数值模拟模型Fig.5 Model of simulation

表2 计算采用的岩体力学参数Table 2 Physical parameters of the rock

4.3 数值计算结果

为充分说明各机组模型的变形及应力特征，分别在3个机组的主厂房、主变室以及尾水调压井中选取了具有代表性的关键点，提取了各关键点部位的位移值。关键点的具体布置如图6所示。

图6 关键点分布图Fig.6 Distribution map of key points

为验证计算结果的正确性，选取了具有代表性的部分监测点的位移现场实测数据与数值计算结果进行对比，分别见图7、8。

分析图7、8可知，现场实测、考虑能量耗散与不考虑能量耗散的三者所得的位移变化规律基本一致：当洞室开挖高程至监测点附近时，测点位移迅速增大，其增长速率随着开挖高程的降低而逐渐平缓，最终经过一段时间的发展之后而逐渐趋于稳定。但在考虑能量耗散作用的情况下，开挖完成后洞室围岩的位移值要比不考虑能量耗散时的大，平均增幅约10%～20%。考虑能量耗散作用的计算结果与现场实测结果吻合得更好。

图7 主厂房下游边墙1 666 m高程处测点位移变化曲线Fig.7 Deformation of the downstream side wall of main power house at elevation 1 666 m

图8 主变室上游边墙1 668 m高程处测点位移变化曲线Fig.8 Deformation of the upstream side wall of main transformer house at elevation 1 668 m

开挖完成后洞室围岩的位移分布趋势是相似的。都表现为：洞室群开挖后，围岩总体朝向开挖临空面变形。边墙部位变形较大，以水平变形为主，拱顶下沉，底板回弹。且主厂房上下游边墙的中部关键点处的位移值较大。

在主厂房与母线洞的连接部位，即5#关键点处产生了较大位移。在1#机组中为11.32 cm的水平位移（见图9(a)），2#机组中为12.09 cm的水平位移（见图9(b)），3#机组中为11.43 cm的水平位移（见图9(c)）。此外，穿过厂房区域的断层f14对洞群稳定性的影响较大。由于，该断层穿过2#、3#机组的主变室部位，造成2#机组主变室下游边墙关键点产生13.61 cm的水平位移值（见图10），3#机组顶拱关键点产生16.02 cm的下沉位移（见图11）。因此，在实际工程中应在断层经过部位加强支护，防止因局部破坏影响整体稳定性。

图9 主厂房下游边墙5#关键点位移变化Fig.9 Deformation of the downstream side wall of main power house KP5

图10 2#机组主变室下游边墙10#关键点位移变化Fig.10 Deformation of the downstream side wall of main transformer house KP10 in unit #2

图11 3#机组主变室顶拱8#关键点位移变化Fig.11 Deformation of the crown of main transformer house KP8 in unit #3

图12为开挖完成后的塑性区三维图。从图可知，开挖完成后各机组剖面洞室周边的塑性区分布趋势一致，均为在洞室的边墙中部位置塑性区深度较大，且劈裂区的分布范围与塑性区相似。考虑能量耗散作用时的范围较不考虑能量耗散时的要大。边墙两侧约1倍洞径范围内处于劈裂破损区，且边墙表面处于受拉状态，易产生劈裂及岩体剥落现象。

图12 开挖完成后的塑性区三维图Fig.12 The 3D plastic zone after excavation

5 结 论

（1）从能量耗散的角度出发，考虑裂纹扩展开裂所产生的能量耗散对岩体本构方程的影响，对FLAC3D原有本构方程进行优化，得出考虑能量耗散作用下的高地应力脆性岩体计算模型。

（2）通过数值模拟计算，对锦屏一级水电站的开挖稳定性进行了分析，分析结果表明，考虑变弹模的弹脆性模型能更真实地反映实际情况；同时指出了工程开挖过程中应重点关注的部位，具有一定的工程指导意义。

（3）为以后进一步研究张开位移提出了一个可行的方法，同时，对FLAC3D的二次开发有一定的借鉴作用。

注：本文相关工作主要是在山东大学期间完成，在中国建筑科学研究院进行了后续的完善和补充工作。

[1]陈梦德. 高地应力地区岩体变形性质探讨[J]. 水电站设计, 1987, (2): 100.

[2]史红光. 二滩水电站地下厂房围岩稳定性因素评价[J].水电站设计, 1996, 15(2): 75－78.

[3]宋宏伟, 贾颖绚, 段艳燕. 开挖中的围岩破裂性质与支护对象研究[J]. 中国矿业大学学报, 2006, (2): 192－193.SONG Hong-wei, JIA Ying-xuan,DUAN Yan-yan. Study on characteristics of rock broken by excavation and roof supporting object[J]. Journal of China University of Mining & Technology, 2006, (2): 192－193.

[4]HIBINO S, MOTOJIMA M. Characteristic behavior of rock mass during excavation of large caverns[C]//Proceedings of the 8th ISRM Congress. Tokyo: A. A.Balkema, 1995: 583－586.

[5]彭琦, 王俤剀, 邓建辉, 等. 地下厂房围岩变形特征分析[J]. 岩石力学与工程学报, 2007, 26(12): 2586－2587.PENG Qi, WANG Di-kai, DENG Jian-hui, et al. Analysis of surrounding rock deformation characteristics in underground powerhouse[J]. Chinese Journal of Rock Mechanics and Engineering, 2007, 26(12): 2586－2587.

[6]江权, 冯夏庭, 陈国庆, 等. 高地应力条件下大型地下洞室群稳定性综合研究[J]. 岩石力学与工程学报, 2008,27(增刊 2): 3768－3777.JIANG Quan, FENG Xia-ting, CHEN Guo-qing, et al.Stability study of large underground caverns under high geostress[J]. Chinese Journal of Rock Mechanics and Engineering, , 2008, 27(Supp.2): 3768－3777.

[7]JING Lan-ru. Formulation of discontinuous deformation analysis(DDA)—An implicit discrere element model for block systems[J]. Engineering Geology, 1998, 49(3－4):371－381.

[8]JIANG Yu-jing, LI Bo, YAMASHITA Yuji. Simulation of cracking near a large underground cavern in a discontinuous rock mass using the expanded distinct element method[J]. International Journal of Rock Mechanics & Mining Sciences, 2009, 46(1): 97－106.

[9]GOLSHANI A, ODA M, OKUI Y, et al. Numerical simulation of the excavation damaged zone around an opening in brittle rock[J]. International Journal of Rock Mechanics & Mining Sciences, 2007, 44(6): 835－845.

[10]EXADAKTYLOST G E, TSOUTRELISJ C E. Pillar failure by axial splitting in brittle rocks[J]. International Journal of Rock Mechanics and Mining Sciences &Geomechanics Abstracts, 1995, 32(6): 551－562.

[11]王川婴, 葛修润, 白世伟. 数字式全景钻孔摄像系统及应用[J]. 岩土力学, 2001, 22(4): 522－525.WANG Chuan-ying, GE Xiu-run, BAI Shi-wei. The digital panoramic borehole camera system and its application[J]. Rock and Soil Mechanics, 2001, 22(4):522－525.

[12]李守巨, 吴金花, 刘迎曦, 等. 某大跨度地下厂房开挖过程变形监测研究[J]. 工程建设, 2008, 40(6): 1－5.LI Shou-ju, WU Jin-hua, LIU Ying-xi, et al. Study on deformation monitoring of large span underground powerhouse during excavation process[J]. Engineering Construction, 2008, 40(6): 1－5.

[13]范天佑. 断裂理论基础[M]. 北京:科学出版社, 2003.

[14]林卓英, 吴玉山, 关玲琍. 岩石在三轴压缩下脆-延性转化的研究[J]. 岩土力学, 1992, 13(2): 45－53.LIN Zhuo-ying, WU Yu-shan, GUAN Ling-li. Study of the inverting brittle and ductile of rock under triaxial compression [J]. Rock and Soil Mechanics, 1992, 13(2):45－53.

[15]陈景涛, 冯夏庭. 高地应力下硬岩的本构模型研究[J].岩土力学, 2007, 28(11): 2271－2278.CHEN Jing-tao, FENG Xia-ting. Study of constitutive model for hard rock under high geostresses[J]. Rock and Soil Mechnics, 2007, 28(11): 2271－2278.

[16]吕森鹏, 陈卫忠, 贾善坡, 等. 脆性岩石破坏试验研究[J].岩石力学与工程学报, 2009, 28(增刊1): 2772－2777.LÜ Sen-peng, CHEN Wei-zhong, JIA Shan-po, et al.Experimental study of brittle rock failure[J]. Chinese Journal of Rock Mechanics and Engineering, 2009,28(Supp.1): 2772－2777.

[17]潘鹏志, 冯夏庭, 周辉. 脆性岩石破裂演化过程的三维细胞自动机模拟[J]. 岩土力学, 2009, 30(5): 1471－1476.PAN Peng-zhi, FENG Xia-ting, ZHOU Hui. Failure evolution processes of brittle rocks using 3D cellular automaton method[J]. Rock and Soil Mechanics, 2009,30(5): 1471－1476.

[18]刘豆豆, 陈卫忠, 杨建平, 等. 脆性岩石卸围压特性试验研究[J]. 岩土力学, 2009, 30(9): 2588－2594.LIU Dou-dou, CHEN Wei-zhong, YANG Jian-ping, et al.Experimental research on strength characteristic of brittle rock unloading confining pressure[J]. Rock and Soil Mechanics, 2009, 30(9): 2588－2594.

[19]刘波, 韩彦辉. FLAC原理、实例与应用指南[M]. 北京:人民交通出版社, 2005.

[20]中国水电顾问集团成都勘测设计研究院.锦屏一级水电站可行性研究报告[R]. 成都: 中国水电顾问集团成都勘测设计研究院, 2008.