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

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

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

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

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

0 引 言

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

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

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

1 計算模型

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

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

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

式中，B為滑行艇艇寬。

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

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

2 數(shù)值計算方法

2.1 計算方法

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

連續(xù)性方程：

動量方程：

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

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

2.2 計算域和邊界條件

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

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

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

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

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

圖3 計算域網(wǎng)格Fig. 3 Grid of computational domain

2.3 航態(tài)估算

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

開始計算時，在不同載荷系數(shù)下設(shè)置的滑行艇初始航態(tài)如圖4 所示。

圖4 滑行艇的航態(tài)估算Fig. 4 Attitude estimation of planing craft

滑行艇的總阻力系數(shù)Ct定義為：

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

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

表2 航態(tài)預(yù)估對阻力預(yù)報的影響Table 2 Influence of trim estimation on resistance prediction

3 計算結(jié)果

3.1 阻力計算結(jié)果

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

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

圖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 壓力分布

滑行艇壓力系數(shù)CP的定義為：

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

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

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

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

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

3.3 自由面波形

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

圖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 水氣分布

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

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

4 結(jié) 論

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

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

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

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

致 謝

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