回轉(zhuǎn)體滑水航行流體動力特性研究

劉富強, 羅凱, 梁紅閣, 黃闖, 耿少航, 董興杰

(1.西北工業(yè)大學(xué) 航海學(xué)院, 陜西 西安 710072; 2.西安慶安電氣控制有限責(zé)任公司, 陜西 西安 710077;3.中國船舶集團有限公司 第705研究所, 陜西 西安 710077)

魚雷具有水下打擊威力大、隱蔽性好等特點，是我國海軍水下攻防的主戰(zhàn)裝備，已成為建設(shè)我國強大海軍的重要力量。然而，傳統(tǒng)魚雷由于水下航行阻力大、平臺空間小及質(zhì)量限制等原因，導(dǎo)致其航程較小，無法有效突破敵方防御圈，發(fā)射后我方平臺易被敵方鎖定，對發(fā)射母艇的威脅較大[1-3]。空中導(dǎo)彈具有航程長、阻力小、殺傷力強等優(yōu)點，然而隨著各國對于空中攔截水平的提升，空中導(dǎo)彈易被敵方鎖定，在未命中目標(biāo)之前將被摧毀[4-5]。而海面低空區(qū)域海況較為復(fù)雜，當(dāng)導(dǎo)彈在近水面掠水航行或者水面滑行時，不易被敵方雷達發(fā)現(xiàn)，極大增強其突防能力，低空近水面區(qū)域的研究現(xiàn)階段已經(jīng)成為了各國研究的重點，因此研究一種水面滑行的攻擊武器非常必要[6-7]。其與水下魚雷相比阻力明顯減小，航程顯著提高，相對于空中導(dǎo)彈，其隱蔽性明顯增強。水面導(dǎo)彈/魚雷武器運動問題可簡化為回轉(zhuǎn)體在水面的滑水問題進行研究。

近年來，國內(nèi)外對于運動體在水面的滑水問題主要集中于對于艦船、滑行艇等的研究。Moctar等[8]對DTC船模在不同速度下的阻力特性進行實驗和數(shù)值仿真研究，數(shù)值計算結(jié)果與試驗結(jié)果吻合性較好；張明霞等[9]基于STAR-CCM+對小水線面三體船數(shù)值仿真，得到不同的水下潛體形狀、側(cè)體位置的阻力變化規(guī)律；Marco等[10]通過數(shù)值仿真階梯形船體在靜水面滑水過程研究船體的沾濕區(qū)域特性以及尾部渦流特性；Amromin[11]通過優(yōu)化船體底部設(shè)計，采用空化的方法為船體減阻提供了解決方案。

國內(nèi)外對小回轉(zhuǎn)體水面滑行的研究較少。Timothy等[12]利用水池進行拖曳實驗，對圓柱棒在多種攻角、不同淹沒深度的滑水問題進行了一系列的實驗研究，獲得流體動力特性的經(jīng)驗解；Waid等[13]對圓柱棒在空化和非空化情況下，不同波形水面滑水工況進行實驗研究，得到空化特性和波面對圓柱棒滑水問題的影響規(guī)律；Serebryakov等[14]對圓柱體沿一定迎角滑水進行了實驗研究，分析了滑水體與型腔表面之間的相互作用過程，研究了深淺度對超空化的影響。

國內(nèi)對于回轉(zhuǎn)體研究主要集中在對于回轉(zhuǎn)體的入水空泡繞流以及水面沖擊特性研究，何春濤等[15]基于VOF(Volume of Fluid)方法求解氣、水兩相流動的RANS方程，并結(jié)合動網(wǎng)格技術(shù)，對回轉(zhuǎn)體垂直入水空泡生成過程，空泡壁面運動特性和空泡表面閉合特性進行了研究；孫凱等[16]采用商用CFD軟件，仿真了平頭回轉(zhuǎn)體模型在亞音速、跨音速以及超音速狀態(tài)下的入水過程，得到了模型在不同速度下的入水阻力特性和頭部直徑、液體可壓縮性及空氣激波對入水過程的影響；張新彬[17]在對仿水黽機器人研究過程中，將機器人水面運動簡化為細長圓柱棒滑水問題，建立了細長圓柱棒與水面接觸相互作用分析模型，研究其運動過程中的力學(xué)特性。

對于大尺度回轉(zhuǎn)體高速滑水問題的實驗研究較為困難，因此提出一種數(shù)值仿真的方法研究回轉(zhuǎn)體滑水問題尤為必要。本文基于STAR-CCM+數(shù)值仿真軟件，選用k-w湍流模型，采用VOF波，構(gòu)建運動體在靜水面滑水?dāng)?shù)值計算模型。對回轉(zhuǎn)體在不同速度、不同淹沒深度滑水工況進行數(shù)值模擬，研究其滑水過程中的流場特性和流體動力特性。

1 數(shù)值仿真模型建立

1.1 數(shù)值方法

本文所有計算均采用西北工業(yè)大學(xué)航海學(xué)院高速航行器創(chuàng)新實驗室的數(shù)值模擬工作站完成，數(shù)值模擬計算軟件版本為STAR-CCM+ 14.02.012。滑水問題屬于氣液兩相問題涉及兩相之間的相互作用, 數(shù)值方法控制方程包括連續(xù)性方程、動量方程和相體積分數(shù)方程[18-19]。在對湍流的非直接數(shù)值模擬中，應(yīng)用最廣泛的是RANS雷諾平均方程，湍流模型采用SSTk-w模型[20]，其基于BSLk-w湍流模型，考慮了湍流剪切力的輸運問題。采用Schnerr & Sauer空化模型[21]模擬回轉(zhuǎn)體水面滑水航行過程中的空化現(xiàn)象。對于滑水問題的求解為了獲得更好的流場特性以及流體動力特性，采用VOF方法模擬氣液兩相界面問題，其可以更好地觀察自由液面的變化[22]。

1.2 VOF波

STAR-CCM+軟件提供VOF波建模，VOF波模型用于模擬輕流體和重流體之間的交界面上的表面重力波。此模型通常與六自由度運動模型一起用于海洋應(yīng)用。STAR-CCM+提供以下 VOF波模型：平波、一階波、五階波、橢圓余弦波、疊加波、不規(guī)則波等。VOF波用于多相模擬，在默認狀態(tài)下，輕流體默認為空氣，重流體默認為水。本文所研究的運動體靜水面滑水問題，采用平波模型構(gòu)建氣液交界面，通過設(shè)置平波水面上的點、水面方向模擬回轉(zhuǎn)體在靜水面滑水航行工況。

1.3 Schnerr & Sauer空化模型

本文采用學(xué)者Schnerr和Sauer[21]在2001年提出 Schnerr & Sauer空化模型，該模型將汽相的體積分數(shù)和單位體積液體所含有的空泡數(shù)量聯(lián)系起來，模型對于相間質(zhì)量傳遞的表達式為：

(1)

(2)

1.4 計算時間步長

本文采用隱式非定常瞬態(tài)求解模型仿真運動體滑水過程[23-24],其計算時間步長的設(shè)置非常關(guān)鍵,設(shè)置時間步長過大,影響仿真結(jié)果的準(zhǔn)確性,設(shè)置時間步長過小,計算較為耗時。因此,一個合適的計算時間步長對于數(shù)值仿真計算同樣重要。

本文研究的運動體滑水問題,計算時間步長選用關(guān)系式[10]如公式(3)所示,最大時間步長為1×10-4s。在本文數(shù)值仿真工況中,運動體速度最大值為21.69 m/s,近壁面網(wǎng)格尺度最小值為3×10-3m,根據(jù)庫朗數(shù)計算公式(4)得最大庫朗數(shù)為0.723,不大于一般系統(tǒng)設(shè)定默認值1.0,選取時間步長合適。

(3)

(4)

式中：V表示運動體的速度；l表示在速度方向上運動體的特征長度；Δt為時間步長；Δx為網(wǎng)格尺度。

2 可行性驗證

2.1 DTC船

Moctar等[8]對DTC船模在6種不同速度下的阻力特性進行實驗,基于本文提出的數(shù)值仿真模型和DTC船模,對6種不同速度下的船模水面航行問題進行數(shù)值模擬。

DTC船模具有嚴(yán)格的對稱性,在滑水計算過程中采用半域計算,減小計算工作量,對稱面附近網(wǎng)格如圖1所示,對其船模附近流域進行多層加密,從而確保仿真計算結(jié)果的準(zhǔn)確性。數(shù)值計算采用k-w湍流模型結(jié)合速度入口實現(xiàn)船模運動,外邊界設(shè)置為速度入口,出口設(shè)置為壓力出口。

圖1 船模對稱面附近網(wǎng)格劃分

本節(jié)對文獻[8]中DTC船模不同速度實驗工況進行數(shù)值仿真,其中DTC船模垂線間長5.976 m,水線寬度0.859 m,靜水面位于龍骨上方0.244 m處,排水量為0.827 m3,速度分別為1.335,1.401,1.469,1.535,1.602,1.668 m/s。

圖2 不同速度模型運動水面示意圖

圖2展示了不同速度下船模穩(wěn)態(tài)航行下的水面示意圖。在DTC船模以不同的速度在水面運動數(shù)值仿真過程中,對其船體受到的阻力進行監(jiān)測,監(jiān)測結(jié)果如表1所示,表中RT表示文獻中試驗得到的總阻力,RF表示文獻中摩擦阻力計算值,RT-S表示數(shù)值仿真過程中船體受到的總阻力,RF-S表示數(shù)值仿真過程中船體受到的摩擦阻力。將不同速度船模數(shù)值仿真模擬結(jié)果的計算值與實驗值對比,總阻力誤差不超過4%,摩擦阻力誤差不超過1.5%,符合工程誤差范圍內(nèi)。

表1 不同速度下船模阻力特性試驗值與仿真值對比

2.2 回轉(zhuǎn)體滑水

Timothy等[12]利用水池進行拖曳實驗,對不同直徑回轉(zhuǎn)體的滑水問題進行了一系列的實驗研究。基于運動體靜水面滑水的數(shù)值仿真模型和文獻[12]中直徑88.9 mm、長度1 422 mm回轉(zhuǎn)體,對回轉(zhuǎn)體以12.19 m/s速度、6°攻角進行一定淹沒深度滑水工況進行數(shù)值仿真。

圖3展示了回轉(zhuǎn)體在沾濕長度為4倍直徑滑水狀態(tài)下的實驗照片和數(shù)值仿真計算結(jié)果。觀察仿真結(jié)果發(fā)現(xiàn)在滑水過程中,運動體尾部飛濺形成大量的水花,與實驗照片吻合性較好。在數(shù)值計算過程中,對回轉(zhuǎn)體滑水的升阻力進行監(jiān)測,該工況下水面滑行實驗和仿真結(jié)果對比如表2所示,數(shù)值模擬升力為19.37 N,阻力為10.21 N。基于公式(5)對升阻力特性進行無量綱化表示,升力系數(shù)為0.032 9,阻力系數(shù)為0.017 3。對比文獻[12]中得到的半域升力系數(shù)0.033 8,阻力系數(shù)0.017 0,升力系數(shù)誤差為-2.66%,阻力系數(shù)誤差為1.76%,升阻力系數(shù)誤差均小于5%,滿足工程誤差要求。

(5)

式中：Fl和Fd分別表示回轉(zhuǎn)體滑水過程中受到的升力和阻力；ρ表示運動環(huán)境介質(zhì)的密度,本文中取水的密度998 kg/m3;v和d分別為回轉(zhuǎn)體的運動速度與直徑。

圖3 滑水實驗照片與數(shù)值仿真對比

表2 回轉(zhuǎn)體特定工況水面滑行實驗仿真結(jié)果對比

本章參考國外運動體在靜水面滑水的文獻,對其實驗工況進行相同工況的數(shù)值仿真計算,將數(shù)值仿真計算結(jié)果與文獻中結(jié)果相對比,發(fā)現(xiàn)數(shù)值模擬流場與實驗照片吻合性較好,且監(jiān)測流體動力與實驗值誤差小于5%,在工程誤差范圍內(nèi)。因此本文所提出的基于STAR-CCM+軟件對運動體水面滑水航行問題數(shù)值模擬計算方法可行,可用于后續(xù)對于回轉(zhuǎn)體在水面滑水航行問題的研究。

3 數(shù)值模擬計算

小運動體低速航行可以依靠實驗進行研究,而對于大尺度、高航速的運動體水面滑行問題實驗較為困難。本章基于前文提出的運動體滑水?dāng)?shù)值計算模型,對回轉(zhuǎn)體在不同工況下的滑水問題進行深入的數(shù)值模擬,研究回轉(zhuǎn)體在不同速度、不同淹沒深度滑水工況下的流場特性和流體動力特性。

3.1 計算物理模型

本章主要對回轉(zhuǎn)體在水面滑水航行不同工況進行數(shù)值仿真計算,其運動體和計算域關(guān)于對稱面對稱,因此在計算過程中采用半域計算,提高計算效率。圖4為計算域尺寸以及邊界條件示意圖,其中回轉(zhuǎn)體直徑為D=88.9 mm,回轉(zhuǎn)體滑水攻角τ為6°,回轉(zhuǎn)體長度為L=1 422 mm,v指示回轉(zhuǎn)體的滑水方向水平向右。考慮滑水運動回轉(zhuǎn)體周圍流場的充分發(fā)展,設(shè)回轉(zhuǎn)體水面沾濕長度Lw為特征長度,回轉(zhuǎn)體迎水面沾濕點距離速度入口流域長度為6Lw,回轉(zhuǎn)體尾部距離壓力出口流域長度為10Lw;高度方向,水面以上空氣流域高度設(shè)定為8D,水面以下水深設(shè)定為12D;流域?qū)挾仍O(shè)定為10D,由于采用半域數(shù)值仿真計算,其實際流域?qū)挾葹?0D。

采用STAR-CCM+軟件進行數(shù)值仿真計算,需要對流域邊界設(shè)定邊界條件,本文研究的回轉(zhuǎn)體滑水航行模型入口采用速度入口,出口采用壓力出口,流域上下界以及寬度方向非對稱面邊界設(shè)定為速度入口,對稱面邊界設(shè)定為對稱平面,運動體即回轉(zhuǎn)體表面設(shè)定為壁面。數(shù)值計算過程中,回轉(zhuǎn)體表面速度設(shè)定為0,通過速度入口給定流域速度實現(xiàn)回轉(zhuǎn)體的相對運動。

圖4 計算域尺寸示意和邊界條件

在利用計算機軟件進行仿真過程中,一般借助網(wǎng)格進行求解。數(shù)值計算模型網(wǎng)格數(shù)量和質(zhì)量嚴(yán)重影響仿真計算的效率和結(jié)果準(zhǔn)確性。在網(wǎng)格邊界層底層網(wǎng)格高度、網(wǎng)格模型、計算域保持不變條件下,對回轉(zhuǎn)體尾部淹沒區(qū)域局部加密得到網(wǎng)格單元總數(shù)目分別為109萬,172萬,198萬,239萬和280萬的5種計算域網(wǎng)格結(jié)果。通過對比回轉(zhuǎn)體12.19 m/s速度、6°攻角特定工況數(shù)值模擬結(jié)果,驗證網(wǎng)格無關(guān)性。

表3為不同網(wǎng)格數(shù)量回轉(zhuǎn)體特定工況水面滑行流體動力特性,與文獻[12]阻力系數(shù)0.017 0、升力系數(shù)0.033 8相比,109萬,172萬和198萬網(wǎng)格數(shù)量模型阻力系數(shù)計算誤差較大,超過20%;而239萬和280萬網(wǎng)格數(shù)量模型流體動力系數(shù)計算誤差明顯較小,均小于3%。考慮網(wǎng)格數(shù)量較大時,增加計算機運行負荷,計算耗時增長。因此選用239萬網(wǎng)格數(shù)量模型用于后續(xù)回轉(zhuǎn)體滑水問題數(shù)值模擬計算。

表3 不同網(wǎng)格數(shù)量回轉(zhuǎn)體特定工況滑水航行流體動力特性

3.2 不同速度滑水

本節(jié)對回轉(zhuǎn)體在不同速度下的滑水工況進行數(shù)值仿真模擬,采用無量綱系數(shù)Cv表征速度見公式(6),采用無量綱系數(shù)λ表征回轉(zhuǎn)體的淹沒深度見公式(7)。回轉(zhuǎn)體滑水特征速度Cv分別取3,8,13,18和23,沾濕特征長度λ取4。

(6)

式中:v為回轉(zhuǎn)體運動速度;d為回轉(zhuǎn)體直徑;g取當(dāng)?shù)刂亓铀俣取?/p>

(7)

式中:Lk為回轉(zhuǎn)體在軸向的沾濕長度;d為回轉(zhuǎn)體的直徑;h為回轉(zhuǎn)體的淹沒深度;τ為回轉(zhuǎn)體的滑水攻角。

圖5 不同速度滑水水面效果圖

本節(jié)通過對不同速度下回轉(zhuǎn)體的滑水工況進行數(shù)值模擬,研究其流場特性和流體動力特性。圖5表征回轉(zhuǎn)體在不同速度下滑水水面效果圖,低速滑水時回轉(zhuǎn)體尾部沾水面積較少,隨著滑水速度的提升,回轉(zhuǎn)體尾部濺起大量水花,回轉(zhuǎn)體尾部發(fā)生大面積沾水,當(dāng)Cv大于13時,明顯觀察到尾部水花飛濺。圖6為回轉(zhuǎn)體在Cv為8,13,23時對稱面速度云圖,回轉(zhuǎn)體在滑水過程中,回轉(zhuǎn)體頭部出現(xiàn)駐點;隨即在周圍展開出現(xiàn)繞流,流域速度最大值為回轉(zhuǎn)體滑水速度的1.2倍左右;回轉(zhuǎn)體尾部在滑水初期局部壓力低于水的飽和蒸氣壓時,發(fā)生空化,同時當(dāng)尾部附近液面壓力高于尾部局部壓力時,液面發(fā)生卷曲,有波浪興起,發(fā)生飛濺。

圖6 不同速度滑水速度云圖

回轉(zhuǎn)體在不同速度下滑水時的流體動力特性如表4所示,阻力和升力采用系數(shù)化表示,其中下標(biāo)f表示摩擦作用下流體動力,p表示壓力作用下的流體動力。圖7為回轉(zhuǎn)體水面滑水過程中受到的壓力和摩擦力在垂直和水平方向的力分解示意圖,水平方向合力為滑水阻力,垂直方向合力為滑水升力。對照圖7分析表4發(fā)現(xiàn),當(dāng)Cv≥8時,滑水阻力中摩擦阻力占總阻力的2/3;滑水升力中壓差升力占比超過97%;而Cv為3時,滑水阻力中壓差阻力明顯高于摩擦阻力,升阻力系數(shù)明顯高于其他速度工況。

圖7 回轉(zhuǎn)體水面滑水力分解示意圖

表4 回轉(zhuǎn)體不同速度滑水流體動力特性

觀察圖5a),Cv=3回轉(zhuǎn)體滑水水面效果圖發(fā)現(xiàn),在回轉(zhuǎn)體以較小速度滑水時,滑水面周圍液滴飛濺較少,高于靜水面區(qū)域回轉(zhuǎn)體幾乎不發(fā)生沾水,相比高速滑水,回轉(zhuǎn)體受到的摩擦力較小;回轉(zhuǎn)體在低速滑水過程中,水下完全沾濕部分提供向上的升力,而在高速滑水過程中由于液滴飛濺在回轉(zhuǎn)體高于靜水面的表面受到向下的壓力產(chǎn)生負升力,因此在回轉(zhuǎn)體低速滑水時升力系數(shù)明顯高于高速滑水時升力系數(shù),尤其壓差作用形成的升力部分。

圖8 不同速度滑水流體動力系數(shù)曲線

圖8為回轉(zhuǎn)體在不同速度滑水時的流體動力系數(shù)曲線。回轉(zhuǎn)體在Cv≥8時升阻力系數(shù)幾乎不變,低速滑水時升力系數(shù)和阻力系數(shù)明顯較高;Cv=3時,阻力系數(shù)高于高速滑水阻力系數(shù)20%,升力系數(shù)是高速滑水升力系數(shù)的3倍,其與回轉(zhuǎn)體低速滑水時流場無明顯液滴飛濺有關(guān);回轉(zhuǎn)體高速滑水時,升阻力特性屬于回轉(zhuǎn)體的固有特性,與滑水速度無關(guān),一般與回轉(zhuǎn)體滑水航行沾濕特性有關(guān)。

3.3 不同淹沒深度滑水

本節(jié)對不同淹沒深度/沾濕長度的回轉(zhuǎn)體滑水工況進行數(shù)值仿真,沾濕特征長度λ取1,2,3和4,回轉(zhuǎn)體滑水速度Cv取13。研究不同淹沒深度回轉(zhuǎn)體滑水過程中的流場特性和流體動力特性。

圖9表征了不同沾濕長度工況下回轉(zhuǎn)體滑水過程水面效果圖,發(fā)現(xiàn)隨著滑水淹沒深度的提高,回轉(zhuǎn)體沾水長度明顯增長,并且在回轉(zhuǎn)體尾部附近水面產(chǎn)生的液滴飛濺效果更加明顯。圖10為λ=4工況回轉(zhuǎn)體滑水液滴飛濺圖,其飛濺長度超過1倍的回轉(zhuǎn)體長度。

圖9 不同沾濕長度滑水水面效果圖

圖10 λ為4工況回轉(zhuǎn)體滑水液滴飛濺圖

觀察圖11不同沾濕長度回轉(zhuǎn)體滑水流域速度云圖發(fā)現(xiàn),在回轉(zhuǎn)體滑水過程中頭部出現(xiàn)駐點,流場速度最大值位于頭部繞流區(qū)域。在4種不同沾濕長度回轉(zhuǎn)體滑水工況中,流場速度最大值幾乎相同,最大值為14.3～14.5 m/s,為流場速度12.19 m/s的1.17～1.20倍,其速度場幾乎不受沾濕長度的影響。

圖11 不同沾濕長度滑水速度云圖

同時在數(shù)值計算仿真的過程中,對回轉(zhuǎn)體滑水的流體動力進行了監(jiān)測,在不同沾濕工況下回轉(zhuǎn)體的升阻力特性如表5所示。回轉(zhuǎn)體在不同沾濕工況滑水流體動力系數(shù)曲線如圖12所示,其升力系數(shù)和阻力系數(shù)均隨著沾濕長度的增加而增大。觀察升阻力系數(shù)曲線的線性特性發(fā)現(xiàn),阻力系數(shù)曲線一次線性特性較強,而升力系數(shù)相比較弱。其主要由于回轉(zhuǎn)體在不同沾濕工況滑水過程中,隨著淹沒深度和沾濕長度的增加,回轉(zhuǎn)體周圍液滴飛濺增強,靜水面以上部分沾水面積增大,隨即受到的壓力作用分布不均勻。而壓差作用力分解主要體現(xiàn)為壓差升力,對于阻力貢獻較小。因此在不同沾濕工況滑水過程中,回轉(zhuǎn)體阻力系數(shù)隨沾濕長度變化線性特性較強,升力系數(shù)隨沾濕長度變化線性特性較弱。

表5 回轉(zhuǎn)體不同沾濕工況滑水流體動力特性

圖12 不同沾濕工況滑水流體動力系數(shù)曲線

4 結(jié) 論

本文基于STAR-CCM+數(shù)值仿真軟件，構(gòu)建運動體在靜水面滑水?dāng)?shù)值計算模型。用數(shù)值仿真的方法對回轉(zhuǎn)體在不同速度和不同淹沒深度工況下的滑水過程進行數(shù)值仿真，得到研究運動體滑水航行流體動力特性的預(yù)報方法，為近水面區(qū)域武器的研究提供理論參考，主要得到以下結(jié)論：

1) 回轉(zhuǎn)體滑水過程中，速度場分布幾乎不受沾濕長度和滑水速度的影響，速度場最大值為滑水速度的1.2倍左右；

2) 回轉(zhuǎn)體在Cv≥8滑水過程中，流體動力系數(shù)幾乎不變，其不受速度大小的影響；而在Cv=3時，流體動力系數(shù)較高，阻力系數(shù)高于高速滑水阻力系數(shù)20%，升力系數(shù)是高速滑水升力系數(shù)的3倍，其與回轉(zhuǎn)體低速滑水時流場無明顯液滴飛濺有關(guān)，造成該現(xiàn)象的根本原因是回轉(zhuǎn)體低速滑水時表面壓力和摩擦力作用分布不同；

3) 回轉(zhuǎn)體不同淹沒深度滑水過程中，其阻力系數(shù)隨淹沒深度增長線性特性較強，升力系數(shù)線性特性較弱，其主要與不同淹沒深度下回轉(zhuǎn)體表面沾濕不均從而導(dǎo)致壓力分布有關(guān)。