孫燁，司先才*，2，馬靖豐，譚俊哲，2，張雪飛，王樹杰，2

1中國海洋大學(xué)工程學(xué)院，山東青島266100

2青島市海洋可再生能源重點實驗室，山東青島266100

0 引 言

在無人帆船實現(xiàn)智能航行的過程中，操縱性對于無人帆船的循跡航行具有一定的影響，是十分重要的航行性能。在船舶操縱性能方面，艉舵發(fā)揮著重要作用，包括轉(zhuǎn)彎能力、初始轉(zhuǎn)彎能力、偏航檢查能力和航向能力。在設(shè)計階段，船—舵系統(tǒng)的水動力分析對準(zhǔn)確預(yù)報船體操縱運(yùn)動具有重要作用。目前，針對船—舵系統(tǒng)操縱性水動力計算的方法主要有經(jīng)驗公式法、系統(tǒng)辨識法、模型試驗法和數(shù)值計算方法，其中基于CFD數(shù)值計算方法預(yù)報船舶操縱水動力的研究較多。Liu與Yasukawa等［1-2］基于 CFD 方法，采用 MMG（Ship manoeuvring mathematical model group）方法建立船舶運(yùn)動數(shù)學(xué)模型，對船—舵輪廓以及船舶操縱性能和KVLCC2模型船操縱性的影響進(jìn)行了研究；Sun等［3］采用CFD方法，對Kriso集裝箱船的船—槳—舵系統(tǒng)進(jìn)行了水動力學(xué)研究，并對模型和全尺寸KCS船的船—槳—舵系統(tǒng)進(jìn)行了瞬態(tài)兩相流計算。Shin等［4］提出了一種用來提高波浪下船舵失速角度的波浪扭曲方向舵，通過CFD方法對其艉舵的水動力性能進(jìn)行了分析，并將結(jié)果與試驗結(jié)果進(jìn)行了比較驗證。許輝和麻紹鈞［5］基于CFD方法計算了不同船速和舵角下的船—舵水動力干擾系數(shù)，并與試驗結(jié)果進(jìn)行了對比驗證分析。馮松波等［6］研究了KVLCC2船在斜航狀態(tài)下船體與艉舵的水動力學(xué)性能，并對CFD方法進(jìn)行了驗證與確認(rèn)。上述研究大多是針對螺旋槳驅(qū)動的船體操縱性水動力進(jìn)行計算分析，對于以風(fēng)帆驅(qū)動的帆船的操縱性研究較少，并且由于無人帆船結(jié)構(gòu)的特殊性，易受到外界環(huán)境的影響而呈現(xiàn)嚴(yán)重的非線性，航向很難控制。因此，準(zhǔn)確預(yù)報無人帆船船體的操縱性能，對安全航行來說就顯得十分重要。

無人帆船船體受風(fēng)帆的作用，會引起不同的漂角，導(dǎo)致船體產(chǎn)生一定的偏航，此即為船體的斜航狀態(tài)。在此狀態(tài)下，操舵可以使船體恢復(fù)原有航向，防止船體出現(xiàn)較大的漂角。為了探討無人帆船的舵力效能對船體操縱性能的影響，本文擬利用Fluent軟件建立船體斜航狀態(tài)下的粘性流場模型。在采用CFD方法模擬船體運(yùn)動時需要使用動網(wǎng)格技術(shù)，會涉及到網(wǎng)格的再生與變形，但當(dāng)船體的偏航角與運(yùn)動幅度過大時，網(wǎng)格易出現(xiàn)負(fù)體積而導(dǎo)致數(shù)值計算發(fā)散或是計算結(jié)果出現(xiàn)錯誤。而動態(tài)重疊網(wǎng)格技術(shù)在處理復(fù)雜曲面離散以及物體的大幅度運(yùn)動方面具有一定優(yōu)勢，處理結(jié)構(gòu)物具有大幅度運(yùn)動的繞流問題時不需要網(wǎng)格再生，可提高動態(tài)網(wǎng)格的處理效率和計算的準(zhǔn)確性［7-9］。因此，本文將結(jié)合重疊網(wǎng)格技術(shù)計算不同漂角與舵角下船體所受到的橫向力和轉(zhuǎn)艏力矩、船舵的反偏航力矩與失速角；同時，利用MMG方法將水動力學(xué)計算結(jié)果代入帆船操縱運(yùn)動模型中，在Matlab軟件中通過編程討論不同舵角與漂角之間的關(guān)系，并利用Simulink模塊綜合分析船體在迎風(fēng)狀態(tài)下進(jìn)行Z字形航行時的操縱性能。

1 無人帆船模型

1.1 模型參數(shù)的選用

本文以自主設(shè)計的無人帆船模型為數(shù)值模擬研究對象。無人帆船船體（附體）與風(fēng)帆的主尺度如表1所示，舵剖面選用NACA 0018翼型。

表1 小型無人帆船模型主尺度Table 1 The main dimensions of small saildrone model

1.2 帆—船—舵系統(tǒng)受力分析

坐標(biāo)系的建立包含大地固定坐標(biāo)系o0-x0y0z0和載體固定坐標(biāo)系o-xyz。形成大地表面固定坐標(biāo)系的x0y0表示船體航行時靜止的水面，z0軸垂直向下。載體固定坐標(biāo)系遵循右手法則，船身軸線為x軸，方向指向船頭為正，y軸與船舷垂直，指向船舷右側(cè)為正，z軸垂直于xy組成的平面。船體重心G在載體固定坐標(biāo)系o-xyz中位于（xG，0，KG-M）點，（x0G，y0G）表示船體重心G在大地固定坐標(biāo)系中的位置。為簡化分析，假設(shè)如下：船體與船舵為剛性體，在低航速下不考慮興波的影響，穩(wěn)心足夠大（橫傾角較?。?，忽略橫搖對船體操縱的耦合作用，并且在此階段只利用xy水平面下船體的三自由度運(yùn)動模型。圖1所示為小型無人帆船船體的受力分析原理圖。圖中：α為船體運(yùn)動時的漂角；ψ為帆船航向角；δ為船體舵角；R為船體受水阻力與橫向力的合力；Ry，Rx分別為由船體及附體引起的橫向力和水阻力；T，N分別為風(fēng)帆產(chǎn)生的推進(jìn)力和橫向力；P為風(fēng)帆橫向力與推力的合力。由于橫傾角的存在，風(fēng)帆產(chǎn)生的推進(jìn)力T與水阻力Rx并沒有作用在同一作用線上，而是產(chǎn)生了一個偏航力矩MT；此外，風(fēng)帆氣動力的作用點與船體水動力的作用點并不在同一鉛垂線上，風(fēng)帆產(chǎn)生的橫向力N與船體以及附體引起的水動力橫向力Ry會產(chǎn)生一個反偏航力矩MN［9］。

圖1 船體受力分析Fig.1 The force analysis of hull

1.3 帆—船—舵運(yùn)動模型

MMG分離建模方法可將作用于船體的外力以及外力矩分別計算并表示出來，其中包含裸船體、船舵與水流，以及風(fēng)與帆相互作用的力及力矩，同時考慮船體航向角與舵角之間的關(guān)系，結(jié)合受力分析，建立無人帆船在靜水中的三自由度數(shù)學(xué)模型為

式中：uG，vG分別為重心處船體前進(jìn)方向速度和橫向方向速度，uG=u，vG=v+xGr，其中u為船中心縱蕩速度，v為船中心橫蕩速度，r為舷搖角速度；Fx，F(xiàn)y，Mz分別為在x，y，z軸附近重心處作用的外力和力矩；m為船舶質(zhì)量；u?G，v?G為對時間的導(dǎo)數(shù)；IzG為船舶在重心附近的慣性矩?？紤]到船體的附加質(zhì)量mx，my以及附加慣性矩JzG，式（1）可表示為

式中：mx，my，mz分別為載體固定坐標(biāo)系下x，y，z軸的附加質(zhì)量；Jz為載體固定坐標(biāo)系下z軸的轉(zhuǎn)動慣量；φ為轉(zhuǎn)艏角度；XH，XF，XR分別為帆船前進(jìn)方向所受到的來自裸船體、風(fēng)帆、舵與海風(fēng)、海流的相互作用力；YH，YF，YR分別為垂直于帆船前進(jìn)方向所受到的來自裸船體、風(fēng)帆、舵的力；NH，NF，NR分別為來自裸船體、風(fēng)帆、舵的舷搖力矩。未知變量u，v，r可通過式（2）求解。式（3）表達(dá)的是船舶重心在大地固定坐標(biāo)系中的位置：

風(fēng)帆的推進(jìn)力XF、橫向力YF以及舷搖力矩NF可表示為

式中：ρ為空氣密度；θ為帆的相對風(fēng)向角；C1，C2分別為風(fēng)帆的升力和阻力系數(shù)；C3，C4分別為風(fēng)帆的最大推力系數(shù)與橫向力系數(shù)；C5為轉(zhuǎn)矩系數(shù)。結(jié)果采用本實驗室項目中風(fēng)帆動力學(xué)的計算值。

以下方程為船—舵橫向力與力矩計算公式。其中，系數(shù)值通過經(jīng)驗公式求解，并與第2節(jié)的CFD數(shù)值計算結(jié)果進(jìn)行對比分析。

式中：tR為船—舵的減額系數(shù)，aH與xH為船—舵的水動力干擾系數(shù)；FN為船—舵的正壓力；xR為舵中心處縱向坐標(biāo)。船體非線性流體動力與力矩的理論計算公式［10］為

2 基于CFD方法的操縱水動力求解

進(jìn)行CFD數(shù)值仿真是為對船—舵系統(tǒng)操縱水動力導(dǎo)數(shù)的求解與恢復(fù)航向效果進(jìn)行判斷，初步確定船—舵的失速角大小以及對船舵的控制參量，以便于后期為船舵的智能控制提供依據(jù)。

2.1 計算方法與數(shù)值求解

應(yīng)用Fluent軟件進(jìn)行流場數(shù)值模擬，船體周圍的三維粘性流場可用雷諾平均的連續(xù)性方程和動量守恒方程描述，選用SST k-ω湍流模型結(jié)合標(biāo)準(zhǔn)壁面函數(shù)模擬邊界層中近壁面附近的流場。壓力速度耦合采用SIMPLEC算法；壓力方程、動量方程和湍流方程均通過二階迎風(fēng)格式離散以保證計算精度。對于兩相流模擬船體自由面繞流問題，水和空氣遵循質(zhì)量守恒，并利用流體體積（VOF）法處理自由表面；船—舵表面采用無滑移壁面條件。如圖2所示，速度入口位于離船艏1.5倍船長處，壓力出口距離船艉部3.5倍船長，底部距離船體水線面1.5倍船長，左、右側(cè)壁距離船側(cè)各1.5倍船長。

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

2.2 重疊網(wǎng)格技術(shù)在三維結(jié)構(gòu)網(wǎng)格中的實現(xiàn)

本文采用重疊網(wǎng)格技術(shù)，分別對流場、船體和船舵這3套結(jié)構(gòu)網(wǎng)格進(jìn)行劃分，流體域作為背景網(wǎng)格，船體和船舵作為重疊網(wǎng)格嵌套在背景網(wǎng)格中。在重疊網(wǎng)格技術(shù)中，固定的背景網(wǎng)格與嵌套網(wǎng)格結(jié)合使用。所謂背景網(wǎng)格，即在流場區(qū)域內(nèi)生成的網(wǎng)格，嵌入背景區(qū)域的重疊區(qū)域生成的網(wǎng)格為重疊網(wǎng)格。嵌套網(wǎng)格與背景網(wǎng)格重疊，并允許嵌套網(wǎng)格在背景網(wǎng)格內(nèi)部自由移動和轉(zhuǎn)動。背景與重疊網(wǎng)格之間的信息交換通過傳輸單元（貢獻(xiàn)單元/受體單元），使用線性或距離加權(quán)內(nèi)插算法完成［11-12］，并對自由表面處進(jìn)行單獨的網(wǎng)格劃分與加密（圖3（a）），有效保證計算要求。

2.3 操縱水動力學(xué)求解結(jié)果

本文以實尺度的船模進(jìn)行數(shù)值模擬，船速U=0.5～1.5 m/s，速度間隔為0.2 m/s，對船體的橫向力與轉(zhuǎn)艏力矩進(jìn)行計算。此外，計算舵角δ的范圍為 0°～40°，舵角間隔為 5°；漂角α范圍為-20°～20°，間隔為5°。

在CFD計算結(jié)果收斂的情況下，計算了船舵的升力系數(shù)Cl與阻力系數(shù)Cd，進(jìn)而求得了船—舵失速角（圖4）。隨后，對不同漂角和舵角下船舶所受橫向力與轉(zhuǎn)艏力矩進(jìn)行求解。為了初步驗證本文數(shù)值方法的可靠性，將CFD計算值與理論計算值進(jìn)行了對比（圖5），結(jié)果顯示吻合度較好。另外，還利用最小二乘法對CFD計算結(jié)果進(jìn)行了曲線擬合（圖6）。船體的位置導(dǎo)數(shù)求解是計算曲線零點斜率，阻力導(dǎo)數(shù)求解是通過對船體進(jìn)行周期性艏搖運(yùn)動，然后將曲線擬合進(jìn)行計算求解（圖7），表2所示為其計算結(jié)果。

圖4 不同舵角下舵的升、阻力系數(shù)Fig.4 The lift and drag coefficients of the rudder under different rudder angles

圖5 不同漂角下船舶所受橫向力系數(shù)和轉(zhuǎn)艏力矩系數(shù)對比（舵角為0°）Fig.5 Comparison between lateral force coefficients and turning torque coefficients of hull under different drift angles（rudder angle of 0 degree）

圖6 最小二乘法擬合結(jié)果圖Fig.6 Fitting results of least square method

圖7 不同風(fēng)向角下船舵角與斜航角的關(guān)系Fig.7 The relationship of the rudder angle and the drift angle under different wind direction angles

表2 船體操縱水動力導(dǎo)數(shù)求解結(jié)果Table 2 Results of hydrodynamic derivatives of ship maneuvering

如圖4所示，當(dāng)舵角為15°時，舵的升力系數(shù)最大，當(dāng)轉(zhuǎn)到20°狀態(tài)時，舵的升力系數(shù)急劇下降。由CFD計算結(jié)果可知，船舵失速角范圍為15°～20°。并且，隨著舵角的不斷增加，舵葉迎流面積增大，舵的阻力系數(shù)呈上升趨勢。

3 操縱性能分析

根據(jù)第1節(jié)對帆船受力的分析以及采用三自由度數(shù)學(xué)模型對帆—船—舵系統(tǒng)進(jìn)行的分析，討論舵角與航向角之間的關(guān)系；通過對船—舵的水動力學(xué)計算結(jié)果進(jìn)行參數(shù)擬合，并整合到帆—船—舵系統(tǒng)的MMG模型中（式2），利用Matlab軟件中的Simulink模塊搭建框圖并對帆—船—舵系統(tǒng)的運(yùn)動進(jìn)行仿真，分析帆船在不同風(fēng)向角下船舵的保持航向效能。由圖7可知，舵角約在10°之后帆—船—舵系統(tǒng)趨于穩(wěn)定，船體以較小的漂角（0°～6°）航行，符合帆船正常的航行姿態(tài)。

為了量化船舵的舵力航向保持能力，進(jìn)行了-10°/10°的 Z字形操縱模擬研究［1］，繪制的舵角δ和航向角ψ的模擬仿真時歷曲線如圖8所示。從模擬結(jié)果中可以看出，船艏轉(zhuǎn)動隨船舵轉(zhuǎn)動的響應(yīng)程度較快。綜上可知，本次船舵的舵效可用于規(guī)定工況下船體的航向控制。

4 結(jié) 語

圖8 航向角隨舵角變化的時歷曲線Fig.8 The time history curves of the variation of the heading angle with the rudder angle

本文以自主設(shè)計的無人帆船的帆—船—舵系統(tǒng)為研究對象，基于Fluent軟件，對小尺度模型下無人帆船的水動力進(jìn)行數(shù)值計算，并與理論計算結(jié)果進(jìn)行了對比驗證，同時，還將計算結(jié)果代入通過受力分析建立的帆—船—舵系統(tǒng)三自由度操縱模型中，對帆船的操縱性進(jìn)行了預(yù)報。研究結(jié)果表明，船—舵在不同工況下其水動力性能有所差異，舵力效能將最終影響到船體操縱的機(jī)動性。由船—舵的水動力性能與帆船操縱模型可知，基于CFD的方法可以應(yīng)用到船體水動力計算中，用于預(yù)報帆船的操縱性運(yùn)動，基于CFD的方法還可以應(yīng)用到帆—船—舵系統(tǒng)在一定漂角、舵角下船體運(yùn)動的水動力計算中。該研究可為無人帆船在結(jié)構(gòu)設(shè)計階段應(yīng)用CFD方法預(yù)報船舶操縱性打下基礎(chǔ)，從而為后期對船舵的智能控制提供依據(jù)。未來，可進(jìn)一步驗證船—舵的水動力學(xué)計算結(jié)果與試驗值，分析預(yù)報無人帆船在復(fù)雜環(huán)境下的操縱性。