考虑节理面法向蠕变的节理岩体蠕变模型
熊良宵，杨林德
(同济大学 地下建筑与工程系 岩土及地下工程教育部重点实验室，上海，200092)
摘 要：采用Fish语言对FLAC3D软件中的Interface单元进行修正，将Interface单元的法向刚度和切向刚度转化为时间的蠕变函数，便可得到在数值分析中考虑节理面法向蠕变的节理岩体蠕变模型。计算结果表明：节理岩体的单轴蠕变量随着节理面间距的增大而减小；当节理面的倾角小于45?时，节理岩体单轴蠕变量随着节理面倾角的增大而增大，但当节理面倾角大于45?时，轴向蠕变量反而随着节理面倾角的增大而减小；若固定节理面的切向蠕变参数，增加其法向蠕变参数，则节理岩体的轴向蠕变量也会随之减小。
关键词：蠕变模型；节理岩体；法向刚度；切向刚度；FLAC3D软件
中图分类号：TU45 文献标识码：A 文章编号：1672-7207(2009)03-0814-08
Creep model for rock mass considering normal creep of rock joint plane
XIONG Liang-xiao, YANG Lin-de
(Key Laboratory of Geotechnical and Underground Engineering of Ministry of Education,
Department of Geotechnical Engineering, Tongji University, Shanghai 200092, China)
Abstract: Revised by Fish language, the normal stiffness and shear stiffness of interface element in FLAC3D were translated into creep functions of time, and the creep model for rock mass was realized in numerical analysis when considering normal creep of rock joint plane. The calculation results show that the uniaxial creep deformation for rock mass decreases with the increase of joint interval, and it increases with the increase of dip angle of joint plane if the angle is smaller than 45?, while the result will be quite on the contrary if the angle is larger than 45?. If the shear creep parameters of rock joint are constant, the uniaxial creep deformation will decrease with the increase of normal creep parameters.
Key words: creep model; rock mass; normal stiffness; shear stiffness; FLAC3D software
节理岩体的蠕变特性一直是岩石力学的重要研究内容，它主要由完整岩石的蠕变和节理面的蠕变组 成。目前，人们对完整岩石和节理面蠕变模型的研究较多[1-8]，而对节理岩体蠕变模型的研究较少。
Curran等[9]采用位移不连续法建立了节理岩体的蠕变模型；刘考学[10]利用对互层岩体或节理岩体的不同介质进行“分离、集中”得到了互层节理岩体的二维和三维黏弹性模型等；杨海天等[11]将受一组正交节理切割的岩体结构视为广义复合材料，用“代表元”的方法，建立了复合蠕变本构方程；杨松林等基于文献[12]中节理岩体的弹性本构关系，在已知岩石和节理蠕变规律的前提下，假定岩石体积变形和节理法向压缩变形为弹性变形，忽略节理的剪胀现象，认为只有岩石畸变和节理剪切滑移与时间有关，推导了含3组相交节理的岩体蠕变模型[13]。
对于深埋隧洞，隧洞开挖造成深部岩体受到加载和卸载的交替作用，节理面的法向压缩变形更不是固定的弹性变形，而是随着时间发生变化。因此，若将节理面的法向压缩变形视为固定的弹性变形，并以此推导节理岩体的蠕变模型，分析节理岩体的蠕变规律，则计算得到的蠕变变形将与实际变形存在偏差。在此，本文作者在文献[12-15]的基础上，通过将FLAC3D软件中Interface单元的法向刚度和切向刚度转化为时间的蠕变函数，以便能得到考虑节理面法向蠕变的节理岩体蠕变模型，并对含1组节理岩体进行蠕变计算，通过解析解和数值解共同分析节理面间距、倾角和蠕变参数对计算结果的影响。
1 节理岩体的蠕变模型
1.1 岩石的蠕变关系
假定岩石材料是各向同性的，岩石体变是弹性的，畸变是黏弹性的，则有[13]：
，
。 (1)
偏应变和偏应力的定义为：
由式(1)和(2)可得到蠕变情况下岩石的应力与应变关系：
岩石的蠕变柔度张量可写为：
由于
，式(4)可简化为：
1.2 节理面的蠕变关系
Huang等[12]给出了单一节理面的弹性本构关系：
假定节理面的变形特性是各向同性的，则有Dss=Dtt，Dsn=Dtn，Dnt=Dns，Dst=Dts，且有[12]：
Dss=1/Knn，Dss=1/Kss。 (8)
式中：Knn为节理面的法向刚度；Kss为节理面的切向刚度。
对于粗糙节理面，发生剪切变形时会引起法向膨胀，即Dns=λDss，Dnt=λDtt。其中，λ为膨胀系数，可表示为：
杨松林等[13]在式(7)的基础上，将变形增量、柔度矩阵和应力增量变为时间的函数，从而给出了单一节理面的蠕变关系：
而杨松林等[13]为了简化，忽略节理面在法向应力作用下的蠕变变形，将法向变形视为固定的弹性变形。实际上，当节理面较软弱，而法向应力较大时，节理面的法向蠕变变形会比较明显。因此，可将法向刚度和切向刚度也表示为时间的函数[14-15]，并代入式(8)可得到：
若节理面的法向和切向蠕变变形规律均符合广义Kelvin模型，则有：
若不考虑节理面的剪胀效应，且认为节理面的法向变形不引起切向变形，节理面内2个相互垂直的剪切变形也互不影响，则有：
基于文献[12]中含M组节理面的节理岩体弹性柔度张量，可得到相应节理面的蠕变柔度张量[13]：
1.3 节理岩体的蠕变关系
节理岩体的蠕变
来自2部分：一是岩石的蠕变
，二是节理面的蠕变
。根据前面的分析，节理岩体的蠕变关系可表示为[13]：
将式(5)和(16)代入式(17)和(18)，可得到考虑节理面法向蠕变的节理岩体蠕变关系。
2 解析解分析
将文献[13]中的3组节理岩体模型简化成含1组节理的岩体模型(见图1)，并分析节理面法向蠕变对节理岩体蠕变变形的影响，以及节理面间距、倾角和蠕变参数对蠕变变形的影响。
图1 岩体的一组节理分布
Fig.1 Rock mass with one set of joints
设单一组节理面的倾角为θ，节理间距为d。在文献[13]的基础上，假定完整岩石和节理面的蠕变规律均符合广义Kelvin模型，蠕变参数具体见表1~3，其中：G1，G2和η分别表示完整岩石的瞬时弹性模量、黏弹性模量和黏滞系数。岩石体积应变相对应的参数J2=0.25×10-3 MPa-1。
表1 完整岩石黏弹性参数
Table 1 Values of viscoelastic model for intact rock
表2 节理面法向黏弹性参数
Table 2 Values of normal viscoelastic model for rock joint
表3 节理面切向黏弹性参数
Table 3 Values of shear viscoelastic model for rock joint
假设不考虑节理面的剪胀效应，且满足式(14a)，(14b)和(14c)，经计算可分别得到完整岩石蠕变、节理面法向蠕变和节理面剪切蠕变引起的岩体轴向蠕变量，从而得到节理岩体总的轴向蠕变量。因此，可设节理面倾角为30?，节理面的间距为2 m，沿Z轴施加的单轴压应力为10 MPa，此时，节理岩体的总轴向蠕变量见图2(a)，节理面切向蠕变和法向蠕变等引起的岩体轴向蠕变量与总轴向蠕变量之比见图2(b)。
(a) 节理岩体的总轴向蠕变量；(b) 蠕变量与节理岩体总轴向蠕变量之比
图2 节理岩体的轴向蠕变量
Fig.2 Uniaxial creep strain of rock mass
由图2可见，节理面法向蠕变引起的岩体轴向蠕变量占岩体总轴向蠕变量的15.5%。当节理面的倾角小于30?时，节理面法向蠕变引起的岩体轴向蠕变量所占的比例将更大，因此，不可忽视节理面的法向蠕变变形。
为了分析节理间距对轴向蠕变变形的影响，对节理面倾角为0?，节理面间距d分别为2，3和4 m时的轴向蠕变进行对比，结果见图3。
间距/m: 1—2; 2—3; 3—4
图3 节理面间距d对轴向蠕变量解析解的影响
Fig.3 Influence of joint plane interval on analytical results of axial creep
由图3可见，随着节理面间距的增大，节理岩体的轴向蠕变量会随之减小。
为了分析节理面倾角对轴向蠕变量的影响，对节理面间距为2 m，节理面倾角为0?，15?，30?，45?，60?，75?和90?时的轴向蠕变进行对比，结果见图4。
由图4可见，当节理面倾角由0?逐渐增加到45?时，节理岩体的轴向蠕变量会随之增大。而当节理面倾角由45?增至90?时，轴向蠕变量反而减小。因此，节理岩体轴向蠕变量并不是随着节理面倾角的增加而一直增加。
倾角/(?): 1—0; 2—15; 3—30; 4—45; 5—60; 6—75; 7—90
图4 节理面倾角对轴向蠕变量解析解的影响
Fig.4 Influence of dip angle of joint plane on analytical results of axial creep
为了分析法向蠕变与剪切蠕变参数对轴向蠕变量的影响，可固定剪切蠕变参数，将法向蠕变参数设为以下几组进行计算：① G1n=4G1s；G2n=4G2s，ηn=4ηs；② G1n=6G1s，G2n=6G2s，ηn=6ηs；③ G1n=8G1s，G2n=8G2s，ηn=8ηs。此时，节理面的倾角设为45?，节理面的间距为2 m，参数变化对节理岩体轴向蠕变量的影响见图5。
由图5可见，固定节理面的切向蠕变参数，增加其法向蠕变参数，则节理岩体的轴向蠕变量也会随之减小。
1—G1n=4G1s, G2n=4G2s, ηn=4ηs; 2—G1n=6G1s, G2n=6G2s, ηn=6ηs; 3—G1n=8G1s, G2n=8G2s, ηn=8ηs
图5 节理面蠕变参数对轴向蠕变量解析解的影响
Fig.5 Influence of creep parameter of rock joint on analytical results of axial creep
3 数值解分析
采用FLAC3D软件对不同类型节理岩体进行流变计算，蠕变本构模型采用Burgers模型，如图6所示。其中：EM和EK分别为Maxwell和Kewin黏弹性模量；ηM和ηK分别为Maxwell和Kewin黏滞系数。
图6 Burgers模型
Fig.6 Burgers model
计算时，不考虑Maxwell黏滞系数，使其退化为广义Kelvin模型。同时，采用Fish语言对Interface单元的法向刚度和切向刚度进行修改，使其与时间的关系符合广义Kelvin模型。在整个计算过程中，法向刚度和切向刚度均会随时间发生变化，其部分Fish语言如下：
def parameter_time
age=age+1e-3
creep_time=age
parameter_kn=1/((1/k_1n)*(1-exp(- (k_1n/v_1n)*
creep_time))+1/m_1n)
parameter_ks=1/((1/k_1s)*(1-exp(- (k_1s/v_1s)*
creep_time))+1/m_1s)
end
其中：m_ln为节理面的法向瞬时弹性模量；k_ln为节理面的法向黏弹性模量；v_ln为节理面的法向黏滞系数；m_1s为节理面的切向瞬时弹性模量；k_1s为节理面的切向黏弹性模量；v_1s为节理面的切向黏滞系数；creep_time为计算时间；parameter_kn和parameter_ks分别为节理面的法向和切向刚度。
对含1组节理面的节理岩体进行计算，长×宽×高为5 m×5 m×10 m，完整岩石的蠕变参数、节理面的法向和切向蠕变参数与表1~3中的参数一致。
为了分析节理面间距对节理岩体轴向蠕变量的影响，可假定节理面的倾角为0?，节理面间距分别为2，3和4 m，计算模型见图7。400 d后，节理面间距不同时节理岩体的轴向位移分布云图见图8；其整个蠕变时间过程中，轴向蠕变量的比较结果见图9。
间距/m: (a) 2; (b) 3; (c) 4
图7 不同节理面间距的计算模型
Fig.7 Computing models with different joint plane intervals
间距/m: (a) 2; (b) 3; (c) 4
图8 节理面间距对轴向位移的影响
Fig.8 Influence of joint plane interval on axial displacement
间距/m: (a) 1~2; (b) 2~3; (c) 3~4
图9 节理面间距对数值解的影响
Fig.9 Influence of joint plane interval on numerical results
由图8和图9可见，随着节理面间距的增大，节理岩体的轴向蠕变量会随之减小。
为了分析节理面倾角对轴向蠕变量的影响，可在节理岩体模型中只设置1条节理面，节理面倾角为30?，45?和60?时的节理岩体计算模型见图10，此时，节理面的间距大于5 m。400 d后，节理面倾角为30?，45?和60?时节理岩体的轴向位移分布云图，见图11；不同倾角对节理岩体轴向蠕变量的影响结果见图12。
倾角/(?): (a) 30; (b) 45; (c) 60
图10 不同节理面倾角的计算模型
Fig.10 Computing model with different dip angles of joint plane
倾角: (a) 30?; (b) 45?; (c) 60?
图11 节理面倾角对轴向位移的影响
Fig.11 Influence dip angle of joint plane on axial displacement
倾角/(?): 1—0; 2—15; 3—30; 4—45; 5—60; 6—75; 7—90
图12 节理面倾角对数值解的影响
Fig.12 Influence of dip angle of joint plane on numerical result
由图11和图12可知，当节理面倾角由0?逐渐增加到45?时，节理岩体的轴向蠕变量会随之增大；而当倾角大于45?时，随着节理面倾角的增加，节理岩体的轴向蠕变量呈减小趋势，这与解析解得到的结论是一致的。
为了分析法向蠕变与剪切蠕变参数对数值解的影响，法向蠕变参数与切向蠕变参数之间的关系与解析解一样，也分为3种情况。采用图10(b)所示的计算模型，400 d后，节理岩体的轴向位移分布云图见图13，轴向蠕变量的比较结果见图14。
(a) G1n=4G1s, G2n=4G2s, ηn=4ηs; (b) G1n=6G1s, G2n=6G2s, ηn=6ηs
图13 节理面蠕变参数对轴向位移的影响
Fig.13 Influence of creep parameter of joint plane on axial displacement
1—G1n=4G1s, G2n=4G2s, ηn=4ηs; 2—G1n=6G1s, G2n=6G2s, ηn=6ηs; 3—G1n=8G1s, G2n=8G2s, ηn=8ηs
图14 节理面倾角对轴向蠕变数值解的影响
Fig.14 Influence of creep parameter of joint plane on numerical results of axial creep
由图11(b)、图13以及由图14可知，当节理面法向蠕变参数分别为切向蠕变参数的4倍、6倍和8倍时，节理岩体的轴向位移和轴向蠕变量依次减小。由FLAC3D计算得到的结果也表明，若固定节理面的切向蠕变参数，则增加其法向蠕变参数，则节理岩体的轴向蠕变量也会随之减小。
综上所述，由FLAC3D软件对含一组节理面的节理岩体进行流变计算时，得到的结论与前面解析解得到的结论基本一致。
4 结 论
a. 若忽略节理面的法向蠕变量，所推导的节理岩体蠕变模型将不能反映节理岩体的真实蠕变规律。
b. 采用Fish语言对FLAC3D软件中的Interface单元进行修正，将其法向刚度和切向刚度转化为时间的蠕变函数，可以得到在数值分析中考虑节理面法向蠕变的节理岩体蠕变模型。计算结果表明：节理岩体的轴向蠕变量随着节理面间距的增大而减小；当节理面倾角小于45?时，节理岩体的轴向蠕变量随着节理面倾角的增大而增大，但当节理面倾角大于45?时，轴向蠕变量反而随着节理面倾角的增大而减小；固定节理面的切向蠕变参数，增加其法向蠕变参数，则节理岩体的轴向蠕变量会随之减小。
参考文献：
[1] 孙 钧. 岩石流变力学及其工程应用研究的若干进展[J]. 岩石力学与工程学报, 2007, 26(6): 1081-1106.
SUN Jun. Rock rheological mechanics and its advance in engineering applications[J]. Chinese Journal of Rock Mechanics and Engineering, 2007, 26(6): 1081-1106.
[2] 徐卫亚, 杨圣奇, 褚卫江. 岩石非线性黏弹塑性流变模型(河海模型)及其应用[J]. 岩石力学与工程学报, 2006, 25(3): 433-447.
XU Wei-ya, YANG Sheng-qi, CHU Wei-jiang. Nonlinear viscoelasto-plastic rheological model (Hohai model) of rock and its engineering application[J]. Chinese Journal of Rock Mechanics and Engineering, 2006, 25(3): 433-447.
[3] 邓荣贵, 周德培, 张倬元, 等. 一种新的岩石流变模型[J]. 岩石力学与工程学报, 2001, 20(6): 780-784.
DENG Rong-gui, ZHOU De-pei, ZHANG Zhuo-yuan, et al. A new rheological model for rocks[J]. Chinese Journal of Rock Mechanics and Engineering, 2001, 20(6): 780-784.
[4] 马明军, 钟时猷. 一个软弱岩石的黏弹塑性流变力学[J]. 中南工业大学学报: 自然科学版, 1990, 21(3): 236-241.
MA Ming-jun, ZHONG Shi-you. A new visco-elastic plastic rheological mechanical model for soft rocks[J]. Journal of Central South University of Technology: Natural Science, 1990, 21(3): 236-241.
[5] 徐卫亚, 杨圣奇. 节理岩石剪切流变特性试验与模型研究[J]. 岩石力学与工程学报, 2005, 24(Suppl.2): 5536-5542.
XU Wei-ya, YANG Sheng-qi. Experiment and modeling investigation on shear rheological property of joint rock[J]. Chinese Journal of Rock Mechanics and Engineering, 2005, 24(Suppl.2): 5536-5542.
[6] 沈明荣, 张学进. 规则齿形结构面剪切蠕变本构方程的参数分析[J]. 地下空间与工程学报, 2006, 2(6): 907-911.
SHEN Ming-rong, ZHANG Xue-jin. Parametric analysis of shear creep constitutive equation of discontinuous plane with regular dentation[J]. Chinese Journal of Underground Space and Engineering, 2006, 2(6): 907-911.
[7] 丁秀丽, 刘 建, 白世伟, 等. 岩体蠕变结构效应的数值模拟研究[J]. 岩石力学与工程学报, 2006, 25(Suppl.2): 3642-3649.
DING Xiu-li, LIU Jian, BAI Shi-wei, et al. Study on numerical simulation of structure effects of rock mass creep[J]. Chinese Journal of Rock Mechanics and Engineering, 2006, 25(Suppl.2): 3642-3649.
[8] 杨圣奇, 徐卫亚, 杨松林. 龙滩水电站泥板岩剪切流变力学特性研究[J]. 岩土力学, 2007, 28(5): 895-902.
YANG Sheng-qi, XU Wei-ya, YANG Song-lin. Investigation on shear rheological mechanical properties of shale in Longtan Hydropower Project[J]. Rock and Soil Mechanics, 2007, 28(5): 895-902.
[9] Curran J H, Crawford A M. A displacement discontinuity approach to modeling the creep behaviour of rock and its discontinuities[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 1983, 7(2): 245-268.
[10] 刘考学. 互层状岩体及节理岩体的黏弹性模型[J]. 陕西机械学院学报, 1990, 6(4): 292-307.
LIU Kao-xue. A viscoelastic model for interlayered rock or joint rock[J]. Journal of Xi’an University of Technology, 1990, 6(4): 292-307.
[11] 杨海天, 邬瑞锋. 双向正交节理岩体的一个复合蠕变模型[J]. 大连理工大学学报, 1991, 31(3): 267-273.
YANG Hai-tian, WU Rui-feng. Composite creep model for rock with bilateral orthogonal joint[J]. Journal of Dalian University of Technology, 1991, 31(3): 267-273.
[12] Huang T H. Elastic moduli for fractured rock mass[J]. Rock Mechanics & Rock Engineering,1995, 28(3): 135-144.
[13] 杨松林, 张建民, 黄启平. 节理岩体蠕变特性研究[J]. 岩土力学, 2004, 25(8): 1225-1228.
YANG Song-lin, ZHANG Jian-min, HUANG Qi-ping. Analysis of creep model of jointed rock[J]. Rock and Soil Mechanics, 2004, 25(8): 1225-1228.
[14] Curran J H, Crawford A M. A comparative study of creep in rock and its discontinuities[C]//21st US Symposium on Rock Mechanics. Rolla, 1980: 596-603.
[15] 林士渊. 规则节理岩体流变模式之研究[D]. 台北: 国立台湾大学土木工程学系, 2002.
LIN Shin-Yuan. The study of rheological model for regularly jointed rock masses[D]. Taipei: Department of Civil Engineering, National Taiwan University, 2002.
收稿日期：2008-06-20；修回日期：2008-09-10
基金项目：国家自然科学基金和雅砻江水电开发联合研究重点资助项目(50639090)
通信作者：熊良宵(1982-)，男，江西九江人，博士研究生，从事地下工程和岩石力学研究；电话：021-65987671；E-mail: xiongliangxiao@163.com
