王威，王聰，李聰慧，宋武超

(哈爾濱工業(yè)大學(xué) 航天學(xué)院，黑龍江 哈爾濱 150001)

0 引言

水下武器借助通氣超空泡的包裹，可以大幅度地減小水下阻力、提高航行速度。超空泡多相流動特性研究始終是各國學(xué)者關(guān)注的熱點(diǎn)問題，其中空泡尾流特性是通氣超空泡多相流動研究的重要方面，空泡尾部的雙渦管泄氣現(xiàn)象多在低弗勞德數(shù)Fr條件下產(chǎn)生，水洞實(shí)驗(yàn)也易于觀察，因此針對雙渦管尾流特性的研究備受各國學(xué)者青睞。

從工程應(yīng)用的角度看，帶有后體的超空泡航行體多相流特性也是需要考慮的問題，袁緒龍等[7]和陳偉政等[8]通過水洞實(shí)驗(yàn)分析了沾濕區(qū)域及重力對空泡形態(tài)的非對稱性影響。航行體后體與流動介質(zhì)相接觸[9]，沾濕區(qū)域是維持航行體受力平衡的關(guān)鍵。沾濕區(qū)域從空泡下部穿出，必然會影響空泡多相尾流結(jié)構(gòu)，這與無后體的自由閉合空泡尾流存在較大差異。Yu等[10]對航行體后體沾濕區(qū)域的升力進(jìn)行了研究，分析了升力突變的原因。Wang等[11]研究了超空泡航行體的尾流結(jié)構(gòu)，采用定常計算方法定性分析了空泡尾部流場壓力和渦量分布特性，指出了沾濕區(qū)域?qū)u量方向的影響，并給出了渦管內(nèi)部流動參數(shù)的變化規(guī)律。

從相關(guān)文獻(xiàn)報道可以看出，對圓盤空化器模型的研究多側(cè)重于空泡尾流結(jié)構(gòu)，對帶有后體的航行體模型研究多側(cè)重于流體動力，但研究后體沾濕區(qū)域?qū)叫畜w空泡尾流結(jié)構(gòu)的影響并不多見。數(shù)值模擬方法是研究超空泡尾流結(jié)構(gòu)的重要方法，其優(yōu)點(diǎn)是易于捕捉實(shí)驗(yàn)方法較難獲得的流場參數(shù)變化特性，分析流場參數(shù)非定常變化。

本文利用數(shù)值模擬方法，在文獻(xiàn)[11]的基礎(chǔ)上，采用非定常計算方法，進(jìn)一步定量分析低Fr條件下航行體攻角逐漸變化時沾濕區(qū)域形成過程中空泡尾流場壓力及渦量分布的演化特性，揭示沾濕區(qū)域?qū)张菸擦鹘Y(jié)構(gòu)的影響機(jī)理，給出沾濕區(qū)域比例與空泡尾流結(jié)構(gòu)的轉(zhuǎn)變關(guān)系。

1 數(shù)學(xué)模型及驗(yàn)證

1.1 航行體模型

根據(jù)低Fr條件(5

航行體模型如圖1所示，該模型由圓盤空化器、圓錐段、圓柱段、圓柱形尾噴管組成。航行體模型圓柱段直徑為D，航行體全長L=12D，錐段長度為L/2，圓盤空化器直徑Dn=0.3D，航行體質(zhì)心位置在xc=6.5D處，坐標(biāo)系原點(diǎn)建立在航行體質(zhì)心。

圖1 計算模型Fig.1 Computation model

1.2 計算流域邊界條件

計算流域的邊界如圖2所示(圖中p∞為遠(yuǎn)場邊界壓力)，計算域?yàn)閳A柱形，直徑為20D，長度為48D，重力方向沿y軸負(fù)向，計算流域采用六面體結(jié)構(gòu)化網(wǎng)格劃分。航行體表面采用壁面邊界，在空化器后部設(shè)有通氣孔，采用質(zhì)量入口條件；流域入口和流域圓柱面邊界都采用速度入口條件，流域出口采用水深5 m壓力出口條件。

圖2 計算域及邊界條件Fig.2 Computational domain and boundary condition

Zhou等[12]研究了無限流域尺度對數(shù)值模擬結(jié)果的影響，指出p∞的等值線應(yīng)在計算流域內(nèi)部。本文利用p∞的等值面判斷三維無限流域數(shù)值模擬對流域的最小尺度要求。從圖2可以發(fā)現(xiàn)，p∞的等值面并未與計算流域邊界交叉，且在計算流域內(nèi)部，表明本文的流域邊界尺度用于模擬無限流域情況下的通氣空泡流動是合理的。

比較不同網(wǎng)格密度條件下航行體全濕狀態(tài)時的流體動力系數(shù)隨網(wǎng)格數(shù)量Ne的變化，結(jié)果如圖3所示。圖3中，CD為阻力系數(shù)，CL為升力系數(shù)。

圖3 不同網(wǎng)格密度下航行體的流體動力系數(shù)Fig.3 Hydrodynamic coefficients of vehicle under different grid densities

從圖3中可以看出：在0.6×106～0.9×106網(wǎng)格密度條件下，航行體的阻力系數(shù)CD和升力系數(shù)CL均隨著網(wǎng)格數(shù)量而變化；在0.9×106～1.2×106的網(wǎng)格密度條件下，航行體的阻力和升力系數(shù)不再隨網(wǎng)格數(shù)量變化。采用0.9×106以上網(wǎng)格數(shù)量進(jìn)行數(shù)值模擬時，計算結(jié)果不會因?yàn)榫W(wǎng)格疏密程度而造成計算誤差，綜合考慮數(shù)值計算速度等因素，本文選取網(wǎng)格數(shù)量Ne=1.0×106的條件進(jìn)行數(shù)值模擬研究。

1.3 模型驗(yàn)證

本文的通氣率系數(shù)為

Cq=Q/(v∞D(zhuǎn)2)，

(1)

式中：Q為氣體體積流量。

在相同的Cq和v∞條件下，將文獻(xiàn)[10]的水洞實(shí)驗(yàn)結(jié)果和本文數(shù)值模擬結(jié)果對比，如圖4和圖5所示。圖4分別給出了數(shù)值模擬和實(shí)驗(yàn)空泡形態(tài)發(fā)展過程，上方是數(shù)值模擬結(jié)果(本文以空氣體積分?jǐn)?shù)為0.5作為空泡界面)，下方是實(shí)驗(yàn)結(jié)果，可見空泡尾部都出現(xiàn)了渦管泄氣的現(xiàn)象。

圖4 數(shù)值模擬結(jié)果和實(shí)驗(yàn)結(jié)果對比Fig.4 Comparison of simulated and experimental results

圖5 數(shù)值模擬和實(shí)驗(yàn)的空泡尺度對比Fig.5 Comparison of simulated and experimental cavity dimensions

空泡尺度的對比結(jié)果如圖5所示。由圖5可以看出，空泡長度Lc和空泡直徑Dc的變化曲線具有很好的一致性，并且數(shù)值模擬結(jié)果和實(shí)驗(yàn)結(jié)果的空泡尺度相差在6%以內(nèi)。圖4和圖5表明，數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果吻合較好，選用的數(shù)學(xué)模型可以很好地預(yù)測通氣空泡多相流動現(xiàn)象。

2 模擬結(jié)果分析與討論

2.1 空泡尾流結(jié)構(gòu)

航行體在不同工況下的空泡形態(tài)如表1所示。

表1 不同工況下航行體的空泡形態(tài)

Tab.1 Supercavity shapes under different working conditions

比較表1中工況1和工況2，發(fā)現(xiàn)：在無沾濕區(qū)域條件下，通氣超空泡的尾流在重力環(huán)境中呈現(xiàn)雙渦管的泄氣形式；工況1僅起對比作用，表明在給定流速條件下，不計重力環(huán)境的壓差因素時，空泡尾部為單一渦管的流動結(jié)構(gòu)。

航行體后體沾濕區(qū)域主要為航行體提供縱平面內(nèi)平衡的升力作用[13-14]。比較表1中工況2、工況3、工況4可以發(fā)現(xiàn)，隨著攻角α的增加，航行體沾濕區(qū)域面積逐步增大，航行體超空泡尾流結(jié)構(gòu)從雙渦發(fā)展到三渦，表明沾濕區(qū)域?qū)张菸擦鹘Y(jié)構(gòu)存在影響。

表1中工況5僅起與工況4的對比作用，給出了不計環(huán)境壓差因素時航行體的超空泡尾部形態(tài)。比較工況4、工況5可以發(fā)現(xiàn)，在相同攻角(α=1.8°)條件下，兩種工況中航行體的超空泡形狀和沾濕區(qū)域都非常相似。二者環(huán)境壓差不同，進(jìn)一步表明航行體沾濕區(qū)域是影響超空泡尾部流場結(jié)構(gòu)的重要因素。

航行體直線運(yùn)動時空泡形態(tài)演化過程如圖6所示。由圖6可見，0.50 s以后，通氣超空泡已經(jīng)包裹整個航行體，形態(tài)尺度趨于穩(wěn)定，空化數(shù)σ≈0.028，此時壓力曲線平直，整個航行體表面壓力分布趨于平穩(wěn)，空泡內(nèi)部壓力差別較小[11]。

圖6 航行體空泡演化Fig.6 Evolution of cavity of vehicle

本文對壓力系數(shù)Cp及渦量Ω分別定義如下：

(2)

(3)

式中：p為流域壓力；Ωx、Ωy、Ωz分別表示x軸、y軸、z軸方向的渦量，渦量的符號按照右手螺旋定則取值；vx、vy、vz分別表示x軸、y軸、z軸方向的流動速度。

空泡在0.30 s形成自由閉合，航行體從0.30 s開始勻速改變攻角至0.34 s結(jié)束，航行體的攻角α從0°增加到1.8°. 航行體尾流場的壓力及渦量變化如圖7所示。

圖7 空泡尾部壓力及渦量云圖Fig.7 Contours of pressure and vorticity at the tail of supercavity

統(tǒng)計圖7中空泡尾部流場最右側(cè)截面在航行體攻角改變過程中空泡尾渦的非定常變化，結(jié)果如表2所示。表2中S0為沾濕區(qū)域面積，S為航行體表面積。

表2 空泡尾流橫截面渦方向的非定常變化

Tab.2 Unsteady variation of vortex direction in cross section of cavity wake

從表2中可以看出：t取值為0.300～0.310 s時，空泡尾部的雙渦管逐漸形成，流域環(huán)境壓差使空泡尾部主要以雙渦管閉合為主，雖然有小部分沾濕區(qū)域影響，但S0/S只有0.87%左右，空泡尾部的渦量方向在此期間保持一致，為逆時針和順時針，如圖7(b)所示；t取值為0.310～0.319 s時，空泡尾部雙渦逐漸消失，隨著S0/S的不斷增大，空泡尾部形成四渦，該情況在圣安東尼瀑布實(shí)驗(yàn)室的研究中也有出現(xiàn)[5]，實(shí)驗(yàn)指出四渦情況發(fā)生在超空泡尾流從回射流到雙渦管之間的轉(zhuǎn)變過程中，由于數(shù)值模擬研究中選擇空氣體積分?jǐn)?shù)不同以及數(shù)值方法在處理細(xì)微氣泡的旋渦結(jié)構(gòu)時還與實(shí)驗(yàn)結(jié)果存在一定差異，故數(shù)值模擬暫時還不能精細(xì)地反映實(shí)驗(yàn)中渦管細(xì)節(jié)；t取值為0.319～0.321 s時，隨著沾濕區(qū)域?qū)翰钣绊懙脑龃螅臏u中的上方雙渦逐漸減小，下方雙渦逐漸增大，渦量云圖如圖7(c)所示；t取值為0.321～0.322 s時，隨著S0/S的進(jìn)一步擴(kuò)大，空泡尾流上方雙渦逐漸被下方雙渦取代；t取值為0.322～0.340 s時，空泡尾部上方的渦管尺度逐漸增大，形成三渦管尾流，上方渦管流域的渦量很小，下方雙渦的渦量保持順時針和逆時針方向，如圖7(d)所示；t取值為0.340～0.440 s時，航行體攻角保持固定，空泡形態(tài)趨于穩(wěn)定，尾渦形態(tài)如圖7(e)所示。

2.2 空泡尾部壓力

航行體尾部渦管附近的流線方向如圖8所示，從圖8(a)和圖8(b)中的流線方向?qū)Ρ瓤梢园l(fā)現(xiàn)，沾濕區(qū)域致使空泡尾部橫截面中渦管附近的流線方向發(fā)生了改變。

圖8 空泡尾部渦管附近的壓力和流線Fig.8 Pressure and streamline near vortex tube at cavity tail

圖9 攻角、沾濕區(qū)域和空泡尾流高壓區(qū)無量綱壓力差系數(shù)隨時間變化曲線Fig.9 Changes of angle-of-attack, wetted area dimensionless pressure difference in the high pressure region of wake flow of cavity with time

由圖9可以看出：t取值為0.300～0.318 s時，受流域壓差影響，上方高壓區(qū)的壓力小于下方，雙渦管之間的流動方向向上，在空泡尾部形成一對中心上卷的雙渦；t取值為0.318～0.322 s時，航行體的沾濕區(qū)域開始逐漸影響尾流壓力分布，使上下方高壓區(qū)逐漸融合，壓力差減??；t取值為0.322～0.340 s時，攻角不斷增加，沾濕區(qū)域繞流的誘導(dǎo)作用對尾流壓力分布的影響進(jìn)一步擴(kuò)大，高壓區(qū)再次分成上下兩個部分，上方高壓區(qū)的壓力大于下方，壓力差值逐漸增加，雙渦管之間的介質(zhì)轉(zhuǎn)為向下流動，在空泡尾部形成一對中心下卷的雙渦；t取值為0.340 s后，航行體攻角不在變化，沾濕區(qū)域趨于穩(wěn)定。

2.3 空泡尾部渦量

為進(jìn)一步分析四渦出現(xiàn)期間沾濕區(qū)域?qū)ξ矞u方向的影響規(guī)律，考慮到流域環(huán)境壓差和沾濕區(qū)域?qū)张菸矞u存在相反方向的影響，可以通過空泡尾部上下流線的速度分布進(jìn)行深入研究[7-8]，如圖10所示?？张菘v向平面上圍繞整個空泡閉合輪廓線的環(huán)量Γ[15]近似表示為

(4)

式中：v為空泡輪廓線速度矢量；s為空泡輪廓線長度矢量；l為空泡輪廓線周長；vx1和vx2分別為空泡尾部上、下表面流線的x軸方向速度。

圖10 空泡尾部上下流線的速度vx1和vx2Fig.10 Velocities vx1 and vx2 of upper and lower streamlines at cavity tail

由于空泡上下表面的速度差異在空泡尾部處較大，用vx1曲線和vx2曲線所圍的面積可以近似地表征空泡尾部的渦量變化，如圖11所示。圖11(a)中空泡剛好閉合在航行體尾部，受流域環(huán)境壓差作用，空泡下表面的速度略大[11]，形成雙渦尾流。圖11(b)中空泡尾流既受環(huán)境壓差影響，又受到不斷擴(kuò)大的沾濕區(qū)域影響，空泡下表面速度vx2受沾濕區(qū)域阻擋逐漸小于vx1，渦量方向(順時針區(qū)域和逆時針區(qū)域的面積之差)變化較快，空泡尾部形成了四渦流動結(jié)構(gòu)。圖11(c)中，空泡尾流隨著沾濕區(qū)域的進(jìn)一步擴(kuò)大，空泡上表面的速度逐漸大于下表面，渦量方向發(fā)生變化，再次形成方向相反的雙渦尾流。

圖11 空泡尾部上下流線的速度比較(g=9.8 m/s2)Fig.11 Comparison of velocities of upper and lower streamlines at cavity tail(g=9.8 m/s2)

進(jìn)一步研究沾濕區(qū)域和流域環(huán)境壓差對空泡尾渦影響的主次關(guān)系，對比分析了不考慮環(huán)境壓差條件下(如工況1，假設(shè)g=0)航行體攻角變化時，沾濕區(qū)域?qū)张菸膊繙u量的影響。如圖12所示，0.300 s時航行體攻角α=0°，空泡尾部不受沾濕區(qū)域影響，航行體空泡尾部上下速度差異較小，空泡橫截面上有4個尾渦。0.340 s時航行體攻角α=1.8°，沾濕面積比例為36%左右，沾濕區(qū)域影響了航行體尾部的流速分布，空泡尾部出現(xiàn)對稱分布的渦量。

圖12 空泡尾部上下流線的速度比較(g=0 m/s2)Fig.12 Comparison of velocities of upper and lower streamlines at cavity tail(g=0 m/s2)

將vx1和vx2沿x軸方向積分，兩條曲線所圍的無量綱面積A為

(5)

用A的統(tǒng)計結(jié)果分析比較考慮沾濕區(qū)域引起的流速差異和環(huán)境壓差引起的流速差異，比較二者對空泡尾部渦量方向的影響關(guān)系，如圖13所示。

圖13 沾濕區(qū)域?qū)张菸矞u的影響規(guī)律Fig.13 Influence law of wetted area on cavity wake vortex

從圖13中可以看出，沾濕區(qū)域影響的曲線和流域環(huán)境壓差影響的曲線相交在S0/S處于7～10%之間：當(dāng)S0/S<7%時，環(huán)境壓差對空泡尾部渦量方向起到主導(dǎo)作用，呈現(xiàn)一對中心上卷的雙渦；當(dāng)7%10%時，沾濕區(qū)域繞流的誘導(dǎo)作用對空泡尾部渦量方向起到主導(dǎo)作用，呈現(xiàn)一對中心下卷的雙渦。另外，由于流動過程中存在多相流之間復(fù)雜的擾動影響，兩條曲線交點(diǎn)處S0/S的數(shù)值是一個范圍。

3 結(jié)論

本文研究了超空泡航行體沾濕區(qū)域?qū)张菸擦鹘Y(jié)構(gòu)的影響，對空泡尾流結(jié)構(gòu)的壓力及渦量的非定常變化進(jìn)行了深入分析，得到了沾濕區(qū)域?qū)张菸擦鹘Y(jié)構(gòu)的影響規(guī)律。主要得出以下結(jié)論：

1)揭示了低弗勞德數(shù)條件下通氣超空泡尾流結(jié)構(gòu)特性。航行體無攻角時，空泡尾部形成一對中心上卷的雙渦；航行體帶有攻角的條件下，空泡尾部形成一對中心下卷的雙渦。

2)給出了沾濕區(qū)域?qū)ξ矞u方向的影響規(guī)律。在考慮重力環(huán)境影響的條件下，沾濕區(qū)域比例小于7%時，環(huán)境壓差是空泡尾渦方向的主導(dǎo)因素；沾濕區(qū)域比例大于10%時，沾濕區(qū)域是空泡尾渦方向的主導(dǎo)因素。

3)指出了空泡的四渦尾流結(jié)構(gòu)產(chǎn)生的原因。沾濕區(qū)域的比例處于7%～10%之間時，在環(huán)境壓差和沾濕區(qū)域的雙重作用下，空泡尾部流場中出現(xiàn)了四渦流動結(jié)構(gòu)。

致謝本文撰寫過程中得到了哈爾濱工業(yè)大學(xué)航天學(xué)院曹登慶教授和徐敏強(qiáng)教授的大力幫助與支持，在此對兩位老師表示衷心的感謝！