考虑初始渗透压力的裂隙岩体本构模型二次开发及其验证

张卫东， 陈士海

(华侨大学 土木工程学院， 福建 厦门 361021)

考虑初始渗透压力的裂隙岩体本构模型二次开发及其验证

张卫东， 陈士海

(华侨大学 土木工程学院， 福建 厦门 361021)

为充分研究含水裂隙岩体的力学特性，在摩尔库伦模型的基础上，充分考虑初始渗透压力的作用.对本构方程中的柔度张量进行适当的修正，生成考虑初始渗透压力的裂隙岩体本构模型.为了增加对比性，按照同样的理论推导方法及二次开发流程，在不考虑初始渗透压力的影响下，成功地生成未考虑初始渗透压力的裂隙岩体本构模型.在此基础上，利用C++语言，将以上两种模型编译成可运行的动态链接库文件.为验证自定义本构模型的准确性，在FLAC3D中建立同一个50 m×50 m×1 m的圆形隧道数值模型进行对比分析.通过分别求解计算3种本构模型,分析各圆形隧道上4个对称点的位移变化情况.结果表明:将初始渗透压力直接引进本构方程中，利用C++二次开发出来的本构模型是成功的.

裂隙岩体; 初始渗透压力； 柔度张量； 本构关系； 二次开发

自然界中的岩体是一种含有初始损伤的介质，它包含着各式各样的非连续面，如断层、节理、裂隙，而含有初始地下水的裂隙岩体在工程的施工中会经常碰到.1999年，朱维申等[1]在Betti能量互易定理和不可逆热力学定律的基础之上，建立了裂隙岩体的三维脆弹性断裂损伤本构模型.2008年，蓝航等[2]在几何损伤力学理论的基础上，建立了节理岩体的采动损伤本构模型.2012年，齐银萍[3]在理论分析本构模型的基础之上，在FLAC3D中编写出了损伤流变本构方程.2014年，文献[4-5]等利用FLAC3D建立了含水裂隙岩体的数值模型.但是，以上研究开发出来的模型大部分是针对无水裂隙岩体，而对于有水裂隙岩体，则是利用模拟软件中相应的渗流模块进行流固耦合计算，而在本构模型当中没有将渗透压力带入到本构方程中进行二次开发.因此，初始含水裂隙岩体工程的模拟分析存在一定的局限性.本文在摩尔库伦模型的基础之上，充分考虑初始渗透压力的作用，并将其直接引进本构方程.

1 含水裂隙岩体的本构关系

从宏观力学效果方面考虑，裂隙岩体的初始损伤及其损伤的演化可用柔度张量的改变进行表示，而裂隙的半径、法向和切向刚度系数、传压和传剪系数、初始渗透压力等都会在不同程度上对柔度张量产生影响.文献[6]通过坐标运算和叠加原理得到了含水裂隙岩体柔度张量的表达式，即

1.1 含水裂隙岩体的初始损伤柔度张量

(2)

式(2)中：KⅠ,KⅡ,KⅢ分别为裂隙尖端的Ⅰ，Ⅱ，Ⅲ型应力强度因子；Ω为裂隙的外表面，积分沿整个裂隙表面进行.

现假定裂隙都为圆形，且半径为a，在查阅应力强度因子手册后，可获得裂隙尖端的Ⅰ，Ⅱ，Ⅲ型应力强度因子[8].

(3)

式(3)中：σ，τ分别为应力张量在裂隙面法向和切向上的投影分量；G1，G2为裂隙形状系数.

在某一应力状态下，裂隙面上实际的法向和切向应力张量是需要进行修正的.为此,引进Cn和Cv,反映裂隙在某一应力状态下，裂隙面传递法向应力和切向应力的能力[9].对于含水裂隙岩体，还需要考虑初始渗透压力(p)的作用，因此,综合以上情况,作用在裂隙面上的有效法向和切向应力分别为

τ′=(1-Cv)τ.

(5)

将式(4)，(5)带入式(3)，化简可得到单位体积岩体内该裂隙产生的附加弹性应变能φd，即

[7]中的方法，求得单位体积岩体内裂隙产生的附加弹性应变能φd.

且有

其中：ni，nj，nk，nl均为裂隙面的法向向量分量.

通过比较式(7)及文献[7]中的表达式，可得由于初始渗透压力存在而产生的附加柔度张量，即

(8)

1.2 裂隙扩展产生的损伤演化柔度张量

1.2.1压剪应力状态下的损伤演化柔度张量 在压剪应力状态下，岩体中的裂隙随着外加作用载荷的增加，开始闭合摩擦滑动，压剪起裂，形成分支张型裂隙，新裂面顺着最大主应力方向不断延伸发展，直至位于分支张型裂隙尖端的微裂隙损伤区相互汇合，导致宏观裂隙的击穿贯通，而使岩体局部破坏[10].

大量的现场勘测结果和理论计算表明，岩体中的裂隙在压剪应力的作用下，将在最大和最小主压应力组成的平面内沿着最大主压应力σ1的方向稳定扩展[11]，其翼形分支裂隙发生扩展.

设远场Cauchy应力张量在裂隙面法向和切向的投影分别为正应力σn和剪应力τn，即

式(9)中：θ为裂隙与最大主压应力的夹角.

当考虑初始渗透压力的影响时，裂隙面上实际传递的法向应力σne和切向应力τne分别为

压剪应力状态下，考虑初始渗透压力的裂隙岩体其裂隙扩展产生的损伤演化柔度张量可根据应变能等效原理和自洽理论推导得出[12]，压剪应力作用下分支裂隙尖端瞬时应力强度因子[13]为

(11)

式(11)中：τe=τne-τs，τs为裂隙面上下表面抵抗外力的能力.

τe=fsσn+Cs表示摩尔库伦准则[14]，fs为裂隙面的摩擦系数(fs=tanφ，φ为岩体介质的内摩擦角)，Cs为裂隙面的黏结力，σn为裂隙面的法向应力，且有σn=σne.

扩展中的翼形分支裂隙逐渐沿平行于最大主压应力的方向稳定扩展，当分支裂隙扩展至KI=KlC时，裂隙停止扩展[15].其中，KI为翼形分支裂隙尖端的应力强度因子；KIC为裂隙岩体的断裂韧度.翼形分支裂隙尖端应力强度因子计算方法采用Kemeny计算模型[13].

1.2.2 拉剪应力状态下的损伤演化柔度张量 在拉剪应力状态下，岩体中的裂隙受拉后,出现张开现象，其裂隙表面的黏结力将变为零,而不能传递拉应力.在最大周向应力准则的基础之上，当分支裂隙尖端瞬时应力强度因子KI0达到岩体的断裂韧度KlC时，裂隙便开始扩展.

拉剪应力状态下，分支裂隙尖端瞬时应力强度因子为

(12)

式(12)中：开裂角θ0可参考文献[12]中的表达式.

在文献[16]的基础之上，翼形分支裂隙的扩展长度，即

可得拉剪应力状态下考虑初始渗透压力的裂隙岩体其裂隙扩展产生的损伤演化柔度张量[13].

2 裂隙岩体自定义本构的二次开发

利用FLAC3D的二次开发平台，可以按照自己的意愿进行本构模型的自定义.自定义本构模型的计算流程图，如图1所示.

图1 自定义本构模型的计算流程图Fig.1 Calculation flow chart of the constitutive model

3 模型的验证及结论

在FLAC3D中,建立的统一圆形隧道模型在水平方向上和竖直方向上的长度都为50 m，进深为1 m，隧道半径为5 m.模型划分单元数为2 500，节点数为5 200.模型中：所有节点y方向上的位移被约束，z=-25 m的底面及x=±25 m的左右两面固定，即模型底面及左右两面上的节点约束其位移.

模型中，节点施加的初始应力条件为

(14)

式(14)中：σx，x,σy，y,σz，z分别表示模型节点受到的x方向、y方向和z方向上的应力.

经测定，在模型底面渗透压力为0.275 MPa，模型顶面为0.225 MPa，且沿z轴呈线性分布；隧道围岩弹性模量为15 MPa，容重为25 kN·m-3，内摩擦角为34°，粘聚力为0.85 MPa，抗拉强度为2 MPa，断裂韧度为1.2 MPa·m-1/2；裂隙特征长度为0.1 m，内摩擦角为11.5°，体积密度为0.33，裂隙面轴向夹角为60°，法向刚度为2.5 MPa·m-1，切向刚度为1.15 MPa.

图2 数值计算模型Fig.2 Numerical calculation model

为验证自定义本构模型的准确性，采用“建立同一数值模型、分别调用3种不同本构模型(即摩尔库伦模型、未考虑初始渗透压力p的裂隙岩体本构模型，以及考虑初始渗透压力p的裂隙岩体本构模型)”的方法进行数值模拟，并分别将以上3个数值模拟记为SZ1，SZ2和SZ3.数值计算模型,如图2所示.

图3 圆形隧道上的对称节点布置Fig.3 Layout of symmetric nodes in circular tunnel

为保证数值模拟中模型所受到的渗透压力条件的统一性，SZ1和SZ2直接通过“INITIAL pp 2.5e5 grad-0.1e4”设定初始的孔隙水压力；SZ3则采用FISH函数及“PROPERTY XXX”命令定义初始的渗透压力.计算求解完成后，统一选取各圆形隧道模型上z方向上位移变化相同的对称点(0,0,5)，(0,0,-5),以及x方向上的位移变化点(5,0,0),(-5,0,0)，如图3所示.

SZ1，SZ2，SZ3求解后，各圆形隧道对称位置上节点的位移(s)变化曲线，如图4所示.图4中:n为时步.由图4可知：各模型节点的位移最终趋于某一个固定值.

SZ3中模型的最大不平衡力记录图，如图5所示.由图5可知：最大不平衡力随着求解时步的增加最终接近于零.因此,SZ1,SZ2,SZ3求解计算到8 000个时步时,系统已经达到了平衡状态.

(a) SZ1节点

(b) SZ2节点 (c) SZ3节点图4 各节点位移变化Fig.4 Nodal displacement

图5 SZ3最大不平衡力记录Fig.5 Maximal unbalanced force of SZ3

由图4(a)可知：各圆形隧道对称位置上的节点位移变化趋势是大致相同的，即节点1 301，3 901，1 353，1；在0～500个时步内，4个对称节点的位移受到外荷载的作用急剧地增加；在500～8 000个时步内，节点3 901的位移在8～9 mm起伏，基本保持不变；节点1 353的位移在0.5 mm左右保持不变；节点1的位移在0.6 mm左右保持不变；在500～4 000个时步内，节点1 301的位移缓慢地增加；在4 000～8 000个时步内，节点1 301的位移在4～5 mm之间呈稳定状态.这说明通过控制数值模拟条件的统一性验证自定义本构模型的方法是成功的.

图4(c)中节点3 901的位移最终趋于稳定值9 mm左右，而图4(a)中节点3 901的位移最终趋于稳定值8 mm左右，这说明采用摩尔库伦模型进行数值模拟计算求解后得到的节点位移偏于保守.

图4(b)中各对称节点最终平衡状态的位移均比图6大，以节点1 353，1为例，图4(b)中位移分别为0.6,0.6 mm，而图4(c)中位移分别为0.4,0.4 mm，这说明当裂隙岩体中存在初始状态下的地下水时，直接通过命令“INITIAL pp”赋予模型初始的孔隙水压力，其求解计算得到的对称节点上的位移比直接通过修正本构模型中刚度矩阵得到对称节点上的位移更大.

4 结论

1) 当裂隙岩体中存在初始地下水时，考虑初始渗透压力的影响是非常有必要的.

2) 在充分考虑初始渗透压力的作用，并将其直接引进本构方程中，利用C++二次开发出来的本构模型是成功的.通过修正摩尔库伦本构模型中的刚度矩阵考虑初始渗透压力的影响会更结合实际，其数值模拟结果相对于直接设置初始孔隙水压力时将更加精确.

3) 当初始条件统一时，采用自定义本构模型数值模拟出来的结果与采用摩尔库伦模型数值模拟出来的结果大同小异，模型中各对称节点的位移变化趋势大致相同，但采用摩尔库伦模型计算得到的结果相对于采用自定义本构模型计算得到的结果偏于保守.

因此，综合以上结论分析，此次二次开发出来的本构模型是成功的.同时，该模型也给一些包含地下水的岩体工程的数值模拟工作提供了一定的参考价值.

参考文献：

[1] 朱维申,张强勇.节理岩体脆弹性断裂损伤模型及其工程应用[J].岩石力学与工程学报,1999,18(3):245-249.

[2] 蓝航,姚建国,张华兴,等.基于FLAC3D的节理岩体采动损伤本构模型的开发及应用[J].岩石力学与工程学报,2008,27(3):572-579.

[3] 齐银萍.裂隙岩体三维损伤流变模型研究及工程应用[D].济南:山东大学,2012：11-40.

[4] 王昆.含裂隙水巷道变形破坏特征研究[D].淮南:安徽理工大学,2014:9-33.

[5] 王昆,赵光明,孟祥瑞，等.含水裂隙岩体本构模型及数值模拟理论研究[J].地下空间与工程学报,2015,11(4):901-908.

[6] 平杨,郑少河,白世伟.考虑渗透压力的裂隙岩体断裂损伤本构模型研究[C]∥中国岩石力学与工程学会第六次学术会议.北京:中国科学技术出版社,2000:134-136.

[7] 郑少河.裂隙岩体渗流场： 损伤场祸合理论研究及应用[D].武汉：中国科学院武汉岩土力学研究所,2000:77-118.

[8] 朱维申,李术才,陈卫忠.节理岩体破坏机理和锚固效应及工程应用[M].北京：科学出版社,2002:53-72.

[9] KEMENY J,COOK N G W.Effective moduli, non-linear deformation and strength of a cracked elastic solid[J].J Rock Mech Min Sci and Geomech Abstr,1986,23(2):107-118.

[10] 易顺民,朱珍德.裂隙岩体损伤力学导论[M].北京:科学出版社,2005:24-36.

[11] 赵延林,曹平,林杭,等.渗透压作用下压剪岩石裂纹流变断裂贯通机制及破坏准则探讨[J].岩土工程学报,2008,30(4):511-517.

[12] 张电吉.裂隙岩质边坡渗流场与应力场耦合分析及工程应用[D]. 武汉:中国科学院武汉岩土力学研究所, 2003:11-21.

[13] 柴红保,曹平,赵延林,等.裂隙岩体损伤演化本构模型的实现及应用[J].岩土工程学报,2010,32(7):1047-1053.

[14] 曹平,蒲成志.单压下有序多裂隙脆性材料破坏机制及其简化模型[J].中国有色金属学报,2011,21(10):2659-2668.

[15] 王涛,韩煊,赵先宇,等.FLAC3D数值模拟方法及工程应用： 深入剖析FLAC3D 5.0[M].北京:中国建筑工业出版社,2015:168-199.

[16] 郑颖人,孔亮.岩土塑性力学[M].北京:中国建筑工业出版社,2010:155-198.

(责任编辑： 陈志贤 英文审校： 方德平)

Secondary Development on Fractured Rock Constitutive Model Considering Initial Seepage Pressure

ZHANG Weidong， CHEN Shihai

(College of Civil Engineering, Huaqiao University, Xiamen 361021, China)

In order to study of mechanical properties of water-bearing fractured rock, based on Mohr coulomb model, considering the role of initial seepage pressure, flexibility tensor of the constitutive equation is appropriately amended, the constitutive model of jointed rock considering initial seepage pressure is obtained. By comparison, the constitutive model of jointed rock is also obtained without considering initial seepage pressure, according to the same theoretical derivation method and secondary development process. The C++program of two models is compiled to verify the accuracy of the models, through building the same 50 m×50 m×1 m numerical model for circular tunnel in the FLAC3D, three kinds of constitutive model and 4 symmetric point displacements are analyzed in circular tunnel. The result shows that constitutive model introducing initial seepage pressure is successful accurate. Keywords：fractured rock； initial seepage pressure; flexibility tensor; constitutive relation; secondary development

10.11830/ISSN.1000-5013.201703007

2016-08-02

陈士海(1964-)，男，教授，博士，主要从事岩土与地下工程防灾减灾，以及工程结构抗震与抗爆的研究.E-mail:cshblast@163.com.

高等学校博士学科点专项基金资助项目(20113718110002); 爆炸冲击防灾减灾国家重点实验室基金资助项目(DPMEIKF201307); 华侨大学科研基金资助项目(13BS402); 华侨大学研究生科研创新能力培育计划资助项目(1400204010)

TU 91

A

1000-5013(2017)03-0319-06