離散裂隙網(wǎng)絡(luò)模擬在地下硐室穩(wěn)定性中的應(yīng)用
——以甘肅北山地下實驗室為例

紀(jì)景仁， 王駒， 唐振平，2， 李亞偉， 凌輝

(1.南華大學(xué)資源環(huán)境與安全工程學(xué)院， 衡陽 421001； 2. 南華大學(xué)稀有金屬礦產(chǎn)開發(fā)與廢物地質(zhì)處置技術(shù)湖南省重點實驗室， 衡陽 421001； 3. 核工業(yè)北京地質(zhì)研究院中核高放廢物地質(zhì)處置評價技術(shù)重點實驗室， 北京 100029； 4. 國家原子能機(jī)構(gòu)高放廢物地質(zhì)處置創(chuàng)新中心， 北京 100029)

隨著中國核能事業(yè)的飛速發(fā)展，高放廢物的處理和處置，逐漸成為一個重大的安全和環(huán)保問題[1]。目前深部地質(zhì)處置是國際公認(rèn)的安全、可行的高放廢物處置方式，采用的是“多重屏障系統(tǒng)”設(shè)計思路，把高放廢物處置在數(shù)百米深的地下圍巖中。一般把廢物體、廢物處置容器和緩沖回填材料稱為“工程屏障”，把周圍的地質(zhì)體稱為“天然屏障”[2]。中國高放廢物地質(zhì)處置遵循“選址-地下實驗室-處置庫”的三步走戰(zhàn)略，現(xiàn)已進(jìn)入地下實驗室研發(fā)階段，確認(rèn)甘肅北山新場花崗巖場址為中國首座高放廢物處置地下實驗室建設(shè)場址[3]，并完成了地下實驗室的初步設(shè)計。地下實驗室已于2021年6月17日正式開工。

對于北山花崗巖這類堅硬的高放廢物處置庫圍巖，巖體中的節(jié)理裂隙是主要的導(dǎo)水通道，核素可能作為溶質(zhì)隨地下水遷移，同時其幾何形態(tài)和組合特征對地下硐室圍巖變形和失穩(wěn)起著重要的控制作用，因此研究巖體裂隙對地下硐室穩(wěn)定性的影響是高放廢物處置工程成功建設(shè)必不可少的一環(huán)。巖體節(jié)理裂隙作為一種地質(zhì)構(gòu)造運(yùn)動的產(chǎn)物，在一定的構(gòu)造應(yīng)力場中生成，所以其空間分布具備一定的統(tǒng)計規(guī)律，研究節(jié)理裂隙幾何特征和分布規(guī)律及構(gòu)建與現(xiàn)場實際圍巖具有相似統(tǒng)計特征的三維節(jié)理網(wǎng)絡(luò)模型，是進(jìn)行節(jié)理裂隙巖體工程力學(xué)行為計算與分析研究的基礎(chǔ)[4]?；谶@種理論思想，在20世紀(jì)80年代，Long等[5]正式提出離散裂隙網(wǎng)絡(luò)(discrete fracture network，DFN)模型。DFN是一種依據(jù)節(jié)理裂隙的方向、位置、密度、大小統(tǒng)計規(guī)律，隨機(jī)生成裂隙圓盤，通過展布于空間中裂隙集團(tuán)來構(gòu)建的錯綜復(fù)雜裂隙網(wǎng)絡(luò)模型，并賦予每個裂隙圓盤相應(yīng)的力學(xué)屬性，從而實現(xiàn)各種地應(yīng)力、地質(zhì)力學(xué)效果的建模和擬合。近年來，隨著計算機(jī)技術(shù)的快速發(fā)展， DFN模型在工程地質(zhì)等領(lǐng)域被廣泛應(yīng)用。張光福等[6]利用3DEC中的DFN技術(shù)，實現(xiàn)了煤巖井壁的仿真模擬，并分析其穩(wěn)定性，應(yīng)用效果較好；羅攀等[7]使用DFN技術(shù)探究了陸相頁巖裂縫的尺寸、密度和角度等儲層參數(shù)對二氧化碳濾失速度的影響；Cacciari等[8]利用地面激光掃描儀掃描隧道圍巖表面結(jié)構(gòu)面自動生成DFN模型，并估計隧道圍巖中結(jié)構(gòu)面的體密度；Wang等[9]基于表征單元體(representative elementary volume，REV)和合成巖體(synthetic rock mass，SRM)建模技術(shù)提出了一種用于節(jié)理巖體開挖響應(yīng)模擬的DFN-DEM多尺度建模方法，并驗證了其適用性。

現(xiàn)以甘肅北山新場地下實驗室場址地區(qū)地表、鉆孔節(jié)理裂隙為研究對象，運(yùn)用數(shù)理統(tǒng)計與空間幾何手段，對節(jié)理裂隙產(chǎn)狀、大小、密度進(jìn)行描述，在此基礎(chǔ)上采用三維離散單元軟件3DEC和Monte Carlo隨機(jī)分析方法，建立離散裂隙網(wǎng)絡(luò)(DFN)，進(jìn)而運(yùn)用離散裂隙網(wǎng)絡(luò)-離散元(DFN-DEM)耦合方法[10]，即用DFN模型對離散元塊體進(jìn)行切割，結(jié)合甘肅北山即將開展建設(shè)地下實驗室工程實踐，構(gòu)建能反映工程巖體裂隙分布特征的等效節(jié)理裂隙巖體計算模型，并進(jìn)行地下硐室的穩(wěn)定性分析，這也是DFN-DEM方法首次運(yùn)用到中國高放廢物地質(zhì)處置地下實驗室硐室穩(wěn)定性的研究，旨在為中國即將開展建設(shè)的高放廢物地質(zhì)處置地下實驗室深部硐室穩(wěn)定性評價提供技術(shù)支撐和方法參考。

1 模擬對象與模擬方法

根據(jù)甘肅北山地下實驗室初步建設(shè)方案，地下實驗室水平實驗主巷道最大埋深可達(dá)地下-560 m。以地下實驗室地下-560 m深度主試驗區(qū)的水平硐室為研究對象，開展地下硐室穩(wěn)定性研究，探索離散裂隙網(wǎng)絡(luò)技術(shù)在北山地下實驗室深部地下硐室穩(wěn)定性分析中的應(yīng)用效果。假設(shè)該水平巷道斷面為圓形，直徑7 m、長度40 m。

采用數(shù)理統(tǒng)計和幾何分析方法，首先對甘肅北山地下實驗室場址地表裂隙和鉆孔裂隙數(shù)據(jù)進(jìn)行規(guī)律統(tǒng)計，獲得研究區(qū)裂隙的產(chǎn)狀、尺寸和密度數(shù)據(jù)，建立離散裂隙網(wǎng)絡(luò)模型；然后構(gòu)建離散裂隙網(wǎng)絡(luò)-離散元等效巖體，獲得能反映北山地下實驗室-560 m深度開挖硐室裂隙圍巖特征的模型；最后開展等效巖體模型的數(shù)值模擬運(yùn)算，得到硐室圍巖的位移特征和應(yīng)力狀態(tài)，實現(xiàn)地下硐室穩(wěn)定性分析與評價。

2 裂隙的幾何描述及統(tǒng)計規(guī)律

選取甘肅北山新場地下實驗室場址地區(qū)，測量3個編號為ZK1、ZK2、ZK7的工程勘察鉆孔的裂隙傾角、傾向、深度數(shù)據(jù)，鉆孔位置如圖1所示；對新場場址勘察孔周邊地表露頭測量得到的4 585條裂隙產(chǎn)狀、跡長信息進(jìn)行幾何特征分析。對上述裂隙統(tǒng)計規(guī)律進(jìn)行總結(jié)，獲取相關(guān)參數(shù)數(shù)據(jù)，為下一步離散裂隙網(wǎng)絡(luò)的生成提供數(shù)據(jù)基礎(chǔ)。

由圖2可知，甘肅北山新場場址勘察孔周邊地表露頭裂隙可分為4個優(yōu)勢產(chǎn)狀組，采用K-means聚類分析法，利用SPSS軟件，分析得到各組優(yōu)勢產(chǎn)狀，并劃分每條裂隙的優(yōu)勢產(chǎn)狀組歸屬，得到勘察孔周邊地表露頭結(jié)構(gòu)面優(yōu)勢產(chǎn)狀組。然后，根據(jù)得到的各組優(yōu)勢結(jié)構(gòu)面，分別對傾向、傾角劃分區(qū)間進(jìn)行概率分布研究，發(fā)現(xiàn)各優(yōu)勢組傾向、傾角符合正態(tài)分布特征。接著，選用地下實驗室場址地表裂隙條數(shù)大于30條，且分布較為均勻的露頭，采用楊春和等[11]提出的圓形窗口法估算裂隙跡長，由頻數(shù)分布發(fā)現(xiàn)各優(yōu)勢組裂隙跡長呈對數(shù)正態(tài)分布。有研究表明裂隙圓盤直徑的分布形式與跡長分布形式相同[12]，因此，圓盤直徑的概率分布形式也確定為對數(shù)正態(tài)分布。而后，采用伍法權(quán)法[13]估算裂隙圓盤直徑分布函數(shù)的均值與標(biāo)準(zhǔn)差。

圖1 鉆孔位置圖Fig.1 Drilling location map

圖2 地表裂隙極點等密度圖Fig.2 Pole and contour plots diagram of surface fractures

最后，因模擬研究對象位于地下560 m深度的空間，故利用3個編號為ZK1、ZK2、ZK7的勘察孔分別在地下500～598、500～568、500～548 m深度使用鉆孔電視裝置測得的共計87條孔壁裂隙數(shù)據(jù)，劃分每條鉆孔裂隙的優(yōu)勢組歸屬，并求得各優(yōu)勢組線密度P10，即單位長度內(nèi)的裂隙數(shù)目。離散裂隙網(wǎng)絡(luò)在達(dá)到鉆孔定義的目標(biāo)裂隙線密度P10時停止。

綜上，通過對甘肅北山新場巖體裂隙特征參數(shù)分析，獲得了裂隙空間各參數(shù)，如表1所示。

3 數(shù)值模型建立

3.1 DFN生成與等效巖體構(gòu)建

確定建模所需裂隙參數(shù)概率模型之后，運(yùn)用3DEC數(shù)值模擬軟件中的離散裂隙網(wǎng)絡(luò)(DFN)功能[14]，基于Monte Carlo隨機(jī)抽樣原則和線性同余法理論，生成能反映研究區(qū)節(jié)理裂隙分布特征的離散裂隙網(wǎng)絡(luò)。根據(jù)裂隙尺寸大小，刪除尺寸較小裂隙，將DFN模型簡化到一個擁有許多原始模型特征的簡單版本，以提高后續(xù)力學(xué)模型計算效率。本文為了模型計算需要，假設(shè)尺寸大于7 m的大尺寸裂隙對硐室圍巖穩(wěn)定性起控制作用。

簡化后的離散裂隙網(wǎng)絡(luò)模型如圖3所示，其空間大小為40 m×40 m×40 m，以正方形形心為坐標(biāo)軸原點。

圖3 DFN模型Fig.3 DFN model

表1 甘肅北山新場離散裂隙網(wǎng)絡(luò)模擬參數(shù)表Table 1 Simulation parameter table of discrete fracture network in Beishan Xinchang， Gansu Province

為了進(jìn)一步模擬實現(xiàn)甘肅北山地下實驗室地下-560 m深度，直徑為7 m的圓形隧道及其所處的節(jié)理裂隙圍巖條件，如圖4所示，建立一個邊長與DFN模型相同，為40 m×40 m×40 m，以正方形形心為坐標(biāo)原點的塊體，并規(guī)定其初始地應(yīng)力最小水平主應(yīng)力方向為X軸方向，最大水平主應(yīng)力方向為Y軸方向，模型高度方向為Z軸方向。同時，在該塊體中開挖一條長40 m，軸向沿Y軸，圓心為坐標(biāo)軸原點，直徑為7 m的水平圓形隧道，然后用離散裂隙網(wǎng)絡(luò)(DFN)模型切割離散元(DEM)塊體，構(gòu)建與實際巖體有相似裂隙分布特征的等效巖體，最終實現(xiàn)對研究區(qū)目標(biāo)巖體的數(shù)值模型構(gòu)建。

圖4 DFN-DEM等效巖體Fig.4 DFN-DEM equivalent rock mass

3.2 模型參數(shù)賦值

3.2.1 初始地應(yīng)力賦值

根據(jù)對甘肅北山新場地下實驗室預(yù)選區(qū)鉆孔地應(yīng)力的實測數(shù)據(jù)分析發(fā)現(xiàn)，研究區(qū)地應(yīng)力的關(guān)系式[15]為

(1)

式(1)中：σH、σh、σv分別為最大主應(yīng)力、最小主應(yīng)力和豎直主應(yīng)力，MPa；H為深度，m。

研究假設(shè)DFN-DEM等效巖體模型中的圓形隧道底板位于地下-560 m深度，以模擬地下實驗室-560 m深度開挖硐室所處的圍巖環(huán)境，該模型地應(yīng)力條件按式(1)賦值。同時，模型的邊界條件為：上邊界為自由邊界，下邊界約束豎向(Z向)位移，左右邊界為側(cè)向約束。

3.2.2 巖石與裂隙物理力學(xué)參數(shù)

在DFN中的裂隙為直徑大小有限的圓盤狀，但是在3DEC中，塊體不能被部分切割，因此在3DEC中所有的裂隙必然都比DFN模型中相應(yīng)的裂隙要大。這個問題可以通過在環(huán)形裂隙的內(nèi)部和外部分別賦予不同的裂隙屬性來解決。如圖5所示，圓盤內(nèi)部為真實裂隙，賦予真實裂隙物理屬性，圓盤外部為虛擬裂隙，賦予與完整巖塊力學(xué)屬性近似一致的虛擬裂隙物理屬性[16]。

在DFN-DEM等效巖體數(shù)值模擬計算中，巖石塊體采用摩爾庫倫本構(gòu)模型，節(jié)理裂隙采用庫倫滑移模型，北山新場花崗巖體物理力學(xué)參數(shù)與裂隙物理參數(shù)取值如表2所示[17-18]。

圖5 裂隙圓盤賦值屬性示意圖Fig.5 Schematic diagram of fissure disc assignment attribute

4 模擬計算結(jié)果

4.1 位移場分析

經(jīng)過3DEC迭代運(yùn)算，甘肅北山地下實驗室圓形開挖巷道的整體圍巖變形運(yùn)動情況如圖6(a)所示。硐室圍巖的整體變形趨勢表現(xiàn)為向臨空面方向移動，頂拱下沉，底板回彈，且頂拱向下的變形量大于底板向上的回彈量；離硐室表面越遠(yuǎn)變形量越小，硐室圍巖變形幾乎關(guān)于硐口左右對稱。硐室頂拱的最大下沉量為5.0 mm，底板的最大的回彈量為3.7 mm，左右邊墻變形為1.0～3.5 mm。

表2 巖體及裂隙物理力學(xué)參數(shù)Table 2 Physical and mechanical parameters of rock mass and fractures

圖6(b)為硐室以軸線為剖面的豎直方向(Z方向)位移云圖，變形經(jīng)300倍放大，可以看出，硐室在與貫穿硐室圍巖的裂隙相交處發(fā)生了輕微錯動，頂拱下沉量與底板回彈量較完整圍巖增大，未貫穿硐室圍巖的裂隙對硐室圍巖的變形量影響較小，且并未擴(kuò)展至貫穿硐室圍巖。

圖7為硐室圍巖變形特征點豎向(Z方向)位移數(shù)據(jù)監(jiān)測曲線，從總體趨勢來看，頂拱位移量大于底板位移量，頂拱位移隨計算時步逐漸變大，頂拱裂隙處圍巖豎向變形可達(dá)4.3 mm；底板位移在前期逐漸增加，但是后期由于底板應(yīng)力重新分布，底板位移變形表現(xiàn)出一定的回彈，底板裂隙處圍巖豎向變形值為2.9 mm；同時可以得知，無論頂拱或底板位置，裂隙與硐室相交的薄弱處位移量大于完整圍巖處位移量，這也與圖6(b)位移云圖表征結(jié)果相符。

圖6 位移云圖Fig.6 Displacement cloud map

圖7 硐室圍巖變形特征監(jiān)測點豎向位移曲線Fig.7 Vertical displacement curve of the monitoring point for the deformation characteristics of the surrounding rock of the chamber

4.2 應(yīng)力場分析

硐室開挖引起圍巖應(yīng)力重分布，硐壁圍巖因開挖卸荷而產(chǎn)生位移變形，應(yīng)力狀態(tài)發(fā)生較大改變[19]。

圖8(a)為硐室圍巖的最大主應(yīng)力云圖，可以看出，硐室周邊圍巖產(chǎn)生較大的切向壓應(yīng)力集中。最大應(yīng)力集中的位置出現(xiàn)在硐室的兩側(cè)邊墻裂隙與硐壁圍巖交錯處，最大主應(yīng)力值可達(dá)48.2 MPa，頂拱、底板為應(yīng)力集中區(qū)域，但是應(yīng)力集中程度與范圍都小于邊墻，其最大主應(yīng)力應(yīng)力為13.0～30.0 MPa；硐室附近應(yīng)力變化顯著，但距硐室5 m之外的區(qū)域，應(yīng)力分布較為均勻。

圖8(b)為硐室沿軸線剖面最大主應(yīng)力云圖，可以得知，頂拱與底板發(fā)生應(yīng)力集中現(xiàn)象，表現(xiàn)為越靠近開挖面應(yīng)力集中程度越大，且硐室與裂隙相交處應(yīng)力發(fā)生突變，應(yīng)力集中現(xiàn)象最為明顯最大主應(yīng)力量值達(dá)30.2 MPa。

5 結(jié)論

(1)通過分析統(tǒng)計研究區(qū)地表裂隙與深部鉆孔裂隙分布規(guī)律，基于3DEC離散元仿真軟件，建立離散裂隙網(wǎng)絡(luò)，構(gòu)建能反映深部地下空間巖體裂隙分布特征的等效巖體模型，進(jìn)行數(shù)值模擬運(yùn)算，為北山地下實驗室深部地下硐室的穩(wěn)定性分析提供一種簡單快速的技術(shù)手段。

圖8 最大主應(yīng)力云圖Fig.8 Maximum principal stress cloud diagram

(2)通過對甘肅北山地下實驗室場址地表裂隙進(jìn)行數(shù)理統(tǒng)計與幾何特征分析，發(fā)現(xiàn)研究區(qū)節(jié)理裂隙可劃分為4個優(yōu)勢組，各組裂隙傾向、傾角均符合對數(shù)正態(tài)分布，裂隙圓盤直徑符合對數(shù)正態(tài)分布規(guī)律；對場址地區(qū)地下500～600 m深度3個勘察孔裂隙數(shù)據(jù)分析可知，地下實驗室-560 m深度主試驗區(qū)附近圍巖裂隙各優(yōu)勢組線密度為0.222 2～0.027 7 條/m。

(3)在現(xiàn)有分析條件下，基于等效巖體模擬結(jié)果表明北山地下實驗室地下硐室開挖后，圍巖位移較小，最大變形量僅為5.0 mm，在地下實驗室圍巖彈性應(yīng)變范圍內(nèi)，圍巖硐室周邊圍巖應(yīng)力集中程度不高，最大主應(yīng)力為48.2 MPa，遠(yuǎn)小于北山花崗巖的單軸抗壓強(qiáng)度141.1 MPa，說明硐室開挖后穩(wěn)定性較好。