小型無人帆船Z形操縱性能預(yù)測分析

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

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

2青島市海洋可再生能源重點(diǎn)實(shí)驗(yàn)室，山東青島266100

0 引 言

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

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

1 無人帆船模型

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

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

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

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

MMG分離建模方法可將作用于船體的外力以及外力矩分別計(jì)算并表示出來，其中包含裸船體、船舵與水流，以及風(fēng)與帆相互作用的力及力矩，同時(shí)考慮船體航向角與舵角之間的關(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為對(duì)時(shí)間的導(dǎo)數(shù)；IzG為船舶在重心附近的慣性矩。考慮到船體的附加質(zhì)量mx，my以及附加慣性矩JzG，式（1）可表示為

式中：mx，my，mz分別為載體固定坐標(biāo)系下x，y，z軸的附加質(zhì)量；Jz為載體固定坐標(biāo)系下z軸的轉(zhuǎn)動(dòng)慣量；φ為轉(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可表示為

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

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

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

2 基于CFD方法的操縱水動(dòng)力求解

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

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

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

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

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

本文采用重疊網(wǎng)格技術(shù)，分別對(duì)流場、船體和船舵這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)部自由移動(dòng)和轉(zhuǎn)動(dòng)。背景與重疊網(wǎng)格之間的信息交換通過傳輸單元（貢獻(xiàn)單元/受體單元），使用線性或距離加權(quán)內(nèi)插算法完成［11-12］，并對(duì)自由表面處進(jìn)行單獨(dú)的網(wǎng)格劃分與加密（圖3（a）），有效保證計(jì)算要求。

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

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

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

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

圖5 不同漂角下船舶所受橫向力系數(shù)和轉(zhuǎn)艏力矩系數(shù)對(duì)比（舵角為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òng)力導(dǎo)數(shù)求解結(jié)果Table 2 Results of hydrodynamic derivatives of ship maneuvering

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

3 操縱性能分析

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

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

4 結(jié) 語

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

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