小水線面雙體船阻力及航態(tài)預(yù)報方法

龔家燁，李云波，常赫斌

(1.哈爾濱工程大學(xué)船舶工程學(xué)院，黑龍江哈爾濱150001;2.中國船舶及海洋工程設(shè)計研究院，上海200011)

小水線面雙體船(small waterplane area twin hull，SWATH)波浪中良好的試航性能使其在軍用和民用上均具有良好的應(yīng)用前景。小水線面船型的不足是濕表面積大導(dǎo)致阻力大，航行姿態(tài)受排水量影響敏感，小的排水量變化會引起較大的吃水和姿態(tài)變化，航行中容易出現(xiàn)埋首現(xiàn)象。這些涉及到船舶安全航行的問題也是船型研發(fā)和設(shè)計中要解決的關(guān)鍵技術(shù)問題。

解決上述問題需要能夠準(zhǔn)確預(yù)報小水線面雙體船的阻力、航行姿態(tài)、姿態(tài)變化對阻力的影響、船體周圍的流線分布情況，進而指導(dǎo)船型優(yōu)化、鰭等附體布置及運動姿態(tài)改善等。而目前對于小水線面雙體船型阻力的預(yù)報方法多為粘流理論、薄船理論以及一些商用軟件。粘流方法難以達到快速預(yù)報的目的，薄船理論則難以對船體周圍興波以及航行姿態(tài)進行模擬，而Rapid等商用軟件常常需要通過升高面元特殊處理來實現(xiàn)自由表面的收斂。因此，探討有效的理論方法來支撐上述研究與預(yù)報在目前小水線面船型研發(fā)設(shè)計中是迫切需要的。

本文在現(xiàn)有小水線面雙體船阻力預(yù)報研究的基礎(chǔ)上［1－4］，基于Dawson方法［5］，針對上述船舶阻力、運動姿態(tài)和流線分布預(yù)報問題，探討理論預(yù)報方法與數(shù)值計算方法［6－12］。采用變化的邊界網(wǎng)格，在船舶運動姿態(tài)變化過程中實時進行濕表面積劃分及迭代，給出計及航行姿態(tài)和興波變化影響的船舶阻力預(yù)報方法;通過不斷迭代給出航行的升沉和縱傾值。通過船體周圍流線模擬和迭代完善，預(yù)報船體近遠(yuǎn)場興波、興波阻力及總阻力。

研究表明，本方法可進行快速的船舶周圍流場、興波、運動姿態(tài)和阻力預(yù)報，可從興波及干擾、流場分布、船舶阻力和運動姿態(tài)等角度進行船體水動力作用原理分析，支撐進行船型優(yōu)化、附體布置，進而研發(fā)設(shè)計出具有較好水動力性能的船型具有重要的理論和實際意義。

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

1.1 基本模型

采用隨船坐標(biāo)系，如圖1所示。坐標(biāo)系原點o位于船舶中心，z軸垂直向上，x軸由船艏指向船艉，y軸指向右舷為正，xy平面與靜水面重合。

圖1 坐標(biāo)系Fig.1 Coordinate system

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

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

流場中有

物面上滿足:

式中:n=nxi+nyj+nzk為物面法向量.

自由面z=ζ(x，y)上滿足:

速度勢滿足遠(yuǎn)處趨于零的遠(yuǎn)方輻射條件。

1.2 流線速度勢自由面邊界條件

考慮到φ關(guān)于靜水面對稱，則

在z=ξ和z=0上對速度勢展開，忽略φ的非線性項，可得z=0上簡化的自由面邊界條件:

假設(shè)▽φ▽φ=φlφl，可得重疊模流線速度勢表示的自由面條件:

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

1.3 奇點分布及數(shù)值計算模型

將船體和自由表面進行數(shù)值離散，取船體網(wǎng)格數(shù)為Nb，自由面網(wǎng)格數(shù)Nf，以四邊形網(wǎng)格為主，船體首尾局部采用三角形網(wǎng)格進行特殊處理，自由面采用貼水線網(wǎng)格，半自由面網(wǎng)格示意圖如圖2。

圖2 小水線面雙體船自由面網(wǎng)格Fig.2 Typical mesh for free surface of one SWATH

離散單元布置源匯后，控制點p處重疊模流場速度勢表示為

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

對式(9)離散，并代入物面不可穿透條件可得p趨近于物面時:

擾動速度勢φ也可用相同方法進行數(shù)值離散。

采用改進后的差分方式求解方程(8)，將流線差分坐標(biāo)變換為關(guān)于橫縱向網(wǎng)格的差分，沿縱向網(wǎng)格采用四點向前差分方式，橫向網(wǎng)格則采用中心差分來計算，因此自由面條件中流線速度勢的二階偏導(dǎo)數(shù)φll、φll的計算，以φll為例表示為

縱向差分中，當(dāng)i=1時，其上游網(wǎng)格i=0為虛單元，仍然計算其受到的單位誘導(dǎo)速度，但其源強為零。橫向差分中，與船體水線相鄰的網(wǎng)格j=1取與其相鄰網(wǎng)格j=2計算兩點差分，與邊界相鄰的網(wǎng)格仍然采用虛單元的方法計算中心差分。

式(12)中(?L/?x)、(?T/?x)通過坐標(biāo)轉(zhuǎn)換得

式中:J為雅各比矩陣。代入求解物面方程(4)和自由面條件(8)可得物面源強σB0、σB1和自由面源強σF。

1.4 興波阻力及深沉、縱傾的計算

船體表面控制點壓力系數(shù)Cpi和船體受到的興波阻力Rw、興波阻力系數(shù)Cw、升力Fz和縱傾力矩Ny表示為

式中:Sw表示船體的濕表面積，xi、zi為控制點坐標(biāo)，xg、zg為船舶的重心坐標(biāo)，nxi、nzi為單元法線方向。

給定時間步長Δt，可以得到船舶下一時刻的吃水St+Δt以及縱傾值Tt+Δt:

對船體水線以下近水面處局部網(wǎng)格進行實時劃分，而不發(fā)生出水或入水的大部分網(wǎng)格始終保持不變。由于小水線面雙體船型支柱相對規(guī)整的船型特點，使得其包括潛體在內(nèi)的絕大多數(shù)網(wǎng)格保持不變，僅僅需要對支柱接近水線部分網(wǎng)格進行重新插值與劃分，這樣既減小了網(wǎng)格重新離散帶來的誤差，也節(jié)省了計算時間。通過反復(fù)迭代直到滿足收斂條件(Fz＜ε、Ny＜ε)即可得到考慮船體姿態(tài)變化影響的興波阻力，進行船體姿態(tài)預(yù)報。

2 數(shù)值計算方法研究

2.1 波面網(wǎng)格劃分及探討

為了在較大弗汝德數(shù)范圍內(nèi)能夠獲得較準(zhǔn)確的興波波形，本研究對網(wǎng)格劃分及形式進行了研討。自由面網(wǎng)格沿縱向從船艏向遠(yuǎn)前方、船艉向后方呈輻射狀，網(wǎng)格長度按比例rx遞增，參考尺度為波長:

式中:Lwl為水線長。

單個波長網(wǎng)格數(shù)隨弗汝德數(shù)增加而減少，但至少保證一個波長30個網(wǎng)格，以免造成內(nèi)存的浪費或計算精度的不足。橫向網(wǎng)格尺度按ry向外遞增。

2.2 自由表面處理及模擬

為了準(zhǔn)確模擬興波，在自由面單元間切向誘導(dǎo)速度處理時，考慮到面源正負(fù)面切向速度連續(xù)性，將控制點上移一定高度以消除奇點影響獲得較好計算結(jié)果。上移距離與縱向網(wǎng)格尺寸成一定比例聯(lián)動變化。研究中發(fā)現(xiàn)，控制點位于單元中心時，單元對本身切向誘導(dǎo)速度貢獻約為零，這在數(shù)值差分時易造成自由面波形的震蕩及不穩(wěn)定，本研究進行了控制點前移處理，進行了不同前移比例對計算結(jié)果及色散影響研究。最后建議自由面所有控制點前移0.1個單元長度可以得到較好的數(shù)值模擬結(jié)果。

2.3 壓力積分影響及出水點的求解

研究發(fā)現(xiàn)，由于壓力積分時控制點壓力小于零的網(wǎng)格為出水網(wǎng)格，忽略其對阻力及航態(tài)的影響。因此，按照傳統(tǒng)的網(wǎng)格劃分方式垂向網(wǎng)格過少，很難較為精確地區(qū)分開船體出水部分，尤其對升力和船體深沉的計算結(jié)果有較大的影響。但是垂向網(wǎng)格過多又增大所需內(nèi)存和計算速度。因此本文探討了求得控制點壓強后，對垂直方向的每一列網(wǎng)格進行細(xì)化，對細(xì)化前所求得的各控制點壓強進行垂向線性回歸，插值得到細(xì)化后網(wǎng)格控制點的壓強及特殊點。通過上述處理可以提高船體運動姿態(tài)及不同航態(tài)下的濕表面積預(yù)報精度。

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

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

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

式中:L為船長，B為船寬，D為吃水，B=L/10，D= L/16。

沿船體基礎(chǔ)網(wǎng)格數(shù)1 584個，自由面2 968個(Fr=0.316)。網(wǎng)格大小隨弗汝德數(shù)變化進行聯(lián)動調(diào)整，網(wǎng)格劃分如圖3、4。

圖3 Wigley船體網(wǎng)格劃分Fig.3 Mesh for Wigley hull

圖4 自由面網(wǎng)格(Fr=0.316)Fig.4 Mesh for free surface of one Wigley(Fr=0.316)

探討了控制點升高高度對數(shù)值計算結(jié)果影響。取Dx為單元長度，共討論了4個方案。方案1～4分別為控制點升高 Dx/20、Dx/40、Dx/80及定值0.000 1的計算。

將方案1～3結(jié)果與試驗和Rapid商業(yè)軟件參考值進行對比，對比結(jié)果和波形模擬見圖5、6。

圖5 興波阻力系數(shù)(方案1～3)Fig.5 Coefficients of wave making ressitance(case 1-case 3)

圖6 方案3波形圖(Fr=0.316)Fig.6 Wave pattern of case 3(Fr=0.316)

圖5可以看出，弗汝德數(shù)低于0.35時，不同方案計算結(jié)果差別不大且與剩余阻力系數(shù)試驗值相符良好。但隨著弗汝德數(shù)增大，方案1、2的數(shù)值計算結(jié)果與試驗結(jié)果偏差逐漸增大。升高比例為Dx/80的方案3與試驗結(jié)果相符良好，且表現(xiàn)出較好的非線性效果。由上述可見，控制點升高比例不宜過大，選取Dx/80可以獲得較好的數(shù)值。從圖6也可以看出，方案3可以較好的模擬船體興波波形，且不會引起自由面的震蕩。

再將方案1～3中的最優(yōu)方案3、4及試驗和Rapid商業(yè)軟件參考值對比，如圖7～10。

圖7 興波阻力系數(shù)(方案3、方案4)Fig.7 Coefficients of wave making ressitance(case 3，case 4)

圖8 升沉曲線(方案3、方案4)Fig.8 Curves of sinkage(case 3、case 4)

圖9 縱傾曲線(方案3、方案4)Fig.9 Curves of trim(case 3，case 4)

圖10 船側(cè)面波形圖(方案3、方案4)(Fr=0.316)Fig.10 Wave profile along ship(case 3、case 4)(Fr=0.316)

圖7～10可見，控制點上移隨網(wǎng)格調(diào)整可提高計算精度，尤其對興波阻力和縱傾預(yù)報尤為明顯。由方案3結(jié)果與商業(yè)軟件非線性計算參考值進行比較也可以看出，采用本研究探討的數(shù)值計算方法，對非線性問題在較大弗汝德數(shù)范圍內(nèi)具有較好的預(yù)報精度。

由圖7可見，弗汝德數(shù)小于0.4范圍內(nèi)，本文方法無論數(shù)值計算結(jié)果還是曲線變化趨勢均與試驗結(jié)果有較好的相符度。由圖10可見，本研究模擬的興波波高與試驗值結(jié)果相符良好，僅船艏相對試驗結(jié)果較低，主要原因可能是忽略了粘性影響及貼船體處的網(wǎng)格劃分處理所致.

3.2 小水線面雙體船數(shù)值計算及結(jié)果分析

采用本研究方法對多型小水線面船型在弗汝德數(shù)0.1～0.45較大的范圍內(nèi)進行了對比分析和驗證。由于數(shù)值偏差和曲線變化趨勢是類似的，因此以某一型為例進行說明。船體和自由面網(wǎng)格數(shù)分別為1 500個和2 800個，船體網(wǎng)格示意圖見圖11。

圖11 某小水線面雙體船網(wǎng)格劃分Fig.11 Mesh for one SWATH

利用該理論方法預(yù)報了該船興波阻力、升沉和縱傾，采用粘流理論預(yù)報了粘性阻力，進而得到總阻力。將本方法所得數(shù)值預(yù)報結(jié)果與不考慮航態(tài)影響的約束模數(shù)值結(jié)果和試驗值進行比較，見圖12～15。

圖12可見，計及航態(tài)的興波阻力系數(shù)計算結(jié)果與試驗所得剩余阻力曲線的趨勢是完全吻合的，峰谷值也相互對應(yīng)。圖13的總阻力系數(shù)對比可以看出，約束模的總阻力系數(shù)計算結(jié)果與試驗值具有較大偏差?？紤]航態(tài)影響的總阻力系數(shù)曲線與在弗汝德數(shù)0.35之前幾乎穿過試驗點，且在峰值點的大小及相位上均與試驗結(jié)果相符良好。由于弗汝德數(shù)0.35相當(dāng)于水線長 100 m的小水線面雙體船約20 kn航速，可滿足多數(shù)該船型阻力預(yù)報和設(shè)計需求。弗汝德數(shù)0.35后理論預(yù)報結(jié)果略小于試驗值，最大偏差在15%。主要原因是高速高階效應(yīng)和粘性效果影響理論預(yù)報結(jié)果。

圖12 計算興波阻力系數(shù)與試驗數(shù)據(jù)比較Fig.12 Comparison of computed wave making resistance coefficient with experimental data

圖13 計算總阻力系數(shù)與試驗數(shù)據(jù)比較Fig.13 Comparison of computed total resistance coefficient with experimental data

圖14 計算縱傾值與試驗數(shù)據(jù)比較Fig.14 Comparison of computed trim with experimental data

從圖14、15可以看出，數(shù)值模擬的升沉、縱傾值與試驗值的變化趨勢相符?？v傾值在弗汝德數(shù)0.3左右有偏差，升沉預(yù)報與試驗結(jié)果符合較好。

該船型在Fr=0.3時產(chǎn)生的興波如圖16、17所示，可以看出興波也很好地滿足了遠(yuǎn)方輻射條件。

圖15 計算升沉值與試驗數(shù)據(jù)比較Fig.15 Comparison of computed Sinkage with experimental data

圖16 波形圖(Fr=0.3)Fig.16 Wave pattern(Fr=0.3)

圖17 波形云圖(Fr=0.3)Fig.17 Cloud picture of wave pattern(Fr=0.3)

4 結(jié)論

1)通過數(shù)值計算方法改善，可以在非全非線性自由面條件下改善數(shù)值預(yù)報結(jié)果，獲得小水線面雙體船阻力及航行姿態(tài)較好的預(yù)報結(jié)果。

2)采用本論文的理論預(yù)報方法，可以克服很多方法在自由面模擬的發(fā)散和震蕩問題，可以快速實現(xiàn)數(shù)值迭代收斂過程，計算結(jié)果穩(wěn)定。

3)利用本文的方法，對小水線面船舶可以較好的預(yù)報船舶航行姿態(tài)以及濕表面積變化，尤其對小水線面船型優(yōu)化和附體設(shè)計具有很好的指導(dǎo)作用。

4)由于采用無升力源匯模型，因此對于雙體之間的興波干擾難以準(zhǔn)確模擬，在今后的工作中需要對預(yù)報所采用的模型改善。

［1］鄒早建，羅青山，史一鳴.小水線面雙體船阻力預(yù)報研究［J］.中國造船，2005，46(1):14－21.

ZOU Zaojian，LUO Qingshan，SHI Yiming.A research on resistance prediction of SWATH ships［J］.Shipbuilding of China，2005，46(1):14－21.

［2］羊少剛，李干洛，余靈.高速雙體船興波阻力的數(shù)值計算［J］.華南理工大學(xué)學(xué)報:自然科學(xué)版，1995，23(10):16－25.

YANG Shaogang，LI Ganluo，YU Ling.Numerical calculation of wavemaking resistance of high－speed catamaran［J］.Journal of South China University of Technology:Natural Science，1995，23(10):16－25.

［3］黃鼎良.小水線面雙體船性能原理［M］.北京:國防工業(yè)出版社，1993.

［4］陳京普，朱德祥，何術(shù)龍.雙體船/三體船興波阻力數(shù)值預(yù)報方法研究［J］.船舶力學(xué)，2006，10(2):23－29.

CHEN Jingpu，ZHU Dexiang，HE Shulong.Research on numerical prediction method for wavemaking resistance of catamaran/trimaran［J］.Journal of Ship Mechanics，2006，10 (2):23－29.

［5］DAWSON C W.A practical computer method for solving shipwave problems［C］//Proceedings of the Second International Conference on Numerical Ship Hydrodynamics.Beckly，1977: 30－38.

［6］RAVEN H C.Variations on a theme by Dawson［C］//17th Symposium on Naval Hydrodynamics.The Hague Netherlands，1988:151－171.

［7］TARAFDER M S，SUZUKI K.Numerical calculation of freesurface potential flow around a ship using the modified Rankine source panel method［J］.Ocean Engineering，2008，35 (5/6):536－544.

［8］RAVEN H C.A solution method for the nonlinear ship wave resistance problem［D］.Delft:Delft University of Technology，1996.

［9］FLODéN D，KIM K，OTTOSSON P.A computational/experimental investigation on resistance and seakeeping characteristics of trimaran configuration in comparison with monohull［C］//The 9th Symposium on Practical Design of Ships and Other Floating Structures.Luebeck－Travemuende，Germany，2004.

［10］荻原誠功.Rankine sourceによる船體まわりの流れの近似計算法［J］.関西造船協(xié)會誌，1983:190.

［11］劉應(yīng)中.船舶興波阻力理論［M］.北京:國防工業(yè)出版社，2003:81－91.

LIU Yingzhong.Theory of ship wave making resistance［M］.Beijing:National Defense Industry Press，2003:81－91.

［12］王中，盧曉平，王瑋.考慮升沉和縱傾的方尾船非線性興波阻力計算［J］.水動力學(xué)研究與進展，2010，25(3): 422－428.

WANG Zhong，LU Xiaoping，WANG Wei.Nonlinear wave making resistance calculation for transom－stern ship with consideration of sinkage and trim［J］.Chinese Journal of Hydrodynamics，2010，25(3):422－428.