基于Fluent-EDEM流固耦合算法的水下臺(tái)階爆破爆堆形成過(guò)程研究*

吳 亮,余 創(chuàng),梁志堅(jiān),段衛(wèi)東,陳 明

(1.武漢科技大學(xué) 理學(xué)院 湖北省智能爆破工程技術(shù)研究中心,武漢 430065;2.武漢大學(xué) 水工巖石力學(xué)教育部重點(diǎn)實(shí)驗(yàn)室,武漢 430072)

眾所周知,巖石的破碎和運(yùn)動(dòng)是巖土爆破中的2個(gè)關(guān)鍵環(huán)節(jié),其中破碎塊體的運(yùn)動(dòng)常被稱拋擲運(yùn)動(dòng)。當(dāng)實(shí)際生產(chǎn)中礦巖與普通巖互為夾層時(shí),礦體運(yùn)動(dòng)結(jié)果對(duì)分離礦石和廢石的影響較大,有時(shí)對(duì)采礦效率的影響遠(yuǎn)大于機(jī)械設(shè)備[1]。水下鉆孔爆破作為工程爆破技術(shù)中的重要組成部分,廣泛應(yīng)用在擴(kuò)展港口、疏浚航道、取水以及搶險(xiǎn)救災(zāi)等工程中。水下與陸地巖石爆破的主要差異在于巖石表面接觸介質(zhì)的不同,直接影響了巖石的爆破破碎以及塊體拋擲運(yùn)動(dòng)效果[2,3]。由于影響水下爆破效果的因素錯(cuò)綜復(fù)雜,很多工程只能憑有關(guān)經(jīng)驗(yàn)公式與試驗(yàn)數(shù)據(jù)估算,往往沒(méi)有達(dá)到破碎與堆積效果;另外,室內(nèi)的模型實(shí)驗(yàn)成本高且周期長(zhǎng),而現(xiàn)場(chǎng)試驗(yàn)無(wú)法及時(shí)獲取水下爆破破碎與堆積效果,且很多工程也不具備現(xiàn)場(chǎng)試驗(yàn)條件。因此,水下爆破破碎與堆積效果的預(yù)測(cè)工作亟待解決。

由于爆破作用過(guò)程的復(fù)雜性,且現(xiàn)有數(shù)學(xué)模型還不完善或存在局限性,目前工程中對(duì)爆破效果的預(yù)測(cè)還主要依靠經(jīng)驗(yàn)判斷,其準(zhǔn)確性及可靠性難以滿足要求[4]。數(shù)十年來(lái),數(shù)值模擬技術(shù)在爆破領(lǐng)域得到了廣泛應(yīng)用,但傳統(tǒng)的連續(xù)介質(zhì)模型很難對(duì)巖石爆炸沖擊破碎后的塊體運(yùn)動(dòng)進(jìn)行模擬。YI采用顆粒模型(PBM)模擬了爆轟過(guò)程以及被爆的巖石等脆性介質(zhì)[5],實(shí)現(xiàn)巖石爆破破碎過(guò)程,突破了連續(xù)介質(zhì)模型難處理大變形的問(wèn)題。PREECE結(jié)合并行計(jì)算技術(shù)采用三維離散元模型分析了數(shù)百個(gè)炮孔[6],離散單元數(shù)超過(guò)100萬(wàn)。安華明等采用有限—離散雜交元模擬了爆破漏斗產(chǎn)生的整個(gè)動(dòng)態(tài)過(guò)程[7]。余德運(yùn)等采用流體與DDA塊體之間網(wǎng)格覆蓋的流固耦合技術(shù)模擬了爆破漏斗形成過(guò)程[8]。蘇都都等采用PFC 2D對(duì)爆破參數(shù)與爆堆形態(tài)的關(guān)系進(jìn)行了研究[4],預(yù)測(cè)了臺(tái)階爆破的爆堆形態(tài)。冷振東等在巖體強(qiáng)度參數(shù)服從Weibull分布規(guī)律的基礎(chǔ)上[9],采用3DEC離散元軟件建立了臺(tái)階爆破的三維臺(tái)階爆破模型,模擬起爆位置對(duì)巖石的動(dòng)態(tài)破碎與拋擲過(guò)程的影響規(guī)律。周旺瀟等和YAN等在考慮不同區(qū)域塊度尺寸分布的基礎(chǔ)上[10,11],采用3DEC離散元模型討論了的爆破漏斗與臺(tái)階爆破爆堆形態(tài)。可見(jiàn),預(yù)測(cè)與控制爆破拋擲運(yùn)動(dòng)的數(shù)值軟件的應(yīng)用極大地提高了生產(chǎn)效率,但水下巖石爆破破碎及其拋擲運(yùn)動(dòng)目前鮮見(jiàn)報(bào)道。

近年來(lái)研究者開(kāi)始采用CFD-DEM方法解決工程中的流固耦合問(wèn)題,例如水下巖塞爆破和深海采礦管道固液混合流體[12,13]。XU等基于Fluent-EDEM耦合方法模擬了管道堵塞過(guò)程[14],通過(guò)考慮連續(xù)流體相和離散顆粒相,初步研究了水合物的堵塞機(jī)理。HE等基于Fluent-EDEM引入滑移網(wǎng)格[15],該模型可以同時(shí)分析大型實(shí)體模型自由移動(dòng)引起的流場(chǎng)和小顆粒的運(yùn)動(dòng)。FU等采用CFD-DEM耦合的計(jì)算方法模擬了滲流力與掌子面巖體顆粒的相互作用[16]。WU等基于DEM中顆粒的接觸模型[17,18],采用Fluent-EDEM耦合方法研究水下臺(tái)階爆破塊體的運(yùn)動(dòng)過(guò)程及爆堆形態(tài)。隨著計(jì)算軟件的開(kāi)發(fā)與應(yīng)用,越來(lái)越多的研究者嘗試采用各種數(shù)值計(jì)算的方法解決水下爆破的問(wèn)題。因此,采用基于UDF開(kāi)發(fā)的Fluent-EDEM耦合方法,考慮顆粒間的連接作用來(lái)研究水下巖石臺(tái)階爆破破碎及塊體運(yùn)動(dòng)的問(wèn)題,期望研究成果能促進(jìn)水下臺(tái)階爆破技術(shù)發(fā)展。

1 計(jì)算方法

1.1 流體基本原理

計(jì)算流體力學(xué)是利用數(shù)學(xué)的計(jì)算方法來(lái)實(shí)現(xiàn)對(duì)流場(chǎng)中的控制方程進(jìn)行計(jì)算,并對(duì)有限元網(wǎng)格節(jié)點(diǎn)進(jìn)行數(shù)值求解。當(dāng)研究對(duì)象為不可壓縮流體時(shí),基本的定律有三個(gè):質(zhì)量守恒定律、動(dòng)量守恒定律、能量守恒定律,如下式(1)和式(2)。在本文中沒(méi)有能量的交換,故不考慮能量守恒方程。

(1)

(2)

1.2 DEM顆粒接觸模型

DEM模型中顆粒受力運(yùn)動(dòng)遵循牛頓第二定律,包括了轉(zhuǎn)動(dòng)和平動(dòng),顆粒的運(yùn)動(dòng)方程(動(dòng)量守恒方程)如下式(3)和(4)所示。

(3)

(4)

式中:v、ω、m分別為顆粒速度、角速度和質(zhì)量,下標(biāo)i,j表示兩個(gè)相互接觸的顆粒;Fij為i和j顆粒間的接觸力;Fg為顆粒間的引力;Ff為流體與顆粒間的相互作用力;I是顆粒轉(zhuǎn)動(dòng)慣量;T是顆粒間的力矩;Tf是流體引起的扭矩。

本文采用Hertz-Mindlin接觸模型計(jì)算DEM顆粒之間的作用力,該模型具備滾動(dòng)摩擦特性,如圖1所示。Hertz-Mindlin接觸模型的計(jì)算表達(dá)式見(jiàn)文獻(xiàn)[19]。

圖1 Hertz-Mindlin接觸模型示意圖Fig. 1 Schematic diagram of Hertz-Mindlin’s contact model

1.3 DEM顆粒平行鍵模型

平行鍵模型(Parallel bond Model,PBM)可以用于模擬微觀結(jié)構(gòu)如何影響宏觀行為,它將固體介質(zhì)離散為具有固定尺寸的膠結(jié)顆粒材料,由于顆粒之間存在空隙,在固體介質(zhì)模型形成平行鍵時(shí),這些顆粒間的空隙可以視為混凝土材料中的自然孔隙或巖石中的自然裂縫,這些裂縫相互作用、合并形成宏觀裂紋。平行鍵模型能表達(dá)顆粒間法向和切向的作用力,如圖2所示,當(dāng)法向或切向應(yīng)力超過(guò)最大許可應(yīng)力時(shí)平行鍵斷裂,顆粒之間的相互作用回歸到Hertz-Mindlin接觸模型。

圖2 平行鍵模型示意圖Fig. 2 Schematic diagram of parallel bond model

在平行鍵生成之后顆粒上的力和扭矩被重置為0,在模型開(kāi)始計(jì)算時(shí)則根據(jù)公式(5)～(8)在每個(gè)時(shí)間步長(zhǎng)中計(jì)算平行鍵受到的力和扭矩的增量。

δFn=-vnSnAδt

(5)

δFt=-vtStAδt

(6)

δMn=-ωnStJδt

(7)

(8)

式中:Fn和Ft為法向力和切向力;Mn和Mt則為法向扭矩和切向扭矩;vn和vt為顆粒的法向和切向速度;ωn和ωt為法向和切向角速度;Sn和St為法向剛度和切向剛度;δt為時(shí)間步長(zhǎng);A是平行鍵截面面積,計(jì)算公式如式(9)所示;J為平行鍵截面極慣性矩,計(jì)算公式見(jiàn)式(10)。

(9)

(10)

式中:RB為平行鍵半徑,是平行鍵合模型輸入?yún)?shù)之一,取值一般大于顆粒半徑。

當(dāng)法向力或切向力超過(guò)其限定值時(shí),顆粒間的平行鍵發(fā)生斷裂

(11)

(12)

式中:σmax、τmax為最大法向應(yīng)力和最大切向應(yīng)力。

2 水下臺(tái)階爆破計(jì)算模型

2.1 計(jì)算模型

參考實(shí)際臺(tái)階爆破工程,數(shù)值模型中上臺(tái)階面距離水面10 m、臺(tái)階面寬6 m,臺(tái)階高度為10 m,炮孔直徑0.105 m,孔深10 m,裝藥長(zhǎng)度為7 m,堵塞3 m,孔間距為4.5 m,最小抵抗線為3 m。本文數(shù)值模型以單段4個(gè)炮孔同時(shí)起爆為對(duì)象進(jìn)行建模,如圖3所示。圖3(a)為Fluent模型,上表面為壓力出口邊界,壓力值設(shè)置為一個(gè)大氣壓。圖3(b)為EDEM模型,顆粒填充區(qū)域?yàn)殚L(zhǎng)23.5 m、寬5.5 m、高10 m的長(zhǎng)方體,顆粒采用直徑0.21 m的球形顆粒,共填充146728顆。

圖3 計(jì)算模型(單位:m)Fig. 3 Calculation models in the numerical simulation(unit:m)

2.2 模型參數(shù)

在EDEM模型中引入平行鍵的模型來(lái)模擬連續(xù)的巖體介質(zhì),平行鍵的模型除了包含摩擦系數(shù)在內(nèi)的基本力學(xué)參數(shù)外還需要額外輸入5個(gè)平行鍵參數(shù)。由于平行鍵的模型能體現(xiàn)微觀上顆粒之間的力學(xué)行為,除了平行鍵半徑取顆粒半徑值外,至今的研究結(jié)果還不能明確地給出實(shí)驗(yàn)中獲得的參數(shù)與平行鍵參數(shù)之間的關(guān)系式,因此為了獲得準(zhǔn)確的數(shù)值計(jì)算參數(shù),需對(duì)平行鍵的參數(shù)進(jìn)行標(biāo)定。許多學(xué)者使用DEM進(jìn)行巴西劈裂數(shù)值計(jì)算[20-22],通過(guò)與實(shí)驗(yàn)結(jié)果對(duì)比,標(biāo)定平行鍵的參數(shù)。

在EDEM中建立直徑50 mm、高100 mm的圓柱模型,投放的球形顆粒直徑均為3 mm,共投放顆粒8473顆,共生成平行鍵29 258個(gè)。分別對(duì)圓柱顆粒模型進(jìn)行巴西劈裂數(shù)值模擬,模型破壞結(jié)果如圖4(a)所示,應(yīng)力-應(yīng)變曲線如圖4(b)。

圖4 巴西劈裂數(shù)值模擬Fig. 4 Numerical simulation of Brazilian splitting

圖4(a)中灰白色為顆粒間的平行鍵,紅色代表已破壞的平行鍵,巴西劈裂模型的破壞模式為拉伸破壞。圖4(b)中三組巴西劈裂實(shí)驗(yàn)曲線來(lái)自文獻(xiàn)[23]的實(shí)驗(yàn)。三組巴西劈裂實(shí)驗(yàn)的平均峰值應(yīng)力為2.46 MPa,而數(shù)值模型的峰值應(yīng)力為2.44 MPa。說(shuō)明基于平行鍵方法建立的材料模型能夠在一定程度上反映實(shí)際工況中巖石、混凝土等脆性材料的力學(xué)行為。

綜上,結(jié)合標(biāo)定后的平行鍵模型參數(shù),將堆積的顆粒建立平行鍵,臺(tái)階模型共生成平行鍵351 460個(gè),具體參數(shù)如下表1所示。需要注意的是,同一組模型中含有平行鍵與無(wú)平行鍵在時(shí)間步長(zhǎng)的選取上存在差異,通常在進(jìn)行含有平行鍵的數(shù)值模型計(jì)算時(shí)選取兩者中時(shí)間步長(zhǎng)較小的一個(gè)[24],本文離散元的計(jì)算步長(zhǎng)為3×10-5s,流體的計(jì)算步長(zhǎng)是離散元的10倍,為3×10-4s。

表1 巖體EDEM模型參數(shù)Table 1 EDEM model parameters of rock mass

3 模擬結(jié)果

3.1 臺(tái)階爆破破碎與塊體運(yùn)動(dòng)過(guò)程

數(shù)值模型中通過(guò)在炮孔位置設(shè)置高壓區(qū)域?qū)崿F(xiàn)爆炸荷載的沖擊作用,塊體顆粒運(yùn)動(dòng)速度與未破壞的平行鍵分布圖見(jiàn)圖5,該圖顯示了在有平行鍵作用下水下臺(tái)階塊體顆粒鼓包拋擲的整個(gè)運(yùn)動(dòng)過(guò)程。從臺(tái)階爆破破碎來(lái)看,初始狀態(tài)顆粒間的平行鍵完整密布,炮孔內(nèi)爆炸荷載作用在炮孔壁后,炮孔壁介質(zhì)顆粒向外傳遞爆炸荷載并向外運(yùn)動(dòng),當(dāng)顆粒間的平行鍵超過(guò)臨界值時(shí)將斷開(kāi),處于炮孔附近的平行鍵大量破裂,結(jié)果顯示平行鍵的斷裂是從炮孔壁逐步向外開(kāi)始擴(kuò)散的。由于臺(tái)階顆粒排列是隨機(jī)的,所以平行鍵斷裂擴(kuò)展到臨水面時(shí)并不是一個(gè)平面,同時(shí)考慮應(yīng)力波的反射作用,最先破裂的部位在炮孔的中部,該部位的塊體顆粒獲得的動(dòng)能較多,所以拋擲的速度也最大,而臺(tái)階底部和上部的顆粒平行鍵的斷裂滯后,特別是臺(tái)階上部以及兩側(cè)殘存有大塊。

圖5 水下臺(tái)階爆破塊體運(yùn)動(dòng)過(guò)程Fig. 5 Movement process of blocks underwater bench blasting block

從塊體運(yùn)動(dòng)來(lái)看,在0.0204 s時(shí)刻,炮孔附近的塊體顆粒速度明顯大于兩側(cè)靠近巖體壁面的顆粒,說(shuō)明鼓包運(yùn)動(dòng)由最小抵抗線方向開(kāi)始。在0.492 s時(shí)刻,表面塊體顆粒拋擲的速度峰值為2.48 m/s,隨著顆粒進(jìn)一步剝離表面,峰值速度持續(xù)增加。在0.918 s時(shí),中部爆破區(qū)域呈現(xiàn)塌陷趨勢(shì),而兩側(cè)未見(jiàn)崩壞,部分顆粒在重力與浮力的作用下進(jìn)行拋擲運(yùn)動(dòng),并在3 s時(shí)刻顆粒拋擲速度峰值達(dá)到4.5 m/s,隨后速度峰值逐漸衰減。在2.004 s時(shí)刻,可以明顯觀察到塊體顆粒群形成了具有一定坡度的堆積體,此時(shí)塊體顆粒的運(yùn)動(dòng)形式滾落為主。隨著塊體顆粒的繼續(xù)滾動(dòng),塊體顆粒在臺(tái)階前表面形成了平緩的爆堆,并伴隨有左側(cè)與堵塞段大塊堆積在上面。

3.2 堆積范圍和分層結(jié)果分析

為了方便觀察臺(tái)階前方塊體顆粒的堆積層次分布規(guī)律,將臺(tái)階前堆積的塊體顆粒群選取并進(jìn)行上色,底部顆粒為藍(lán)色、中部顆粒為紅色、上部顆粒為綠色,見(jiàn)圖6(a)。圖6(b)為水下臺(tái)階爆堆中間截面最終堆積形狀的側(cè)視圖。通過(guò)側(cè)視圖測(cè)出爆堆堆積形狀的休止角約為20°。爆后塊體顆粒的分層延續(xù)了其在原臺(tái)階模型中的分布形式,塊體爆堆中藍(lán)色塊體顆粒位于最下層、紅色塊體顆粒位于中間、綠色塊體顆粒位于上層,且中部紅色塊體顆粒運(yùn)動(dòng)距離最遠(yuǎn),中部紅色塊體顆粒覆蓋了大部分的底部藍(lán)色塊體顆粒,主要原因是塊體顆粒堆積分層與其初始速度和位置高度有關(guān),臺(tái)階中部塊體顆粒相對(duì)位置較高,在爆炸荷載作用下水平運(yùn)動(dòng)距離最遠(yuǎn),而上部堵塞部分的塊體顆粒受爆炸荷載作用小,其主要在重力與浮力作用下翻滾覆蓋在中部塊體顆粒層之上。

圖6 塊體顆粒堆積與分層效果Fig. 6 Range of accumulation and layered accumulation of blocks

3.3 塊體運(yùn)動(dòng)速度和流體流場(chǎng)分析

為更具體地分析塊體顆粒的堆積規(guī)律,在臺(tái)階前表面中軸線上選取6個(gè)測(cè)點(diǎn),測(cè)點(diǎn)間距為2 m,每個(gè)測(cè)點(diǎn)均包含若干塊體顆粒,見(jiàn)圖3(b)所示,繪制顆粒的平均速度隨時(shí)間變化的時(shí)程圖,見(jiàn)圖7,該速度方向垂直于臺(tái)階前表面。在炸藥反向起爆之后,炸藥起爆點(diǎn)位置處的2#和3#測(cè)點(diǎn)水平速度峰值最大,分別為3.4 m/s和3.0 m/s;位于底部的1#測(cè)點(diǎn)受到顆粒下沉碰撞影響,其峰值速度分別為2.1 m/s;臺(tái)階堵塞段的5#和6#測(cè)點(diǎn)的水平速度峰值曲線加速段斜率較低且到達(dá)最大速度峰值出現(xiàn)較晚,說(shuō)明臺(tái)階堵塞段塊體顆粒受爆炸荷載作用的強(qiáng)度較小,另外也說(shuō)明該部分塊體滾動(dòng)滑落運(yùn)動(dòng)明顯。隨著塊體顆粒的運(yùn)動(dòng),顆粒的速度受水的阻力而衰減,在鼓包運(yùn)動(dòng)結(jié)束之后,塊體顆粒以滾動(dòng)滑落的運(yùn)動(dòng)形式為主。結(jié)合爆堆分層規(guī)律,并觀察速度時(shí)程曲線后半段可知,塊體顆粒的滾動(dòng)滑落速度與其所處位置有關(guān),位于底部的1#測(cè)點(diǎn)的塊體顆粒由于受到底面的摩擦阻礙以及上層顆粒的覆蓋作用,其塊體顆粒的滑動(dòng)速度較小,而4#～6#測(cè)點(diǎn)的滾動(dòng)滑落速度隨著高度的增加而衰減較慢。

圖7 測(cè)點(diǎn)水平速度峰值時(shí)程曲線Fig. 7 Time history curve of horizontal velocity peak at measuring points

圖8為臺(tái)階模型中軸面的流線分布圖。起爆初期爆炸能量主要是破碎巖石,水體無(wú)明顯流速,隨著塊體顆粒的鼓包拋擲運(yùn)動(dòng),臺(tái)階前表面附近的流體流線速度最大。在2.0 s時(shí)刻臺(tái)階上部塊體顆粒潰散滑落,導(dǎo)致了其表面附近出現(xiàn)流線回流。隨著塊體顆粒的滾動(dòng)滑落形成具有平緩坡度的爆堆,流線的回流形式進(jìn)一步擴(kuò)大,但是流線的峰值速度在逐漸衰減。當(dāng)塊體顆粒開(kāi)始滾動(dòng)滑落時(shí),流體流線則主要受到顆粒的運(yùn)動(dòng)形式影響,流線峰值速度區(qū)域均出現(xiàn)在塊體顆粒附近。

圖8 流線分布圖Fig. 8 Streamline distribution diagram

4 結(jié)論

基于Fluent-EDEM耦合方法并采用平行鍵與Hertz-Mindlin接觸模型模擬了水下巖體破碎和運(yùn)動(dòng)的力學(xué)作用機(jī)理,討論了爆破塊體拋擲速度、塊體分層堆積和爆區(qū)流場(chǎng)特點(diǎn),揭示了水下巖體爆破破碎塊體堆積的運(yùn)動(dòng)特點(diǎn)。初步結(jié)論如下:

(1)水下臺(tái)階爆破塊體顆粒堆積分層與其初始速度和位置高度有關(guān),其本質(zhì)是爆破能量密度分布與塊體重力勢(shì)能的外在表現(xiàn)。

(2)起爆初期爆炸能量主要是破碎巖石,水體無(wú)明顯流速,隨著塊體顆粒的鼓包拋擲運(yùn)動(dòng),最小抵抗線方向流體流線速度最大。當(dāng)塊體顆粒開(kāi)始滾動(dòng)滑落時(shí),流線則主要受到顆粒的運(yùn)動(dòng)形式影響。

(3)采用平行鍵與Hertz-Mindlin接觸模型能反映巖體爆破破碎過(guò)程,揭示大塊產(chǎn)生的機(jī)理與范圍。計(jì)算結(jié)果表明FLUENT-EDEM流固耦合的數(shù)值模型能較準(zhǔn)確地模擬水下臺(tái)階爆破巖體破碎與塊體運(yùn)動(dòng)問(wèn)題,可為水下巖體爆破工程提供新的研究手段。