魯建華,周東輝,賈會霞

(1.浙江工業(yè)職業(yè)技術(shù)學(xué)院, 浙江 紹興 312000;2 浙江理工大學(xué) 機(jī)械與自動控制學(xué)院, 杭州 310018)

0 引言

當(dāng)水下航行體高速運(yùn)動時,其表面所承受的壓力會降低到水的飽和蒸汽壓,從而導(dǎo)致航行體表面附近的水汽化形成超空泡,超空泡可使航行體受到的阻力最大減少90%左右[1]。由于超空泡顯著的減阻作用,超空泡流動在高速魚雷、射彈、反潛火箭等水中兵器設(shè)計中有著重要的應(yīng)用[2]。在實(shí)際應(yīng)用中,由于海洋水域環(huán)境和發(fā)射方式的不同,還會涉及超空泡航行體在淺水區(qū)運(yùn)動,此種情況存在超空泡演化和自由面波動的耦合效應(yīng),與深水位(無限水域)下的超空泡有顯著差別。近些年,相關(guān)的研究成果在國防軍事工程和海洋運(yùn)輸工業(yè)領(lǐng)域有著重要的應(yīng)用。

國內(nèi)外學(xué)者對近自由面的空化問題開展了一些研究,其中對自由面與空化水翼的相互影響進(jìn)行了較多研究,目的是開發(fā)出能高速航行的艦船。Dawson等[3]進(jìn)行了近自由面楔型水翼的空化流動研究,獲得了空泡的發(fā)展過程。Faltinsen等[4]通過非線性分析研究了水翼在重力作用下自由液面附近的空化流動,研究表明自由面會影響水翼的升阻系數(shù)。陳鑫等[5]開展了淺水中二維對稱楔的空化問題數(shù)值模擬研究,研究表明空化數(shù)越小,自由液面的變形范圍越大。陳慶任等[6]采用非定常面元法數(shù)值分析了近自由面三維振動水翼的水動力特性,研究了不同浸深、航速和振動頻率的影響。Zhou等[7]基于OpenFOAM平臺開出了用于求解近自由面空化問題的多相流求解器,可以較好的捕捉自由面的變化。王一偉等[8]通過實(shí)驗(yàn)和數(shù)值模擬研究了近自由面射彈的云狀空化問題,獲得了空泡的演化特性,分析了自由面對再進(jìn)入射流和空泡脫落的影響。徐暢等[9]利用實(shí)驗(yàn)和數(shù)值模擬方法研究了近自由面條件下高速射彈誘導(dǎo)的通氣超空泡流動,分析了超空泡的形成機(jī)制和減阻效果。施紅輝等[10]利用輕氣泡發(fā)射技術(shù)開展了近自由面超空泡射彈實(shí)驗(yàn),獲得了超空泡的演化特性和射彈的運(yùn)動特性,分析了自由面/超空泡的相互作用。張亞濤等[11]利用VOF數(shù)值模擬方法開展了近自由面射彈的超空泡流動特性研究,分析了水深對超空泡形態(tài)的影響。

目前近自由面的空化問題國內(nèi)外的研究主要針對單個航行體,而在一些情況下需要在淺水區(qū)對水面目標(biāo)發(fā)射超空泡射彈或者超空泡魚雷進(jìn)行飽和攻擊,這必然涉及到超空泡航行體的串列運(yùn)動,例如對于大型水面艦艇,往往需要魚雷群進(jìn)行飽和打擊。因此需要弄清超空泡流場之間的相互作用規(guī)律以及自由面的影響,這些問題關(guān)系到超空泡武器串列運(yùn)動的彈道穩(wěn)定性。因此本文基于求解N-S方程的有限體積法,利用VOF的數(shù)值模擬方法計算了近自由面串列航行體的超空泡流場,獲得了自由液面影響下串列航行體的超空泡演化特性及流場特性。

1 控制方程及數(shù)值方法

1.1 控制方程

本文采用流體體積函數(shù)(VOF)模型來捕獲相間界面,混合相的連續(xù)方程和動量方程分別為

(1)

(2)

其中:xi、xj分別是笛卡爾坐標(biāo)分量;ui、uj分別為笛卡爾坐標(biāo)系中的速度分量;P為流場的壓力;SM為附加的源相;ρm、μm分別為混合相的密度和動力粘度,滿足如下的方程

ρm=αvρv+αgρg+α1ρ1

(3)

μm=αvμv+αgμg+α1μ1

(4)

其中:ρv、ρg、ρ1分別為水蒸氣、空氣、水的密度;μv、μg、μ1分別為水蒸氣、空氣、水的動力粘度;αv、αg、α1分別為水蒸氣、空氣、水的動力粘度,并滿足αv+αg+α1=1。

采用RNGk-ε湍流模型[12]對流體的控制方程進(jìn)行封閉求解,空化是流線強(qiáng)烈彎曲導(dǎo)致的,該模型可以更好地處理高彎曲流線及高應(yīng)變率的流動。流動中的空化問題采用Schnerr-Sauer空化模型[13]求解。水蒸氣體積分?jǐn)?shù)的一般方程為

(5)

其中:Re、Rc分別為水蒸氣的蒸發(fā)速率和冷凝速率。

當(dāng)Pv≥P時

(6)

當(dāng)Pv≤P時

(7)

1.2 數(shù)值方法和邊界條件

航行體采用的物理模型為直徑D=6 mm,長度L=48 mm的柱體,如圖1所示。計算域和邊界條件設(shè)置如圖2所示,整個計算域的尺寸為2 500 mm×960 mm,前發(fā)航行體距離來流入口為480 mm,兩航行體之間的距離為480 mm。定義航行體中心線距離自由面的高度為航行體浸沒深度h。計算域左側(cè)邊界設(shè)定為速度入口,來流速度大小為60 m/s,右側(cè)邊界設(shè)置為壓力出口,上邊界設(shè)置為Symmetry,下邊界設(shè)為壁面條件。計算域的網(wǎng)格采用結(jié)構(gòu)網(wǎng)格進(jìn)行劃分,在航行體周圍和自由面附近進(jìn)行加密處理,如圖3所示。計算方法中,壓力與速度之間的耦合求解采用 Coupled算法,各項(xiàng)體積率離散采用Modified HRIC格式。

圖1 航行體的物理模型

圖2 計算域和邊界條件設(shè)置示意圖

圖3 網(wǎng)格劃分示意圖

1.3 數(shù)值方法驗(yàn)證

首先進(jìn)行了網(wǎng)格無關(guān)性驗(yàn)證,建立3種不同數(shù)量的網(wǎng)格進(jìn)行計算結(jié)果對比,網(wǎng)格數(shù)分別為85.8萬(case1)、119.3萬(case2)、247.1萬(case3)。圖4為不同網(wǎng)格密度下前發(fā)航行體超空泡輪廓的計算結(jié)果。由圖4可知,隨著網(wǎng)格數(shù)量的增加,在case2和case3網(wǎng)格密度條件下前發(fā)航行體的超空泡輪廓已經(jīng)基本一致了。綜合考慮計算效率和計算結(jié)果的準(zhǔn)確性,計算時選擇了case2的網(wǎng)格。

圖4 t=4 ms,不同網(wǎng)格數(shù)量條件下前發(fā)航行體的超空泡輪廓對比

由于目前還沒有近自由面串列航行體實(shí)驗(yàn)的相關(guān)文獻(xiàn),故選擇了文獻(xiàn)[14]中的近自由面單個航行體實(shí)驗(yàn)進(jìn)行數(shù)值模擬方法有效性驗(yàn)證。實(shí)驗(yàn)中航行體的浸沒深度h為18 mm,速度為38.6 m/s。圖5航行體超空泡形態(tài)的數(shù)值模擬結(jié)果和實(shí)驗(yàn)結(jié)果的對比圖。從圖5可以看出,兩者的超空泡形態(tài)有較好的一致性。由于超空泡上表面的水層很薄,實(shí)驗(yàn)照片不能明顯的辨認(rèn)出空泡上表面的輪廓,因此只給出了超空泡下表面輪廓曲線的定量對比,如圖6所示。由圖6可知,整體上兩者的超空泡下表面輪廓一致性較好,但是在空泡中部截面附近和空泡尾部附近實(shí)驗(yàn)所得的結(jié)果稍大于數(shù)值模擬結(jié)果,但誤差在可接受范圍內(nèi)。誤差的主要原因:一是實(shí)驗(yàn)存在測量偏差,主要為實(shí)驗(yàn)拍攝空泡時存在光線折射;二是由于超空泡尾部閉合區(qū)域是充滿蒸汽、液滴和旋渦的多相流湍流區(qū),還涉及自由面的作用,非常復(fù)雜,數(shù)值模擬很難精確模擬超空泡尾部滯止及崩潰。

圖6 數(shù)值模擬獲得的超空泡下表面輪廓與實(shí)驗(yàn)結(jié)果的對比

2 結(jié)果與討論

2.1 超空泡的演化特性分析

圖7給出了串列航行體在不同浸沒深度下的水相圖。從圖7中可以看出,4種浸沒深度情況下串列航行體的超空泡發(fā)展規(guī)律是相同的,即隨著時間的推移,前發(fā)航行體與后發(fā)航行體的超空泡流場產(chǎn)生干擾,發(fā)生了2個超空泡融合、分離、尾空泡潰滅等流動現(xiàn)象,最后導(dǎo)致后發(fā)航行體進(jìn)入前發(fā)航行體的超空泡內(nèi)部。通過對比還可以發(fā)現(xiàn),串列航行體離自由面越近,自由面的作用效果越強(qiáng),超空泡形態(tài)的上下對稱性越差,并且超空泡發(fā)生融合的時間越早,當(dāng)h=30 mm時,串列航行體的超空泡尾部發(fā)生了自由面上方空氣侵入現(xiàn)象,其他3種浸沒深度的工況沒有發(fā)生這種現(xiàn)象。

(8)

其中:cx0取值為0.827;k為經(jīng)驗(yàn)常數(shù),其取值范圍一般為0.9～1.0;σ為空化數(shù)。

超空泡外形輪廓公式為

(9)

圖8表示的為t=4 ms時,不同浸沒深度的前發(fā)航行體的超空泡前沿輪廓與Logvinovich半經(jīng)驗(yàn)公式的對比。隨著浸沒深度的減少,自由面作用越明顯,超空泡上表面輪廓厚度越大,與Logvinovich半經(jīng)驗(yàn)公式計算出的空泡輪廓相差越大。在h=300 mm時,超空泡上下輪廓基本對稱,與Logvinovich半經(jīng)驗(yàn)公式計算的空泡輪廓一致性較好。

圖8 不同浸沒深度的前發(fā)航行體超空泡前沿輪廓對比

2.2 流場特性分析

圖10為不同浸沒深度下串列航行體頭部中心點(diǎn)壓力變化曲線,P為流體的絕對壓強(qiáng),P0為一個標(biāo)準(zhǔn)大氣壓。從圖中可以看出,浸沒深度較小時,兩發(fā)航行體頭部中心點(diǎn)壓力有一定程度的減小。對比前發(fā)航行體和后發(fā)航行體的頭部壓力曲線,在t=6 ms之后,后發(fā)航行體頭部中心點(diǎn)的壓力明顯小于前發(fā)航行體的;隨著時間的推移,后發(fā)航行體頭部中心點(diǎn)壓力逐漸降到最低,約為0.05個大氣壓,這表明此時后發(fā)航行體受到的壓差阻力幾乎為0,其原因?yàn)楹蟀l(fā)航行體的流場受到前發(fā)航行體超空泡流場的干擾,后發(fā)航行體頭部的流體介質(zhì)逐漸由液相水變成水蒸氣相,水蒸氣的密度為液態(tài)水密度的0.056%,相差了3個數(shù)量級。

圖10 不同浸沒深度條件下串列航行體頭部中心點(diǎn)壓力變化曲線

圖11給出了t=6 ms時,不同浸沒深度下串列航行體的壓力云圖及其徑向壓力分布,H1為點(diǎn)到射彈頭部中心點(diǎn)前方2*D處位置的豎直距離(豎直向上為負(fù),豎直向下為正)。

圖11 t=6 ms時,不同浸沒深度條件下串列航行體的壓力云圖及其徑向壓力分布

由圖11可知,航行體頭部前方為高壓區(qū),其他部分被一個低壓場包裹,低壓場壓力約為飽和蒸汽壓3 540 Pa。隨著浸沒深度的減小,航行體頭部前方的壓力沿徑向分布越來越不對稱,這表明航行體距離自由面越近,航行體上下表面水層對壓力場的約束程度差異越大。

圖12為h=30 mm時,12 ms時刻串列航行體超空泡附近的流線分布。從圖中可知,后發(fā)航行體超空泡尾部水層破碎,空泡內(nèi)部壓強(qiáng)小,外界壓強(qiáng)大,在內(nèi)外壓差作用下自由面上方的空氣被迅速吸入空泡內(nèi)部,其速度最大可達(dá)300 m/s。

圖12 h=30 mm時,12 ms時刻的流線圖

3 結(jié)論

本文中采用數(shù)值模擬方法對近自由面串列航行體的多相流動特性進(jìn)行了研究,獲得的主要結(jié)論如下:

1) 近自由面串列雙航行體的超空泡演化經(jīng)歷了2個空泡的獨(dú)立發(fā)展、空泡彼此的相互融合與分離、尾空泡潰滅等過程,最終導(dǎo)致后發(fā)航行體進(jìn)入了前發(fā)航行體的超空泡內(nèi)部,使后發(fā)航行體受到的阻力減小,表明航行體的串列布置提升了后發(fā)航行體的減阻性能。

2) 航行體距離自由面高度越近,其超空泡形態(tài)上下對稱性越差,與Logvinovich半經(jīng)驗(yàn)公式計算的空泡輪廓偏差越大,并且自由面波動興起的波浪特征高度越大。

3) 浸沒深度較小時,自由面的作用明顯,航行體前端高壓區(qū)壓力分布沿航行體中心線不對稱;浸沒深度較大時,自由液面的作用大大減弱,航行體前端高壓區(qū)壓力分布基本對稱。