基于有限元的巖石爆炸波傳播分析

周升 鄧康宇

(1、湖南有色新田嶺鎢業有限公司,湖南 長沙 423000 2、中國礦業大學深部巖土力學與地下工程國家重點實驗室,江蘇 徐州 221116)

1 概述

在礦山開采中,爆破技術常常被用于破碎巖石。在巖石爆破過程中,鉆孔處會承受載荷[1-2]。了解爆破時波在隧道結構中的傳播過程是保證爆破操作安全可靠的關鍵。本文采用COMSOL 有限元軟件,對巖石表面上的短時載荷引起的巖體中的波傳播進行了瞬態模擬分析。

2 數值模擬

2.1 模型建立

圖1 有限元幾何模型

本文使用了長方體模型幾何。模型的四個側壁中有兩個是對稱平面。另外兩個側面的作用是截斷計算域,在該切面上,巖石的深度尺寸遠大于切向尺寸。長方體的尺寸及材料參數采用了參考文獻1 提供的數據,爆破發生在兩個平面與底面相交的位置。幾何模型如圖1 所示,其中寬深高都為0.24m。采用映射和掃掠方式來劃分網格,網格總共包含13824 個單元,最小質量為1.0,平均質量為1.0。

巖石材料為花崗巖,固體模型為各向同性,默認為線彈性材料,坐標系為全局坐標系,其中具體材料參數如表1 所示。

表1 材料屬性

爆炸采用的固體力學方程為：

其中,u 為位移,Fv為體積力,∈為應變,C 為材料的各向屬性矩陣。

結構瞬態特性包含慣性項,力矩計算參考點x、y 和z 方向都設置為0,左前方和右前方兩個邊界為對稱,位移場采用線性單元來離散化,因變量位移場分量分別設置為u、v 和w。左后方和右后方兩個邊界為低反射邊界條件,上邊界為自由邊界,在底部邊界施加壓力脈沖載荷大小如下：

圖2 底部施加載荷函數曲線

載荷表示巖石內部靠近表面的地方發生了爆炸,載荷函數曲線如圖2 所示[3-4]。對截斷邊界分別采用低反射邊界條件和無反射邊界條件進行建模,阻尼類型為P 波和S 波。默認情況下,低反射邊界條件會從相鄰域中獲取材料數據,嘗試為壓力波和剪切波創建完美的阻抗匹配,因此：

其中n 和t 分別邊界處的單位法向和切向矢量,cp和cs是材料中壓力波和剪切波的速度。當波方向接近壁的法向時這種方法效果最佳。

研究類型采用低反射邊界條件和不采用低反射邊界條件,研究總時長為150ms,時間步長為2ms,采用默認的PARDISO直接求解器和全耦合求解器來對模型進行求解計算。PARDISO直接求解器預排序算法為自動,調度方法為自動,分別勾選“行預牌序”、“重用預排序”、“Bunch-Kaufman 主元”、“多線程前推和后溯求解”。主元擾動設置為1×10-9,勾選“用于集群的并行直接稀疏求解器”,核外模式選擇“自動”,核外的內存分數設置為0.99,核內內存法選擇“自動”,最小核內存為512MB,總內存使用比例為0.8,內部內存使用因子為3。檢查誤差估計設置為自動,誤差估計因子設置為1,勾選迭代求精,最大網格細化數位15,誤差率范圍為0.5。

全耦合求解器的非線性方法為恒定牛頓,阻尼系數設置為1,雅克比矩陣更新設置為最小,終止技術為容差判據,最大迭代次數設置為4,容差因子設置為1,終止準則為解。

2.2 結果分析

圖3 彈性波開始傳播時長方體中的應力

針對長方體中波傳播的瞬態研究,我們設置了150ms 的時間間隔。圖3 展示了波的典型傳播模式。該圖是巖石四分體爆炸波的題應力張量和Z 分量位移等值面圖。從圖3 中我們可以發現,隨著時間的推移,爆炸波傳播范圍越來越廣,波形形狀類似于鐘罩型,并且,爆炸波傳播的越遠,能量越弱,體應力張量越小,最大爆炸位移,在底邊界的有效爆炸半徑內。

圖4 在使用和不使用低反射邊界條件的情況下,巖石上表面的垂直位移

3 結論

3.1 巖石爆破一般取四分之一的對稱模型,然后可以利用鏡像還原出原始巖石爆破結果,這樣做可以減小模擬計算量,更快得到需要的結果。

3.2 模擬大型振動結構中波的傳播是一項極具挑戰性的任務,工程師必須在減小計算域尺寸與減少表面邊界的反射之間做出平衡。利用COMSOL 軟件的低反射邊界條件,我們可以輕松將計算域減小到合理大小,同時保證仿真結果的精確度[5-6]。

3.3 低反射邊界條件默認采用了相鄰域的材料數據,從而創建對壓力波和剪切波的完美阻抗匹配。這種方法最適用于波方向接近于壁法向的情況。本文介紹的內容只是使用低反射邊界條件簡化建模過程的一個案例,后續讀者可以在本文模型的基礎上對其他仿真研究進行優化[7]。

3.4 求解器的選擇,在有限元數值模擬中非常重要,其中線性問題和非線性問題的求解器又差別很大,對于多物理場的耦合問題,建議采用全耦合處理器。直接求解器一般收斂性比較高,但是計算時間長,耗費內存。非線性問題,一般采用迭代求解器,其對初值選擇,要求很高。