張藝瀚 王 平 蔡新功

（中國船舶及海洋工程設(shè)計研究院 上海200011）

考慮粘性的三體船非線性橫搖預(yù)報研究

張藝瀚 王 平 蔡新功

（中國船舶及海洋工程設(shè)計研究院 上海200011）

三體船主體與側(cè)體之間的流動干擾以及舭龍骨等附體造成的漩渦分離，使其在波浪中運動時的橫搖非線性效應(yīng)更加強烈而復(fù)雜。應(yīng)用CFD方法模擬某典型三體船在靜水中的衰減和強迫橫搖運動，得到非線性橫搖阻尼并分析其變化規(guī)律。應(yīng)用等效線性化原理修正并求解基于三維線性勢流理論建立的橫搖運動方程。研究結(jié)果與模型試驗比較后得出，采用粘性流CFD方法并結(jié)合三維勢流理論可較準(zhǔn)確地預(yù)報三體船的非線性橫搖阻尼和波浪中的橫搖運動。

三體船；非線性橫搖；三維勢流理論；粘流理論；計算流體動力學(xué)

引 言

三體船是一種新型的水面平臺，水下部分由主體和兩個側(cè)體共三個細長船體組成，兩個側(cè)體通過連接橋與中體連接成一體。三體船的船型特點使其具有良好的穩(wěn)性、阻力和耐波性等性能。橫搖運動性能是耐波性研究中的重點，對船舶的適居性、使用性、安全性等都有重要的影響。目前，國內(nèi)外單體船的橫搖運動預(yù)報主要包括二維切片法、三維頻域法、三維時域法等，其中對于非線性橫搖一直是難以完善的問題［1］。

目前，三體船的橫搖運動主要依靠船模試驗，國內(nèi)外針對三體船的橫搖運動的理論預(yù)報研究較少。Zhang J W研究了三體船非線性橫搖阻尼的處理方法，通過試驗和經(jīng)驗公式兩種方法計算非線性橫搖阻尼［2］，并對橫搖阻尼的非線性成分進行了系統(tǒng)分析［3］。李培勇等通過三體船模型的靜水橫搖試驗和規(guī)則波中的橫搖運動試驗，得到側(cè)體橫向位置以及舭龍骨的存在對三體船橫搖運動的影響［4-5］。對于橫搖運動的非線性計算方法，基于段文洋（Duan Wenyang）在二維半理論時域格林函數(shù)上的研究［6］，項久洋［7］運用三維時域Rankine源方法研究船型要素對三體船耐波性的影響，李輝等［8］研究基于等效線性原理的橫搖阻尼勢流計算法，針對不同船型指出非線性橫搖阻尼經(jīng)驗公式的適用性。

然而，這些理論預(yù)報結(jié)果無法準(zhǔn)確考慮三體船非線性效應(yīng)及流體粘性影響，與試驗值存在一定差距。

CFD方法可以在橫搖阻尼中有效考慮流體粘性的影響，Yu Yi Hsiang等［9］對四個不同的船體剖面進行自由橫搖衰減和強迫橫搖運動模擬，給出不同剖面不同橫搖幅值下的橫搖力矩時歷和橫搖阻尼系數(shù)，黃昊等［10］對系列60船型進行強迫橫搖數(shù)值模擬，提出通過強迫橫搖得到橫搖阻尼系數(shù)的方法，朱仁傳等［11-12］也應(yīng)用CFD方法進行單體船的橫搖運動數(shù)值計算，但針對三體船的三維CFD橫搖數(shù)值計算，目前仍沒有深入研究。

本文提出勢流理論結(jié)合粘性修正的方法，應(yīng)用CFD方法計算三體船的非線性橫搖阻尼，修正三維勢流理論建立的橫搖運動方程中，試圖求解三體船在波浪中的橫搖運動響應(yīng)。

1 CFD法求解靜水非線性橫搖阻尼

1.1 計算模型建立

應(yīng)用GAMBIT前處理工具建立三體船三維非結(jié)構(gòu)網(wǎng)格模型，在FLUENT平臺應(yīng)用UDF編程的動網(wǎng)格技術(shù)和滑移網(wǎng)格技術(shù)模擬靜水橫搖衰減運動和固定頻率下強迫橫搖運動，計算采用有限體積法建立離散控制方程，應(yīng)用VOF多相流模型，湍流模式采用標(biāo)準(zhǔn)k-ω模式。計算模型和計算域網(wǎng)格見圖1和圖2。

圖1 三體船計算模型

圖2 三體船計算網(wǎng)格

1.2 強迫橫搖運動模擬

橫搖強迫運動指強迫計算模型以一定的頻率和幅值進行橫搖運動，分析橫搖阻尼系數(shù)和頻率以及它們與幅值之間的關(guān)系。目前國內(nèi)外對于三體船型的二維剖面強迫橫搖運動研究較多，本文采用三維計算模型可以真實反映三體船橫搖阻尼中側(cè)體與主體干擾以及舭龍骨等附體產(chǎn)生的非線性部分。

采用滑移網(wǎng)格的方法，在模擬過程中不僅不需要進行網(wǎng)格重畫，而且可以有較大的時間步長，大大提高計算精度和效率。橫搖強迫運動采用UDF工具進行編程，使用FLUENT自帶的DEFINE_CG_ MOTION宏進行定義，模擬橫搖運動見式（1）：

圖3 強迫橫搖運動動畫捕捉

其中1.57 rad/s頻率下得到的橫搖力矩曲線如圖4所示。可見，隨著搖幅和頻率的增加，橫搖阻尼力矩逐漸增加，與文獻中的單體船具有相似規(guī)律。

圖4w=1.57 rad/s頻率下的橫搖力矩時歷曲線

根據(jù)上述橫搖力矩曲線可得到線性化橫搖阻尼系數(shù)：

當(dāng)橫搖角度為0°，將橫搖角速度最大時的瞬時橫搖力矩代入式（2）便可求出線性阻尼力矩，無因次化后得到：

此線性阻尼系數(shù)可視為已考慮了漩渦阻尼、流體粘性等非線性因素的等效線性化阻尼。計算模型的無因次阻尼系數(shù)見表1?？梢?，橫搖阻尼系數(shù)隨著波幅和頻率的增加而增大。

表1 三體船三維模型無因次橫搖阻尼系數(shù)

1.3 橫搖衰減運動模擬

船舶的靜水橫搖衰減試驗通常用于測定線性和平方橫搖阻尼系數(shù)，CFD方法可有效模擬三體船的靜水橫搖試驗，得到衰減曲線。

以上述計算模型為例，計算零航速下初始橫搖角為5°、9°和15°時的橫搖衰減運動，結(jié)果如圖5所示。

圖5 不同初始角下的橫搖衰減曲線

可見，不同初始角下的橫搖固有周期略有差距，初始角越大、衰減幅度越大，體現(xiàn)了大幅橫搖下的橫搖的非線性效應(yīng)。

基于Fround［13］及Haddara［14］等人的研究和目前普遍采用的形式，本文對阻尼力矩采用的是線性項加平方項的關(guān)系，因此其對應(yīng)關(guān)系如下：

式中：a、b為衰減系數(shù)；A、B為阻尼力矩系數(shù)；分別為橫搖角速度、衰減曲線相鄰兩振幅平均值、衰減曲線相鄰兩振幅之差。

根據(jù)衰減曲線得到橫搖角的消滅曲線并通過最小二乘法擬合得到衰減系數(shù)，由式（5）得到橫搖阻尼力矩系數(shù)（見表2）。

表2 橫搖阻尼

可見，隨著初始角度的增加，消滅系數(shù)a逐漸增加，而系數(shù)b呈減少的趨勢，符合理論計算規(guī)律［15］。

式中：D為船的質(zhì)量；h為船的橫穩(wěn)心高；為船自由橫搖的近似固有圓頻率。

1.4 橫搖阻尼比較分析

通過強迫橫搖運動結(jié)果可直接計算三體船在特定頻率和特定搖幅下的線性化橫搖阻尼，近似船舶在規(guī)則波中橫搖時的非線性橫搖阻尼。表3、表4為CFD計算值與試驗值的比較，其中試驗值采用衰減試驗得到的線性和平方阻尼，利用等效線性化方法［15］得到橫搖阻尼?？梢?，CFD方法與試驗值較為接近，誤差在5%以內(nèi)。

表3 零航速下強迫運動9°橫搖阻尼系數(shù)比較

表4 10 kn航速下強迫運動9°橫搖阻尼系數(shù)比較

CFD橫搖衰減運動模擬了船模水池靜水橫搖試驗，可實現(xiàn)橫搖線性阻尼和平方阻尼的求解，表5—表7為CFD計算值與船模靜水橫搖衰減試驗結(jié)果的比較。可見，隨著航速的增加，橫搖阻尼系數(shù)增加，且非線性橫搖阻尼的比例減少，數(shù)值計算結(jié)果與試驗值較為接近，誤差在5%以內(nèi)。CFD方法可以在沒有條件進行水池船模試驗時，在一定工程誤差范圍內(nèi)預(yù)報三體船的橫搖阻尼。

表5 零航速衰減運動橫搖阻尼系數(shù)比較

表6 10 kn航速衰減運動橫搖阻尼系數(shù)比較

表7 18 kn航速衰減運動橫搖阻尼系數(shù)比較

2 波浪中橫搖運動預(yù)報

三維線性勢流理論是預(yù)報船舶在波浪中運動的重要方法。針對三體船型，勢流理論建立的橫搖運動方程無法準(zhǔn)確考慮片體主體流動干擾、附體漩渦分離等橫搖非線性因素，而直接求解非線性橫搖運動方程在數(shù)學(xué)上較為困難。本文在三體船的非線性橫搖阻尼求解的基礎(chǔ)上，提出采用勢流理論結(jié)合CFD粘性修正的方法近似求解非線性橫搖運動方程，獲得不同航速規(guī)則波中的橫搖運動響應(yīng)。

2.1 線性運動方程

根據(jù)三維勢流理論，船舶在波浪中的運動總方程為：

水動力系數(shù)由式（8）求解：

2.2 非線性橫搖修正

式（6）中，橫搖線性阻尼系數(shù)為B44，其可通過式（8）求得。考慮非線性橫搖阻尼后，采用線性加平方阻尼力矩形式：

將上式直接代入運動方程后的非線性微分方程求解困難，本文采用等效線性化原理得到線性化橫搖阻尼力矩。

計算關(guān)鍵在于線性阻尼系數(shù)A和平方阻尼系數(shù)B。本文提出粘流結(jié)合勢流的方法，即采用CFD數(shù)值模擬得到線性和平方阻尼系數(shù)，等價線性化后計入橫搖運動方程中，取代原來的線性橫搖阻尼系數(shù)B44，這樣既包含興波阻尼，又包含摩擦阻尼、漩渦阻尼及流體干擾等非線性因素，在線性橫搖運動方程中體現(xiàn)非線性橫搖阻尼力矩的影響。另外，本文在橫搖方程求解時同時考慮恢復(fù)力矩的非線性，將非線性項展開，略去高階項得到線性化的恢復(fù)力矩系數(shù)。

2.3 實例計算與分析

根據(jù)上述勢流和粘流相結(jié)合的方法，針對上述三體實船計算其在規(guī)則波中的橫搖運動響應(yīng)。應(yīng)用三次樣條曲線建立三維面元水動力計算模型，如圖6所示。

圖6 目標(biāo)三體船水動力計算模型

CFD計算結(jié)果與試驗值、線性勢流理論值及船模衰減試驗修正值進行比較，其中船模衰減試驗修正值是指采用船模衰減試驗得到的橫搖阻尼系數(shù)代替數(shù)值計算結(jié)果進行橫搖運動方程的非線性修正。

在橫浪90°和斜浪下，計算三體船在0 kn、10 kn、18 kn下的橫搖運動頻響函數(shù)見圖7。其中規(guī)則波參數(shù)為：單位波幅、浪向角90°、斜浪45°、波浪圓頻率0.1～2.0 rad/s。

圖7 橫搖頻域響應(yīng)

從橫浪和斜浪規(guī)則波橫搖運動的預(yù)報情況來看，橫搖阻尼系數(shù)對響應(yīng)幅值影響較大。采用本文勢流粘流組合方法得到的非線性橫搖阻尼修正勢流理論橫搖方程的方法得到的結(jié)果與試驗值符合較好，并與通過船模試驗獲得阻尼系數(shù)修正方法結(jié)果較為接近。另外，如果不對線性橫搖阻尼進行非線性修正，則得到的響應(yīng)幅值普遍高于試驗值。

橫搖幅值在波浪周期接近船舶橫搖固有周期時達到峰值，符合理論規(guī)律。橫浪中的預(yù)報值高于斜浪中預(yù)報值，且隨著航速增加，橫搖阻尼增大，橫搖峰值逐漸減少，線性結(jié)果與非線性修正后結(jié)果之間的差距也在減少。這說明零速橫搖時非線性因素成分較大，隨著航速增加，線性阻尼占主導(dǎo)因素。

3 結(jié) 論

（1）三體船型較為特殊，橫搖阻尼的求解必須綜合考慮主體與側(cè)體的流體干擾、漩渦產(chǎn)生和分離以及附體阻尼等現(xiàn)象引起的粘流非線性效應(yīng)，而CFD數(shù)值模擬方法可較準(zhǔn)確求解三體船的橫搖運動，并考慮上述非線性。其中，強迫運動可以模擬三體船在固定頻率和搖幅下的橫搖運動，三維衰減運動可以模擬靜水橫搖運動，從而得到非線性橫搖阻尼系數(shù)。

（2）本文提出的勢流粘流組合方法可以較準(zhǔn)確地預(yù)報三體船在規(guī)則波中橫搖運動，既體現(xiàn)了勢流理論在求解方程的優(yōu)勢，又能考慮非線性效應(yīng)，比單純的三維勢流理論和切片理論更加精確，因此可作為三體船在波浪中橫搖運動預(yù)報方法在工程中廣泛應(yīng)用。另外，橫搖阻尼的非線性修正對縱搖、垂蕩等預(yù)報不會產(chǎn)生影響，也適合基于勢流理論的六自由度運動預(yù)報。

［ 1 ］劉宇辰，任毅，蔡新功，等.三體船橫搖運動試驗研究［J］.船舶，2011（2）：20-23.

［ 2 ］Zhang J W.Roll damping characteristics of a trimaran displacement ship［J］.International Shipbuilding Progress.1999，448：445-472.

［ 3 ］Zhang J W.Investigation into roll damping of trimaran ships［J］.NARG Report， 1996，1026： 225-276.

［ 4 ］李培勇，裘泳銘，顧敏童，等.超細長三體船耐波性試驗研究［J］.海洋工程，2002（4）：1-10.

［ 5 ］李培勇，馮鐵城，裘泳銘.三體船橫搖運動［J］.中國造船，2003（1）： 24-30.

［ 6 ］Duan Wenyang， Hudson D A， Price W G.Theoretical prediction of the motions of fast displacement vessels in long-crested head seas［C］//The 3rd International Conference for High Performance Marine Vehicles，Shanghai， China， 2000： 116-122.

［ 7 ］項久洋.船型要素對三體船耐波性和波浪載荷影響的數(shù)值計算［D］.武漢：武漢理工大學(xué)，2008：25-45.

［ 8 ］李輝，張藝瀚，任慧龍.非線性橫搖阻尼對波浪扭矩和運動的統(tǒng)一修正［J］.哈爾濱工程大學(xué)學(xué)報，2013（7）：1-6.

［ 9 ］Yu Yi Hsiang， Spyros A Kinnas.Prediction of hydrodynamic performance for various ship-shaped hulls under excessive roll motions using an unsteady Navier-Stokes solver［C］//Proceedings of the 18th IOPEC， Vancouver BC，Canada， 2008： 265-286.

［10］黃昊，郭海強，朱仁傳，等.粘流中船舶橫搖阻尼計算［J］.船舶力學(xué)，2008（4）：568-573.

［11］朱仁傳，楊春蕾，繆國平，等.粘性流中三維船體橫搖運動數(shù)值模擬及阻尼計算［C］//第十一界全國水動力學(xué)學(xué)術(shù)會議暨第二十四屆全國水動力學(xué)研討會周培源誕辰110周年紀(jì)念大會文集（下冊）.中國，2012：903-912.

［12］朱仁傳，郭海強，繆國平，等.一種基于CFD理論船舶附加質(zhì)量與阻尼的計算方法［J］.上海交通大學(xué)學(xué)報，2009（2）：198-203.

［13］Froude J E.Remarks on Mr Scott-Russell's paper on rollin［M］.The Papers of William Froude.1955： 112-118.

［14］Haddara M R， Bennett P A.A Study of the angle dependence of roll damping moment［J］.Ocean Engineering，1989（4）：26-32.

［15］李積德.船舶耐波性［M］.第2版.哈爾濱：哈爾濱工程大學(xué)出版社，1991：1-3.

［16］戴仰山，沈進威，宋競正.船舶波浪載荷［M］.第1版.北京：國防工業(yè)出版社，2007：1-4.

Nonlinear rolling prediction for trimaran ships considering viscous effect

ZHANG Yi-han WANG Ping CAI Xin-gong
(Marine Design & Research Institute of China, Shanghai 200011, China)

The nonlinear effect is more intense and complex for the trimaran rolling in waves due to the disturbance between the main hull and side hulls, and the vortex shedding that is caused by the appendages as bilge keels.The roll decay and forced rolling in the calm water for a typical trimaran are simulated by CFD methods.The nonlinear roll damping is then obtained for the analysis of the roll damping rule.The trimaran rolling equation which is constructed by the three-dimensional (3D) potential flow theory is modified and solved by applying the equalization linearised theory.By comparison with the model test, the results show that the nonlinear roll damping and rolling in waves can be accurately predicted through the viscous flow CFD method combined with 3D potential flow theory by comparison with the model test results.

trimaran; nonlinear rolling; 3D potential flow theory; viscous flow theory; CFD(computation fluid dynamics)

U661.32

A

1001-9855（2016）06-0012-07

2015-09-09；

2016-07-27

張藝瀚（1987-），男，工程師，研究方向：艦船總體設(shè)計。王 平（1965-），男，研究員，研究方向：艦船總體設(shè)計。蔡新功（1974-），男，研究員，研究方向：艦船總體設(shè)計。

10.19423/j.cnki.31-1561/u.2016.06.012