基于CFD的排水型三體無人艇阻力優化研究

2023-05-10 04:41:12孫文祺姜兆禎王文龍

艦船科學技術 2023年7期

孫文祺，張奇，姜兆禎，王文龍

(海軍潛艇學院，山東青島 266199)

0 引言

三體船是一種新興的多體船型，兼具單體船和雙體船的優點，具有較優的耐波性和快速性，是當前船舶領域的研究熱點。多體船的片體對水體擾動強，與單體船型相比流場更加復雜，因此數值計算的難度更大。

目前多體船的參考規范和船型資料較少，傳統的經驗公式對多體船難以實現準確的阻力預報，但是研究表明借助計算流體動力學(computational fluid dynamics,CFD)技術可以有效提升預報準確度。張磊等[1]詳細講解了應用CFD技術預報船體阻力與運動響應的流程。陳康等[2]將仿真結果與實驗結果比對，驗證了仿真結果的可信度。謝云平等[3]利用CFD 方法計算某三體船模的阻力，經試驗驗證計算方法的準確性后，對不同側體型線方案進行計算和比對。張明霞等[4]基于試驗數據對比分析了仿真中船體周圍網格加密尺寸、邊界層首層網格節點高度、湍流模型等因素對計算準確度的影響。周廣利等[5]研究了三體船側體排水量占總排水量比例對船體阻力的影響，并通過試驗驗證了數值仿真的準確性。付麗寧等[6]針對某雙體船提出了多種船型優化方案，并使用CFD商業軟件對阻力優化效果進行仿真評估。李志君等[7]從主側體主尺度和主側體布局2個角度出發，對船體阻力進行仿真和優化。金夢顯等[8]針對某三體船設計了多種側體布局方案，并根據阻力仿真計算進行選優，結果表明側體位于船尾對阻力性能有利。

為準確預報項目研發的無人艇阻力，探討中低速三體船的水動力特性，為船體的快速性優化積累經驗，本文利用STAR-CCM+軟件對該無人艇進行流體仿真，考慮側體布置和側體型線2種因素設置不同的計算方案，對比其阻力計算結果并分析流場特征，為后續的船體研究提供參考。

1 數值計算方法

1.1 數學模型

在粘性流動的范疇內，湍流是一種非線性的復雜流體運動，具有不規則脈動、三維有旋等特征。目前湍流的數值模擬方法主要有直接數值模擬(direct numerical simulation,DNS)、雷諾平均法(Reynolds average Navier-Stokes,RANS)和大渦模擬(large eddy simulation,LES)，其中RANS方法可以模擬高雷諾數流動，計算量較小，能對絕大多數工程問題進行求解，因此得到廣泛應用。本文采用RANS方法模擬數值拖曳水池，假設流體為不可壓縮粘性流，平均流動的連續方程和雷諾方程分別如下：

自由液面的處理采用VOF多相流模型，該模型根據不同位置處每個物相的體積分數追蹤氣液交界面的相位置和相分布。某一相i的體積分數定義如下：

式中：Vi為網格單元中相i的體積；V為單元總體積；N為相數。

根據體積分數可以確定某網格單元中流體的狀態：

根據質量守恒原理得出相的連續性方程如下：

1.2 數值拖曳水池

計算模型為穿浪三體船型，采用深V 型舭部和方尾設計，建模比例為1∶1，側體布置于主體兩側中部偏后位置，如圖1所示。以主側體的中縱剖面為基準定義片體的橫向跨距為a，以主側體尾柱間距表示側體縱向位置為b。原船型全長4.75 m，主體寬(圖1 中c)0.52 m，吃水0.36 m，排水量380 kg，側體長2.5 m，寬(圖1中d)0.3 m。船體模型由CAID軟件建立，經ANSA 進行面網格優化后可直接導入STAR-CCM+進行計算域劃分。船體靜水力參數由Maxsurf 軟件進行核算。

圖1 片體布局示意圖Fig.1 Hull layout diagram

計算域為長方形，如圖2所示。計算模型具有對稱性，為減少計算量采用半船計算。L為船總長，上游距船首1L，下游距船尾4L，側壁距對稱面1L，上下邊界分別距水線面1L和2L。上游為入口并設置流體速度分量，下游為壓力出口，側壁和對稱面均設置為對稱邊界條件，船殼設置為無滑移壁面條件。使用STAR-CCM+的Surface Remesher，Automatic Surface Repair，Prism Layer Mesher，Trimmed Cell Mesher 生成面網格和體網格。在船體舷側表面設置面網格加密，對船體附近、自由液面附近、船行波區域設置3層體網格過渡加密，以捕捉船體周圍流場和尾流信息。根據文獻[1，9，10]的結論調整船體舷側邊界層厚度，使數值模擬所關注區域的y+值保持在200萬左右。

圖2 計算域Fig.2 Computational domain

2 計算方案設置及結果分析

2.1 單個片體阻力分析

根據無人艇的設計指標指定水流速度為2.572 m/s，原船型在額定載況下的主體水線長約4.7 m，設計航速下Fr=0.38，傅汝德數較小說明該無人艇屬于中低速船。通常船體阻力可分為摩擦阻力和壓阻力，為了確定該船在設計航速下各阻力成分的占比，針對主體、側體、三體船3種工況進行計算。分別建立主體、側體和三體船的計算模型，并通過仿真軟件監測船體的摩擦阻力、壓阻力和總阻力，繪制阻力隨物理時間變化的曲線如圖3所示，計算結果穩定無發散。

圖3 阻力曲線Fig.3 Drag curves

各計算方案依次編號為A1～A5，匯總的計算結果如表1所示。可以看出三體船組合前后同一片體的摩擦阻力基本不變，組合后片體間興波擾動會導致主體壓阻力減小和側體壓阻力增大，且兩者變化幅度相近。對單個片體的阻力成分進行分析后可以發現，主體的壓阻力約占其總阻力的52.7%，側體的壓阻力約占53.0%，兩片體的壓阻力占比均較大。

表1 各片體的阻力計算結果Tab.1 The resistance calculation results of each hull

2.2 側體縱向位置對阻力的影響

根據已有的三體船阻力研究，側體縱向位置是影響阻力的重要因素[11]。本文針對側體縱向布局設置了多組計算方案，各方案的b值及計算結果如表2所示。b為負值表示側體尾部在主體尾部之后。隨著側體位置從后向前移動，船體阻力先增大后減小，最小阻力值出現在方案B1處，與原船型相比阻力降低約6.7%。從圖4可以看出，隨著側體位置后移船體摩擦阻力曲線基本保持水平，說明側體縱向位置變動對摩擦阻力影響不明顯；壓阻力與總阻力的變化趨勢相同，說明總阻力的變化主要受壓阻力的影響。

表2 不同側體縱向位置的阻力計算結果Tab.2 Resistance calculation results for different longitudinal positions of the side-hull

圖4 b 值-阻力曲線Fig.4 b value-resistancecurve

從圖5可以看出，B1～B3船體周圍流場趨向復雜，液面等高線更密集，波形更復雜，片體間的興波擾動更明顯，同時波峰升高。上述流場變化導致船體興波阻力增大。

圖5 不同b 值的興波高程圖Fig.5 Wave elevation map with different b values

圖6為各方案的艇底壓力分布圖。可以看出從B1～B3隨著側體位置前移，艇尾底部區域的壓力逐漸減小，艇首底部區域的壓力逐漸增大，主體的首尾壓力差增大。B3和B4的艇底壓力分布較為接近。從B4～B6艇尾底部壓力逐漸增大，艇首底部壓力逐漸減小，首尾壓力差減小。為對比艇底壓力的縱向分布，對B1，B3，B6分別沿140水線和240水線提取壓力值，匯總并繪制曲線如圖7所示。其中橫坐標的零點位于尾柱，向船首方向為負。可以看出兩圖的壓力變化規律大致相同，B3比B1的艇尾壓力小，而B6的艇中后部壓力較低，尾端壓力明顯高于B1和B3，艇首壓力略高于前兩者。

圖6 船底壓力縱向分布曲線Fig.6 Longitudinal distribution curves of bottom pressure

圖7 不同側體型寬的興波高程圖Fig.7 Wave elevation map with different side-hull beam

2.3 側體型寬對阻力的影響

原船型的側體寬0.296 m，長寬比約為8.4，半進流角約為8.2°，在保持側體總長和三體船排水量不變的基礎上，減小側體型寬，以增大側體長寬比并減小側體半進流角，改善側體周圍流場以減小壓阻。根據外觀、穩性等方面的設計要求，設置了3種不同的方案，具體側體船型參數如表3所示。

表3 側體船型參數Tab.3 Side-hull form parameters

各方案的計算結果匯總如表4 所示。可以看出b值為0 m 時側體型寬對船體阻力的影響很小，其中C4，C5的壓阻力結果相近，C6的壓阻力較低，但是吃水增大導致摩擦阻力升高，因此總阻力基本不變。b值為0.4 m 時三種船型的阻力差距稍明顯，C3相對于C1阻力降低1.3%，其中壓阻力降低約4.5%，摩擦阻力升高約2.2%。繪制船體興波高程圖，取C1和C3方案的興波高程作對比如圖7所示。可以看到C3的舷側興波減小，等高線變疏，主體尾部興波幅度也有些許降低，“雞尾波”波峰高度由0.0952 m 降至0.0925 m。