李井煜，盧曉平

海軍工程大學(xué)艦船工程系，湖北武漢430033

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

李井煜，盧曉平

海軍工程大學(xué)艦船工程系，湖北武漢430033

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

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

0 引 言

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

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

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

1 數(shù)值計(jì)算

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

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

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

式中：?為哈密頓算子；Φ為流場(chǎng)總速度勢(shì)；φ為船體擾動(dòng)速度勢(shì)；n為船體法線方向向量。

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

式中：ρ為流體密度；mi為點(diǎn)源的源強(qiáng)；qi為點(diǎn)源處的速度。

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

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

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

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

式中：(x，y，z)為場(chǎng)點(diǎn)坐標(biāo)；(ξ，η，ζ)為源點(diǎn)坐標(biāo)；r1，r2，K，K0，θ均為計(jì)算的過(guò)程量，其中：

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

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

根據(jù)文獻(xiàn)［9］，式（8）的第1項(xiàng)為0，故興波阻力的表達(dá)式可進(jìn)一步簡(jiǎn)化為

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

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

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

2 數(shù)值驗(yàn)證

2.1 Wigley船型興波阻力計(jì)算分析

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

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

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

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

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

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

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

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

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

2.2 低速段誤差分析

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

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

主要來(lái)自如下項(xiàng)，為

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

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

圖5 不同速度段的振蕩項(xiàng)QFig.5 The oscillating termQin different speed range

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

3 結(jié) 論

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

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

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

3）在低速段，采用復(fù)合算法計(jì)算的精度較低，主要是因?yàn)樵摲椒ㄖ械腜函數(shù)和Q函數(shù)包含了振蕩項(xiàng)，故處理振蕩函數(shù)的方法仍需進(jìn)一步研究。

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

［2］徐剛.不規(guī)則波中浮體二階水動(dòng)力時(shí)域數(shù)值模擬［D］.哈爾濱：哈爾濱工程大學(xué)，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］陳紀(jì)康，段文洋，朱鑫.三維泰勒展開(kāi)邊界元方法及其數(shù)值驗(yàn)證［J］.水動(dòng)力學(xué)研究與進(jìn)展（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］李井煜，盧曉平，趙鵬偉.直接邊界元法解勢(shì)流速度場(chǎng)問(wèn)題［J］.中國(guó)艦船研究，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］韓玉超，盧曉平，王中.無(wú)限區(qū)域二維勢(shì)流直接邊界元法精度分析［J］.中國(guó)艦船研究，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］ 劉應(yīng)中.船舶興波阻力理論［M］.北京：國(guó)防工業(yè)出版社，2003.

［8］王獻(xiàn)孚，周樹(shù)信，陳澤梁.計(jì)算船舶流體力學(xué)［M］.上海：上海交通大學(xué)出版社，1992.

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

［10］熊鷹，盧曉平.艦船快速性計(jì)算圖冊(cè)［M］.武漢：海軍工程學(xué)院，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ù)求解興波阻力的復(fù)合算法［J］.中國(guó)艦船研究，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ò)出版時(shí)間：

時(shí)間：2017-11-28 11:10

國(guó)家部委基金資助項(xiàng)目

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

E-mail：935228691@qq.com

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

E-mail：luxiaoping100@163.com