胡曉，郜冶，熊永亮

(哈爾濱工程大學(xué)航天與建筑工程學(xué)院，黑龍江哈爾濱150001)

航行體在高速入水和水下高速運(yùn)動過程中，會產(chǎn)生空泡(或超空泡)，空泡的外形和尺寸會對航行體的水動力特性產(chǎn)生重大影響，從而影響實際的水彈道以及對航行體的控制［1-2］.

要維持理想的流動體制，對超空泡外形的控制就非常必要，目前國內(nèi)公開的文獻(xiàn)中主要是通過改變空化器線性和外形參數(shù)來控制空泡形態(tài)［3-5］，Savchenko Yu.N對錐形空化器和Σ形空化器開展了研究，給出了它們的阻力特性與形狀參數(shù)和空化數(shù)的函數(shù)關(guān)系式，同時還對Σ形空化器的穩(wěn)定性進(jìn)行了深入的研究分析［6］;或者通過采用非零攻角的空化器來實現(xiàn)非對稱超空泡流動模式，從而達(dá)到控制水下高速航行體運(yùn)動姿態(tài)和受力的目的，Waid［7］和Semenneko［8］對帶攻角的空化器形成的空泡形態(tài)進(jìn)行了研究，黃海龍等［9］通過數(shù)值模擬的方法詳細(xì)分析了不同攻角的圓盤空化器生成自然超空泡的形態(tài)特征及升阻力特性.

當(dāng)航行體運(yùn)動速度增加并接近蒸汽空泡狀態(tài)時，可以通過改變空化器阻力系數(shù)cx0來控制超空泡參數(shù)，隨著空化器阻力系數(shù)的變化，超空泡尺寸也會相應(yīng)改變.

本文對不同中心錐角的變阻力空化器在空泡外形的控制能力進(jìn)行了數(shù)值研究，比較了不同中心錐角空化器的阻力系數(shù)和空泡尺寸響應(yīng)速度;并提出了一種新的空泡控制方案，數(shù)值模擬結(jié)果也證明了該方案在阻力控制方面的靈敏性和可操作性.

1 控制方程與數(shù)值模型

1.1 控制方程

利用 FLUENT6.2中的 Mixture多相流模型［10］，對混合相的連續(xù)性方程和動量方程進(jìn)行求解.

混合相連續(xù)性方程:

式中:xi為笛卡爾坐標(biāo)向量，ρm為混合介質(zhì)密度:

式中:αk、ρk分別表示第k相的體積分?jǐn)?shù)和密度，ui

為混合物的質(zhì)量平均速度:

式中:uk為第k相的速度矢量.

混合相動量方程:

式中:p為壓力，μt為湍流粘性系數(shù)，μm為混合粘性系數(shù)，其表達(dá)式如下:

式中:μk為第k相的粘度.

湍流模型采用標(biāo)準(zhǔn)k-ε兩方程模型，近壁區(qū)域的流動采用標(biāo)準(zhǔn)壁面函數(shù)法處理.

1.2 自然空化模型

FLUENT6.2中的自然空化模型［11］模擬的空泡流包括水和水蒸汽介質(zhì)，空泡流建模采用混合均質(zhì)流理論，整個混合物允許對流，蒸汽相傳輸方程如下:

式中:f為蒸汽相質(zhì)量分?jǐn)?shù)，Re和Rc是蒸汽產(chǎn)生和潰滅項:

式中:下標(biāo)l和v分別表示液相和汽相，σ為液體的表面張力，Ce和Cc為經(jīng)驗常數(shù)，其值分別為0.02和 0.01.

1.3 動網(wǎng)格技術(shù)

采用動態(tài)分層法對變形區(qū)域網(wǎng)格進(jìn)行更新，與動邊界相鄰的網(wǎng)格層(圖1中的層j)依據(jù)其網(wǎng)格高度h被切割或與下一網(wǎng)格層(圖1中的層i)合并，具體法則如下:

hmin＞(1+αs)hideal，j層網(wǎng)格被分割;hmin＜αchideal，j層網(wǎng)格被合并(潰滅).

其中，hmin為網(wǎng)格層j的最小網(wǎng)格高度，hideal為理想網(wǎng)格高度，αs、αc分別為網(wǎng)格層分割和潰滅因子.

圖1 動態(tài)分層法Fig.1 Dynamic layering method

2 數(shù)值計算

2.1 傳統(tǒng)變阻力空化器

圖2是文獻(xiàn)［12-13］中提到的變阻力空化器.空化器由帶尖銳前緣的外腔體1，以及安裝于腔體內(nèi)部，與1有著相同直徑Dn的中心錐形元件2組成.通過改變外體1和中心元件2的相對位置，可以使空化器在不改變直徑的情況下阻力發(fā)生變化.

圖2 傳統(tǒng)變阻力空化器示意圖Fig.2 Scheme of the traditional variab le-drag cavitator

試驗［6］研究表明，空化器阻力系數(shù)cx和工作行程x的變化幅度均隨中心元件錐角2β的減小而增大，可以看出，中心錐角2β=60°時，在工作行程內(nèi)可得到阻力系數(shù)幾乎3倍的變化.

由于以上試驗是基于穩(wěn)態(tài)過程進(jìn)行的，并不能反映超空泡外形隨阻力變化的瞬態(tài)響應(yīng)，本文展開了基于動網(wǎng)格技術(shù)的空泡流瞬態(tài)數(shù)值計算.圖3為中心元件錐角2β=90°時變阻力空化器附近的網(wǎng)格劃分.

計算中設(shè)置入口速度為30 m/s，空化數(shù)為0.174.初始時刻，中心元件最左端與外體尖緣部分位于同一豎直平面(x=h)，首先進(jìn)行自然空泡流的穩(wěn)態(tài)計算，并以此作為動態(tài)計算的初始流場，之后中心元件在12 s內(nèi)以uc=0.1m/s的速度往前移動，其中空泡面定義為蒸汽相體積分?jǐn)?shù)為0.5的等值面.

圖3 90°錐角變阻力空化器附近網(wǎng)格Fig.3 Grid around the variable-drag cavitator，2β =90°

2.2 變阻力空化器瞬態(tài)空泡流模擬結(jié)果分析

2.2.1 不同中心錐角空化器的空泡流計算結(jié)果

圖4給出了2β=60°、90°的變阻力空化器在不同工作行程x處的空泡外形.可以看出:1)雖然空化器一直向前運(yùn)動，但在中心元件完全露出來之前(0≤x≤h)，空泡長度幾乎沒有變化;2)隨著中心元件的繼續(xù)運(yùn)動，其尖緣部分開始空化，空泡長度也逐漸減小，達(dá)到穩(wěn)定后不再變化;3)中心錐角越大，穩(wěn)定狀態(tài)所對應(yīng)的工作行程越大:2β=90°時，空泡長度在x=-0.8Dn之后便不再變化，而且由于空化器阻力在x=-0.3Dn處的突然降低，在原來較大的空泡與新生成的較小空泡之間形成了明顯的階梯式邊界(見圖4(b))，這與文獻(xiàn)［6］中觀察到的空泡形態(tài)相符.

圖4 不同工作行程處的空泡外形Fig.4 Cavity shapes at different working stokes

圖5給出了60°、90°和120°中心錐角空化器的阻力系數(shù)計算結(jié)果(其中阻力系數(shù)Cd=F/0.5ρv2S，F(xiàn)為阻力，S為空化器的最大橫截面積).圖5(a)為90°中心錐角空化器的阻力系數(shù)計算值與試驗數(shù)據(jù)的對比，由曲線可以看出，模擬值與試驗值變化趨勢大體一致，與模擬結(jié)果相比，穩(wěn)態(tài)試驗中的阻力隨工作行程的增加，其增長速度更為緩慢.由圖5(b)可以看出，中心元件的運(yùn)動使空化器的阻力系數(shù)發(fā)生了明顯變化，其變化的幅度隨中心錐角的減小而增大，其中60°錐角的變阻力空化器在工作行程內(nèi)可得到阻力系數(shù)約2.4倍的變化，這與試驗數(shù)據(jù)是相吻合的;阻力系數(shù)在初始時刻最大，其最大值與圓盤空化器的阻力系數(shù)相接近;隨著中心錐角的減小，空化器對應(yīng)的阻力系數(shù)下限值也有所降低.

圖5 阻力系數(shù)隨工作行程的變化曲線Fig.5 Curves of drag coefficients of the cavitator on working stokes

圖6為無量綱空泡尺寸的變化曲線，由圖6知，空泡長度Lc和最大直徑Dc的變化趨勢基本一致:即在x=0.0Dn處空泡尺寸由最大值開始減小，x=-0.80Dn處，空泡長度及最大直徑降至最低，之后不再變化;Lc、Dc的變化幅度隨錐角的減小而增大:2β=60°時，工作行程內(nèi)可得到空泡長度1.15倍和最大直徑1.29倍的變化.不同中心錐角的變阻力空化器可獲得的最大空泡尺寸與中心元件形狀無關(guān)，當(dāng)工作行程值x接近中心元件高度h時(見圖2(a))，空泡尺寸達(dá)到最大.

結(jié)合阻力系數(shù)和空泡尺寸的變化曲線可以看出，兩者變化并不同步，空泡長度及最大直徑變化要滯后于阻力的變化，初始時刻(x=h)，阻力系數(shù)和空泡尺寸取得最大值，之后隨著工作行程值x的逐漸減小，阻力系數(shù)開始減小，在x=-0.50Dn達(dá)到最小，而空泡長度及最大直徑在x=0.0Dn處才開始變化，在x=-0.80Dn處降至最低.

圖6 不同錐角變阻力空化器的空泡尺寸變化曲線Fig.6 Curves of variance of cavity size for the variable-drag cavitatorswith different cone angles

2.2.2 中心元件速度對模擬結(jié)果的影響

為研究中心元件的運(yùn)動速度對空泡流場結(jié)構(gòu)的影響，本文選取中心錐角2β=60°的變阻力空化器，對其中心元件速度uc=1 m/s時的瞬態(tài)空泡流進(jìn)行了模擬計算.

圖7為不同uc下的阻力系數(shù)Cd曲線，由圖知2種速度下的阻力系數(shù)變化趨勢基本一致:即隨著工作行程值x的減小，阻力系數(shù)逐漸下降，降至最低值后不再變化，中心元件速度為1 m/s時，計算的阻力系數(shù)值要略大些.

圖7 阻力系數(shù)隨工作行程變化曲線，2β=60°Fig.7 Curves of drag coefficient on working stoke，2β =60°

空泡尺寸隨工作行程變化曲線如圖8，隨著工作行程值x的減小，空泡長度Lc和最大直徑Dc并沒有馬上減小，而是經(jīng)過一段時間的遲滯才開始變化;由于與速度為0.1 m/s的計算結(jié)果相比，中心元件速度為1 m/s時得到的阻力系數(shù)值要稍微大些，故該速度下計算的空泡尺寸也有小幅增加，由于遲滯效應(yīng)，1 m/s速度下計算的空泡長度在x=0.466Dn時開始增加，直到工作行程值x=-0.80Dn，2種中心元件運(yùn)動速度下的空泡長度不再變化;而空泡最大直徑的增加與阻力系數(shù)同步.

圖8 不同中心元件速度的空泡尺寸隨工作行程變化曲線，2β=60°Fig.8 Curves of cavity size on working stoke at different velocity magnitudes of central element，2β =60°

結(jié)合圖7、8可以看出，中心元件運(yùn)動速度地大小對變阻力空化器的空泡流模擬結(jié)果有一定影響，隨著速度的增大，空泡尺寸及阻力系數(shù)值均有所增加，其中空泡長度的增加幅度最大達(dá)到了10%，而空泡最大直徑和阻力系數(shù)受速度影響較小.

2.3 改進(jìn)的可變側(cè)向力空化器

2.3.1 控制方案與內(nèi)部結(jié)構(gòu)

上述變阻力空化器雖然可以比較快速的改變航行體阻力和空泡尺寸，然而要想獲得良好的控制效果，需要其控制機(jī)構(gòu)能夠輕易提供數(shù)十萬牛頓的力，這一控制方法很難實施.考慮到空化器側(cè)面的壓力較低，而且通過增加空化器有效面積同樣可以改變航行體的受力，本文提出了如圖9所示的可變側(cè)向力空化器的控制方案:在錐形空化器內(nèi)部安裝一外徑小于Dn的可動元件3，通過可動元件3的橫向運(yùn)動來控制元件2的位置，元件2、3均設(shè)計成4瓣，控制元件2內(nèi)部共有8對弧形鋼片4，其中每瓣的兩端均有一對弧形鋼片，可動元件3由錐形可動部分和圓柱操縱桿連接而成，4瓣圓柱操縱桿分別控制4瓣圓錐部分的橫向位移，從而改變控制元件的周向位移hn值，使前端空化器1的有效工作面積變化，進(jìn)而改變空泡尺寸;當(dāng)可動元件任意兩側(cè)位移不一致時，控制元件具有不對稱性，前端錐形空化器產(chǎn)生不對稱受力，使空泡形態(tài)發(fā)生變化，可以有效的控制航行體運(yùn)動姿態(tài)，有關(guān)可變側(cè)向力空化器的具體內(nèi)容可參考文獻(xiàn)［14］.

圖9 可變側(cè)向力空化器的控制方案Fig.9 Control scheme of the variable-lateral force cavitator

2.3.2 空泡流模擬結(jié)果

圖10、11給出了45°錐角空化器的控制元件2上半周部分突出1/10Dn時的空泡流計算結(jié)果，由于元件2的不對稱性，在空化器上產(chǎn)生豎直向下的升力，導(dǎo)致空泡上下截面不對稱(如圖10所示).3種空化數(shù)下得到的空化器阻力系數(shù)分別為0.620、0.553和0.544，相比于圓錐空化器，可變側(cè)向力空化器在相同空化數(shù)下的阻力系數(shù)增加了近60%，證實了其在阻力控制方面的靈敏性.

圖10 可變側(cè)向力空化器的空泡形態(tài)Fig.10 Cavity shapes of the variable-lateral force cavitator

空化數(shù)為0.190時的空化器前端縱剖面的壓力系數(shù)分布如圖11所示，可以看出，凸出的元件2導(dǎo)致了不對稱的壓力分布，從而在空化器上產(chǎn)生了相當(dāng)于30%阻力的升力，升力系數(shù)達(dá)到0.177;此外壓力系數(shù)等值線顯示，控制元件2外表面的壓力非常低，這保證了元件的可操作性.

圖11 可變側(cè)向力空化器前端壓力系數(shù)等值線，σ=0.190Fig.11 Contours of pressure coefficient in front of the variable-lateral force cavitator，σ =0.190

3 結(jié)論

本文利用FLUENT6.2的自然空化模型和動網(wǎng)格技術(shù)對水下超空泡航行體的瞬態(tài)空泡流進(jìn)行了數(shù)值模擬，得到了以下幾點結(jié)論:

1)傳統(tǒng)變阻力空化器的阻力響應(yīng)速度較快，空泡外形的變化則要滯后于阻力變化;隨著工作行程值的減小，不同錐角的變阻力空化器的阻力系數(shù)和空泡尺寸也逐漸減小;

2)中心錐角越小，對應(yīng)的空化器工作行程及阻力系數(shù)(空泡尺寸)的變化幅度也就越大，2β=60°時在工作行程內(nèi)可得到阻力系數(shù)幾乎3倍的變化;

3)中心元件運(yùn)動速度對模擬結(jié)果也有一定影響，速度值從0.1m/s增至1m/s時，空泡尺寸有小幅增加;

4)數(shù)值模擬結(jié)果顯示可變側(cè)向力空化器對阻力和升力的控制非常靈敏.

［1］MAY A.Water entry and cavity-running behavior ofmissles［R］.AD-A020259，1975.

［2］WOLFE wP.Drag of revolution in supercavitating flow［J］.Journal of Fluids Engineering，1989，112:23-30.

［3］賈力平，張嘉鐘，于開平，等.空化器線形與超空泡減阻效果關(guān)系研究［J］.船舶工程，2006，28(2):20-23.JIA Lingping，ZHANG Jiazhong，YU Kaiping，et al.Research on the relation of cavitator profile and supercavity dragreduction［J］.ShipEngineering，2006，28(2):20-23.

［4］顧建農(nóng)，張志宏，高永琪，等.充氣頭型對超空泡軸對稱體阻力特性影響的試驗研究［J］.兵工學(xué)報，2004，25(6):766-769.GU Jiangnong，ZHANG Zhihong，GAO Yongqi，et al.An experimental study about the influence of ventilated headforms on the drag characters of a super cavitation axisymmetric body［J］.Acta Armamentarll，2004，25(6):766-769.

［5］賈力平，王聰，于開平，等.空化器參數(shù)對通氣超空泡形態(tài)影響的實驗研究［J］.工程力學(xué)，2007，24(3):159-164.JIA Liping，WANG Cong，YU Kaiping，etal.Experimental investigation of cavitator parameters effecting on ventilated supercavity shape［J］.Engineering Mechanics，2007，24(3):159-164.

［6］SAVCHENKO Y N.Control of supercavitation flowand stability of supercavitating motion of bodies［C］//The RTO AVT Lectures Series on"Supercavitating Flows"，RTO Lecture Series 005.Ottawa:St.Joseph，2002:313-341.

［7］WAID R L.Cavity shapes for circular disks at angles of attack［R］.California Institute of Technology.Report No.E-73.4，1975.

［8］SEMENENKO V N.Artificial supercavitation physics and calculation［C］//RTO/AVT lecture series on"supercavitating flows"，RTO Lecture Series 005，Ottawa:St.joseph，2002:195-227.

［9］黃海龍，王聰，黃文虎，等.變攻角圓盤空化器生成自然超空泡數(shù)值模擬［J］.北京理工大學(xué)學(xué)報，2008，28(2):100-103.HUANG Hailong，WANG Cong，HUANG Wenhu，et al.Numerical simulation on natural supercavity based on diskcavitator with variable attack angle［J］.Transactions of Beijing Institute of Technology，2008，28(2):100-103.

［10］MANNINEN M，TAIVASSALO V，KALLIO S.On the mixturemodel for multiphase flow［M］.Espoo:Technical Research Center of Finland，VTT Publications，1996:3-67.

［11］SINGHALA K，ATHAVALEmM，LIH Y，etal.Mathematical basis and validation of the full cavitation model［J］.Journal of Fluids Engineering，2002，124:617-624.

［12］FU Huiping，LU Chuanjing，LIJie.Numercial research on drag reduction characteristics of supercavitating body of revolution［J］.Journal of ShipMechanics，2004，8(3):1-8.

［13］SAVCHENKO Y N，VLASENKO Y D，SEMENENKO V N.Experimental study of high-speed cavitation flows［J］.Jof Fluid Mechanics Research，1998，26(3):103-111(in Russian).

［14］哈爾濱工程大學(xué).可變側(cè)向力空化器［P］.中國:20091007，2007.3，2011，4.21.