張 飛，張東賓，王宇偉

(中國礦業(yè)大學(北京)力學與建筑工程學院，北京 100083)

0 引 言

在地礦開采過程中，往往對巖體施加爆炸、沖擊等強動荷載作用使其破碎，而巖體中富集的結(jié)構(gòu)面影響其動態(tài)力學行為，給工程活動的順利進行帶來了諸多困擾。巖體中的結(jié)構(gòu)面形式多樣，如孔洞、節(jié)理、裂隙等，其幾何特征對巖體的力學特性影響也各有不同。眾多學者對含結(jié)構(gòu)面巖體的動態(tài)力學特性開展了廣泛的研究工作。

譚卓英等[1]對含圓形孔洞結(jié)構(gòu)面巖體的沖擊動態(tài)響應(yīng)進行了實驗?zāi)M研究，結(jié)果表明孔洞結(jié)構(gòu)對巖體沖擊破壞形式、峰值應(yīng)力水平等有一定影響。朱哲明等[2]對含缺陷巖體的爆炸動態(tài)響應(yīng)進行了研究，發(fā)現(xiàn)不同缺陷種類的巖體動態(tài)強度和破壞范圍有所差別。楊仁樹等[3]通過實驗分析了含貫通裂紋的類巖石介質(zhì)動態(tài)破裂行為和預(yù)應(yīng)力場中含層理類巖石介質(zhì)爆破破裂的動態(tài)響應(yīng)。王蒙等[4]基于SCSCC試件對巖石復合型裂紋的動態(tài)擴展特性進行了實驗研究。YAVUZ等[5]對含多裂紋無限大平板進行分析，得出裂紋間距和長度對應(yīng)力強度因子的影響規(guī)律，并給出應(yīng)力強度因子計算式。GREGOIRE等[6]在進行沖擊試驗時發(fā)現(xiàn)動態(tài)裂紋擴展時有止裂又繼續(xù)擴展的現(xiàn)象。ZHU等[7]對巷道圍巖的動態(tài)響應(yīng)問題進行了數(shù)值模擬和實驗研究，并分析了應(yīng)力波作用下的動態(tài)響應(yīng)及加載方式對動態(tài)響應(yīng)的影響。AVACHAT等[8]對巖體的動態(tài)斷裂行為引入高速攝影法進行了影像學分析。FUNATSU等[9]在分析巖體斷裂破壞時，考慮了壓力和溫度等因素對砂巖動態(tài)破壞的影響。魏晨慧等[10]和謝冰等[11]對節(jié)理角度和間距變化時的巖體爆炸動態(tài)響應(yīng)進行了數(shù)值模擬研究。以上諸多研究從多個角度對巖體的動態(tài)力學行為進行了研究，分析了動態(tài)應(yīng)力強度因子、裂紋擴展行為等動力學特性。本文基于ABAQUS數(shù)值分析平臺，建立模型，分析巖體結(jié)構(gòu)面幾何特征變化時的沖擊動態(tài)響應(yīng)，以期找到結(jié)構(gòu)面幾何特征對巖體沖擊動態(tài)響應(yīng)的影響。

1 應(yīng)力波衰減規(guī)律

應(yīng)力波在巖體等非連續(xù)性介質(zhì)中傳播時，往往受地質(zhì)條件、結(jié)構(gòu)面幾何特征等因素影響，其反射透射及衰減狀況尤為復雜，眾多學者提出了一些理論模型，如接觸界面模型、完好黏結(jié)界面模型、位移不連續(xù)體模型等?；谶@些理論模型，應(yīng)力波經(jīng)過非連續(xù)界面時的傳播規(guī)律了較為完善的理論支撐，如李業(yè)學等[12]推導了應(yīng)力波穿越節(jié)理時反射透射系數(shù)FR1、FT1的解析解，其受節(jié)理等效剛度、波阻抗、頻率等因素影響，見式(1)和式(2)。

(1)

(2)

應(yīng)力波在巖體中傳播時，受到阻滯作用發(fā)生衰減，攜帶能量耗散；巖體某一質(zhì)點位移表示為式(3)。

u(x,t)=Aei(ωt-Ux)

(3)

其中：

(4)

(5)

(6)

式中：A為質(zhì)點振幅；λ為波長；G為剪切模量；η為黏滯系數(shù)；ρ為密度；ω為應(yīng)力波頻率。

聯(lián)立式(4)～式(6)可得式(7)。

U=β-iα

(7)

其中：

α=

β=

將式(7)代入式(3)可得式(8)。

u(x,t)=Ae-αxei(ωt-βx)

(8)

由此可知，應(yīng)力波在巖體中傳播時以式(8)的函數(shù)關(guān)系衰減，并且與巖體的黏滯系數(shù)、剪切模量、密度、應(yīng)力波頻率等因素有關(guān)。

2 模型建立

巖體本構(gòu)模型選擇德魯克-普拉格準則(以下簡稱“D-P準則”)作為失效準則。D-P準則同時量化了體積應(yīng)力、剪應(yīng)力和中間主應(yīng)力對巖石彈塑性強度的影響，能理想地反映巖石的力學行為。

為分析巖體結(jié)構(gòu)面幾何特征不同時沖擊荷載作用下的動態(tài)響應(yīng)，建立A組、B組、C組三組模型。平面應(yīng)力狀態(tài)如圖1所示，外觀尺寸150 mm×50 mm，上部中點受沖擊荷載，速度30 m/s，下部及左右部設(shè)置為約束邊界，x、y、z方向位移為0；網(wǎng)格劃分時，單元形狀為四邊形，采用進階算法(advancing front)生成均勻網(wǎng)格，單元屬性為四結(jié)點雙線性平面應(yīng)力四邊形單元；為便于分析，取分析點及分析路徑，以A組模型為例，模型A1結(jié)構(gòu)面左側(cè)端部分析點記為A1L，右側(cè)端部分析點記為A1R，下部中間分析點記為A1D，左右兩側(cè)以端部為起始點豎直向下方向上分別取分析路徑記為path1和path2，其他模型同理。

圖1 模型示意圖Fig.1 Diagram of model

選用ABAQUS/Explicit，其本質(zhì)上是基于動力學方程求解，對于爆炸、沖擊的數(shù)值模擬具有很好的仿真效果。

3 結(jié)果分析

3.1 結(jié)構(gòu)面厚度變化時的動態(tài)響應(yīng)

A組模型的結(jié)構(gòu)面厚度值分別為0.1 mm、0.2 mm、0.3 mm、0.4 mm，截取模型中間部分如圖2所示。

圖2 A組模型示意圖Fig.2 Diagram of group A model

四個分析點A1D、A2D、A3D、A4D的Mises應(yīng)力值隨時間變化如圖3所示。以沖擊塊體接觸模型為零時刻，在70 μs左右時，Mises應(yīng)力值出現(xiàn)峰值，分別為12.52 MPa、9.89 MPa、9.80 MPa、9.34 MPa。 模型A1結(jié)構(gòu)面厚度為0.1 mm，分析點A1D的Mises應(yīng)力峰值高于同組其他3個模型，但其他3個模型此時峰值差別不大，這可能是結(jié)構(gòu)存在某一特定厚度使應(yīng)力波沒有穿透結(jié)構(gòu)面所導致。隨著時間的延后，分析點A1D的Mises應(yīng)力值震蕩下降，在200 μs之前，其值基本上大于其他3個模型分析點，200 μs左右時，其他3分析點應(yīng)力峰值高于分析點A1D，這是因為應(yīng)力波發(fā)生繞射，導致結(jié)構(gòu)面下部分析點處的應(yīng)力值又有所上升，200 μs之后，Mises應(yīng)力值震蕩下降。

圖3 下部分析點應(yīng)力值Fig.3 Stress of down monitoring point

以模型A1為例，在分析路徑path1上，其Mises應(yīng)力值隨距離衰減曲線如圖4所示。擬合函數(shù)為f(x)=ax-b，三個時刻的擬合函數(shù)關(guān)系式分別為：f1(x)=28.13x-0.43，f2(x)=23.72x-0.43，f3(x)=18.06x-0.40，符合應(yīng)力波衰減規(guī)律。分析路徑path2情況與此類似。

圖4 應(yīng)力值衰減曲線Fig.4 Curve fitting of stress attenuation

模型A1應(yīng)力云圖變化如圖5所示。60 μs時結(jié)構(gòu)面兩側(cè)以端部為中心，并主要沿下側(cè)、外側(cè)逐漸擴散，可以觀察到結(jié)構(gòu)面下部區(qū)域云圖顏色偏冷色，這說明應(yīng)力波繞射、透射到此區(qū)域較少。120 μs左右，應(yīng)力云圖持續(xù)保持以結(jié)構(gòu)面兩端為中心，主要向下側(cè)、外側(cè)擴展，向結(jié)構(gòu)面下部中間集中程度并不大，這正是所取分析點Mises應(yīng)力值震蕩變化的直觀現(xiàn)象，其根本原因是應(yīng)力波的持續(xù)傳播，之所以出現(xiàn)不只一次的峰值，是因為結(jié)構(gòu)面尖端集聚能量釋放的結(jié)果，180 μs左右也基本呈現(xiàn)此現(xiàn)象，云圖顏色逐漸轉(zhuǎn)冷色調(diào)。 實際巖體承受強動荷載時，弱面處出現(xiàn)應(yīng)力集中，集聚較大應(yīng)變能，并常從此起裂破壞。

圖5 應(yīng)力云圖變化圖Fig.5 Changing diagram of stress cloud

3.2 結(jié)構(gòu)面長度變化時的動態(tài)響應(yīng)

B組四個模型結(jié)構(gòu)面長度分別為10 mm、20 mm、30 mm、40 mm，截取模型中間部分如圖6所示。

圖6 B組模型示意圖Fig.6 Diagram of group B model

四個分析點B1D、B2D、B3D、B4D的Mises應(yīng)力值隨時間變化見圖7。分析點B1D和分析點B2D變化規(guī)律相似，先震蕩上升，200 μs左右達到峰值，分別為18.02 MPa和16.15 MPa，后者較小，隨后震蕩下降，這是因為隨著結(jié)構(gòu)面長度的增加，應(yīng)力波較難以繞射到結(jié)構(gòu)面下部區(qū)域。同樣，分析點B3D的峰值和整體應(yīng)力值水平更低。但隨著結(jié)構(gòu)面長度的進一步增加，情況有所不同，模型B4的結(jié)構(gòu)面長度為40 mm，分析點B4D在100 μs左右時達到峰值19.23 MPa，隨后有很短暫的下降后又快速上升，在160 μs左右時再次達到峰值18.86 MPa，隨后震蕩下降。

圖7 下部分析點應(yīng)力值Fig.7 Stress of down monitoring point

以模型B4為例，在分析路徑path1上，其Mises應(yīng)力值隨距離衰減曲線如圖8所示。擬合函數(shù)為f(x)=ax-b，三個時刻的擬合函數(shù)關(guān)系式分別為：f1(x)=50.44x-0.59、f2(x)=42.47x-0.60、f3(x)=46.88x-0.60，符合應(yīng)力波衰減規(guī)律。分析路徑path2情況與此類似。

圖8 應(yīng)力值衰減曲線Fig.8 Curve fitting of stress attenuation

Mises應(yīng)力云圖變化如圖9所示，90 μs左右，結(jié)構(gòu)面下部開始出現(xiàn)“水滴”狀顏色轉(zhuǎn)暖色區(qū)域，隨后快速消失，在150 μs左右時再次出現(xiàn)此區(qū)域；同時可以觀察到結(jié)構(gòu)面?zhèn)让婧拖聜?cè)面接觸，這說明結(jié)構(gòu)面的應(yīng)變積累，結(jié)構(gòu)面的豎向位移值超過了結(jié)構(gòu)面厚度，相當于節(jié)理的等效剛度降低，透射系數(shù)增加，在模型建立階段結(jié)構(gòu)面的上下側(cè)表面相互作用設(shè)置為硬接觸，上表面將應(yīng)力傳遞到下表面，應(yīng)力波發(fā)生了透射現(xiàn)象。

圖9 應(yīng)力云圖變化圖Fig.9 Changing diagram of stress cloud

3.3 結(jié)構(gòu)面曲率變化時的動態(tài)響應(yīng)

C組模型結(jié)構(gòu)面變化情況如圖10所示。為直觀表示，圖中數(shù)值為結(jié)構(gòu)面的曲率半徑大??；由曲率半徑r和曲率k之間的關(guān)系式k=1/r可知兩者呈反比，即曲率越來越小，依次為500、333、250、200。

圖10 C組模型示意圖Fig.10 Diagram of group C model

取結(jié)構(gòu)面左側(cè)端部四個分析點C1L、C2L、C3L、C4L，其Mises應(yīng)力值隨時間變化情況如圖11所示。當結(jié)構(gòu)面曲率越來越小時，Mises應(yīng)力值在60 μs之前的增加速度越來越大，結(jié)構(gòu)面曲率減小，結(jié)構(gòu)面趨于平直，應(yīng)力波傳播到端部的距離縮短，因此端部的應(yīng)力值會在較短時間內(nèi)開始快速增加。分析點C1L在120 μs左右時第一次達到較大值為35.78 MPa，隨后震蕩下降；分析點C2L在60 μs左右時首次達到峰值42.57 MPa，隨后震蕩下降，但隨著結(jié)構(gòu)面曲率的進一步減小；分析點C3L在50 μs左右時便首次出現(xiàn)峰值，大小為42.61 MPa；分析點C4L的峰值出現(xiàn)時間更早，基本上在50 μs左右，這和結(jié)構(gòu)面幾何特征是平直時的端部分析點出Mises應(yīng)力值峰值出現(xiàn)時間基本一致，這說明隨著曲率減小，結(jié)構(gòu)面趨于平直，其端部Mises應(yīng)力值在應(yīng)力波傳播初期有著相似的變化情況。分析點C3L在首次峰值出現(xiàn)后有短暫時間的緩慢下降，趨勢較為平直，這可能是時間間隔不夠短引起的誤差，其下降趨勢應(yīng)該是震蕩下降，分析點C4L出現(xiàn)不只一次峰值后下降。

圖11 左側(cè)分析點應(yīng)力值Fig.11 Stress of left monitoring point

應(yīng)力波在傳播過程中攜帶的能量主要轉(zhuǎn)化為應(yīng)變能儲存在介質(zhì)內(nèi)部，取分析點處的應(yīng)變能隨時間變化情況如圖12所示。結(jié)構(gòu)面左側(cè)端部的應(yīng)變能隨時間變化趨勢基本和Mises應(yīng)力值隨時間變化趨勢保持協(xié)調(diào)。分析C1L在100 μs左右時應(yīng)變能積累到最大值，隨后趨于下降；分析點C2L的應(yīng)變能在60 μs左右時達到峰值，隨后震蕩下降；分析點C3L和C4L的應(yīng)變能在50 μs左右快速上升到峰值，隨后趨于較平緩下降，在200 μs左右時快速降低；可見隨著結(jié)構(gòu)面曲率的減小，結(jié)構(gòu)面端部積累應(yīng)變能的情況有所不同，曲率大時，也即結(jié)構(gòu)面越來越“圓”，此時應(yīng)力波從結(jié)構(gòu)面弧所對應(yīng)的弦的中垂線方向入射，結(jié)構(gòu)面反射較多的應(yīng)力波，應(yīng)力波攜帶的能量也就不容易在此聚集，隨著結(jié)構(gòu)面趨于平直，應(yīng)力波反射情況趨同于直線型結(jié)構(gòu)面。

圖12 左側(cè)分析點應(yīng)變能Fig.12 Strain energy of left monitoring point

4 結(jié) 論

1) 在以結(jié)構(gòu)面端部為起始點豎直向下的分析路徑上，Mises應(yīng)力值以函數(shù)關(guān)系式f(x)=ax-b衰減。

2) 結(jié)構(gòu)面厚度增大，應(yīng)力波難以發(fā)生透射，下部分析點應(yīng)力值水平降低；結(jié)構(gòu)面長度增加，上下表面容易發(fā)生閉合現(xiàn)象，應(yīng)力波發(fā)生透射，下部分析點應(yīng)力值會增大。

3) 結(jié)構(gòu)面曲率減小，即結(jié)構(gòu)面趨于“平直”時，結(jié)構(gòu)面端部應(yīng)力峰值出現(xiàn)時間早，說明更容易積累應(yīng)變能；應(yīng)變能與應(yīng)力值變化情況保持協(xié)調(diào)。