分段黏結(jié)非連續(xù)變形分析方法及其在砂巖破裂分析中的應(yīng)用

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

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

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

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

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

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

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

1 基于分段式黏結(jié)模型的DDA 方法

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

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

1.1.1數(shù)學原理與控制方程

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

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

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

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

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

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

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

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

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

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

1.2 分段式塊體黏結(jié)模型

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

1.2.1分段式黏結(jié)模型概述

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

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

1.2.2黏結(jié)本構(gòu)關(guān)系及分段式表征

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

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

基于摩爾庫倫強度準則判別黏結(jié)單元的法向拉伸破壞和切向剪切破壞。若黏結(jié)單元平均正應(yīng)力超過單元抗拉強度值，黏結(jié)發(fā)生拉伸斷裂，黏結(jié)單元失效，兩側(cè)塊體由黏結(jié)狀態(tài)變?yōu)榧兘佑|狀態(tài)。

式中，Tb為黏結(jié)單元的抗拉強度。

若黏結(jié)單元平均剪應(yīng)力超過剪切極限值，則黏結(jié)單元發(fā)生剪切斷裂，黏結(jié)單元失效，同時兩側(cè)塊體由黏結(jié)狀態(tài)變?yōu)榧兘佑|狀態(tài)。

式中，cb為邊-邊黏結(jié)單元的黏聚力；φb為邊-邊黏結(jié)單元的內(nèi)摩擦角。

1.2.3分段式黏結(jié)模型的力學求解

黏結(jié)單元的本構(gòu)模型采用罰函數(shù)方法計算，邊-邊黏結(jié)單元的法向約束和切向約束的能量泛函可表示為

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

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

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

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

2 模型連續(xù)變形與應(yīng)力校驗

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

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

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，內(nèi)摩擦角φ=42°。在梁的左端約束其水平位移和垂直位移，右端只約束其垂直位移，分別模擬簡支梁中的鉸支座和滾動支座。在梁的中間處施加Fp=20 kN 的集中荷載，梁中軸線中心點處垂直變形w的解析解為

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

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

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

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

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

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

2.2 試樣單軸壓縮應(yīng)力驗證

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

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

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

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

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

2.3 巴西圓盤應(yīng)力分布驗證

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

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

式中，α為弧度。

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

圖7 巴西劈裂監(jiān)測點應(yīng)力與相對誤差Fig.7 Stress measurement of track point in Brazilian splitting test and relative errors

3 砂巖準靜態(tài)變形及強度特性

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

表1 模擬參數(shù)Table 1 Simulation parameters

3.1 單軸壓縮強度實驗?zāi)M

選擇對楊圣奇等[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 含孔洞巖樣數(shù)值模型Fig.8 Numerical models of rock samples with holes

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

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

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

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

圖10 應(yīng)力應(yīng)變曲線對比Fig.10 Comparisons of stress-strain curves

3.2 巴西圓盤劈裂實驗?zāi)M

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

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

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

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

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

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

4 結(jié)論

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

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

(3) 針對煤礦頂?shù)装宄Ｓ龅降纳皫r，進行了含孔洞砂巖試樣、含預(yù)制裂紋圓盤試樣的加載模擬，發(fā)現(xiàn)采用分段式黏結(jié)模型的DDA 方法可以準確捕捉裂紋的萌生位置、裂紋擴展路徑和多裂紋作用下的砂巖破壞過程，驗證了模型在砂巖破壞分析中的良好適用性。