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

0 引 言

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

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

2016年底，在全市如火如荼、积极争创全国“双拥模范城”的进程中，在娄烦县对口扶贫的杏花岭区委统战部副部长王业发，被组织调任区双拥办主任。他工作履历比较丰富，公安干警出身，当过宗教局长，一身正气、两袖清风。生活中为人热情，说话做事雷厉风行，粗中有细。新的起点，新的征程。从此，王业发从零开始，以抓铁留痕的决心，踏石留印的力度，系统学习了解志愿兵的历史背景，查询退役军人的各种政策文件，吃准吃透文件精神，边学习边与老兵们接触，以人性化、情感化的工作方式化解老兵们的不安情绪，与老兵们无障碍沟通，深入持久地抓好“双拥”工作。

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

主尺度 LOA/m LWL/m B/m T/m CB实船模型33.2 3.32 30.5 3.05 9 0.9 3.25 0.325 0.649 0.649

2 计算方法

2.1 控制方程

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

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

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

应用水波理论中船行波波高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所示。

式中其 RNG k-ε湍流模型常数见表 2。

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

C1ε C2ε σk σε Cμ β η0 1.42 1.68 1.0 1.3 0.084 5 0.012 4.38

2.2 边界层理论

本文考虑边界层对数值计算结果的影响，对应实船设计航速，船模航速为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为总阻力相对误差。

努力改善教育支出的区域不平衡，加大对中西部教育的投入力度，缩小地区之间的差异；同时加大对农村的教育投入力度，缩小城乡差距．

3 计算模型建立与设置

3.1 模型及域选择

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

1.借助学分银行，实现继续教育与各类培训、学历教育的衔接与沟通。银行是个金融术语，指以开展信用经营为主，为资金的借贷双方提供中介平台的机构。“学分银行”是借鉴了金融领域的银行术语在教育领域的运用，指学习者把在学校教育、职业技能培训机构等多渠道获得学习成果进行学分化转换，以量化的学分记录在其学习档案中，并存在一个类似于银行数据交换中心的系统中，供学校、企业、政府等有关机构去查询其学习成果等。当学习者的学分积累够一定的学分量，达到某专业毕业所需的学分后，可以申请该专业的学历毕业证书，亦可通过积累一定量的学分申请某些职业技能证书。

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

设计水线长（m）水线以上高度（m）3.05 4.1 13.35 4 6 2.5进流段长度（m）去流段长度（m）水深（m）半宽（m）

3.2 网格划分

2.城镇住房供给调节机制：由计划机制向市场机制演进。改革开放40年来中国住房供给调节机制演进的最大特征是由单位分配到市场购买的结构性变迁。与此同时，国家持续加大保障房的供应。总的来看，尽管改革开放40年里中国城镇住房供给调节机制市场化改革成就斐然，但同样也带来了负面效应，表现为商品房供给被住房投资者和投机者引导，炒房时有发生，偏离了住房供给是为了满足民生需求的根本目标。这表明，新时代中国城镇住房供给调节机制改革应当是在市场调节机制基础上，加大政府调节，使城镇住房供给调节机制在一定范围内发挥作用，范围之外(譬如保障房、自建房、共建房)构建政府主导、社会参与和民众配合的城镇住房供给调节机制。

图1 船体模型 Fig.1 Model ship

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

4 数值计算结果及分析

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

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

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

△y 2 mm 3 mm 4 mm 5 mm 6 mm y+Rf Rr Error 110.960 13.472 24.955 2.65%166.440 13.688 23.519 0.61%221.920 13.731 22.315 3.70%277.400 13.762 21.484 5.84%332.881 13.802 21.240 6.38%

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

Fn Rt△y=4 mm y+ Rf Rr Error y+ Rf Rr Error△y=3 mm 0.237 0.274 0.283 0.292 0.311 0.329 27.283 37.434 40.163 43.533 50.314 58.579 146.327 166.440 171.425 176.394 186.286 196.119 10.471 13.688 14.494 15.335 17.156 19.028 16.512 23.519 25.148 27.381 31.091 37.707 0.79%0.61%1.29%1.87%2.14%3.14%195.102 221.920 228.566 235.192 248.381 261.493 10.523 13.731 14.590 15.407 17.218 19.104 15.774 22.315 23.909 26.232 30.331 36.091 3.61%3.70%4.14%4.35%5.49%5.77%

图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 波形及自由面分析

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

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

5 结 论

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

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

该矿体隐伏于上棋盘村西北约810 m，为(构造)蚀变石英脉型辉钼矿，矿体呈灰白色—褐黄色—褐红色，矿体外推长度200 m，厚度1.07 m。矿石钼品位为0.030%。

（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.