近物面爆炸氣泡潰滅與射流載荷特性研究

古 濱， 劉亮濤， 魏銘利， 邱清水， 姚熊亮， 王志凱

(1.西華大學(xué)土木建筑與環(huán)境學(xué)院， 成都 610039；2.哈爾濱工程大學(xué)船舶工程學(xué)院， 哈爾濱 150001)

引 言

水下爆炸指炸藥、魚雷、炸彈或核彈等在水中的爆炸。爆炸后在水中形成向四周擴(kuò)展并不斷減弱的沖擊波。爆炸產(chǎn)物形成的氣泡在水中膨脹，然后回縮，進(jìn)行振蕩并不斷上浮，同時向四周發(fā)出二次壓力脈沖。當(dāng)沖擊波遇到物體時發(fā)生反射、折射和繞射，物體在沖擊波和二次壓力脈沖的作用下發(fā)生位移、變形或破壞。沖擊波到達(dá)水面和氣球突出水面后，在水面激起表面波。水下爆炸已有了較為系統(tǒng)的研究，并取得了較豐碩的成果，如古濱等[1-2]通過ALE方法結(jié)合數(shù)值模擬研究了水下沖擊波作用下雙層殼結(jié)構(gòu)的響應(yīng)特征，姜忠濤等[3]研究了近場水下爆炸氣泡射流載荷沖擊船體外板的動響應(yīng)。

水下氣泡在不同邊界條件附近會呈現(xiàn)出不同的運(yùn)動特性，由于氣泡周圍邊界條件的不對稱，會使氣泡周圍壓力分布朝某個方向產(chǎn)生梯度，在坍塌過程中氣泡將會朝壓力梯度大的方向產(chǎn)生射流，并穿透氣泡。根據(jù)邊界類型的不同，射流產(chǎn)生的方向會發(fā)生改變，如在自由液面附近，射流遠(yuǎn)離自由液面，而壁面附近的氣泡將會朝向壁面產(chǎn)生射流，這種射流現(xiàn)象最早由Naude和Ellis(1961)[4]、Benjamin和Ellis(1966)[5]等在實(shí)驗(yàn)中觀測到，并被認(rèn)為是結(jié)構(gòu)空化剝蝕的主要原因。對于近物面的氣泡，除了距離對氣泡運(yùn)動特性影響顯著之外，初始參數(shù)的不同，氣泡的運(yùn)動形式和射流載荷特性也會呈現(xiàn)出不同的特征。對于爆炸氣泡等尺度較大的氣泡，浮力對氣泡的作用較為顯著。而對于空化等尺度較小的氣泡，粘性對氣泡的作用逐漸凸顯。對于存在自由液面影響的氣泡，由于物面和自由液面的耦合作用，氣泡射流又會有不一樣的表現(xiàn)形式。

盡管氣泡在不同邊界附近的運(yùn)動特性已經(jīng)被很多學(xué)者深入地研究，如王加夏等[6]通過電火花氣泡裝置和高速攝像設(shè)備對氣泡與附近剛性/柔性邊界間的耦合現(xiàn)象進(jìn)行了實(shí)驗(yàn)觀測和研究；賀銘等[7]利用歐拉有限元數(shù)值模型對近場水下爆炸氣泡與雙層破口結(jié)構(gòu)的相互作用機(jī)理進(jìn)行了研究，主要分析了艙室涌流及流場演化等物理屬性；李梅[8]等研究了雙爆源氣泡與水面的相互作用過程；朱錫、張阿漫等[9-13]都開展過小當(dāng)量真實(shí)裝藥水下爆炸產(chǎn)生氣泡的實(shí)驗(yàn)，研究氣泡脈動周期和最大半徑等參數(shù)，得到了氣泡全周期的運(yùn)動形態(tài)。但是，由于很多方面的局限性，目前對于近物面，尤其一些極近物面不同參數(shù)下的氣泡坍塌運(yùn)動和射流載荷特性的研究并不系統(tǒng)。

為此，本文基于不可壓縮流體控制方程建立氣泡脈動及射流載荷特性的計(jì)算模型。通過分析不同物面深度下射流特征參數(shù)分別與藥包質(zhì)量和無量綱距離參數(shù)的變化趨勢，開展了考慮水下爆炸氣泡射流速度、射流寬度和射流高度等方面的較為系統(tǒng)的研究，基于數(shù)值模擬與模型實(shí)驗(yàn)相結(jié)合的方法，驗(yàn)證了計(jì)算方法的有效性，并提出了3種有效的數(shù)值預(yù)報模型。

1 計(jì)算模型

1.1控制方程

假設(shè)流體絕熱，氣泡內(nèi)部和外界流場不存在質(zhì)量和熱交換。氣泡在外界流場的速度在大部分時間內(nèi)都要遠(yuǎn)小于聲速，因此外界流場的可壓縮性忽略不計(jì)。氣泡運(yùn)動之前，假設(shè)流體靜止，隨后氣泡的運(yùn)動被認(rèn)為是軸對稱的，如圖1所示，因此流體的運(yùn)動由不可壓縮、粘性的Navier-Stokes方程控制，柱坐標(biāo)下的Navier-Stokes方程有如下表現(xiàn)形式：

(1)

(2)

(3)

式中，

其中:x和y分別為軸對稱坐標(biāo)軸r和z；u和v分別表示坐標(biāo)軸x和y方向的速度；ρ1表示流體密度；p表示流體壓力；gy表示沿y方向的重力加速度；υ表示運(yùn)動粘性系數(shù)。當(dāng)式(1)～式(3)中的粘性項(xiàng)不考慮時，方程變?yōu)椴豢蓧嚎sEuler方程。

圖1 軸對稱坐標(biāo)示意圖

1.2氣體狀態(tài)方程

如果直接采用流體控制方程求解氣泡內(nèi)部可壓氣體和外部不可壓液體，界面處的數(shù)值不穩(wěn)定將是個巨大的挑戰(zhàn)。為了避免出現(xiàn)這種情況，通常假設(shè)氣泡內(nèi)部壓力均勻分布，把壓力作為外界流場的壓力邊界條件。這種假設(shè)對氣泡的運(yùn)動影響較小，因此是合理的，此時氣泡內(nèi)部的氣體的狀態(tài)可以用理想氣體狀態(tài)方程表示。氣泡內(nèi)部壓力由常值蒸汽壓力和隨體積變化的非冷凝氣體壓力組成，根據(jù)絕熱假設(shè)，氣泡內(nèi)部的理想氣體狀態(tài)方程可以轉(zhuǎn)化成另一種表現(xiàn)形式：

(4)

其中：pb表示氣泡內(nèi)部壓力；pv表示飽和蒸汽壓力；p0表示氣泡初始壓強(qiáng)；V0對應(yīng)的非冷凝氣體壓力；Vt表示氣泡在任一時刻的體積；?表示氣體的比熱。

1.3邊界條件

在求解氣泡外界流場之前，需要確定流場的邊界條件。垂直于對稱軸和剛性壁面的速度分量滿足不可穿透條件，可以寫成：

v·n=0

(5)

其中，v表示流場的速度矢量，n表示邊界指向流場的單位法向速度。

無窮遠(yuǎn)處的壓力邊界條件由靜水壓力表示：

pout=p

(6)

其中，pout表示無窮遠(yuǎn)處的靜水壓力。

由于氣泡內(nèi)部氣體的運(yùn)動狀態(tài)不予以直接求解，因此需要把壓力當(dāng)成外界流場的壓力邊界條件。在氣泡界面處，切向壓力條件和法向壓力條件需要同時滿足。由于氣泡內(nèi)部氣體的粘性極小，通常忽略不計(jì)，因此切向壓力邊界條件可以寫成：

μ·τ·δ·n=0

(7)

其中：τ表示界面的單位切向量;δ表示偏應(yīng)力張量。

法向壓力邊界條件需要滿足：

(8)

其中，pl表示在外界流場一側(cè)邊界處的壓力，σ表示表面張力系數(shù)；κ表示氣泡表面的局部曲率，可以通過Young-Laplace方程進(jìn)行求解：

(9)

其中，R1和R2分別表示沿x和y軸的主曲率半徑。對于軸對稱模型氣泡表面質(zhì)點(diǎn)的位置x(s)和y(s)可以當(dāng)成弧長s的函數(shù)，因此主曲率半徑R1和R2可以表示成：

(10)

(11)

2 計(jì)算結(jié)果驗(yàn)證與分析

2.1氣泡潰滅運(yùn)動過程驗(yàn)證

實(shí)驗(yàn)驗(yàn)證的工況見表1。這里采用相同的工況對邊界積分方程進(jìn)行理論驗(yàn)證，理論驗(yàn)證主要將邊界積分方程計(jì)算的結(jié)果和Rayleigh-Plesset方程的理論值進(jìn)行對比，包括半徑曲線和氣泡表面節(jié)點(diǎn)速度。在這之前，先對邊界積分方程進(jìn)行收斂性驗(yàn)證，確定氣泡表面節(jié)點(diǎn)數(shù)量N和Δφ速度勢變化值。氣泡的初始半徑R0=0.0204，初始內(nèi)部壓力P0=66.928 MPa。圖2給出了不同速度勢變化值Δφ下的氣泡射流時歷曲線。當(dāng)Δφ取值較大時，射流速度不穩(wěn)定，出現(xiàn)明顯的振蕩現(xiàn)象。當(dāng)Δφ取值減小時，射流速度逐漸收斂。當(dāng)Δφ等于0.01時，射流速度基本收斂。因此在隨后的計(jì)算中，將速度勢變化值Δφ取為0.01。圖3中給出了氣泡表面不同節(jié)點(diǎn)數(shù)量的射流速度收斂曲線，當(dāng)節(jié)點(diǎn)速度大于60時基本收斂，隨后的計(jì)算將節(jié)點(diǎn)數(shù)量N取為60。圖4和圖5分別給出邊界積分方程數(shù)值解和Rayleigh-Plesset方程理論解兩者關(guān)于氣泡表面節(jié)點(diǎn)速度和氣泡半徑的對比曲線。由圖4和圖5可見，無論是節(jié)點(diǎn)速度曲線還是半徑曲線，邊界積分方程的數(shù)值結(jié)果和理論解都吻合得較好。

表1水下爆炸實(shí)驗(yàn)初始參數(shù)

圖2 速度勢變化值Δφ收斂性分析

圖3 氣泡表面節(jié)點(diǎn)數(shù)量收斂性分析

圖4 邊界積分方程數(shù)值解和理論解關(guān)于氣泡節(jié)點(diǎn)速度對比

采用邊界積分方程模擬的數(shù)值結(jié)果與實(shí)驗(yàn)對比情況如圖6所示，數(shù)值結(jié)果中的箭頭表示速度方向，云圖表示流場壓力。由圖6可知，除了環(huán)狀階段的凸起，數(shù)值結(jié)果與實(shí)驗(yàn)結(jié)果吻合良好。由于凸起主要由氣液混合組成，里面包含了很多破碎的小氣泡，因此很難采用邊界積分方程對氣液混合物進(jìn)行模擬，通常在計(jì)算過程中會忽略不計(jì)。因此采用邊界積分方程來研究水下爆炸氣泡射流特征參數(shù)是合適的。

圖5 邊界積分方程數(shù)值解和理論解關(guān)于氣泡半徑對比

圖6 水下爆炸氣泡數(shù)值與實(shí)驗(yàn)對比

2.2氣泡射流載荷計(jì)算

2.2.1 射流速度模型

氣泡射流與無量綱距離參數(shù)的關(guān)系可以采用多項(xiàng)式來表示，通過比較一系列不同的數(shù)學(xué)模型，最后選用以下數(shù)學(xué)模型來描述射流速度V與藥包質(zhì)量和無量綱距離參數(shù)之間的關(guān)系：

(12)

其中，ai(i=1,2,3,4,5)表示未知系數(shù)。不同物面深度下未知系數(shù)ai的取值情況見表2。

表2不同物面深度下的射流速度系數(shù)

圖7給出了物面深度等于50 m和250 m時射流速度數(shù)學(xué)模型與邊界元數(shù)值結(jié)果的整體對比。從圖中可以看出，邊界元數(shù)值結(jié)果散落在射流速度數(shù)學(xué)模型周圍，通過計(jì)算，兩者之間的最大誤差不超過2%。

圖7 射流速度數(shù)學(xué)模型與邊界元數(shù)值結(jié)果的整體對比

為了進(jìn)一步驗(yàn)證式(12)中的射流速度數(shù)學(xué)模型和表2中系數(shù)的準(zhǔn)確性，對射流速度數(shù)學(xué)模型與邊界元數(shù)值結(jié)果的進(jìn)行局部對比分析。圖8給出了物面深度h等于50 m和250 m時射流速度數(shù)學(xué)模型與邊界元數(shù)值結(jié)果隨著藥包質(zhì)量變化的局部對比，其中無量綱距離參數(shù)分別等于0.65、1.0和1.5。無論從曲線的變化趨勢上，還是從結(jié)果的吻合程度上來看，射流速度數(shù)學(xué)模型在藥包質(zhì)量這個變化方向上的結(jié)果和邊界元數(shù)值結(jié)果吻合較為一致。圖9給出了物面深度h=50 m和250 m時射流速度數(shù)學(xué)模型與邊界元數(shù)值結(jié)果隨著無量綱距離參數(shù)變化的局部對比，其中藥包質(zhì)量W分別等于50 kg, 500 kg和1000 kg。從圖9中可以看出，邊界元的數(shù)值結(jié)果與射流速度數(shù)學(xué)模型計(jì)算的結(jié)果吻合較為一致。即使射流速度在無量綱距離參數(shù)γ=0.7存在拐點(diǎn)，射流速度數(shù)學(xué)模型也能對邊界元數(shù)值結(jié)果進(jìn)行完美表示，如圖9(b)所示。

圖8 射流速度數(shù)學(xué)模型與邊界元數(shù)值結(jié)果隨著藥包質(zhì)量變化的局部對比

圖9 射流速度數(shù)學(xué)模型與邊界元數(shù)值結(jié)果隨著無量綱距離參數(shù)變化的局部對比

2.2.2 射流寬度模型

射流寬度與無量綱距離參數(shù)的關(guān)系采用多項(xiàng)式進(jìn)行表示。通過對多種數(shù)學(xué)函數(shù)模型進(jìn)行比較，選取以下數(shù)學(xué)模型來描述射流寬度D與藥包質(zhì)量和無量綱距離參數(shù)之間的關(guān)系：

肉牛養(yǎng)殖中，胃腸道疾病是十分常見的一類疾病。引起肉牛胃腸道疾病的致病原因多樣，包括致病原因和飼養(yǎng)管理。在具體診治中，應(yīng)該結(jié)合肉牛臨床癥狀，明確致病原因，采取針對性措施進(jìn)行治療，避免隨意用藥，增加藥物耐藥性，降低治療效果。

(13)

其中，ai(i=1,2,3,4,5)表示未知系數(shù)。不同物面深度下未知系數(shù)ai的取值情況見表3。

圖10給出了物面深度分別等于50 m和250 m時，射流寬度數(shù)學(xué)模型與邊界元數(shù)值結(jié)果的整體對比。從兩者相對位置來看，邊界元數(shù)值結(jié)果的散點(diǎn)幾乎都座落在數(shù)學(xué)模型曲面上，最大誤差只有1.4%?？梢姡淞鲗挾葦?shù)學(xué)模型結(jié)合相應(yīng)的物面深度系數(shù)能較精確地表示射流寬度。

通過圖10中數(shù)學(xué)模型曲面上顏色的分布情況，能很直觀地反映出射流寬度隨著藥包質(zhì)量和無量綱距離參數(shù)的變化趨勢。射流寬度沿著無量綱距離參數(shù)這個方向的變化不是很明顯，較大的射流寬分布在藥包質(zhì)量較大的紅色區(qū)域。

表3不同物面深度下的射流寬度系數(shù)

圖10 射流寬度數(shù)學(xué)模型與邊界元數(shù)值結(jié)果的整體對比

2.2.3 射流高度模型

射流高度與無量綱距離參數(shù)的關(guān)系也采用多項(xiàng)式進(jìn)行描述。通過多個函數(shù)對比和篩選，選擇以下數(shù)學(xué)模型來描述射流高度L與藥包質(zhì)量和無量綱距離參數(shù)之間的關(guān)系：

(14)

其中，ai(i=1,2,3,4,5)表示未知系數(shù)。不同物面深度下未知系數(shù)ai的取值情況見表4。

表4不同物面深度下的射流高度系數(shù)

圖11給出了物面深度分別等于50 m和250 m時，射流高度數(shù)學(xué)模型與邊界元數(shù)值結(jié)果沿著藥包質(zhì)量和無量綱距離參數(shù)這兩個方向的整體對比情況。從圖11可以看出，射流高度數(shù)學(xué)模型的曲面結(jié)果和邊界元數(shù)值結(jié)果吻合得較為一致，通過分析，兩者之間最大誤差不超過3.1%。通過射流高度數(shù)學(xué)模型曲面顏色的分布情況，能很直觀的看出射流高度隨著藥包質(zhì)量和無量綱距離參數(shù)的變化趨勢。隨著藥包質(zhì)量的增加，每一藥包質(zhì)量下對應(yīng)的最大射流高度逐漸朝著無量綱距離參數(shù)增大方向偏移。

圖11 時射流高度數(shù)學(xué)模型與邊界元數(shù)值結(jié)果的整體對比

2.2.4 射流載荷模型

由于氣泡射流載荷的復(fù)雜性，目前關(guān)于氣泡射流載荷的表現(xiàn)形式并沒有明確的定論。但是，國內(nèi)外關(guān)于水射流的砰擊特性卻做了較為系統(tǒng)的研究，并取得了較豐碩的成果[14-16]，因此這一節(jié)將在射流速度數(shù)學(xué)模型、射流寬度數(shù)學(xué)模型、射流高度數(shù)學(xué)模型和水射流理論的基礎(chǔ)上，研究描述射流載荷變化趨勢的數(shù)學(xué)模型，旨在為近物面水下爆炸氣泡射流載荷工程預(yù)報提供一定的指導(dǎo)。

當(dāng)水射流以較高的速度砰擊到結(jié)構(gòu)表面時，由于速度較快，流體會被瞬間壓縮產(chǎn)生沖擊波，如圖12所示，因此高速水射流砰擊平板那一瞬時的壓力通常認(rèn)為是水錘壓力。在一維假設(shè)下，基于可壓縮流體力學(xué)的質(zhì)量和動量定理，水射流作用在剛性板表面時，在砰擊區(qū)域中心點(diǎn)處的初始壓力Pm可以表示如下：

Pm=ρ1cV

(15)

其中：ρ1表示流體密度；c表示沖擊波傳播速度；V表示水射流速度。

圖12 水射流砰擊過程

水錘壓力持續(xù)的時間很短，其持續(xù)的時間為稀疏波從水射流邊界傳播到中心的時間(圖12)，因此水錘壓力作用的時間τ可以表示為：

(16)

其中，D為水射流寬度。

水錘壓力釋放之后，砰擊區(qū)域中心點(diǎn)處的壓力逐漸趨于穩(wěn)定，即進(jìn)入駐壓階段(圖12)，該階段的水動壓力大小為：

(17)

當(dāng)考慮材料變形對水射流砰擊的影響時，需要對公式(15)中的壓力進(jìn)行修正，修正后砰擊區(qū)域中心點(diǎn)處的初始壓力可以表示為：

(18)

其中：ρ1和ρ2表示流體和材料的密度；c1和c2表示沖擊波在流體和材料中的傳播速度。

圖13給出較典型的水射流以570 m/s沖擊平板時中心點(diǎn)處的壓力曲線。壓力曲線的峰值壓力為850 MPa，數(shù)值與實(shí)驗(yàn)水錘壓力幾乎相等(未超過10%)。水錘壓力持續(xù)的時間約為1 μs與式(16)計(jì)算的較為一致，這里水射流的寬度為3 mm。隨后，壓力進(jìn)入駐壓階段，該階段的壓力約為式(17)計(jì)算的壓力值。接下來，在圖13水射流實(shí)驗(yàn)載荷的基礎(chǔ)上，建立水射流載荷數(shù)學(xué)模型。通過對實(shí)驗(yàn)數(shù)據(jù)分析可知，稀疏波從水射流邊緣傳播到中心這個階段，水錘壓力主要以指數(shù)形式進(jìn)行衰減，而在穩(wěn)態(tài)階段，壓力趨于平穩(wěn)，因此建立以下數(shù)學(xué)模型來描述水射流載荷的變化趨勢：

P=Aexp(-Bt)+C

(19)

其中：A,B和C為未知系數(shù)；通過初始時刻壓力峰值Pm，τ時刻的壓力Pτ和無窮遠(yuǎn)時刻的壓力Pe可以對未知系數(shù)進(jìn)行求解，然后再將系數(shù)代入就可以得到水射流載荷隨時間變化的表達(dá)式：

(20)

圖13 水射流載荷數(shù)學(xué)模型與實(shí)驗(yàn)數(shù)據(jù)對比

3 結(jié)束語

針對水下爆炸氣泡射流速度、射流高度、射流寬度和射流載荷缺乏系統(tǒng)性研究的這一問題，本文基于邊界積分方程研究了水下爆炸氣泡射流特性與藥包質(zhì)量、無量綱距離參數(shù)和物面深度之間的關(guān)系。通過分析不同物面深度下射流特征參數(shù)分別與藥包質(zhì)量和無量綱距離參數(shù)的變化趨勢，建立了射流速度、射流寬度和射流高度3種數(shù)學(xué)模型，并驗(yàn)證了它們的有效性。爆炸氣泡射流速度隨著藥包質(zhì)量的增加而減小，變化趨勢符合冪函數(shù)Wa(0