基于Fluent-EDEM流固耦合算法的水下臺階爆破爆堆形成過程研究*

吳 亮,余 創,梁志堅,段衛東,陳 明

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

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

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

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

1 計算方法

1.1 流體基本原理

計算流體力學是利用數學的計算方法來實現對流場中的控制方程進行計算,并對有限元網格節點進行數值求解。當研究對象為不可壓縮流體時,基本的定律有三個:質量守恒定律、動量守恒定律、能量守恒定律,如下式(1)和式(2)。在本文中沒有能量的交換,故不考慮能量守恒方程。

(1)

(2)

1.2 DEM顆粒接觸模型

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

(3)

(4)

式中:v、ω、m分別為顆粒速度、角速度和質量,下標i,j表示兩個相互接觸的顆粒;Fij為i和j顆粒間的接觸力;Fg為顆粒間的引力;Ff為流體與顆粒間的相互作用力;I是顆粒轉動慣量;T是顆粒間的力矩;Tf是流體引起的扭矩。

本文采用Hertz-Mindlin接觸模型計算DEM顆粒之間的作用力,該模型具備滾動摩擦特性,如圖1所示。Hertz-Mindlin接觸模型的計算表達式見文獻[19]。

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

1.3 DEM顆粒平行鍵模型

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

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

在平行鍵生成之后顆粒上的力和扭矩被重置為0,在模型開始計算時則根據公式(5)～(8)在每個時間步長中計算平行鍵受到的力和扭矩的增量。

δFn=-vnSnAδt

(5)

δFt=-vtStAδt

(6)

δMn=-ωnStJδt

(7)

(8)

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

(9)

(10)

式中:RB為平行鍵半徑,是平行鍵合模型輸入參數之一,取值一般大于顆粒半徑。

當法向力或切向力超過其限定值時,顆粒間的平行鍵發生斷裂

(11)

(12)

式中:σmax、τmax為最大法向應力和最大切向應力。

2 水下臺階爆破計算模型

2.1 計算模型

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

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

2.2 模型參數

在EDEM模型中引入平行鍵的模型來模擬連續的巖體介質,平行鍵的模型除了包含摩擦系數在內的基本力學參數外還需要額外輸入5個平行鍵參數。由于平行鍵的模型能體現微觀上顆粒之間的力學行為,除了平行鍵半徑取顆粒半徑值外,至今的研究結果還不能明確地給出實驗中獲得的參數與平行鍵參數之間的關系式,因此為了獲得準確的數值計算參數,需對平行鍵的參數進行標定。許多學者使用DEM進行巴西劈裂數值計算[20-22],通過與實驗結果對比,標定平行鍵的參數。

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

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

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

綜上,結合標定后的平行鍵模型參數,將堆積的顆粒建立平行鍵,臺階模型共生成平行鍵351 460個,具體參數如下表1所示。需要注意的是,同一組模型中含有平行鍵與無平行鍵在時間步長的選取上存在差異,通常在進行含有平行鍵的數值模型計算時選取兩者中時間步長較小的一個[24],本文離散元的計算步長為3×10-5s,流體的計算步長是離散元的10倍,為3×10-4s。

表1 巖體EDEM模型參數Table 1 EDEM model parameters of rock mass

3 模擬結果

3.1 臺階爆破破碎與塊體運動過程

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

圖5 水下臺階爆破塊體運動過程Fig. 5 Movement process of blocks underwater bench blasting block

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

3.2 堆積范圍和分層結果分析

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

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

3.3 塊體運動速度和流體流場分析

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

圖7 測點水平速度峰值時程曲線Fig. 7 Time history curve of horizontal velocity peak at measuring points

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

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

4 結論

基于Fluent-EDEM耦合方法并采用平行鍵與Hertz-Mindlin接觸模型模擬了水下巖體破碎和運動的力學作用機理,討論了爆破塊體拋擲速度、塊體分層堆積和爆區流場特點,揭示了水下巖體爆破破碎塊體堆積的運動特點。初步結論如下:

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

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

(3)采用平行鍵與Hertz-Mindlin接觸模型能反映巖體爆破破碎過程,揭示大塊產生的機理與范圍。計算結果表明FLUENT-EDEM流固耦合的數值模型能較準確地模擬水下臺階爆破巖體破碎與塊體運動問題,可為水下巖體爆破工程提供新的研究手段。