駱 健， 章楊松，*， 李曉昭

(1. 南京理工大學理學院， 江蘇 南京 210094; 2. 南京大學地球科學與工程學院， 江蘇 南京 210093)

0 引言

巖體在形成和演化的過程中，由于受到各種地質(zhì)作用和構(gòu)造力的影響，會形成大小不一、形態(tài)復雜、分布錯綜復雜的巖體結(jié)構(gòu)面，并與隧道開挖形成的臨空面組成不同規(guī)模的巖石塊體，這些塊體的失穩(wěn)或垮落會破壞巖體整體的穩(wěn)定性。因此,在設(shè)計階段，根據(jù)形成的結(jié)構(gòu)面信息分析塊體的穩(wěn)定性，對優(yōu)化設(shè)計和開挖支護方案具有重要的指導意義。

由石根華[1]提出的關(guān)鍵塊體理論和P.M.Warburton[2]提出的塊體可移動性的矢量法是巖石塊體穩(wěn)定性分析方法的重大突破，為用數(shù)值方法研究裂隙巖體穩(wěn)定奠定了重要基礎(chǔ)。研究裂隙巖體穩(wěn)定的方法主要有連續(xù)介質(zhì)法、非連續(xù)介質(zhì)法和連續(xù)-非連續(xù)介質(zhì)的混合法等[3]。黃達等[4]運用數(shù)值計算軟件FLAC3D，通過連續(xù)介質(zhì)法考慮塊體的變形以及塑性區(qū)分布特征來評價塊體的穩(wěn)定性。張雨霆等[5-6]提出復雜巖石塊體識別的單元重構(gòu)-聚合法，并進一步采用連續(xù)介質(zhì)法，運用數(shù)值分析軟件FLAC3D中的Interface界面單元模擬巖石塊體穩(wěn)定性。Wu等[7]提供了一種DDA(非連續(xù)變形分析)的離散化方法來計算塊體穩(wěn)定性，此方法能夠模擬出塊體的掉落行為。

連續(xù)介質(zhì)法的缺點是不能模擬塊體的掉落行為；在模擬具有膠結(jié)結(jié)構(gòu)面的塊體失穩(wěn)時，目前常用的非連續(xù)介質(zhì)法無法分析其膠結(jié)面的漸進損傷過程及狀況。針對以上問題，本文采用顯式動力有限元模擬裂隙巖體的連續(xù)-非連續(xù)變形以及含膠結(jié)的結(jié)構(gòu)面組成塊體的局部失穩(wěn)過程。選取甘肅北山某典型花崗巖研究區(qū)域測量的結(jié)構(gòu)面數(shù)據(jù)以及搜索出的多種塊體，運用有限元軟件ABAQUS中的黏聚力模型(cohesive zone model，CZM)模擬具有填充和膠結(jié)的結(jié)構(gòu)面的抗拉強度和黏聚強度。根據(jù)測量的結(jié)構(gòu)面數(shù)據(jù)信息，分4種方案對膠結(jié)結(jié)構(gòu)面的抗拉強度和黏聚強度進行賦值；運用Matlab程序?qū)λ淼蓝词夷Ｐ偷膇np文件進行修改，從而可以對洞室開挖邊界施加相應的節(jié)點荷載，通過對反加于圍巖邊界荷載的逐漸釋放來分析初始地應力和開挖誘導的二次地應力對塊體穩(wěn)定的影響。結(jié)合以上幾點，綜合考慮初始地應力、開挖誘導的二次地應力、結(jié)構(gòu)面的黏聚強度和抗拉強度，以及它們之間的摩擦應力等因素來研究塊體的穩(wěn)定性。

1 黏聚力模型

(a) 法向

(b) 切向圖1 雙線性張力-位移關(guān)系Fig. 1 Bi-linear tension-displacement relationship

1.1 黏聚力單元損傷原理

在黏聚力單元損傷之前，認為其應力應變關(guān)系符合線彈性準則，如式(1)所示。2個相鄰單元之間用類似膠材料的黏聚力單元連接，該黏聚力單元很薄，有時考慮到實際的工程問題，黏聚力單元厚度常常設(shè)置成0。

(1)

式中：t為單元應變矢量；tn、ts、tt分別為單元法向和2個切向的牽引力；K為單元剛度矩陣；ε為應變矢量；εn、εs、εt分別為單元法向和2個切向的應變。

矢量ε表達式如下：

εn=dn/To,εs=ds/To,εt=dt/To

。

(2)

式中：dn、ds、dt分別為單元法向和2個切向變形黏聚力單元分開的距離；To為黏聚力單元最初的厚度。

利用名義應力比的二次交互函數(shù)準則來確定膠結(jié)結(jié)構(gòu)面是否會斷裂，當臨界應力比平方和等于1時，認為損傷完全。此準則用公式表達如下[9]：

(3)

1.2 損傷模式

研究表明，裂紋的擴展是一個復合型損傷演化過程，包括法向和2個切向的擴展。B-K準則[10]被廣泛應用于預測裂縫擴展，其表達式見式(4)。該式用于當巖體完全沿第1、2剪切方向變形時，2個方向的臨界斷裂能相同的情況。

(4)

1.3 損傷演化

模擬材料的斷裂損傷演化通常采用剛度退化準則，該準則表達式如下[11]：

(5)

(6)

零厚度的黏聚力單元模型如圖2所示。圖中的中間層表示具有零厚度的黏聚力單元層，當塊體受到外力作用時，中間的黏聚力層會逐漸損傷；當其應力或應變滿足定義的初始臨界損傷準則時，則開始退化，黏聚力單元逐漸刪除，單元間的黏聚力消失，而上下面之間的摩擦力從分析步開始只要面面接觸就會一直起作用。根據(jù)摩爾-庫侖理論，黏聚力單元層的強度τ=c+σtanφ，其中c為結(jié)構(gòu)面的黏聚力強度，σtanφ為結(jié)構(gòu)面之間的摩擦強度。

圖2 零厚度的黏聚力單元模型Fig. 2 Cohesive element model of zero thickness

2 圍巖塊體的穩(wěn)定性分析

在甘肅北山某花崗巖典型區(qū)域研究范圍內(nèi)進行塊體的搜索，得到的塊體如圖3所示。其中，1、2、4號塊體是四面體，3號塊體是五面體，且由幾何位置初步判斷，1、2、3號塊體可能出現(xiàn)失穩(wěn)情況，4號塊體是穩(wěn)定塊體。

1—4表示塊體編號。

圖3搜索出的塊體示意圖
Fig. 3 Sketch of searched blocks

2.1 模型建立

根據(jù)得到的塊體結(jié)構(gòu)面數(shù)據(jù)，以花崗巖典型區(qū)域地下開挖隧道為例，對開挖范圍內(nèi)搜索出的塊體的穩(wěn)定性進行研究。開挖隧道斷面由邊墻與拱頂組成，可以分別定義為直線段和圓弧段?？紤]到模擬隧道開挖僅對一定范圍內(nèi)的巖體產(chǎn)生影響，根據(jù)工程經(jīng)驗，模型邊界距開挖洞口邊的距離超過5倍洞口直徑時，開挖產(chǎn)生的影響就會變得很小。在本次模擬中，考慮到洞口直徑為5 m左右，模型尺寸可取為40 m×60 m×60 m (長×寬×高)。根據(jù)北山巖體條件及核廢料地下處置的基本要求，計算模型上表面位于地下600 m埋深處，同時上表面的豎向應力按實測應力取值； 下表面位于660 m埋深處。模型的邊界條件為： 上邊界為自由邊界，下邊界約束豎向(z向)位移，左右邊界為側(cè)向約束。單元網(wǎng)格的尺寸從模型邊界到孔洞壁附近逐漸被細化，如圖4所示。

圖4 隧道模型的網(wǎng)格劃分Fig. 4 Meshing of tunnel model

2.2 材料及結(jié)構(gòu)面參數(shù)

模型頂部為應力邊界條件，應力值為上覆巖石的質(zhì)量，依據(jù)趙宏剛等[12]得出的巖石埋深與地應力的關(guān)系，取垂直主應力為16.2 MPa，x方向的側(cè)壓力系數(shù)為0.92，y方向的側(cè)壓力系數(shù)為0.65。根據(jù)北山某典型花崗巖采石場地區(qū)的統(tǒng)計資料，分別取表1所示的巖體物理力學參數(shù)，按照等效連續(xù)介質(zhì)進行計算。

表1 巖體的物理力學參數(shù)Table 1 Physico-mechanical parameters of rock mass

結(jié)構(gòu)面的強度包含抗剪強度和抗拉強度2部分，結(jié)合搜索出塊體的結(jié)構(gòu)面特征以及相應的資料[13]，參照GB/T 50218—2014《工程巖體分級標準》[14]，得出其抗拉強度在0～0.04 MPa，相應的黏聚強度在0～0.12 MPa，內(nèi)摩擦角在19°～29°。

2.3 模型分析步

本次模擬分為2個分析步，每個分析步與時間步相對應。第1個分析步是初始地應力平衡，時間步是0～0.1； 第2個分析步是模擬隧道開挖，時間步是0.1～1.4。為了分析開挖誘導的二次地應力的影響，把初始平衡步得到的開挖邊界應力加到隧道邊界上，如圖5所示。隨著時間步的增加，開挖邊界的節(jié)點力逐漸釋放，當時間步增加到1時，節(jié)點力釋放完成。釋放過程中，由于中間巖體被挖出，原有的地應力平衡被打破，將由此產(chǎn)生應力二次重分布； 同時，原有的結(jié)構(gòu)面與新產(chǎn)生的臨空面形成關(guān)鍵塊體，在地應力以及自身重力的作用下，此塊體可能會出現(xiàn)失穩(wěn)。

圖5 隧道邊界節(jié)點力Fig. 5 Joint force of tunnel boundary

2.4 模擬結(jié)果

考慮到不同位置的塊體可能的失穩(wěn)情況，分4種方案進行模擬,具體參數(shù)見表2。整個方案中，4號塊體始終保持穩(wěn)定。為了僅考慮摩擦力的影響，方案1把抗拉強度和黏聚強度賦予一個很小的數(shù)值(1×10-9MPa)，相應的摩擦因數(shù)見表2； 方案2、3、4分別考慮不同位置塊體可能的失穩(wěn)情況，通過調(diào)整不同塊體結(jié)構(gòu)面的抗拉強度和黏聚強度，確定影響塊體穩(wěn)定的主要因素。

表2 塊體結(jié)構(gòu)面參數(shù)的4種方案Table 2 Four schemes of block structural surface parameters

1)方案1模型不同時間步的位移云圖如圖6所示。在抗拉強度和黏聚強度接近于0，僅摩擦應力起作用的情況下，1、2、3號塊體幾乎同時掉落。圖7示出開挖時間步為0.1～1.4時的摩擦應力變化情況。由圖7可知，摩擦應力在起始階段逐漸增大，達到峰值后逐漸降低并趨于0。對應圖6來看，趨于0即塊體已經(jīng)完全脫離巖體結(jié)構(gòu)面并向下滑落，頂部3號塊體大約在時間步為0.5時最先脫離，另外2個塊體大約在時間步為0.6時脫離。2號塊體的摩擦應力始終相對較小，因此若想讓2號塊體保持穩(wěn)定，提高塊體結(jié)構(gòu)面的抗拉強度是最有效的辦法。對于3號塊體，由于其位置的特殊性而最容易失穩(wěn)，在接下來的探討中，發(fā)現(xiàn)即使賦予該結(jié)構(gòu)面參數(shù)區(qū)間內(nèi)最大的抗拉強度、黏聚強度和摩擦因數(shù)，該塊體還是無法保持穩(wěn)定。

2)方案2模型的位移云圖如圖8所示。由圖8可知： 開挖時間步為0.1時，3號塊體已經(jīng)失穩(wěn)，而2號塊體此時還處于穩(wěn)定狀態(tài)； 隨著開挖進行，二次應力逐漸釋放，在時間步為1時，應力釋放完成，2號塊體已經(jīng)失穩(wěn)，1號塊體始終保持穩(wěn)定狀態(tài)，塊體的位移曲線如圖9所示。下面主要探討3號塊體結(jié)構(gòu)面的損傷情況，黏聚力單元的損傷情況如圖10所示。圖10(a)示出時間步為0.1時，3號塊體的結(jié)構(gòu)面接近損傷完全(時間步為1時損傷破壞)，黏聚力單元的底邊邊角損傷最為嚴重； 時間步為0.3時，底邊邊角和兩側(cè)面部分已經(jīng)完全損傷，并失效刪除； 時間步為0.4時，已經(jīng)完全失效刪除。

(a) 時間步為0.3

(b) 時間步為1圖6 方案1不同時間步的位移云圖(單位： m)Fig. 6 Displacement nephograms of block with different time steps in scheme 1 (unit: m)

圖7 方案1不同塊體的摩擦應力的時間步曲線Fig. 7 Curves of relationship between time steps and friction stresses of different blocks in scheme 1

(a) 時間步為0.1

(b) 時間步為1圖8 方案2不同時間步的位移云圖(單位： m)Fig. 8 Displacement nephograms of block with different time steps in scheme 2 (unit: m)

圖9 方案2塊體位移隨時間步的變化情況Fig. 9 Curves of relationship between time steps and friction stresses of different blocks in scheme 2

3)方案3模型的位移云圖如圖11所示。由圖11可知: 降低1號塊體結(jié)構(gòu)面的摩擦因數(shù)而黏聚力不變，塊體在時間步為1時出現(xiàn)失穩(wěn)并逐漸滑落，此時由于2號塊體結(jié)構(gòu)面的黏聚力增加到0.1 MPa，摩擦因數(shù)為0.45，整個分析步過程中，2號塊體都保持穩(wěn)定； 3號塊體此時已被賦予結(jié)構(gòu)面參數(shù)區(qū)間的最大抗拉強度、黏聚強度和摩擦因數(shù)，但依然失穩(wěn)。方案3塊體位移曲線如圖12所示，1號塊體結(jié)構(gòu)面的損傷情況如圖13所示。

(a) 時間步為0.1

(b) 時間步為0.3圖10 3號塊體結(jié)構(gòu)面隨時間步變化的損傷情況Fig. 10 Damage of structural surfaces of block 3 variating with time steps

(a) 時間步為0.3

(b) 時間步為1圖11 方案3不同時間步的位移云圖(單位： m)Fig. 11 Displacement nephograms of block with different time steps in scheme 3 (unit: m)

4)方案4位移云圖如14所示。把2號塊體結(jié)構(gòu)面的摩擦因數(shù)降到最小值0.34，此時黏聚力沒有變，2號塊體在時間步為0.8時出現(xiàn)失穩(wěn)，到時間步為1時出現(xiàn)明顯的滑落情況。值得注意的是，計算過程中1號塊體結(jié)構(gòu)面的摩擦因數(shù)由方案3的0.4提高到0.42，而此時的1號塊體卻比方案3滑落的早。由此可見，塊體間的滑落是相互影響的，一個塊體的滑落勢必會影響其他塊體。方案4塊體的位移曲線如圖15所示。

圖12 方案3塊體位移隨時間步的變化情況Fig. 12 Curves of relationship between time steps and friction stresses of different blocks in scheme 3

(a) 時間步為0.3

(b) 時間步為0.9圖13 1號塊體結(jié)構(gòu)面隨時間步變化的損傷情況Fig. 13 Damage of structural surfaces of block 1 variating with time steps

(a) 時間步為0.5

(b) 時間步為1圖14 方案4不同時間步的位移云圖(單位： m)Fig. 14 Displacement nephograms of block with different time steps in scheme 4 (unit: m)

圖15 方案4塊體位移隨時間步的變化情況Fig. 15 Curves of relationship between time steps and friction stresses of different blocks in scheme 4

3 結(jié)論與討論

本文運用黏聚力模型，結(jié)合甘肅北山某花崗巖典型區(qū)域搜索出的塊體，通過考慮具有填充和膠結(jié)的結(jié)構(gòu)面的抗拉強度、黏聚強度以及摩擦應力等因素，分析不同位置塊體的穩(wěn)定情況。得到主要結(jié)論如下。

1)采用黏聚力單元不僅能模擬僅有摩擦應力作用的結(jié)構(gòu)面情況，還能較好地模擬具有填充和膠結(jié)的結(jié)構(gòu)面的抗拉強度和黏聚強度，以及結(jié)構(gòu)面的滑移和脫開行為，可作為塊體穩(wěn)定分析計算的工具。

2)實際隧道開挖過程中，頂部的塊體最容易失穩(wěn)，因此在為隧道設(shè)計襯砌時應適當加強襯砌頂部的強度。此外，抗拉強度、黏聚強度和摩擦應力共同構(gòu)成了結(jié)構(gòu)面的強度，不同位置的塊體滑落是相互影響的，一個塊體的滑落勢必會影響其他塊體的穩(wěn)定性。

3)采用黏聚力單元模擬膠結(jié)結(jié)構(gòu)面的抗拉強度和黏聚強度，通過把初始平衡步得到的開挖邊界應力反加到隧道邊界上，并在開挖步逐漸釋放邊界應力，從而可以表征出黏聚力單元隨二次地應力釋放而逐漸損傷的過程，這為模擬結(jié)構(gòu)面損傷隨時間的漸進性發(fā)展提供了一個新思路。

本文的研究尚有許多不足之處，使用的黏聚力損傷模型為常用的雙線性函數(shù)，與實際膠結(jié)結(jié)構(gòu)面隨外力具體損傷的契合度也有待改進。此外，模擬含有巖橋的塊體穩(wěn)定也是以后研究的重點與難點。