李海東，齊江輝，吳述慶

(武漢第二船舶設計研究所， 武漢 430064)

1 引言

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

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

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

本文基于VOF多相流模型，采用Rayleigh-Plesset方程的Singhal全空化模型模擬水下航行體超空泡流動。采用重疊網格方法實現航行體運動及橫流擾動，模擬不同航行速度及橫流擾動大小下空泡特性及航行體水動力特性。對比航行過程中橫流擾動及擾動消失后的流場特性，為復雜工況的超空泡航行體設計提供了參考。

2 數值仿真模型建立

2.1 數值方法

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

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

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

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

2.2 數值模型及網格劃分

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

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

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

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

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

圖3 網格劃分結果示意圖Fig.3 Results of meshing

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

圖4 網格無關性驗證曲線Fig.4 Grid independence verification

2.3 數值方法驗證

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

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

圖6 數值仿真空泡輪廓與實驗結果曲線Fig.6 Comparison of numerical simulation cavitation contour and experimental results

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

3 結果及分析

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

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

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

3.2 橫流擾動下的空泡變化分析

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

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

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

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

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

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

3.3 橫流速度對超空泡的影響

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

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

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

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

4 結論

1) 建立了水下航行體超空泡數值計算模型，驗證了計算模型的精度及可靠性;進行了網格不確定度分析，驗證了數值計算網格無關性。

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

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