航行器水面滑行流場特性數值仿真

劉富強,孫 元,王廣平,王雪峰

(中國船舶集團有限公司 第705 研究所,陜西 西安,710077)

0 引言

魚雷作為我國海軍水下攻防的主戰裝備,因其具有水下隱蔽性好、打擊威力大等特點,已成為海軍的重要打擊力量。傳統魚雷由于水下航行阻力大、質量限制等原因,導致其航程較小,無法有效突破敵方防御[1-2]。空中導彈具有阻力小、航程長等優點,然而隨著世界各國空中攔截水平的提升,導彈極易被敵方雷達鎖定,在未命中目標之前即可被摧毀[3-4]。海面低空區域海況較為復雜,當武器在水面滑行或近水面掠水航行時,不易被敵方雷達發現,極大增強其突防能力,低空近水面區域的武器研究已成為各國研究的熱點,因此研究一種水面滑行的攻擊武器非常必要[5-6]。水面導彈/魚雷武器運動問題可被簡化為航行器水面滑行問題進行研究。

近年來,國內外學者針對艦船、滑行艇等航行器水面滑行問題開展了多項研究。Moctar 等[7]對杜伊斯堡測試船型(Duisburg test case,DTC)在不同速度下的阻力特性進行試驗和數值仿真研究,數值計算結果與試驗結果吻合性較好;De Marco等[8]通過數值仿真階梯形船體在靜水面的滑水過程,研究了船體的沾水區域特性以及尾部渦流特性。此外,國內外學者對小尺寸圓柱棒水面滑行問題進行了相關研究。Timothy等[9]利用水池進行拖曳試驗,對圓柱棒在多種攻角和不同淹沒深度下的滑水問題進行了一系列試驗研究,獲得流體動力的經驗解;張珂等[10]對回轉體在水面滑行問題進行數值仿真,研究了滑行過程中回轉體沾濕區域的流場特性以及空泡特性;張新彬[11]將機器人水面運動簡化為細長圓柱棒滑水問題,建立了細長圓柱棒與水面接觸相互作用分析模型,研究其運動過程中的力學特性;劉富強等[12]通過對回轉體在不同速度、不同淹沒深度下水面自由航行過程進行數值仿真,探究其水面滑行流場特性和流體動力特性。

上述研究主要是針對低速航行器,而對于大尺度航行器高速水面滑行問題的研究較少。對大尺度航行器水面滑行問題進行試驗研究較為困難,因此提出一種數值仿真方法研究航行器水面滑行尤為必要。文中基于STAR-CCM+數值仿真軟件,選用剪切應力傳輸(shear stress transport,SST)k-ω湍流模型,采用流體體積(volume of fluid,VOF)波和多重參考系(multiple reference frame,MRF)模型運動參考系,構建了航行器水面滑行數值仿真模型。對航行器在不同速度水面滑行工況進行數值仿真,研究其滑行過程中的流場特性和流體動力特性。

1 數值仿真模型

1.1 數值方法

水面滑行問題屬于氣液兩相問題,涉及兩相之間的相互作用,在湍流的非直接數值仿真中,應用最廣泛的是時均性質的RANS 方程,控制方程包括連續性方程、動量方程、能量方程和相體積方程[13-14]。湍流模型采用SSTk-ω模型[15],其基于Baselinek-ω(BSL)湍流模型,考慮了湍流剪切力的輸運問題。采用Schnerr &Sauer[16]空化模型仿真航行器在水面滑行過程中的空化現象。采用MRF運動參考系技術[14,17-18]將航行器運動規律以平移速度的方式賦予參考系,然后把相對速度代入控制方程進行計算。對于滑水問題的求解,為了獲得更好的流場特性以及流體動力特性,采用VOF方法仿真氣液兩相界面問題,可以更好的觀察自由液面的變化[19]。

1.2 SST k-ω 湍流模型

SSTk-ω湍流模型基于Baselinek-ω湍流模型,考慮了湍流剪切力的輸運問題,通過對渦粘性方程設置限制器,對逆壓梯度條件下流動分離的發生和發展具有較高精度的預測。Menter[15]認為BSL 模型的缺陷主要是由對渦粘性的預測不當引起的,可通過給渦粘性設置限制器的方法解決,并且由此提出了SSTk-ω湍流模型。

1.3 Schnerr &Sauer 空化模型

Schnerr 和Sauer[16]在2001 年提出了 Schnerr &Sauer 空化模型,該模型將汽相的體積分數和單位體積液體所含有的空泡數量聯系起來,模型對于相間質量傳遞的表達式為

1.4 MRF 運動參考系

MRF 模型是學者Issa 等[18]于1994 年提出的一種對于多區域計算較為簡單的定常計算模型。定義坐標系ox′y′z′固連于繞流物體,并且相對于慣性坐標系oxeyeze(以地面坐標系為例),坐標系ox′y′z′的原點位置矢量為r0,平移速度為vt,角速度為 ω。對于計算域內的任一流體微元,假設其相對于坐標系ox′y′z′原點位置為r,則該點在慣性坐標系oxeyeze中的速度可表示為

式中:v為 流體微元在慣性參考系oxeyeze中的絕對速度;vr為流體微元在運動參考系ox′y′z′的相對速度。

采用運動參考系方法,將繞流物體的運動規律以平移速度或旋轉速度的方式賦予參考系,然后把相對速度代入控制方程進行計算。相對速度形式的動量守恒方程需要添加額外的體積力,主要用于描述旋轉和加速等非慣性運動,包括科氏加速度、向心加速度等。運動參考系中流體質點受到的體積力為

式中:α和a分別為運動參考系相對于慣性參考系的旋轉加速度和平移加速度;等號右邊5 個加數項依次為科氏加速度、向心加速度、角加速度、平移加速度和重力加速度。

1.5 VOF 波

STAR-CCM+軟件可提供VOF 波建模,VOF波模型用于仿真輕流體和重流體之間交界面上的表面重力波。此模型通常與6 自由度運動模型一起用于海洋應用。STAR-CCM+提供以下VOF 波模型:平波、1 階波、5 階波、橢圓余弦波、疊加波和不規則波等。VOF 波用于多相仿真,在默認狀態下,輕流體默認為空氣,重流體默認為水。文中所研究的靜水面滑行問題,采用平波模型構建氣液交界面,通過設置平波水面上的點,在水面方向仿真航行器水面滑行航行工況。

2 可行性驗證

2.1 回轉體水面滑行

Timothy 等[9]利用水池進行拖曳試驗,對直徑為88.9 mm 和165.1 mm 圓柱棒的滑水問題進行了一系列試驗研究。基于運動體靜水面滑行的數值仿真模型和文獻[9]中直徑88.9 mm、長度142 2 mm的圓柱棒模型,對圓柱棒以12.19 m/s 水平速度、6°攻角水面滑行工況進行數值仿真。

圖1 展示了圓柱棒在沾濕長度為4 倍直徑滑水狀態下的試驗照片和數值仿真計算結果。觀察仿真結果發現,在滑水過程中,運動體尾部飛濺產生大量水花,與試驗照片吻合性較好。在數值計算過程中,對圓柱棒滑水的升阻力進行監測,該工況下數值計算升力為19.37 N,阻力為10.21 N。

圖1 滑水試驗照片與數值仿真結果Fig.1 Planing test photo and numerical simulation result

基于式(5)[9]對升阻力特性進行無量綱化表示,升力系數CL為0.032 9,阻力系數CD為0.017 3。對比文獻中得到的升力系數0.033 8,阻力系數0.017,升力系數誤差為?2.5%,阻力系數誤差為2.28%,均小于5%,滿足工程誤差要求。

式中:FL和FD分別表示運動體航行過程中受到的升力和阻力;ρ表示運動環境介質的密度,文中取水的密度ρ=998 kg/m3;v和d分別為運動體的運動速度和直徑。

2.2 航行器水下航行

文中所研究的航行器是以某533 mm 口徑魚雷為原型的去鰭舵簡化模型,下文對該簡化模型在水下20 m/s 零攻角航行工況進行數值仿真。圖2為魚雷簡化模型在水下航行的流場圖,觀察壓力云圖可知,魚雷頭部壓力達到最高0.3 MPa,全域最低壓力為19 182 Pa,明顯高于水的飽和蒸氣壓3 170 Pa,在該工況下不發生空化;由圖2(b)速度云圖可見,魚雷頭部形成駐點,速度為零,流域最高速度為22.99 m/s,為航行器航行速度的1.15 倍。魚雷為回轉體模型,在不考慮重力情況下流場上下對稱分布。

圖2 航行器20 m/s 時水下航行流場Fig.2 Underwater flow field around the vehicle navigating at 20 m/s

在航行器水下航行數值仿真計算過程中,對航行器航行阻力進行監測,航行器以20 m/s 在水下穩態航行時的半域阻力為2 372 N,則全域阻力為4 744 N,根據式(6)計算阻力系數為0.106 3。由文獻[20]查詢得該魚雷X 型鰭舵布局阻力系數為0.141,其中單一鰭舵阻力系數為0.009,則魚雷光體阻力系數約為0.105,數值仿真計算阻力系數與文獻值誤差為1.24%,小于5%,處于工程應用誤差范圍內,該數值計算結果可信。

其中,A為航行器最大橫截面積。

參考經典文獻中運動體在靜水面滑行工況以及某533 mm 口徑魚雷水下航行工況,對其進行相同工況的數值仿真計算,將數值仿真計算結果與文獻中結果相對比,發現數值仿真流場與試驗照片吻合性較好,且監測流體動力數值與試驗值誤差小于5%,在工程誤差范圍內。因此文中所提出的基于STAR-CCM+對航行器水面滑水航行問題的數值仿真計算方法可行,可用于后續對航行器水面滑水航行問題的研究。

3 數值仿真計算

基于前文提出的航行器水面滑行數值計算模型對航行器在不同速度工況下的滑水特性進行數值仿真,研究航行器在不同速度水面滑行的流場特性和流體動力特性,仿真速度分別為10,20,30,40,50,60,70,80 和90 m/s 共9 組工況。

3.1 計算物理模型

文中對航行器在不同速度水面滑行工況進行數值仿真,其運動體與計算域關于航行器縱平面對稱,因此在計算過程中采用半域計算,以提高計算效率。圖3 為計算域尺寸以及邊界條件示意圖,航行器直徑D為533 mm,滑水航行攻角取11.5°,航行器在水面滑行初始狀態,水面恰好淹沒航行器尾端,特征沾濕長度L為2.54 m,v指示航行器水面滑行方向水平向左。

圖3 計算域尺寸及邊界示意圖Fig.3 Diagram of calculation domain size and boundaries

考慮到航行器在水面滑行數值仿真計算過程中流場的充分發展,航行器前端流域長度取2.5L,后端流域長度取5L;高度方向,水平面以上空氣部分和水平面以下均取10D;流域寬度方向設置為10D,由于采用對稱半域計算,實際流域寬度為20D。

采用STAR-CCM+軟件進行數值仿真計算,需要對流域邊界進行設定,文中所構建的航行器水面滑行數值計算模型入口采用速度入口,出口采用壓力出口,流域上下界以及寬度方向非對稱面邊界設定為壓力出口,對稱面邊界設定為對稱平面,運動體即航行器表面設定為壁面。

在利用計算機軟件進行仿真過程中,一般借助網格進行求解。數值計算模型網格的數量和質量會嚴重影響仿真計算的效率和結果準確性。文中采用切割體網格和棱柱層網格對流域網格進行劃分,采用表面重構對網格與結構體銜接處進行優化。同時針對航行器尾端沾水重點區域采用體積控制法進行局部體加密。

在網格模型和計算域保持不變的條件下,采用體形狀內不同網格加密尺度對航行器尾部淹沒區域進行局部加密,得到網格單元總數目分別為237 萬,344 萬,496 萬,581 萬和674 萬的5 種計算域網格結果。通過對航行器90 m/s 速度、11.5°攻角、水面恰好淹沒航行器尾端特定工況下的數值仿真結果進行對比,驗證網格數量無關性。

表1 為不同網格數量航行器特定工況下水面滑行流體動力特性,流體動力系數計算參考式(6),237 萬和344 萬網格數量模型阻力系數和升力系數計算值明顯高于其他網格數量模型計算值;581 萬和674 萬網格數量模型流體動力系數計算值幾乎相同。考慮網格數量較大時,會增加計算機運行負荷,計算耗時增長。因此選用581 萬網格數量模型用于后續航行器水面滑行問題數值仿真計算。圖4 為581 萬網格數量模型航行器尾端網格局部放大圖,航行器尾端網格平均尺寸為5 mm,網格質量較好。

表1 不同網格數量下水面滑行流體動力特性Table 1 Hydrodynamic characteristics of the vehicle planing under different mesh numbers

圖4 航行器尾部網格局部Fig.4 Local mesh at the tail of the vehicle

3.2 水面滑行流場特性

圖5 為航行器不同速度滑行工況水面氣液交界面及空化效果圖,圖中黃色為航行器本體,綠色為航行器在水面滑行時體積分數為0.5 的液面,用于表示氣液交界面,紅色橢圓框內為尾部空化區域。圖6 為航行器在不同速度水面滑行工況下的密度云圖,觀察發現航行器在水面滑行過程中尾部發生沾濕,興起波浪,液體飛濺形成水花,并且隨著航行速度的提高,興起水花長度明顯增長,波浪浪高明顯增高。

圖5 不同速度下航行器滑行氣液交界面及空化效果圖Fig.5 Diagram of gas-liquid interfaces and cavitation effects of the vehicle planning at different speeds

圖6 不同速度下航行器水面滑行密度云圖Fig.6 Density contours of the vehicle planning at different speeds

觀察航行器在滑水速度為10 和20 m/s 時水面效果圖(圖5)發現,航行器尾部完全被液面包圍,幾乎不發生空化,無明顯空泡;當滑行速度高于30 m/s時航行器尾部出現空化泡,并且隨著航行速度的進一步增大,航行器尾部空化區域逐漸增大。對比圖6 水面滑行密度云圖,在10 m/s 水面滑行速度時,航行器尾部幾乎沒有空化發生;在50 m/s 水面滑行速度時,航行器尾部有明顯的空泡發生,且空泡閉合;在90 m/s 水面滑行速度時,航行器尾部空泡與大氣貫通,空泡發生破裂。對比航行器不同速度水面滑行尾部液體飛濺情況發現,隨著航行速度的增加,液滴飛濺越高,從最初的小液滴演變為大的尾流波浪。

圖7 為航行器在不同速度水面滑行工況下對稱面的壓力云圖分布,觀察不同速度下水面滑行最高壓力發現:航行器在水面滑水航行過程中,壓力最高點位于航行器沾水最前端,此處由于水面滑行沖擊作用形成高壓區,航行器滑行速度越高,局部最高壓力越大。由圖可見,航行器水面滑行過程中的低壓區主要存在于2 個位置,其一為航行器尾部空化區域,航行器尾部速度略高于航行速度,此處壓力較低,易發生空化;其二為在航行器貼近尾部尾流區域形成的錐狀低壓區。在航行器水面滑行工況數值仿真計算過程中,水的飽和蒸氣壓為3 170 Pa,航行器10 m/s 工況壓力場最低壓力為99 186 Pa,壓力遠高于水的飽和蒸氣壓,不發生空化;而30 m/s 航行器滑水工況壓力場最低壓力為3 169.6 Pa,航行器尾部開始空化,與水面滑行密度場以及空化效果圖顯示結果一致。

圖7 不同速度下航行器水面滑行壓力云圖Fig.7 Pressure contours of the vehicle planing at different speeds

然而觀察航行器在以速度為40 和50 m/s 工況下空化效果圖發現,航行器貼近尾部尾流區域有錐狀封閉區域,且觀察50 m/s 水面滑行壓力云圖發現在該處壓力明顯小于水的飽和蒸氣壓。圖8 為航行器以50 m/s 水面滑行過程中尾部流線局部放大圖,對照該速度下的水面效果、壓力場及流線圖發現,在航行器尾部出現錐狀低壓區,低壓區內主要包含空化的水蒸氣,其密閉于氣泡內未與大氣相通,泡內出現繞流現象。當航行器水面滑行速度繼續提高,在高速滑行工況中,尾流呈波浪狀,原低壓區與空氣相通,空泡潰滅,低壓區消失。觀察航行器在70 和90 m/s 水面滑行工況的壓力場發現最低壓力低于水的飽和蒸氣壓,并且水面滑行速度越高,低壓區范圍越大,航行器尾部空化越明顯。

圖8 速度50 m/s 時航行器水面滑行尾部流線Fig.8 Streamline diagram of the vehicle planing at 50 m/s

3.3 水面滑行流體動力特性

在對航行器不同速度下水面滑行過程進行數值仿真的過程中,對其穩定狀態下水面滑行流體動力進行監測,監測數據如表2所示。

表2 不同速度下航行器水面滑行流體動力特性Table 2 Hydrodynamic characteristics of the vehicle planing at different speeds

圖9 為航行器在不同速度下水面滑行流體動力系數曲線,對比發現航行器在不同速度下的流體動力特性明顯不同。當航行器水面滑行速度為20,30,40 和50 m/s 時,阻力系數明顯較高,升力系數為負值。其他滑行速度工況下,航行器阻力系數基本相同,升力系數為正值。

圖9 航行器水面滑行流體動力特性曲線Fig.9 Curves of hydrodynamic characteristics of vehicle planing

由圖5～圖7 可見,當滑行速度為20,30,40 和50 m/s 時,航行器尾部由于液滴飛濺,沾濕面積明顯大于其他滑行速度工況。圖10 為航行器在水面滑行速度為10,50 和90 m/s 工況下的沾濕效果圖。觀察發現,航行器在水面滑行速度為50 m/s時,表面沾濕區域明顯大于10 和90 m/s 工況下的表面沾濕面積。圖11 為航行器在水面滑行過程中受到的壓力和摩擦力在垂直和水平方向的力分解示意圖,水平方向合力為水面滑行阻力,垂直方向合力為水面滑行升力。從圖中可以看出,在航行器水面滑行過程中,航行阻力主要由摩擦阻力構成,航行升力主要由壓差升力構成。因此航行器在20～50 m/s 速度區間內水面滑行阻力系數明顯較高。

圖10 航行器水面滑行沾濕效果圖Fig.10 Diagram of wetting effect of vehicle planing

圖11 航行器水面滑行力分解示意圖Fig.11 Diagram of force resolution for vehicle planing

當航行器水面滑行速度為20,30,40 和50 m/s時,升力系數為負值,由圖5 可見,此時航行器尾部空泡逐漸生成。圖12 為航行器在50 m/s 速度下水面滑行過程中距離尾部1 倍直徑處橫截面密度云圖,可以明顯看到在航行器尾部周圍有明顯的空泡存在,空泡幾乎完全包裹著航行器尾部。航行器在水面滑行的過程中,尾端發生空化,形成閉合空泡,空泡內壓力接近水的飽和壓力。而航行器尾端上部處于沾濕狀態或裸露狀態,其表面壓力接近于大氣壓。水的飽和壓力遠小于大氣壓力,約為大氣壓力的1/30。因此,在航行器尾端壓差作用效果為向下的負升力。

圖12 速度50 m/s 時航行器水面滑行截面密度云圖Fig.12 Cross-sectional density contour of the vehicle planning at 50 m/s

圖13 為航行器水面滑行速度為30,40 和50 m/s 時表面壓力云圖,航行器表面最低壓力點位于尾端下部,約為水的飽和壓力,而航行器尾端上方壓力為大氣壓。隨著航行器滑行速度的提高,航行器尾端下部低壓區面積逐漸增大,航行器尾端負升力顯著增強。

圖13 不同速度下航行器表面壓力云圖Fig.13 Surface pressure contours of the vehicle at different speeds

然而隨著航行器水面滑行速度的提高,空泡潰滅,與大氣相通,泡內壓力恢復至大氣壓,航行器尾部負升力消失。由航行器在70,80 和90 m/s 水面滑行工況下的升力系數可知,隨著速度的提高,升力系數略微增大,這與航行器在水面高速航行過程中尾部上端液滴飛濺減少有關。

4 結束語

文中基于STAR-CCM+數值仿真軟件,構建了航行器水面滑行的數值仿真模型,對航行器在不同速度下的水面滑行工況進行仿真,通過研究航行器水面滑行過程中的流場特性和流體動力特性,獲得一種研究航行器水面滑行流場特性的預報方法,主要得到以下結論。

1) 航行器在水面滑行過程中,低速工況下幾乎不發生空化,當速度高于30 m/s 時,在航行器尾端發生空化。此時空泡內壓力低于航行器尾端沾濕面壓力,空泡發生形變,液面向航行器尾部卷曲,形成飛濺。

2) 航行器在低速水面滑行初始空化的過程中,尾端形成封閉空泡,泡內為低壓區,泡內出現繞流,此時航行器的升力為負值,待滑行速度提高,空泡潰滅與大氣連通,升力值明顯提高。

3) 航行器以不同速度水面滑行時的流場明顯不同,導致不同速度下升力系數與阻力系數差異較大,升力甚至出現負值,這主要與不同速度下航行器尾端空化效果不同造成沾濕和表面壓力分布存在差異有關。

文中對航行器水面航行流場特性研究采用簡化航行器定常研究,在后續研究中將會針對帶附體復雜航行器水面動態航行過程進行更深層次研究,包括航行器的跨介質水面航行、多自由度運動特性等,以期服務于航行器水面滑行工程應用。