基于RANS法和边界层理论预报三维船体阻力

胡俊明，李铁骊，林 焰，徐雪锋，徐利刚

（1.大连理工大学 船舶工程学院，辽宁 大连 116024；2.海军装备部直属工作部，北京 100841；3.海军驻无锡地区军事代表室，江苏 无锡 214061）

0 引 言

船舶阻力问题是船舶行业重要的研究课题，如何快速精确预报船体阻力成为问题研究的关键。模型试验法在船舶阻力预报中一直占据着主导地位，其预报结果精确可靠，实用性强，但其经济性较差且无法完全模拟复杂的实际海况，存在尺度效应[1]。理论计算方法由来已久，由于船体复杂的三维曲面，水的粘性和自由面的非线性，使理论求解相当困难，目前求解摩擦阻力沿用傅汝德的相当平板假定，而兴波阻力理论由1898年Michell薄船理论的出现，奠定线性兴波理论的基础，之后的慢船理论Dawson法的产生，借助摄动理论使线性兴波向非线性兴波发展成为一种可能[2-4]。随着技术进步，计算机性能更新和完善，计算流体力学（CFD）的兴起，使数值求解复杂曲面粘性自由面绕流场成为现实，计及波浪破碎、涡分离、抨击等非线性现象，且计算结果与船模试验结果一致，真实可靠，经济性好，可用范围广，成为研究船舶阻力性能的重要手段[5-7]。

本文应用二因次RANS法研究船舶阻力，基于粘性流理论N-S方程求解自由面数值绕流场，重点研究了船体自由面绕流模型的网格划分方法，应用ICEM网格处理软件处理内外域交接面，实现交界面的网格重构，以提高内外域不同网格类型之间的数据交换效率，缩短计算时间，同时应用UDF控制压强出口边界处压强随水深变化以提高数值计算稳定性，提高计算效率。数值计算中基于边界层理论考虑边界层第一层网格高度对数值计算结果的影响，重在研究其对船体各阻力分量的影响。基于此方法所得数值结果与实验流体力学法（EFD）做比较分析，通过其计算精度、船体阻力各分量值和兴波波形的比较得到其变化规律，结果表明，此方法预报船体阻力其数值精度满足工程需要，具有较强的实用性。

1 船型特征及模型试验

文中以总长为33.2 m的渔船为实船原型，进行船模试验，实船与模型的缩尺比为10，且实船与模型在压载状态下主尺度如表1所示，表中LOA为总长，LWL为设计水线长，B为船宽，T为设计吃水，CB为方形系数。应用三因次EFD法对船模进行低速流试验，根据傅汝德数Fn=0.1～0.2范围内的试验结果，采用普鲁哈斯卡法确定形状因子1+K的值，应用1957-ITTC平板摩擦阻力公式，计算船模平板摩擦阻力Rfm，根据测定的船模总阻力Rtm，计算船模剩余阻力Rrm=Rtm-Rfm、粘压阻力Rpvm=KRfm、兴波阻力及各阻力系数[8]。

表1 船型主尺度Tab.1 Principal dimensions of the real ship and the model

2 计算方法

2.1 控制方程

RANS二因次法基于粘性流理论N-S方程求解自由面船体阻力，考虑自由面非线性且流体不可压缩，其流体的连续方程和雷诺平均N-S方程如下[9]：

式中：ρ是流体密度，t是时间，υ是流体运动粘性系数，p是流体压力，fi是质量力，ui为时均速度；ui′为脉动速度为雷诺应力项。雷诺应力对于精确求解丰满船型船体阻力至关重要，为使RANS方程封闭，应用RNG k-ε湍流模型，其输运方程为[10]：

表2 RNG k-ε湍流模型常数Tab.2 RNG k-εturbulence model constant

2.2 边界层理论

边界层影响船体阻力的数值精度，特别是边界层第一层网格高度的选择和确定，决定数值计算结果的偏差大小。本文采用RNG k-ε模型，应用壁面函数处理粘性底层和过渡层，用一维数学模型代替求解三维N-S方程，以降低计算资源的使用，同时选择合理边界层第一层网格高度使其网格节点落于对数律区域内，避免网格节点落在粘性底层而导致数值计算结果的偏差。其边界层第一层网格高度无量纲计算公式如下[11]：

式中：υ为流体的运动粘性系数，Δy为第一层网格高度，ρ为流体的密度，K和E为常数[9-10]。

3 计算模型建立与设置

3.1 模型及域选择

因船型对称且直航，受力对称，文中所建模型为半船模型，应用Auto CAD可视化绘图软件和ICEM网格处理软件联合建立自由面网格绕流计算模型。因自由面兴波对数值流场存在影响，为提高计算精度和效率应合理选择计算域的大小，采用Hirt和Nichols提出的多相流模型方法（VOF）追踪自由液面[12]。

应用水波理论中船行波波高ZA、波长λw与波速的关系选择计算域的大小，其绕流计算模型为长方体。由公式确定水线以上计算域高度大于船行波最大波高，基于计算船行波波长，由于自由面影响，波面在船体前缘隆起，取进流段长度大于λw，尾部采用压强出口以提高计算效率和数值结果的稳定，其长度应考虑首尾波系的影响，为实现其充分发展，文中去流段长度取为13.35 m，约大于6倍最大船行波波长。根据开尔文波系尽可能减小宽度区域边界对船行波的反射，由计算域长度和开尔文角选取合适的区域宽度，计算域深度大于λw时认为船行波对绕流场影响甚微[13]。由模型航速范围1.3～1.8 m/s，估算船行波波高范围0.043 1～0.082 5 m，波长范围1.081～2.073 m，其自由面网格绕流模型计算域尺寸划分如表3所示。

表3 模型计算域划分Tab.3 Size of the computational domain

3.2 网格划分

网格计算模型考虑了边界层和船体线型对流场的影响，文中借助网格处理软件ICEM实现内外域交界面网格重构以合成整体域，以提高内外域不同网格类型之间的数据交换效率，缩短计算时间，在船艏艉区域进行网格加密处理，网格尺寸取为0.003L，船体表面最大网格尺寸为0.01L，体网格由船体表面按照一定的增长率向外扩展[14]，内域增长率为1.1，在外域中，采用分块结构化网格，增长率为1.3，最大网格尺寸0.08L，考虑边界层影响，在船体表面划分8层边界层，取边界层为3 mm，增长率为1.1，在水线面附近进行网格加密，水面第一层高度为0.005L，增长率为1.2。其半船体模型如图1所示，整个流场网格和局部边界层网格划分如图2所示。

图1 船体模型Fig.1 Model ship

图2 网格划分图Fig.2 Mesh information

4 数值计算结果及分析

4.1 边界层对数值结果的影响

本文考虑边界层对数值计算结果的影响，对应实船设计航速，船模航速为1.5 m/s，傅汝德数Fn=0.274，边界层船体表面第一层网格高度分别取为2-6 mm，其数值计算结果如表4。由表中可以看出边界层第一层网格高度Δy充分影响数值结果精度，因其第一层网格高度节点未落入粘性底层，随着Δy值的减小，其摩擦阻力值略微减小，但影响甚小，其总阻力值增加，数值结果和试验值总体相差不大，同时基于粘性Δy值的变化改变了流线的形状，改变了船体压力分布，影响兴波阻力，直接影响剩余阻力值，对设计航速，当Δy=3 mm时，其数值计算结果与模型试验结果的误差最小，达到较好的吻合状态。此外，本文分别应用船体表面边界层第一层网格高度3 mm和4 mm，对船体模型进行各航速数值模拟计算，其计算结果详见表5。由表5中可以看出，第一层网格高度Δy=3 mm的数值结果的精度优于Δy=4 mm，随值增大，总阻力数值在高傅汝德数下其误差有增大趋势。结合表4和5，可以得到，对应任何一个航速模态，存在较佳的第一层网格高度Δy和值，使数值结果和试验结果达到较佳吻合状态，表中Rf为摩擦阻力、Rr为剩余阻力、Rt为模型试验总阻力，单位为N，error为总阻力相对误差。

表4 边界层对设计航速数值结果的影响Tab.4 Influence of boundary layer to numerical results

表5 边界层对各航速数值结果的影响Tab.5 Influence of boundary layer to numerical results under different speed

图3 摩擦和剩余阻力系数曲线图Fig.3 Curve of friction resistance coefficient and residual resistance coefficient

图4 粘压和兴波阻力系数曲线图Fig.4 Curve of viscous pressure resistance coefficient and wave resistance coefficient

4.2 阻力分析

应用边界层船体表面第一层网格高度3 mm，考虑带自由面粘性流计算，数值模拟不同傅汝德数下的船体阻力，所得结果与EFD法做分析比较，通过基于CFD二因次法所得摩擦阻力系数Cfm0与1957-ITTC平板摩擦阻力系数Cfm随雷诺数的变化曲线图3（a）中可以看出基于CFD带自由面粘性流得到的摩擦阻力较平板摩擦阻力略大，由图3（b）中可以看出基于CFD二因次法得到剩余阻力Crm0较EFD法剩余阻力Crm略小。基于三因次EFD法得到粘压阻力系数、兴波阻力系数的变化曲线如图4所示，Cpvm随雷诺数增加而减小，兴波阻力系数Cwm随傅汝德数增加而增大，与实际相符。

图5 不同航速下的自由面波形图Fig.5 Contour of wave pattern at various speed

4.3 波形及自由面分析

自由面波形反应了兴波特性，为了分析研究其变化规律，本文基于Tecplot技术处理RANS法所得自由面波形与EFD法做分析比较，取三个航速点（Fn=0.274、0.292和0.329），其自由面波形和等值线图如5所示。由图5中可以看出随着傅汝德数增加，船体首部波高等值线数值增大，波形等值线覆盖区域增加，船中及尾部区域规律变化愈发明显，船两侧受扰动区域增加，兴波区域向外扩散，波峰数值增加且范围扩大，船首尾肩部开尔文波系明显且波峰后移，其兴波波形完全符合开尔文波系特点，与EFD法兴波图形变化趋势吻合一致，这为船型优化提供了强有力的证据。

5 结 论

本文基于RANS法对总长为33.2 m带球鼻艏新型渔船进行船体阻力数值模拟，考虑粘性流自由面效应及边界层的影响，结合EFD方法对数值结果和兴波波形做分析比较，得到以下结论：

（1）应用ICEM处理网格绕流计算模型，实现内外域交接面网格重构，以提高内外域不同网格类型之间的数据交换效率，同时应用UDF控制压强出口边界处压强随水深变化以提高数值计算稳定性，缩短计算时间。

（2）基于边界层理论研究粘性流带自由面船体绕流问题，表明边界层对数值计算结果精度存在重要影响，边界层第一层网格高度节点未落入粘性底层时，在相同航速下，摩擦阻力值随边界层第一层网格高度的增加而略微增大，剩余阻力值却明显减少，总阻力在高傅汝德数下其数值误差有增大趋势，且边界层第一层网格高度相同时，随傅汝德数增加其船体总阻力误差增大。

（3）基于边界层数值结果分析得到对应任何一个航速模态，存在较佳的第一层网格高度Δy和y+值，使数值结果和EFD试验结果达到较佳的吻合状态。其RANS法数值结果得到的总阻力、摩擦阻力和剩余阻力系数与EFD法吻合较好，表明此方法用于船舶阻力性能预报的合理性、可行性和可靠性。

（4）自由面波形反应兴波特性，随傅汝德数增大船体首部波高越发明显，船体两侧受扰动区域增加，船中及尾部区域兴波向外扩散，船首尾肩部开尔文波系明显且波峰后移，其兴波图形变化趋势与EFD法吻合一致，为船型优化提供有力参考。

[1]倪崇本,朱仁传,缪国平,范佘明.一种基于CFD的船舶总阻力预报方法[J].水动力学研究与进展,2010,25(5):579-586.Ni Chongben,Zhu Renchuan,Miao Guoping,Fan Sheming.A method for ship resistance prediction based on CFD computation[J].Chinese Journal of Hydrodynamics,2010,25(5):579-585.

[2]Dawson C W.A practical computer method for solving ship-wave problems[C].In:Proceeding of Second International Conference on Numerical Ship Hydrodynamics,1977:30-38.

[3]Suzukik,Iokamorin.Studies on minimization of wave making resistance based on Rankine source method[J].Kansai Society of Naval Architecture in Japan,1999,185:9-19.

[4]Raven H C.A practical nonlinear method for calculating ship wave making and wave resistance[C]//Preprints 19th Symp.on Naval Hydrodynamics.Seoul,Korea,1992:60-75.

[5]刘应中,张怀新,李谊乐,缪国平.21世纪的船舶性能计算和RANS方程[J].船舶力学,2001,5(5):66-84.Liu Yingzhong,Zhang Huaixin,Li Yile,Miao Guoping.Calculation of the ship performance and solving of RANSequations in the 21st century[J].Journal of Ship Mechanics,2001,5(5):66-84.

[6]Choi J E,Min K S,Kim J H,Lee SB,Seo H W.Resistance and propulsion characteristics of various commercial ships based on CFD results[J].Ocean Engineering,2010,37:549-566.

[7]Ciortan C,Wanderley J,Guedes Soares C.Turbulent free-surface flow around a wigley hull using the slightly compressible flow formulation[J].Ocean Engineering,2007,34:1383-1392.

[8]盛振邦,刘应中.船舶原理[M].上海:上海交通大学出版社,2004.Shen Zhenbang,Liu Yingzhong.Principles of Ship.[M].Shanghai:SJTU Press,2004.

[9]Versteeg H K,Malalasekera W.An introduction to computation fluid dynamics:The Finite Volume Method[M].Wigley,New York,1995.

[10]Yakhot V,Orzag SA.Renormalization group analysis of turbulence:Basic theory[J].Scient Comput,1986(1):3-11.

[11]Schlichting H.Boundary Layer Theory[M].McGrawhill,New York,1979.

[12]Hirt C W,Nichols B D.Volume of Fluid(VOF)Method for the dynamics of free boundary[J].Journal of Computational Physics,1981(39):201-225.

[13]顾溟宇,王文华,王言英.基于FLUENT的大型集装箱船航速预报方法研究[J].船舶工程,2012,34(2):28-31.Gu Mingyu,Wang Wenhua,Wang Yanying.Studies of FLUENT-based speed prediction method for large containers[J].Ship Engineering,2012,34(2):28-31.

[14]朱德祥,张志荣,吴乘胜,赵 峰.船舶CFD不确定度分析及ITTC临时规程的初步应用[J].水动力学研究与进展,2007,22(3):363-370.Zhu Dexiang,Zhang Zhirong,Wu Chengsheng,Zhao Feng.Uncertainty analysis in ship CFD and the primary application of ITTCprocedures[J].Journal of Hydrodynamics,2007,22(3):363-370.