基于節(jié)理巖體損傷本構(gòu)的洞室位移反分析研究與應(yīng)用

楊云浩，徐衛(wèi)亞

(河海大學(xué) 巖土力學(xué)與堤壩工程教育部重點(diǎn)實(shí)驗(yàn)室 巖土工程研究所，江蘇 南京，210098)

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

1 節(jié)理巖體損傷本構(gòu)在 FLAC3D的實(shí)現(xiàn)

1.1 節(jié)理巖體的等效柔度張量

從損傷力學(xué)的觀點(diǎn)來(lái)看，節(jié)理裂隙可視為巖體的損傷，由此產(chǎn)生的附加柔度張量表示為Cd，則節(jié)理巖體的等效柔度張量Ce-d為[8,12]：

式中：Ce為巖石的柔度張量；E0為巖石的彈性模量；M為裂隙組數(shù)；體積密度和統(tǒng)計(jì)特征尺寸第k組裂隙面法向的方向余弦k組裂隙的傳壓和傳剪系數(shù)和G2為裂隙形狀系數(shù)，當(dāng)將裂隙形狀簡(jiǎn)化為圓盤(pán)時(shí)，

由式(2)可見(jiàn)：節(jié)理的存在使巖體的力學(xué)性質(zhì)表現(xiàn)為明顯的各向異性(彈性常數(shù)個(gè)數(shù)大于 2)。在以FLAC3D內(nèi)置應(yīng)變軟化本構(gòu)模型為基礎(chǔ)開(kāi)發(fā)節(jié)理巖體損傷本構(gòu)模型時(shí)，應(yīng)首先從 FLAC3D所采用的最一般的彈塑性本構(gòu)增量計(jì)算格式入手對(duì)該模型進(jìn)行修改，使其可適用于各向異性的巖體。

1.2 FLAC3D內(nèi)置應(yīng)變軟化本構(gòu)模型計(jì)算格式的修改

FLAC3D中的彈塑性本構(gòu)模型都遵循同樣的增量計(jì)算算法格式，即：首先由 t時(shí)刻的應(yīng)力以及時(shí)步Δt內(nèi)產(chǎn)生的總應(yīng)變?cè)隽?，依?jù)彈性剛度張量求得時(shí)步Δt內(nèi)的應(yīng)力增量以及t+Δt時(shí)刻的應(yīng)力，這一步稱為彈性應(yīng)力猜想，然后由屈服準(zhǔn)則判斷是否出現(xiàn)塑性屈服，若是，則按照關(guān)聯(lián)或非關(guān)聯(lián)塑性流動(dòng)法則對(duì)應(yīng)力進(jìn)行修正以獲取計(jì)算單元真正的應(yīng)力狀態(tài)[13]。若以和別表示n維廣義應(yīng)力矢量的第i個(gè)分量在計(jì)算時(shí)刻t的值及t+Δt時(shí)刻的彈性猜想值，則有[13]

式中： Δ ε1,… ,Δ εn為Δt時(shí)段內(nèi)廣義應(yīng)變?cè)隽渴噶康膎個(gè)分量； Si(Δ ε1,… ,Δ εn)為一線性函數(shù)，由彈性應(yīng)力增量與應(yīng)變?cè)隽筷P(guān)系確定。

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

式中：f(·)和f*(·)分別為屈服準(zhǔn)則函數(shù)及其去掉常數(shù)項(xiàng)以后的部分；g為塑性勢(shì)函數(shù)。

在 FLAC3D的應(yīng)變軟化本構(gòu)模型計(jì)算格式中，屈服準(zhǔn)則函數(shù)以主應(yīng)力形式表示，上述廣義應(yīng)力和應(yīng)變相應(yīng)地變?yōu)橹鲬?yīng)力和主應(yīng)變(此時(shí)，n=3)。在各向同性情況下，應(yīng)力主向和應(yīng)變主向是一致的，所以式(4a)中的 S1(·)可按下式計(jì)算(S2(·), S3(·)與此類似)

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

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

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

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

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

1.3 多步開(kāi)挖問(wèn)題的等效柔度張量更新方法

對(duì)于多步開(kāi)挖的地下洞室工程，每一步開(kāi)挖都造成圍巖應(yīng)力狀態(tài)的重新調(diào)整，從而可能導(dǎo)致分支裂隙的起裂和擴(kuò)展(即損傷的演化)。為在等效柔度張量中考慮分支裂隙的貢獻(xiàn)，將原始裂隙與分支裂隙構(gòu)成的拐折裂隙等效為一條直裂隙(這一簡(jiǎn)化處理雖欠嚴(yán)謹(jǐn)?shù)珶o(wú)疑是合理的)進(jìn)行處理。等效直裂隙的平均半徑a*和裂隙面法向n*采用下式計(jì)算

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

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

將a*和n*的計(jì)算過(guò)程編制為FISH程序，在每一步開(kāi)挖模擬結(jié)束后調(diào)用該程序，則各計(jì)算單元的等效柔度張量即可更新并被用于下一步開(kāi)挖模擬計(jì)算。

1.4 節(jié)理巖體損傷本構(gòu)模型與常規(guī)彈塑性模型的比較

為對(duì)比所實(shí)現(xiàn)的節(jié)理巖體彈塑性損傷本構(gòu)模型與常規(guī)彈塑性模型，進(jìn)行單軸壓縮數(shù)值試驗(yàn)。模擬試樣為圓柱形，加載采用應(yīng)力控制，逐級(jí)施加荷載至塑性屈服出現(xiàn)。模擬試樣的力學(xué)參數(shù)取為：體積模量14.0 GPa，剪切模量10.5 GPa，黏聚力10.5 MPa，內(nèi)摩擦角 50°，抗拉強(qiáng)度 1.0 MPa，斷裂韌度 1.8 MPa·m1/2。應(yīng)用節(jié)理巖體損傷本構(gòu)模型的試樣除采用相同的力學(xué)參數(shù)之外，還包括其內(nèi)斷續(xù)節(jié)理(僅考慮一組節(jié)理裂隙)的統(tǒng)計(jì)幾何參數(shù)：裂隙平均半徑為2 cm(或4 cm)、體積密度0.03條/cm3、裂隙面法向矢量(0.68, 0.27, 0.68)。數(shù)值試驗(yàn)得到的軸向應(yīng)變-應(yīng)力曲線如圖1所示。

圖1 單軸壓縮數(shù)值試驗(yàn)軸向應(yīng)變-應(yīng)力關(guān)系Fig.1 Axial strain versus axial stress from uniaxial compression test

圖1中，“彈塑性損傷模型-1”和“彈塑性損傷模型-2”分別對(duì)應(yīng)裂隙半徑為2 cm和4 cm時(shí)的情況。由圖1可見(jiàn)：彈塑性損傷模型的曲線斜率略大于常規(guī)模型的曲線斜率，且表現(xiàn)為隨荷載的增加而增大(荷載較小時(shí)增幅不明顯，荷載較大時(shí)較為顯著)；與常規(guī)彈塑性模型相比，損傷模型在較低的應(yīng)力下進(jìn)入塑性屈服，裂隙半徑越大(即等效柔度越大)，則相應(yīng)的屈服點(diǎn)應(yīng)力越低。數(shù)值試驗(yàn)表明：所實(shí)現(xiàn)的彈塑性損傷本構(gòu)模型能正確地描述含斷續(xù)節(jié)理巖體在初始損傷及變載作用下的損傷演化所導(dǎo)致的力學(xué)性質(zhì)的弱化現(xiàn)象。

2 支持向量機(jī)與PSO-DE優(yōu)化算法相結(jié)合的位移反分析流程

基于位移監(jiān)測(cè)信息的力學(xué)參數(shù)反分析目標(biāo)函數(shù)一般采用如下形式

式中：u和u分別為計(jì)算位移向量和監(jiān)測(cè)位移向量，由各測(cè)點(diǎn)計(jì)算(監(jiān)測(cè))位移值構(gòu)成；u由有限元(有限差分)計(jì)算求取。

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

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

3 工程實(shí)例應(yīng)用

3.1 工程概況與地質(zhì)條件

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

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

表1 廠房區(qū)節(jié)理分組及其幾何特征參數(shù)Table 1 Characteristic geometric parameters of joint sets in powerhouse area

3.2 開(kāi)挖支護(hù)及位移監(jiān)測(cè)

主廠房共分11層開(kāi)挖施工，分層高度5～8 m(見(jiàn)圖3)。為減小爆破對(duì)邊墻圍巖的擾動(dòng)，每一分層采用分塊開(kāi)挖的方式，首先開(kāi)挖中心區(qū)塊然后擴(kuò)挖上下游邊墻2 m厚的保護(hù)區(qū)塊，每層開(kāi)挖結(jié)束后立即進(jìn)行噴混凝土和系統(tǒng)錨桿支護(hù)。為監(jiān)測(cè)主廠房的圍巖變形，沿其軸向共布置有4個(gè)位移監(jiān)測(cè)斷面，隨開(kāi)挖的進(jìn)行逐步安設(shè)多點(diǎn)位移計(jì)，多點(diǎn)位移計(jì)測(cè)孔布置4個(gè)錨點(diǎn)，錨點(diǎn)間距不一，最深錨點(diǎn)距孔口24 m。2-2斷面變形監(jiān)測(cè)布置如圖3所示。由于有相當(dāng)數(shù)量的測(cè)孔量測(cè)數(shù)據(jù)不能很好地反映圍巖的變形規(guī)律，因而真正用于反分析的測(cè)孔只有4個(gè)：1-1斷面的MZC1-5；2-2斷面的MZC2-4和MZC2-5；3-3斷面的MZC3-4。提取的量測(cè)數(shù)據(jù)為測(cè)孔安裝之后至廠房開(kāi)挖完畢這一期間所發(fā)生的孔口相對(duì)孔底的位移，即增量相對(duì)位移。

3.3 數(shù)值模擬模型與位移反分析

數(shù)值模擬模型是進(jìn)行參數(shù)反分析的關(guān)鍵，在分析地勘資料的基礎(chǔ)上，建立了廠房洞群的三維數(shù)值模型。模型長(zhǎng)480 m，寬280 m，模型底部標(biāo)高EL230.0，向上延伸至地表。為實(shí)現(xiàn)后續(xù)的洞群穩(wěn)定性分析的目的，所建洞群除主廠房、主變室及母線洞之外，還包括引水隧洞和尾水隧洞。模型內(nèi)包含了洞群穿越的所有地層(P2w5，P2c，…，T1yn1)，以及2條層間錯(cuò)動(dòng)帶(fj1，fj2)，洞群與地層的空間關(guān)系如圖4(a)所示。模型計(jì)算網(wǎng)格劃分如圖4(b)所示。模型初始地應(yīng)力場(chǎng)采用側(cè)壓系數(shù)法擬合獲得，模擬開(kāi)挖步完全按實(shí)際分層開(kāi)挖情況來(lái)確定。對(duì)主廠房所處的T1y22地層應(yīng)用所開(kāi)發(fā)的節(jié)理巖體損傷本構(gòu)模型，所需的力學(xué)參數(shù)有彈性模量、泊松比、黏聚力初值和殘余值、內(nèi)摩擦角初值和殘余值、巖石斷裂韌度，此外還有裂隙的幾何參數(shù)(其值見(jiàn)表1)。其余地層則應(yīng)用莫爾-庫(kù)倫模型。

圖3 典型斷面多點(diǎn)位移計(jì)布設(shè)及主廠房開(kāi)挖分層Fig.3 Excavation stages and layout of extensometers of powerhouse cavern

圖4 三維數(shù)值計(jì)算模型Fig.4 3D numerical model

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

在參數(shù)敏感性分析基礎(chǔ)上，由均勻設(shè)計(jì)確定待反分析參數(shù)的取值方案，代入數(shù)值模擬模型計(jì)算后提取對(duì)應(yīng) 4個(gè)測(cè)孔的增量相對(duì)位移，形成“參數(shù)-增量相對(duì)位移”樣本集供訓(xùn)練 SVR使用。在正式反分析之前，首先以一組假想量測(cè)位移檢驗(yàn)了前述位移反分析方法的可靠性，結(jié)果表明：敏感性最強(qiáng)的 E_T1y22反分析值與假想真實(shí)值的相對(duì)誤差為 0.87%～3.56%，的誤差較大，總體而言，反分析方法可靠。以真實(shí)量測(cè)位移反分析參數(shù)的過(guò)程見(jiàn)圖5，由此得到的E_T1y22，GPa，1.56 MPa，0.31 MPa，23.37 GPa 和 11.8 GPa。反分析所得參數(shù)對(duì)應(yīng)的計(jì)算增量相對(duì)位移見(jiàn)表 2，與量測(cè)值相比，最大和最小相對(duì)誤差分別為 16.56%和7.54%，圖 5(a)所示為誤差平方和最終收斂至 0.098 mm2，而實(shí)際誤差平方和為0.152 mm2，二者比較接近，因此反分析成果是令人滿意的。為進(jìn)行對(duì)比，所有地層均應(yīng)用莫爾-庫(kù)倫本構(gòu)模型，對(duì)參數(shù)進(jìn)行了反分析，得到的T1y22地層等效彈性模量和黏聚力分別為25.38 GPa和1.30 MPa，二者對(duì)比可見(jiàn)：采用考慮損傷的本構(gòu)對(duì)應(yīng)的參數(shù)值明顯大于等效參數(shù)值，這是因?yàn)閾p傷模型中的彈性模量為巖石的彈性模量，該值理所當(dāng)然的應(yīng)比巖體的要大，由此可以近似認(rèn)為開(kāi)挖擾動(dòng)后圍巖的彈性模量平均降低20.13%，靠近開(kāi)挖邊界處圍巖的等效彈模降低值應(yīng)高于該值，而在遠(yuǎn)離開(kāi)挖邊界處，則應(yīng)遠(yuǎn)小于該值。

圖5 均方誤差收斂曲線及參數(shù)識(shí)別過(guò)程曲線Fig.5 Processes of SSE convergence and parameters identification

表2 測(cè)孔孔口處計(jì)算增量相對(duì)位移值與實(shí)測(cè)值對(duì)比Table 2 Incremental displacements from back analysis and field measurement

將反分析所得參數(shù)代入數(shù)值計(jì)算模型進(jìn)行思林主廠房及相鄰洞室分步開(kāi)挖正演分析，得到的圍巖內(nèi)裂隙相對(duì)擴(kuò)展長(zhǎng)度比動(dòng)態(tài)變化的模擬結(jié)果如圖6所示。

圖6 裂隙相對(duì)擴(kuò)展長(zhǎng)度比隨開(kāi)挖過(guò)程的動(dòng)態(tài)變化Fig.6 Evolution of Ex_r during excavation

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

4 結(jié)論

(1) 對(duì)開(kāi)挖于含成組展布的且相對(duì)結(jié)構(gòu)尺寸為小量的斷續(xù)節(jié)理裂隙巖體內(nèi)的地下洞室，開(kāi)挖擾動(dòng)區(qū)的形成與裂隙的擴(kuò)展演化密切相關(guān)，只有采用節(jié)理巖體損傷本構(gòu)模型才能夠較好地描述這一實(shí)際情況?；贔LAC3D內(nèi)置應(yīng)變軟化模型所開(kāi)發(fā)的節(jié)理巖體損傷本構(gòu)模型充分考慮節(jié)理巖體的各向異性特征，可將其用于考慮擾動(dòng)區(qū)的位移反分析。

(2) 采用節(jié)理巖體損傷本構(gòu)的數(shù)值計(jì)算雖然不能給出一個(gè)直觀的圍巖體內(nèi)漸變參數(shù)場(chǎng)，但是由于每個(gè)單元體的裂隙擴(kuò)展長(zhǎng)度不同，所以單元的剛度張量也不同，這相當(dāng)于間接地給出一個(gè)漸變參數(shù)場(chǎng)。

(3) 節(jié)理巖體損傷本構(gòu)模型參數(shù)較多，若所有參數(shù)都進(jìn)行反分析，將導(dǎo)致反分析解的不唯一性問(wèn)題更加突出，因此，在反分析過(guò)程中通過(guò)參數(shù)敏感性分析減少待反分析參數(shù)的個(gè)數(shù)是必要的。此外，作為替代數(shù)值計(jì)算的“參數(shù)-位移”非線性映射建模工具和正演反分析尋優(yōu)控制算法都對(duì)反分析解的可靠性有決定性影響，采用SVR和PSO-DE優(yōu)化算法進(jìn)行反分析可在一定程度上提高反分析解的可靠性。

[1] 楊林德. 巖土工程問(wèn)題的反演理論與工程實(shí)踐[M]. 北京: 科學(xué)出版社, 1996: 101-106.YANG Lin-de. Theory of back analysis in geotechnical engineering and its application[M]. Beijing: Science Press, 1996:101-106.

[2] 侯哲生, 韓文峰, 李曉. 金川二礦區(qū)巷道圍巖力學(xué)參數(shù)空間效應(yīng)數(shù)值分析[J]. 巖石力學(xué)與工程學(xué)報(bào), 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] 李寧, 段小強(qiáng), 陳方方, 等. 圍巖松動(dòng)圈的彈塑性位移反分析方法探索[J]. 巖石力學(xué)與工程學(xué)報(bào), 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] 高廣運(yùn), 馮世進(jìn), 李鴻博. 卵礫石層中隧道考慮松動(dòng)圈的位移反分析應(yīng)用[J]. 地下空間與工程學(xué)報(bào), 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ìn)兵, 鄧建輝. 高地應(yīng)力條件下大型地下廠房松動(dòng)區(qū)變化規(guī)律及參數(shù)反演[J]. 巖土力學(xué), 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] 倪紹虎, 肖明. 基于圍巖松動(dòng)圈的地下工程參數(shù)場(chǎng)位移反分析[J]. 巖石力學(xué)與工程學(xué)報(bào), 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] 江權(quán). 高地應(yīng)力硬巖彈脆塑性劣化本構(gòu)模型與大型地下洞室群圍巖穩(wěn)定性分析[D]. 武漢: 中科院武漢巖土力學(xué)研究所,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] 孫衛(wèi)軍, 周維垣. 裂隙巖體彈塑性-損傷本構(gòu)模型[J]. 巖石力學(xué)與工程學(xué)報(bào), 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] 周維垣, 楊延毅. 節(jié)理巖體的損傷斷裂模型及驗(yàn)證[J]. 巖石力學(xué)與工程學(xué)報(bào), 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] 藍(lán)航. 節(jié)理巖體采動(dòng)損傷本構(gòu)模型及其在露井聯(lián)采工程中的應(yīng)用[D]. 北京: 煤炭科學(xué)研究總院, 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] 柴紅保, 曹平, 趙延林, 等. 裂隙巖體損傷演化本構(gòu)模型的實(shí)現(xiàn)及應(yīng)用[J]. 巖土工程學(xué)報(bào), 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] 朱維申, 李術(shù)才, 陳衛(wèi)忠, 等. 節(jié)理巖體破壞機(jī)理和錨固效應(yīng)及工程應(yīng)用[M]. 北京: 科學(xué)出版社, 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] 吳群英, 林亮. 應(yīng)用數(shù)理統(tǒng)計(jì)[M]. 天津: 天津大學(xué)出版社,2004: 149-160.WU Qun-ying, LIN Liang. Applied mathematical statistics[M].Tianjin: Tianjin University Press, 2004: 149-160.

[17] 李江騰, 古德生, 曹平, 等. 巖石斷裂韌度與抗壓強(qiáng)度的相關(guān)規(guī)律[J]. 中南大學(xué)學(xué)報(bào): 自然科學(xué)版, 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.