水下高速航行體并聯(lián)出水過程的空化特性研究

周東輝，施紅輝，魯建華，薛明瑞

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

0 引 言

并聯(lián)出水問題來源于水下彈道導(dǎo)彈的齊射出水、水下槍炮多管并射形成彈幕攔截或攻擊低空飛行目標(biāo)等海戰(zhàn)背景，具有很強的工程應(yīng)用價值。當(dāng)采用并聯(lián)發(fā)射方式潛射導(dǎo)彈/射彈時，可以顯著地增加武器的突防概率和毀傷效果。并聯(lián)出水過程中不僅涉及了流體介質(zhì)的突變，氣相、固體、液相的三相耦合流動，還存在航行體之間的流場彼此干擾，導(dǎo)致航行體的流場和運動特性非常復(fù)雜。

導(dǎo)彈/射彈以較高速度水下發(fā)射時，其表面發(fā)生空化形成空泡，空泡的演化過程及形態(tài)關(guān)系到水下航行體的受力情況，進而影響到航行體的彈道特性。因此，空泡演化特性的研究一直是該領(lǐng)域的研究重點。在國外，LOGVINOVICH 等[1]提出的空泡截面獨立膨脹原理，可用于計算空泡的形狀。SAVCHENKO 等[2]提出了計算超空泡輪廓的半經(jīng)驗公式。WAUGH 等[3]通過實驗研究了發(fā)射角度和空泡對導(dǎo)彈出水姿態(tài)的擾動影響，給出了出水空泡的形態(tài)。NGUYEN 等[4]利用數(shù)值模擬方法研究了射彈勻速出水過程的超空泡流動，給出了射彈出水過程的阻力變化和超空泡形態(tài)。在國內(nèi)，賈會霞等[5,6]利用高速攝像技術(shù)開展了超空泡射彈出水的實驗研究，獲得了出水過程中空泡的演化過程和射彈的速度變化，給出了不同模型和參數(shù)下超空泡出水形態(tài)尺寸的變化規(guī)律。施紅輝等[7]采用實驗和數(shù)值模擬相結(jié)合的方法研究了超空泡射彈出水的現(xiàn)象，分析了射彈頭部形狀對射彈運動、自由面破碎、超空泡形狀的影響。別慶等[8]基于Mixture 多相流模型，采用動網(wǎng)格技術(shù)和6DOF 技術(shù)建立了潛射導(dǎo)彈水下運動模型，獲得了變速狀態(tài)下導(dǎo)彈的空泡特性。魯傳敬等[9]基于多相流模型，采用動網(wǎng)格技術(shù)數(shù)值模擬了細長體出水過程。魏海鵬等[10]比較了不同多相流模型在航行體出水過程中對流體界面捕獲、壓力計算方面的差異，給出了VOF 模型和Mixture 模型的適用范圍。顏開等[11]對出水空泡的行為進行了很多數(shù)值計算，給出了空泡周圍速度場和壓力場的變化規(guī)律。ZHANG 等[12]數(shù)值模擬了出水空泡潰滅的過程，獲得了出水速度對空泡形態(tài)及潰滅壓力的影響。陳瑛等[13]利用大渦模擬的數(shù)值方法對導(dǎo)彈出水過程中空泡脫離和潰滅進行了研究。

針對并聯(lián)出水問題，MNASRI 等[14]利用VOF 模型和動網(wǎng)格技術(shù)進行了雙圓柱體低速并聯(lián)出水的二維數(shù)值模擬，分析了出水過程中自由面的變形以及流場的相互干擾。盧佳興等[15]進行了回轉(zhuǎn)體齊射出水實驗，分析了回轉(zhuǎn)體間的流體動力干擾和艇速對其影響。畢鳳陽等[16]建立了多細長體水下齊射多相流動與多體運動耦合數(shù)值模擬的計算模型，分析了典型工況下的齊射擾動特性。

綜上所述，現(xiàn)有的研究主要集中在單航行體出水過程的空泡演化、流體動力以及彈道特性等，但對高速航行體并聯(lián)出水的超空泡流動特性研究還偏少。本文基于VOF 多相流模型，引入Schnerr-Sauer 空化模型、SST κ -ω湍流模型、6DOF 剛體運動模型，利用重疊網(wǎng)格技術(shù)建立高速航行體并聯(lián)出水的數(shù)值模型，模擬發(fā)射速度為200 m/s 時，單航行體出水以及不同并聯(lián)間距下雙航行體出水的自由運動過程，獲得了航行體的超空泡演化特性和運動特性，可為并聯(lián)超空泡射彈出水的工程應(yīng)用提供參考。

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

1.1 控制方程

采用VOF 多相流模型描述水、水蒸汽、空氣構(gòu)成的多相流動系統(tǒng)，該模型將三相當(dāng)作密度可變的單一介質(zhì)的混合相來處理，各相共享同一壓力差和速度場?；旌舷嗟倪B續(xù)性方程和動量方程分別為：

式中： ρm為混合相密度， ρm=αvρv+αgρg+α1ρ1， αv，αg， α1分 別為水蒸氣、空氣、水的體積分數(shù)；ui為i軸方向的速度；P為壓力； μm為混合相動力粘度，μm=αvμv+αgμg+α1μ1， μv、 μg、 μ1分別為水蒸汽、空氣、水的動力粘度；uj為j軸方向的速度；gi為重力加速度的分量，F(xiàn)i為體積力的分量。

采用SST κ -ω湍流模型[17]對RANS 方程提供湍流封閉。空化現(xiàn)象采用Schnerr-Sauer 空化模型描述[10]，其控制方程如下：

1.2 邊界條件、網(wǎng)格劃分及數(shù)值方法

數(shù)值模型中采用的航行體模型為截錐型射彈，如圖1 所示。該彈由1 個圓臺和1 個圓柱組成，質(zhì)量為2.98 g，圓柱段直徑D=6 mm，錐段長L0=18 mm，射彈全長L=48 mm，空化器直徑D0=3 mm。

圖1 幾何模型示意圖Fig. 1 Schematic diagram of geometric model

計算域和邊界條件設(shè)置如圖2 所示，計算域為長方體，高780 mm，長210 mm，寬120 mm。計算域的底部為壓力入口，計算域的側(cè)面和上部為壓力出口，采用用戶自定義場函數(shù)定義邊界面上的壓力?？諝庥蚋?00 mm，水域高480 mm，航行體頭部距離自由液面的高為400 mm，航行體A 和航行體B 的表面均設(shè)定為無滑移壁面。上述計算模型中壓力與速度耦合的求解采用C o u p l e d 算法，壓力場和空間離散采用PRESTO!格式，體積率離散采用Modified HRIC。定義兩航行體軸線之間的距離為并聯(lián)間距△d。航行體的質(zhì)心在y方向的位移為縱向位移Sy，在y方向的速度為縱向速度Vy，取豎直向上為正，定義航行體無量綱縱向位移

圖2 計算域及邊界條件設(shè)置Fig. 2 Computational domain and boundary condition setting

在數(shù)值模擬中重疊網(wǎng)格技術(shù)是研究流場中多物體運動最為直接有效的方法之一，因此本文選擇重疊網(wǎng)格技術(shù)實現(xiàn)2 個航行體并聯(lián)出水的數(shù)值模擬。計算區(qū)域分為背景網(wǎng)格區(qū)域和子網(wǎng)格區(qū)域，由于航行體A 和航行體B 分別獨立運動，因此需要劃分2 個子網(wǎng)格區(qū)域。圖3 為計算域的網(wǎng)格劃分示意圖，均采用結(jié)構(gòu)化網(wǎng)格進行劃分，全計算域的網(wǎng)格總數(shù)約為216.5 萬。背景網(wǎng)格區(qū)域中的航行體運動區(qū)域和自由液面附近進行局部加密，以便精確捕捉空泡界面和自由液面的變形，外部背景區(qū)域網(wǎng)格較疏，加快收斂速度。2 個子網(wǎng)格區(qū)域大小相同，均為包裹航行體的圓柱，其長度為10D，直徑為2D，并且采用相同的網(wǎng)格劃分方式，對航行體壁面附近區(qū)域進行網(wǎng)格加密。計算中，航行體的運動通過6DOF 剛體運動模型求解。

圖3 網(wǎng)格劃分示意圖Fig. 3 Schematic diagram of meshing

1.3 數(shù)值方法驗證

根據(jù)文獻[7]中射彈模型2 的實驗數(shù)據(jù)，對水下超空泡射彈出水進行數(shù)值方法驗證。圖4 給出了出水超空泡演化過程的實驗和數(shù)值模擬結(jié)果，相鄰圖片的時間間隔為Δt=0.25 ms。通過對比可以發(fā)現(xiàn)，數(shù)值模擬較好地模擬出了實驗中的出水超空泡演化過程。為了進一步驗證數(shù)值模擬方法有效性，將射彈的無量綱縱向位移的數(shù)值結(jié)果與實驗數(shù)據(jù)進行對比，如圖5 所示?？梢钥闯?，射彈無量綱縱向位移的數(shù)值結(jié)果與實驗數(shù)據(jù)具有較好的一致性，最大誤差約為3.9%。通過上述比較驗證了本文數(shù)值模擬方法的有效性。

圖4 射彈出水超空泡的演化過程對比（Δt=0.25 ms）Fig. 4 Comparisons of water-exit supercavity evolutions of the projectile between numerical simulation and experiment(Δt=0.25 ms)

圖5 關(guān)于射彈無量綱縱向位移的數(shù)值模擬和實驗結(jié)果[7]對比Fig. 5 Comparisons of dimensionless longitudinal displacements of the projectile between numerical simulation and experiment

2 結(jié)果與討論

本文以單航行體出水為參照工況，考慮并聯(lián)間距為2D、3D、4D、5D四種工況，分別對不同工況開展數(shù)值模擬研究，計算工況如表1 所示。表中VA0，VB0，V0分別為航行體A、航行體B、單航行體出水時的初速度，方向為豎直向上。

表1 計算工況表Tab. 1 Numerical simulation cases

2.1 并聯(lián)出水過程的超空泡演化特性

圖6（a）給出了航行體單獨出水過程的超空泡演化過程?？梢钥闯?，航行體在水下運動至出水過程中經(jīng)歷了形成自然超空泡、超空泡隨動、超空泡局部潰滅等。在水下運動階段，由于航行體的運動速度較高，導(dǎo)致航行體肩部位置的壓力降低到水的飽和蒸汽壓以下，此時液相水會不斷汽化形成空泡，隨著空泡不斷的發(fā)展形成超空泡，超空泡的形狀為細長的橢球體，航行體被包裹在超空泡內(nèi)部，只有頭部沾濕，如圖6（a）中t=1 ms 所示。隨著航行體繼續(xù)向上運動，超空泡跟隨航行體向上運動。航行體穿越水面時，被航行體帶入空氣中的空泡在大氣壓力的作用下發(fā)生潰滅，如圖6（a）t=2.8 ms 所示。由于航行體出水速度較高，水面以下的空泡來不及隨航行體出水而被剝離在自由液面以下，如圖6（a）中t=3 ms 所示。當(dāng)航行體完全進入空氣之后，出現(xiàn)了向上運動的水的噴濺，如圖6（a）中t=3.4 ms 所示。6（b）給出了并聯(lián)間距△d= 4D時，航行體并聯(lián)出水過程中超空泡演化過程。為了便于分析，以t=2 ms 的圖像為例，定義2 個航行體相隔較近的一側(cè)為內(nèi)側(cè)，相隔較遠的一側(cè)為外側(cè)?？梢钥闯?，2 個航行體都各自形成了超空泡，航行體被包裹在超空泡內(nèi)部。由于受到相鄰航行體的限制，超空泡內(nèi)側(cè)的擴張受到了抑制，而超空泡的外側(cè)自由擴展，導(dǎo)致雙航行體超空泡外側(cè)輪廓的曲率要大于內(nèi)側(cè)輪廓，形成了一對非對稱超空泡，但雙超空泡形態(tài)在空間上呈現(xiàn)出較好的鏡面對稱特征，如圖6（b）中t= 1~1.8 ms 所示。隨著雙航行體繼續(xù)向上運動，超空泡繼續(xù)發(fā)展，其直徑和長度增加。在穿越水面階段，并聯(lián)出水與單獨出水相似，出現(xiàn)水面以上的空泡潰滅、還未潰滅的空泡被自由液面阻隔等行為，如圖6（b）t=2.6~3 ms 所示。另外還可以看出，超空泡的后半段以并聯(lián)航行體的中軸線為基準(zhǔn)相互靠攏。在雙航行體完全出水后，噴濺的水花是不對稱的，外側(cè)噴濺高度大于內(nèi)側(cè)，如圖6（b）中t= 3.4 ms 所示。

圖6 單獨出水和并聯(lián)出水的超空泡演化Fig. 6 Supercavity evolutions of the single vehicle exiting water and two vehicles exiting water in parallel

圖7為t=2 ms 時，不同并聯(lián)間距下的出水超空泡形態(tài)對比?？梢钥闯?，隨著航行體并聯(lián)間距的減小，超空泡尾部向內(nèi)側(cè)擴張的程度逐漸增加。當(dāng)并聯(lián)間距△d為2D和3D時，雙超空泡尾部壁面融合形成“U”形空泡，而△d為3D和4D的工況沒有出現(xiàn)雙超空泡的尾部融合。其原因主要是由于隨著雙航行體內(nèi)側(cè)流體域變小，該區(qū)域接受雙航行體傳遞動能的水流體減少，引起每個流體質(zhì)點獲得的動能平均值增加，根據(jù)伯努利定理，內(nèi)側(cè)區(qū)域的平均壓力必然降低，增強了超空泡向內(nèi)側(cè)擴張的速度，從而導(dǎo)致雙空泡在內(nèi)側(cè)壁面發(fā)生了融合。

圖7 t=2 ms 時，不同并聯(lián)間距的出水超空泡形態(tài)對比Fig. 7 Comparisons of the shape of supercavitation with different parallel distances when t = 2 ms

2.2 航行體并聯(lián)出水過程中的運動特性

定義航行體的偏航角θ為當(dāng)前航行體軸線與y軸的夾角，逆時針旋轉(zhuǎn)方向為偏航角的正方向。圖8 給出了并聯(lián)間距△d= 4D時，雙航行體偏航角隨時間變化的曲線。從圖中可知，并聯(lián)出水時，航行體A 向外側(cè)進行偏轉(zhuǎn)，航行體B 向另一外側(cè)進行偏轉(zhuǎn)，即2 個航行體的運動向著兩者頭部遠離、尾部靠近的方向偏轉(zhuǎn)，從整體上看航行體A 和航線體B 的偏轉(zhuǎn)運動表現(xiàn)為良好的對稱性。航行體偏轉(zhuǎn)的原因主要是航行體誘導(dǎo)的超空泡沿軸線不對稱，導(dǎo)致航行體受到了不對稱的水動力。在此種工況下，隨著航行體運動時間的增加，偏航角也越來越大，在時間為3.4 ms 時，偏航角達到了3.06○。

圖8 并聯(lián)間距△d = 4D 時，雙航行體的偏航角隨時間變化Fig. 8 Variations of yaw angle of two vehicles with time when parallel distance △d = 4D

鑒于兩航行體運動的對稱性，將基于航行體A 進行并聯(lián)出水過程的運動分析。圖9 給出了不同工況下航行體的縱向位移、縱向速度隨時間的變化?？梢钥闯觯? 種工況下航行體的速度變化趨勢相似，在運動初始階段，航行體的速度衰減很快，這是由于超空泡尚未形成，受到水的阻力較大；隨著航行體產(chǎn)生穩(wěn)定的超空泡，減小了受到的水的阻力，速度衰減趨勢逐漸變緩；2.6 ms 之后，航行體的速度衰減變得很小，原因為2.6 ms 時，航行體開始穿越水面進入空氣中，在穿越過程中，航行體頭部進入空氣中，尾部在空泡中，頭部流體介質(zhì)由液相水變成氣相，密度突變，從而導(dǎo)致航行體受到的壓差阻力變得很小，當(dāng)航行體完全進入空氣中，航行體只受到量級較小的空氣阻力。通過對比縱向速度和無量綱縱向位移的放大圖可以發(fā)現(xiàn)，雖然各工況下航行體縱向速度和無量綱縱向位移相差較小，但是仍有區(qū)別。對于并聯(lián)出水，隨著并聯(lián)間距的減小，航行體的縱向速度衰減略微增快，無量綱縱向位移變小。航行體單獨出水的縱向速度和縱向無量綱位移大于并聯(lián)出水的，這表明并聯(lián)出水對超空泡航行體的減阻產(chǎn)生不利影響。

圖9 不同工況下航行體無量綱縱向位移、縱向速度隨時間的變化Fig. 9 Variations of dimensionless longitudinal displacement and longitudinal velocity of vehicle with time under different cases

圖10 為不同工況下航行體出水過程的偏航角變化。從圖中可知，航行體單獨出水時，航行體幾乎沒有發(fā)生偏轉(zhuǎn)，偏航角的變化近似為 0 度附近的一條直線。對于并聯(lián)出水，并聯(lián)間距△d為2D和3D時，航行體的偏航角隨時間的增加先增加后減小，其原因為航行體向外側(cè)偏轉(zhuǎn)，當(dāng)偏航角達到某一臨界值時，本文工況條件下約為3.4°，航行體的尾部刺穿超空泡內(nèi)測壁面沾濕，產(chǎn)生了橫向沾濕面作用的均布力，均布力產(chǎn)生向內(nèi)側(cè)的偏轉(zhuǎn)力矩，此力矩大于航行體受到的向外側(cè)偏轉(zhuǎn)的力矩，因此航行體開始向內(nèi)側(cè)偏轉(zhuǎn)，對應(yīng)偏航角減小?！鱠為4D和5D時，隨著并聯(lián)間距的增大，航行體受到相鄰航行體流場干擾的影響減弱，航行體向外側(cè)的偏航角減小，一直到航行體完全出水時偏航角還沒有達到臨界值，所以航行體只發(fā)生了向外側(cè)偏轉(zhuǎn)。因此為確保航行體出水彈道的穩(wěn)定性，需要合理控制并聯(lián)間距。

圖10 不同工況下航行體出水的偏航角變化Fig. 10 Variations of the yaw angle of vehicles exiting water under different cases

3 結(jié) 語

本文采用數(shù)值模擬方法對高速航行體并聯(lián)出水問題進行研究，獲得主要結(jié)論如下：

1）并聯(lián)出水過程中，航行體周圍流場彼此發(fā)生干擾，出水超空泡內(nèi)側(cè)擴張受到抑制，導(dǎo)致超空泡外側(cè)輪廓的曲率要大于內(nèi)側(cè)輪廓。而且當(dāng)并聯(lián)間距較小時，2 個出水超空泡的尾部發(fā)生融合。

2）相比航行體單獨出水，并聯(lián)出水過程中航行體的縱向速度衰減較快，使航行體的減阻性能小幅下降。

3）航行體并聯(lián)出水受到不對稱的水動力作用，航行體的彈道發(fā)生偏轉(zhuǎn)，當(dāng)并聯(lián)間距較小時，航行體先發(fā)生向外側(cè)偏轉(zhuǎn)之后轉(zhuǎn)向內(nèi)側(cè)偏轉(zhuǎn)，當(dāng)并聯(lián)間距較大時，航行體只發(fā)生向外側(cè)的偏轉(zhuǎn)。