分段黏結非連續變形分析方法及其在砂巖破裂分析中的應用

鄭 飛，鄧慶龍，李 芷，王培瑞，靳 陸，焦玉勇

(中國地質大學(武漢) 工程學院,湖北 武漢 430074)

深部煤炭開采中，高地溫、高地壓、高巖溶水壓及開采擾動影響更為顯著，巖石力學問題也更為復雜[1]。沖擊地壓、煤與瓦斯突出災害與巖石的多尺度斷裂行為密切相關，在實驗室尺度研究煤層頂底板巖層代表性巖石在動靜荷載作用下的變形、斷裂、強度等特性及演化規律，對于認識上述災害的發生規律具有重要的意義[2-5]。

大量試驗及相關研究表明，巖石的變形及強度特性與細觀尺度的顆粒、微裂紋等特征具有較強的相關關系。精細刻畫巖石細觀結構，考慮顆粒及微裂紋等巖石細觀特征，進行宏觀變形及強度特性的分析是重要的研究課題。在巖石的室內實驗方面，聲發射[6]、CT 掃描[7]、3D 打印[8]等先進的技術得到運用，推動了對巖石強度特性、動力學特性的研究。

數值模擬技術由于其靈活性、定量性、過程可視性等特點，也是巖石多尺度斷裂分析研究的重要手段。分析巖石斷裂過程的數值方法，通常可劃分為基于連續介質假定的方法、基于非連續介質假定的方法、以及混合連續與非連續框架進行求解的方法。基于連續介質假定的方法包括擴展有限元(XFEM)[9]、巖石破裂過程分析方法(RFPA)[10]、邊界單元方法[11]等；非連續類數值模擬方法主要包括基于顯式接觸力求解的離散單元法(Discrete Element Method,DEM)和基于隱式位移求解的非連續變形分析方法(Discontinuous Deformation Analysis，DDA)[12]；混合連續與非連續框架求解的方法包括FDEM[13-14]、CDEM[15]、拉格朗日元-離散元耦合法[16]等方法。

非連續類數值計算方法可以考慮較大尺度斷層、節理、層理等不連續面的張開、閉合、滑移，也可以模擬較小尺度的巖石斷裂問題，其中DDA 方法在接觸的數學計算、塊體的積分、計算收斂性方面具有嚴密的數學基礎。不同于采用顯式時間積分計算的DEM，DDA 可視為基于Newmark 隱式積分和接觸開閉迭代求解的隱式非連續計算方法，在離散塊體系統大變形、大位移計算中可以保持時間步取值的無條件穩定。DDA 方法最初用于解決裂隙切割的塊體系統的靜態穩定性及動態響應問題。KE[17]提出采用子塊體劃分DDA 方法模擬裂隙巖體變形問題，焦玉勇等[18]則通過引入節理約束的三角形子塊體剖分及黏結單元來模擬節理巖體的漸進破裂過程，甯尤軍等[19]和ZHENG 等[20]相繼對黏結單元的破壞準則、單元積分模型進行修正和優化。面對深部開采中與巖層破斷、非連續大變形相關的巖石力學問題，DDA 方法具有天然的優勢。許多學者利用DDA 方法研究了深部沿空掘巷過程[21]、煤層開采導致的地表沉降[22]、煤礦巖體大變形和開挖過程[23]、礦山陷落柱[24]等問題。

非連續類數值模擬方法已成為分析巖石斷裂問題的重要手段，傳統子塊體DDA 方法基于邊-邊接觸關系和法向、切向強度參數的修正模擬斷裂過程，對復雜巖石材料變形破壞過程的適用性和計算精度有所欠缺?；陔x散表征黏結單元力學破壞過程的樸素思想，筆者提出了一種分段式黏結模型，提升傳統黏結模型在模擬復雜本構關系時的適用性，更準確地模擬巖石動靜荷載下的斷裂問題。

1 基于分段式黏結模型的DDA 方法

1.1 子塊體剖分DDA 方法的基本原理

巖體可以看作是由節理、裂隙等不連續面切割形成的斷續介質系統?？紤]不連續面的在幾何方面的有限延展性和在力學方面的非連續變形特性，采用考慮非連續界面約束的子塊體幾何剖分和黏結/接觸單元進行力學近似求解，是常用的非連續計算方法[17]。子塊體剖分后，黏結單元使用由巖石性質等效的虛擬節理參數，接觸單元則根據真實節理性質進行賦值[18]。

1.1.1數學原理與控制方程

離散塊體系統可通過連接或接觸關系形成統一的分析對象。塊體中任意一點p(x,y)的位移向量u={u,v}T可表示為u=Td，其中，T為位置矩陣；d為表征塊體剛體位移和應變自由度變量的向量[12]：

式中，u0、v0為塊體沿x軸和y軸的平動位移量；r0為塊體在xoy平面的旋轉角(逆時針為正)；εx、εy、γxy分別為x軸和y軸方向的應變和剪切應變。

p點的位置矩陣T可由式(2)計算[12]：

式中，x0、y0為塊體的形心點坐標分量。

塊體之間的接觸/黏結關系可由多種方式進行求解，例如罰函數法、拉格朗日乘子法等。由于不增加額外的自由度，以及對于隱式求解中矩陣對陣正定性質的滿足，罰函數法使用更為廣泛?？紤]所有塊體之間的接觸/黏結關系以后構造總勢能泛函，基于最小勢能原理可獲得塊體系統的動力平衡方程[20]：

式中，M為質量矩陣；C為阻尼矩陣；KE為彈性剛度矩陣；KB為黏結單元剛度矩陣；KC為塊體接觸矩陣；F為系統外力向量；D為塊體的廣義位移自由度向量，其上1 個點表示對時間的一階導數，2 個點表示對時間的二階導數。

1.1.2時間積分方案與隱式求解模式

求解動力學方程，Newmark 積分法是常用的一種積分方案。DDA 方法中原始的常加速度積分方案可視為Newmark 積分方法的一個隱式求解方案的特例[25]。從tn時刻到tn+1時刻，間隔為Δt的求解步中，常加速度積分方案可表示為

將式(5)代入動力平衡方程，可獲得每一時步的位移增量Dn+1的求解公式：

由此可見，每一個計算步中需要求解一個整體矩陣，可以證明該矩陣滿足對稱正定特征[12]。

1.2 分段式塊體黏結模型

基于可變形塊體及黏結/接觸的方法可實現巖石材料的連續變形及非連續破壞過程分析。常用的子塊體DDA 方法中，每個子塊體均為常應變彈性單元，采用考慮拉伸破壞的修正摩爾-庫倫準則，破裂發生在子塊體邊界上，不考慮塊體內部發生破壞。黏結模型可視為一種特殊的接觸模型，可模擬黏結單元兩側塊體的連續變形、黏結單元的損傷及失效?；诹P函數的處理方法，在子塊體間引入黏結模型，可以同時模擬分析域的連續變形、漸進破裂、以及破裂后的動力學過程。隨著對巖體及巖石材料精細力學模擬的發展，傳統的黏結模型需要更新以更好的還原巖體結構/巖石材料的原始幾何形態及物理狀態，尤其是對于巖石裂紋面的非線性破裂過程，因此本研究提出分段式黏結單元模型。

1.2.1分段式黏結模型概述

基于離散化表征的思想，筆者提出了一種分段式單元黏結模型，以更準確的表征塊體間黏結的逐步損傷失效過程。如圖1 所示，黏結模型在幾何上從屬于2 個離散可變形塊體的2 條重疊邊，在力學上通過本構關系定義2 個邊的相對位移與作用力的關系。幾何方面，長度為lb的一對黏結邊可離散化為長度均等的nb個線段，在每對線段的中心點沿著邊的法向和切向可施加一對罰彈簧進行黏結對邊的位移約束。分段式黏結模型主要包括以下參數：黏結分段數量nb，黏結長度lb，法向黏結剛度切向黏結剛度和可定義為

其中，Ec、h分別為塊體的特征彈性模量和特征長度；αb為黏結剛度比，表征法向黏結剛度與塊體彈性模量和特征長度的比例關系；βtn為切向剛度比，表征切向黏結剛度與法向黏結剛度的比值?；诜侄尾呗?，黏結單元內非線性應力分布、非線性變形-應力關系可以等效為分段線性的模式進行表征和求解，而每個分段的應力、變形由設置在分段形心點處的黏結彈簧模擬。

1.2.2黏結本構關系及分段式表征

邊-邊黏結單元的平均正應力σm可由式(9)獲得：

邊-邊黏結單元的平均剪應力τm可由式(10)獲得：

基于摩爾庫倫強度準則判別黏結單元的法向拉伸破壞和切向剪切破壞。若黏結單元平均正應力超過單元抗拉強度值，黏結發生拉伸斷裂，黏結單元失效，兩側塊體由黏結狀態變為純接觸狀態。

式中，Tb為黏結單元的抗拉強度。

若黏結單元平均剪應力超過剪切極限值，則黏結單元發生剪切斷裂，黏結單元失效，同時兩側塊體由黏結狀態變為純接觸狀態。

式中，cb為邊-邊黏結單元的黏聚力；φb為邊-邊黏結單元的內摩擦角。

1.2.3分段式黏結模型的力學求解

黏結單元的本構模型采用罰函數方法計算，邊-邊黏結單元的法向約束和切向約束的能量泛函可表示為

式中，p0、q0為時步初始時刻分段si中線點的空間位置向量；ni和ti為當前分段的法向和切向單位向量；up和uq為當前時步p點和q點的位移增量向量，可由式u=Td求得。

以第si分段為例，根據能量泛函最小化原理，法向約束能量泛函對未知量D的二階導數項和一階導數項分別以子矩陣和子向量的形式加入總剛度矩陣和總廣義力向量：

同理，切向約束能量泛函對未知量D的二階導數項和一階導數項也分別加入總剛度矩陣和總廣義力向量：

黏結單元各分段按上述公式求解并疊加可得到整個邊-邊黏結單元模型在DDA 計算框架中的罰函數求解公式。

2 模型連續變形與應力校驗

準確模擬材料在連續變形階段的變形與應力是模擬巖石變形和破壞的基礎。對于分段黏結模型，需要確定其幾何參數(分段數目、黏結單元尺寸)、變形參數(黏結剛度比、切向剛度比)、強度參數(抗拉強度、黏聚力、內摩擦角)。

采用簡支梁中心受壓、單軸壓縮和巴西劈裂這3個實驗，通過對比模型的解析解與模擬結果，對分段式黏結模型的幾何參數、變形參數進行校驗并提出合適的取值范圍或推薦值。

2.1 簡支梁中心受壓變形驗證

本文選取簡支梁模型如圖2 所示，長度l=1 m，高度H=0.1 m，寬度取單位值b=1 m。簡支梁模型的密度ρ=2 500 kg/m3，彈性模量E=50 GPa，泊松比μ=0.2，抗拉強度T=12 MPa，黏聚力c=63 MPa，內摩擦角φ=42°。在梁的左端約束其水平位移和垂直位移，右端只約束其垂直位移，分別模擬簡支梁中的鉸支座和滾動支座。在梁的中間處施加Fp=20 kN 的集中荷載，梁中軸線中心點處垂直變形w的解析解為

圖2 不同塊體尺寸的簡支梁模型Fig.2 Illustrations of simplely supported beam models with different block sizes

其中，I為梁的截面慣矩，計算公式為

可求得梁中心點處垂直位移的理論值為-0.1 mm。

在本算例模擬中，考慮了4 種塊體(黏結單元)尺寸(15.0、10.0、7.5、5.0 mm)以及6 種黏結剛度比(5、10、25、50、75、100)，黏結單元分段數目取2，切向黏結剛度比取1.0，設置梁的中軸線中心點為監測點。監測點得到的垂直位移如圖3(a)所示，模擬結果與理論結果的相對誤差如圖3(b)所示，可以看出，隨著黏結剛度比增加，監測點的位移逐漸接近理論值，且隨著劃分塊體尺寸的增大，監測點的位移越來越接近理論值。當黏結剛度比大于25 時，所有模型的相對誤差都達到了10%以下，達到可以接受的精度。另外，控制其他變量相同，取黏結單元分段數為2、3、5 時的計算結果表明結果相對偏差遠小于1%，取分段數為2 即可滿足要求。

圖3 簡支梁監測點的垂直位移與相對誤差Fig.3 Simulated vertical displacement and relative error of the track point in simply supported beam example

并且當黏結剛度比比大于50 后，相對誤差的變化趨于穩定，且削弱了塊體尺寸的影響。故建議黏結剛度比應在[50,100]為宜。選取黏結剛度比為50，塊體尺寸為7.5 mm 的簡支梁模型，在梁中心處分別施加1Fp至5Fp的力，在線彈性階段模型的相對誤差保持在4.09%，不發生顯著變化，表明未發生破壞階段變形計算的穩定性。

2.2 試樣單軸壓縮應力驗證

巖石單軸壓縮試驗是測定巖石強度的基本試驗，其測定的單軸抗壓強度在地下工程中應用廣泛，具有重要的工程意義。筆者建立長100 mm、寬50 mm 的矩形作為標準試樣單軸壓縮二維模型，模型的物理力學參數與簡支梁一致。如圖4 所示，在模型上邊界施加0.1 MPa 的均布載荷，下邊界固定，模型中心處設置一監測點。對于均勻試樣，在軸向應力應在垂直方向上均勻分布。

圖4 不同塊體尺寸單軸壓縮模型Fig.4 Illustrations of uniaxial compression models with different block sizes

在本算例模擬中，考慮了4 種塊體尺寸(10.0、7.5、5.0、2.5 mm)以及6 種黏結剛度比(5、10、25、50、75、100)，黏結單元分段數目取2，切向黏結剛度比取1.0，模型中心的應力如圖5(a)所示，模擬結果與理論結果的相對誤差如圖5(b)所示。由圖5 可知，所有模型的相對誤差均低于2%，表明分段式黏結模型在連續變形階段的應力分布精度足夠高。且隨著黏結剛度比的增加，中心點的軸向應力越接近理論值，當黏結剛度比大于25 時，相對誤差低于0.5%?？刂破渌兞肯嗤?，取黏結單元分段數為2、3、5 時的計算結果表明結果相對偏差遠小于1%，取分段數為2 即可滿足要求。選取黏結剛度比為50，塊體尺寸為5 mm 的單軸壓縮模型，在模型上邊界施加0.1～1.0 MPa 的分布力，在線彈性階段模型的相對誤差保持在0.13%，不發生顯著變化，表明未發生破壞階段應力計算的穩定性。

圖5 單軸壓縮監測點的軸向應力與相對誤差Fig.5 Simulated stress and relative error of the track point in uniaxial compression test

結合簡支梁中心受壓與單軸壓縮的模擬結果進行分析可知，在幾何參數方面，塊體尺寸、黏結單元數目對于算例中連續變形、及應力計算精度的影響不顯著；在變形參數方面，當黏結剛度比≥50 時，可以保證變形與應力的相對誤差都滿足精度要求，同時減小了塊體尺寸的影響，因此在后續模擬中將采用黏結剛度比50 進行實驗；而剪切剛度比值設為1.0，可滿足算例中對于計算精度的要求。對比原有黏結模型[18-19]，分段式黏結模型中黏結剛度比可以克服對于不同塊體尺寸下取不同剛度的復雜取值問題，根據經驗在[50,100]內可獲得合適的精度要求。

2.3 巴西圓盤應力分布驗證

巴西劈裂試驗是一種廣泛使用的間接測量巖石抗拉強度的方法。本例中采用平臺巴西劈裂試驗，巴西圓盤模型如圖6 所示，半徑R=0.025 m，厚度取單位值w=1 m，物理力學參數設置與簡支梁一致，黏結剛度比取50，塊體尺寸取1～2 mm。模擬時將圓盤的下邊界固定，在上邊界加載平臺上假定試驗機施加了一集中力F=10 kN，平臺加載角2α=12°，則加載平臺上的均布荷載P為

圖6 巴西劈裂模型Fig.6 An illustration of the Brazilian splitting model

式中，α為弧度。

平臺巴西劈裂中圓盤各個點的x方向應力與y方向應力的解析解參考黃耀光等[26]的計算公式。為驗證應力計算的準確性，在圓盤中心沿垂直方向均勻設置9 個監測點。監測點x方向應力與y方向應力的模擬結果與解析解的對比如圖7(a)所示，相對誤差如圖7(b)所示。由圖7 可以看出，監測點應力的模擬結果與解析解基本吻合。y方向應力的相對誤差均低于10%，x方向應力的相對誤差表現為越遠離端面精度越高，在靠近端面時會表現出較大的誤差，由于數值求解時模型與邊界條件與解析求解時并不完全一致，其結果符合圣維南原理。算例進一步驗證了分段式黏結模型的準確性。

圖7 巴西劈裂監測點應力與相對誤差Fig.7 Stress measurement of track point in Brazilian splitting test and relative errors

3 砂巖準靜態變形及強度特性

在天然的巖體中，總是存在著許多裂隙、節理和斷層等不連續面，這些缺陷很大程度的影響了巖體的強度特性以及斷裂模式，所以探究巖石在含預制缺陷時的強度特性與裂紋擴展規律對地下工程圍巖破壞問題有重要價值。筆者選取煤礦開采等地下工程中最為常見的巖石種類之一——砂巖為研究對象，設置了2 組模擬實驗，探究了缺陷對砂巖強度特性和變形的影響，表1 總結了2 組實驗所使用到的參數。

表1 模擬參數Table 1 Simulation parameters

3.1 單軸壓縮強度實驗模擬

選擇對楊圣奇等[27]的含孔洞的砂巖單軸壓縮實驗進行模擬，探究孔洞對單軸壓縮強度特性以及裂紋萌生的影響。該砂巖試樣致密均勻，平均密度為2 620 kg/m3，高度為120 mm，寬度為60 mm，厚度為30 mm。選擇長度為120 mm 寬度為60 mm 的矩形截面作為單軸壓縮模擬二維模型，分別建立孔洞直徑d=0(完整巖樣)、5、10、15 mm 的模型，在孔洞周圍劃分更小的塊體以更好的捕捉裂紋萌生過程，如圖8 所示。單軸壓縮模擬過程中，模型的下邊界固定，上邊界以0.5 mm/s 的速度垂直向下加載。

圖8 含孔洞巖樣數值模型Fig.8 Numerical models of rock samples with holes

進行數值模擬實驗前，需對模型的子塊體彈性參數(彈性模量、泊松比)、黏結剛度參數(黏結剛度比，經驗取值50；切向與法向剛度比，取值1.0)、黏結強度參數(抗拉強度、黏聚力、內摩擦角)進行校準，采用“試錯法”對以上參數進行標定。

筆者通過進行大量的標準巖樣單軸壓縮和巴西劈裂試驗模擬尋找參數規律，確定的標定步驟為：①通過模型的單軸壓縮試驗來校準彈性模量、泊松比，根據巖樣的真實參數進行賦值并進行±5%范圍內調整；②通過模型的巴西劈裂試驗來校準模型的抗拉強度，賦予模型真實巖樣的抗拉強度并進行±5%范圍內調整；③通過單軸抗壓強度來校驗黏聚力與內摩擦角，對比強度和破壞模式確定最終取值。結合該實驗完整巖樣單軸壓縮曲線，標定結果見表1，單軸壓縮模擬的應力應變響應曲線與實驗曲線的對比如圖9 所示，可以看出2 條曲線峰值強度與總體趨勢基本吻合。

圖9 完整巖樣實驗與模擬應力應變曲線對比Fig.9 Comparisons of experimental and simulated stress-strain curves for uniaxial compression of intact rock sample

使用同一套參數進行含孔洞砂巖的模擬實驗，得到的應力應變響應曲線與裂紋萌生狀態如圖10 所示。由圖10 可知，在含孔洞模型中，拉伸裂紋最先在孔洞中心附近上下同時萌生，同時，強度特性表現為隨孔洞直徑的增加，單軸壓縮峰值強度與峰值應變均出現明顯的下降，這些現象均與實驗室實驗觀察到的一致。對比圖9、圖10 可以看出，使用分段式黏結模型模擬可以很好的捕捉到試樣的抗壓強度與裂紋萌生點，最大的差別在于加載的初始階段，這是由于實驗室實驗所使用的試樣內部存在一些微裂隙，在位移加載初期表現為壓實階段，而模擬時是一個完整的連續體，并沒有這個階段。

圖10 應力應變曲線對比Fig.10 Comparisons of stress-strain curves

3.2 巴西圓盤劈裂實驗模擬

本算例探究砂巖巴西圓盤在單裂紋、多裂紋條件下的裂紋擴展模式，對ZHAO 等[28]的實驗進行了數值模擬，使用的參數取值見表1。

單裂紋巴西圓盤的模型如圖11(a)所示，圓盤的直徑為62 mm，裂隙的長度為20 mm，寬度為1 mm。裂隙與水平方向分別呈0°、30°、60°、90°，塊體尺寸為1～3 mm，在有裂紋的地方劃分更小的子塊體以更好地捕捉裂紋萌生和裂紋擴展的過程。將模型的下邊界固定，在上邊界中以0.5 mm/s 的速度垂直向下進行位移加載模擬試樣巴西劈裂過程。由圖11(b)、(c)可知，當裂紋傾角為0°時，裂紋從預制裂紋的中部上下同時萌生，并沿著垂直方向擴展，直至貫穿。對于其他3 種裂紋傾角模型，裂紋都是從預制裂紋的尖端開始萌生，并沿著加載方向繼續擴展，直至完全貫穿。

圖11 單裂紋巴西圓盤模型及其裂紋萌生和破壞模式Fig.11 Brazilian disk models with single crack and their failure mode

含有多條裂紋的巴西圓盤模型如圖12(a)所示，圓盤的直徑、裂紋的大小、子塊體剖分策略、加載方式均與單裂紋模型一致，多裂紋構成的裂紋系統有如圖①、②、③、④四種形式。由圖12(b)、(c)可以看出，對于模型①和模型②，裂紋從預制裂紋L1 和L2 的尖端同時萌生，靠近試樣端面的2 個尖端從裂紋萌生開始，不斷沿著加載方向擴展直至貫穿，另外2 個尖端裂紋則不斷沿著靠近模型中軸線的方向擴展。對于模型③，裂紋從預制裂紋L2 的上下尖端同時萌生，并沿著加載方向擴展直至貫通破壞，裂紋擴展模式與單裂紋巴西劈裂類似。對于模型④，裂紋最先從預制裂紋L1 和L3 的上下中部萌生，緊接著預制裂紋L2 的上下尖端也同時出現裂紋，并且逐漸向預制裂紋L1和L3 合并，最后一起上下貫穿。

圖12 含有多條裂紋的巴西圓盤及其裂紋萌生和破壞模式Fig.12 Brazilian discs containing multiple cracks and their pattern of crack initiation and crack propagation

將模擬結果與室內實驗結果的進行了比較可以發現，在預制單裂紋和預制多裂紋的巴西圓盤在裂紋萌生、裂紋擴展以及最終的破壞模式上，取得了很高的一致性。表明了分布式黏結模型在模擬裂紋行為方面能夠較為準確的預測裂紋萌生位置與擴展軌跡。

4 結論

(1) 提出一種分段式邊-邊黏結模型對黏結效應進行等效計算，推導了分段式式黏結模型的求解公式，擴充了基于子塊體剖分的非連續變形分析計算框架，提高了黏結模型在面對復雜黏結本構行為的適用性。

(2) 通過簡支梁中心受壓、單軸壓縮、巴西圓盤受壓3 個算例驗證了邊-邊黏結模型在連續變形和應力計算中的準確性，對于黏結模型中的重要參數如黏結剛度比的取值提出了合理建議范圍。

(3) 針對煤礦頂底板常遇到的砂巖，進行了含孔洞砂巖試樣、含預制裂紋圓盤試樣的加載模擬，發現采用分段式黏結模型的DDA 方法可以準確捕捉裂紋的萌生位置、裂紋擴展路徑和多裂紋作用下的砂巖破壞過程，驗證了模型在砂巖破壞分析中的良好適用性。