周文哲，徐立君，郭錦澤，劉籌資，程永舟,3，鄭長(zhǎng)青

(1.長(zhǎng)沙理工大學(xué) 水利與環(huán)境工程學(xué)院，長(zhǎng)沙 410114；2.湖南水利水電職業(yè)技術(shù)學(xué)院，長(zhǎng)沙 410100；3.水沙科學(xué)與水災(zāi)害防治湖南省重點(diǎn)實(shí)驗(yàn)室，長(zhǎng)沙 410114；4.珠海爆破新技術(shù)開發(fā)公司，珠海 519099)

水下爆破在現(xiàn)階段是一種快速有效的水下施工技術(shù)，在水下作業(yè)的很多方面都有著不可替代的作用，有著廣泛的應(yīng)用[1]。水下爆炸形成的氣泡脈動(dòng)與自由液面相互作用產(chǎn)生形態(tài)不同的水冢[2]，水冢在重力的作用下回落到水面會(huì)對(duì)周圍的結(jié)構(gòu)等物體產(chǎn)生較大的沖擊。當(dāng)炸藥的裝藥量過(guò)大，形成的水下波動(dòng)會(huì)對(duì)周圍的水生動(dòng)植物產(chǎn)生沖擊[3]，還會(huì)產(chǎn)生較大的水面波動(dòng)，威脅周圍船舶、建筑安全。因此研究水下爆炸氣泡運(yùn)動(dòng)對(duì)水面形態(tài)的影響和水面興波規(guī)律，對(duì)指導(dǎo)水下爆炸實(shí)施時(shí)周圍建筑物的安全防護(hù)具有重要意義。

在水下爆炸氣泡運(yùn)動(dòng)對(duì)自由液面影響方面，Szymezak等基于BEM方法[4]，對(duì)水下爆炸氣泡運(yùn)動(dòng)進(jìn)行模擬，將水面的變化的過(guò)程分為空穴形成、空穴塌落、氣泡上浮出水、氣泡射流和水柱回落五個(gè)階段；張阿漫等基于BEM方法[5]，研究在不同距離、浮力、強(qiáng)度參數(shù)下氣泡的脈動(dòng)過(guò)程，對(duì)比各過(guò)程下的自由面的水?，F(xiàn)象。李健等采用等效初始?xì)馀萏娲ㄋ幍姆绞絒6,7]，研究了近自由液面時(shí)的氣泡膨脹過(guò)程和相對(duì)應(yīng)的自由液面變化過(guò)程，研究了距離參數(shù)對(duì)自由面水冢變化的影響規(guī)律。Li等[8]、陳瑩玉等利用物質(zhì)歐拉求解方法及流固耦合技術(shù)[9]，研究氣泡運(yùn)動(dòng)形態(tài)對(duì)水冢寬度與高都的影響。李帥等對(duì)皇冠型水冢的形成過(guò)程進(jìn)行了模擬[10]，研究了水深距離、韋博數(shù)對(duì)皇冠型水冢特征參數(shù)的影響。Safiyari等基于SPH方法[11,12]，對(duì)爆炸氣泡與自由液面的相互作用過(guò)程簡(jiǎn)化，將自由液面變化過(guò)程簡(jiǎn)化為空腔的擴(kuò)張、坍縮，進(jìn)而引起水面波動(dòng)的過(guò)程，研究爆炸后自由液面的波動(dòng)變化。李汪諱等將爆炸后自由液面的變化簡(jiǎn)化為錐形水柱的沖擊水面的過(guò)程[13]。Wang等將氣泡運(yùn)動(dòng)在水面附近形成的水冢簡(jiǎn)化為水下空腔和水面水柱兩種形態(tài)[14]，對(duì)水面變化形態(tài)進(jìn)行研究。

前人在研究氣泡與自由液面相互作用時(shí)，大多考慮水深距離參數(shù)的對(duì)氣泡脈動(dòng)以及水冢形態(tài)的影響，較少考慮裝藥量的不同；在自由液面形態(tài)變化方面，主要關(guān)注于氣泡運(yùn)動(dòng)對(duì)水冢形態(tài)的影響，在對(duì)水冢產(chǎn)生后的水面波動(dòng)過(guò)程研究較少，對(duì)水面興波過(guò)程的研究基本建立在水冢形態(tài)簡(jiǎn)化條件下。在實(shí)際工程中，水下爆破后水面將會(huì)產(chǎn)生一個(gè)持續(xù)的波動(dòng)過(guò)程，且裝藥量的變化對(duì)水面運(yùn)動(dòng)形態(tài)的影響較大?；诖?，采用ANSYS非線性有限元軟件，基于ALE算法，結(jié)合水冢形成過(guò)程及其特征幾何參數(shù)、氣泡運(yùn)動(dòng)特性與水面波動(dòng)過(guò)程，就不同裝藥量條件下，水下爆炸氣泡運(yùn)動(dòng)對(duì)水面形態(tài)的影響機(jī)理開展研究。

1 水下爆炸氣泡運(yùn)動(dòng)數(shù)值模擬

1.1 基于ALE方法的氣泡運(yùn)動(dòng)模擬

采用ALE(Arbitrary Lagrangian-Eulerian)算法對(duì)氣泡運(yùn)動(dòng)進(jìn)行三維數(shù)值模擬，ALE算法具有在保持物質(zhì)界面的同時(shí)防止網(wǎng)格過(guò)度變形造成計(jì)算失效的特點(diǎn)，該算法兼顧Lagrange算法與Euler算法的優(yōu)點(diǎn)，在模擬爆炸產(chǎn)生氣體的大變形現(xiàn)象的同時(shí)，保持對(duì)不同物質(zhì)之間界面的捕捉，能較好地模擬水下爆炸過(guò)程中流體的大范圍變形與氣泡、水流的高速流動(dòng)。

在ANSYS軟件中進(jìn)行數(shù)值計(jì)算時(shí)，為了提高計(jì)算效率同時(shí)保證模擬結(jié)果的正確性，在氣泡最大半徑1.5倍范圍內(nèi)采用與裝藥半徑尺寸相同的均勻網(wǎng)格，其他區(qū)域采用由均勻網(wǎng)格區(qū)域至流域外邊界逐漸漸變的網(wǎng)格。為了兼顧計(jì)算效率，采用三維1/4模型，在對(duì)稱面上施加對(duì)稱邊界條件，在空氣域和水域外邊界處設(shè)置pressure outflow邊界單元層，以模擬無(wú)限邊界條件，模型的x/y方向和裝藥下方的水域深度應(yīng)至少為氣泡脈動(dòng)最大半徑的4倍，空氣層的高度滿足水冢中心水柱向上發(fā)展的要求。幾何模型如圖1。

炸藥使用*MAT_HIGH_EXPLOSIVE_BURN材料模型，采用JWL狀態(tài)方程

(1)

式中：p為爆轟產(chǎn)物壓強(qiáng),Pa；V相對(duì)體積,V=v/v0，v0為初始相對(duì)體積；e為單位體積內(nèi)能,J/m3；A、B、R1、R2、ω為狀態(tài)方程參數(shù)，見表1。

圖 1 模型示意圖Fig. 1 Schematic diagram of the model

水使用*MAT_NULL材料模型，采用Grüneisen狀態(tài)方程

(γ0+αμ)E

(2)

式中：p為壓強(qiáng),Pa；ρ0為水體密度,kg/m3；C為水中聲速；E單位體積內(nèi)能,J/m3；V為相對(duì)體積；γ0、S1、S2、S3為GRUNEISEN狀態(tài)方程參數(shù)，α為γ0的一階體積系數(shù)；Pc為截?cái)鄩簭?qiáng)，Pa，水的截?cái)鄩簭?qiáng)為飽和蒸氣壓；η水的粘滯系數(shù)，Pa·S。水的材料方程與狀態(tài)方程參數(shù)見表2。

空氣使用*MAT_NULL材料模型，采用線性多項(xiàng)式狀態(tài)方程

p=C0+C1μ+C2μ2+C3μ3+(C4+C5μ+C6μ2)E

(3)

(4)

表 1 TNT炸藥的材料方程和狀態(tài)方程參數(shù)取值[15,16]

表 2 水的材料方程和狀態(tài)方程參數(shù)取值[15,16]

對(duì)于理想氣體有下列關(guān)系式

C0=C1=C2=C3=C6=0,C4=C5=γ-1

(5)

(6)

式中：p為壓強(qiáng),Pa；γ為單位熱值率，對(duì)于理想氣體γ=1.4；E0為單位體積初始內(nèi)能；C0、C1、C2、C3、C4、C5、C6為線性多項(xiàng)式方程參數(shù)?？諝獾牟牧戏匠膛c狀態(tài)方程參數(shù)見表3。

表 3 空氣的材料方程和狀態(tài)方程參數(shù)取值[15,16]

1.2 氣泡運(yùn)動(dòng)模型的驗(yàn)證

為了驗(yàn)證氣泡運(yùn)動(dòng)模型的有效性，設(shè)置與文獻(xiàn)[17]實(shí)驗(yàn)相同的工況，使用5.2 g TNT的球形裝藥(與實(shí)驗(yàn)中4 g PENT炸藥等效)，裝藥半徑為0.009 m，模型計(jì)算尺寸為1 m×1 m×4 m，漸變網(wǎng)格劃分尺寸最小為0.009 m，最大為0.1 m，均勻網(wǎng)格處網(wǎng)格尺寸為0.009 m。如圖2所示，就文獻(xiàn)中近水面條件(裝藥位置位于水面以下0.13 m)下的氣泡運(yùn)動(dòng)趨勢(shì)相對(duì)比，發(fā)現(xiàn)實(shí)驗(yàn)結(jié)果和數(shù)值模擬結(jié)果兩者在氣泡脈動(dòng)的模擬結(jié)果，一致性較好。

對(duì)于TNT炸藥水下爆炸氣泡脈動(dòng)的第一周期及氣泡半徑計(jì)算，已經(jīng)形成了能夠應(yīng)用于工程實(shí)際的經(jīng)驗(yàn)公式[17]。如表4所示，通將數(shù)值模擬結(jié)果與經(jīng)驗(yàn)公式計(jì)算值進(jìn)行對(duì)比，發(fā)現(xiàn)氣泡半徑及氣泡周期變化誤差在10%以內(nèi)，表明數(shù)值模擬結(jié)果可靠。

圖 2 近自由液面條件下氣泡運(yùn)動(dòng)趨勢(shì)對(duì)比Fig. 2 Comparison of bubble movement trends near free surface of water

(7)

(8)

式中:R為氣泡第一次脈動(dòng)的氣泡最大半徑,m；T為脈動(dòng)周期,s；W為裝藥當(dāng)量,kg；D為裝藥位置處水深,m。

表 4 氣泡最大半徑與周期誤差計(jì)算結(jié)果

2 氣泡運(yùn)動(dòng)對(duì)水面形態(tài)的影響

2.1 不同類型水冢形成過(guò)程與水面形態(tài)

不同水?，F(xiàn)象的發(fā)生取決于起爆點(diǎn)水深距離參數(shù)γf=h/Rmax[18]，式中：h為炸藥裝藥位置距水面的距離,m；Rmax為氣泡脈動(dòng)過(guò)程的最大半徑,m。那立民等通過(guò)大量工況計(jì)算[19]，總結(jié)出水冢形態(tài)與距離參數(shù)γf的關(guān)系，將水冢形態(tài)劃分為零碎型(0≤γf≤0.06)、飛濺型(0.06≤γf≤0.03)、酒杯型(0.3≤γf≤0.7)、皇冠型(0.7≤γf≤1.4)等。為研究不同水深距離參數(shù)下水冢的形成過(guò)程，選用5 kg球型TNT炸藥，控制炸藥質(zhì)量不變，改變炸藥中心至水面距離，對(duì)四種特征水冢的形成進(jìn)行模擬。模擬參數(shù)與結(jié)果如表5。

表 5 不同水深距離參數(shù)γf下的水冢形態(tài)

如圖3所示，5 kg球形TNT炸藥在水下爆炸的模擬結(jié)果表明，隨著水深距離參數(shù)的不斷增大，依次產(chǎn)生了零碎型、飛濺型、酒杯型及皇冠型等不同類型的水冢，且水冢所產(chǎn)生的空腔最大半徑也在逐步增大，在酒杯型水冢空腔最大半徑達(dá)到最大。

2.2 不同裝藥量對(duì)水冢形態(tài)的影響

皇冠型水冢形態(tài)較為穩(wěn)定，以形成皇冠型水冢的特征水深距離參數(shù)γf=1作為基準(zhǔn)參數(shù)，對(duì)不同裝藥量下水冢形態(tài)發(fā)展進(jìn)行研究。表6所示，在0.5 kg、5 kg、50 kg裝藥量條件下，爆炸氣泡脈動(dòng)第一周期時(shí)間分別為186 ms、318 ms、665 ms。氣泡第一周期時(shí)間與裝藥水深相關(guān)[17]，在同一距離水深參數(shù)條件下，不同裝藥量產(chǎn)生的氣泡半徑大小不同，需要的裝藥水深不同，氣泡第一周期脈動(dòng)時(shí)長(zhǎng)不同。隨著裝藥量增大，氣泡脈動(dòng)第一周期時(shí)長(zhǎng)增加。

圖 3 不同水深距離參數(shù)γf下的水冢形態(tài)Fig. 3 Water mound morphology under different water depth and distance parameters γf

表 6 不同裝藥量下氣泡特征

為了對(duì)水冢形成過(guò)程進(jìn)行直觀對(duì)比，以當(dāng)前時(shí)刻/氣泡脈動(dòng)第一周期時(shí)長(zhǎng)對(duì)時(shí)間單位進(jìn)行歸一化處理。5 kg球形TNT炸藥，在水深距離參數(shù)γf=0.998條件下將產(chǎn)生皇冠型水冢，如圖4所示。水冢形態(tài)發(fā)展過(guò)程為：0.483時(shí)刻，氣泡膨脹至最大半徑，水面被氣泡頂起；0.9時(shí)刻，由于氣泡坍塌，形成向下氣泡射流，并向上形成水射流，與氣泡頂起的水體一起組成了中心水柱；1.317時(shí)刻，隨著氣泡射流發(fā)展的同時(shí)，向上水射流也在不斷發(fā)展，并促進(jìn)了中心水柱的發(fā)展。環(huán)形水體也在氣泡膨脹的推動(dòng)下開始突出水面向上發(fā)展，并與中心水柱一起構(gòu)成了皇冠型水冢；2.4時(shí)刻，氣泡經(jīng)過(guò)多次脈動(dòng)分裂成兩個(gè)部分，水冢也在多次氣泡脈動(dòng)作用下形成了位置不同的環(huán)狀水體；2.65時(shí)刻后，水面水柱達(dá)到最大高度，此時(shí)四周的水冢已經(jīng)在重力作用下下墜，沖擊水面并形成波動(dòng)向外傳播。

圖5、圖6分別為0.5 kg、50 kg裝藥量下水冢發(fā)展過(guò)程。在歸一化1時(shí)刻，自由液面均形成中心水柱(圖5d、圖6d)，隨后中心水柱與環(huán)形水體發(fā)展成為冠狀結(jié)構(gòu)(圖5e、圖6e)，該過(guò)程與5 kg裝藥量下皇冠型水冢的形成過(guò)程類似。對(duì)比圖4、圖5、圖6可知，在相同水深距離參數(shù)條件下，裝藥量的增大并不會(huì)改變水冢的基本形態(tài)。

圖 4 5 kg裝藥量，γf=0.998條件下水冢形態(tài)發(fā)展過(guò)程(a)～(i)為歸一化時(shí)間Fig. 4 The development process of water spike morphology under the condition of γf=0.998 and charge amount of 5 kg((a)～(i) are normalized time)

圖 5 0.5 kg裝藥量，γf=1.003條件下水冢形態(tài)發(fā)展過(guò)程(a)～(i)為歸一化時(shí)間Fig. 5 The development process of water spike morphology under the condition of γf=1.003 and charge amount of 0.5 kg((a)～(i) are normalized time)

圖 6 50 kg裝藥量，γf=0.999條件下水冢形態(tài)發(fā)展過(guò)程(a)～(i)為歸一化時(shí)間Fig. 6 The development process of water spike morphology under the condition of γf=0.999 and charge amount of 50 kg((a)～(i) are normalized time)

對(duì)比三種裝藥量下的水冢形成過(guò)程。在氣泡脈動(dòng)第一周期收縮階段，0.5 kg、5 kg裝藥量下形成的收縮氣泡呈現(xiàn)果仁狀，產(chǎn)生向下的射流，同時(shí)水面形成向上的水射流，對(duì)應(yīng)的歸一化時(shí)刻分別為0.935 s(圖5c)、0.900(圖4c)；50 kg裝藥量下在氣泡收縮過(guò)程中，在氣泡的上表面形成了向下射流，下表面形成了向上射流，對(duì)應(yīng)時(shí)刻為0.937 s(圖6c)。在水冢冠狀結(jié)構(gòu)發(fā)展過(guò)程中，0.5 kg、5 kg裝藥量下形成多個(gè)不同位置的環(huán)狀水體，對(duì)應(yīng)時(shí)刻為2.871 s(圖5f)、2.400(圖4f)；50 kg裝藥量下形成的環(huán)狀結(jié)構(gòu)較為破碎，對(duì)應(yīng)時(shí)刻為2.424 s(圖6f)。

2.3 不同裝藥量對(duì)氣泡運(yùn)動(dòng)的影響

對(duì)比圖4、圖5、圖6，裝藥量增大使爆炸形成的氣泡體積增大，氣泡運(yùn)動(dòng)趨勢(shì)發(fā)生改變。對(duì)不同裝藥量下的氣泡運(yùn)動(dòng)特征進(jìn)行對(duì)比分析，圖7(a)為在歸一化時(shí)間下的氣泡時(shí)程變化曲線，氣泡最大半徑隨著裝藥增大而增大。圖7(b)為在歸一化時(shí)間下氣泡的無(wú)量綱位移(氣泡中心位移/氣泡最大半徑)歷程曲線，在第一周期膨脹階段，三種裝藥量條件下的氣泡中心均向自由面移動(dòng)。在第一周期收縮階段，0.5 kg裝藥量下氣泡中心的向遠(yuǎn)離自由面方向移動(dòng)；5 kg裝藥量下氣泡中心位移基本穩(wěn)定在初始位置附近；50 kg裝藥量下氣泡中心始終向自由面移動(dòng)。在張阿漫等對(duì)近自由面氣泡運(yùn)動(dòng)規(guī)律研究中[5]，無(wú)量綱水深距離γf=1的條件下，氣泡中心位移變化與浮力參數(shù)相關(guān)?？梢娫谙嗤罹嚯x參數(shù)下，裝藥量的增大會(huì)使氣泡最大半徑增大，引起浮力參數(shù)的變化，從而產(chǎn)生不同的射流情況。在0.5 kg、5 kg裝藥量情況產(chǎn)生的氣泡體積較小，自由面的Bjerknes作用力占主導(dǎo)地位，氣泡在收縮階段形成遠(yuǎn)離自由面方向的射流；而在50 kg裝藥量情況下，氣泡體積較大受到的浮力作用增強(qiáng)，在收縮階段氣泡表面和底端均形成射流。

圖 7 不同裝藥量氣泡運(yùn)動(dòng)特征對(duì)比Fig. 7 Comparison of bubble motion characteristics with different charges

2.4 不同裝藥量下的水冢特征參數(shù)及水面波動(dòng)對(duì)比

冠型水冢特征寬度為氣泡脈動(dòng)第二周期收縮至最小時(shí)刻形成的水冢底部寬度，特征高度為中心水柱最高點(diǎn)距自由液面高程的距離。對(duì)0.5 kg、5 kg、50 kg三種裝藥量下的計(jì)算結(jié)果進(jìn)行了比較，如圖8所示。在0.5 kg時(shí)水冢寬度為2.428 m，中心水柱高度為2.315 m；5 kg時(shí)水冢寬度為4.978 m，中心水柱高度為4.437 m；50 kg時(shí)水冢寬度為11.538 m，中心水柱高度為10.179 m。隨著裝藥量的不斷增大，水冢寬度、中心水柱高度逐漸增大。

圖 8 不同裝藥量下水冢特征尺寸對(duì)比Fig. 8 Comparison of characteristic sizes of water mounds with different charges

在最大氣泡半徑2倍位置處設(shè)置監(jiān)測(cè)點(diǎn)，對(duì)0.5 kg、5 kg、50 kg裝藥量下的水面波動(dòng)特性進(jìn)行分析，如圖9所示。0.5 kg、5 kg、50 kg裝藥量在監(jiān)測(cè)點(diǎn)處形成的最大水面波動(dòng)高度分別為40.978 cm、53.636 cm、153.684 cm，隨著裝藥量的增加，產(chǎn)生的水面波動(dòng)在增大。對(duì)于50 kg裝藥量，由于氣泡的上浮效應(yīng)，能夠托動(dòng)更多的水體上升，沖擊水面形成的水面波動(dòng)最大高度也相對(duì)較大。在2時(shí)刻以內(nèi)，水面波動(dòng)由爆炸初期氣泡頂起的水體產(chǎn)生，氣泡半徑將影響水面波動(dòng)的大小。隨著中心水柱以及環(huán)形水體的發(fā)展，監(jiān)測(cè)點(diǎn)處自由液面先上升后下降，并隨著中心水柱及環(huán)形水體的下墜沖擊，出現(xiàn)波動(dòng)狀況。從圖9中可看出，隨著裝藥量增大，達(dá)到最大波動(dòng)的時(shí)間在不斷提前。經(jīng)計(jì)算，水面波動(dòng)從水域中心傳播到2倍氣泡半徑處的平均速度分別為1.985 m/s、2.894 m/s、3.533 m/s，由于水面中心水柱和環(huán)狀水體所具有的勢(shì)能大小不同，從而影響了水面波動(dòng)到達(dá)測(cè)點(diǎn)的時(shí)間。隨著裝藥量的增加，水面波動(dòng)高度在逐漸增大，且水面波動(dòng)達(dá)到最大的時(shí)刻也在不斷提前。

3 結(jié)論

通過(guò)ANSYS非線性有限元程序?qū)λ卤馀葸\(yùn)動(dòng)過(guò)程進(jìn)行了模擬，研究水下爆炸氣泡運(yùn)動(dòng)對(duì)自由液面形態(tài)變化的影響，得到如下結(jié)論：

圖 9 不同裝藥量條件下水面波動(dòng)歷時(shí)曲線Fig. 9 The water surface fluctuation duration curve under the condition of different charge amount

(1)相同裝藥量下，隨著水深距離參數(shù)的增大，自由液面由破碎狀態(tài)逐漸形成為具有氣泡脈動(dòng)周期的皇冠型水冢，爆炸所產(chǎn)生的空腔半徑在逐步增大。

(2)在相同水深距離參數(shù)下，裝藥量大小變化對(duì)水?；拘螒B(tài)影響較小，水冢類型主要受水深距離參數(shù)影響。爆炸產(chǎn)生的氣泡最大半徑受裝藥量影響，不同半徑氣泡在第一周期收縮階段產(chǎn)生的射流方向不同，氣泡中心點(diǎn)在豎直方向無(wú)量綱位移趨勢(shì)不同，隨著裝藥量增大，逐漸向自由液面方向移動(dòng)。

(3)隨著裝藥量增加，水冢特征寬度、特征高度的變化呈增大趨勢(shì)。爆炸產(chǎn)生的水面波動(dòng)高度與氣泡半徑大小相關(guān)，水面波動(dòng)速度與波動(dòng)高度成正比。