結合Savitsky 方法和重疊網格技術的滑行艇阻力數值計算與分析

郭軍，扈喆，朱子文，陳作鋼，崔連正，李貴斌

1 集美大學 輪機工程學院，福建 廈門 361021

2 上海交通大學 海洋工程國家重點實驗室，上海 200240

3 上海交通大學 船舶海洋與建筑工程學院，上海 200240

4 中國船舶及海洋工程設計研究院 噴水推進技術重點實驗室，上海 200011

0 引 言

滑行艇是依靠流體的動力升力將艇體托出水面而使其處于滑行狀態的一種高速艇，具有航速高、機動性好等特點。憑借這些優越的性能，滑行艇已被應用于導彈艇、魚雷艇、偵察艇、交通艇和游艇等船型，在軍用和民用領域具有廣闊的應用前景。滑行艇在整個航行過程中航態變化大、阻力變化劇烈，如何快速、準確地預報滑行艇的阻力性能，是研究人員普遍關心的問題。

早期的滑行艇阻力性能預報方法主要有半理論半經驗公式方法（如Savitsky 法[1]等）、根據系列滑行艇試驗資料整理得出的圖譜方法，以及滑行艇阻力模型試驗。其中，半經驗半理論方法和圖譜方法的預報精度受限于計算船型與母型船的相似程度，而滑行艇的模型試驗雖然能夠直觀地反映出艇體的受力情況以及流場特性，但成本高、耗時長。隨著計算流體動力學（CFD）計算精度的不斷提高，采用CFD 技術不僅可以預報滑行艇的受力及運動情況，還可以捕捉艇體周圍的流場細節，是研究滑行艇水動力問題的一種經濟、有效的方法。近年來，越來越多的研究人員采用CFD 技術對滑行艇進行數值模擬，并取得了一定的研究成果。Brizzolara 和Serra[2]對棱柱型滑行艇流場的CFD 數值模擬精度進行了研究，結果表明CFD 技術可以用于滑行艇水動力問題的分析與計算。曹洪建[3]基于FLUENT 軟件預報了滑行艇在固定姿態下的實船阻力，結果顯示計算誤差約為10%，認為采用FLUENT 軟件對滑行艇進行阻力預報可行，但對升力的預報則不理想。Ghadimi 等[4]采用RANS 方法對滑行艇的姿態進行了數值模擬，結果顯示RANS 方法可以用于滑行艇的初步設計。馬偉佳等[5]采用混合網格對滑行艇的阻力進行了數值模擬，結果顯示數值計算結果與試驗結果有一定的誤差。蔣一[6]對某超高速三體滑行艇的快速性進行了數值計算，結果顯示在速度較低時阻力的計算值與試驗值吻合較好，但隨著航速的增加，誤差會逐漸增大，在最大航速下誤差約為30%。Lotfi 等[7]基于 CFX 軟件對某斷級滑行艇的水動力性能進行了研究，結果顯示其阻力預報精度約為5%。Frisk 等[8]采用FLUENT 軟件和STAR-CCM+軟件對滑行艇的阻力與姿態進行了數值計算，結果顯示滑行狀態下數值預報的阻力偏低，在半滑行狀態下傾角的誤差達32%。De Marco 等[9]采用重疊網格技術和變形網格技術對滑行艇的水動力進行了分析，發現重疊網格技術表現較好。孫華偉等[10]采用STAR-CCM+軟件探討了船體表面第1 層網格節點高度、船體表面網格尺度以及網格節點分布系數這3 個因素對棱柱型滑行艇阻力計算精度、收斂速度及計算穩定性的影響。邵文勃等[11]采用STAR-CCM+軟件探究了時間步長分析與近壁面網格劃分對滑行艇阻力計算精度的影響，結果顯示時間步長和近壁面網格劃分均會對滑行艇的摩擦阻力與艇底的水氣分布造成影響。魏子凡等[12]采用NUMECA軟件對翼滑艇進行了數值模擬，發現加水翼和防飛濺條能改善滑行艇的阻力性能。丁江明等[13]研究了在RANS 計算中網格類型、網格尺寸等因素對滑行艇艇底水氣分布、興波形狀、艇體姿態以及阻力的影響，并與船模阻力試驗結果進行了比較。易文彬等[14]采用FLUENT，CFX 以及 STAR-CCM+軟件探究了時間步長及網格數對預報結果的影響，發現采用CFX 和STAR-CCM+軟件模擬的滑行艇艇底表面有部分水體積分數約為0.5～0.8，艇底的水體積分數預報偏低，艇底水氣分布異常，并認為較小的時間步長以及較細致的網格分布能夠改善艇底水氣分布的效果。李屺楠等[15]采用重新建模法對滑行艇繞流場進行了數值模擬，結果顯示重新建模方法相比重疊網格方法其計算精度和計算效率更高。王慧等[16]采用重疊網格技術對滑行艇的興波進行了模擬，結果顯示采用CFD 方法能夠準確模擬滑行艇在靜水中的航行姿態。

綜上所述，雖然現在采用CFD 技術可以預報滑行艇的阻力性能，但在半滑行狀態和滑行狀態下，滑行艇艇體與流體之間存在劇烈的相互作用，在數值模擬中常存在艇底水氣分布異常以及噴濺區域較難真實模擬的現象，且當滑行艇處于高速航行狀態時其阻力預報與模型試驗之間仍存在較大的誤差，因此需進一步研究半滑行狀態和滑行狀態下的滑行艇阻力高精度數值預報方法。為此，本文擬基于CFD 方法，首先采用Savitsky方法對滑行艇的航行姿態進行預估，再結合重疊網格技術對滑行艇的阻力性能進行預報，以提供一種準確、有效的滑行艇阻力數值預報方法，并分析滑行艇航態預估對阻力預報的影響；然后，在3 種載荷工況及6 個航速下對滑行艇的阻力性能進行數值預報，并與試驗結果進行對比；最后，分析滑行艇的流場特性。

1 計算模型

本文將以Fridsma[17]于1969 年進行的滑行艇系列模型試驗中的一種模型為研究對象。該艇長1.143 m，寬0.228 6 m，型深0.143 m，排水量7.26 kg，底部斜升角為10°。滑行艇的艇體型線和幾何模型如圖1 所示。沿艇長方向，將滑行艇等分為20 站，從第2 站到20 站之間的型線相同，第0 站到第2 站之間的型線如圖1（a）所示，圖中數值表示型線圖的站位。

圖1 滑行艇模型Fig. 1 Geometry model of planing craft

滑行艇的無量綱速度為體積弗勞德數Fr?：

式中，B為滑行艇艇寬。

滑行艇的不同載荷工況如表1 所示。

表1 滑行艇的載荷工況Table 1 Load conditions of planing craft

2 數值計算方法

2.1 計算方法

對滑行艇水動力性能的數值模擬采用商業軟件STAR-CCM+，不可壓縮流體的湍流基本方程由連續性方程和動量方程組成[18]。

連續性方程：

動量方程：

本文的湍流模型采用SSTk-ω 湍流模型，該模型常被用于船舶水動力性能的數值計算，其表達形式為[18]：

在滑行艇阻力數值計算中，采用流體體積（VOF）[19]方法和高分辨率界面捕捉格式（high resolution interface capturing，HRIC）對自由液面進行捕捉，求解器采用基于分離流的黏性求解器，壓力與速度的耦合運用SIMPLE 法，對流項離散格式為二階迎風格式。采用動態流體?固體相互作用模型（dynamic fluid-body interaction，DFBI）模擬滑行艇的六自由度運動，在計算中具有升沉和縱傾的自由度，滑行艇的釋放時間為0.5 s，緩沖時間為2.5 s。采用隱式不定常格式，國際拖曳水池ITTC 推薦的船舶阻力數值計算時間步長 ?t與垂線間長LPP及航速U的關系為 ?t=(0.005～0.01)LPP/U[20]，本文針對滑行艇阻力計算采用的時間步長為0.005LPP/U。

2.2 計算域和邊界條件

由于模型幾何及其繞流場關于中縱剖面基本對稱，故數值模擬對象為滑行艇的半個繞流場。數值計算的計算域和邊界條件如圖2 所示。計算域由靜止的背景區域和隨艇體運動的重疊區域組成，采用重疊網格法處理滑行艇高速航行時艇體姿態變化較大的問題。計算域的長度、寬度和高度分別為13Lpp，5Lpp和5Lpp，速度入口位于艇前2Lpp處，計算域頂部距離滑行艇1.0Lpp。將計算域入口、側面、頂部和底部設為速度入口邊界，計算域的出口設為壓力出口邊界，艇體表面設為無滑移壁面邊界。在計算域的出口前設置消波阻尼以消除傳至出口處的船行波，為避免波的反射，阻尼長度為。

圖2 計算域和邊界條件Fig. 2 Computational domain and boundary conditions

網格類型為切割體非結構化網格，在滑行艇周圍和自由面等區域進行了網格細化。在滑行艇周圍設置3 個開爾文波系形狀的網格加密區以進行局部加密，在設計水線處設置3 個長方體的網格加密區以對自由面網格進行局部加密，在背景區域與重疊區域之間設置1 個網格加密區以進行網格過渡，以使背景區域和重疊區域的網格尺寸處于相同量級。滑行艇的表面網格尺寸為1.75%LPP，滑行艇邊界層網格使用棱柱形網格進行劃分，共15 層，增長率為1.2，邊界層的厚度δ為：

式中：LW為滑行艇的平均浸濕長度；Re為雷諾數。

計算域的自由面網格和滑行艇表面網格的分布如圖3 所示。背景區域和重疊區域的網格量分別為90 萬和141 萬，最終的計算域網格總數為231 萬。

圖3 計算域網格Fig. 3 Grid of computational domain

2.3 航態估算

滑行艇在排水型狀態下的艇體姿態與靜浮時差別不大，而半滑行狀態和滑行狀態下的艇體姿態則與靜浮時差別較大。由于重疊區域內的網格是隨滑行艇姿態的變化而做整體運動，此時艇體周圍及自由面等區域需要更大的網格加密區，但這樣會極大地增加網格數量，降低計算效率。因此在開始計算時，調整滑行艇的航態以使最終計算的航態接近于初始航態，這樣不僅能降低網格數量，提高計算效率，還能改善計算的收斂性。此外，由于網格的方向與來流方向較為接近，故還能改善艇底的水氣分布異常現象[23]，從而提高計算精度。本文將采用半理論半經驗公式Savitsky方法[24-25]對滑行艇的航態進行估算。

開始計算時，在不同載荷系數下設置的滑行艇初始航態如圖4 所示。

圖4 滑行艇的航態估算Fig. 4 Attitude estimation of planing craft

滑行艇的總阻力系數Ct定義為：

式中：Rt為滑行艇總阻力；Awet為滑行艇的動態濕表面積。

在載荷系數C?= 0.608 及航速U= 4 m/s（Fr?=2.904）的工況下，對滑行艇進行數值計算，分析航態預估對滑行艇阻力預報的影響，計算結果與試驗結果的比較如表2（表中，σ 為滑行艇的深沉，Lwet為浸濕長度）所示。由表可見，與單純采用重疊網格相比，結合Savitsky 經驗公式進行航態預估后，滑行艇航態、阻力和浸濕長度的預報精度均有所提高。

表2 航態預估對阻力預報的影響Table 2 Influence of trim estimation on resistance prediction

3 計算結果

3.1 阻力計算結果

在3 種載荷工況及6 個航速（U= 1，2，3，4，5，6 m/s）下，對滑行艇的阻力性能進行數值計算，并與阻力試驗結果進行比較，以驗證數值方法的準確性及有效性。

滑行艇升沉計算結果與試驗結果的對比如圖5 所示。由圖5 可見，升沉的計算結果稍微偏小，3 種載荷工況下升沉的平均誤差分別為0.12%，0.31%和?5.06%。滑行艇縱傾角計算結果與試驗結果的對比如圖6 所示。由圖6 可見，縱傾角的計算結果稍微偏小，3 種載荷工況下縱傾角的平均誤差分別為?0.62%，?2.11% 和?5.86%。滑行艇總阻力計算結果與試驗結果的對比如圖7所示。由圖7 可見，阻力的數值計算結果偏大，3 種載荷工況下阻力的平均誤差分別為9.03%，4.09%和6.93%。滑行艇平均浸濕長度計算結果與試驗結果的對比如圖8 所示。由圖8 可見，3 種載荷工況下浸濕長度的平均誤差分別為?0.63%，1.23% 和1.61%。整體上，滑行艇的升沉、縱傾角、總阻力和平均浸濕長度的計算結果均與試驗值吻合良好，說明本文所提的數值計算方法有效，數值計算結果準確可靠，可用于半滑行狀態和滑行狀態下滑行艇的阻力預報。

圖5 升沉計算值與試驗值對比Fig. 5 Comparison of heaving between CFD and EFD results

圖6 縱傾計算值與試驗值對比Fig. 6 Comparison of trim between CFD and EFD results

圖8 平均浸濕長度計算值與試驗值對比Fig. 8 Comparison of average wet length between CFD and EFD results

3.2 壓力分布

滑行艇壓力系數CP的定義為：

式中：P為滑行艇表面的絕對壓力；P0為參考點壓力（1 個大氣壓）。

當航速U= 4 m/s 時，在不同的載荷工況下，滑行艇艇體表面的壓力分布規律如圖9 所示，其中圖9（a）為滑行艇龍骨線的壓力分布圖，圖9（b）為滑行艇艇體表面壓力分布云圖。圖中，x為沿船長方向的位置。從中可以看出，龍骨線的壓力系數在與水接觸的位置顯著增加，在駐點處時達到峰值，然后減小并達到穩定，在尾部逐漸減小至負值，至艉封板時達到最小；當C?= 0.304，0.608 時，壓力系數的峰值比較接近，導致滑行艇的縱傾角也比較接近；當C?= 0.912時，壓力系數的峰值增加，壓力中心位置前移，導致滑行艇的縱傾角增加。

圖9 不同載荷工況下的壓力分布Fig. 9 Pressure distribution at different load conditions

當載荷系數C?= 0.912 時，不同航速下滑行艇艇體表面的壓力分布規律如圖10 所示。從中可以看出，當Fr?=1.357 時，滑行艇處于半滑行狀態的初始階段，在噴濺區域其壓力系數是先快速增加至最大，隨后迅速減小至常壓，然后緩慢增加，在艉部又迅速減小至負值；當Fr?= 2.714 時，滑行艇處于半滑行狀態的結束階段，此時龍骨線壓力分布的變化規律與滑行狀態的變化規律基本相同；隨著航速的增加，駐線壓力系數的峰值減小，駐點線逐漸向艉部方向移動，駐點線與中縱剖面的夾角逐漸減小。

圖10 不同航速下的壓力分布Fig. 10 Pressure distribution at different speeds

3.3 自由面波形

當載荷系數C?= 0.912 時，不同航速下滑行艇中縱剖面的自由面波高對比如圖11（圖中，z為自由面的波高位置）所示。滑行艇的自由面波形如圖12 所示，其中左圖為整體波形，右圖為噴濺區域波形。由圖可見，隨著航速的增加，船行波的角度逐漸減小，艏部興波的高度逐漸降低，艇后雞尾狀水丘的高度逐漸增大，且逐漸向后移動，艇后“空穴”的深度也逐漸減小，而“空穴”的長度則逐漸增大；當Fr?= 1.357時，在滑行艇艏部發生了明顯的噴濺現象，水流沖擊滑行艇艏部后向前翻卷，噴濺的浸濕面積較大；當Fr?= 2.714 時，在艇體底部發生了噴濺現象，噴濺區域的浸濕形狀接近于三角形；當Fr?= 4.071 時，噴濺區域的浸濕形狀已完全為三角形。

圖11 不同航速下y = 0 剖面的自由面波高對比Fig. 11 Wave pattern of y = 0 at different speeds

圖12 不同航速下的自由波面分布Fig. 12 Free wave surface distribution at different speeds

3.4 水氣分布

當載荷系數C?= 0.912 時，不同航速下滑行艇表面的水氣分布如圖13 所示。由圖可見，滑行艇艇底與水接觸部分水相的體積分數接近于1.0，說明本文艇底水氣分布預報正常，未出現艇底水氣分布模擬異常的現象，這從另一個角度也說明了本文滑行艇阻力數值模擬的準確性。

圖13 不同航速下的滑行艇表面的水氣分布Fig. 13 Volume fraction of water of planing craft at different speeds

4 結 論

本文在求解RANS 方程的基礎上，首先采用Savitsky 方法對滑行艇的航行姿態進行了預估，然后結合重疊網格技術對滑行艇的阻力性能進行預報，并提供了一種滑行艇阻力高精度數值預報方法，接著，在此基礎上對3 種載荷工況下滑行艇的阻力性能進行數值預報，并與試驗結果進行了對比，最后，分析了滑行艇的流場特性。通過研究，主要得出以下結論：

1） 隨著載荷系數的增加，龍骨線壓力系數的峰值增加，壓力中心位置前移，滑行艇的縱傾角增加。

2） 隨著航速的增加，龍骨線壓力系數的峰值減小，駐點線逐漸后移，駐點線與中縱剖面的夾角減小；隨著航速的增加，滑行艇船行波的角度逐漸減小，艏部興波逐漸減弱，艉部雞尾流的高度逐漸增大，艇后“空穴”的深度逐漸減小，但長度逐漸增大。

3） 本文模擬的艇底水氣分布正常，能夠模擬出滑行艇的噴濺現象，且滑行艇阻力數值計算結果與試驗結果吻合良好，這些均說明本文數值方法準確可靠，可為滑行艇水動力數值研究提供有力的技術支撐。

致 謝

本文得到了中國船舶及海洋工程設計研究院噴水推進技術重點實驗室的支持，在此深表謝意！