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

(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]等在實驗中觀測到，并被認(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)行了實驗觀測和研究；賀銘等[7]利用歐拉有限元數(shù)值模型對近場水下爆炸氣泡與雙層破口結(jié)構(gòu)的相互作用機(jī)理進(jìn)行了研究，主要分析了艙室涌流及流場演化等物理屬性；李梅[8]等研究了雙爆源氣泡與水面的相互作用過程；朱錫、張阿漫等[9-13]都開展過小當(dāng)量真實裝藥水下爆炸產(chǎn)生氣泡的實驗，研究氣泡脈動周期和最大半徑等參數(shù)，得到了氣泡全周期的運(yùn)動形態(tài)。但是，由于很多方面的局限性，目前對于近物面，尤其一些極近物面不同參數(shù)下的氣泡坍塌運(yùn)動和射流載荷特性的研究并不系統(tǒng)。

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

1 計算模型

1.1控制方程

假設(shè)流體絕熱，氣泡內(nèi)部和外界流場不存在質(zhì)量和熱交換。氣泡在外界流場的速度在大部分時間內(nèi)都要遠(yuǎn)小于聲速，因此外界流場的可壓縮性忽略不計。氣泡運(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)中的粘性項不考慮時，方程變?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)部氣體的粘性極小，通常忽略不計，因此切向壓力邊界條件可以寫成：

μ·τ·δ·n=0

(7)

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

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

(8)

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

(9)

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

(10)

(11)

2 計算結(jié)果驗證與分析

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

實驗驗證的工況見表1。這里采用相同的工況對邊界積分方程進(jìn)行理論驗證，理論驗證主要將邊界積分方程計算的結(jié)果和Rayleigh-Plesset方程的理論值進(jìn)行對比，包括半徑曲線和氣泡表面節(jié)點速度。在這之前，先對邊界積分方程進(jìn)行收斂性驗證，確定氣泡表面節(jié)點數(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時，射流速度基本收斂。因此在隨后的計算中，將速度勢變化值Δφ取為0.01。圖3中給出了氣泡表面不同節(jié)點數(shù)量的射流速度收斂曲線，當(dāng)節(jié)點速度大于60時基本收斂，隨后的計算將節(jié)點數(shù)量N取為60。圖4和圖5分別給出邊界積分方程數(shù)值解和Rayleigh-Plesset方程理論解兩者關(guān)于氣泡表面節(jié)點速度和氣泡半徑的對比曲線。由圖4和圖5可見，無論是節(jié)點速度曲線還是半徑曲線，邊界積分方程的數(shù)值結(jié)果和理論解都吻合得較好。

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

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

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

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

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

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

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

2.2氣泡射流載荷計算

2.2.1 射流速度模型

氣泡射流與無量綱距離參數(shù)的關(guān)系可以采用多項式來表示，通過比較一系列不同的數(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é)模型周圍，通過計算，兩者之間的最大誤差不超過2%。

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

為了進(jì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é)模型計算的結(jié)果吻合較為一致。即使射流速度在無量綱距離參數(shù)γ=0.7存在拐點，射流速度數(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)系采用多項式進(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é)果的散點幾乎都座落在數(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)系也采用多項式進(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ū)域中心點處的初始壓力Pm可以表示如下：

Pm=ρ1cV

(15)

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

圖12 水射流砰擊過程

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

(16)

其中，D為水射流寬度。

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

(17)

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

(18)

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

圖13給出較典型的水射流以570 m/s沖擊平板時中心點處的壓力曲線。壓力曲線的峰值壓力為850 MPa，數(shù)值與實驗水錘壓力幾乎相等(未超過10%)。水錘壓力持續(xù)的時間約為1 μs與式(16)計算的較為一致，這里水射流的寬度為3 mm。隨后，壓力進(jìn)入駐壓階段，該階段的壓力約為式(17)計算的壓力值。接下來，在圖13水射流實驗載荷的基礎(chǔ)上，建立水射流載荷數(shù)學(xué)模型。通過對實驗數(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ù)據(jù)對比

3 結(jié)束語

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