基于約束收斂法的Hoek-Brown 巖體巷道大變形分析

劉洪磊，張權云，關 凱，朱萬成，于慶磊

（東北大學資源與土木工程學院，遼寧 沈陽 110819）

0 引言

近年來，在巷道開挖過程中，伴有襯里開裂甚至坍塌的大變形，對巷道施工過程中的安全構成了重大威脅，特別是在高地應力條件下，深部軟巖巷道的大變形如幫壁膨出、底鼓等問題愈發突出，因此，從工程實際方面，需合理預測大變形條件下圓巷圍巖力學響應，為圍巖穩定性評估提供理論參考。而在理論分析方面，約束收斂法的原理是根據圍巖特征曲線與支護特征曲線相交，得到的交點即為圍巖與支護相互作用達到平衡時的應力與位移，由于巷道支護的目的是限制巷道位移量，因此結合工程經驗，可為深部軟巖巷道初期支護設計提供理論參考。而圍巖特征曲線可通過圓形巷道的彈塑性解得到，因此合理預測巷道變形和圍巖塑性區對支護優化設計和圍巖結構穩定性有重要影響。同時采用不同的理論分析與不同破壞準則計算的巷道圍巖應力狀態與位移分布可能會存在不同程度的差異。選擇合適的屈服準則和適用于深部實際巷道工程施工的分析方法進行深部軟巖巷道的變形機理研究對于深部巷道支護設計具有重要意義。

長期以來，國內外已有眾多學者對圍巖彈塑性分析開展了大量有關小變形理論分析與大變形理論分析的研究。許多學者基于小變形理論，采用約束收斂法對圍巖的變形特征展開研究[1-6]，盡管小變形理論很簡單，但是由于在建立控制方程時忽略了物質點的變化，因此在預測巷道施工后的巖體行為方面存在缺陷。特別是在高地應力條件下，可能會高估圍巖變形能力，有時甚至會導致預測的變形大于開挖半徑[7]。目前，地下金屬礦山的開采深度逐漸增加，而小變形理論通常適用于淺層、硬巖的環境條件，但深部巖石地下工程具有高地應力、高溫和高滲透水壓等復雜的環境條件，巖土材料具有幾何非線性特征，因此，小變形理論體系在深部巖石地下工程中不再適用。

近年來，針對上述缺陷，對深部巷道開挖的圍巖大變形分析的研究越來越多，例如VRAKAS[8]提出了基于Mohr-Coulomb 準則（以下簡稱“MC 準則”）和Hoek-Brown 準則（以下簡稱“HB 準則”）的應變軟化巖體在大變形理論條件下的開挖圓形巷道的解析解，但其基于理想彈塑性模型，未考慮圍巖應變軟化模型。ZHANG 等[9]提出了應變軟化巖體圍巖特征曲線的大變形解析解，發現大變形理論得出的位移、軟化區半徑均小于小變形解析解，且其差異隨著彈性模量的減小而增大。關凱等[10-11]在VRAKAS 研究的基礎上，考慮了巖體應變軟化行為，對MC 準則巖體中開挖圓形巷道后引起的圍巖響應進行研究。相關研究大多基于圍巖理想彈塑性或彈-脆-塑性巖體力學行為進行研究，且求得的解析解結構較復雜，而巖石的破壞過程實際上是強度與變形參數逐漸劣化直至消失的過程，因而無法準確反映深部巷道圍巖變形特征。

同時，針對地下金屬礦山，當前眾多學者對深部圓巷圍巖洞壁變形的理論解析解分析多采用MC 準則或HB 屈服準則[11-13]，但是由于地下礦山復雜的節理巖體條件，基于考慮地質強度指標和節理賦存狀態的HB 強度準則預測其變形和破壞特性具有適用性。在VRAKAS[8]研究基礎上，利用塑性區半徑對物質點當前坐標及位移進行無量綱轉換，并考慮強度與變形參數劣化的過程，基于HB 屈服準則提出了一種簡單的深部軟巖大變形預測方法，并以理論方法和前常銅鐵礦巷道（矽卡巖、破碎、膨脹）監測數據進行有效驗證，在驗證的基礎上進行了參數分析，并基于不同變形理論條件下圍巖特征曲線區別，針對軟弱破碎巖體提出支護建議。

1 理論介紹和理論方法簡述

1.1 理論介紹

在連續、均勻、各向同性、初始彈性的巖體中開挖一個半徑為R0的圓形巷道，如圖1 所示，且其在橫斷面上所受靜水應力為σ0，巖體彈塑性交界處即塑性區半徑為Rp，巷道內壁受內部支護力為σa。

對于以上的軸對稱平面應變問題，當前狀態下的應力平衡方程見式（1）。

式中：σr為巖體的徑向應力；σt為巖體的切向應力。

但需注意在大變形分析時，上述應力平衡方程應建立在巷道的已變形構型上，而變形的描述應該在圖1 中的圍巖穩定后的巷道半徑處。在該問題中對數應變可以準確地描述應變-位移關系，因此在當前極坐標條件下，計算徑向應變與切向應變的公式見式（2）和式（3）[8]。

當內部支護力逐漸減小時，圍巖發生屈服，隨著巷道開挖，圍巖逐漸到達彈性極限狀態，應力逐步達到巖石破壞的極限應力。而兩種最常用的巖體破壞極限主應力間的關系有非線性HB 屈服準則和線性MC 屈服準則，此處采用的是HB 屈服準則，計算見式（4）。

式中：σ1為巖體的最大主應力，因本文研究軸對稱問題，其值等于切向應力 σt，MPa；σ3為巖體的最小主應力，等于徑向應力 σr，MPa；σci為完整巖體的單軸抗壓強度，MPa；m為巖石的無量綱經驗參數，反映巖石的軟硬程度；s、a為無量綱材料常數。

在巖體發生屈服后，巖體強度突然下降，并遵循屈服后的軟化行為。根據雙線性函數計算應變軟化巖體的強度和變形參數，見式（5）和式（6）。

在大變形分析中，有限變形理論給出的應力、應變在Eulerian 和Lagrange 兩種坐標系下的轉換關系有利于非線性方程的求解。同時，通過有限變形理論分析大變形問題時，采用基于變形率張量加性分解的次彈塑性模型進行求解，可得出類似于小變形理論的本構方程及塑性流動法則，見式（7）。需注意，上述方程均針對變形后的當前構型，且后文中的計算過程中應用的本構方程與塑性流動法則參考式（7）[11]。

式中：C為四階彈性張量；σ為Cauchy 應力張量；Q為塑性勢函數；ε為應變張量；λ為塑性乘子。

相較于之前的小變形分析[1-3]，首先采用有限變形次彈塑性理論，不同于建立在變形前的初始構型上的小變形分析，應力邊界和平衡方程均建立在已變形結構上；隨后定義無量綱當前位置與無量綱位移，見式（8）和式（9），方便對數應變、圓環內任一點物質點處的半徑與位移的無量綱變換，簡化后續理論推導。

1.2 方法簡述

利用有限差分方法，將塑性區（包括塑性軟化區和殘余區）劃分為一組同心圓環（圖2），其中，ρ（j）和ρ（j+1）為第j環內外邊界的無量綱半徑；j=1和j=n+1分別為對應于彈塑性邊界和巷道洞壁；εrRp和 εtRp分別為彈塑性邊界處的徑向應變和切向應變。因此，σrRp、σt和 εrRp、εtRp等于初始值 σr1、σt1和 εr1、εt1。假設每個塑性區圓環的徑向應力增量 Δσr都是恒定的，可通過式（8）計算。

圖2 塑性區劃分圓環Fig.2 Concentric annulus in plastic zone

式中：n為劃分的同心圓環總數；σr（j）為在r=r（j）（j=2,3,...,n+1）處的徑向應力。

同時，為了便于后續理論推導和數值編程，定義了物質點無量綱當前位置及無量綱位移，見式（9）和式（10）。

當第j個塑性圓環足夠小時，可假設無窮小圓環內巖體內徑、外徑受力近似與脆塑性巖體響應相同，如圖3 所示。據此進行解析解推導時可將彈脆塑性巖體內徑、外徑所受徑向應力 σa和 σRp分別用 和代替，彈塑性交界物質點處半徑Rp和巷道開挖變形后半徑a分別用 ρ（j）和 ρ（j+1）代替，且第j個塑性圓環可視為均質圓環，即第j環內任意物質點的力學性質均看作與 ρ=ρ（j）處的受力情況相同。

圖3 第j 個圓環受力分布圖Fig.3 Stress distribution in j annulus

根據式（8）中的塑性圓環各物質點分布規律，ρ（j+1）處物質點的徑向應力計算見式（11）。

若劃分的圓環數目足夠多，即n值足夠大，則第j個塑性圓環內巖體的徑向和切向應力的解析式可參考VRAKAS 等研究的HB 巖體的理想彈塑性應力解析表達式[8]求得。

2 理論驗證

VRAKAS 等[8,14]推導出基于HB 準則和MC 準則，圓形巷道軟弱圍巖的理想彈塑性大變形理論解析解，而基于以上研究，提出了大變形理論下，基于HB 破壞準則的巷道應變軟化數值解法。為了能與VRAKAS等的理想彈塑性結果進行對比，采用臨界塑性內變量值為0.000 01 與10 000，使研究退化為彈脆塑性問題和理想彈塑性問題。相關計算參數[5]：r0=2、μ=0.3、mp=1.7、sp=0.003 9、mr=0.85、sr=0.001 9、E=5.7 GPa、Ψp=0、Ψr=0、σcip=30 MPa、σcir=25 MPa、ap=0.55、ar=0.6、σ0=15 MPa。計算得出與VRAKAS 等[14]彈脆塑性和理想彈塑性對比結果，如圖4 和圖5 所示（VRAKAS結果用折線圖表示，本文研究結果用散點圖表示）。由圖4 和圖5 可知，結果能夠很好地吻合，表明提出的大變形解析解能有效預測圍巖擠壓響應。

圖4 理想彈塑性巖體結果對比Fig.4 Comparison of ideal elastic-plastic rock mass results

圖5 彈脆塑性巖體結果對比Fig.5 Comparison of elastic brittle plastic rock mass results

3 前常銅鐵礦巷道變形預測應用

為便于分析，以前常銅鐵礦-480 m 中段試驗巷道為例進行擠壓巷道圍巖變形預測應用探究，通過巖體基本參數計算理論推導條件下的松動圈、收斂位移以及損失位移大小，進而與監測數據對比來驗證提出模型和方法的有效性與準確性。根據圖6 所示近似方形的巷道斷面尺寸為2.5 m×2.5 m，若將巷道視為等效圓形，可以得到巷道初始設計等效半徑為1.25 m。

圖6 -480 m 試驗巷道初始設計斷面Fig.6 Initial design cross-section of -480 m test tunnel

根據過往的研究[15-16]可知，巖體的峰后本構模型取決于巖石類型及地質強度指標GSI：對于較完整脆性巖體（GSI＞75），巖體力學行為一般表現為彈脆塑性行為；對于中等節理巖體（25＜GSI＜75），巖體一般表現為應變軟化行為；對于破碎巖體（GSI＜25），巖體本構模型通常表現為理想彈塑性本構模型。由于試驗巷道圍巖主要為矽卡巖，根據有關矽卡巖強度和膨脹變形機制的研究可知，試驗巷道圍巖的矽卡巖大部分屬于軟巖，因此可采用應變軟化模型（圖7）研究矽卡巖巷道開挖的力學響應。根據ZHANG 等[6]研究，此時臨界塑性內變量ηc的取值為0 與無窮大之間（在分析中以ηc=0.114 1 進行等效計算）。

圖7 應變軟化模型Fig.7 Strain softening model

前常銅鐵礦-480 m 中段試驗巷道斷面尺寸、原巖應力及矽卡巖力學參數見表1。值得注意的是，在分析時假設施加的支護措施相對于變形壓力可以忽略不計，即巷道開挖后支護結構均發生完全毀壞σσ=0；此外，靜水原巖應力大小采用自重應力γh進行計算，其中，γ和h分別為巖體容重和埋深，即σ0為12.48 MPa。同時，為了確保數值結果的精確度在可控范圍內，彈性區物質點數ne以及塑性區物質點數np均取100，計算容差Tol=1e-7。

表1 -480 m 中段試驗巷道相關計算參數Table 1 Calculation parameters related to the -480 m middle section test tunnel

于慶磊等[17]、褚吉祥等[18]對-480 m 巷道分別進行了矽卡巖松動圈和巷道最大收斂變形的現場測量，為了驗證提出方法的準確性，將現場監測結果值與提出的理論模型數值結果進行對比驗證。巷道初始設計半徑為a0=1.25 m，表2 為計算結果與矽卡巖巷道現場監測（探測）結果的對比。

表2 巷道開挖大變形理論與現場監測結果對比Table 2 Comparison of the tunnel excavation based on the large deformation theory and field monitoring results

由表2 可知，基于理論分析得到的塑性區大?。?.14 m）位于現場監測松動圈范圍（2.5～3.8 m）之內，且平均相對誤差僅為0.32%。而在巷道收斂變形預測方面，理論計算結果（76.00 mm）遠大于現場監測結果（48.00 mm），且平均相對誤差高達58.3%；由于理論計算得到的巷道收斂變形為巷道開挖擾動引起的總變形量，而現場監測的收斂變形僅是在多點位移計安裝之后的部分位移，使得在設備安裝之前發生的圍巖位移監測不到（圖8），這部分位移也稱損失位移，其可通過超前鉆孔測量或基于巷道縱向變形曲線等理論方法獲得，因此造成收斂變形的理論計算結果遠大于現場監測值。

圖8 損失位移概念圖Fig.8 Loss displacement

根據表2 可以反演出多點位移計安裝前圍巖的損失位移為28 mm。由此可知，考慮到多點位移計監測收斂變形的局限性以及基于理論模型得到的巷道塑性區與實測松動圈范圍基本一致，可以認為研究提出的理論方法是正確有效的，能較好地預測大變形應變軟化圍巖巷道的開挖響應。

4 巖體參數分析

4.1 強度參數a 對廣義HB 巖體彈塑性解的影響

研究提出的有限差分方法中所有的強度參數均可通過塑性內變量求出（式（6））。但是HOEK 等[16]研究得出，在HB 破壞準則中，強度參數a不再是固定的0.5，而是隨著地質條件變化，經驗關系見式（12）。

式中，GSI為反映巖體斷裂面狀況的地質強度指數，GSI的范圍為10～100，故a在0.5～0.6 之間變化。因此，為了正確使用該準則，需分析強度參數a是如何影響彈塑性行為的。

采用以下參數研究強度參數a對圍巖特征曲線和應力位移分布曲線的影響：r0=2 m、μ=0.3、mp=1.7、sp=0.003 9、mr=0.85、sr=0.001 9、E=5.7 GPa、Ψp=0°、Ψr=0°、σcip=30 MPa、σcir=25 MPa、ap=0.55、ar=0.6、ηc=0.004、σ0=15 MPa。圖9 為假設的不同強度參數a條件下的四種應力及位移分布圖。其中三條曲線假設強度參數a的殘余值和峰值分別為0.50、0.55、0.60，而其中一條曲線假設參數a隨著臨界塑性內變量 ηc從0.50 到0.60 變化，其他輸入數據保持不變。由圖9 可知，當內部支護力 σa很小時，圖中的曲線差別較大，且ap與ar不等時，曲線更接近ap=0.60 時的圍巖特征曲線，這說明峰值強度參數對結果影響更大；而a值較大時，塑性區半徑也大，如圖10 所示，說明峰值強度參數對塑性區半徑也有一定影響。從圖9 還可以看出，a值越大洞壁收斂位移也會越大，而由圖11 可知，ap=ar=0.60 是其等于0.50 位移的三倍，故HB 準則中的強度參數a值對巖體開挖響應的影響相對較大，同時a值越大巖體更易破壞。

圖9 參數a 對應力及位移分布的影響Fig.9 Influence of the parameter a on the stress and displacement distribution

圖10 洞壁收斂位移、塑性區半徑與強度參數a 的關系曲線Fig.10 Relationship curves of cave wall convergence displacement,plastic zone radius and strength parameter a

圖11 參數m 對圍巖特征曲線的影響Fig.11 Influence of parameter m on the characteristic curve of surrounding rock

4.2 強度參數m、s 的影響

HOEK 等[16]提出強度參數m、s也與地質強度指標GSI相關，可通過經驗公式計算，見式（13）和式（14）。

式中：D為應力釋放或爆炸破壞對巖體擾動程度的系數；mi為組成巖體的完整巖塊的Hoek-Brown 常數；m、s為巖體的強度常數，其中，m為巖石的軟硬程度，其取值范圍為0.000 000 1（嚴重擾動巖體）～25（理想完整的堅硬巖體）；s為巖體破碎程度，其取值范圍為0（極破碎巖體）～1（理想完整巖體）。地質強度指標GSI與巖體的表面粗糙性和表面風化程度、結構特征等特征有關。因此，強度參數m和s能夠綜合反映巖體結構、巖礦狀態和巖體不連續面質量，為了準確分析巖體質量與強度參數的變化關系，從而正確使用該準則，還需分析強度參數m、s是如何影響彈塑性行為。

為了探討強度參數m對圍巖特征曲線的影響，繪制參數m對圍巖特征曲線的影響曲線圖，如圖11所示。由圖11 可知，隨著強度參數的增大，相同內部支護力條件下的位移對應減小，說明m越大巖體強度越高，而且圍巖特征曲線的增幅與強度參數的增幅不同，相同內部支護力條件下的位移增幅明顯更大，強度參數在0.6、1.3、2.0 值時隨著強度參數的增大，位移逐漸減小且減小速度快，而在內部支護力值大于8 時，位移與塑性區半徑增加速度逐漸平緩。因m反映巖石的軟硬程度，m越大巖體越完整堅硬、巖體內的節理裂隙越少，以上變化規律說明軟巖對圍巖特征曲線、位移及塑性區半徑的影響較大，而硬巖對位移及塑性區半徑幾乎沒有影響，即相對較硬的巖體經開挖產生的洞壁收斂位移與塑性區半徑相對較小。

同時強度參數s對圍巖特征曲線、位移及塑性區半徑的影響和強度參數m的影響類似，如圖12 和圖13 所示。由圖12 和圖13 可知，相對破碎巖體產生的位移及塑性區范圍比相對完整巖體的位移及塑性區范圍更大。同時相比于強度參數m，巖體力學行為對強度參數s更敏感。

圖12 強度參數m、s 對塑性區半徑的影響曲線Fig.12 Influence of parameter m、s on the plastic zone

圖13 參數s 對圍巖特征曲線的影響Fig.13 Influence of parameter s on the ground curve

4.3 臨界塑性內變量的影響

理想彈塑性和彈脆塑性是巖體力學行為的兩種特殊情況，通過調節臨界塑性內變量的數值大小，使得巖體基本力學參數在峰值參數與殘余參數之間變化，從而達到兩種特殊狀態。為了直觀體現不同臨界塑性內變量對巖體巷道開挖響應的影響，選取巖體基本力學參 數：r0=3 m、μ=0.25、mp=2、sp=0.004、mr=0.6、sr=0.002、E=5.7 GPa、Ψp=1.55°、Ψr=1.55°、σcip=30 MPa、σcir=25 MPa、ap=0.506、ar=0.530、ηc=0.01、σ0=15 MPa、σa=0.75 MPa，繪制了應力分布、洞壁收斂位移與臨界塑性內變量的關系曲線，如圖14 和圖15所示。由圖14 可知，隨著臨界塑性內變量的增大，應力、塑性區半徑以及收斂位移均在不斷減小，因此不同的臨界塑性內變量值將導致巷道圍巖塑性區應力分布的不同。由圖15 可知，在應變軟化階段，臨界塑性內變量的輕微降低將會引起較大巷道圍巖變形及洞壁位移增量，同時，臨界塑性內變量約大于0.1 或臨界塑性內變量約小于0.001 時，洞壁收斂位移趨于平緩，此時臨界塑性內變量幾乎沒有影響，說明此時比較接近巖體力學行為的兩種特殊情況——彈脆塑性與理想彈塑性本構模型。

圖14 臨界塑性內變量對應力分布的影響Fig.14 Influence of critical plastic internal variables on stress distribution

圖15 洞壁收斂位移與臨界塑性內變量的關系圖Fig.15 Relationship between wall convergence displacement and critical plastic internal variables

4.4 小變形和大變形理論下的圍巖響應

在圍巖與支護耦合作用中，收斂約束法是以巷道開挖后圍巖彈塑性理論為基礎、參考現場監測數據、結合工程經驗判斷的支護設計方法，其主要由圍巖特征曲線和支護特征曲線組成[11,19]。為了探討不同理論方法對工程施工過程中巷道圍巖變形預測及支護受力分析的影響，基于約束收斂法，根據不同分析方法即大變形與小變形理論及破壞準則，研究不同條件下的圍巖特征曲線的區別。

根據約束收斂法，圖16 為采用不同變形理論及破壞準則條件下的圍巖-支護相互作用曲線。當在大變形理論分析中，采用MC 破壞準則時，圍巖特征曲線與支護特征曲線相交于圖中C 點處，表明支護與圍巖平衡時的支護所受壓力為σaC；如果在大變形分析中，采用HB 破壞準則，則此時圍巖特征曲線在彈性階段的應力與位移與MC 準則的基本一致，而在塑性階段，相同洞壁收斂位移條件下，支護壓力反而較小，與支護特征曲線相交于D 點；同時當采用小變形理論分析MC 巖體圓形巷道開挖時的圍巖響應，圍巖與支護特征曲線在A 點相交，此時會高估巷道圍巖的擠壓勢；而當采用小變形理論分析HB 巖體巷道開挖時，兩曲線相交于B 點，由圖16 可知，對于分析方法，大變形理論明顯比小變形理論結果更準確預測圍巖變形，而針對破壞準則分析，HB 破壞準則明顯比MC 破壞準則更準確。因HB 準則預測節理巖體更準確，因此基于HB 準則的大變形計算方法更準確，但在深部軟巖巷道開挖的高地應力賦存條件下，為了防止因提供的支護力不足導致圍巖-支護整體結構失穩，采用基于MC 準則的大變形理論進行支護設計既節約成本又可為支護結構破壞預留一定壓力。

圖16 擠壓巖體圍巖-支護相互作用曲線Fig.16 Interaction curve of surrounding rock-support of extruded rock mass

5 結論

在VRAKAS 等[8,14]研究的基礎上，采用Hoek-Brown 破壞準則，提出了針對圓形巷道軟弱圍巖應變軟化條件下的大變形理論方法，研究了m、s、a等不同強度參數對圍巖響應的影響，并基于約束收斂法通過大變形與小變形理論分析結果對比，為支護設計提供分析手段和建議，得到結論如下所述。

1）根據平衡方程、Hoek-Brown 破壞準則等提出相對簡單的應變軟化大變形分析方法。

2）通過理想彈塑性和彈脆塑性理論結果對比以及工程實例中的松動圈以及收斂位移大小對比驗證本文提出模型和方法的有效性與準確性。

3）研究了不同強度參數及理論方法對圍巖響應的影響：隨著強度參數m、s逐漸增加，巷道圍巖的位移及應力逐漸減小，而強度參數a的影響相反。大變形與小變形理論方法條件下的圍巖特征曲線對比，說明采用大變形分析方法，所得支護壓力與收斂位移相對更小，結果更準確；而MC 準則比HB 準則所得洞壁收斂位移更小，結果更準確，但在高原巖應力條件下，則很可能由于支護提供的壓力不足以抵抗圍巖破壞從而導致支護結構發生破壞或巷道失穩，故基于MC 準則的大變形理論預測的支護效果相對HB 準則更偏于保守。