Maxsurf Motions模塊在船舶耐波性上的應(yīng)用與驗證

張大朋，嚴(yán)謹(jǐn)，趙博文，朱克強

(1.廣東海洋大學(xué)船舶與海運學(xué)院，湛江 524088；2.浙江大學(xué)海洋學(xué)院，舟山 316021; 3.寧波大學(xué)海運學(xué)院，寧波 315211)

海浪是作用在海洋運動物體上的最主要的外力[1]。船舶作為海洋中最主要的運動物體，其在海浪中的運動問題，是船舶耐波性與適航性的基礎(chǔ)。準(zhǔn)確預(yù)報和正確了解船舶在風(fēng)浪中的水動力性能，有利于設(shè)計和建造耐波性能優(yōu)良的船舶，同時也有利于在船舶運營過程中采取適當(dāng)措施以規(guī)避傾覆等海難事故的發(fā)生[2]。

隨著高性能計算技術(shù)的迅速發(fā)展，船舶在波浪中的運動理論也越來越成熟[3]。各種理論和計算方法開始應(yīng)用于船舶運動問題的研究，包括切片理論、高速細長體理論、三維頻域理論和三維時域理論、非線性理論等[4-6]。船舶在波浪上運動的理論和計算方法從二維切片近似方法發(fā)展至三維邊界元法、從零航速發(fā)展至有航速、從頻域法發(fā)展至?xí)r域法，取得了一系列重要的進展[7-10]。這些方法都是應(yīng)用勢流理論，并且自由面條件進行了線性化近似。大多數(shù)方法都是應(yīng)用面元法進行的。從頻域計算轉(zhuǎn)向時域計算是船舶運動的線性理論的另一個重要發(fā)展趨勢，這也反映了工程設(shè)計中對一些非線性問題的需求。時域計算方法較頻域法在處理一些帶非線性問題上有明顯的適用性。目前時域方法中一種是應(yīng)用頻域計算的結(jié)果，通過相應(yīng)的變換確定輻射力的時延函數(shù)和波浪干擾力的脈沖響應(yīng)函數(shù)。應(yīng)用這些函數(shù)，通過對運動時歷或來波的卷積計算相應(yīng)的輻射力或干擾力。最后通過對時域微分方程積分解出運動。

計算機的高速發(fā)展為船舶在波浪中的運動響應(yīng)預(yù)報問題提供了一條新的途徑，計算流體動力學(xué)(computational fluid dynamics，CFD)方法[11-13]。應(yīng)用CFD技術(shù)進行船舶與海洋工程結(jié)構(gòu)物流動數(shù)值模擬和水動力性能預(yù)報的優(yōu)越性和應(yīng)用前景日益顯現(xiàn)。耐波性中涉及的甲板上浪、船體砰擊和載液船體的液艙晃蕩等強非線性問題己經(jīng)可以通過將求解非定常RANS方程進行研究，在一定程度上可以解決部分工程實際問題[14]。近年來，基于黏流理論的CFD軟件如FLUENT、STAR-CCM+、OpenFOAM等得到了快速地發(fā)展，已經(jīng)被廣泛應(yīng)用于船舶耐波性及水動力計算的各類問題，為眾多學(xué)者的研究提供了強有力的計算工具。然而，目前來看理論分析和數(shù)值計算仍不能完全替代模型實驗，最后通常還是需要用模型實驗的結(jié)果來驗證理論計算是否正確。此外，在非線性的運動計算分析、結(jié)果預(yù)報中還是有很多問題存在，需要進一步改進的理論部分還有很多。

目前，應(yīng)用勢流理論計算方法快速地預(yù)估某一船舶在指定海況中的運動特征，在設(shè)計方案的比較階段用理論與計算相結(jié)合的手段對其耐波性能進行比較和選擇，已經(jīng)成為船舶設(shè)計流程中的常規(guī)環(huán)節(jié)[15]。Maxsurf Motions是一款基于勢流理論的船體耐波性預(yù)測程序，目前在船體耐波性計算的研究上已得到較多的應(yīng)用，但其計算精度仍需要進一步討論。采用Maxsurf Motions程序?qū)?種國際標(biāo)準(zhǔn)船模進行耐波性計算，驗證Maxsurf Motions中的切片理論和面元法的可靠性，探討Motions程序在船舶耐波性分析中的可行性，以期對Maxsurf Motions的應(yīng)用起到一定的借鑒和指導(dǎo)作用。

1 基本理論

Maxsurf Motions程序有兩種方法用于計算船舶在波浪中的運動響應(yīng)，分別是切片理論和面元法。切片理論用于計算船舶的垂蕩和縱搖響應(yīng)。橫搖響應(yīng)通過線性橫搖阻尼理論計算。面元法是一種一階衍射/輻射水動力分析方法，其中使用了基于面元的恒定邊界元法(boundary element method，BEM)，能夠計算流體動力附加質(zhì)量和阻尼、一階波浪激振力和力矩、平均漂移力和力矩以及船舶濕表面上的壓力。

1.1 切片理論

使用切片理論計算船舶垂蕩和縱搖的耦合運動時，包含以下基本假設(shè)：船長遠大于船寬和吃水，船寬遠小于波長；船體為剛性，舷側(cè)直立；航速適中且保持不變；相比于波幅，船體運動很小且呈線性；水深遠大于波長；船體對入射波沒有影響(Froude-Kriloff假設(shè))。

船舶的垂蕩、縱搖和橫搖運動本質(zhì)上都是振蕩的，這是由于這些運動中均包含由浮力變化產(chǎn)生的恢復(fù)力。船舶在波浪作用下的運動可視為阻尼-彈簧-質(zhì)量系統(tǒng)。

垂蕩和縱搖的耦合運動中，垂蕩方程為

C35η5=F3eiωt

(1)

縱搖方程為

(2)

當(dāng)船舶在波浪中航行時，除了會產(chǎn)生與航行速度和水深相關(guān)的船行波外，還會產(chǎn)生另一種與入射波浪誘導(dǎo)船體運動相關(guān)的波。這兩種波均會消耗船體的能量且互相疊加，因此船舶在波浪中航行時所消耗的能量要高于靜水中，其中額外的能量損失即為所謂的波浪增阻。Maxsurf Motions程序的波浪增阻計算方法分別為：Gerritsma&Beukelman法、Salvesen法和Havelock法。

Gerritsma&Beukelman法所描述的波浪增阻RAW與船舶的相對垂直速度有關(guān)，計算公式為

(3)

Salvesen法中，波浪增阻的計算公式為

(4)

ib33(ξ)]}dξ

(5)

(6)

(7)

Havelock法中，波浪增阻的計算公式為

(8)

式(8)中：ε為運動與相應(yīng)激振力或力矩的相位差；下角標(biāo)3表示垂蕩；下角標(biāo)5表示縱搖運動。

1.2 面元法

面元法將物體表面離散，用一個平面或曲面代替原來的物面(稱為面元)，在該面元上布置流動的奇點如源、渦、偶極子及其組合進行求解。面元法中引入一個旋轉(zhuǎn)中心，如圖1所示。物體遭遇波浪后通過求解繞旋轉(zhuǎn)中心上的簡諧運動，進而求解船體在波浪中的運動問題。

圖1 面元法旋轉(zhuǎn)中心

假定波高和陡度較小，同時流體不可壓縮、無黏、無旋，因此可以使用線性波理論。流場用速度勢的梯度來描述，速度勢由Laplace方程控制，同時應(yīng)滿足適當(dāng)?shù)倪吔鐥l件。對于簡諧運動，速度勢可表示為

Φ(x;t)=Re[Φ(x)e-iωt]

(9)

式(9)中：Re為雷諾數(shù)。

在整個流體運動域內(nèi)，必須滿足Laplace方程：

?2Φ(x)=0,x∈Ω

(10)

自由水面邊界條件為

(11)

物面邊界條件為

(12)

海底邊界條件為

(13)

式中：?為拉普拉斯算子；Ω為實數(shù)集；g為重力加速度；n為海底的外法線方向；Un為速度；h為水深。

在線性理論的框架內(nèi)，速度勢Φ被細分為入射勢ΦI、衍射勢ΦD和輻射勢Φj。

(14)

式(14)中:xj為第j個運動模式的復(fù)振幅。

入射波的速度勢定義為

(15)

式(15)中：h為水深；β為波浪入射角。

假設(shè)浮體是剛性的，并且在平靜的水中處于穩(wěn)定的平衡狀態(tài)。考慮作用在浮體表面上的水動力，平臺運動方程為

(Cij+Kij)]=Fi

(16)

2 船舶耐波性計算與驗證

2.1 Series-60船舶

Series-60是一艘不帶球鼻艏的貨船船型，是國際公認(rèn)的用于驗證的標(biāo)準(zhǔn)船模，有較豐富的實驗數(shù)據(jù)和文獻資料。Series-60船模的方形系數(shù)Cb=0.68，縮尺比λ=3.93。幾何模型如圖2所示。幾何參數(shù)如表1所示。

圖2 Series-60船舶幾何模型

表1 Series-60船舶幾何參數(shù)

驗證工況為迎浪，分別計算兩個航速Fn=0.198和Fn=0.297。將垂蕩和縱搖運動的響應(yīng)幅值算子(response amplitude operators，RAOs)與實驗和SWAN(simulating waves near shore)程序形成對比，繪制在圖3和圖4中。其中，實驗結(jié)果是由Gerritsma等[16]開展的實驗，即G-B實驗，MIT-SWAN是Sclavounos等[17]開發(fā)的一款船舶非線性時域分析程序。

圖3 Fn=0.198的RAOs曲線

圖4 Fn=0.297的RAOs曲線

當(dāng)Fn=0.198時，垂蕩響應(yīng)中，Motions計算的垂蕩幅值響應(yīng)算子(response amplitude operator，RAO)精度不如MIT-SWAN，但成功預(yù)測到了處于共振峰值時的遭遇頻率范圍，Motions對峰值的預(yù)測誤差約15%，SWAN的誤差約2%。縱搖響應(yīng)中，Motions對峰值的預(yù)測誤差約-16%，SWAN的誤差約8%。當(dāng)Fn=0.297時，Motions對縱搖響應(yīng)的計算精度要優(yōu)于垂蕩響應(yīng)。Motions能合理地預(yù)測Series-60船舶的運動響應(yīng)。當(dāng)Series-60船處于垂蕩和縱搖的共振頻率范圍時，此時船體運動最復(fù)雜，Motions的計算誤差較大。當(dāng)遭遇頻率較低或較高時，Motions的計算結(jié)果與實驗和SWAN程序較吻合。Motions的優(yōu)勢在于能夠在短時間內(nèi)(約幾十秒)預(yù)測多個頻率的RAOs，而相同的工況SWAN程序需要超過8 h的時間。

2.2 SA船舶

SA(SA指在特殊區(qū)域special area內(nèi)航行的船舶)船舶是一艘不帶有球鼻艏和方艉的快速貨船，船型瘦長，符合切片理論的假設(shè)。SA船舶的幾何模型如圖5所示。幾何參數(shù)如表2所示。

表2 SA船舶幾何參數(shù)

圖5 SA船舶幾何模型

驗證工況為迎浪，計算4個航速Fn=0.15、0.2、0.25、0.3，對應(yīng)5.8、7.73、9.67、11.6 m/s。驗證對象主要是垂蕩和縱搖運動的RAOs曲線、運動相位滯后以及附加阻力，與Journee[18]在2001年公布的實驗數(shù)據(jù)以及Geritsma[19]的實驗數(shù)據(jù)形成對比。相位滯后是指垂蕩和縱搖運動與相應(yīng)激振力/力矩的相位差。

通過圖6和圖7可以看出，基于切片理論的SA船舶的運動響應(yīng)與相位滯后總體上與實驗基本吻合，但在某些頻率點出現(xiàn)較大的誤差。Fn=0.15時，縱搖運動的RAOs曲線比垂蕩運動更貼合實驗，而縱搖運動的相位滯后卻與實驗差距較大，尤其是當(dāng)船長波長比大于1(也就是波長小于船長)時，實驗的縱搖運動相位滯后與計算值相差一個負(fù)號。同樣的情況發(fā)生在Fn=0.3。由于實驗數(shù)據(jù)不足，實驗結(jié)果存在相當(dāng)大的分散性，且在兩種航速下實驗都沒有很好地捕捉到響應(yīng)的峰值，因此很難說明計算結(jié)果與實驗結(jié)果不吻合的原因是切片理論的精度不準(zhǔn)確。然而大部分計算結(jié)果顯示，當(dāng)波長較大時，SA船舶的垂蕩和縱搖耦合運動響應(yīng)與實驗吻合良好，一定程度上可以證明對于SA船舶的運動計算，Motions程序在大部分工況下具有較強的精度。

L為船長；λ為縮尺比

圖7 Fn=0.3的SA船體運動響應(yīng)

Raw為附加阻力；ρ為水的密度；g為重力；ζ為波高；B為船寬；L為船舶水線長

觀察圖8可知，無論是低頻段還是高頻段，Havelock法的計算結(jié)果均比較符合實驗情況，尤其是在航速較大的情況下，高頻段的計算結(jié)果與實驗數(shù)據(jù)十分貼合。Salvesen法明顯不適用于SA船舶的附加阻力計算，無論是曲線趨勢還是數(shù)據(jù)點，均與實驗相差甚遠。

值得注意的是，無論是哪一種計算方法，都無法做到在整個波頻內(nèi)與實驗完全吻合。這其中有尺度效應(yīng)的影響。SA船舶是實尺度，而實驗是模型尺寸。在拖曳水池中精準(zhǔn)地測量出船舶的附加阻力是十分困難的，且從模型尺度到實尺度的換算存在修正系數(shù)。SA船舶的修正系數(shù)到目前為止仍是未知數(shù)。綜合以上原因，可以認(rèn)為Motions模塊中的附加阻力計算方法具備一定的計算精度，能夠應(yīng)用于船舶早期設(shè)計階段和附加阻力的估算。

2.3 Wigley船舶

Wigley船型滿足切片理論和Michelle興波阻力積分公式的假設(shè)，型值由式(19)可得

(19)

式(19)中：Lpp為垂線間長；B為最大船寬；d為吃水；x、y、z為空間三維點坐標(biāo)。

d/Lpp=0.062 5；B/Lpp=0.1。幾何模型和參數(shù)分別如圖9和表3所示。

表3 Wigley船舶幾何參數(shù)

圖9 Wigley船舶幾何模型

驗證工況為迎浪，計算3個航速Fn=0.2、0.3、0.4，對應(yīng)1.09、1.64、2.19 m/s，驗證對象主要是附加質(zhì)量、阻尼系數(shù)和波浪力。

圖10～圖12為基于Motions模塊得到的不同航速下Wigley船體的附加質(zhì)量系數(shù)隨遭遇頻率的變化，并與Journée[20]的實驗結(jié)果進行對比。Wigley船體在波浪中運動時，會促使周圍流體的速度發(fā)生變化。由于流體存在慣性，流體將阻礙船體運動并反作用于船體上，其作用力稱為附加慣性力，把該附加慣性力假想為在物體上附加了一部分質(zhì)量，這部分質(zhì)量稱為附加質(zhì)量。附加質(zhì)量與物體本身的形狀及運動方向有關(guān)。

由圖10～圖12可知，不同航速下Journee的實驗結(jié)果均勻地分散在Motions曲線附近，說明基于切片理論的Motions模塊在不同航速下計算船舶附加質(zhì)量時均能得到與實驗較為接近的結(jié)果。其中，中高頻段的吻合度較好，部分低頻段的吻合度較差，說明Motions模塊對中高頻的計算精度較高。3種航速下Motions計算的附加質(zhì)量系數(shù)差距不大，說明航速效應(yīng)對附加質(zhì)量的計算影響較低。4個附加質(zhì)量系數(shù)中，A33和在3個航速下的計算數(shù)據(jù)幾乎一致，說明因垂蕩產(chǎn)生的垂蕩附加質(zhì)量不隨航速的變化而發(fā)生明顯變化。A35和A53的結(jié)果相差一個負(fù)號，說明因縱搖產(chǎn)生的垂蕩慣性力和因垂蕩產(chǎn)生的縱搖慣性力大小相同，方向相反。A55是因縱搖產(chǎn)生的縱搖附加質(zhì)量，航速的增大會使縱搖附加質(zhì)量略微增大。4個附加質(zhì)量系數(shù)隨遭遇頻率的變化是一致的：在中低頻段，附加質(zhì)量系數(shù)隨遭遇頻率的增大快速減小或增大，在高頻段，附加質(zhì)量系數(shù)曲線幾乎趨于一條直線，說明此時的遭遇頻率基本不會對附加質(zhì)量造成影響。

V為船舶排水體積

圖11 Fn=0.3時Wigley附加質(zhì)量隨遭遇頻率的變化

圖12 Fn=0.4時Wigley附加質(zhì)量隨遭遇頻率的變化

圖13～圖15是基于Motions模塊得到的不同航速下Wigley船體的阻尼系數(shù)隨遭遇頻率的變化，并與Journee的實驗結(jié)果進行對比。在勢流理論框架下已經(jīng)假設(shè)流體無旋、無黏、不可壓，因此這里的阻尼不是粘性阻尼，而是物體在波浪運動中所做的功。由于物體在含有自由表面的流體中運動，勢必會興波，而興起波浪必然會帶走能量，這個能量就是物體運動做的功，體現(xiàn)為物體受到的阻力，大小與物體運動的速度相關(guān)。

圖13 Fn=0.2時Wigley阻尼系數(shù)隨遭遇頻率的變化

圖14 Fn=0.3時Wigley阻尼系數(shù)隨遭遇頻率的變化

圖15 Fn=0.4時Wigley阻尼系數(shù)隨遭遇頻率的變化

相比于附加質(zhì)量，Motions對阻尼系數(shù)的計算精度明顯降低。除了B33之外，其余3個阻尼系數(shù)均與計算值出現(xiàn)了較大的差別。對于B35和B55，中低頻段的計算值偏大，高頻段的計算值偏小，B53的現(xiàn)象與B35相反。這3個系數(shù)均與縱搖模態(tài)相關(guān)。而且也看出，隨著航速的增大，阻尼系數(shù)與實驗值的偏離度增大，航速效應(yīng)對阻尼系數(shù)的影響要大于對附加質(zhì)量的影響，尤其是對縱搖模態(tài)的阻尼系數(shù)的影響最為明顯。

2.4 KVLCC2船舶

KVLCC2是一艘由韓國KRISO設(shè)計的現(xiàn)代油船，具有詳細的試驗數(shù)據(jù)。KVLCC2船模帶有球鼻艏和U形船尾框架線，屬于肥大型船。通過KVLCC2船舶驗證面元法對斜浪中的運動響應(yīng)與漂移力的計算情況。KVLCC2實物模型如圖16[21]所示，幾何參數(shù)如表4所示。

表4 KVLCC2船舶幾何參數(shù)

圖16 KVLCC2實物模型[21]

驗證的工況共有2個，分別是迎浪和艏斜浪135°，航速為0。與切片方法不同的是，面元法的分析需要網(wǎng)格的支持。根據(jù)KVLCC2的船體特征，船首、球鼻艏和船尾的曲率較大，應(yīng)該生成較小尺寸的網(wǎng)格以捕捉到船體曲面，平行中體處的曲率變化較小，可以將網(wǎng)格尺寸適當(dāng)調(diào)大。

KVLCC2各曲面的網(wǎng)格尺寸如圖17所示。網(wǎng)格劃分結(jié)果如圖18所示。圖19為面元法對180°浪向下的KVLCC2運動響應(yīng)計算結(jié)果，并與CFD結(jié)果以及三維勢流理論結(jié)果[22]形成對比。對于垂蕩運動，面元法的計算結(jié)果與CFD和三維勢流理論吻合良好。誤差較大的地方發(fā)生在1<λ/L<2的區(qū)間，也就是波長為1倍和2倍船長的工況。對于縱搖運動，面元法結(jié)果的總體趨勢與CFD及三維勢流理論一致。在長波工況下誤差較大。

Mesh用于指定表面是否自動網(wǎng)格化，勾選表示該表面進行網(wǎng)格的劃分；Type表示網(wǎng)格類型，目前只有Tri.一種，即三角形網(wǎng)格；Min.edge length用于指定三角形網(wǎng)格的最小邊長；Max.edge length用于指定三角形網(wǎng)格的最大邊長；Patch egde angle用于指定網(wǎng)格的接觸角度；Node merge tolerance指相對容差，用于指示兩個相鄰的網(wǎng)格節(jié)點是否由于過于靠近而合并，該值應(yīng)在0.001～0.05；Geometric tolerance指幾何公差，用于指定網(wǎng)格邊可以偏離與其關(guān)聯(lián)的曲面的距離，該值應(yīng)在0.001～0.05，較低的值將在高曲面曲率區(qū)域創(chuàng)建較小的三角形；Surface Mesh Parameters為表面網(wǎng)格參數(shù)；Surface Name為表面名稱；Characteristic Leng為特征長度

Zero pt.為零點；AP為艉垂線；FP為艏垂線；Baseline為基線；CF為漂心；CG為重心；CB為浮心

圖19 180°浪向下的響應(yīng)幅值算子

對于在波浪中做搖蕩運動的船舶來說，作用在KVLCC2的外力包括入射力、繞射力、輻射力和靜水回復(fù)力。CFD結(jié)果基于黏流理論，考慮了黏性阻尼的影響。三維勢流理論在求解邊界積分方程時，用無航速的脈動源格林函數(shù)代替難以計算的移動源格林函數(shù)，在理論上與有航速船舶的水動力問題并不完全吻合。因此對前進速度效應(yīng)的準(zhǔn)確預(yù)測是格林函數(shù)的主要局限。面元法計算船舶在波浪中的運動響應(yīng)時，將船體周圍流場速度勢分解為波浪入射勢、自身繞射勢、相鄰船體產(chǎn)生的繞射勢、自身輻射勢和相鄰船體運動產(chǎn)生的輻射勢，通過源匯分布法進行求解。相比CFD方法，面元法能夠大幅減少計算成本；相比三維勢流理論，面元法通過對網(wǎng)格單元進行壓力積分得到二階漂移力和力矩，計算內(nèi)容更加豐富。

3 結(jié)論

采用Maxsurf Motions模塊對Series-60船舶、SA船舶、Wigley船型以及KVLCC2船舶的耐波性進行了驗證，將RAOs曲線、運動相位滯后、波浪增阻、附加質(zhì)量以及阻尼系數(shù)與實驗值及其他方法所得到的值進行對比。得到如下結(jié)論。

(1)Motions對縱搖響應(yīng)的計算精度要優(yōu)于垂蕩響應(yīng)。當(dāng)船舶處于垂蕩和縱搖的共振頻率范圍時，船體運動最復(fù)雜，Motions的計算誤差較大。當(dāng)遭遇頻率較低或較高時，Motions的計算結(jié)果與實驗值較吻合。

(2)附加質(zhì)量的計算中，Motions對中高波頻的計算精度較高，對部分低波頻段的計算精度較差。阻尼系數(shù)的計算中，隨著航速的增大，Motions的計算值與實驗值的偏離度增大，航速效應(yīng)對阻尼系數(shù)的影響要大于對附加質(zhì)量的影響，尤其是對縱搖模態(tài)的阻尼系數(shù)的影響最為明顯。

(3)Motions包含兩種耐波性計算方法，分別是切片理論和面元法。Motions對耐波性的計算精度本質(zhì)上是切片理論和面元法的計算精度。切片理論主要計算垂蕩和縱搖，使用條件是船體盡量細長。面元法可以計算6個自由度的運動，但只能針對零航速。