橫流擾動(dòng)下的魚雷超空泡及水動(dòng)力特性研究

李海東，齊江輝，吳述慶

(武漢第二船舶設(shè)計(jì)研究所， 武漢 430064)

1 引言

在高速航行水下航行體周圍由于壓力降低至汽化壓力以下，導(dǎo)致航行體周圍產(chǎn)生爆發(fā)式的汽化，汽化形成的小氣泡逐漸融合形成大的氣泡包裹住航行體，這稱作超空化現(xiàn)象。超空泡魚雷就是利用這一原理，由于空氣的密度遠(yuǎn)遠(yuǎn)小于水的密度，因此被空泡包裹住的航行體航行阻力相比于水中會(huì)大大降低，極大提升了航行體的航行速度。超空泡航行是水下航行體高速航行的主要狀態(tài)，因此受到了廣泛的研究和關(guān)注。

國(guó)內(nèi)外對(duì)于超空泡現(xiàn)象開展了很多研究工作，相關(guān)研究方法較為成熟。在試驗(yàn)方面，主要是開展水洞試驗(yàn)測(cè)量航行體的空泡特性及水動(dòng)力性能。張珂等[1]在多功能高速空泡水洞中開展超空泡航行器尾部滑行流場(chǎng)特性試驗(yàn)，并開展了數(shù)值仿真計(jì)算，驗(yàn)證了數(shù)值仿真方法的有效性。栗夫園等[2-3]通過試驗(yàn)研究了錐形空化器的空化特性。時(shí)素果等[4]開展超空泡航行體自由航行試驗(yàn)，觀察到航行體航行過程中出現(xiàn)的周期性尾拍現(xiàn)象。Savchenko[5-6]開展大量的航行體約束模型試驗(yàn)，得到了很多空泡及水動(dòng)力變化規(guī)律。在數(shù)值模擬方面，陳偉善等[7]研究了空化器形狀對(duì)超空泡射彈尾拍運(yùn)動(dòng)的影響，具體分析了3種頭型射彈的尾拍特性。孟慶昌[8]對(duì)不同彈頭形狀的高速射彈垂直入水過程的空泡形態(tài)及水動(dòng)力特性進(jìn)行了研究。劉富強(qiáng)等[9]對(duì)并列超空泡射彈的彈道特性進(jìn)行了數(shù)值分析，研究了并列間距對(duì)彈體運(yùn)動(dòng)穩(wěn)定性的影響。

在公開發(fā)表的研究文獻(xiàn)中，對(duì)于橫流擾動(dòng)下的超空泡流動(dòng)研究較少。范春永等[10]對(duì)水下航行體在側(cè)方來流作用下的超空泡現(xiàn)象進(jìn)行了研究，研究發(fā)現(xiàn)相同側(cè)方流動(dòng)速度下，對(duì)流速度越大，超空泡受沖擊越弱，航行體阻力變化較小。其主要研究了不同側(cè)方來流擾動(dòng)下的空泡特性，并未對(duì)航行體的水動(dòng)力特性進(jìn)行分析。橫流擾動(dòng)下的超空泡水動(dòng)力特性變化較大，涉及空泡的多階段特性，因此對(duì)橫流擾動(dòng)下的超空泡問題研究很有必要。

本文基于VOF多相流模型，采用Rayleigh-Plesset方程的Singhal全空化模型模擬水下航行體超空泡流動(dòng)。采用重疊網(wǎng)格方法實(shí)現(xiàn)航行體運(yùn)動(dòng)及橫流擾動(dòng)，模擬不同航行速度及橫流擾動(dòng)大小下空泡特性及航行體水動(dòng)力特性。對(duì)比航行過程中橫流擾動(dòng)及擾動(dòng)消失后的流場(chǎng)特性，為復(fù)雜工況的超空泡航行體設(shè)計(jì)提供了參考。

2 數(shù)值仿真模型建立

2.1 數(shù)值方法

超空泡航行涉及水-汽-液三相流動(dòng)，為求解得到超空泡形態(tài)及航行體水動(dòng)力性能，采用VOF多相流模型可以精準(zhǔn)捕捉到三相交界面的變化。超空泡流動(dòng)中選用Shih[10]提出的標(biāo)準(zhǔn)k-ε兩方程模型，近壁面采用標(biāo)準(zhǔn)壁面函數(shù)處理方法，可以準(zhǔn)確地模擬粘性流動(dòng)中的各項(xiàng)不均勻性，同時(shí)具有較高的數(shù)值精度。

空化模型選取基于Rayleigh-Plesset方程的Schnerr-Sauer空化模型，該模型是Rayleigh-Plesset方程的簡(jiǎn)化形式，忽略了空泡生成過程中加速及表面張力的影響，對(duì)于大多數(shù)的工程應(yīng)用而言，這種簡(jiǎn)化的影響是可以忽略的[11]。

為方便地實(shí)現(xiàn)對(duì)穩(wěn)定航行體突加橫流待空泡穩(wěn)定后橫流消失這一過程的模擬，采用重疊網(wǎng)格方法實(shí)現(xiàn)航行體的運(yùn)動(dòng)及橫流擾動(dòng)，將計(jì)算域劃分為重疊網(wǎng)格區(qū)域和背景網(wǎng)格區(qū)域。重疊網(wǎng)格區(qū)域與背景網(wǎng)格區(qū)域之間通過重疊網(wǎng)格邊界條件實(shí)現(xiàn)數(shù)值傳輸，重疊網(wǎng)格的應(yīng)用文獻(xiàn)[12-14]中已經(jīng)有詳細(xì)的介紹，此處不再累述。

流體屬性為不可壓縮流體，黃闖等[15-16]研究顯示當(dāng)航行體航速小于900 m/s時(shí)，流體的可壓縮性對(duì)航行體附近的空泡影響很小，因此本文采用不可壓縮流體，飽和蒸氣壓為3 540 Pa。

2.2 數(shù)值模型及網(wǎng)格劃分

本文參考“風(fēng)雪”魚雷，對(duì)魚雷外形進(jìn)行簡(jiǎn)化，結(jié)合Serebryakov給出的無尾翼魚雷模型，設(shè)計(jì)出一種無尾翼超空泡魚雷。無尾翼超空泡魚雷由空化器、圓錐段、平行段及尾噴管組成，其外形及尺寸參數(shù)如圖1所示，計(jì)算中模型采用1∶1比例進(jìn)行數(shù)值分析。

圖1 超空泡魚雷示意圖Fig.1 Supercavitating torpedo shape and size parameters

為避免邊界可能對(duì)超空泡魚雷的流場(chǎng)產(chǎn)生影響，計(jì)算域要設(shè)置足夠大，由于魚雷需要考慮橫向運(yùn)動(dòng)，不能采用對(duì)稱域計(jì)算，因此本文計(jì)算域選為圓柱體，圓柱體直徑為魚雷最大直徑的50倍，圓柱體長(zhǎng)度為30倍魚雷長(zhǎng)度。計(jì)算域分為重疊網(wǎng)格區(qū)域和背景網(wǎng)格區(qū)域，計(jì)算域尺寸及邊界條件設(shè)置如圖2所示。計(jì)算中通過給定重疊網(wǎng)格區(qū)域的運(yùn)動(dòng)速度V0實(shí)現(xiàn)魚雷的水下航行模擬。

圖2 計(jì)算域尺寸及邊界條件示意圖Fig.2 Computational domain size and boundary conditions

網(wǎng)格劃分方法采用切割體網(wǎng)格，重疊網(wǎng)格區(qū)域設(shè)置網(wǎng)格基本尺寸為20 mm，目標(biāo)表面尺寸及最小表面尺寸均為2 mm，魚雷表面網(wǎng)格尺寸為1 mm，魚雷表面附近劃分邊界層，邊界層厚度為1 mm，劃分6層邊界層，邊界層增長(zhǎng)率為1.05。背景網(wǎng)格基本尺寸為0.5 m，最小網(wǎng)格尺寸為0.05 m。采用重疊網(wǎng)格方法計(jì)算時(shí)，計(jì)算結(jié)果的精度很大程度取決于重疊網(wǎng)格與背景網(wǎng)格之間數(shù)據(jù)傳遞的精度，因此需要對(duì)背景網(wǎng)格進(jìn)行網(wǎng)格加密。加密區(qū)域?yàn)楸尘熬W(wǎng)格與重疊網(wǎng)格重疊的區(qū)域，加密區(qū)域一般稍大于重疊網(wǎng)格區(qū)域。加密區(qū)域采用各向同性網(wǎng)格，即網(wǎng)格3個(gè)方向尺寸相同均設(shè)置為2 mm。最終網(wǎng)格劃分結(jié)果如圖3所示，網(wǎng)格數(shù)量約為337萬。

圖3 網(wǎng)格劃分結(jié)果示意圖Fig.3 Results of meshing

在采用CFD方法進(jìn)行數(shù)值模擬時(shí)，需要進(jìn)行網(wǎng)格無關(guān)性分析[17]，采用同樣的網(wǎng)格劃分比例，僅改變基本網(wǎng)格尺寸，建立3種不同尺度網(wǎng)格數(shù)量分別為206萬、337萬、506萬，在速度為100 m/s無橫流擾動(dòng)工況下進(jìn)行數(shù)值模擬，所得到空泡對(duì)比如圖4所示。從圖4中可以看出，337萬網(wǎng)格和506萬網(wǎng)格計(jì)算結(jié)果基本一致，空泡輪廓最大相對(duì)偏差為0.89%，而206萬網(wǎng)格則差別稍大，與337萬網(wǎng)格計(jì)算空泡輪廓最大相對(duì)偏差為3.1%。為保證較高計(jì)算精度的同時(shí)提高計(jì)算效率，本文選取337萬網(wǎng)格劃分方案。

圖4 網(wǎng)格無關(guān)性驗(yàn)證曲線Fig.4 Grid independence verification

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

為驗(yàn)證本文數(shù)值方法的準(zhǔn)確性，采用Hrubes[15]實(shí)驗(yàn)中相同的工況進(jìn)行數(shù)值模擬，數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果，如圖5、圖6所示。從圖6中可以看出，本文數(shù)值模擬的空泡形態(tài)、大小與Hrubes實(shí)驗(yàn)基本一致，空泡輪廓與試驗(yàn)值最大相對(duì)誤差約為4.3%，得到的空泡形態(tài)曲線也與試驗(yàn)結(jié)果吻合良好，證明本文所建數(shù)值模型在進(jìn)行超空泡模擬時(shí)有良好的精度。

圖5 數(shù)值仿真空泡形態(tài)與實(shí)驗(yàn)結(jié)果圖Fig.5 Comparison of numerical simulation cavitation shape and experimental results

圖6 數(shù)值仿真空泡輪廓與實(shí)驗(yàn)結(jié)果曲線Fig.6 Comparison of numerical simulation cavitation contour and experimental results

本文研究橫流擾動(dòng)下的超空泡魚雷空泡及水動(dòng)力性能，通過給定重疊網(wǎng)格域的速度模擬魚雷的水下航行，待魚雷空泡形態(tài)及水動(dòng)力特性穩(wěn)定后，給重疊網(wǎng)格區(qū)域施加橫向運(yùn)動(dòng)模擬突加的橫流擾動(dòng)，待魚雷空泡形態(tài)及水動(dòng)力特性再次穩(wěn)定后，撤去重疊網(wǎng)格的橫向運(yùn)動(dòng)分量，繼續(xù)觀察魚雷的空泡及水動(dòng)力性能。

3 結(jié)果及分析

3.1 空泡形態(tài)分析

魚雷速度為V=100 m/s時(shí)超空泡流動(dòng)如圖7所示。圖7中坐標(biāo)原點(diǎn)位置為魚雷艉部中心點(diǎn)，從圖7可以看出，該魚雷超空泡尾部形態(tài)不規(guī)則且有一定程度的向內(nèi)凹陷，這是由于魚雷模型艉部為斷截面圓柱，使得流動(dòng)在艉部產(chǎn)生分離，因此會(huì)產(chǎn)生超空泡在尾部向內(nèi)凹陷。同時(shí)可以看到，尾部超空泡非完全對(duì)稱，這是由于計(jì)算中考慮了重力的影響，因此超空泡會(huì)產(chǎn)生一定程度的彎曲，尾部空泡長(zhǎng)度有所不同。

圖7 空泡形態(tài)仿真圖Fig.7 Calculation results of cavitation morphology

3.2 橫流擾動(dòng)下的空泡變化分析

在魚雷速度為100 m/s時(shí)且超空泡形態(tài)穩(wěn)定后，給魚雷橫向方向突加流動(dòng)，速度大小為20 m/s，方向?yàn)?Y方向，觀察超空泡形態(tài)在達(dá)到穩(wěn)定前的形狀，如圖8所示。Δt為以施加橫向流動(dòng)為起點(diǎn)的時(shí)間變化，從圖8(a)可以看出，在施加橫向流動(dòng)后很短的時(shí)間內(nèi)(Δt=0.001 s)，超空泡形態(tài)已經(jīng)發(fā)生了明顯的變化，超空泡尾部發(fā)生了嚴(yán)重的流動(dòng)分離現(xiàn)象，圓錐段及平行段超空泡發(fā)生嚴(yán)重的不對(duì)稱性。同時(shí)魚雷側(cè)面有較大的沾濕面積，超空泡長(zhǎng)度暫未發(fā)生明顯變化。由圖8(b)可以看出，圓錐段及平行段超空泡-Y方向半徑明顯增大，這是由于橫向流動(dòng)擾動(dòng)相當(dāng)于魚雷產(chǎn)生橫向速度，魚雷橫向速度與水平速度疊加，使得超空泡沿疊加速度方向延伸。魚雷表面沾濕面積分為圓錐段部分及平行段部分，兩部分沾濕面積被魚雷肩部分隔開來，同時(shí)超空泡長(zhǎng)度有明顯的減小。圖8(c)為超空泡穩(wěn)定過渡段，此時(shí)圓錐段及平行段超空泡已基本穩(wěn)定，尾部超空泡還存在一定程度的流動(dòng)分離。圖8(d)為超空泡穩(wěn)定段，此時(shí)超空泡形態(tài)已基本穩(wěn)定，超空泡不再沿魚雷中線對(duì)稱，-Y方向超空泡半徑相對(duì)無橫流擾動(dòng)時(shí)有明顯增大，尾部超空泡沿著-Y方向有較大程度的延伸。從圖8可以看出，在增加橫向流擾動(dòng)后，魚雷超空泡長(zhǎng)度先減小后逐漸增加至穩(wěn)定值，平行段超空泡半徑在-Y方向明顯增加，但在魚雷Y方向表面有較大的沾濕面積。

圖8 橫流擾動(dòng)下超空泡形狀仿真圖Fig.8 Changes of supercavitation morphology under cross-flow disturbance

在疊加橫流擾動(dòng)下魚雷超空泡穩(wěn)定后，分析橫流突然消失后魚雷超空泡形狀如圖9所示。從圖9(a)可以看出，在橫流突然消失后很短時(shí)間內(nèi)超空泡形態(tài)發(fā)生了明顯的變化，超空泡形態(tài)極為不規(guī)整。在隨后的時(shí)間，超空泡逐漸完全包裹魚雷彈體，同時(shí)尾部的超空泡長(zhǎng)度逐漸變小，超空泡整體呈現(xiàn)軸對(duì)稱形式，如圖9(b)所示。超空泡尾部凹陷位置基本不變，同時(shí)尾部超空泡繼續(xù)延伸直至形成穩(wěn)定的超空泡，此時(shí)超空泡形態(tài)及尺寸與未加橫流擾動(dòng)時(shí)的超空泡穩(wěn)定形態(tài)基本一致。

圖9 橫流擾動(dòng)消失后超空泡形狀仿真圖Fig.9 Changes in supercavitation shape after the cross-flow disturbance disappears

表1給出了無橫流擾動(dòng)階段、橫流擾動(dòng)階段、橫流擾動(dòng)消失階段的魚雷超空泡形態(tài)及阻力變化，由表1可以看出：橫流擾動(dòng)下超空泡長(zhǎng)度稍有增大，超空泡+Y方向尺寸減小，-Y方向尺寸約增大為原來的1.5倍。在橫流擾動(dòng)消失后，超空泡尺寸逐漸恢復(fù)，但超空泡長(zhǎng)度以及Y方向尺寸仍比無橫流穩(wěn)定超空泡稍大，這是由于空泡發(fā)展過程中的時(shí)滯性所致。橫流擾動(dòng)使得魚雷航行阻力約增大為原來的1.3倍，待橫流擾動(dòng)消失后，阻力恢復(fù)為無橫流穩(wěn)定值。

表1 橫流擾動(dòng)下的超空泡形態(tài)及阻力

3.3 橫流速度對(duì)超空泡的影響

魚雷速度為100 m/s，施加橫流速度(Vt)分別為魚雷速度的1%、10%、15%、20%和30%，計(jì)算不同橫流速度下魚雷超空泡形態(tài)變化及水動(dòng)力特性。圖10(a)為未加橫流擾動(dòng)時(shí)的穩(wěn)態(tài)超空泡形態(tài)，圖10(b)～圖10(f)為不同橫流速度時(shí)的穩(wěn)態(tài)超空泡。從圖10中可以看出，當(dāng)橫流速度很小時(shí)(Vt/V=1%)，魚雷超空泡形態(tài)基本沒有變化。魚雷超空泡形態(tài)變化隨著橫流擾動(dòng)速度的增大而變大。Vt/V=1%時(shí)，魚雷體表面已經(jīng)開始出現(xiàn)沾濕面積，隨著橫流速度的增大沾濕面積也隨之增大。

圖10 不同橫流擾動(dòng)速度下超空泡形態(tài)圖Fig.10 Supercavitation morphology under different

圖11為不同橫流擾動(dòng)速度時(shí)的魚雷航行阻力曲線，無橫流擾動(dòng)時(shí)魚雷航行阻力約為365.5 N。橫流速度較小(Vt/V=1%)時(shí)，魚雷航行阻力變化較小，約為383.4 N。隨著橫流擾動(dòng)速度的增大，魚雷航行阻力增量變大，橫流速度Vt/V=20%時(shí)，魚雷航行阻力約為718.8 N，為無橫流擾動(dòng)時(shí)的2倍。同時(shí)可以看出，隨著橫流擾動(dòng)速度的增大，魚雷航行阻力曲線震蕩特性變得劇烈，這是由于橫流使得魚雷周圍流場(chǎng)湍動(dòng)加劇，流場(chǎng)變得更加不穩(wěn)定。

圖11 不同橫流擾動(dòng)速度下魚雷航行體阻力曲線Fig.11 Resistance of torpedo body under different cross-flow disturbance speeds

4 結(jié)論

1) 建立了水下航行體超空泡數(shù)值計(jì)算模型，驗(yàn)證了計(jì)算模型的精度及可靠性;進(jìn)行了網(wǎng)格不確定度分析，驗(yàn)證了數(shù)值計(jì)算網(wǎng)格無關(guān)性。

2) 基于重疊網(wǎng)格方法，模擬橫流擾動(dòng)下魚雷超空泡形態(tài)及水動(dòng)力特性，橫流擾動(dòng)使得魚雷超空泡不再具有對(duì)稱性，尾部流場(chǎng)流動(dòng)分離加劇，超空泡不再完全包裹魚雷表面，魚雷表面產(chǎn)生較多沾濕面積，航行阻力增加。

3) 在魚雷航行速度一定時(shí)，橫流擾動(dòng)速度越大，魚雷超空泡變形越嚴(yán)重甚至發(fā)生脫落，魚雷表面沾濕面積變大，航行阻力增量越大；在橫流擾動(dòng)速度達(dá)到魚雷航行速度的20%時(shí)，航行阻力大幅增大，約為無橫流擾動(dòng)時(shí)的2倍。