基于Rankine源和Kelvin源格林函数求解兴波阻力的复合算法

李井煜，卢晓平430033海军工程大学 舰船工程系，湖北 武汉

摘 要：［目的］运用边界元法计算船舶兴波阻力基本上是先求解船体附近的速度分布，然后采用伯努利方程进行压力积分，其计算过程复杂，且误差非常大。［方法］提出一种可快速计算船舶兴波阻力的复合算法，利用Rankine Lagally Kelvin源格林函数求解船体表面源强，结合 定理进行受力计算，并基于 源格林函数求解船舶兴Wigley波阻力。运用该算法对 船的兴波阻力进行计算。［结果］计算结果表明，所用算法相较于运用线性兴波阻Kelvin力中的薄船理论得到的结果精度更高，而且与完全使用 源格林函数的算法相比效率也更高。［结论］所用算法可在计算兴波阻力时作为精度与效率之间的一种折中方法。关键词：线性兴波阻力；格林函数；Rankine源；Kelvin源；边界元法中图分类号：U661.33 文献标志码：A DOI：10.3969/j.issn.1673-3185.2017.06.001

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

Abstract：［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. Key words：linear wave resistance；Green function；Rankine source；Kelvin source；Boundary Element Method（BEM）

0引言

Boundary Element Method，BEM）边界元法（已被广泛应用于船舶兴波阻力的数值计算中，该方法可分为间接边界元法（分布源方法，或称间接法）和直接边界元法（源—偶混合分布法，或称直接法）。间接法可直接、方便地求解得到流场中的速度，且能推导并用于流体力学中常用到的源汇强度概念，因此在各种商用软件得到了广泛应用［1］。这两种方法的计算效率相差无几，但对计算精度而言，直接法略高于间接法。近年来，间接法的发展遇到了一些瓶颈。例如，文献［2］提到间接法在求解非光滑边界处的切向诱导速度时计算精度极差，而对于非线性问题，其计算精度很低，甚至导致发散，故近年来有学者提出发展直接法来提高计算精度［3-6］。无论是使用直接法还是间接法，都涉及到格林函数（基本解）的选择问题。Kelvin源格林函数可满足线性自由面兴波条件，在基于此格林函数的船舶兴波阻力计算中，可以不需要在流场自由液面布置源汇。Rankine源格林函数形式简单，但在基于此格林函数的船舶兴波阻力计算中，需要在所有物面布置源汇，故仅能近似地在有限区域求解兴波问题［7-8］。实际上，若仅关注船舶航行时的兴波阻力，计算求解出物面上的源汇分布是关键。Lagally本文将基于 定理和线性兴波阻力理Rankine论，提出一种新复合算法。该算法运用 源Kelvin格林函数计算船体表面的源强密度，采用源格林函数计算兴波阻力。与线性兴波理论的薄船理论相比，复合算法的计算精度较高，而与完全Kelvin使用 源格林函数的方法相比，复合算法的计算效率更高。本文所用算法可作为兴波阻力计算时精度与效率之间的一种折中方法。

1 数值计算

本文假定船舶在无限深广的自由水面上，以航速U在x轴负向做匀速直线运动。流场满足均匀、不可压缩和无粘、无旋流动条件，且船体兴波属于微振幅波。取固结在船上的随船坐标系，如1图 所示。 控制方程为不可压缩、流动的连续性方程，边1）~界条件为物面不可穿透条件，表达式如式（式（3）所示。 Ñ2Φ =0 （1） Φ(x y z) = Ux + φ(x y z) （2） ¶Φ =0 （3） ¶n式中： Ñ为哈密顿算子；Φ 为流场总速度势；φ为

船体扰动速度势；n为船体法线方向向量。由拉盖尔定理可知，流体对物体的作用力F可认为是多个点源i对其产生作用力的总和，即F = 4πρå miqi （4） i式中：ρ为流体密度；mi 为点源的源强；qi 为点源处的速度。针对上述船舶兴波问题，本文认为流场由物面分布的点源诱导产生，各单元的源强 σi dA（σ i为各单元的源强密度，A为面积）将受到沿航行方向的阻力影响，故船体兴波阻力 R 被认为是作用w于整个分布的源汇上的作用力积分之和，即R =- 4πρ  σiVdA =- 4πρ  σi ¶Φ （5） dA ¶x w A A式中，V为单元内流场沿x轴方向的速度分量。根据势流理论，船体扰动速度势 φ由船体表

面分布的源汇诱导产生，为

（6）

φ =-  σi dA × Gs A式中，G 为单位源强产生的速度势，或称格林函s数，下标 s表示某面元上的格林函数。在兴波问2 Kelvin题中，常见的 种格林函数包括 源格林函Rankine Kelvin数和 源格林函数。 源格林函数满足线性自由面条件，表达式如式（7）所示。-ππ G = 1 - 1 - 2K0 2 sec2 θdθ× s r r π 1 2 2 cos[K (x - ξ )cos θ]cos[K ( y - η)sin θ]dK+ 0¥ K (z + ξ ) e K - K0 sec2θ

π  K sec2 θ(z + ξ) 2K 22 0 e sin[k 0(x - ξ )sec θ]× 0 -π cos[k 0( y - η)sec2 θ sin θ]× （7） cos[K 0( y - η)sec2 θ sin θ]sec2 θdθ式中：(x y z) 为场点坐标；(ξ η ζ )为源点坐标； r1 ，r ，K ，K ，θ 均为计算的过程量，其中： 2 0 g K = 0 2 U r1 = [(x - ξ )2 + ( y - η)2 + (z - ζ )2]1 2 r = [(x - ξ )2 + ( y - η)2 + (z + ζ )2]1 2 2 Kelvin式中， g 为过程量。相比于 源格林函数，

Rankine源格林函数不满足自由面兴波条件等有限域内的边界条件，但可简化为： G = 1 s r

r= (x - ξ )2 + ( y - η)2 + (z - ζ )2另一方面，将式（2）和式（6）代 式（5），可入 得到如下兴波阻力表达式：

R =- 4πρU  σi (ξ η ζ )dA(ξ η ζ )w A

¶φ 8 4πρ  σi (ξ η ζ) dA(ξ η ζ) （ ） ¶x A根据文献［9］，式（8）的第1 0，故兴波阻项为力的表达式可进一步简化为¶φ （9） R =- 4πρ  σi (ξ η ζ ) dA(ξ η ζ) ¶x w A将式（7）代入式（9）可 4得到 个积分项。因积分区域是关于点源和场点对称，而根据文献［9］的4分析，只有第 项对兴波阻力起作用，故其表达式为：

π

10  R = 8πρK 02 2π (P + Q 2)sec3θ dθ （ ） 2 w - 2 Q}= cos}[K0 P  σi sec2 θ(x cos θ + y sin θ)]× sin A 11 exp(K z sec2 θ)dA （ ） 0式中，P 函数和 Q 函数为计算过程量，故只需求

解得到船体物面上的源强密度分布，即可求解兴波阻力。根据线性兴波理论的薄船理论，船体表面的源强密度直接与船体形状相关。但是，根据Neumann-Kelvin Kelvin理论或其他使用 源格林函数的非线性兴波阻力理论，计算源强分布的过程十分复杂。从整体上看，后者的精度要高于前者。在物理直观上，前者的精度可以通过简化的Rankine线性兴波理论得到提高，即通过使用 源格林函数求解无限域中的叠模来确定各单元源强密度的分布。综上所述，本文拟提出一种求解兴波阻力的Rankine快速复合算法。该算法使用 源格林函数求解船体表面各单元的源强密度，结合基于Kelvin源格林函数的式（10）和式（11）求解兴波阻力。本文所用算法实际上可在计算兴波阻力时作为精度与效率之间的一种折中方法。其中，使用Rankine源格林函数求解源强密度可以使用直接边界元法和间接边界元法求解。根据直接边界元法，求解源强密度的公式为¶ ¶φ(N) C(M )φ(M ) = 1 AN φ(N ) (1) - 1 dA 4π ¶nN r r ¶nN N

（12）

2 数值验证 2.1 Wigley船型兴波阻力计算分析

Taylor Wigley图谱［10］适在 用范围内选取一种1船型作为第 个算例用以分析和验证本文所用复1合算法的有效性。该船型参数如表 所示。 本文对船体的网格划分使用三角形网格。船40×4×2体划分为 个网格单元。考虑到上叠模部640 2分的网格单元，单元总数为 个，单元信息如图3 1 Wigley所示。图 为第 个算例中 船型的兴波阻力计算结果比较。 3 Rankine由图 可以看出，结合了叠模技术、源格林函数法、Lagally定理和直接边界元法（复合Michell算法）的计算结果和 薄船线性兴波阻力理

Taylor论的计算结果相比更接近于 图谱计算值， Fr>0.35尤其是在 的高速阶段。同时，还发现在低速阶段存在直接边界元法计算不准确的现象。下文将对该部分的误差进行详细分析。2 Wigley第 个算例为另一种 船型。该模型曾605 =5m，在 所进行过阻力试验，模型参数如下：L B=0.4 m，T=0.178 m。通过Prohaska方法确定模型形状因子，进而可近似得到兴波阻力的试验值。4 2 Wigley图 所示为第 个算例中 船型的兴波阻力试验结果与计算结果的比较。 4由图 可以看出，基于复合算法的直接边界Michell元法和间接边界元法的计算精度要高于薄船线性兴波阻力理论的精度，尤其是在中、高速段。对比直接边界元法和间接边界元法，发现直接边界元法的计算精度比间接边界元法的稍高。Fr<0.2但正如前文所述，本文所用算法在 的低速2段计算误差较大。表 所示为各算法精度的定量分析结果。

2.2 低速段误差分析

式（11）可知，本文所用复合算法的误差分析 主要来自如下项，为cos}[K0 sec2 θ(x cos θ + y sin θ)]× exp(K z sec2 θ) sin 0当速度较低时， K 很大，这将导致此项沿 x 0

轴方向产生高频振荡，指数项在很大程度上增加5所示为某一水平线上式（11）中了振荡幅值。图的 Q 函数随 x轴振荡的趋势，其中纵坐标表示π J =  2π sin[K0 sec2θ(x cos θ + y sin θ)]· - 2 exp(K z sec2 θ)dθ 0由 Q 函数的主要积分部分交换积分次序消去 Q

后组成。 针对这部分误差，一方面可以通过增大 x 轴方向的网格划分密度加以改进，另一方面需要进一步研究振荡函数的积分算法来估算误差。

3结论

本文提出了一种求解船舶兴波阻力的复合算法，介绍了复合算法的原理，并验证了该方法。通

过数值算例，得到如下结论： 1）总体上，无论是使用基于直接边界元法还是使用间接边界元法的复合算法来求解船舶兴波Michell阻力，与 积分公式等传统的数值计算方法NeumannKelvin相比都具有较高的精度，同时与基于理论的求解方法相比复合算法的效率更高，故可在兴波阻力计算时作为精度与效率之间的一种折中方法。2 ）直接边界元法的求解精度略高于间接边界元法。3）在低速段，采用复合算法计算的精度较低，主要是因为该方法中的 P 函数和 Q 函数包含了振荡项，故处理振荡函数的方法仍需进一步研究。

参考文献：

1］ . ［ 戴遗山，段文洋 船舶在波浪中运动的势流理论［M］. 北京：国防工业出版社，2008. 2］ 徐刚. ［ 不规则波中浮体二阶水动力时域数值模拟［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］ 陈纪康，段文洋，朱鑫. 三维泰勒展开边界元方法 及其数值验证［ 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］ 韩玉超，卢晓平，王中. 无限区域二维势流直接边析［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］. ［ 北京：国防工业出版社，2003. 8］ 王献孚，周树信，陈泽梁. 计算船舶流体力学［M］. ［ 上海：上海交通大学出版社，1992.

9］ 程天柱，石仲堃. ［ 兴波阻力理论及其在船型设计中的应用［M］. 武汉：华中工学院出版社，1987. 10］ 熊鹰，卢晓平. M］. ［ 舰船快速性计算图册［ 武汉：

海军工程学院，1992.

图3 1 Wigley第 个算例中 船型兴波阻力计算结果比较Fig.3 Results comparison of the Wigley ship form in case 1

图4 2 Wigley第 个算例中 船型兴波阻力计算结果比较Fig.4 Results comparison of the Wigley ship form in case 2

2 Wigle图 船型网格Fig.2 Panel mesh of the Wigley ship form

图1 坐标系Fig.1 Coordinate system