基于重疊網格的航行體水動力特性計算

韓守根，楊 巖，劉丙鑫，陶 鋼

(1.南京理工大學 能源與動力工程學院，南京 210094；2.南京航空航天大學 直升機傳動技術重點實驗室，南京 210016)

1 引言

潛射航行體的發射是一個涉及多因素干擾，目前為應對各種復雜環境需求，航行體朝著變深度發射、簡化結構、提高發射可靠性等趨勢發展，其中發射可靠性一直是其重點研究方向之一[1]。以往的研究往往只關注水下過程[2-4]，忽略了出水過程的阻力動態演變。因此出水過程的阻力來源對于水下航行體的減阻設計是非常必要的。利用實驗方法與重疊網格技術結合，揭示了航行體出水過程中水冢現象、浮力變化以及附連水質量變化對于水動力及水動力系數影響的內在規律，為工程設計提供了指導。

2 數學模型

2.1 VOF模型

VOF模型中動量方程為求解整個區域內統一的動量方程，作為結果的速度場由各相共享[5]，即：

ρg+F

(1)

相間界面跟蹤通過求解一項或多項容積比率的連續方程來完成。對第q相方程為：

(2)

2.2 水動力系數

水動力系數是水下航行體結構設計的重要參數，一般而言，水動力系數大多通過水動力測量試驗和理論方法計算相結合獲得，由于航行體水動力系數很多，想要全部辨識是不切實際的，本文選取最常用的無量綱化水動力系數來建立數學模型[6]，有：

Cd=F/(0.5ρV2S)

(3)

式(3)中：ρ為水的密度；V為水深z處航行體最大縱向速度，即為質心速度在縱向方向的分量；S為航行體特征橫截面積；L為航行體軸向長度。

航行體質心的流體力合力為：

(4)

3 實驗方案

實驗前航行體與發射平臺在室內進行吊裝，吊裝過程如圖1所示，吊裝完成后,將航行體和發射平臺移至湖內水下，豎直固定在湖內水下。

圖1 航行體發射平臺室內吊裝工作示意圖

為獲得航行體表面壓力，實驗時在航行體表面布置了多個壓力監測點，測點1為航行體頭部駐點，測點2分布于頭部母線段，測點3分布在柱段，測點4位于尾段。航行體幾何尺寸及測點具體位置如圖2所示。

圖2 航行體幾何尺寸及測點位置

試驗當天無風，湖面基本無浪，航行體從水下15m處垂直起射直至射出水面。其水下彈射出筒過程如圖3所示，射出水面過程如圖4所示。

圖3 航行體水下垂直彈射過程

圖4 航行體出水過程

實驗中航行體從初始體到完全出水等過程時序如圖5所示。航行體尾部離開發射筒時，高壓發射氣體在沖出筒口的過程中將會對流場產生了劇烈擾動，形成筒口效應[7]，這將會使得此時的數據誤差較大、可信度低，數據分析時要避免這一段的筒口效應影響，因此本文數值模擬重點為模擬筒口效應消失到航行體完全出水這一段過程，并對這一過程的流場以及水動力系數進行分析討論。

圖5 出水過程時序圖

4 數值模擬

采用重疊網格技術，利用k-epsilon(2eqn)湍流模型、VOF模型計算當前時刻瞬態流場，通過對航行體表面壓力積分，獲得當前時刻航行體所受水動力，在此基礎上計算航行體的水動力系數等; 利用迭代獲得下一時刻的航行體邊界和網格，進而求解下一時刻的瞬態流場，接著計算下一時刻航行體的加速度、速度和位移等得到航行體運動特性。上述過程不斷循環，直至航行體完全出水為止。

為簡化模型和計算，做如下假設：

1)流體常物性，不可壓縮；

2)流體在航行體運動之前是靜止的；

3)航行器是剛體，其外形關于xOz、yOz平面對稱；

4)不考慮其偏航運動；

5)不考慮空泡的影響。

4.1 流體控制方程

流體力學的控制方程為3個基本的守恒方程，即連續性方程(質量守恒方程)、動量守恒方程(運動方程)和能量守恒方程，其中能量守恒方程是包含有熱交換的流動問題必須滿足的基本方程，因為本文所做數值計算中涉及到的流體(海水)是不可壓縮流體，故將其忽略，數值計算時只需要考慮質量守恒方程和動量守恒方程即可。

考慮到本文中的流體處于湍流狀態，所以還需遵守附加的湍流輸運方程。Menter[8]和Rattanasiri[9]等人分別將k-ω和k-ε模型應用于近壁面和遠流場區。本文應用的湍流模型為SSTk-ω湍流輸運方程。

非定常N-S方程[10]：

(5)

湍流動能方程[11]：

(6)

湍流動能耗散率方程：

(7)

式(5)～(6)中：μ為流體粘度；σk為湍流動能普朗特數；σω為湍流耗散普朗特數；Gk和Gω分別為k和ω的耗散；Sω為湍流耗散率的源項；Dω為交叉擴散項；μt為湍流黏度。

4.2 幾何模型及邊界和初始條件

將整個計算域劃分為2個區域：一個是數值波浪水池區域(背景區域)；另一個為航行體區域(子區域)。航行體區域設置相對比較簡單，對于包裹體(航行體的重疊區域)的外部邊界設置為重疊網格邊界。由于沒有波浪和海流，因此不需要考慮波浪和海流的進出口，邊界設置相對簡單，波浪水池四周設置為對稱面，底部設置為滑移壁面,頂部設置為壓力出口。邊界條件設定如圖6所示。對于模擬輸入的初始條件，我們采用了實驗測定的0 s時的速度值，即：Vz0=25 m/s。

圖6 邊界條件設定

4.3 重疊網格模型

重疊網格的劃分主要分為2個區域，即背景區域和航行體區域(子區域)。背景區域在自由表面附近劃分很細的網格，用來捕捉航行體穿越自由液面時自由液面的細節變化，在豎直方向劃分出與子區域尺寸一致的切割體網格，用來更好地與子網格裝配。背景計算域網格劃分完成后，對子區域進行網格劃分，并對其近壁面網格細化。子網格與背景網格重疊，子網格通過插值向背景網格傳遞計算結果，加快計算效率和提高計算精度[12]。網格劃分如圖7所示。

圖7 背景區域和航行體區域網格劃分

5 結果與分析

將實驗數據與仿真計算的結果整理，繪制速度曲線如圖8所示，航行體表面壓力曲線圖如圖9所示。從速度曲線可知，航行體在水下速度下降顯著，水下0.88 s內速度實驗值由25 m/s下降至9.5 m/s，計算值由25 m/s下降至9.7 m/s。從壓力曲線分布可知，頭部駐點壓強隨著航行體速度減小而降低；頭錐段壓力，因為繞流作用和環境壓力的下降而呈現出較穩定的下降趨勢，0.88 s后頭部出水，駐點壓強呈現明顯的下降趨勢，出水前達到最大壓力為0.3 MPa。仿真計算結果與實驗吻合較好，驗證了用重疊網格方法處理航行體運動的準確性和可靠性。

圖8 航行體運動速度曲線

圖9 航行體表面壓力曲線

取航行體在水下0.7 s截面壓力云圖如圖10所示，由圖10可以看出，航行體運動對整個流場壓力的分布影響很小，其主要影響集中在航行體周圍，可以看出明顯的高壓和低壓區。其中最明顯的高壓區包括航行體頂部高壓區，達到0.386 MPa，是阻力的主要來源，低壓區主要表現在航行體肩部附近的低壓區和尾端，分別達到0.184 MPa和0.132 MPa。

圖10 0.7 s時航行體表面壓力云圖

航行體出水過程數值仿真如圖11所示，其中t0指的是航行體頭部出水時刻；t1是肩部出水時刻；t2是因為肩部到主段之間承接處出水時刻；t3指的是航行體尾部剛離開自由液面時刻；t4時刻拐點是航行體出水后完全擺脫附連水時刻。由圖11可知，航行體出水產生水冢，還會穿越水面而形成水冢破碎[13]，之后隨著航行體出水的距離增加，水冢還會增長一段時間，在達到最大水冢之后，水面開始沿著航行體表面回落，完全出水后僅有少量附連水附著于航行體表面，附連水引起的粘滯力是航行體完全出水后主要的流體阻力來源[14]。

圖11 航行體出水過程數值仿真

同時，繪制水動力曲線如圖12所示，水動力系數曲線如圖13所示。其中縱軸代表力的大小和方向。由圖12、圖13可以看出，從航行體頭部到完全出水過程中，航行體水動力及系數在t0時刻達到最大值，分別為0.65 kN和0.184。在整個過程中，水動力及力系數主要是沿Z軸方向變化明顯，X方向波動很小，基本可以忽略。在航行體出水過程中，水動力和力系數波動變化拐點主要出現在整個航行體出水過程，水下段水動力和力系數變化相對穩定。但從筒口效應到出水過程中，水動力及系數在0 s最大，這主要是因為水的靜壓所導致的。

圖12 水動力曲線

圖13 水動力系數曲線

結合圖11，分析航行體出水過程水力和力系數變化曲線上各個拐點出現的原因。t0時刻出現的拐點是因為航行體頭部出水導致的，這是航行體出水過程的開始時刻，水動力以及水動力系數達到最大值(從頭部接觸水面到完全出水)；t1時刻出現的拐點是由于肩部出水導致受力發生變化；t2時刻拐點是因為肩部到尾部之間有段承接處幾何形狀改變導致的；t3時刻拐點是由于航行體尾部出水，這是航行體出水過程的結束時刻。此時航行體已經出水，所受浮力變為零，出水時刻速度較低，其空氣阻力又特別小，可以忽略不計，但此時航行體水動力不為零，這是由于航行體出水時其表面所附帶的水沾濕產生的粘滯力導致的，粘滯力是航行體剛出水后主要的流體阻力來源[15]。t4時刻拐點是因為航行體出水后完全擺脫水的影響，而空氣阻力又特別小，可以忽略不計，因此所受到的流體力及力系數趨于零。

6 結論

1)航行體在水下時頂部壓力是航行體阻力的主要來源力。水動力及系數呈現較為穩定的上升趨勢，而速度下降顯著，計算值由25 m/s下降9.7 m/s。

2)出水過程水動力及力系數主要受水冢現象、浮力變化以及附連水質量變化影響而呈現波動并逐漸趨向于零。

3)粘滯力是航行體完全出水后主要的流體阻力來源。