應力波在各向異性介質中傳播及損傷數值模擬★

鞠冬冬，王 清，陳 鵬，胡少斌，蔡余康

(1.南京市江北新區公共工程建設中心,江蘇 南京 210098; 2.河海大學土木與交通學院,江蘇 南京 210098)

在人類的生活中經常有諸如沖擊、爆炸、碰撞等強荷載作用下材料發生高應變的現象。比如地震、洪水、交通事故這樣的災難,還有礦山開采、鉆爆法開挖、結構抗震等造福民生的技術方法。這些現象中諸如地震、礦山開采以及鉆爆法開挖都可納入巖石力學范疇,更具體地講,屬于巖石動力學范疇,即應變率在1×10-1～1×104s-1之間的情況[1]。而巖石動力學以應力波作為其理論基礎,主要研究涵蓋應力波在巖體中的傳播與衰減[2-4],應力波與節理裂隙的相互作用[5-7],應力波與結構面的透、折、反射關系[8-10],應力波入射角[11-12]等多個方面。除此之外,巖石動力學所研究的另一個基本問題是巖石或巖體受沖擊后的動態響應,包括變形、損傷、斷裂的結果以及演化過程。目前,對應力波傳播及巖石的動態響應的研究多種多樣,從一開始常規的應力波的傳播和衰減規律[13],到后來探究循環荷載下的巖石力學響應[14-15],進一步地研究其裂紋擴展和損傷演化的規律[16-18],并考慮不同的工況如高地應力爆破和水下深孔爆破等。模擬研究的對象也從各向同性巖石擴展到正交各向異性[19]和各向異性巖石。具體到工程應用領域,如典型的巖石動力學應用范疇——礦山開采及鉆爆法破巖。這兩個領域經典的核心施工流程是鉆孔、裝藥以及爆破。然而從目前來看,還很少有人基于應力波對該過程的巖體,尤其是各向異性巖體損傷演化以及能量耗散進行研究,從而確定更合理的布孔間距和布孔方式。本文將采用數值模擬的手段對無限元邊界下各向異性巖體在孔周脈沖波作用情況下巖體的損傷演化展開研究,并根據損傷程度來確定各向異性條件下合理的布孔間距和布孔角。

1 數值建模

本章將基于COMSOL MULTIPHYSICS建立圍壓狀態下爆破花崗巖層的數值模擬模型,并給出本模型的相關控制方程、幾何建模、網格劃分、邊界條件,以及模型所涉及的材料參數。

1.1 控制方程

本模型采用相場損傷法來描述線彈性巖石的損傷應變軟化過程,其本構方程采用胡克定律:

σ=D×ε

(1)

其中,σ為彈性應力;D為考慮損傷的彈性矩陣;ε為彈性應變。

本文研究應力波在各向異性介質中的傳播,故彈性矩陣D自定義為:

(2)

由于考慮了應變軟化過程,引入關于損傷的裂紋相場因子,其相場損傷演化函數fd為:

fd=1-(1-φ)2

(3)

其中,φ為裂紋相場因子,與fd間存在弱函數關系約束,弱表達式為:

(4)

其中,Dd為裂紋驅動力;test()為試函數;l為長度尺寸;φx和φy分別為φ在x和y方向上的導數。同時,為確定固體的應力應變,還需要加入平衡微分方程,考慮慣性項后為:

(5)

其中,ρs為固體密度;x為固體位移;f為單位體積力;t為時間。此外,還需通過變形協調方程來根據應變確定固體位移,表達式為:

(6)

1.2 幾何模型

本模型模擬爆破孔處由爆破或高能氣體壓裂產生的沖擊應力波的巖層中的傳播,故采用長寬為8 m×4 m的模型尺寸,取模型尺寸寬度的1/8為無限元域層以表示有限域內的無限空間,其物理寬度為8 944 m,爆破孔直徑為160 mm,間距為Δd,爆破孔圓心相連的軸線與豎向軸的夾角為α,所加荷載為指數函數與正弦函數的乘積形式,其表達式為:

(7)

其中,a為放大系數,取10;P0為作用力大小取7 000 kN;γ為衰減率,取1.86;R為爆破孔直徑;l為爆破孔深度;t0為加載時間,取2.7×10-8s。

具體模型如圖1所示。

在二維平面的網格劃分當中,主要有三角形網格單元和四邊形網格單元兩種,相比于三角形網格單元,四邊形網格具有更為平順,更適用于高度非線性問題,同等密度下網格數更少,收斂性更強等優點[20-22],因此本模型中,周邊無限元域采用四邊形網格,中間采用三角形網格,在接近圓孔處進行加密,如圖2所示,網格數共11 336,一定程度上保證計算結果的精度。

本模型邊界的水平載荷與豎向荷載皆取ρ=5 MPa,所采用的花崗巖參數如表1所示。

表1 花崗巖參數表

2 模型結果與分析

本章展現了無限元邊界下應力波在花崗巖層中的傳播過程以及伴隨發生的損傷演化過程。通過對爆破孔間距和夾角α進行參數化掃描得到了不同工況下的模擬結果。同時令損傷大于0.9的面積作為有效破壞面積,將有效破壞面積作為爆破效果的評價指標,展示各向異性條件下孔間距和軸線轉角對應力波沖擊破壞巖石的影響。

2.1 應力波在各向異性巖層中的傳播及損傷規律

本文研究爆破孔間距和軸線轉角對應力波傳播以及損傷演化的影響,因此對間距和轉角進行參數化掃描,即在完全相同的模型條件及求解器設置下計算僅改變指定參數下的所有結果,具體模擬組設置如表2所示,其中峰值應力考慮到爆破時的放大效應,設為435 MPa。

表2 模擬組參數

本文將A組作為對照組,其應力波傳播過程如圖3所示,可以看出應力波的橫向傳播速度顯著大于豎向傳播速度,由彈性波波速計算公式v=√(E/P)可得彈性波在水平方向的相速度為7 817.4 m/s,在豎直方向上的相速度為4 792.0 m/s,符合圖3中t=2×10-5s的結果。同時可以觀察到水平方向的應力遠大于豎直方向上的應力,聯系各向同性條件下的計算結果(如圖4所示)可知,這是由彈性模量在方向上的差異性所導致的。

A組的損傷演化過程如圖5所示,由于本模型采用的是理想的線彈性模型,沒有考慮到巖石的脆性,因此裂隙區較不明顯,對于破碎區,將損傷面積(d>0.9)作為等效的破碎面積,其等效直徑為爆破孔的7倍,符合其他學者得出的破碎區為孔徑2倍～10倍的結論[23],從圖5中還可以看出兩個爆破孔的連線中點處損傷較中點兩側更高,從時間在0.296 ms的損傷云圖中可以發現中點處更早發生破碎,這是由于應力波在中點處產生了疊加,這點在下文關于孔間距的影響中會進一步說明。

2.2 軸線轉角對應力波傳播及損傷的影響

在各向異性的均質巖石中,應力波向不同的方向傳播會產生有規律性的差異,從各組的應力結果中(如圖6所示)可以看出不同軸線轉角下的峰值應力相對近似,且由于距離水平邊界和豎直邊界的距離不等導致沒有明顯規律。

而觀察模型中心點的峰值應力曲線(如圖7所示)可以觀察到不同角度下第一個波峰的峰值和出現時間不同,其出現時間隨著轉角的減小逐個延后,而峰值應力也逐漸減小。但之后的波峰由于次波的干擾并未展現清晰的規律,其中D組相對其他組較大是由于傳播過程中產生的次波共振恰好處于中心位置。

如上所述本文將損傷d>0.9的區域作為破碎區,并得出A組—E組的損傷破碎面積,可以看出在不同轉角下,一開始的破碎面積擴展是一致的,這時雙孔處加載的應力波使各自臨近處發生破碎,但尚未傳播到鄰孔處,但在一定時間后,應力波傳播到鄰孔位置,對彼此的應力場造成擾動,同時由于各向異性條件下在不同方向上應力波的傳播速率和衰減率的差異導致破碎區的擴展產生了較為復雜的變化,但從8×10-4s的結果可以看出破碎區面積隨著角度的減小而增大,值得注意的是該時間點B組面積較C組更大(見圖8),然而從該點的切線斜率可以判斷在更大的時間尺度上C組的最終損傷面積更大。本文認為是角度更小的情況下,雙孔間彈性模量更小,其間得到損傷演化速率更快,兩孔間更容易貫通。

2.3 孔間距對應力波傳播及損傷的影響

爆破孔間距向來是爆破孔布置方案中的關鍵參數,對爆破效果有著重要影響。由圖9可以看出F,G兩組峰值應力較大,且出現多個波峰,而H,I兩組應力較小,且波動不明顯。表明孔間距較小時產生的峰值應力較大,且由于距離較近,相應傳播時間更短,因此應力波的疊加效應更加明顯。而從損傷演化曲線(如圖10所示)可以看出破碎區大致隨著距離的增加而增加,本文認為其原因為雙孔應力波的相對沖擊增強了能量的衰減,導致近孔間距爆破的能量利用率降低。

3 結語

本文通過數值模擬對應力波在各向異性巖層中傳播展開研究與分析,得到如下結論:

1)各向異性條件下,應力波在彈性模量更大的方向上傳播速率更快,且應力幅值較高。2)當軸線轉角不同時,雙孔中心點處應力的第一個波峰的出現時間隨著轉角的減小逐個延后,且峰值應力也逐漸減小。同時破碎區面積隨著角度的減小而增大。3)孔間距較小時產生的峰值應力較大,應力波的疊加效應更加明顯。破碎區面積大致隨著距離的增加而增加。