小水線面雙體船阻力及航態預報方法

龔家燁，李云波，常赫斌

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

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

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

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

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

1 數學模型及數值方法

1.1 基本模型

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

圖1 坐標系Fig.1 Coordinate system

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

基于上述假設，可知φ?φ。

流場中有

物面上滿足:

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

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

速度勢滿足遠處趨于零的遠方輻射條件。

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

考慮到φ關于靜水面對稱，則

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

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

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

1.3 奇點分布及數值計算模型

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

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

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

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

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

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

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

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

式(12)中(?L/?x)、(?T/?x)通過坐標轉換得

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

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

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

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

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

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

2 數值計算方法研究

2.1 波面網格劃分及探討

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

式中:Lwl為水線長。

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

2.2 自由表面處理及模擬

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

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

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

3 算例及結果分析

3.1 排水式船體數值計算及驗證

采用Wigley船型，對不同的數值處理進行探討，從而得出較優的參數設置。Wigley船型其數學表達式為

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

沿船體基礎網格數1 584個，自由面2 968個(Fr=0.316)。網格大小隨弗汝德數變化進行聯動調整，網格劃分如圖3、4。

圖3 Wigley船體網格劃分Fig.3 Mesh for Wigley hull

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

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

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

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

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

圖7 興波阻力系數(方案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 船側面波形圖(方案3、方案4)(Fr=0.316)Fig.10 Wave profile along ship(case 3、case 4)(Fr=0.316)

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

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

3.2 小水線面雙體船數值計算及結果分析

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

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

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

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

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

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

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

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

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

圖15 計算升沉值與試驗數據比較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 結論

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

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

3)利用本文的方法，對小水線面船舶可以較好的預報船舶航行姿態以及濕表面積變化，尤其對小水線面船型優化和附體設計具有很好的指導作用。

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

［1］鄒早建，羅青山，史一鳴.小水線面雙體船阻力預報研究［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］羊少剛，李干洛，余靈.高速雙體船興波阻力的數值計算［J］.華南理工大學學報:自然科學版，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］.北京:國防工業出版社，1993.

［4］陳京普，朱德祥，何術龍.雙體船/三體船興波阻力數值預報方法研究［J］.船舶力學，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］.関西造船協會誌，1983:190.

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

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

［12］王中，盧曉平，王瑋.考慮升沉和縱傾的方尾船非線性興波阻力計算［J］.水動力學研究與進展，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.