基于节理岩体损伤本构的洞室位移反分析研究与应用

杨云浩，徐卫亚

(河海大学 岩土力学与堤坝工程教育部重点实验室 岩土工程研究所，江苏 南京，210098)

摘要：基于对洞室位移反分析中开挖扰动区处理方法的总结分析，探讨将节理岩体损伤本构应用于反分析研究中的实现方法。通过FLAC^3D预留的二次开发接口，修改其内置应变软化模型的数值计算格式，即对位于屈服面以外的应力点采用修正应力张量分量的方式替代主应力修正的方式，使修改后的本构模型可适用于各向异性的节理岩体，同时结合FISH语言编制的节理裂隙几何参数随开挖过程的演化程序实现对节理岩体等效柔度张量的更新，从而完成节理岩体损伤本构的FLAC^3D二次开发。采用基于支持向量回归机与PSO-DE优化算法的增量位移反分析方法，应用该本构模型对思林水电站地下厂房围岩参数进行反分析。研究结果表明：对于开挖于含成组展布且相对结构尺寸为小量的断续节理岩体内的地下洞室位移反分析问题，将节理岩体损伤本构引入其中是一个很好的尝试，扰动区内围岩的等效弹模平均降低约为20.13%，计算增量位移值与量测值的相对误差在7.54%~16.56%，反分析成果是合理和可靠的。

关键词：节理岩体；开挖扰动区；损伤；位移反分析；FLAC^3D自定义本构

中图分类号：TU923         文献标志码：A         文章编号：1672-7207(2012)07-2723-10

Displacement-based back-analysis using damage mechanics constitutive model for jointed rock mass and its application in underground opening

YANG Yun-hao, XU Wei-ya

(Key Laboratory of Ministry of Education for Geomechanics and Embankment Engineering,

Research Institute of Geotechnical Engineering, Hohai University, Nanjing 210098, China)

Abstract: Based on the overview of back analysis methods for excavation disturbance zone (EDZ), the implementation of employing damage mechanics constitutive model for jointed rock mass (JRM) in displacement-based back analysis is discussed. Following the guide of writing constitutive model in FLAC^3D, the damage mechanic constitutive model for JRM is developed successfully by modifying the basic strain-softening model. Specifically, in order to make the strain-softening model capable of handling the case when elastic flexibility tensor is anisotropic such as the equivalent elastic compliance tensor of JRM, a correction is made to the six components of stress tensor instead of the three principal stresses for the stress point located above the yield surface. Meanwhile, a FISH program is written to update characteristic parameters of joint sets, so the equivalent elastic compliance tensor of each zone in numerical model can be updated accordingly. Applying the damage mechanics constitutive model to numerical simulation of underground opening excavation of Silin hydropower station, the back-analysis of rock mass deformation and strength parameters is conducted using the displacement monitoring data from multi-point extensometers. The result shows that for underground opening excavated within jointed rock mass containing a considerable number of intermittent joints with stable trending, employing the damage mechanics constitutive model for JRM in back-analysis is a good attempt, the influence of EDZ can be taken into consideration reasonably, the average decrease of equivalent elastic modulus of rock mass in EDZ is around 20.13%, the fractional errors of calculated incremental displacements are in the range of 7.54%-16.56%, and the result of back-analysis is reliable.

Key words: jointed rock mass; excavation disturbance zone; damage; displacement-based back-analysis; FLAC^3D user-defined constitutive model

地下工程的开挖实际上是对处于某种应力状态的岩体进行人为的卸荷与扰动。卸荷导致围岩应力的重分布，使某些区域的岩体进入塑性工作状态，塑性变形可能随开挖的持续进行而不断累积，从而导致岩体发生明显的劣化而形成一定深度的扰动区。鉴于这一点，在洞室位移反分析中，必须对扰动区进行适当处理才有可能使所得结果更接近实际。考虑扰动区的位移反分析方法大多都基于杨林德等^[1]较早提出的双区域介质模型。双区域介质模型因其简单实用和可操作性强的特点得到了广泛的应用^[2-5]，但其不足之处在于未考虑扰动区内参数是渐变的且对于三维反分析问题存在数值建模上的困难，为此，倪绍虎等^[6]提出基于“参数场”概念的反分析方法，该方法引入三维弹塑性损伤有限元分析，将围岩扰动与损伤程度联系起来，损伤程度以损伤降低系数来量化，扰动区内岩体参数值假设为未扰动前参数值与损伤降低系数的乘积，由此给出扰动区内的参数场，该方法相对比较合理。此外，江权^[7]基于CWFS模型提出高地应力硬岩弹脆塑性劣化本构模型，该模型将扰动区内围岩力学参数看作等效塑性剪应变的函数，因而也是渐变的。该本构模型可用于考虑扰动区的位移反分析，但由于参数的增加，反分析的复杂度和不唯一性也相应增大。对于在含成组展布的、相对结构尺寸为小量的断续节理裂隙岩体内开挖的地下洞室，若不考虑爆破震动影响，则可以认为围岩扰动区主要是因围岩内部应力重分布导致原生节理裂隙(或裂纹)的进一步扩展所造成的，相应地，考虑扰动区的位移反分析应采用含节理岩体损伤本构模型较为适宜。孙卫军等^[8]采用二阶对称损伤张量考虑裂隙的力学影响，基于损伤力学原理建立这种节理岩体的本构模型，但是二阶损伤张量在节理密度较大时往往会出现损伤变量大于1的不合理情形。周维垣等^[9]采用四阶损伤张量，从自一致理论和即时模量概念出发推导节理岩体等效柔度张量，并由节理的扩展机理建立损伤张量演化方程，提出了相应的弹塑性损伤本构模型。蓝航^[10]和柴红保等^[11]均基于FLAC^3D的莫尔-库伦模型，实现了节理岩体损伤本构的开发。FLAC^3D内置莫尔-库伦模型只适用于各向同性岩土体，但节理岩体具有明显的各向异性，其等效

柔度张量反映了这一力学特征，这一点在文献[10-11]的研究中并未得到足够重视。在此，本文作者基于FLAC^3D内置应变软化本构模型，通过修改本构模型的数值计算格式使其可以适用于各向异性的节理岩体，同时采用FISH语言编制裂隙扩展(即损伤)演化程序，从而完成节理岩体损伤本构模型的开发。在此基础上，采用支持向量回归机(SVR)与粒子群-差异演化杂交优化算法(PSO-DE)相结合的位移反分析方法，以思林水电站地下厂房开挖期位移监测数据为依托，开展考虑开挖扰动劣化的位移反分析研究。

1  节理岩体损伤本构在FLAC3D的实现

1.1  节理岩体的等效柔度张量

从损伤力学的观点来看，节理裂隙可视为岩体的损伤，由此产生的附加柔度张量表示为C^d，则节理岩体的等效柔度张量C^e-d为^{[8, 12]}：

               (1)

    (2)

式中：C^e为岩石的柔度张量；E₀为岩石的弹性模量；M为裂隙组数；和分别为第k组裂隙的平均体积密度和统计特征尺寸；，，和为第k组裂隙面法向的方向余弦；和分别为第k组裂隙的传压和传剪系数，，；G₁和G₂为裂隙形状系数，当将裂隙形状简化为圆盘时，

，。

由式(2)可见：节理的存在使岩体的力学性质表现为明显的各向异性(弹性常数个数大于2)。在以FLAC^3D内置应变软化本构模型为基础开发节理岩体损伤本构模型时，应首先从FLAC^3D所采用的最一般的弹塑性本构增量计算格式入手对该模型进行修改，使其可适用于各向异性的岩体。

1.2  FLAC^3D内置应变软化本构模型计算格式的修改

FLAC^3D中的弹塑性本构模型都遵循同样的增量计算算法格式，即：首先由t时刻的应力以及时步?t内产生的总应变增量，依据弹性刚度张量求得时步?t内的应力增量以及t+?t时刻的应力，这一步称为弹性应力猜想，然后由屈服准则判断是否出现塑性屈服，若是，则按照关联或非关联塑性流动法则对应力进行修正以获取计算单元真正的应力状态^[13]。若以和分别表示n维广义应力矢量的第i个分量在计算时刻t的值及t+?t时刻的弹性猜想值，则有^[13]

 (i=1, …, n)      (3)

式中：为?t时段内广义应变增量矢量的n个分量；为一线性函数，由弹性应力增量与应变增量关系确定。

如果应力点位于屈服面以外，则需将其拉回至屈服面上，相应的应力修正式为

 (i=1, …, n)    (4a)

  (4b)

式中：和分别为屈服准则函数及其去掉常数项以后的部分；g为塑性势函数。

在FLAC^3D的应变软化本构模型计算格式中，屈服准则函数以主应力形式表示，上述广义应力和应变相应地变为主应力和主应变(此时，n=3)。在各向同性情况下，应力主向和应变主向是一致的，所以式(4a)中的可按下式计算(, 与此类似)

     (5)

式中：K和G分别为体积模量和剪切模量。将式(5)代入式(4a)和(4b)，即可得主应力修正值。在各向异性情况下，正应力不仅产生正应变，而且还会产生剪应变，此时式(5)便不成立了，因此必须以应力张量分量(i, j=1, 2, 3)修正的方式替代3个主应力修正的方式，相应地，(i=1, …, 6)的计算式变为

  (6)

式中：为刚度张量的矩阵形式。FLAC^3D的应变软化模型是基于莫尔-库伦模型的(区别仅在于模型的黏聚力、内摩擦角、剪胀角和抗拉强度随塑性应变的增加而软化)，其屈服准则仍然是带拉伸截止限的复合莫尔-库伦屈服准则，模型的塑性势函数g有2个，即剪切和拉伸塑性势函数g^s和g^t。下面求g^s和g^t对应力张量分量的偏导。首先根据主应力与应力张量不变量的关系，将带拉伸截止限的莫尔-库伦屈服准则的剪切和拉伸屈服函数和以及g^s和g^t表达为下述形式：

 (7)

式中：I₁，J₂，J₃和θ_σ分别为应力张量第一不变量、应力偏张量的第二、第三不变量和应力洛德角；c，，，分别为黏聚力、内摩擦角、剪胀角和抗拉强度。由式(7)可得g^s和g^t对的偏导为

，，    (8)

J₂，J₃对的偏导容易求出，此处不赘述。将式(8)代入式(6)，再在式(4a)和(4b)中应用式(6)和(7)，便得到应力张量分量的修正式。

为检验上述做法的正确性，对FLAC^3D应变软化本构程序的应力修正代码段按上述方式进行了修改，应用修改后的应变软化模型和原模型对各向同性材料进行单轴和三轴压缩的数值模拟试验。对比表明：由于应力修正方式的差异使得二者的计算结果并不完全相同，但差别是微小的(相对误差在5%以内)，从而验证了以应力张量分量修正代替主应力修正的可行性。这样，就可以在修改后的应变软化模型中应用前述节理岩体等效刚度张量计算弹性应力猜想值并进行正确的塑性修正，从而实现节理岩体损伤本构模型的开发。

1.3  多步开挖问题的等效柔度张量更新方法

对于多步开挖的地下洞室工程，每一步开挖都造成围岩应力状态的重新调整，从而可能导致分支裂隙的起裂和扩展(即损伤的演化)。为在等效柔度张量中考虑分支裂隙的贡献，将原始裂隙与分支裂隙构成的拐折裂隙等效为一条直裂隙(这一简化处理虽欠严谨但无疑是合理的)进行处理。等效直裂隙的平均半径a^*和裂隙面法向n^*采用下式计算

，          (9)

式中：a和n分别为原始裂隙的平均半径和裂隙面法向矢量；?a和?n与分支裂隙的起裂角和扩展长度l有关(具体计算式见文献[9])。l采用Ashby等^[14]提出的支裂纹起裂准则(见式(10a))与扩展长度计算式(见式(10b))确定。

     (10a)

       (10b)

式中：L=l/a；；为摩擦系数；K_IC为岩石Ⅰ型断裂韧度。当单元的应力状态满足式(10a)时，则采用迭代计算确定l。

将a^*和n^*的计算过程编制为FISH程序，在每一步开挖模拟结束后调用该程序，则各计算单元的等效柔度张量即可更新并被用于下一步开挖模拟计算。

1.4  节理岩体损伤本构模型与常规弹塑性模型的  比较

为对比所实现的节理岩体弹塑性损伤本构模型与常规弹塑性模型，进行单轴压缩数值试验。模拟试样为圆柱形，加载采用应力控制，逐级施加荷载至塑性屈服出现。模拟试样的力学参数取为：体积模量14.0 GPa，剪切模量10.5 GPa，黏聚力10.5 MPa，内摩擦角50°，抗拉强度1.0 MPa，断裂韧度1.8 MPa·m^1/2。应用节理岩体损伤本构模型的试样除采用相同的力学参数之外，还包括其内断续节理(仅考虑一组节理裂隙)的统计几何参数：裂隙平均半径为2 cm(或4 cm)、体积密度0.03条/cm³、裂隙面法向矢量(0.68, 0.27, 0.68)。数值试验得到的轴向应变-应力曲线如图1所示。

图1  单轴压缩数值试验轴向应变-应力关系

Fig.1  Axial strain versus axial stress from uniaxial compression test

图1中，“弹塑性损伤模型-1”和“弹塑性损伤模  型-2”分别对应裂隙半径为2 cm和4 cm时的情况。由图1可见：弹塑性损伤模型的曲线斜率略大于常规模型的曲线斜率，且表现为随荷载的增加而增大(荷载较小时增幅不明显，荷载较大时较为显著)；与常规弹塑性模型相比，损伤模型在较低的应力下进入塑性屈服，裂隙半径越大(即等效柔度越大)，则相应的屈服点应力越低。数值试验表明：所实现的弹塑性损伤本构模型能正确地描述含断续节理岩体在初始损伤及变载作用下的损伤演化所导致的力学性质的弱化现象。

2  支持向量机与PSO-DE优化算法相结合的位移反分析流程

基于位移监测信息的力学参数反分析目标函数一般采用如下形式

         (11)

式中：u和分别为计算位移向量和监测位移向量，由各测点计算(监测)位移值构成；u由有限元(有限差分)计算求取。

为提高参数反分析的效率和可靠性，在反分析流程中(见图2)，采用支持向量回归机(SVR)建立“参   数-位移”非线性映射关系以代替正演过程中频繁的数值计算工作，参数寻优过程则由PSO-DE杂交优化算法进行控制，该算法由Mahamed等^[15]提出，本文作者对其进行了适当改进后将其应用于位移反分析问题求解。

图2中，G_max为粒子种群最大演化代数；为第G+1代第i个粒子位置向量的第j个分量；和分别为第G代第i1和i2个粒子位置向量的第j个分量；为第G代第i3个粒  子的自身最佳位置向量的第j个分量，；，和分别为第G代第i个粒子的自身和邻居群最佳位置向量的第j个分量；；为交叉率，取值在[0,1]内；为第G代第i个粒子的适应度；，和为相关参数。

3  工程实例应用

3.1  工程概况与地质条件

思林水电站位于贵州乌江中游河段，电站装机4台，总装机容量1 000 MW。地下主厂房洞室全长177.8 m，最大开挖宽度28.4 m，开挖总高度74.7 m，最大开挖断面面积1 653 m²。主厂房轴线方位NW0°，轴线地面高程482~553 m，埋深90~120 m。据思林水电站地质可研编修报告，厂房区位于塘头向斜倒转翼，地层产状N35°~42°E/NW∠70°。主厂房与主变室洞周以Ⅱ和Ⅲ类围岩为主，岩性以三叠系夜郎组T₁y²²中厚、厚层灰岩和白云质灰岩为主。厂房区内无大的断层，但有2条规模较大的层间错动带，此外还有3组构造裂隙(见表1，表中的裂隙平均半径及体积密度是将裂隙面简化为圆盘状，按文献[10]采用的公式计算得来)，3组裂隙中以第Ⅰ组最为发育，占工程区统计量的45%；Ⅱ组次之，占30%；Ⅲ组为25%。厂房洞群区空间地应力量测表明：厂区最大主应力平均值约17 MPa，方位N60°W左右，倾角小于10°，符合测区区域构造应力场分布规律。

图2  基于SVR与PSO-DE的位移反分析流程

Fig.2  Flow chart of back-anaylsis based on SVR and PSO-DE algorithm

表1  厂房区节理分组及其几何特征参数

Table 1  Characteristic geometric parameters of joint sets in powerhouse area

3.2  开挖支护及位移监测

主厂房共分11层开挖施工，分层高度5~8 m(见图3)。为减小爆破对边墙围岩的扰动，每一分层采用分块开挖的方式，首先开挖中心区块然后扩挖上下游边墙2 m厚的保护区块，每层开挖结束后立即进行喷混凝土和系统锚杆支护。为监测主厂房的围岩变形，沿其轴向共布置有4个位移监测断面，随开挖的进行逐步安设多点位移计，多点位移计测孔布置4个锚点，锚点间距不一，最深锚点距孔口24 m。2-2断面变形监测布置如图3所示。由于有相当数量的测孔量测数据不能很好地反映围岩的变形规律，因而真正用于反分析的测孔只有4个：1-1断面的MZC1-5；2-2断面的MZC2-4和MZC2-5；3-3断面的MZC3-4。提取的量测数据为测孔安装之后至厂房开挖完毕这一期间所发生的孔口相对孔底的位移，即增量相对位移。

3.3  数值模拟模型与位移反分析

数值模拟模型是进行参数反分析的关键，在分析地勘资料的基础上，建立了厂房洞群的三维数值模型。模型长480 m，宽280 m，模型底部标高EL230.0，向上延伸至地表。为实现后续的洞群稳定性分析的目的，所建洞群除主厂房、主变室及母线洞之外，还包括引水隧洞和尾水隧洞。模型内包含了洞群穿越的所有地层(P₂w⁵，P₂c，…，T₁yn¹)，以及2条层间错动带(fj1，fj2)，洞群与地层的空间关系如图4(a)所示。模型计算网格划分如图4(b)所示。模型初始地应力场采用侧压系数法拟合获得，模拟开挖步完全按实际分层开挖情况来确定。对主厂房所处的T₁y²²地层应用所开发的节理岩体损伤本构模型，所需的力学参数有弹性模量、泊松比、黏聚力初值和残余值、内摩擦角初值和残余值、岩石断裂韧度，此外还有裂隙的几何参数(其值见表1)。其余地层则应用莫尔-库伦模型。

图3  典型断面多点位移计布设及主厂房开挖分层

Fig.3  Excavation stages and layout of extensometers of powerhouse cavern

图4  三维数值计算模型

Fig.4  3D numerical model

由于反分析所选用的变形监测测孔均布置在主厂房内，而主厂房完全处于T₁y²²地层内(洞群区厚度最大的地层)，其上、下游相邻地层分别为T₁y²¹²、T₁y²³¹及T₁y²³²，按地质报告，T₁y²³²和T₁y²¹²在力学性质上可视为同一地层，故反分析时只考虑3个地层的力学参数。但若3个地层的变形及强度参数均进行反分析，则将导致参数过多而使反分析无法实施，为此，首先采用正交试验方差分析筛选对位移有显著影响(即敏感性强)的参数。为尽量减少正交试验方案总数，3个地层的变形参数只考虑弹性模量，这样共有11个参数参与筛选。选用L₁₂(2¹¹)正交表设计参数取值组合，代入数值模拟模型进行正演计算求得对应4个变形监测测孔处的增量位移数据后，即可采用方差分析确定敏感性强的参数，具体分析方法详见文献[16]。需说明的是，在正交表中未安排误差列，在方差分析时将极差相对很小的列作为误差列处理。综合对应各测孔试验数据的方差分析结果，最终确定的待反分析参数为：T₁y²²地层弹性模量E_T₁y²²、初始黏聚力c_i_T₁y²²、残余黏聚力c_r_T₁y²²，T₁y²³¹地层弹性模量E_T₁y²³¹，T₁y²³²(T₁y²¹²)地层弹性模量E_T₁y²³²。参照地质报告，将E_T₁y²²，E_T₁y²³¹，E_T₁y²³²，c_i_T₁y²²和c_r_T₁y²²的反分析取值区间分别设定为20~35 GPa、7~13 GPa、18~24 GPa、1.0~1.6 MPa和0.10~0.58 MPa。其余参数不进行反分析，取地质报告建议值，T₁y²²地层白云质灰岩的断裂韧度K_IC则依据文献[17]提出的岩石Ⅰ型断裂韧度K_IC与抗压强度σ_c的线性经验关系，确定为1.88 MPa·m^1/2。

在参数敏感性分析基础上，由均匀设计确定待反分析参数的取值方案，代入数值模拟模型计算后提取对应4个测孔的增量相对位移，形成“参数-增量相对位移”样本集供训练SVR使用。在正式反分析之前，首先以一组假想量测位移检验了前述位移反分析方法的可靠性，结果表明：敏感性最强的E_T₁y²²反分析值与假想真实值的相对误差为0.87%~3.56%，c_i_T₁y²²的误差为0.63%~6.25%，E_T₁y²³²的误差为0.40%~7.11%，而敏感性较弱的c_r_T₁y²²和E_T₁y²³¹的误差较大，总体而言，反分析方法可靠。以真实量测位移反分析参数的过程见图5，由此得到的E_T₁y²²，c_i_T₁y²²，c_r_T₁y²²，E_T₁y²³²和E_T₁y²³¹值分别为：31.78 GPa，1.56 MPa，0.31 MPa，23.37 GPa和11.8 GPa。反分析所得参数对应的计算增量相对位移见表2，与量测值相比，最大和最小相对误差分别为16.56%和7.54%，图5(a)所示为误差平方和最终收敛至0.098 mm²，而实际误差平方和为0.152 mm²，二者比较接近，因此反分析成果是令人满意的。为进行对比，所有地层均应用莫尔-库伦本构模型，对参数进行了反分析，得到的T₁y²²地层等效弹性模量和黏聚力分别为25.38 GPa和1.30 MPa，二者对比可见：采用考虑损伤的本构对应的参数值明显大于等效参数值，这是因为损伤模型中的弹性模量为岩石的弹性模量，该值理所当然的应比岩体的要大，由此可以近似认为开挖扰动后围岩的弹性模量平均降低20.13%，靠近开挖边界处围岩的等效弹模降低值应高于该值，而在远离开挖边界处，则应远小于该值。

图5  均方误差收敛曲线及参数识别过程曲线

Fig.5  Processes of SSE convergence and parameters identification

表2   测孔孔口处计算增量相对位移值与实测值对比

Table 2   Incremental displacements from back analysis and field measurement

将反分析所得参数代入数值计算模型进行思林主厂房及相邻洞室分步开挖正演分析，得到的围岩内裂隙相对扩展长度比Ex_r ()动态变化的模拟结果如图6所示。

图6  裂隙相对扩展长度比随开挖过程的动态变化

Fig.6  Evolution of Ex_r during excavation

由图6可见：第3步开挖结束后，发生裂隙扩展的区域范围较小且扩展长度比也较小；第8步开挖结束后，发生裂隙扩展的区域显著扩大且扩展长度比有所上升，在洞周应力集中部位(如顶拱)和卸荷松弛部位(如边墙)均出现裂隙扩展；第8步开挖后主厂房大规模开挖完成，转入小规模槽挖阶段，在此期间内，发生裂隙扩展的区域范围增加不大。进一步将塑性区与裂隙相对扩展长度比等值线图相对照，发现二者有较好的一致性，裂隙扩展区与塑性区范围基本相当。

4  结论

(1) 对开挖于含成组展布的且相对结构尺寸为小量的断续节理裂隙岩体内的地下洞室，开挖扰动区的形成与裂隙的扩展演化密切相关，只有采用节理岩体损伤本构模型才能够较好地描述这一实际情况。基于FLAC^3D内置应变软化模型所开发的节理岩体损伤本构模型充分考虑节理岩体的各向异性特征，可将其用于考虑扰动区的位移反分析。

(2) 采用节理岩体损伤本构的数值计算虽然不能给出一个直观的围岩体内渐变参数场，但是由于每个单元体的裂隙扩展长度不同，所以单元的刚度张量也不同，这相当于间接地给出一个渐变参数场。

(3) 节理岩体损伤本构模型参数较多，若所有参数都进行反分析，将导致反分析解的不唯一性问题更加突出，因此，在反分析过程中通过参数敏感性分析减少待反分析参数的个数是必要的。此外，作为替代数值计算的“参数-位移”非线性映射建模工具和正演反分析寻优控制算法都对反分析解的可靠性有决定性影响，采用SVR和PSO-DE优化算法进行反分析可在一定程度上提高反分析解的可靠性。

参考文献：

[1] 杨林德. 岩土工程问题的反演理论与工程实践[M]. 北京: 科学出版社, 1996: 101-106.
YANG Lin-de. Theory of back analysis in geotechnical engineering and its application[M]. Beijing: Science Press, 1996: 101-106.

[2] 侯哲生, 韩文峰, 李晓. 金川二矿区巷道围岩力学参数空间效应数值分析[J]. 岩石力学与工程学报, 2004, 23(1): 64-68.
HOU Zhe-sheng, HAN Wen-feng, LI Xiao. Numerical analysis on spatial effect of mechanical parameter of surrounding rocks in Jinchuan depositⅡ[J]. Chinese Journal of Rock Mechanics and Engineering, 2004, 23(1): 64-68.

[3] 李宁, 段小强, 陈方方, 等. 围岩松动圈的弹塑性位移反分析方法探索[J]. 岩石力学与工程学报, 2006, 25(7): 1304-1308.

LI Ning, DUAN Xiao-qiang, CHEN Fang-fang, et al. A back analysis method for elasto-plastic displacement of broken rock zone around tunnel[J]. Chinese Journal of Rock Mechanics and Engineering, 2006, 25(7): 1304-1308.

[4] 高广运, 冯世进, 李鸿博. 卵砾石层中隧道考虑松动圈的位移反分析应用[J]. 地下空间与工程学报, 2008, 4(1): 57-61.

GAO Guang-yun, FENG Shi-jin, LI Hong-bo. Application of displacement back analysis for tunnel in gravel-cobble stratum with considering the loosing zone[J]. Chinese Journal of Underground Space and Engineering, 2008, 4(1): 57-61.

[5] 魏进兵, 邓建辉. 高地应力条件下大型地下厂房松动区变化规律及参数反演[J]. 岩土力学, 2010, 31(增刊1): 330-336.

WEI Jin-bing, DENG Jian-hui. Variation of excavation damaged zone and back analysis of large scale underground powerhouse with high geostress[J]. Rock and Soil Mechanics, 2010, 31(Supp.1): 330-336.

[6] 倪绍虎, 肖明. 基于围岩松动圈的地下工程参数场位移反分析[J]. 岩石力学与工程学报, 2009, 28(7): 1439-1446.

NI Shao-hu, XIAO Ming. Displacement back analysis of parameter field in underground engineering based on excavation damaged zone of surrounding rock[J]. Chinese Journal of Rock Mechanics and Engineering, 2009, 28(7): 1439-1446.

[7] 江权. 高地应力硬岩弹脆塑性劣化本构模型与大型地下洞室群围岩稳定性分析[D]. 武汉: 中科院武汉岩土力学研究所, 2007: 13-34.

JIANG Quan. Study on model and stability of surrounding rock of large underground caverns under high geo-stress condition[D]. Institute of Rock and Soil Mechanics of Chinese Academy of Science, 2007: 13-34.

[8] 孙卫军, 周维垣. 裂隙岩体弹塑性-损伤本构模型[J]. 岩石力学与工程学报, 1990, 9(2): 108-119.

SUN Wei-jun, ZHOU Wei-yuan. An elasto-plastic damage mechanics constitutive model for jointed rockmass[J]. Chinese Journal of Rock Mechanics and Engineering, 1990, 9(2): 108-119.

[9] 周维垣, 杨延毅. 节理岩体的损伤断裂模型及验证[J]. 岩石力学与工程学报, 1991, 10(1): 43-54.

ZHOU Wei-yuan, YANG Yan-yi. A damage fracture mechanics model for jointed rock masses and its verification[J]. Chinese Journal of Rock Mechanics and Engineering, 1991, 10(1): 43-54.

[10] 蓝航. 节理岩体采动损伤本构模型及其在露井联采工程中的应用[D]. 北京: 煤炭科学研究总院, 2007: 11-48.

LAN Hang. Constitutive model of joint rock mass damage due to mining and its application in engineering of open-underground combined mining[D]. Beijing: China Coal Research Institute, 2007: 11-48.

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

CHAI Hong-bao, CAO Ping, ZHAO Yan-lin, et al. Implementation and application of constitutive model for damage evolution of fractured rock mass[J]. Chinese Journal of Geotechnical Engineering, 2010, 32(7): 1047-1053.

[12] 朱维申, 李术才, 陈卫忠, 等. 节理岩体破坏机理和锚固效应及工程应用[M]. 北京: 科学出版社, 2002: 53-72.
ZHU Wei-shen, LI Shu-cai, CHEN Wei-zhong, et al. Failure mechanism for jointed rock mass, effect of anchoring and engineering applications[M]. Beijing: Science Press, 2002: 53-72.

[13] Itasca Software Company. Theory and back ground, constitutive models: Theory and implementation[M]. Minneapolis Minnesota USA: Itasca Software Company, 2003: 1-49.

[14] Ashby M F, Hallam S D. The failure of brittle solids containing small cracks under compressive stress states[J]. Acta Metal, 1986, 34(3): 497-510.

[15] Mahamed G H, Omran A P, Engelbrecht A S. Bare bones differential evolution[J]. European Journal of Operational Research, 2009, 196: 128-139.

[16] 吴群英, 林亮. 应用数理统计[M]. 天津: 天津大学出版社, 2004: 149-160.
WU Qun-ying, LIN Liang. Applied mathematical statistics[M]. Tianjin: Tianjin University Press, 2004: 149-160.

[17] 李江腾, 古德生, 曹平, 等. 岩石断裂韧度与抗压强度的相关规律[J]. 中南大学学报: 自然科学版, 2009, 40(6): 1695-1699.

LI Jiang-teng, GU De-sheng, CAO Ping, et al. Interrelated law between mode (Ⅰ): Fracture toughness and compression strength of rock[J]. Journal of Central South University: Science and Technology, 2009, 40(6): 1695-1699.

(编辑 何运斌)

收稿日期：2011-06-14；修回日期：2011-08-22

基金项目：国家重点基础研究发展计划(“973”计划)项目(2011CB013504)；国家自然科学基金资助项目(50979030)

通信作者：杨云浩(1977-)，男，山西黎城人，博士研究生，从事岩石力学与工程方面的研究工作；电话：13585389375；E-mail: haoyunyang2001@163.com