劉楊，龔家燁，劉家瑞，陳京普

(1.中國船舶科學研究中心上海分部，上海 200011；2.哈爾濱工程大學，哈爾濱 150001；3.中國船舶科學研究中心，江蘇 無錫 214000)

?

三體船舶阻力及航態(tài)預報方法

劉楊1，龔家燁2，劉家瑞3，陳京普1

(1.中國船舶科學研究中心上海分部，上海 200011；2.哈爾濱工程大學，哈爾濱 150001；3.中國船舶科學研究中心，江蘇 無錫 214000)

針對三體船型的自身特點，基于勢流理論，在Dawson方法的基礎(chǔ)上對數(shù)值處理進行改進，并運用于三體船舶的阻力預報。通過近水面網(wǎng)格的快速劃分與迭代計算，實現(xiàn)對航行姿態(tài)的預報。根據(jù)計算所得的船側(cè)興波，對不同航行姿態(tài)下的濕表面積進行計算，提高總阻力預報的精度。以三體Wigley船型為例，將理論計算結(jié)果和試驗值及商業(yè)軟件參考值進行對比。證明了理論預報方法的可行性與有效性，改進后的數(shù)值方法結(jié)果準確，可很好地避免自由面的震蕩。

三體船；Dawson方法；興波阻力；航行姿態(tài)；濕表面積

目前對于三體船型阻力的預報方法主要分為粘流理論以及勢流2種。粘流理論主要通過一些商用軟件(CFX、Fluent等)直接對總阻力進行計算，但是并未對興波阻力的產(chǎn)生機理進行探討，不利于針對興波阻力的優(yōu)化，同時在多種船型、布局優(yōu)化時需要大量的計算時間。勢流方面主要運用有Noblesse細長理論、薄船理論、Michel積分等方法，運用勢流理論可以達到快速預報的目的，但傳統(tǒng)方法并不能同時針對興波阻力、航行姿態(tài)、近場興波等多方面同時進行計算。為此，考慮在現(xiàn)有三體船阻力預報研究[1-5]基礎(chǔ)上，對傳統(tǒng)Dawson方法[6]進行改進，探討考慮航行姿態(tài)變化的興波阻力預報方法，并通過對濕表面積的計算，提高總阻力預報的精度。通過對數(shù)值計算方法[7-12]的對比分析，根據(jù)三體船型自身特點對網(wǎng)格的處理進行改進，采用控制點上移等方法實現(xiàn)在較高航速下得到光順的自由面波形。運用對近水面網(wǎng)格實時迭代重新劃分與表面壓力回歸分析的方法提高了船體受力與濕表面積計算的精度。

1 數(shù)學模型及數(shù)值方法

1.1 坐標系選取

采用隨船坐標系，見圖1。坐標系原點取船舶中心，z軸垂直向上，x軸由船首指向船尾，y軸指向右舷為正，xy平面與靜水面重合。

圖1 隨船坐標系

1.2 基本方程

假設(shè)流體理想，流動定常，無限水深，船舶興波為微幅波。流場速度勢ψ分解為無自由面影響的重疊模速度勢φ和考慮自由液面作用的擾動速度勢φ，即

(1)

(2)

基于上述假設(shè)，φ?φ。

速度勢滿足以下定解條件

1)流場中。

(3)

2)物面上。

(4)

3)自由面上。

(5)

式中，n為物面法向量。

同時，速度勢滿足無窮遠處趨于零的遠方輻射條件。

自由面上的φ關(guān)于靜水面對稱，則有

(6)

先后在在z=ξ和z=0上將速度勢進行泰勒展開，忽略φ的非線性項及高階項，并結(jié)合流線假設(shè)▽φ▽φ=φlφl，可得簡化的自由面邊界條件

(7)

式中，l為重疊模流線方向。

1.3 數(shù)值離散

將船體表面及自由面進行數(shù)值離散，取船體網(wǎng)格數(shù)為NB，自由面網(wǎng)格數(shù)Nf，以四邊形網(wǎng)格為主，部分不規(guī)則區(qū)域采用三角形網(wǎng)格，自由面網(wǎng)格劃分采用貼水線網(wǎng)格，對于片體間網(wǎng)格采用均勻大小網(wǎng)格，船體外側(cè)橫向及縱向網(wǎng)格則按一定比例增大。

對計算邊界進行數(shù)值離散后，任意控制點p處重疊模流場速度勢勢表示為

(8)

式中，SB表示物面；σB0為物面控制點處源強；rpq為任意源點q到控制點p的距離。

由船體與自由面相互作用產(chǎn)生的擾動速度勢表示為

(9)

對于自由面條件的求解，是該方法求解過程中的關(guān)鍵一步，尤其對于高性能船舶，自由面方程中差分格式尤為重要。本方法采用改進后的差分格式，將流線差分坐標變換為關(guān)于橫縱向網(wǎng)格的差分。因此自由面條件中流線速度勢的二階偏導數(shù)φll，φll表示為

(10)

縱向差分統(tǒng)一采用四點向前差分，對于自由面前緣單元，采用虛單元的方法進行計算；橫向差分計算時，對于片體與主體間網(wǎng)格采用中心差分，而對于船體外側(cè)的網(wǎng)格則采用三點向前差分進行計算。具體差分格式如下。

(11)

式中，(?L/?x)，(?T/?x)通過坐標轉(zhuǎn)換得，轉(zhuǎn)換矩陣如下

(12)

求解后得到物面源強σB1和源強σF。

1.4 深沉、縱傾的計算

船體表面控制點壓力系數(shù)為

(13)

式中，zi為控制點到靜水面間距離。

船體受到的升力和縱傾力矩可表示為

(14)

(15)

式中，xi,zi為控制點坐標；xg,zg為船舶的重心坐標；nxi,nzi為單元法線方向。

根據(jù)縱傾力矩及升力確定縱傾、升沉，快速劃分近水面局部網(wǎng)格，進行快速迭代計算，忽略大部分網(wǎng)格重新計算的內(nèi)存浪費，給定時間步長和收斂條件值，迭代得到定常狀態(tài)下的深沉縱、傾值。

1.5 數(shù)值處理

根據(jù)興波波長的不同采取不同的自由面網(wǎng)格大小。隨著弗汝德數(shù)增加網(wǎng)格數(shù)減小，但至少保證船長方向有25個網(wǎng)格，在滿足精度的基礎(chǔ)上，加快計算速度。經(jīng)過分析計算，與船體相連的第一層網(wǎng)格的寬度取0.15Lwl，片體間網(wǎng)格采用均布網(wǎng)格，并進行一定程度的加密；船體外側(cè)的網(wǎng)格，橫向、縱向增長比例rx、ry均取1.03。

由于三體船舶片體之間相互興波干擾較大，單元對于坐落于單元中心的控制點的切向誘導速度為零，在自由面條件向前差分過程中容易產(chǎn)生自由面震蕩。因此，在處理自由面控制點分布時，將所有控制點根據(jù)單元網(wǎng)格長度，進行向前移動，并且根據(jù)不同航速及區(qū)域的單元長度進行控制點上移，以消除震蕩。

采用船體表面壓力積分的方式求解船舶阻力，并基于粘流理論計算粘性阻力，得到船舶總阻力。在表面壓力積分和濕表面積計算時，結(jié)合直舷假設(shè)，根據(jù)船側(cè)興波對船體濕表面進行重新離散與細化，并對垂向一列單元的壓強進行線性回歸，進而通過垂向積分得到船體的興波阻力及濕表面積。這樣既可以提高阻力和濕表面積的計算精度，也無需增加計算過程中的物面網(wǎng)格數(shù)量，節(jié)省了計算時間。

2 算例及結(jié)果分析

2.1 排水式船體數(shù)值計算及驗證

采用三體Wigley船型，對不同的數(shù)值處理進行探討，從而得出較優(yōu)的參數(shù)設(shè)置。Wigley船型其數(shù)學表達式為

(16)

式中，L為船長；B為船寬；D為吃水，主體L/B=12.5,L/D=28.1；側(cè)體L1/B1=12.5,L1/D1=28.1。船體基礎(chǔ)網(wǎng)格數(shù)660個，自由面2 376個(Fr=0.25)。網(wǎng)格數(shù)隨弗勞德數(shù)變化聯(lián)動調(diào)整，網(wǎng)格劃分如圖2、3。

圖2 Wigley船體網(wǎng)格劃分

圖3 自由面網(wǎng)格(Fr=0.25)

探討了在三體船舶阻力及航態(tài)數(shù)值模擬時，控制點升高高度對數(shù)值計算結(jié)果影響。取Dx為單元長度，計算方案如表1。

表1 控制點上移方案

將方案1、方案2、方案3計算結(jié)果進行對比，阻力、縱傾及升沉的對比結(jié)果和波形模擬見圖4～6。方案3的計算所得波形見圖7。

圖4 計算興波阻力系數(shù)同試驗剩余阻力系數(shù)比較(方案1、2、3)

圖5 縱傾曲線(方案1、2、3)

圖6 升沉曲線(方案1、2、3)

通過圖4的對比結(jié)果可以看出，理論計算所得興波阻力系數(shù)與試驗所得剩余阻力系數(shù)趨勢及峰谷值吻合良好，而同圖5、圖6中深沉縱傾的對比可以明顯看出控制點上移比例為Dx/150的方案3計算精度最高。對比結(jié)果表明，控制點上移方案3可以在消除自由面震蕩的前提下保證計算精度。圖5表明，方案3計算所得的自由面波形也與實際符合良好的，并且光順。

由上述對比可知，方案3優(yōu)于方案1、方案2，因此繼續(xù)將方案3同約束模型方案4及試驗和商業(yè)軟件參考值結(jié)果對比，如圖8～10，方案3計算所得的濕表面積變化見圖11。

圖8 計算興波阻力系數(shù)同試驗剩余阻力系數(shù)比較(方案3、4)

圖9 升沉曲線(方案3、4)

圖10 縱傾曲線(方案3、4)

圖11 濕表面積變化曲線(方案3)

圖8表明，由于忽略了升沉、縱傾的影響，方案4的計算結(jié)果明顯小于程序和商業(yè)軟件結(jié)果。圖9、10表明，總體上本方法計算結(jié)果并不亞于商業(yè)軟件的非線性結(jié)果，且本方法可以再任意航速下得到光滑而完整的自由面興波，不存在自由面發(fā)散的問題。根據(jù)圖9的深沉結(jié)果，本方法計算結(jié)果優(yōu)于參考結(jié)果，與試驗值也較接近，這對于粘性類流體力的估算是很有價值的。但是在圖10中可以看出，在弗勞德數(shù)0.5左右，理論預報方法所得的縱傾角相對較大，其主要原因時忽略了自由面高階項的作用，對自由面興波造成一定的影響，也影響了壓力分布的精度。圖11的濕表面積變化曲線則可以看出，隨著航速的增大，三體船舶的濕表面積變化比例較大，因此在總阻力估算時，準確計算濕表面的影響是十分必要的。

2.2 不同片體布局的數(shù)值模擬

運用上述方法及優(yōu)化后的數(shù)值參數(shù)，對相同三體Wigley船型不同布局下的三體船舶進行數(shù)值模擬和計算。三體船片體布局示意圖見圖12，具體方案見表2。

圖12 片體布局

表2 片體布局方案

表2中，L為三體船主體的水線長度，a為片體中心縱向位置與船長的關(guān)系。片體中心縱向位置與主體相同時a為零，片體向船艉移動為正。

首先對方案1與方案3的計算結(jié)果與試驗值，結(jié)果見圖13～16。

圖13 計算興波阻力系數(shù)計算值與試驗剩余阻力系數(shù)比較(方案1)

圖14 興波阻力系數(shù)(方案3)

圖15 升沉曲線(方案1)

圖16 升沉曲線(方案3)

從圖13～16可以看出，計算興波阻力系數(shù)與試驗所得剩余阻力系數(shù)吻合良好，升沉值也吻合較好，可以看出本方法對于不同片體布局的三體船舶的阻力及航態(tài)變化可以準確進行預報，對于片體位置變化所造成的影響是比較敏感的，且在較大的弗汝德數(shù)范圍內(nèi)能保證精度，這對于高速三體船舶的船型優(yōu)化及片體布局優(yōu)化是有較大實際意義的。

方案1～方案4中不同布局的三體船舶計算值對比見圖17～20。

圖17 計算興波阻力系數(shù)同試驗剩余阻力系數(shù)比較(方案1～方案4)

圖18 升沉對比

圖19 縱傾對比

圖20 濕表面積變化對比

對比結(jié)果中可以明顯看出，片體后移的方案3、方案4隨著弗汝德數(shù)的增加阻力性能優(yōu)于方案1、方案2，而且在深沉和濕表面積上也小于方案1、方案2，這對于減小粘性阻力也是有益的，計算結(jié)果與實際情況是相符的。理論方法可以很好地對不同片體布局的三體船舶進行預報。

方案1、方案4的興波波形圖見圖21、22。

圖21 方案1波形圖(Fr=0.496)

圖22 方案4波形圖(Fr=0.496)

圖21、22中所反映出的興波情況與計算結(jié)果時相符，在Fr=0.496時，前述計算結(jié)果方案4的興波阻力系數(shù)也是大大小于方案1的，而且即使在弗汝德數(shù)較大的情況下，本方法所用的數(shù)值處理仍然可以得到較良好的波形，未發(fā)生自由面震蕩的情況。

3 結(jié)論

1)通過數(shù)值計算方法改善，可以在非全非線性自由面條件下改善數(shù)值預報結(jié)果，克服很多方法在自由面模擬的發(fā)散問題，獲得三體船阻力及航行姿態(tài)較好的預報結(jié)果。

2)本方法可在三體船設(shè)計階段為片體布局及船型初步優(yōu)化提供技術(shù)支撐。

3)由于三體船側(cè)體與主體與片體之間存在較強的興波干擾問題，因此傳統(tǒng)的模型很難準確地對非對稱流場進行模擬，在后續(xù)的工作中，應(yīng)當考慮如何采用有升力模型對非對稱有升力流場進行數(shù)值模擬，進而得到更準確的流場分布及數(shù)值結(jié)果。

[1] 何術(shù)龍,李百齊,程明道,等.三體船船型分析及興波干擾的模型試驗研究[J].水動力學研究與進展:A輯,2006,21(1):122-129.

[2] 吳廣懷,吳培德,蔣耀軍,等.基于興波阻力的三體船片體位置快速優(yōu)化方法[J].船舶力學,2005,9(4):1-8.

[3] 酈云,盧曉平.三體船阻力性能研究[J].船舶力學,2007,11(2):191-198.

[4] 蔡新功,王平,謝小敏.三體船方案優(yōu)化布局的阻力計算與試驗研究[J].水動力學研究與進展:A輯,2007,22(2):202-207.

[5] 陳京普，朱德祥，何術(shù)龍.雙體船/三體船興波阻力數(shù)值預報方法研究[J].船舶力學,2006,10(2)：23-29.

[6] DAWSON C W. A practical computer method for solving ship-wave problems[C]. The second International Conference on Ship Hydrodynamics.Beckly,1977:30-38.

[7] RAVEN H C. Variations on a theme by Dawson[C]. 17th Symposium on Naval Hydrodynamics. The Hague NetheRlands,1988:151-171.

[8] Md Shahjada Tarafder. Numerical calculation of free-surf-ace potential flow around a ship using the modified Rankine source method[J]. Ocean engineering,2008,35:536-544

[9] RAVEN H C. A solution method for the nonlinear ship wave resistance problem[D]. Delft, Netherlands, Delft University of technology,1996.

[10] 賀五洲，程明道，俞漢祥.方艉艦船及加裝尾板的興波阻力計算[J].中國造船,1998(2)：13-18.

[11] 劉應(yīng)中.船舶興波阻力理論[M].北京：國防工業(yè)出版社，2003.

[12] 盧曉平，王中，孫永華，等.Rankine源Dawson型方法求解三體船興波阻力[J].華中科技大學學報(自然科學版),2008,36(11):103-107.

Research on Prediction Method for Resistance and Hull Gesture of Trimaran

LIU Yang1, GONG Jia-ye2, LIU Jia-rui3, CHEN Jing-pu1

(1.China Ship Scientific Research Center, Shanghai 200011, China;2.Harbin Engineering University, Harbin 150001, China;3.China Ship Scientific Research Center, Wuxi Jiangsu 214082, China)

In light of characteristics of trimaran, some improvements were made on the numerical treatment of Dawson’s method based on potential flow theory to predict the resistance. By quick re-mesh of the grids near the water-plane in iterative calculation, the prediction of hull gesture was realized. On the basis of the calculated wave profile along the ship, the wetted surface area of trimaran in different hull gestures was calculated, and the prediction accuracy of total resistance was improved. With the example of tri-Wigley hull trimaran, theoretical calculation result was compared with the experimental and numerical results. The improved numerical method was proved to be feasible and effective, which can prevent the oscillation of the free surface.

trimaran; Dawson’s method; wave making resistance; hull gesture; wetted surface

10.3963/j.issn.1671-7953.2017.04.001

2016-09-12

劉楊(1990—)，男，碩士，助理工程師

研究方向：船舶阻力與推進技術(shù)

U661.3

A

1671-7953(2017)04-0001-06

修回日期：2016-10-11