遍布節理模型在層狀巖體模擬中的適用性研究

黃娟 ，周世杰 ，賈朝軍 ?，宋銀濤 ，張建

（1.中南大學 土木工程學院，湖南 長沙 410075；2.中鐵五局集團有限公司，湖南 衡陽 420002）

層狀巖體是巖土與地下工程建設中經常遇到的一類巖體，在自然界中廣泛分布.長期地質構造作用下所形成的層理面使巖體在強度和變形等方面都表現出明顯的各向異性，這對隧道開挖時巖體錨固［1］、襯砌開裂、仰拱隆起［2］等工程問題有著顯著的影響.因此，層狀巖體的力學行為與響應機制研究具有重要意義與研究價值.

近年來，隨著材料本構不斷完善以及計算機技術的更新迭代，越來越多的數值模擬技術用于巖體力學特性的研究，為室內試驗或現場測試的局限性提供了補充和解決辦法.王培濤等［3］應用顆粒流軟件PFC2D研究了不同層理角度的黑云變粒巖的強度特性.Singh 等［4］通過UDEC 探究了節理巖體在單軸加載條件下產生高側向應變的原因.劉愛華等［5］采用有限元軟件ANSYS 模擬了不同層面傾角的巖體抗拉、抗壓試驗下的破壞形態.此外一部分學者還將有限元法［6-7］、有限差分法［8］、離散元法［9］、有限-離散元法［10］、真實破裂過程分析方法［11］等數值方法應用于模擬層狀圍巖地下洞室的變形和破壞機理等方面.

雖然數值模擬方法繁多，但相比之下，采用離散元法能夠最有效地描述層狀巖體等不連續材料的力學性能［12］.然而，考慮到離散元法計算的效率，若要模擬全部的節理或層理構造以進行某些大型地下工程的開挖掘進是不太可取的［13］.近年來，有學者研發了高效顆粒離散元軟件MatDEM［14］，但顆粒離散元軟件很大程度上依賴于本構參數的準確標定，且該軟件暫未廣泛應用于層狀巖體模擬之中.為了避免這些限制，通常可以采用FLAC3D 中的遍布節理模型來表示一些層狀各向異性巖體.例如，蔣青青等［15］采用FLAC3D 內置的Ubiquitous-Joint 模型分析了層狀巖質邊坡開挖過程中層理傾角和傾向與安全系數之間的關系.朱澤奇等［16］、周鵬發等［17］采用改進的Ubiquitous-Joint 模型模擬了層狀圍巖地下洞室開挖時的變形和破壞.Sainsbury 等［18］針對巖體中普遍存在的隨機節理，通過建立與主節理或層理方向正交的節理集，并提出遍布節理模型參數修正準則，較好地描述了自然界中各向異性巖體強度和變形特性.

然而，目前遍布節理模型中參數的物理意義不夠明確，不能由試驗結果直接獲取.由于遍布節理模型沒有考慮節理間距和節理剛度，如果直接將其材料參數與離散元模型的巖塊和結構面參數賦值一致，模擬結果不能真實地反映實際工程中的巖體強度或變形.因此，需要對遍布節理模型的參數進行校準修正，以便為工程設計或施工提供有意義的參考.

本文通過總結部分學者對層狀各向異性巖體的研究，在Sainsbury 研究的基礎上，分別采用FLAC3D中的Ubiquitous-Joint 模型和Softening-ubiquitous 模型（考慮應變軟化的Ubiquitous-Joint 模型，以下簡稱Subiquitous 模型）以及塊體離散元軟件3DEC 對層狀巖體的力學特征進行模擬并作對比分析.基于校準后的Subiquitous 模型，通過分析新華山隧道開挖和支護過程，揭示層理對圍巖變形和破壞特征的影響，驗證遍布節理模型的適用性.

1 （應變軟化）遍布節理模型

在FLAC3D中，遍布節理模型有Ubiquitous-Joint和Subiquitous 模型兩種.Ubiquitous-Joint 模型［19］對應于摩爾-庫侖模型，即在摩爾-庫侖體中加入節理面，該節理面也服從摩爾-庫侖屈服準則，使材料表現出強度各向異性.Ubiquitous-Joint 模型同時考慮了巖石基質和節理的物理力學屬性，必須在模型的指定區域內同時賦予基質和節理的參數.

節理面的破壞包括拉伸和剪切破壞，如圖1 所示，其中剪切破壞包絡線AB表示為fs=0：

圖1 節理面破壞準則Fig.1 Joint failure criterion

拉伸破壞包絡線BC表示為ft=0：

式中：?j、cj和分別為節理面的內摩擦角、黏聚力和抗拉強度；σ3′3′為節理面上的正應力.

該模型的計算公式中未涉及節理的間距、長度以及巖層的彎曲剛度等.如果不對相應的力學參數進行校準，可能會得到錯誤的巖體強度和變形響應.

Subiquitous 模型［19］是廣義的Ubiquitous-Joint 模型，該模型中基質和節理強度符合雙線性摩爾-庫侖準則，且允許材料基質和節理的強度發生硬化或軟化.Subiquitous 模型和Ubiquitous-Joint 模型都是先根據摩爾-庫侖準則檢測基質的屈服，并進行相應的塑性修正，然后分析在新的應力狀態下節理面上的破壞，在材料未達到極限強度前力學行為一致.

在遍布節理模型中，弱面對巖體強度的影響通常與Jaeger 提出的單弱面理論［20］進行比較.單弱面理論指出，當1 -tan?tanβ>0 時，若滿足式（3）則會發生結構面的剪切破壞.

式中：c、?分別為結構面的黏聚力和內摩擦角；β為結構面的傾角.當1 -tan?tanβ<0 時，巖體不會沿結構面破壞，只會發生基質的破壞.故該理論只允許出現沿結構面的剪切滑移破壞和基質的破壞兩種破壞模式.

圖2 為Ubiquitous-Joint 模型［19］和Jaeger 單弱面理論的巖體承載強度與結構面傾角的關系的對比，可以看出兩者緊密匹配.圖2 中?w為結構的內摩擦角，其中曲線為帶有“肩部”的“U”形曲線.當β

圖2 Ubiquitous-Joint模型三軸抗壓強度值與Jaeger解析解的比較Fig.2 Comparison of triaxial compressive strength values-Ubiquitous-Joint model versus analytical solution

圖3 為部分已有的層狀巖體三軸壓縮試驗研究［21-24］，由圖3 可知，巖體的強度隨著層面傾角連續變化，這一特征也得到了許多研究人員的驗證.而單弱面理論不能充分描述自然存在的層狀巖體的各向異性.遍布節理模型也存在同樣的局限性，故需要進一步探討其適用性.

圖3 層狀巖體三軸壓縮強度隨傾角變化特性Fig.3 Variation of triaxial compressive strength of layered rock mass with bedding angle

2 遍布節理模型與離散元模型的對比

為了探討遍布節理模型對層狀巖體模擬的有效性，針對已有的層狀頁巖三軸壓縮試驗結果，采用FLAC3D建立與試樣同等規模的數值模型，用其內置的Ubiquitous-Joint 模型和Subiquitous 模型進行分析計算，并與3DEC的模擬結果作比較.

2.1 三軸壓縮試驗模擬中的比較

2.1.1 塊體離散元方法和Ubiquitous-Joint模型

為了研究層狀巖體的強度和變形特性，參考頁巖［22］的三軸壓縮試驗數據（如圖4 所示），使用3DEC建立了直徑50 mm、高100 mm 的標準圓柱體模型.層理傾角分別設置為0°、15°、30°、45°、60°、75°和90°，層厚5 mm，巖體參數標定結果見表1.同時基于Ubiquitous-Joint 本構模型建立了類似的FLAC3D 模型，將表1 中的巖體參數直接用作模型中巖石基質和節理的參數輸入，3DEC 模型和Ubiquitous-Joint 模型的強度響應如圖4所示.

表1 巖石和節理力學參數Tab.1 Rock and joint mechanical properties

圖4 離散元和Ubiquitous-Joint模型強度各向異性曲線Fig.4 Anisotropic strength curves of discontinuum and Ubiquitous-Joint model

正如預期，離散元模型隨著β角的增大而遵循連續變化的強度曲線.其與室內試驗不同傾角下的峰值強度相對誤差小于8%，結果基本吻合.Ubiquitous-Joint模型在β角小于15°時其強度不受節理的影響，與室內試驗結果相差超過20%，這種“U”形強度曲線上的肩部清楚地表明了模型的局限性.

提取較為典型的層理傾角為60°時巖石破壞模式的試驗結果與模擬結果，如圖5 所示.可知此時巖石表現為沿層理面的滑移破壞，其中從離散元模型結果可以看到層理面的錯動，與試驗結果一致.而Ubiquitous-Joint 模型顯示大量的節理剪切破壞，但無法得知具體的破裂面位置和破裂形態.

圖5 層理傾角60°時巖石破壞模式Fig.5 Failure mode of rock with bedding angle of 60°

圖6表示了不同傾角下離散元模型和Ubiquitous-Joint 模型的彈性模量和應力-應變曲線.由圖6 可知，Ubiquitous-Joint 模型沒有體現出峰后的應變軟化行為.當直接在模型中采用3DEC 巖石塊體的剛度參數時，所得到的彈性模量明顯高于3DEC的模擬結果，這是Ubiquitous-Joint 模型未考慮節理剛度和節理間距導致的，在實際工程中要特別注意這一點.

圖6 兩種模型的彈性模量和應力-應變曲線Fig.6 Elastic modulus and stress-strain curve of discontinuum and Ubiquitous-Joint model

因此，建議不要將3DEC中的巖石塊體和節理參數直接用作Ubiquitous-Joint 模型的參數，為使其產生有意義的結果，需要對巖石基質和節理參數進行校準，以匹配離散元模型的結果.以下將對此進行探討.

2.1.2 考慮應變軟化的Subiquitous模型參數校準

與Ubiquitous-Joint 模型相比，Subiquitous模型在校準巖石基質和節理參數方面提供了更大的靈活性.通過雙線性軟化關系，可以更好地表示層狀巖體的強度和變形特性.其參數校準準則如下［18］：

1）將離散元模擬結果視為實際層狀巖體的各向異性行為.

2）節理黏聚力和內摩擦角的初始值不變，巖體達到峰值后，節理黏聚力與巖體基質黏聚力以相同的速率軟化至0.

3）校準巖石基質的強度和變形響應，以補償節理剛度和節理間距參數的缺失.

β在0°和90°的情況下，試樣的峰值強度取決于巖石基質的黏聚力和內摩擦角，這些參數對應于β=0°時的離散元模型的強度響應進行校準.巖體基質和節理黏聚力的軟化速率參考離散元模型的結果.

在整個校準過程中，強度和剛度參數以及試樣的破壞過程都得以考慮.比較離散元模型和Subiquitous 模型的破壞模式，將其分為劈裂張拉破壞（β為90°時）、剪切滑移破壞（β為60°時）和復合破壞（β為30°時）.通過監測加載過程中基質的屈服和節理的滑移剪切，可以揭示試樣的破壞機制.

前文中的三軸壓縮數值試驗已用Subiquitous 模型重建，采用經過校準的參數，具體取值見表2.

表2 校準的Subiquitous模型力學參數Tab.2 Calibrated mechanical properties of Subiquitous model

離散元模型和Subiquitous 模型的各向異性“U”形曲線如圖7 所示，并與開始的Ubiquitous-Joint 模型的結果進行了比較.經過校準后的Subiquitous 模型隨著β角的增大同樣遵循連續變化的強度曲線，與離散元模型的結果更加貼切.圖8 顯示了Subiquitous 模型在不同層理傾角下的應力-應變曲線和彈性模量的變化，都與離散元模型更緊密地匹配.

圖7 離散元和Subiquitous模型強度各向異性曲線Fig.7 Anisotropic strength curves of discontinuum and Subiquitous model

圖8 彈性模量變化曲線及不同角度下的應力-應變曲線Fig.8 Elastic modulus and stress-strain response of discontinuum and Subiquitous model

2.2 二維圓形隧道開挖分析

為了驗證2.1 節中開發的校準后的Subiquitous模型在工程中的應用效果，建立了一個圓形隧道模型，研究隧道開挖后的力學響應，該模型是在不考慮重力加速度的各向同性應力場中模擬的.為了比較模擬效果，建立了巖層厚度0.5 m 的3DEC 模型和等效的Ubiquitous-Joint 模型.模型參數取值見表3 和表4.圖9 比較了隧道開挖完成時每個模型的塑性區、位移和最大主應力.

表3 模型材料參數Tab.3 Details of model material

表4 校準的Subiquitous模型力學參數Tab.4 Calibrated mechanical properties of Subiquitous model

圖9 5 m直徑的圓形隧道開挖后的塑性區、位移和最大主應力Fig.9 Plastic zones，displacement and major principal stress around 5 m-diameter tunnel

3DEC 模型中顯示隧道側壁中有少量巖石塊體的拉伸破壞，層理的剪切滑移破壞主要在洞頂和底部沿垂直層理方向延伸約2 m.位移場分布明顯受到了層理的影響，最大位移約為20 mm.圖10顯示的是3DEC 模型放大20 倍的變形狀況，在隧道頂部和底部可以清楚地看到巖層的彎曲.

圖10 3DEC中顯示的巖層彎曲變形（放大20倍）Fig.10 Bending deformation of bedding rock sown in 3DEC（magnified 20 times）

Ubiquitous-Joint 模型（直接對巖石基質和節理賦予和3DEC中塊體和節理相同的參數）中沒有顯示出基質的屈服破壞，而主要為節理的滑移和拉伸破壞，在隧道頂部和底部延伸約4 m，模型的最大位移明顯小于3DEC 模型.節理拉伸破壞導致區域周圍出現顯著的應力重分布，其破壞機制是因為遍布節理模型的公式中沒有考慮巖層的彎曲剛度.經過校準的Subiquitous 模型的隧道側壁上也有少量的基質拉伸破壞，使得隧道周圍出現更具有代表性的應力重分布，其位移場也更接近3DEC模型.

3 工程驗證

為了更好地研究Subiquitous 模型在實際工程中的使用性能，以新華山隧道為例，探討隧道開挖以及在支護結構作用下圍巖的變形和破壞特性，并通過現場實測數據驗證模型的可靠性.

3.1 工程概況和工程地質

新華山隧道位于湖南省張家界市和湖南省湘西州永順縣境內.該隧道為單洞雙線隧道，起止里程為DK26+104.00-DK32+034.49，全長5 930.49 m，最大埋深約為383 m，開挖高度和寬度分別為12.64 m 和14.96 m.

新華山隧道所處地貌為剝蝕低山地貌，地勢較起伏，山坡自然坡度一般為30°～70°.隧道穿越地層受區域構造影響嚴重，節理裂隙發育、巖體破碎.本文以新華山隧道進口段DK26+490斷面附近為研究對象.根據前期地質勘查資料，新華山隧道圍巖主要為層狀特征較明顯的炭質頁巖，由于其所具有的各向異性和開挖后風化較快等特殊工程特性，使得隧道的開挖引起軟弱圍巖向洞內發生不均勻對稱的變形.

3.2 模型建立

根據縱斷面圖可以發現，所模擬區段的埋深約110 m，運用FLAC3D 建立如圖11 所示模型.為降低模型中的單元數量，僅在模型中創建部分上覆巖體，并通過在地層上表面施加荷載模擬其余上覆巖體的自重應力.設定模型x、y、z三個方向上的尺寸分別為100 m、50 m 和100 m，采用位移邊界條件，除上表面外，其余5 個邊界面約束法向位移.模型中，巖體層理傾角采用現場調查得到的層理傾角，即為75°.隧道采用三臺階法開挖，模擬區段并未施做二次襯砌，故模型中支護體系僅包括錨桿和初期支護，相關力學計算參數根據支護設計方案確定（見表5）.采用3DEC 建立同等規模的模型，根據現場測試以及《鐵路隧道設計規范》（TB 10003―2016）取得如表6所示參數.其中節理剛度參數參考文獻［25］，并執行2.1節的校準程序取得Subiquitous 模型的參數，如表7所示.

表5 支護結構計算參數Tab.5 Parameters for the support system

表6 巖石和層理面力學參數Tab.6 Rock and bedding plane mechanical properties

表7 校準的Subiquitous模型力學參數Tab.7 Calibrated mechanical properties of Subiquitous model

圖11 數值模型及細部構造（單位：m）Fig.11 Numerical model and detailed construction（unit：m）

3.3 數值模擬結果與分析

根據上述參數和模型，計算得到隧道中部橫截面處（Y=25 m）開挖并施加初期支護后的圍巖變形和塊體塑性區以及節理塑性區情況如圖12 所示，從中可以看出：

圖12 離散元和校準的Subiquitous模型位移和塑性區對比Fig.12 Comparison of displacement and plastic zones of discontinuum and calibrated Subiquitous model

1）兩種模擬方法的圍巖變形和塑性區響應非常接近，說明經過校準后的Subiquitous 模型能夠較好地體現層狀巖體的力學特性.

2）隧道開挖完成后實測拱頂沉降和水平收斂分別為259.9 mm 和173.5 mm，而3DEC 中拱頂沉降和水平收斂分別為243.9 mm 和160.2 mm，與實測值的差異分別為-5.4% 和-7.6%，FLAC3D 中分別為282.5 mm 和189.2 mm，與實測值的差異分別為8.7%和9.0%，差異性較小，表明建立的模型能夠較好地反映新華山隧道開挖及初期支護后的變形情況.

3）受層理的影響，拱頂和拱底都朝著層理傾角方向產生位移梯度，圍巖位移場呈現出顯著的不對稱性，這也與現場觀察到的非對稱變形情況相符合.FLAC3D 中邊墻附近圍巖位移比3DEC 稍大，是因為Subiquitous 模型無法表示完整巖層的屈曲變形，而巖層的厚度對圍巖位移有重要影響.

4）隧道開挖擾動作用下，圍巖產生了節理剪切破壞、節理張拉破壞、巖石基質剪切破壞和巖石基質張拉破壞4 種類型的破壞，主要處于節理和基質的剪切破壞狀態，且大部分圍巖體同時出現了多種破壞模式.圍巖塑性區分布也顯示為極不對稱性，圍巖深部的塑性區主要集中在左拱腳和右拱肩.3DEC中少量的節理張拉破壞沿洞周分布，FLAC3D中節理張拉破壞更少，這也與Subiquitous 模型無法解釋巖層間距有關.基質的張拉破壞只出現在隧道底部，拱頂的塑性區范圍都很小.

圖13給出了3DEC和FLAC3D模型（與實際掘進過程一致）Y=25 m 斷面處的拱頂沉降監測曲線與現場監測數據的比較，可以發現，3DEC 與實測數據更為接近，而FLAC3D 中采用Subiquitous 模型的計算結果也能較好地吻合.綜合以上分析，校準的Subiquitous模型在工程中有較好的實用性，能為相應工程提供參考，且層理的存在對圍巖的變形和破壞有重要影響.此外，就計算效率而言，兩者計算時長相差30～40倍.

圖13 實測和模擬的拱頂沉降（Y=25 m）Fig.13 Measured and simulated vault settlement（Y=25 m）

3.4 層理傾角對隧道開挖的影響

當隧道施工穿越炭質頁巖地層時，圍巖和支護結構的變形很可能因圍巖層理傾角的變化產生顯性差異.為分析層理傾角對圍巖和支護結構變形的影響，用FLAC3D 依次建立層理傾角為0°、30°、45°、60°和90°等5種工況的仿真模型.采用校準的Subiquitous模型，除傾角外其余參數保持不變.計算得到巖體圍巖和支護結構的變形以及圍巖塑性區分布，如圖14所示，提取各個角度下拱頂的沉降得圖15 所示曲線.由圖14、圖15可知：

圖15 拱頂沉降隨層理傾角的變化Fig.15 Vault settlement varies with bedding angles

1）圍巖和支護結構的變形顯著受到層理傾角的影響.層理傾角從0°到90°變化過程中，初期支護拱頂沉降呈現出倒“V”形變化，即先增大后減小，45°時達到最大值.位移場分布隨著傾角改變，只有0°和90°時存在對稱性.

2）隧道開挖引起的塑性區形狀和范圍與層理傾角密切相關.0°時塑性區范圍最小，當層理傾角小于30°時，圍巖深部塑性區沿垂直于層理方向發展；而當傾角為75°～90°時，深部塑性區主要沿層理方向發展；層理傾角為45°～60°時，塑性區呈現出“X”形狀，且范圍較大，與前文所述巖體在45°～60°時強度較低相對應，表明該傾角范圍內易使隧道圍巖產生破壞.

4 結論

本文通過對比分析遍布節理模型與離散元模型在層狀巖體三軸壓縮以及層狀圍巖隧道開挖應用中的模擬效果，探討采用等效參數的遍布節理模型在層狀巖體模擬中的適用性，得出以下結論：

1）離散元模型能夠更好地體現層狀巖體的變形和破壞特性，但若考慮計算效率，更適合于描述中小尺度層狀巖體力學性質；而遍布節理模型由于其本身對節理裂隙考慮的不足，在模擬層狀巖體時，需要對部分參數（彈性模量、泊松比以及巖石基質的黏聚力和內摩擦角）進行修正才能用于工程分析，且更適用于大尺度工程巖體的力學行為研究.

2）對于新華山隧道工程而言，兩種模型在網格單元劃分接近的情況下，計算效率相差30～40 倍，而校準的遍布節理模型得到的圍巖位移與實測結果分別相差8.7%和9.0%，差異性較小，表明該模型兼顧效率的情況下準確度良好.

3）層理弱面的抗剪強度和抗拉強度較低，故層狀巖體在工程擾動的情況下，容易造成層理剪切滑移破壞以及張拉破壞，在工程中要重點關注.