基于節理巖體損傷本構的洞室位移反分析研究與應用

楊云浩，徐衛亞

(河海大學 巖土力學與堤壩工程教育部重點實驗室 巖土工程研究所，江蘇 南京，210098)

地下工程的開挖實際上是對處于某種應力狀態的巖體進行人為的卸荷與擾動。卸荷導致圍巖應力的重分布，使某些區域的巖體進入塑性工作狀態，塑性變形可能隨開挖的持續進行而不斷累積，從而導致巖體發生明顯的劣化而形成一定深度的擾動區。鑒于這一點，在洞室位移反分析中，必須對擾動區進行適當處理才有可能使所得結果更接近實際。考慮擾動區的位移反分析方法大多都基于楊林德等[1]較早提出的雙區域介質模型。雙區域介質模型因其簡單實用和可操作性強的特點得到了廣泛的應用[2-5]，但其不足之處在于未考慮擾動區內參數是漸變的且對于三維反分析問題存在數值建模上的困難，為此，倪紹虎等[6]提出基于“參數場”概念的反分析方法，該方法引入三維彈塑性損傷有限元分析，將圍巖擾動與損傷程度聯系起來，損傷程度以損傷降低系數來量化，擾動區內巖體參數值假設為未擾動前參數值與損傷降低系數的乘積，由此給出擾動區內的參數場，該方法相對比較合理。此外，江權[7]基于CWFS模型提出高地應力硬巖彈脆塑性劣化本構模型，該模型將擾動區內圍巖力學參數看作等效塑性剪應變的函數，因而也是漸變的。該本構模型可用于考慮擾動區的位移反分析，但由于參數的增加，反分析的復雜度和不唯一性也相應增大。對于在含成組展布的、相對結構尺寸為小量的斷續節理裂隙巖體內開挖的地下洞室，若不考慮爆破震動影響，則可以認為圍巖擾動區主要是因圍巖內部應力重分布導致原生節理裂隙(或裂紋)的進一步擴展所造成的，相應地，考慮擾動區的位移反分析應采用含節理巖體損傷本構模型較為適宜。孫衛軍等[8]采用二階對稱損傷張量考慮裂隙的力學影響，基于損傷力學原理建立這種節理巖體的本構模型，但是二階損傷張量在節理密度較大時往往會出現損傷變量大于 1的不合理情形。周維垣等[9]采用四階損傷張量，從自一致理論和即時模量概念出發推導節理巖體等效柔度張量，并由節理的擴展機理建立損傷張量演化方程，提出了相應的彈塑性損傷本構模型。藍航[10]和柴紅保等[11]均基于FLAC3D的莫爾-庫倫模型，實現了節理巖體損傷本構的開發。FLAC3D內置莫爾-庫倫模型只適用于各向同性巖土體，但節理巖體具有明顯的各向異性，其等效柔度張量反映了這一力學特征，這一點在文獻[10-11]的研究中并未得到足夠重視。在此，本文作者基于FLAC3D內置應變軟化本構模型，通過修改本構模型的數值計算格式使其可以適用于各向異性的節理巖體，同時采用FISH語言編制裂隙擴展(即損傷)演化程序，從而完成節理巖體損傷本構模型的開發。在此基礎上，采用支持向量回歸機(SVR)與粒子群-差異演化雜交優化算法(PSO-DE)相結合的位移反分析方法，以思林水電站地下廠房開挖期位移監測數據為依托，開展考慮開挖擾動劣化的位移反分析研究。

1 節理巖體損傷本構在 FLAC3D的實現

1.1 節理巖體的等效柔度張量

從損傷力學的觀點來看，節理裂隙可視為巖體的損傷，由此產生的附加柔度張量表示為Cd，則節理巖體的等效柔度張量Ce-d為[8,12]：

式中：Ce為巖石的柔度張量；E0為巖石的彈性模量；M為裂隙組數；體積密度和統計特征尺寸第k組裂隙面法向的方向余弦k組裂隙的傳壓和傳剪系數和G2為裂隙形狀系數，當將裂隙形狀簡化為圓盤時，

由式(2)可見：節理的存在使巖體的力學性質表現為明顯的各向異性(彈性常數個數大于 2)。在以FLAC3D內置應變軟化本構模型為基礎開發節理巖體損傷本構模型時，應首先從 FLAC3D所采用的最一般的彈塑性本構增量計算格式入手對該模型進行修改，使其可適用于各向異性的巖體。

1.2 FLAC3D內置應變軟化本構模型計算格式的修改

FLAC3D中的彈塑性本構模型都遵循同樣的增量計算算法格式，即：首先由 t時刻的應力以及時步Δt內產生的總應變增量，依據彈性剛度張量求得時步Δt內的應力增量以及t+Δt時刻的應力，這一步稱為彈性應力猜想，然后由屈服準則判斷是否出現塑性屈服，若是，則按照關聯或非關聯塑性流動法則對應力進行修正以獲取計算單元真正的應力狀態[13]。若以和別表示n維廣義應力矢量的第i個分量在計算時刻t的值及t+Δt時刻的彈性猜想值，則有[13]

式中： Δ ε1,… ,Δ εn為Δt時段內廣義應變增量矢量的n個分量； Si(Δ ε1,… ,Δ εn)為一線性函數，由彈性應力增量與應變增量關系確定。

如果應力點 σiI( i = 1 ,… ,n )位于屈服面以外，則需將其拉回至屈服面上，相應的應力修正式為

式中：f(·)和f*(·)分別為屈服準則函數及其去掉常數項以后的部分；g為塑性勢函數。

在 FLAC3D的應變軟化本構模型計算格式中，屈服準則函數以主應力形式表示，上述廣義應力和應變相應地變為主應力和主應變(此時，n=3)。在各向同性情況下，應力主向和應變主向是一致的，所以式(4a)中的 S1(·)可按下式計算(S2(·), S3(·)與此類似)

式中：K和 G分別為體積模量和剪切模量。將式(5)代入式(4a)和(4b)，即可得主應力修正值。在各向異性情況下，正應力不僅產生正應變，而且還會產生剪應變，此時式(5)便不成立了，因此必須以應力張量分量σij(i, j=1, 2, 3)修正的方式替代3個主應力修正的方式，相應地， Si(·)(i=1, …, 6)的計算式變為

式中： ( sij)6×6為剛度張量的矩陣形式。FLAC3D的應變軟化模型是基于莫爾-庫倫模型的(區別僅在于模型的黏聚力、內摩擦角、剪脹角和抗拉強度隨塑性應變的增加而軟化)，其屈服準則仍然是帶拉伸截止限的復合莫爾-庫倫屈服準則，模型的塑性勢函數g有2個，即剪切和拉伸塑性勢函數gs和gt。下面求gs和gt對應力張量分量的偏導。首先根據主應力與應力張量不變量的關系，將帶拉伸截止限的莫爾-庫倫屈服準則的剪切和拉伸屈服函數 fs和 ft以及gs和gt表達為下述形式：

式中：I1，J2，J3和 θσ分別為應力張量第一不變量、應力偏張量的第二、第三不變量和應力洛德角；c，φ，ψ，tσ分別為黏聚力、內摩擦角、剪脹角和抗拉強度。由式(7)可得gs和gt對ijσ的偏導為

J2，J3對ijσ的偏導容易求出，此處不贅述。將式(8)代入式(6)，再在式(4a)和(4b)中應用式(6)和(7)，便得到應力張量分量ijσ的修正式。

為檢驗上述做法的正確性，對 FLAC3D應變軟化本構程序的應力修正代碼段按上述方式進行了修改，應用修改后的應變軟化模型和原模型對各向同性材料進行單軸和三軸壓縮的數值模擬試驗。對比表明：由于應力修正方式的差異使得二者的計算結果并不完全相同，但差別是微小的(相對誤差在5%以內)，從而驗證了以應力張量分量修正代替主應力修正的可行性。這樣，就可以在修改后的應變軟化模型中應用前述節理巖體等效剛度張量計算彈性應力猜想值并進行正確的塑性修正，從而實現節理巖體損傷本構模型的開發。

1.3 多步開挖問題的等效柔度張量更新方法

對于多步開挖的地下洞室工程，每一步開挖都造成圍巖應力狀態的重新調整，從而可能導致分支裂隙的起裂和擴展(即損傷的演化)。為在等效柔度張量中考慮分支裂隙的貢獻，將原始裂隙與分支裂隙構成的拐折裂隙等效為一條直裂隙(這一簡化處理雖欠嚴謹但無疑是合理的)進行處理。等效直裂隙的平均半徑a*和裂隙面法向n*采用下式計算

式中：a和n分別為原始裂隙的平均半徑和裂隙面法向矢量；Δa和Δn與分支裂隙的起裂角和擴展長度 l有關(具體計算式見文獻[9])。l采用 Ashby等[14]提出的支裂紋起裂準則(見式(10a))與擴展長度計算式(見式(10b))確定。

式中：L=l/a；13/σσλ=；μ為摩擦系數；KIC為巖石Ⅰ型斷裂韌度。當單元的應力狀態滿足式(10a)時，則采用迭代計算確定l。

將a*和n*的計算過程編制為FISH程序，在每一步開挖模擬結束后調用該程序，則各計算單元的等效柔度張量即可更新并被用于下一步開挖模擬計算。

1.4 節理巖體損傷本構模型與常規彈塑性模型的比較

為對比所實現的節理巖體彈塑性損傷本構模型與常規彈塑性模型，進行單軸壓縮數值試驗。模擬試樣為圓柱形，加載采用應力控制，逐級施加荷載至塑性屈服出現。模擬試樣的力學參數取為：體積模量14.0 GPa，剪切模量10.5 GPa，黏聚力10.5 MPa，內摩擦角 50°，抗拉強度 1.0 MPa，斷裂韌度 1.8 MPa·m1/2。應用節理巖體損傷本構模型的試樣除采用相同的力學參數之外，還包括其內斷續節理(僅考慮一組節理裂隙)的統計幾何參數：裂隙平均半徑為2 cm(或4 cm)、體積密度0.03條/cm3、裂隙面法向矢量(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優化算法相結合的位移反分析流程

基于位移監測信息的力學參數反分析目標函數一般采用如下形式

式中：u和u分別為計算位移向量和監測位移向量，由各測點計算(監測)位移值構成；u由有限元(有限差分)計算求取。

為提高參數反分析的效率和可靠性，在反分析流程中(見圖 2)，采用支持向量回歸機(SVR)建立“參數-位移”非線性映射關系以代替正演過程中頻繁的數值計算工作，參數尋優過程則由PSO-DE雜交優化算法進行控制，該算法由 Mahamed等[15]提出，本文作者對其進行了適當改進后將其應用于位移反分析問題求解。

圖 2中，Gmax為粒子種群最大演化代數；G+1代第i個粒子位置向量的第j個分量；為第G代第i1和i2個粒子位置向量的第j個分量； pi3,j(G)為第G代第i3個粒子的自身最佳位置向量的第j個分量， i ≠ i1 ≠ i2 ≠i 3；和pgi,j(G )分別為第G代第i個粒子的自身和鄰居群最佳位置向量的第j個分量；uj,u1,j,u2,j～U(0,1)；RC( G,i)為交叉率，取值在[0,1]內； fi(G)為第G代第i個粒子的適應度；α，β和γ為相關參數。

3 工程實例應用

3.1 工程概況與地質條件

圖2 基于SVR與PSO-DE的位移反分析流程Fig.2 Flow chart of back-anaylsis based on SVR and PSO-DE algorithm

思林水電站位于貴州烏江中游河段，電站裝機 4臺，總裝機容量1 000 MW。地下主廠房洞室全長177.8 m，最大開挖寬度28.4 m，開挖總高度74.7 m，最大開挖斷面面積1 653 m2。主廠房軸線方位NW0°，軸線地面高程482～553 m，埋深90～120 m。據思林水電站地質可研編修報告，廠房區位于塘頭向斜倒轉翼，地層產狀N35°～42°E/NW∠70°。主廠房與主變室洞周以Ⅱ和Ⅲ類圍巖為主，巖性以三疊系夜郎組 T1y22中厚、厚層灰巖和白云質灰巖為主。廠房區內無大的斷層，但有2條規模較大的層間錯動帶，此外還有3組構造裂隙(見表1，表中的裂隙平均半徑及體積密度是將裂隙面簡化為圓盤狀，按文獻[10]采用的公式計算得來)，3組裂隙中以第Ⅰ組最為發育，占工程區統計量的45%；Ⅱ組次之，占30%；Ⅲ組為25%。廠房洞群區空間地應力量測表明：廠區最大主應力平均值約17 MPa，方位N60°W左右，傾角小于10°，符合測區區域構造應力場分布規律。

表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，向上延伸至地表。為實現后續的洞群穩定性分析的目的，所建洞群除主廠房、主變室及母線洞之外，還包括引水隧洞和尾水隧洞。模型內包含了洞群穿越的所有地層(P2w5，P2c，…，T1yn1)，以及2條層間錯動帶(fj1，fj2)，洞群與地層的空間關系如圖4(a)所示。模型計算網格劃分如圖4(b)所示。模型初始地應力場采用側壓系數法擬合獲得，模擬開挖步完全按實際分層開挖情況來確定。對主廠房所處的T1y22地層應用所開發的節理巖體損傷本構模型，所需的力學參數有彈性模量、泊松比、黏聚力初值和殘余值、內摩擦角初值和殘余值、巖石斷裂韌度，此外還有裂隙的幾何參數(其值見表1)。其余地層則應用莫爾-庫倫模型。

圖3 典型斷面多點位移計布設及主廠房開挖分層Fig.3 Excavation stages and layout of extensometers of powerhouse cavern

圖4 三維數值計算模型Fig.4 3D numerical model

由于反分析所選用的變形監測測孔均布置在主廠房內，而主廠房完全處于地層內(洞群區厚度最大的地層)，其上、下游相鄰地層分別為按地質報告，學性質上可視為同一地層，故反分析時只考慮3個地層的力學參數。但若3個地層的變形及強度參數均進行反分析，則將導致參數過多而使反分析無法實施，為此，首先采用正交試驗方差分析篩選對位移有顯著影響(即敏感性強)的參數。為盡量減少正交試驗方案總數，3個地層的變形參數只考慮彈性模量，這樣共有11個參數參與篩選。選用L12(211)正交表設計參數取值組合，代入數值模擬模型進行正演計算求得對應4個變形監測測孔處的增量位移數據后，即可采用方差分析確定敏感性強的參數，具體分析方法詳見文獻[16]。需說明的是，在正交表中未安排誤差列，在方差分析時將極差相對很小的列作為誤差列處理。綜合對應各測孔試驗數據的方差分析結果，最終確定的待反分析參數為：的反分析取值區間分別設定為20～35 GPa、7～13 GPa、18～24 GPa、1.0～1.6 MPa和 0.10～0.58 MPa。其余參數不進行反分析，取地質報告建議值，T1y22地層白云質灰巖的斷裂韌度 KIC則依據文獻[17]提出的巖石Ⅰ型斷裂韌度KIC與抗壓強度σc的線性經驗關系，確定為1.88 MPa·m1/2。

在參數敏感性分析基礎上，由均勻設計確定待反分析參數的取值方案，代入數值模擬模型計算后提取對應 4個測孔的增量相對位移，形成“參數-增量相對位移”樣本集供訓練 SVR使用。在正式反分析之前，首先以一組假想量測位移檢驗了前述位移反分析方法的可靠性，結果表明：敏感性最強的 E_T1y22反分析值與假想真實值的相對誤差為 0.87%～3.56%，的誤差較大，總體而言，反分析方法可靠。以真實量測位移反分析參數的過程見圖5，由此得到的E_T1y22，GPa，1.56 MPa，0.31 MPa，23.37 GPa 和 11.8 GPa。反分析所得參數對應的計算增量相對位移見表 2，與量測值相比，最大和最小相對誤差分別為 16.56%和7.54%，圖 5(a)所示為誤差平方和最終收斂至 0.098 mm2，而實際誤差平方和為0.152 mm2，二者比較接近，因此反分析成果是令人滿意的。為進行對比，所有地層均應用莫爾-庫倫本構模型，對參數進行了反分析，得到的T1y22地層等效彈性模量和黏聚力分別為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

將反分析所得參數代入數值計算模型進行思林主廠房及相鄰洞室分步開挖正演分析，得到的圍巖內裂隙相對擴展長度比動態變化的模擬結果如圖6所示。

圖6 裂隙相對擴展長度比隨開挖過程的動態變化Fig.6 Evolution of Ex_r during excavation

由圖6可見：第3步開挖結束后，發生裂隙擴展的區域范圍較小且擴展長度比也較小；第8步開挖結束后，發生裂隙擴展的區域顯著擴大且擴展長度比有所上升，在洞周應力集中部位(如頂拱)和卸荷松弛部位(如邊墻)均出現裂隙擴展；第 8步開挖后主廠房大規模開挖完成，轉入小規模槽挖階段，在此期間內，發生裂隙擴展的區域范圍增加不大。進一步將塑性區與裂隙相對擴展長度比等值線圖相對照，發現二者有較好的一致性，裂隙擴展區與塑性區范圍基本相當。

4 結論

(1) 對開挖于含成組展布的且相對結構尺寸為小量的斷續節理裂隙巖體內的地下洞室，開挖擾動區的形成與裂隙的擴展演化密切相關，只有采用節理巖體損傷本構模型才能夠較好地描述這一實際情況。基于FLAC3D內置應變軟化模型所開發的節理巖體損傷本構模型充分考慮節理巖體的各向異性特征，可將其用于考慮擾動區的位移反分析。

(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.