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

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

0 引言

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

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

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

1 黏聚力模型

(a) 法向

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

1.1 黏聚力單元損傷原理

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

(1)

式中：t為單元應(yīng)變矢量；tn、ts、tt分別為單元法向和2個(gè)切向的牽引力；K為單元?jiǎng)偠染仃嚕沪艦閼?yīng)變矢量；εn、εs、εt分別為單元法向和2個(gè)切向的應(yīng)變。

矢量ε表達(dá)式如下：

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

。

(2)

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

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

(3)

1.2 損傷模式

研究表明，裂紋的擴(kuò)展是一個(gè)復(fù)合型損傷演化過(guò)程，包括法向和2個(gè)切向的擴(kuò)展。B-K準(zhǔn)則[10]被廣泛應(yīng)用于預(yù)測(cè)裂縫擴(kuò)展，其表達(dá)式見(jiàn)式(4)。該式用于當(dāng)巖體完全沿第1、2剪切方向變形時(shí)，2個(gè)方向的臨界斷裂能相同的情況。

(4)

1.3 損傷演化

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

(5)

(6)

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

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

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

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

1—4表示塊體編號(hào)。

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

2.1 模型建立

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

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

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

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

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

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

2.3 模型分析步

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

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

2.4 模擬結(jié)果

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

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

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

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

(a) 時(shí)間步為0.3

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

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

(a) 時(shí)間步為0.1

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

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

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

(a) 時(shí)間步為0.1

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

(a) 時(shí)間步為0.3

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

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

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

(a) 時(shí)間步為0.3

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

(a) 時(shí)間步為0.5

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

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

3 結(jié)論與討論

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

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

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

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

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