刘杨，龚家烨，刘家瑞，陈京普

(1.中国船舶科学研究中心上海分部，上海 200011；2.哈尔滨工程大学，哈尔滨 150001；3.中国船舶科学研究中心，江苏 无锡 214000)



三体船舶阻力及航态预报方法

刘杨1，龚家烨2，刘家瑞3，陈京普1

(1.中国船舶科学研究中心上海分部，上海 200011；2.哈尔滨工程大学，哈尔滨 150001；3.中国船舶科学研究中心，江苏 无锡 214000)

针对三体船型的自身特点，基于势流理论，在Dawson方法的基础上对数值处理进行改进，并运用于三体船舶的阻力预报。通过近水面网格的快速划分与迭代计算，实现对航行姿态的预报。根据计算所得的船侧兴波，对不同航行姿态下的湿表面积进行计算，提高总阻力预报的精度。以三体Wigley船型为例，将理论计算结果和试验值及商业软件参考值进行对比。证明了理论预报方法的可行性与有效性，改进后的数值方法结果准确，可很好地避免自由面的震荡。

三体船；Dawson方法；兴波阻力；航行姿态；湿表面积

目前对于三体船型阻力的预报方法主要分为粘流理论以及势流2种。粘流理论主要通过一些商用软件(CFX、Fluent等)直接对总阻力进行计算，但是并未对兴波阻力的产生机理进行探讨，不利于针对兴波阻力的优化，同时在多种船型、布局优化时需要大量的计算时间。势流方面主要运用有Noblesse细长理论、薄船理论、Michel积分等方法，运用势流理论可以达到快速预报的目的，但传统方法并不能同时针对兴波阻力、航行姿态、近场兴波等多方面同时进行计算。为此，考虑在现有三体船阻力预报研究[1-5]基础上，对传统Dawson方法[6]进行改进，探讨考虑航行姿态变化的兴波阻力预报方法，并通过对湿表面积的计算，提高总阻力预报的精度。通过对数值计算方法[7-12]的对比分析，根据三体船型自身特点对网格的处理进行改进，采用控制点上移等方法实现在较高航速下得到光顺的自由面波形。运用对近水面网格实时迭代重新划分与表面压力回归分析的方法提高了船体受力与湿表面积计算的精度。

1 数学模型及数值方法

1.1 坐标系选取

采用随船坐标系，见图1。坐标系原点取船舶中心，z轴垂直向上，x轴由船首指向船尾，y轴指向右舷为正，xy平面与静水面重合。

图1 随船坐标系

1.2 基本方程

假设流体理想，流动定常，无限水深，船舶兴波为微幅波。流场速度势ψ分解为无自由面影响的重叠模速度势φ和考虑自由液面作用的扰动速度势φ，即

(1)

(2)

基于上述假设，φ≪φ。

速度势满足以下定解条件

1)流场中。

(3)

2)物面上。

(4)

3)自由面上。

(5)

式中，n为物面法向量。

同时，速度势满足无穷远处趋于零的远方辐射条件。

自由面上的φ关于静水面对称，则有

(6)

先后在在z=ξ和z=0上将速度势进行泰勒展开，忽略φ的非线性项及高阶项，并结合流线假设▽φ▽φ=φlφl，可得简化的自由面边界条件

(7)

式中，l为重叠模流线方向。

1.3 数值离散

将船体表面及自由面进行数值离散，取船体网格数为NB，自由面网格数Nf，以四边形网格为主，部分不规则区域采用三角形网格，自由面网格划分采用贴水线网格，对于片体间网格采用均匀大小网格，船体外侧横向及纵向网格则按一定比例增大。

对计算边界进行数值离散后，任意控制点p处重叠模流场速度势势表示为

(8)

式中，SB表示物面；σB0为物面控制点处源强；rpq为任意源点q到控制点p的距离。

由船体与自由面相互作用产生的扰动速度势表示为

(9)

对于自由面条件的求解，是该方法求解过程中的关键一步，尤其对于高性能船舶，自由面方程中差分格式尤为重要。本方法采用改进后的差分格式，将流线差分坐标变换为关于横纵向网格的差分。因此自由面条件中流线速度势的二阶偏导数φll，φll表示为

(10)

纵向差分统一采用四点向前差分，对于自由面前缘单元，采用虚单元的方法进行计算；横向差分计算时，对于片体与主体间网格采用中心差分，而对于船体外侧的网格则采用三点向前差分进行计算。具体差分格式如下。

(11)

式中，(∂L/∂x)，(∂T/∂x)通过坐标转换得，转换矩阵如下

(12)

求解后得到物面源强σB1和源强σF。

1.4 深沉、纵倾的计算

船体表面控制点压力系数为

(13)

式中，zi为控制点到静水面间距离。

船体受到的升力和纵倾力矩可表示为

(14)

(15)

式中，xi,zi为控制点坐标；xg,zg为船舶的重心坐标；nxi,nzi为单元法线方向。

根据纵倾力矩及升力确定纵倾、升沉，快速划分近水面局部网格，进行快速迭代计算，忽略大部分网格重新计算的内存浪费，给定时间步长和收敛条件值，迭代得到定常状态下的深沉纵、倾值。

1.5 数值处理

根据兴波波长的不同采取不同的自由面网格大小。随着弗汝德数增加网格数减小，但至少保证船长方向有25个网格，在满足精度的基础上，加快计算速度。经过分析计算，与船体相连的第一层网格的宽度取0.15Lwl，片体间网格采用均布网格，并进行一定程度的加密；船体外侧的网格，横向、纵向增长比例rx、ry均取1.03。

由于三体船舶片体之间相互兴波干扰较大，单元对于坐落于单元中心的控制点的切向诱导速度为零，在自由面条件向前差分过程中容易产生自由面震荡。因此，在处理自由面控制点分布时，将所有控制点根据单元网格长度，进行向前移动，并且根据不同航速及区域的单元长度进行控制点上移，以消除震荡。

采用船体表面压力积分的方式求解船舶阻力，并基于粘流理论计算粘性阻力，得到船舶总阻力。在表面压力积分和湿表面积计算时，结合直舷假设，根据船侧兴波对船体湿表面进行重新离散与细化，并对垂向一列单元的压强进行线性回归，进而通过垂向积分得到船体的兴波阻力及湿表面积。这样既可以提高阻力和湿表面积的计算精度，也无需增加计算过程中的物面网格数量，节省了计算时间。

2 算例及结果分析

2.1 排水式船体数值计算及验证

采用三体Wigley船型，对不同的数值处理进行探讨，从而得出较优的参数设置。Wigley船型其数学表达式为

(16)

式中，L为船长；B为船宽；D为吃水，主体L/B=12.5,L/D=28.1；侧体L1/B1=12.5,L1/D1=28.1。船体基础网格数660个，自由面2 376个(Fr=0.25)。网格数随弗劳德数变化联动调整，网格划分如图2、3。

图2 Wigley船体网格划分

图3 自由面网格(Fr=0.25)

探讨了在三体船舶阻力及航态数值模拟时，控制点升高高度对数值计算结果影响。取Dx为单元长度，计算方案如表1。

表1 控制点上移方案

将方案1、方案2、方案3计算结果进行对比，阻力、纵倾及升沉的对比结果和波形模拟见图4～6。方案3的计算所得波形见图7。

图4 计算兴波阻力系数同试验剩余阻力系数比较(方案1、2、3)

图5 纵倾曲线(方案1、2、3)

图6 升沉曲线(方案1、2、3)

通过图4的对比结果可以看出，理论计算所得兴波阻力系数与试验所得剩余阻力系数趋势及峰谷值吻合良好，而同图5、图6中深沉纵倾的对比可以明显看出控制点上移比例为Dx/150的方案3计算精度最高。对比结果表明，控制点上移方案3可以在消除自由面震荡的前提下保证计算精度。图5表明，方案3计算所得的自由面波形也与实际符合良好的，并且光顺。

由上述对比可知，方案3优于方案1、方案2，因此继续将方案3同约束模型方案4及试验和商业软件参考值结果对比，如图8～10，方案3计算所得的湿表面积变化见图11。

图8 计算兴波阻力系数同试验剩余阻力系数比较(方案3、4)

图9 升沉曲线(方案3、4)

图10 纵倾曲线(方案3、4)

图11 湿表面积变化曲线(方案3)

图8表明，由于忽略了升沉、纵倾的影响，方案4的计算结果明显小于程序和商业软件结果。图9、10表明，总体上本方法计算结果并不亚于商业软件的非线性结果，且本方法可以再任意航速下得到光滑而完整的自由面兴波，不存在自由面发散的问题。根据图9的深沉结果，本方法计算结果优于参考结果，与试验值也较接近，这对于粘性类流体力的估算是很有价值的。但是在图10中可以看出，在弗劳德数0.5左右，理论预报方法所得的纵倾角相对较大，其主要原因时忽略了自由面高阶项的作用，对自由面兴波造成一定的影响，也影响了压力分布的精度。图11的湿表面积变化曲线则可以看出，随着航速的增大，三体船舶的湿表面积变化比例较大，因此在总阻力估算时，准确计算湿表面的影响是十分必要的。

2.2 不同片体布局的数值模拟

运用上述方法及优化后的数值参数，对相同三体Wigley船型不同布局下的三体船舶进行数值模拟和计算。三体船片体布局示意图见图12，具体方案见表2。

图12 片体布局

表2 片体布局方案

表2中，L为三体船主体的水线长度，a为片体中心纵向位置与船长的关系。片体中心纵向位置与主体相同时a为零，片体向船艉移动为正。

首先对方案1与方案3的计算结果与试验值，结果见图13～16。

图13 计算兴波阻力系数计算值与试验剩余阻力系数比较(方案1)

图14 兴波阻力系数(方案3)

图15 升沉曲线(方案1)

图16 升沉曲线(方案3)

从图13～16可以看出，计算兴波阻力系数与试验所得剩余阻力系数吻合良好，升沉值也吻合较好，可以看出本方法对于不同片体布局的三体船舶的阻力及航态变化可以准确进行预报，对于片体位置变化所造成的影响是比较敏感的，且在较大的弗汝德数范围内能保证精度，这对于高速三体船舶的船型优化及片体布局优化是有较大实际意义的。

方案1～方案4中不同布局的三体船舶计算值对比见图17～20。

图17 计算兴波阻力系数同试验剩余阻力系数比较(方案1～方案4)

图18 升沉对比

图19 纵倾对比

图20 湿表面积变化对比

对比结果中可以明显看出，片体后移的方案3、方案4随着弗汝德数的增加阻力性能优于方案1、方案2，而且在深沉和湿表面积上也小于方案1、方案2，这对于减小粘性阻力也是有益的，计算结果与实际情况是相符的。理论方法可以很好地对不同片体布局的三体船舶进行预报。

方案1、方案4的兴波波形图见图21、22。

图21 方案1波形图(Fr=0.496)

图22 方案4波形图(Fr=0.496)

图21、22中所反映出的兴波情况与计算结果时相符，在Fr=0.496时，前述计算结果方案4的兴波阻力系数也是大大小于方案1的，而且即使在弗汝德数较大的情况下，本方法所用的数值处理仍然可以得到较良好的波形，未发生自由面震荡的情况。

3 结论

1)通过数值计算方法改善，可以在非全非线性自由面条件下改善数值预报结果，克服很多方法在自由面模拟的发散问题，获得三体船阻力及航行姿态较好的预报结果。

2)本方法可在三体船设计阶段为片体布局及船型初步优化提供技术支撑。

3)由于三体船侧体与主体与片体之间存在较强的兴波干扰问题，因此传统的模型很难准确地对非对称流场进行模拟，在后续的工作中，应当考虑如何采用有升力模型对非对称有升力流场进行数值模拟，进而得到更准确的流场分布及数值结果。

[1] 何术龙,李百齐,程明道,等.三体船船型分析及兴波干扰的模型试验研究[J].水动力学研究与进展:A辑,2006,21(1):122-129.

[2] 吴广怀,吴培德,蒋耀军,等.基于兴波阻力的三体船片体位置快速优化方法[J].船舶力学,2005,9(4):1-8.

[3] 郦云,卢晓平.三体船阻力性能研究[J].船舶力学,2007,11(2):191-198.

[4] 蔡新功,王平,谢小敏.三体船方案优化布局的阻力计算与试验研究[J].水动力学研究与进展:A辑,2007,22(2):202-207.

[5] 陈京普，朱德祥，何术龙.双体船/三体船兴波阻力数值预报方法研究[J].船舶力学,2006,10(2)：23-29.

[6] DAWSON C W. A practical computer method for solving ship-wave problems[C]. The second International Conference on Ship Hydrodynamics.Beckly,1977:30-38.

[7] RAVEN H C. Variations on a theme by Dawson[C]. 17th Symposium on Naval Hydrodynamics. The Hague NetheRlands,1988:151-171.

[8] Md Shahjada Tarafder. Numerical calculation of free-surf-ace potential flow around a ship using the modified Rankine source method[J]. Ocean engineering,2008,35:536-544

[9] RAVEN H C. A solution method for the nonlinear ship wave resistance problem[D]. Delft, Netherlands, Delft University of technology,1996.

[10] 贺五洲，程明道，俞汉祥.方艉舰船及加装尾板的兴波阻力计算[J].中国造船,1998(2)：13-18.

[11] 刘应中.船舶兴波阻力理论[M].北京：国防工业出版社，2003.

[12] 卢晓平，王中，孙永华，等.Rankine源Dawson型方法求解三体船兴波阻力[J].华中科技大学学报(自然科学版),2008,36(11):103-107.

Research on Prediction Method for Resistance and Hull Gesture of Trimaran

LIU Yang1, GONG Jia-ye2, LIU Jia-rui3, CHEN Jing-pu1

(1.China Ship Scientific Research Center, Shanghai 200011, China;2.Harbin Engineering University, Harbin 150001, China;3.China Ship Scientific Research Center, Wuxi Jiangsu 214082, China)

In light of characteristics of trimaran, some improvements were made on the numerical treatment of Dawson’s method based on potential flow theory to predict the resistance. By quick re-mesh of the grids near the water-plane in iterative calculation, the prediction of hull gesture was realized. On the basis of the calculated wave profile along the ship, the wetted surface area of trimaran in different hull gestures was calculated, and the prediction accuracy of total resistance was improved. With the example of tri-Wigley hull trimaran, theoretical calculation result was compared with the experimental and numerical results. The improved numerical method was proved to be feasible and effective, which can prevent the oscillation of the free surface.

trimaran; Dawson’s method; wave making resistance; hull gesture; wetted surface

10.3963/j.issn.1671-7953.2017.04.001

2016-09-12

刘杨(1990—)，男，硕士，助理工程师

研究方向：船舶阻力与推进技术

U661.3

A

1671-7953(2017)04-0001-06

修回日期：2016-10-11