于開平，張 廣，鄒 望，李振旺

（哈爾濱工業(yè)大學 航天學院，哈爾濱 150001）

超空泡航行體縱向平面動力學行為的CFD分析

于開平，張 廣，鄒 望，李振旺

（哈爾濱工業(yè)大學 航天學院，哈爾濱 150001）

超空泡航行體穩(wěn)定性及彈道特性的仿真研究是設(shè)計超空泡航行體控制和制導系統(tǒng)的基礎(chǔ)。文章基于多相流URANS（Unsteady Reynolds Averaged Navier-Stokes）方程和剛體縱向平面運動學方程，建立了三維超空泡航行體動力學仿真模型，對縱向平面內(nèi)超空泡航行體巡航狀態(tài)的動力學行為進行了仿真研究。仿真結(jié)果表明，超空泡航行體在無控條件下仍具有一定的穩(wěn)定性；尾部控制面可以有效縮短航行體達到穩(wěn)定所需要的時間，并且能使航行體穩(wěn)定于空泡中心構(gòu)成理想的穩(wěn)定模式。

流體力學；超空泡航行體；動力學行為；CFD

1 引 言

超空泡技術(shù)為大幅提高水下航行體速度帶來了光明的前景，但是由于超空泡航行體系統(tǒng)本身的復雜性，使得該項技術(shù)在實際的應用中面臨著諸多挑戰(zhàn)，其中超空泡航行體的穩(wěn)定性和控制問題尤為突出[1-2]。

近些年來，國內(nèi)外學者基于表征空泡形態(tài)的數(shù)學模型[3-4]，通過相應模型的簡化和假設(shè)，構(gòu)建了系列超空泡航行體的動力學模型并探討了相關(guān)控制方法。所研究的問題涉及單自由度模型，三自由度模型和六自由度模型。Kirschner[5]和Kulkarni[6]基于簡單的單自由度模型，研究了無控超空泡射彈的動力學行為，得到了航行體水平航行的一些典型特征參數(shù)。

文獻[7-8]考慮了空泡的記憶效應以及航行體尾部滑行的非線性作用，建立了超空泡航行體非線性動力學模型。Goel[9]對一個有固定空泡外形的六自由度航行體進行了研究，對模型進行了基于小擾動的線性化后分析了穩(wěn)定性，得到了縱向和橫向無控運動的穩(wěn)定性結(jié)論。從以上的研究工作來看，在超空泡航行體動力學建模過程中，空泡形態(tài)和相應的流體動力項主要采用經(jīng)典的勢流理論經(jīng)驗公式構(gòu)建，對模型進行了不同程度的線性化處理，這樣導致的結(jié)果是難以精確預測航行體與空泡之間復雜的相互作用，影響數(shù)值計算精度。

隨著計算機水平的不斷提高，計算流體動力學（CFD）方法已經(jīng)成為預測和揭示一些復雜流動現(xiàn)象的有效手段。本文在文獻[10-12]工作的基礎(chǔ)上，通過對商用CFD軟件CFX進行二次開發(fā)，耦合多相流URANS方程和剛體縱向平面運動學方程，結(jié)合動網(wǎng)格技術(shù)，建立了三維超空泡航行體縱向平面數(shù)值計算模型。數(shù)值模型綜合考慮了航行體和航行體控制面與空泡表面復雜的非線性力作用，空泡形態(tài)的三維瞬態(tài)特性及湍流影響，對超空泡航行體巡航階段的動力學行為及彈道特性進行數(shù)值研究。

2 數(shù)值計算模型

2.1 基本控制方程

本文研究自然空泡流動問題，采用均質(zhì)平衡流模型，假設(shè)汽水兩相具有相同的速度場和壓力場，兩相之間不存在相對速度。

連續(xù)性方程：

動量方程：

其中：ρm，μm分別為按體積分數(shù)確定的混合物的密度和動力粘度，u為速度，SM為源項。

2.2 湍流模型［10］

針對本文所涉及的問題，在穩(wěn)態(tài)數(shù)值仿真中使用SST（Shear Stress Transport）湍流模型，而非穩(wěn)態(tài)仿真采用能夠精確預測湍流影響的LES模型，分別簡要介紹如下：

基于Baseline（BSL）k-ω模型的SST湍流模型解決了湍流剪切應力的傳輸問題，同時又對逆壓梯度下產(chǎn)生的分離流具有較高精度的預測，基本方程如下：

其中：ρ為密度、U為速度矢量、Pk表示湍流的生產(chǎn)率，其它參量的具體表達式本文不再贅述。

LES（Large Eddy Simulations）模型的基本思想是認為湍流流場是由不同尺度的渦旋組成，通過濾波函數(shù)把流場瞬時變量分為大尺度渦運動和小尺度渦運動。大尺度渦通過直接求解瞬態(tài)N-S方程模擬，不直接模擬小尺度渦，小尺度渦對大尺度渦的影響采用亞格子模型（Subgrid-scale model）進行模擬。

大渦的運動基本方程：

在LES方法中通過使用濾波函數(shù)，對于任意流場變量φ可以劃分為兩個部分，即濾波后的大尺度平均分量和通過亞格子模型表示的小尺度分量φ′，濾波后的變量可以通過下式得到：

其中：V表示流動區(qū)域，x為濾波后大尺度空間上的空間坐標，y為實際流動區(qū)域中的空間坐標，G(x,y ）為濾波器函數(shù)且滿足：

假設(shè)濾波過程和求導過程可以交換，用濾波函數(shù)處理瞬時N-S方程和連續(xù)方程得到：

2.3 空化模型[10]

忽略熱傳輸和非平衡相變效應，空化流動中蒸汽相體積分數(shù)輸運方程為：

本文采用Rayleigh-Plesset氣泡動力學方程來描述液體中氣泡的生長和潰滅過程，上式源項可分別表示為：

式中：pv為氣泡內(nèi)壓力，p 為環(huán)境壓力，αv是蒸汽的體積分數(shù)，rnuc是成核點體積分數(shù)，RB=1×10-6m，rnuc=5×10-4，F(xiàn)e=50，F(xiàn)c=0.01。

2.4 縱向平面剛體運動學方程

選取地面坐標系OXYZ并建立原點位于航行體質(zhì)心的彈體坐標系oxyz來描述超空泡航行體運動規(guī)律，如圖1所示。本文假定航行體為質(zhì)量恒定的剛體，建立在航行體質(zhì)心的縱向平面運動學方程組如下：

圖1 坐標系示意圖Fig.1 The sketch map of reference frame

式中：m為航行體質(zhì)量，u，v為航行體速度在體坐標系x，y方向上的分量，r為航行體系繞z軸的角速度，F(xiàn)x，F(xiàn)y為航行體所受合外力在體坐標系x，y方向上的分量，Iz，Mz分別為航行體相對z軸的轉(zhuǎn)動慣量和合外力力矩，Vx，Vy為航行體在地面坐標系下的速度分量，θ為航行體俯仰角。為了研究巡航段航行體運動規(guī)律，進一步假定航行體推力能夠保證航行體維持x方向速度恒定。以上運動方程經(jīng)過時間離散化后利用CFX內(nèi)部的CEL語言嵌入數(shù)值計算模型。

3 仿真模型及邊界條件設(shè)置

數(shù)值模擬采用兩種超空泡航行體模型，模型布局參考文獻[13]，如圖2所示。航行體主要由頭部圓盤空化器、錐柱段彈身、楔形尾翼和尾噴管組成。尾翼的作用類似于常規(guī)水下航行體的十字形舵控制面，尾翼穿透空泡與流體相互作用而產(chǎn)生用于平衡航行體重力的升力及轉(zhuǎn)動力矩。兩種模型的基本參數(shù)詳見表1。計算域及邊界條件設(shè)置如圖3所示，整個計算域均采用六面體結(jié)構(gòu)化網(wǎng)格離散，網(wǎng)格單元總數(shù)為1 121 500，如圖4所示。

圖2 計算模型Fig.2 Computation models

圖3 計算域及邊界條件Fig.3 Computational region and boundary conditions

表1 航行體模型基本參數(shù)Tab.1 Vehicle model parameters

計算域的進口和出口分別定義為速度入口和壓力出口，外邊界采用速度入口條件。本文建立的數(shù)值模型將相對運動原理和非慣性坐標系動力學方程相互結(jié)合。針對航行體縱向平面的兩個平動自由度采用相對運動原理進行處理，平動速度及加速度的變化通過運動學方程計算后施加到相應進口邊界。而對與航行體姿態(tài)的變化則完全采用動量矩方程控制動網(wǎng)格運動求解，通過這樣的設(shè)置可以大大縮減完全使用動網(wǎng)格進行數(shù)值求解的計算量。

圖4 航行體模型計算網(wǎng)格Fig.4 Computational grids around vehicle model

圖5 穩(wěn)態(tài)仿真計算結(jié)果Fig.5 Steady-state simulation results

4 計算結(jié)果與討論

首先，通過調(diào)整流場參數(shù)得到能夠覆蓋整個航行體的超空泡穩(wěn)態(tài)解，將其作為縱向平面彈道模擬的初始解。圖4給出了水深H=5 m、航行體水平速度V0=100 m/s，自然空化數(shù)σv=0.024條件下的穩(wěn)態(tài)初始解。這里自然空化數(shù)定義為：σv=2（ P∞- Pv）/ρV20，其中P∞為環(huán)境壓力、Pv=3 170 Pa為25℃水的飽和蒸汽壓力、ρ為水的密度。

圖6給出了超空泡航行體模型1各項參數(shù)的無控彈道仿真結(jié)果，仿真時間為T=1 s，時間步長 Δt=0.001 s。 圖 6（a），圖 6（b）分別為航行體垂向速度和深度的變化曲線。由于超空泡航行體被氣相介質(zhì)包絡(luò)基本喪失浮力作用，初始階段在重力的作用產(chǎn)生了垂直向下的速度分量，航行體小幅下沉，但持續(xù)時間較短。航行體下沉導致其尾部表面與空泡下表面發(fā)生碰撞，產(chǎn)生沾濕區(qū)，在來流的高速沖擊作用下，沾濕區(qū)域產(chǎn)生“滑行力”的作用。在初始階段尾部滑行力隨著垂向速度的減小迅速增大，在垂向速度達到最小值的同時滑行力也達到峰值，如圖6（d）所示，可見航行體垂向速度對空泡尾部滑行力即沾濕區(qū)域的變化影響顯著。

受尾部滑行力的作用航行體垂向速度方向迅速改變并進一步增加，最后逐漸穩(wěn)定在Vy=4 ms。整個仿真過程中航行體向上爬升約3.5 m。圖6（c）給出航行體俯仰角的變化曲線，可以看到航行體姿態(tài)不斷變化，俯仰角在仿真過程前期持續(xù)增大而后逐漸穩(wěn)定于θ≈3.44°。各項仿真參數(shù)變化在t=0.6 s左右均趨于穩(wěn)定，航行體最后處于穩(wěn)定滑行狀態(tài)如圖7所示。這同文獻[14]所指出的超空泡尾部滑行穩(wěn)定狀態(tài)相同，可見無尾翼的超空泡航行體在無控條件下依然可以保持一定的穩(wěn)定性。

圖8給出了帶尾翼超空泡航行體模型2的無控彈道仿真結(jié)果。可以明顯發(fā)現(xiàn)模型2的仿真結(jié)果與模型1呈現(xiàn)較大差異，尾翼對超空泡航行體的彈道特性產(chǎn)生顯著影響。由于航行體自身重力、尾部控制面升力及俯仰力矩的作用，航行體從一開始便逐漸下沉且垂向速度不斷減小，航行體的姿態(tài)由初始的水平直航狀態(tài)逐漸過渡到以俯仰角θ≈-2.86°的下潛狀態(tài)。模型2仿真的最終結(jié)果如圖9所示，航行體完全位于空泡內(nèi)部中心位置，沒有出現(xiàn)類似于模型1的尾部滑行狀態(tài)。可見由于航行體尾部控制面所產(chǎn)生的升力代替了模型1中的滑行力作用從而構(gòu)成了超空泡航行體比較理想的穩(wěn)定航行方式。

圖6 航行體模1無控彈道仿真結(jié)果Fig.6 The simulation results of uncontrolled trajectory of model 1

圖7 航行體模型1尾部滑行模式Fig.7 Planing model of vehicle model-1

圖8 模型2無控彈道仿真結(jié)果及模型阻力變化曲線Fig.8 The simulation results of uncontrolled trajectory of model 2 and resistance curves

對比圖6和圖8，航行體從初始狀態(tài)到穩(wěn)定狀態(tài)各項參數(shù)變化的歷程可以發(fā)現(xiàn)，具有尾部控制面的航行體能夠更平穩(wěn)、更快地達到穩(wěn)定航行狀態(tài)。所不足的是尾翼的存在大大增加了航行體的流動阻力，從圖8（d）可見模型2的阻力約為模型1阻力的1.52倍，這就意味著需要消耗更多的推進能量來維持帶有尾翼控制面航行體的穩(wěn)定巡航狀態(tài)。

圖9 航行體模型2的穩(wěn)定狀態(tài)Fig.9 Steady state of model 2

5 結(jié) 論

本文耦合計算了多相流體URANS方程和剛體運動學方程，結(jié)合動網(wǎng)格技術(shù)建立了三維超空泡航行體縱向平面無控彈道數(shù)值模型。針對兩種不同的航行體模型進行了無控彈道仿真，分析了超空泡航行體縱向平面的動力學行為，得到的主要結(jié)論如下：

（1）垂直于航行體軸向的速度分量導會致空泡尾部偏移變形進而對尾部滑行力產(chǎn)生顯著影響。

（2）在無控條件下超空泡航行體仍具有一定的運動穩(wěn)定性，本文通過CFD仿真得到了Savchenko所提出的超空泡航行體尾部滑行穩(wěn)定運動方式。

（3）尾部控制面可以有效縮短航行體達到穩(wěn)定所需要的時間，并且能使航行體穩(wěn)定于空泡中心構(gòu)成理想的穩(wěn)定模式。

本文建立的數(shù)值模型為縱向平面超空泡航行體彈道仿真提供了相應的數(shù)值研究手段，據(jù)此可進一步開展大量的數(shù)值模擬研究，系統(tǒng)的研究超空泡航行體彈道變化規(guī)律及各種外在因素對彈道特性的影響。

[1]David R Stinebring,Robert B Cook,John E Dzielski,Robert F Kunz.High-speed supercavitating vehicles[C]//AIAA Guidance,Navigation,and Control Conference and Exhibit.Colorado,August 2006.

[2]Kirschner I N,Kring D C,Stokes A W,Fine N E,Uhlman J S.Control strategies for supercavitating vehicles[J].Journal of Vibration and Control,2002,8(2):219-242.

[3]Reichardt H.The laws of cavitation bubbles as axially symmetrical bodies in a flow[R].Ministry of Aircraft Productuin(Great Britian),Reports and Translations,1946(766):322-326.

[4]Vasin A D.The principle of independence of the cavity sections expansionas the basis for investigation on cavitation flows[Z].VKI Special Course on Supercavitating Flows,Brussels:RTO-AVT and VKI,2001:RTO-EN-010(8).

[5]Kulkarni S,Pratap R.Studies on dynamics of a supercavitating projectile[J].Applied Mathematical Modeling,2000,24(2):113-129.

[6]Kirschner I N,Ulman J S.Simplified dynamical system analysis of supercavitating high-speed bodies[C]//In Fifth International Symposium on Cavitation(CAV2003).Osaka,Japan,2003.

[7]Geol A.Robust control of supercavitating vehicles in the presence of dynamic and uncertain cavity[D].Florida:University Of Florida,2005.

[8]王京華,魏英杰,于開平,王 聰,黃文虎,呂 瑞.基于空泡記憶效應的水下超空泡航行體建模與控制[J].振動與沖擊,2010,29(8):160-163.

[9]Vanek B,Bokor J,Balas G J,Arndt R E.Longitudinal motion control of a high-speed supercavitation vehicle[J].Journal of Vibration and Control,2007,13(2):159-184.

[10]周景軍,于開平.空化器傾斜角對超空泡流影響的三維數(shù)值仿真研究[J].船舶力學,2011,15(1):74-80.

Zhou Jingjun,Yu Kaiping.Three dimensional numerical simulation on the influence of cavitator inclination angle to supercavity flow[J].Journal of Ship Mechanics,2011,15(1):74-80.

[11]Yu Kaiping,Zhou Jingjun,Min Jingxin,Zhang Guang.A contribution to study on the lift of ventilated supercavitating vehicle with low Froude number[J].Journal of Fluids Engineering,2010,132:111303-1～111303-7.

[12]Zhou Jingjun,Yu Kaiping,Zhang Guang.Research on the process of supercavity development and the evaluation on two methods of studying the planing state[J].Journal of Ship Mechanics,2011,15(3):200-206.

[13]Pan Zhancheng,Lu Chuanjing,Chen Ying,Hu Shiliang.Numerical study of periodically forced-pitching of a supercavitating vehicle[C]//9th International Conference on Hydrodynamics.Shanghai,China,2010:856-861.

[14]Savchenko Y N.Control of supercavitating flow and stability of supercavitating motion of bodies[R].RTO AVT lecture series on supercavitating flows,Von Karman Institute,Brussels Belgium,2001:12-16.

CFD analysis of the dynamic behavior of supercavitating vehicle in the longitudinal plane

YU Kai-ping,ZHANG Guang,ZOU Wang,LI Zhen-wang

(School of Astronautics,Harbin Institute of Technology,Harbin 150001,China)

Numerical simulation research on the stability of supercavitating vehicle and trajectory is the basis of designing control and guidance system of supercavitating vehicle.In this paper,based on the multiphase flow URANS(Unsteady Reynolds Averaged Navier-Stokes)equations and the kinematical equation of the rigid body in the longitudinal plane,the three-dimensional dynamics simulation model of supercavitating vehicle is established.On this basis,the dynamic behavior of supercavitating vehicle is simulated in the cruising state.The results show that the supercavitating vehicle without control may still be stable to some extent,and fins can effectively decrease the time which the vehicles take to be stable and make the vehicle locate in supercavity as a stable motion mode.

hydromechanics;supercavitating vehicle;dynamic behavior;CFD

U661.2+2 TV131.2

A

10.3969/j.issn.1007-7294.2014.04.003

1007-7294（2014）04-0370-07

2013-11-10

國家自然科學基金重點項目（10832007）

于開平（1968-），男，哈爾濱工業(yè)大學航天學院教授，博士生導師；

張 廣（1983-），男，哈爾濱工業(yè)大學航天學院博士研究生，Email:zhangguang925@163.com。