榮吉利,王 超,趙 瑞,李海旭,嚴 昊
(1.北京理工大學 宇航學院力學系,北京100081;2.中國船舶工業系統工程研究院,北京100036)
艦載直升機在執行海上補給、搜救、反潛等任務中表現出巨大的優勢,因此越來越多的艦載直升機被裝備于非航空型艦船上。但與陸基直升機不同,由于艦船運動以及艉部流場結構復雜等因素,使直升機起降受到嚴重的影響[1-3]。當氣流流經艦船會形成湍流邊界層,隨后會與飛行甲板前方的機庫壁面發生碰撞,由于幾何形狀的突然變化,會導致流動發生膨脹,此時湍流邊界層在機庫的分離點處分離,形成自由剪切層,在機庫后方的飛行甲板處形成巨大的回流區,該回流區包含了隨時間變化的渦結構,且具有極強的非定常特性[4-6]。
相關學者使用計算流體力學(CFD)對影響艦船艉流場的因素進行了數值模擬研究。黃斌等人[7]通過建立一個基于NS 方程的艦載直升機著艦域流場分析方法,探究了機庫開合狀態對艦船艉部著艦域的影響,研究結果表明打開機庫門更有利于直升機安全著艦;洪偉宏等人[8]以LHA 兩棲攻擊艦作為參考通過CFD 定常計算研究發現,艦船的上層建筑是誘發艦船艦面空氣湍流的因素之一,而島式上層建筑靠前的艦船可以改善后部甲板的流場環境;但國外相關的計算流體力學技術研究表明,時間精確的非定常計算能夠精確捕捉由艦船船體結構產生的非穩態湍流氣流[9-11]。趙瑞等人[12]通過高精度的DES 方法對SFS 與SFS2 等船型展開非定常數值模擬,得到了詳細的船體流場渦脫落和時間精確的流場數據,并從數值上探究了有無上層建筑與機庫門開合狀態對艉流場的影響;劉長猛等人[13]使用FLUENT 的非穩態標準k-ε 模型進行計算,發現機庫的形狀與尺寸會對飛行甲板周圍渦旋的強度和位置產生影響。
由此可見,艦船艉部空氣流場的變化取決于諸多因素,其中類似后臺階類型的艦船艉部會因為機庫外形尺寸的不同對其流場產生一定的影響。本文使用基于一方程SA(Spalart-Allmaras)模型[14]的脫體渦模擬(Detached-Eddy Simulation,DES)[15]對典型護衛艦船型進行數值模擬,研究了不同機庫外形對艉流場的影響,并且根據非定常計算結果推導了預測艉流場回流區范圍的經驗公式。
1.1.1 Spalart-Allmaras 湍流模型
基于一方程SA 模型的DES 方法,其湍流模型求解有關渦粘性變量的偏微分方程為:


1.1.2 DES 模型
DES 方法是將大渦模擬(LES)方法和雷諾平均(RANS)方法相結合的湍流模擬方法,通過對Spalart-Allmaras RANS 湍流模型進行修正,使在離壁面近處的區域采用RANS 方法模擬邊界層內小尺度渦的運動,而在遠離壁面的區域處采用LES 方法模擬,對小尺度渦采用亞格子模型模擬,對大尺度渦進行直接模擬。通過結合LES 方法和RANS 方法的優點,DES 方法可以有效處理邊界層內小尺度渦的脈動,同時也降低了計算資源和計算時間[16]。
在DES 方法中,SA 模型中的長度尺度,即流場中任意點到最近物面的距離d,用d~來代替,其中d~表達為

在當前計算中,Δ=max(Δx,Δy,Δz)為網格中心到相鄰單元中心距離中最大的一個,CDES=0.65。其他細節可以參考文獻[15]。
1.1.3 計算設置
本文通過CFD++軟件,采用有限體積法對控制方程進行離散,并基于具有空間二階精度的多尺度TVD(Total Variation Diminishing)限制器。時間推進采用具有二階精度的雙時間步長(Runge-Kutta)全隱式算法。定常計算采用一方程SA 模型并將其計算的結果作為非定常湍流計算的初始流場。非定常計算3 000 步,其中時間步長Δt*=1.88×10-2s。
為了研究不同機庫外形對艦船艉部流場的影響,并且減少計算資源的消耗,本文使用三種簡化的護衛艦CFD 模型,如圖1所示,分別命名為A 船、B 船與C 船。三種簡化船型在外形和尺寸上均保持一致,僅對艉部機庫的外形進行改動。其中,A 船機庫橫截面為“口”字型,機庫高度6.5 m;B 船機庫橫截面為“凸”字型,最大機庫高度6.5 m,最小機庫高度4.5 m;C 船機庫橫截面為“口”字型,機庫高度較低,為4.5 m。

圖1 三種簡化船型A 船(左)、B 船(中)與C 船(右)及其尺寸細節Fig.1 Dimensions of three simplified frigate ships
整套網格為圓柱形區域,如圖2所示,船體模型被放置于圓柱形區域的中央,半徑r=4.5ls,高d=0.75ls,其中ls為船長。本文對三種簡化模型采用混合網格進行劃分,該混合網格由兩部分組成:一是整個流場的非結構網格,船體表面利用三角形單元形成的40 層棱柱結構(垂直于壁面的第一層網格高度為5 mm,保證y+=o(10),增長率為1.3)來保證粘性層流動的捕捉;二是艉部飛行甲板上方區域的結構網格,根據Spalart 的“重點區域”理論[17],為了能充分地發揮DES 方法的優點,網格單元應該盡可能地小并且各向同性,故在艉流場重點關注區域處使用結構網格能更好地控制網格的精細程度。網格細節圖如圖3所示。A 船、B 船與C 船的網格量分別為1 278 萬、1 240 萬和1 217 萬。

圖3 網格細節Fig.3 Details of simplified frigate ships’ grid
為了準確地研究艦船艉流場特性并驗證CFD 的計算效果,國際項目TTCP[18]得到了簡化的護衛艦船型SFS2,其模型及尺寸如圖4所示。加拿大國家研究委員會對兩個簡化模型都做過一系列的風洞實驗,得到了飛行甲板及前部上層建筑中特定位置的平均速度場及湍流數據[19-20]。
為了驗證本文CFD 方法的有效性,以SFS2 模型作為標準算例,針對0°風向角下的流場進行了計算。圖5給出了艉部甲板上方直線的時均速度分布曲線與風洞實測數據的對比,該直線與機庫等高度,位于甲板y-z 平面內正中間位置,橫坐標為坐標值與甲板寬度d 的比值,縱坐標為速度分量u、v、w與來流速度U∞的比值。從圖5 中可以看出本文的計算結果總體上與風洞實測數據相吻合,這表明本文采用的計算方法是有效的。

圖4 SFS2 模型及尺寸細節Fig.4 The geometry of SFS2 and its dimensions

圖5 0°風向角計算結果對比Fig.5 Headwind mean velocities normalized by U∞
機庫的幾何形狀是影響艉流場回流區的重要因素之一,保證船體整體尺寸外形一致,探究不同機庫縱截面外形對艉流場的影響。

圖6 0°風向角下y/b=0 截面處速度云圖與流線圖Fig.6 Streamlines and velocity magnitude contours for a headwind
圖6 為0°風向角下y/b=0 截面的速度流線圖與云圖。從圖中可以看出,A 船回流區縱向長度約為2.05 倍機庫高度,B 船回流區縱向長度約為1.95 倍機庫高度,C 船回流區縱向長度約為2.27 倍機庫高度?;亓鲄^范圍差距不大,直升機著艦域均在回流區內。但從速度云圖可以明顯看出,B 船飛行甲板上方的下洗速度要明顯小于A 船;在飛行甲板長度相同的情況下,A 船與B 船回流區縱向長度占據了整個甲板的約74%,而C 船約為57%;從速度云圖也可以看出C 船回流區更靠近機庫,直升機著艦域位于回流區邊緣處。

圖7 右舷45°風向角下橫截面合速度梯度云圖Fig.7 Contours of velocity magnitude for the Green 45° case
圖7給出了右舷45°風向角下橫截面x/l=0,25%,50%,75%,100%處的合速度梯度云圖??梢钥吹?,在甲板的右側出現了y 方向的“回流區”,A 船速度梯度的變化程度明顯大于B 船,特別是在近機庫區域處(x/l=25%與x/l=50%),而相同位置處C 船速度梯度的變化程度要小于A 船。從由高到低的速度梯度變化中可以得到機庫上方邊界產生的剪切層位置,同位置處速度梯度的大幅變化會對直升機安全著艦產生一定的影響。

圖8 0°風向角下俯視截面z/h=1.15 處合速度梯度云圖(上:時均結果,下:瞬時結果)Fig.8 Contours of mean (top) and instantaneous (bottom) velocity magnitude for a headwind,plotted at z/h=1.15 above the deck
圖8 與圖9 分別為0°與右舷45°風向角下俯視截面z/h=1.15 處合速度的梯度云圖。由時均結果可以看出:0°風向角下時均流場結構較為整齊劃一,渦結構簡單,流動會在建筑的邊緣產生分離,導致了自由剪切層;B 船機庫兩側較低高度的建筑外形減弱了后方艉流場的速度分離,與A 船相比速度梯度變化較為平緩;45°風向角下,A 船艉部右舷處的速度梯度變化程度明顯大于B 船,且影響范圍更廣。由瞬時結果可以看出:0°風向角下合速度在背風處以及艉部飛行甲板位置處尾跡的不對稱性表明了船體的幾何建筑構型會出現相干的渦脫落特性。0°與45°風向角下,A 船與C 船飛行甲板上空的流動均較為紊亂,且渦區的覆蓋范圍更廣,而該位置一般為直升機在艉部飛行甲板上方高空盤旋的高度,說明等高機庫的外形使得艉流場具有相對紊亂的湍流結構,會對直升機在飛行甲板上方安全盤旋造成不利的影響。


圖9 右舷45°風向角下俯視截面z/h=1.15處合速度梯度云圖(上:時均結果,下:瞬時結果)Fig.9 Contours of mean (top) and instantaneous (bottom) velocity magnitude for Green 45° case,plotted at z/h=1.15 above the deck
圖10 與圖11 分別為0°與右舷45°風向角下艉部飛行甲板上方中心直線的下洗速度分布曲線,該直線的位置選擇基于直升機在甲板上空懸停時的位置。由結果可見,0°風向角下,A 船最大下洗速度為-0.18U∞,B 船為-0.15U∞,C 船為-0.14U∞,最大下洗速度的位置均約與機庫同高。45°風向角時下洗速度有更為明顯的梯度變化,約在0.36 倍機庫高度處達到極值,A 船最大下洗速度為-0.20U∞,而B 船為-0.15U∞,C 船為-0.12U∞。由此可見,在不同風向角下B 船的最大下洗速度均小于A 船,且具有較低的機庫高度的C 船回流區關鍵區域處的下洗速度明顯減弱。

圖10 0°風向角下甲板上方中心直線下洗速度分布曲線Fig.10 The downwash velocity component of the flight deck center for a headwind case

圖11 45°風向角下甲板上方中心直線下洗速度分布曲線Fig.11 The downwash velocity component of the f light deck center for Green 45° case
0°與右舷45°風向角下監測點的瞬時速度時域特性如圖12 與圖13所示,監測點位于飛行甲板正后方(x/l=0.8,y/b=0,z/H=0.6),該位置約為直升機進艦的途經點??梢钥闯?,瞬時速度分量隨著時間發生周期性的脈動,表明了流動的非定常特性。0°風向角下,A 船垂向速度最大振幅達到了-0.4U∞~+0.35U∞,B 船為-0.25U∞~+0.2U∞,C 船為-0.05U∞~+0.1U∞,A 船速度震蕩劇烈且幅值要遠遠大于B 船與C 船,且約為B 船的2 倍、C 船的4 倍;相同的趨勢也存在于45°風向角下,此時A 船垂向速度的最大振幅達到了-0.6U∞~+0.8U∞,B 船為-0.05U∞~+0.6U∞,C 船為-0.05U∞~+0.5U∞。因此,A 船甲板后方的艉流場速度變化極快,會對直升機進艦的操控性帶來嚴重的影響。
圖14 為0°風向角下某瞬時時刻流場的渦量等值面圖。渦量的判別依據為λ2準則,λ2為S2+ Ω2的2階特征值,其中S 為應變率幅值,是速度梯度的對稱張量;Ω 為渦量幅值,是速度梯度的反對稱張量。從圖中可以明顯地看出,因艦艏上層建筑而產生的渦結構強度隨著氣流向艉部流動而逐漸減小,但是遇到艉部機庫后渦結構強度又會加強。A 船與C 船整體的渦結構較為相似,而C 船在甲板處的渦結構強度有所減弱。相比擁有與A 船同高機庫外形的B 船,B 船因其兩側較低的機庫高度使其在甲板處渦結構加強的程度也有所減弱。

圖12 0°風向角下監測點瞬時速度時域特性Fig.12 Time domain characteristics of spectra point for a headwind case

圖13 45°風向角下監測點瞬時速度時域特性Fig.13 Time domain characteristics of spectra point for Green 45° case

圖14 0°風向角下某瞬時時刻流場渦結構Fig.14 Visualization of time-accurate vortex in the airwake (flooded by the magnitude of downwash flow)
Schulman[21]總結并推導了0°風向角下二維后臺階問題中回流區范圍的經驗公式。后臺階的建筑外形會影響流場回流區的整體結構,建筑外形包括其高度(H),迎風寬度(W)以及氣流流經臺面的長度(L)。Wilson[22]給出了評價后臺階的長度尺度R 為

式中:BS為較小的迎風高度H,BL為較大的迎風寬度W。
當L>0.9R 時,流動會在建筑邊緣分離后形成自由剪切層,隨著流動的發展會在某一位置與后方壁面發生再附,此時回流區的最大高度為HR=H。Fackrell[23]給出了再附點與建筑背風面之間的距離,即回流區縱向范圍LR為

用機庫高度H 對其進行無量綱化,得LR,H為

由A 船與C 船的機庫外形尺寸:
HA=6.5 m,HC=4.5 m,WA=WC=16 m,LA=LC=18 m
可得A 船與C 船的長度尺度RA=8.76 m 和RC=6.88 m,由
LA=18 m>0.9RA=7.88 m,LC=18 m>0.9Rc=6.19 m
可知當氣流流經A 船與C 船的機庫后會與甲板產生再附,其位置為

這與前文所述結論完全吻合。
Schulman[21]給出了二維后臺階問題中回流區垂向范圍變化公式:

對于艦船艉流場回流區問題,由公式(3)可以確定不同機庫外形的長度尺度,即機庫外形的幾何因子R。本文根據艦船特殊的建筑構型與公式(5),并對非定常數值仿真計算結果進行擬合,得到了確定等高機庫艉部不同橫向位置處回流區縱向長度的經驗公式LR* 和不同縱向位置處回流區垂向高度的經驗公式zH:

式中:下標H 是用機庫高度無量綱化的值。對于公式(8),當y<0 時取其絕對值。

圖15 A 船經驗公式驗證Fig.15 Verification of empirical equation for ship A
通過經驗公式LR*和zH可以對等高機庫外形的艉流場回流區范圍進行預估。圖15 與圖16 分別是對A 船與C 船艉流場回流區范圍經驗公式的驗證,由圖可以看出本文所推導的經驗公式能夠準確預測回流區的范圍。

圖16 C 船經驗公式驗證Fig.16 Verification of empirical equation for ship C
本文采用基于SA 湍流模型的DES 方法對艦船艉流場進行了非定常數值模擬,詳細分析了0°與右舷45°風向角下不同機庫外形尺寸對艉部流場的影響,得到的主要結論如下:
(1)雖然具有等高機庫外形的A 船與兩側具有較低機庫高度外形的B 船在艉部產生的回流區范圍大體一致,但是和B 船相比,A 船更不利于直升機在艉部飛行甲板上方安全起降與懸停,具體表現為:A 船艉流場速度梯度變化程度大,飛行甲板上空流動更為紊亂,渦區的覆蓋范圍更廣且機庫后方渦結構強度更強,飛行甲板中心直線上的下洗速度比B 船分別大3%(0°風向角)與5%(右舷45°風向角),飛行甲板艉部上方的垂向速度振幅過大,相當于B 船的2 倍;
(2)降低機庫高度會使回流區位置相對前移,因此通過合理設計機庫高度可以使回流區有效避開直升機的著艦域。較低的機庫高度使得艉部速度梯度變化趨于平緩,且大幅度降低了關鍵位置處的下洗速度和艉部渦結構的強度,同時提升了直升機進艦的穩定性;
(3)通過經驗公式可以快速預測艉流場回流區的范圍與位置,為直升機著艦提供相應的參考信息,增加其著艦的安全性。