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

周升 鄧康宇

(1、湖南有色新田嶺鎢業(yè)有限公司,湖南 長(zhǎng)沙 423000 2、中國(guó)礦業(yè)大學(xué)深部巖土力學(xué)與地下工程國(guó)家重點(diǎn)實(shí)驗(yàn)室,江蘇 徐州 221116)

1 概述

在礦山開(kāi)采中,爆破技術(shù)常常被用于破碎巖石。在巖石爆破過(guò)程中,鉆孔處會(huì)承受載荷[1-2]。了解爆破時(shí)波在隧道結(jié)構(gòu)中的傳播過(guò)程是保證爆破操作安全可靠的關(guān)鍵。本文采用COMSOL 有限元軟件,對(duì)巖石表面上的短時(shí)載荷引起的巖體中的波傳播進(jìn)行了瞬態(tài)模擬分析。

2 數(shù)值模擬

2.1 模型建立

圖1 有限元幾何模型

本文使用了長(zhǎng)方體模型幾何。模型的四個(gè)側(cè)壁中有兩個(gè)是對(duì)稱平面。另外兩個(gè)側(cè)面的作用是截?cái)嘤?jì)算域,在該切面上,巖石的深度尺寸遠(yuǎn)大于切向尺寸。長(zhǎng)方體的尺寸及材料參數(shù)采用了參考文獻(xiàn)1 提供的數(shù)據(jù),爆破發(fā)生在兩個(gè)平面與底面相交的位置。幾何模型如圖1 所示,其中寬深高都為0.24m。采用映射和掃掠方式來(lái)劃分網(wǎng)格,網(wǎng)格總共包含13824 個(gè)單元,最小質(zhì)量為1.0,平均質(zhì)量為1.0。

巖石材料為花崗巖,固體模型為各向同性,默認(rèn)為線彈性材料,坐標(biāo)系為全局坐標(biāo)系,其中具體材料參數(shù)如表1 所示。

表1 材料屬性

爆炸采用的固體力學(xué)方程為：

其中,u 為位移,Fv為體積力,∈為應(yīng)變,C 為材料的各向?qū)傩跃仃嚒?/p>

結(jié)構(gòu)瞬態(tài)特性包含慣性項(xiàng),力矩計(jì)算參考點(diǎn)x、y 和z 方向都設(shè)置為0,左前方和右前方兩個(gè)邊界為對(duì)稱,位移場(chǎng)采用線性單元來(lái)離散化,因變量位移場(chǎng)分量分別設(shè)置為u、v 和w。左后方和右后方兩個(gè)邊界為低反射邊界條件,上邊界為自由邊界,在底部邊界施加壓力脈沖載荷大小如下：

圖2 底部施加載荷函數(shù)曲線

載荷表示巖石內(nèi)部靠近表面的地方發(fā)生了爆炸,載荷函數(shù)曲線如圖2 所示[3-4]。對(duì)截?cái)噙吔绶謩e采用低反射邊界條件和無(wú)反射邊界條件進(jìn)行建模,阻尼類型為P 波和S 波。默認(rèn)情況下,低反射邊界條件會(huì)從相鄰域中獲取材料數(shù)據(jù),嘗試為壓力波和剪切波創(chuàng)建完美的阻抗匹配,因此：

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

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

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

2.2 結(jié)果分析

圖3 彈性波開(kāi)始傳播時(shí)長(zhǎng)方體中的應(yīng)力

針對(duì)長(zhǎng)方體中波傳播的瞬態(tài)研究,我們?cè)O(shè)置了150ms 的時(shí)間間隔。圖3 展示了波的典型傳播模式。該圖是巖石四分體爆炸波的題應(yīng)力張量和Z 分量位移等值面圖。從圖3 中我們可以發(fā)現(xiàn),隨著時(shí)間的推移,爆炸波傳播范圍越來(lái)越廣,波形形狀類似于鐘罩型,并且,爆炸波傳播的越遠(yuǎn),能量越弱,體應(yīng)力張量越小,最大爆炸位移,在底邊界的有效爆炸半徑內(nèi)。

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

3 結(jié)論

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

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

3.3 低反射邊界條件默認(rèn)采用了相鄰域的材料數(shù)據(jù),從而創(chuàng)建對(duì)壓力波和剪切波的完美阻抗匹配。這種方法最適用于波方向接近于壁法向的情況。本文介紹的內(nèi)容只是使用低反射邊界條件簡(jiǎn)化建模過(guò)程的一個(gè)案例,后續(xù)讀者可以在本文模型的基礎(chǔ)上對(duì)其他仿真研究進(jìn)行優(yōu)化[7]。

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