李井煜，盧曉平

海軍工程大學艦船工程系，湖北武漢430033

基于Rankine源和Kelvin源格林函數(shù)求解興波阻力的復合算法

李井煜，盧曉平

海軍工程大學艦船工程系，湖北武漢430033

［目的］運用邊界元法計算船舶興波阻力基本上是先求解船體附近的速度分布，然后采用伯努利方程進行壓力積分，其計算過程復雜，且誤差非常大。［方法］提出一種可快速計算船舶興波阻力的復合算法，利用Rankine源格林函數(shù)求解船體表面源強，結(jié)合Lagally定理進行受力計算，并基于Kelvin源格林函數(shù)求解船舶興波阻力。運用該算法對Wigley船的興波阻力進行計算。［結(jié)果］計算結(jié)果表明，所用算法相較于運用線性興波阻力中的薄船理論得到的結(jié)果精度更高，而且與完全使用Kelvin源格林函數(shù)的算法相比效率也更高。［結(jié)論］所用算法可在計算興波阻力時作為精度與效率之間的一種折中方法。

線性興波阻力；格林函數(shù)；Rankine源；Kelvin源；邊界元法

0 引 言

邊界元法（Boundary Element Method，BEM）已被廣泛應用于船舶興波阻力的數(shù)值計算中，該方法可分為間接邊界元法（分布源方法，或稱間接法）和直接邊界元法（源—偶混合分布法，或稱直接法）。間接法可直接、方便地求解得到流場中的速度，且能推導并用于流體力學中常用到的源匯強度概念，因此在各種商用軟件得到了廣泛應用［1］。這兩種方法的計算效率相差無幾，但對計算精度而言，直接法略高于間接法。近年來，間接法的發(fā)展遇到了一些瓶頸。例如，文獻［2］提到間接法在求解非光滑邊界處的切向誘導速度時計算精度極差，而對于非線性問題，其計算精度很低，甚至導致發(fā)散，故近年來有學者提出發(fā)展直接法來提高計算精度［3-6］。

無論是使用直接法還是間接法，都涉及到格林函數(shù)（基本解）的選擇問題。Kelvin源格林函數(shù)可滿足線性自由面興波條件，在基于此格林函數(shù)的船舶興波阻力計算中，可以不需要在流場自由液面布置源匯。Rankine源格林函數(shù)形式簡單，但在基于此格林函數(shù)的船舶興波阻力計算中，需要在所有物面布置源匯，故僅能近似地在有限區(qū)域求解興波問題［7-8］。實際上，若僅關(guān)注船舶航行時的興波阻力，計算求解出物面上的源匯分布是關(guān)鍵。

本文將基于Lagally定理和線性興波阻力理論，提出一種新復合算法。該算法運用Rankine源格林函數(shù)計算船體表面的源強密度，采用Kelvin源格林函數(shù)計算興波阻力。與線性興波理論的薄船理論相比，復合算法的計算精度較高，而與完全使用Kelvin源格林函數(shù)的方法相比，復合算法的計算效率更高。本文所用算法可作為興波阻力計算時精度與效率之間的一種折中方法。

1 數(shù)值計算

本文假定船舶在無限深廣的自由水面上，以航速U在x軸負向做勻速直線運動。流場滿足均勻、不可壓縮和無粘、無旋流動條件，且船體興波屬于微振幅波。取固結(jié)在船上的隨船坐標系，如圖1所示。

圖1 坐標系Fig.1 Coordinate system

控制方程為不可壓縮、流動的連續(xù)性方程，邊界條件為物面不可穿透條件，表達式如式（1）～式（3）所示。

式中：?為哈密頓算子；Φ為流場總速度勢；φ為船體擾動速度勢；n為船體法線方向向量。

由拉蓋爾定理可知，流體對物體的作用力F可認為是多個點源i對其產(chǎn)生作用力的總和，即

式中：ρ為流體密度；mi為點源的源強；qi為點源處的速度。

針對上述船舶興波問題，本文認為流場由物面分布的點源誘導產(chǎn)生，各單元的源強σidA（σi為各單元的源強密度，A為面積）將受到沿航行方向的阻力影響，故船體興波阻力Rw被認為是作用于整個分布的源匯上的作用力積分之和，即

式中，V為單元內(nèi)流場沿x軸方向的速度分量。

根據(jù)勢流理論，船體擾動速度勢φ由船體表面分布的源匯誘導產(chǎn)生，為

式中，Gs為單位源強產(chǎn)生的速度勢，或稱格林函數(shù)，下標s表示某面元上的格林函數(shù)。在興波問題中，常見的2種格林函數(shù)包括Kelvin源格林函數(shù)和Rankine源格林函數(shù)。Kelvin源格林函數(shù)滿足線性自由面條件，表達式如式（7）所示。

式中：(x，y，z)為場點坐標；(ξ，η，ζ)為源點坐標；r1，r2，K，K0，θ均為計算的過程量，其中：

式中，g為過程量。相比于Kelvin源格林函數(shù)，Rankine源格林函數(shù)不滿足自由面興波條件等有限域內(nèi)的邊界條件，但可簡化為：

另一方面，將式（2）和式（6）代入式（5），可得到如下興波阻力表達式：

根據(jù)文獻［9］，式（8）的第1項為0，故興波阻力的表達式可進一步簡化為

將式（7）代入式（9）可得到4個積分項。因積分區(qū)域是關(guān)于點源和場點對稱，而根據(jù)文獻［9］的分析，只有第4項對興波阻力起作用，故其表達式為：

式中，P函數(shù)和Q函數(shù)為計算過程量，故只需求解得到船體物面上的源強密度分布，即可求解興波阻力。根據(jù)線性興波理論的薄船理論，船體表面的源強密度直接與船體形狀相關(guān)。但是，根據(jù)Neumann-Kelvin理論或其他使用Kelvin源格林函數(shù)的非線性興波阻力理論，計算源強分布的過程十分復雜。從整體上看，后者的精度要高于前者。在物理直觀上，前者的精度可以通過簡化的線性興波理論得到提高，即通過使用Rankine源格林函數(shù)求解無限域中的疊模來確定各單元源強密度的分布。

綜上所述，本文擬提出一種求解興波阻力的快速復合算法。該算法使用Rankine源格林函數(shù)求解船體表面各單元的源強密度，結(jié)合基于Kelvin源格林函數(shù)的式（10）和式（11）求解興波阻力。本文所用算法實際上可在計算興波阻力時作為精度與效率之間的一種折中方法。其中，使用Rankine源格林函數(shù)求解源強密度可以使用直接邊界元法和間接邊界元法求解。根據(jù)直接邊界元法，求解源強密度的公式為

2 數(shù)值驗證

2.1 Wigley船型興波阻力計算分析

在Taylor圖譜［10］適用范圍內(nèi)選取一種Wigley船型作為第1個算例用以分析和驗證本文所用復合算法的有效性。該船型參數(shù)如表1所示。

表1 Wigley船型參數(shù)Table 1 Parameters of Wigley ship form

本文對船體的網(wǎng)格劃分使用三角形網(wǎng)格。船體劃分為40×4×2個網(wǎng)格單元?？紤]到上疊模部分的網(wǎng)格單元，單元總數(shù)為640個，單元信息如圖2所示。圖3為第1個算例中Wigley船型的興波阻力計算結(jié)果比較。

圖2 Wigle船型網(wǎng)格Fig.2 Panel mesh of the Wigley ship form

圖3 第1個算例中Wigley船型興波阻力計算結(jié)果比較Fig.3 Results comparison of the Wigley ship form in case 1

由圖3可以看出，結(jié)合了疊模技術(shù)、Rankine源格林函數(shù)法、Lagally定理和直接邊界元法（復合算法）的計算結(jié)果和Michell薄船線性興波阻力理論的計算結(jié)果相比更接近于Taylor圖譜計算值，尤其是在Fr＞0.35的高速階段。同時，還發(fā)現(xiàn)在低速階段存在直接邊界元法計算不準確的現(xiàn)象。下文將對該部分的誤差進行詳細分析。

第2個算例為另一種Wigley船型。該模型曾在605所進行過阻力試驗，模型參數(shù)如下：L=5 m，B=0.4 m，T=0.178 m。通過Prohaska方法確定模型形狀因子，進而可近似得到興波阻力的試驗值。圖4所示為第2個算例中Wigley船型的興波阻力試驗結(jié)果與計算結(jié)果的比較。

圖4 第2個算例中Wigley船型興波阻力計算結(jié)果比較Fig.4 Results comparison of the Wigley ship form in case 2

由圖4可以看出，基于復合算法的直接邊界元法和間接邊界元法的計算精度要高于Michell薄船線性興波阻力理論的精度，尤其是在中、高速段。對比直接邊界元法和間接邊界元法，發(fā)現(xiàn)直接邊界元法的計算精度比間接邊界元法的稍高。但正如前文所述，本文所用算法在Fr＜0.2的低速段計算誤差較大。表2所示為各算法精度的定量分析結(jié)果。

2.2 低速段誤差分析

分析式（11）可知，本文所用復合算法的誤差

表2 不同算法的相對誤差比較Table 2 Comparison of relative error of different methods

主要來自如下項，為

當速度較低時，K0很大，這將導致此項沿x軸方向產(chǎn)生高頻振蕩，指數(shù)項在很大程度上增加了振蕩幅值。圖5所示為某一水平線上式（11）中的Q函數(shù)隨x軸振蕩的趨勢，其中縱坐標表示

由Q函數(shù)的主要積分部分交換積分次序消去Q后組成。

圖5 不同速度段的振蕩項QFig.5 The oscillating termQin different speed range

針對這部分誤差，一方面可以通過增大x軸方向的網(wǎng)格劃分密度加以改進，另一方面需要進一步研究振蕩函數(shù)的積分算法來估算誤差。

3 結(jié) 論

本文提出了一種求解船舶興波阻力的復合算法，介紹了復合算法的原理，并驗證了該方法。通過數(shù)值算例，得到如下結(jié)論：

1）總體上，無論是使用基于直接邊界元法還是使用間接邊界元法的復合算法來求解船舶興波阻力，與Michell積分公式等傳統(tǒng)的數(shù)值計算方法相比都具有較高的精度，同時與基于Neumann-Kelvin理論的求解方法相比復合算法的效率更高，故可在興波阻力計算時作為精度與效率之間的一種折中方法。

2）直接邊界元法的求解精度略高于間接邊界元法。

3）在低速段，采用復合算法計算的精度較低，主要是因為該方法中的P函數(shù)和Q函數(shù)包含了振蕩項，故處理振蕩函數(shù)的方法仍需進一步研究。

［1］戴遺山，段文洋.船舶在波浪中運動的勢流理論［M］.北京：國防工業(yè)出版社，2008.

［2］徐剛.不規(guī)則波中浮體二階水動力時域數(shù)值模擬［D］.哈爾濱：哈爾濱工程大學，2010.XU G.Time-domain simulation of second-order hydro?dynamic force on floating bodies in irregular waves［D］.Harbin：Harbin Engineering University，2010（in Chinese）.

［3］MIAO Q M，CHWANG A T.Ship waves by direct boundary element method［J］.Journal of Ship Mechan?ics，2003，7（6）：27-36.

［4］陳紀康，段文洋，朱鑫.三維泰勒展開邊界元方法及其數(shù)值驗證［J］.水動力學研究與進展（A輯），2013，28（4）：482-485.CHEN J K，DUAN W Y，ZHU X.Three-dimension Taylor expansion boundary element method and it's val?idation［J］.Chinese Journal of Hydrodynamics（Ser.A），2013，28（4）：482-485（in Chinese）.

［5］李井煜，盧曉平，趙鵬偉.直接邊界元法解勢流速度場問題［J］.中國艦船研究，2015，10（1）：68-75.LI J Y，LU X P，ZHAO P W.Direct boundary element method for the problem of potential flow velocity field［J］.Chinese Journal of Ship Research，2015，10（1）：68-75（in Chinese）.

［6］韓玉超，盧曉平，王中.無限區(qū)域二維勢流直接邊界元法精度分析［J］.中國艦船研究，2015，10（4）：39-45，54.HAN Y C，LU X P，WANG Z.Precision analysis of the two-dimensional potential flow problem in an infi?nite region with the direct boundary element method［J］.Chinese Journal of Ship Research，2015，10（4）：39-45，54（in Chinese）.

［7］ 劉應中.船舶興波阻力理論［M］.北京：國防工業(yè)出版社，2003.

［8］王獻孚，周樹信，陳澤梁.計算船舶流體力學［M］.上海：上海交通大學出版社，1992.

［9］程天柱，石仲堃.興波阻力理論及其在船型設(shè)計中的應用［M］.武漢：華中工學院出版社，1987.

［10］熊鷹，盧曉平.艦船快速性計算圖冊［M］.武漢：海軍工程學院，1992.

Wave resistance calculation method combining Green functions based on Rankine and Kelvin source

LI Jingyu，LU Xiaoping
Department of Naval Architecture Engineering，Naval University of Engineering，Wuhan 430033，China

［Ojectives］At present,the Boundary Element Method（BEM） of wave-making resistance mostly uses a model in which the velocity distribution near the hull is solved first，and the pressure integral is then calculated using the Bernoulli equation.However，the process of this model of wave-making resistance is complex and has low accuracy.［Methods］To address this problem，the present paper deduces a compound method for the quick calculation of ship wave resistance using the Rankine source Green function to solve the hull surface's source density，and combining the Lagally theorem concerning source point force calculation based on the Kelvin source Green function so as to solve the wave resistance.A case for the Wigley model is given.［Results］The results show that in contrast to the thin ship method of the linear wave resistance theorem，this method has higher precision，and in contrast to the method which completely uses the Kelvin source Green function，this method has better computational efficiency.［Conclusions］In general，the algorithm in this paper provides a compromise between precision and efficiency in wave-making resistance calculation.

linear wave resistance；Green function；Rankine source；Kelvin source；Boundary Element Method（BEM）

U661.33

A

10.3969/j.issn.1673-3185.2017.06.001

http://kns.cnki.net/kcms/detail/42.1755.TJ.20171128.1110.022.html期刊網(wǎng)址：www.ship-research.com

李井煜，盧曉平.基于Rankine源和Kelvin源格林函數(shù)求解興波阻力的復合算法［J］.中國艦船研究，2017，12（6）：1-5.

LI J Y，LU X P.Wave resistance calculation method combining Green functions based on Rankine and Kelvin source［J］.Chinese Journal of Ship Research，2017，12（6）：1-5.

2017-04-06 < class="emphasis_bold"> 網(wǎng)絡(luò)出版時間：

時間：2017-11-28 11:10

國家部委基金資助項目

李井煜，男，1990年生，博士生，研究方向：艦船流體動力性能。

E-mail：935228691@qq.com

盧曉平（通信作者），男，1957年生，博士，教授。研究方向：艦船流體動力性能。

E-mail：luxiaoping100@163.com