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

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

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

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

數(shù)值模擬技術(shù)由于其靈活性、定量性、過程可視性等特點(diǎn)，也是巖石多尺度斷裂分析研究的重要手段。分析巖石斷裂過程的數(shù)值方法，通?？蓜澐譃榛谶B續(xù)介質(zhì)假定的方法、基于非連續(xù)介質(zhì)假定的方法、以及混合連續(xù)與非連續(xù)框架進(jìn)行求解的方法?；谶B續(xù)介質(zhì)假定的方法包括擴(kuò)展有限元(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ù)值計(jì)算方法可以考慮較大尺度斷層、節(jié)理、層理等不連續(xù)面的張開、閉合、滑移，也可以模擬較小尺度的巖石斷裂問題，其中DDA 方法在接觸的數(shù)學(xué)計(jì)算、塊體的積分、計(jì)算收斂性方面具有嚴(yán)密的數(shù)學(xué)基礎(chǔ)。不同于采用顯式時(shí)間積分計(jì)算的DEM，DDA 可視為基于Newmark 隱式積分和接觸開閉迭代求解的隱式非連續(xù)計(jì)算方法，在離散塊體系統(tǒng)大變形、大位移計(jì)算中可以保持時(shí)間步取值的無條件穩(wěn)定。DDA 方法最初用于解決裂隙切割的塊體系統(tǒng)的靜態(tài)穩(wěn)定性及動(dòng)態(tài)響應(yīng)問題。KE[17]提出采用子塊體劃分DDA 方法模擬裂隙巖體變形問題，焦玉勇等[18]則通過引入節(jié)理約束的三角形子塊體剖分及黏結(jié)單元來模擬節(jié)理巖體的漸進(jìn)破裂過程，甯尤軍等[19]和ZHENG 等[20]相繼對(duì)黏結(jié)單元的破壞準(zhǔn)則、單元積分模型進(jìn)行修正和優(yōu)化。面對(duì)深部開采中與巖層破斷、非連續(xù)大變形相關(guān)的巖石力學(xué)問題，DDA 方法具有天然的優(yōu)勢。許多學(xué)者利用DDA 方法研究了深部沿空掘巷過程[21]、煤層開采導(dǎo)致的地表沉降[22]、煤礦巖體大變形和開挖過程[23]、礦山陷落柱[24]等問題。

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

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

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

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

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

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

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

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

式中，x0、y0為塊體的形心點(diǎn)坐標(biāo)分量。

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

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

1.1.2時(shí)間積分方案與隱式求解模式

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.1 簡支梁中心受壓變形驗(yàn)證

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

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

其中，I為梁的截面慣矩，計(jì)算公式為

可求得梁中心點(diǎn)處垂直位移的理論值為-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è)置梁的中軸線中心點(diǎn)為監(jiān)測點(diǎn)。監(jiān)測點(diǎn)得到的垂直位移如圖3(a)所示，模擬結(jié)果與理論結(jié)果的相對(duì)誤差如圖3(b)所示，可以看出，隨著黏結(jié)剛度比增加，監(jiān)測點(diǎn)的位移逐漸接近理論值，且隨著劃分塊體尺寸的增大，監(jiān)測點(diǎn)的位移越來越接近理論值。當(dāng)黏結(jié)剛度比大于25 時(shí)，所有模型的相對(duì)誤差都達(dá)到了10%以下，達(dá)到可以接受的精度。另外，控制其他變量相同，取黏結(jié)單元分段數(shù)為2、3、5 時(shí)的計(jì)算結(jié)果表明結(jié)果相對(duì)偏差遠(yuǎn)小于1%，取分段數(shù)為2 即可滿足要求。

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

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

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

巖石單軸壓縮試驗(yàn)是測定巖石強(qiáng)度的基本試驗(yàn)，其測定的單軸抗壓強(qiáng)度在地下工程中應(yīng)用廣泛，具有重要的工程意義。筆者建立長100 mm、寬50 mm 的矩形作為標(biāo)準(zhǔn)試樣單軸壓縮二維模型，模型的物理力學(xué)參數(shù)與簡支梁一致。如圖4 所示，在模型上邊界施加0.1 MPa 的均布載荷，下邊界固定，模型中心處設(shè)置一監(jiān)測點(diǎn)。對(duì)于均勻試樣，在軸向應(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é)果的相對(duì)誤差如圖5(b)所示。由圖5 可知，所有模型的相對(duì)誤差均低于2%，表明分段式黏結(jié)模型在連續(xù)變形階段的應(yīng)力分布精度足夠高。且隨著黏結(jié)剛度比的增加，中心點(diǎn)的軸向應(yīng)力越接近理論值，當(dāng)黏結(jié)剛度比大于25 時(shí)，相對(duì)誤差低于0.5%。控制其他變量相同，取黏結(jié)單元分段數(shù)為2、3、5 時(shí)的計(jì)算結(jié)果表明結(jié)果相對(duì)偏差遠(yuǎn)小于1%，取分段數(shù)為2 即可滿足要求。選取黏結(jié)剛度比為50，塊體尺寸為5 mm 的單軸壓縮模型，在模型上邊界施加0.1～1.0 MPa 的分布力，在線彈性階段模型的相對(duì)誤差保持在0.13%，不發(fā)生顯著變化，表明未發(fā)生破壞階段應(yīng)力計(jì)算的穩(wěn)定性。

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

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

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

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

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

式中，α為弧度。

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

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

3 砂巖準(zhǔn)靜態(tài)變形及強(qiáng)度特性

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

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

3.1 單軸壓縮強(qiáng)度實(shí)驗(yàn)?zāi)M

選擇對(duì)楊圣奇等[27]的含孔洞的砂巖單軸壓縮實(shí)驗(yàn)進(jìn)行模擬，探究孔洞對(duì)單軸壓縮強(qiáng)度特性以及裂紋萌生的影響。該砂巖試樣致密均勻，平均密度為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

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

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

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

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

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

3.2 巴西圓盤劈裂實(shí)驗(yàn)?zāi)M

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

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

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

含有多條裂紋的巴西圓盤模型如圖12(a)所示，圓盤的直徑、裂紋的大小、子塊體剖分策略、加載方式均與單裂紋模型一致，多裂紋構(gòu)成的裂紋系統(tǒng)有如圖①、②、③、④四種形式。由圖12(b)、(c)可以看出，對(duì)于模型①和模型②，裂紋從預(yù)制裂紋L1 和L2 的尖端同時(shí)萌生，靠近試樣端面的2 個(gè)尖端從裂紋萌生開始，不斷沿著加載方向擴(kuò)展直至貫穿，另外2 個(gè)尖端裂紋則不斷沿著靠近模型中軸線的方向擴(kuò)展。對(duì)于模型③，裂紋從預(yù)制裂紋L2 的上下尖端同時(shí)萌生，并沿著加載方向擴(kuò)展直至貫通破壞，裂紋擴(kuò)展模式與單裂紋巴西劈裂類似。對(duì)于模型④，裂紋最先從預(yù)制裂紋L1 和L3 的上下中部萌生，緊接著預(yù)制裂紋L2 的上下尖端也同時(shí)出現(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)實(shí)驗(yàn)結(jié)果的進(jìn)行了比較可以發(fā)現(xiàn)，在預(yù)制單裂紋和預(yù)制多裂紋的巴西圓盤在裂紋萌生、裂紋擴(kuò)展以及最終的破壞模式上，取得了很高的一致性。表明了分布式黏結(jié)模型在模擬裂紋行為方面能夠較為準(zhǔn)確的預(yù)測裂紋萌生位置與擴(kuò)展軌跡。

4 結(jié)論

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

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

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