離散裂隙網絡模擬在地下硐室穩定性中的應用
——以甘肅北山地下實驗室為例

紀景仁， 王駒， 唐振平，2， 李亞偉， 凌輝

(1.南華大學資源環境與安全工程學院， 衡陽 421001； 2. 南華大學稀有金屬礦產開發與廢物地質處置技術湖南省重點實驗室， 衡陽 421001； 3. 核工業北京地質研究院中核高放廢物地質處置評價技術重點實驗室， 北京 100029； 4. 國家原子能機構高放廢物地質處置創新中心， 北京 100029)

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

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

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

1 模擬對象與模擬方法

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

采用數理統計和幾何分析方法，首先對甘肅北山地下實驗室場址地表裂隙和鉆孔裂隙數據進行規律統計，獲得研究區裂隙的產狀、尺寸和密度數據，建立離散裂隙網絡模型；然后構建離散裂隙網絡-離散元等效巖體，獲得能反映北山地下實驗室-560 m深度開挖硐室裂隙圍巖特征的模型；最后開展等效巖體模型的數值模擬運算，得到硐室圍巖的位移特征和應力狀態，實現地下硐室穩定性分析與評價。

2 裂隙的幾何描述及統計規律

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

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

圖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條孔壁裂隙數據，劃分每條鉆孔裂隙的優勢組歸屬，并求得各優勢組線密度P10，即單位長度內的裂隙數目。離散裂隙網絡在達到鉆孔定義的目標裂隙線密度P10時停止。

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

3 數值模型建立

3.1 DFN生成與等效巖體構建

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

簡化后的離散裂隙網絡模型如圖3所示，其空間大小為40 m×40 m×40 m，以正方形形心為坐標軸原點。

圖3 DFN模型Fig.3 DFN model

表1 甘肅北山新場離散裂隙網絡模擬參數表Table 1 Simulation parameter table of discrete fracture network in Beishan Xinchang， Gansu Province

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

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

3.2 模型參數賦值

3.2.1 初始地應力賦值

根據對甘肅北山新場地下實驗室預選區鉆孔地應力的實測數據分析發現，研究區地應力的關系式[15]為

(1)

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

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

3.2.2 巖石與裂隙物理力學參數

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

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

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

4 模擬計算結果

4.1 位移場分析

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

表2 巖體及裂隙物理力學參數Table 2 Physical and mechanical parameters of rock mass and fractures

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

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

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

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

4.2 應力場分析

硐室開挖引起圍巖應力重分布，硐壁圍巖因開挖卸荷而產生位移變形，應力狀態發生較大改變[19]。

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

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

5 結論

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

圖8 最大主應力云圖Fig.8 Maximum principal stress cloud diagram

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

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