不同圍壓條件下砂巖動態剪切特性數值模擬分析

朱看遠，趙格立，徐 穎，江麗媛

(1.天津大學水利工程仿真與安全國家重點實驗室，天津 300072；2.天津大學建筑工程學院，天津 300350)

巖石的剪切破壞是工程建設中常見的破壞模式，因此準確測定巖石的剪切強度是至關重要的[1-3]。室內實驗中，直剪和沖剪實驗是測定靜態剪切強度常用的實驗方法[1]，由于沖剪實驗具有試樣制作簡單，適用于高強度固體等優勢，因此被廣泛用于巖石的剪切強度研究[4]。沖剪實驗是用來測量圓盤狀巖石的剪切強度而專門設計的，最早采用的是圓形沖頭，但彎曲應力會對巖石造成額外損傷使得測定結果不精確[5-6]。為減少彎曲應力的影響，Stacey[7]采用方形沖頭進行沖剪實驗，從而提高了剪切強度測定的精確性。該方法經過多年的發展和優化，已經成為測量脆性材料靜態剪切強度的常用方法[8-9]。在工程實踐中，圍巖易受到工程爆破、地震等動態載荷的影響，因此對巖石動態剪切特性的研究尤為重要[10-11]。基于此，Huang等[1]通過改進的分離式霍普金森壓桿(SHPB)裝置，實現了對無圍壓條件下砂巖的動態剪切強度特性的研究，獲得了砂巖剪切強度與加載率的關系。Du等[12]采用SHPB對傾斜試樣進行動態沖擊，計算了巖石動態剪切強度的率相關性。為了考慮深部圍壓的影響，Xu等[13]采用圍壓缸與SHPB裝置結合，實現了0～20 MPa圍壓條件下砂巖動態剪切強度特性的研究。綜上所述，目前對巖石剪切強度的研究主要集中在靜態實驗條件下，對于巖石動態剪切特性的研究也進行不少的嘗試，但大多集中在無圍壓或圍壓較低的實驗條件下。這與深部巖石所處的地應力條件相差甚大。因此，有必要對深部高地應力環境下巖石動態剪切特性進行研究。

數值模擬能夠實現室內實驗不易達到的條件，對實驗進行驗證與補充，因此被廣泛用于巖石力學特性研究。張酈榮等[14]基于有限元方法，模擬了松軟土的剪切特性，數值結果和實驗呈現出一致性；楊開新等[15]基于三維離散元模擬了圍壓狀態下的砂巖的靜態剪切破裂特性；蔡國慶等[16]建立三維離散元模型，模擬了不同圍壓與含水率下的土體的剪切特性；梁冰寒等[17]基于離散元軟件，模擬了圍壓對堆石料剪切破壞帶的影響。

強度準則的選取是巖石剪切破壞模擬的重要步驟。Mohr-Coulomb準則(M-C)是最常用的強度準則，但M-C屈服面在偏平面上為六角形，塑性分析時，因塑性流動方向不唯一會引起收斂困難[18]。因此學者對于模型準則的選取進行相關研究，岑奪豐等[19]采用Hoek-Brown強度準則替代M-C強度準則，有效地模擬了砂巖的靜態拉剪破壞特征。在M-C準則的基礎上提出的Drucker-Prager(D-P)準則，避免了塑性分析時收斂困難的問題，并且能較好地反映巖石的強度特性，但其參數的選取會影響結果的準確性[20]。劉金龍[21]、郝明輝等[22]基于數值軟件與理論分析，對D-P準則的相關特性參數以及模型的適用性展開了討論，得到了適用于剪切破壞擬合的參數值選取范圍?；谠撊≈捣秶?，張強等[23]采用ABAQUS軟件，將D-P模型應用于鉆頭破巖數值模擬，得到了影響鉆頭破巖效果的因素。Genikomsou[24]基于ABAQUS三維建模，采用D-P準則對不同形狀的混凝土板進行沖剪效應模擬，并和實驗數據對比，較好地重現了準靜態下混凝土的沖剪破壞響應?？梢姡珼-P塑性模型在模擬巖石材料破壞時，也具有一定的優勢。

綜上，目前對巖石的剪切特性的數值研究以靜態居多，而砂巖的動態剪切特性的數值研究較少，圍壓條件下的數值計算更是鮮為報道。基于此，本文采用ABAQUS有限元軟件和D-P塑性模型，模擬不同圍壓下砂巖的動態剪切特性。以引文[13]中砂巖試樣在不同圍壓下的實驗數據為基礎，創建并驗證了模型的合理性和適用性。并對高圍壓下的砂巖動態沖剪特性進行模擬，研究其動態沖剪特性，以期為深部巖石動態剪切特性及工程結構設計研究提供參考。

1 創建模型

在ABAQUS中，模擬材料的剪切失效主要有兩個階段：損傷開裂階段和損傷演化階段。損傷開裂階段由斷裂應變來控制，還可設置應變率與斷裂應變的關系。損傷演化是指當單元達到某一指標后，剛度降低的過程。當剛度降低至0時，產生裂紋。一般選用斷裂能作為判斷依據，對于巖石或者混凝土，斷裂能一般設定為40～120 N/m[25]，因此本文選取該范圍作為設定。但剪切損傷準則并不能單獨使用，因其須先對受損單元進行斷裂應變的判別，才能進行損傷演化。而斷裂應變的判別標準為最大塑性主應變，因此要與塑性模型配合使用[26]。

1.1 Druker-Prager塑性模型

本文模擬的是砂巖的動態沖剪實驗，采用ABAQUS/Explicit分析模塊，D-P塑性模型的屈服面在子午面上所提供的線性、雙曲線和指數函數模型中，只有線性D-P模型才能適用于Explicit分析模塊[25-26]。線性D-P模型在子午面上的屈服軌跡如圖1所示。

屈服函數為

F=t-ptanβ-d=0

(1)

線性D-P模型在π平面上的屈服面不是圓形，非圓形的屈服面可以真實地反映不同的三軸拉伸和壓縮屈服強度，π平面上的塑性流動以及內摩擦角和剪脹角。線性D-P模型在π平面上的屈服面如圖2所示。

d為屈服面在p-t應力空間上的截距，是另一種形式的黏聚力，由于本實驗是剪切實驗，所以可根據剪切強度τ來定義，如下式所示：

(2)

在D-P線性模型中塑性流動規則的表達式為

G=t-tanψ

(3)

式中：G為塑性流動勢；ψ為p-t應力空間上的剪脹角，當ψ=0，塑性變形時的材料體積不發生變化。

1.2 率相關性

在動態實驗中，應變率對材料的強度有很大的影響，稱之為率相關性。線性D-P模型的硬化規律支持材料屈服強度的率相關性，為了反映應變率對巖石強度的影響，需設置等效應力的率相關參數。通過調試線性D-P模型中的等效應力，得到的透射波的波形數據與實驗數據相同，計算巖石模型的剪切應變率，從而得到剪切應變率與等效應力之間的關系，剪切應變率與等效應力之間的關系用對數函數進行擬合。

(4)

擬合得a=-156.266 37；b=-36.610 67；c則隨著圍壓的增大呈線性變化。

c=9.671σ+111.618 33

(5)

式中：σ為圍壓，MPa。

通過以上的分析，建立了圍壓-剪切應變率-等效應力的關系。巖石在不同的圍壓條件下，計算得到不同剪切應變率下所對應的等效應力，再由剪切損傷模型對發生塑性變形的單元進行刪除，產生裂紋，從而達到沖剪的目的，由透射桿的波形數據可分析得到巖石的動態剪切強度。

1.3 沖剪實驗模型

本文以引文[13]中使用的動態加載設備及試樣作為建模基礎。其所采用的SHPB加載系統主要由撞擊桿、入射桿、透射桿和圍壓油缸等裝置組成。撞擊桿、入射桿和透射桿的材質為鋼，入射桿的應變片距離試樣1.040 m，透射桿的應變片距離試樣0.503 m。實驗裝置[13]的各參數如表1所示。應力波由入射桿和透射桿上的應變片捕捉。入射波εi，反射波εr，透射波εt。通過這三組波形，可以計算試樣兩端的動態力[28]：

p1(t)=EA[εi(t)+εr(t)]

(6)

p2(t)=EA[εt(t)]

(7)

式中：E為桿的彈性模量；A為桿的橫截面面積。當試樣處于動態力平衡狀態時(p1=p2)，動態實驗中的慣性效應可以忽略[29]，此時可以計算試樣所受的動態剪切力：

(8)

式中：p(t)=p1(t)=p2(t)；D為入射桿的直徑；B為試樣的厚度。τ(t)的最大值則為動態剪切強度。

在ABAQUS創建入射桿、透射桿、后支撐和巖樣并對其進行裝配。沖剪模擬需要模擬材料的失效問題，顯示分析中可以很好地解決收斂問題，因此選用動態顯示分析方法作為分析步。根據波速以及模型尺寸，分析時長為700 μs。為了體現單元的失效和刪除，對輸出結果的狀態(STATUS)進行設置。

SHPB沖剪實驗中接觸部分主要為兩部分，分別為入射桿與巖樣的接觸面以及透射桿支撐架與巖樣的接觸面。接觸屬性設置為法向硬接觸和切向無摩擦，法向硬接觸即為兩部件之間的網格不能相互嵌入，切向無摩擦即為兩部件在切向光滑接觸。此外，桿件的剛度比巖樣大得多，因此將桿件端面設置為主面、巖樣為從面。

模擬中直接在入射桿的入射端面施加隨時間變化的脈沖荷載，即在入射桿端面導入文獻[13]的原始波形數據。為使脈沖荷載更平滑，選擇幅值中的平滑分析曲線，作為荷載隨時間變化的曲線。

在0 MPa圍壓條件下對巖樣不設置圍壓，而對于有圍壓的條件，應對巖樣施加預定的圍壓。動態分析前采用靜力步分析巖樣在圍壓加載條件下的應力場，并將結果導入到動態分析模型中。模型網格劃分及其剖視圖，如圖3所示。

2 數值模擬方法的可靠性驗證

2.1 力平衡與加載率的對比驗證

為了驗證模型的合理性與適用性，在0～20 MPa圍壓下，采用與實驗相同的脈沖荷載來模擬砂巖的沖剪實驗，將所得結果與相同條件下的實驗數據進行對比。在SHPB實驗中，進行動態力平衡，用于消除慣性效應[2, 30]，因此對沖剪模型進行力平衡模擬，試驗與模擬力平衡對比如圖4、圖5所示。由力平衡圖可以看出，無圍壓條件下實驗與模擬到達峰值點的時間(實驗87 μs,模擬83 μs)，以及對應的峰值大小(實驗35 kN，數值38 kN)(見圖4)基本一致。20 MPa圍壓條件下，二者到達峰值點的時間(實驗90 μs,模擬85 μs)，以及對應的峰值大小(實驗80 kN，數值82 kN)(見圖5)也基本吻合。模擬的力平衡效果與實驗中力平衡相差不大，這表明采用D-P塑性模型在力平衡的模擬上是可行的。

模擬力平衡后，便可根據式(8)計算試樣的沖剪應力。數值模擬的結果也同樣出現了明顯的峰值點，對應著動態沖剪強度與剪切破壞的開始。沖剪實驗的加載率由破壞開始前剪應力的斜率確定[31]。本文通過最小二乘法對所得到的剪應力時程曲線擬合，得到了加載率的數值。并將模擬的結果與實驗進行對比，結果如圖6、圖7所示。無圍壓條件下，實驗與模擬到達峰值的時間(實驗90 μs，模擬88 μs)，對應的峰值強度(實驗36 MPa，模擬37.5 MPa)，以及由最小二乘法擬合得到的加載率(實驗779.3 GPa/s，模擬779.2 GPa/s)，在數值上基本一致(見圖6)。在圍壓20 MPa條件下，實驗與模擬到達峰值的時間(實驗90 μs，模擬95 μs)，對應的峰值強度(實驗92 MPa，模擬93.5 MPa)，以及由最小二乘法擬合得到的加載率(實驗3 087.5 GPa/s，模擬3 083.5 GPa/s)(見圖7)，可見具有較高的吻合度。通過對0 MPa和20 MPa圍壓下所得到的剪應力時程曲線分析可知，實驗結果與模擬結果在加載率等方面基本一致，即采用D-P塑性模型進行數值模擬能夠重現沖剪實驗的剪應力時程曲線，得到與實驗結果基本一致的加載率。

2.2 動態剪切強度的對比驗證

通過對動態力平衡和加載率的對比分析可知，采用D-P塑性模型能夠實現動態力平衡且加載率與實驗呈現較高的吻合度。本文將模擬的不同圍壓和加載率下巖樣的動態沖剪實驗的數據，按照式(8)進行了計算。將得到的結果進行了擬合和實驗結果進行對比驗證，如圖8所示。結果表明：在圍壓一定的條件下，隨著加載率的升高其動態剪切強度逐漸增大；在加載率相差不大的條件下，動態剪切強度會隨著圍壓的增大而逐漸增大。文獻[13]對所得的實驗結果進行擬合，其數據較為均勻地分布在趨勢線兩側。將模擬所得的剪切強度結果與實驗對比，模擬計算所得的數據大致分布于文獻[13]所擬合的直線兩側，與實驗的結果呈現出較高的一致性。這表明D-P塑性模型適用于動態沖剪實驗的模擬，能夠重現動態沖剪實驗，表現出巖樣沖剪實驗的剪切強度特性。

2.3 巖樣破裂的對比驗證

試樣的破壞模式能直觀的展示其破裂機制，引文[13]對沖剪破壞的試樣進行CT掃描。為了進一步驗證模型的合理性，本文將數值模擬巖樣的破壞結果與實驗進行對比，如圖9所示。數值模擬所呈現的試樣破壞與實驗具有較高的吻合程度。數值模擬上也呈現與實驗類似的結果：圍壓一定時，隨著加載率的增大，試樣破壞產生的裂紋越多；在加載率相差不大的情況下，圍壓增大裂紋的數量減少。文獻[13]沖剪實驗過程中由于圍壓缸的存在，無法對沖剪實驗破壞的全過程進行觀測。數值模擬可以直觀的展示沖剪實驗破壞的全部過程，對研究沖剪破壞的機制具有一定參考價值。采用ABAQUS軟件對砂巖動態沖剪實驗模擬的破壞過程如圖10所示，通過試樣的剖切圖展示了試樣裂紋的產生以及生長過程。將實驗無法觀測的裂紋產生以及生長過程，通過模擬呈現出來，具有較高的參考意義。

表2 a～h工況匯總

3 預測性研究

文獻[13]未對砂巖施加更高的圍壓，但現如今國內礦物的開采深度已經達到2 000 m，國外礦產開采深度甚至更深[32]。在硐室開挖施工中，工程爆破所產生的動態剪切破壞是不可忽視的?；诖?，本文對30、40 MPa圍壓下巖石的動態剪切強度特性進行探究，為工程施工提供一定的參考價值。模擬沖剪破壞時，對模型輸入不同振幅與脈寬的正弦波，實現巖樣的沖剪實驗的數值模擬。動態剪切強度的變化規律如圖11所示。

將模擬所得的數據點進行線性擬合，數值大小關系上滿足圖11所擬合的直線方程，從圖中可以看出兩組圍壓下的數據都均勻地分布在直線兩側。在加載率相差不大的情況下，隨著圍壓的增大巖樣的動態剪切強度增大；在圍壓一定的情況下，巖樣的動態剪切強度隨著加載率的增加而增大。以30 MPa圍壓，2 350 GPa/s加載率和40 MPa圍壓，2 450 GPa/s加載率為例對比可以看出：在加載率相差不大的條件下，巖樣的破壞程度隨著圍壓的增大而減少。這一結果也與低圍壓狀態下的實驗與數值結果具有一定的相似性。

4 結論

1)通過模擬結果與實驗結果對比可知，模擬所得的加載率、剪切強度和試樣的破壞模式與實驗結果基本一致。表明基于ABAQUS/explicit平臺，采用剪切損傷模型和D-P塑性模型作為材料的本構，能夠重現沖剪實驗現象，反映巖石的動態剪切強度與破壞特征。

2)文獻[13]進行沖剪實驗時，由于圍壓缸的存在不能實現對巖石破壞過程的觀測。而數值模擬能夠實現巖石的沖剪破壞過程的可視化，可以觀測巖石的破裂過程，對研究巖石的破裂特性有一定價值。

3)基于模型對0～20 MPa圍壓下模擬與實驗對比呈現良好的吻合性，本文選擇不同振幅與脈寬的正弦波形對30、40 MPa圍壓下試樣的沖剪破壞進行了模擬預測。結果表明：隨著圍壓的增大，巖樣的破壞程度減弱；加載率相差不大的條件下，圍壓的增大巖樣的動態剪切強度增大；圍壓一定的條件下，巖樣的動態剪切強度隨著加載率的增大而增大。該結果也與引文[13]低圍壓狀態下的破壞規律具有一定的相似性，在高地應力施工時，可以對動態剪切強度進行預測估計，具有一定的參考價值。