基于重叠网格M型船水气二相流数值模拟研究

余泽爽 毛筱菲

(武汉理工大学交通学院1) 武汉 430063) (武汉理工大学教育部高性能船重点实验室2) 武汉 430063)

基于重叠网格M型船水气二相流数值模拟研究

余泽爽1,2)毛筱菲1,2)

(武汉理工大学交通学院1)武汉 430063) (武汉理工大学教育部高性能船重点实验室2)武汉 430063)

M型船是一种新型滑行艇，同时具备普通滑行艇的阻力特性、高速多体船的耐波性能，以及微气泡减阻船优异的阻力性能.为了实现对M型船不同航行阶段中复杂变航态过程的动态捕捉、对阻力的精确预报，以及对非排水航行阶段复杂的水气二相流场的模拟，基于粘性CFD软件STAR-CCM+的重叠网格技术，对M型船在静水中从排水、过渡到滑行航行阶段进行数值模拟，从船体表面网格大小讨论了网格独立性，结合计算资源以确定最佳的网格方案；将阻力和纵倾角的计算值与试验值进行对比，从湍流模型与y+值以及时间步长独立性讨论了不同航行阶段的模拟需求；从水气二相流的模拟讨论了M型船在不同航行阶段的减阻原理.

M型船；重叠网格；水气二相流；湍流模型；时间步长适应性

0 引 言

M型船是为避免游船在高速航行时产生的兴波对河岸的冲刷破坏而发展的新型消波船型，因其横剖面形状为M形而得名，由此还衍生出W型船、双M型船等船型.M型船因其独特的几何形状和优秀的流体动力性能，综合了普通滑行艇、高速多体船和微气泡减阻船等高性能船的优点，能够充分利用滑行效应、浮力、表面效应，以及首波与空气捕捉等原理，产生高效的浮力、水气动升力来支持艇重，从而降低船体的总阻力、提高航速.M型船由中心排水体(CDS)、滑行涵道(PT)和刚性侧裙(RS)组成.

双M槽道滑行艇(短剑Stiletto)在美国进行了试航和验收，并交付近海特殊作战部队，此外还是首次采用了碳纤维湿法缠绕环氧复合材料的整体成型技术的新型高速艇[1].迟玉鹏等[2]通过系列船模试验，对槽道艇和单体V型艇阻力和耐波性能进行了对比分析，发现槽道艇具有横摇小、喷溅少和垂向加速度小的优点，并提出了阻力最优的重心纵向位置.苏永昌等[3]通过系列试验，对槽道滑行艇的基本原理、运动特性及其变化规律进行了研究，并重点讨论了船体型线、槽道尺寸和形状，以及重心纵向位置等影响因素.陈辉[4]基于Fluent对M型船静水阻力进行数值模拟，并与试验结果进行对比，在非滑行阶段(Fr▽<3)较为精确，但当Fr▽>3时数值计算误差较大.唐建飞等[5-6]通过系列船模试验，讨论了排水量和重心纵向位置对M船型阻力和航态的影响，提出了该船型存在两个阻力峰值的特性，同时还比较了M船型与普通槽道滑行艇的阻力与耐波性能，提出M船型优良的水动力性能.Ghassabzadeh等[7]采用动网格技术(dynamic mesh)，考虑了自由面水气二相流场，对槽道多体船不同航行阶段进行了模拟，实现了对高速运动的模拟和流场信息的捕捉.邹劲等[8]基于STAR-CCM+，讨论了船表网格尺寸、近船面网格节点分布、船体周围加密区网格尺寸以及变形网格技术等对三体滑行艇的阻力预报、收敛速度以及运动稳定性的影响.

综上，M型船数值模拟研究的难点在于：①航行过程中大幅度的航态变化，对网格要求极高，普通网格已无法满足模拟需求；②M型船航速变化范围很大，对不同航行阶段，有不同的模拟需求；③滑行涵道中水气两相流场的模拟，对网格、壁面处理和湍流模型提出更严格的需求.

1 理论背景

1.1 控制方程

流体的流动需满足质量守恒定律，对于不可压缩流体，即有连续性方程：

(1)

式中：u，v，w分别为x，y，z方向上的分量.

流体流动还应满足动量守恒定律，即由Naviver-Stokes提出的N-S方程：

(2)

Su=Fx+sx=0

Sv=Fy+sy=0

Sw=Fz+sz=ρg

(3)

直接求解N-S方程过于困难，而雷诺平均的RANS方程则提供了便利.RANS方程是平均变量的控制方程，基于湍流模式假定，将高频脉动的流场变量分解为时均量与脉动量之和.

1.2 自由面捕捉

自由面捕捉采用的VOF方法，根据网格单元中不同流体的体积分数的变化来模拟自由面的位置和形状，流体的体积分数即为流体体积与网格单元体积之比：

(4)

式中：f1，f2分别为空气相和水相的体积分数，其控制方程为

(5)

以体积分数f表示时，f=1即为水相，f=0即为空气相，f=0.5即为自由面.

1.3 湍流模型

由于湍流现象的高度复杂性，所以至今未有一种通用的方法能够精确模拟所有流动问题中的湍流现象.雷诺平均的N-S方程虽然解决了方程的变量问题，但是计算中涉及的湍流问题，都是以湍流模型的形式添加到方程中的[9].

Star-CCM+中提供了k-ε模型、SSTk-ε模型、Spalart-Allmaras模型和雷诺应力湍流模型，本文使用和讨论了前两个模型.

1.4 边界层的处理

粘性流体绕流静止的壁面时，由于物面为不可滑移面，故在紧靠物面的薄层内，流速由物面上的零增加到与来流类同，这种在高雷诺数下，在紧靠物面处形成的薄层即为边界层.在边界层中，速度梯度很大，湍流增强，因此壁面处网格的处理对捕捉该处流场的物理特性格外重要.

在使用湍流模型时，一般采用壁面函数法对壁面处网格做特殊处理，.通常采用无因次距离y+来定义船体表面第一层网格厚度，可采用经验公式进行估算

(6)

式中：Δy为船体表面第一层网格厚度；L为湿体长度.

2 数值模拟与试验验证

2.1 模型试验

本文采用的模型和实验数据均来自文献[5-6]，图1为该M型船横剖面图、使用CATIA三维建模软件建立的M型船三维模型、实验用船模图.M型船主要参数及实验数据见表1～2.

图1 M型船

表1 M型船主要参数

表2 M型船模型阻力试验数据

2.2 计算设置

图2为计算域和网格区域的划分情况，分别表示了重叠网格区域、船表周围加密区、自由面加密区以及为了捕捉船后兴波的兴波加密区.

图2 区域划分情况

2.3 网格独立性研究

由重叠网格特点，重叠网格区域网格尺寸与船体表面网格一致，背景网格为重叠网格的四倍，由于重叠网格与背景网格之间采用线性插值，因此船表周围加密区网格为重叠网格两倍，自由面加密区网格(Z向)为重叠网格的80%，船后兴波加密区(XY向)为重叠网格的2倍.表3为网格数量随船表网格尺寸变化得到的不同网格方案.

表3 网格数量随船表网格尺寸变化

选取k-ε湍流模型，边界层以固定壁厚给定，计算时，保持y+值为30～200，时间步长取为0.001 s，时间离散为二阶形式.图3为不同网格方案计算的阻力与纵倾角随网格数量变化的趋势与试验值的对比.

图3 不同网格方案计算值随网格因子的变化趋势

网格因子RGi为单元网格体积之比，即网格数量反比的三次方根，其表达式为

RGi=(Ni/N1)1/3

(7)

式中：Ni为方案i的网格数量.

由图3可知，计算结果随网格数量增加呈现单调趋势，说明网格收敛性得到验证.通过比较数值变化趋势，综合考虑计算资源，认定方案3作为后续研究的网格方案，因为在此方案基础上，增加网格数量对计算结果的影响并不显著，即网格独立性得到满足.

2.4 湍流模型与y+值适应性研究

2.4.1k-ε湍流模型

选用k-ε湍流模型时，分别取y+为50，100和200，根据式(6)估算得到船表第一层网格高度Δy，取边界层数为5，厚度增长率为1.3，则可得到边界层厚度δ.表4为采用k-ε湍流模型时，不同航速下阻力和纵倾角随y+值变化的计算值与试验值的误差.

表4 y+值对阻力和纵倾角的影响

由表4可知，在排水航行阶段(Fr▽<1)，y+值对计算结果的影响并不大，与实验结果的误差均保持在5%以内，横向比较可知选用y+=200为最佳；在过渡航行阶段(13)，误差则存在极大不确定性，可见在该航行阶段，k-ε湍流模型已无法满足模拟需求，因此需要考虑SSTk-ε湍流模型.

2.4.2SSTk-ε湍流模型

船舶在流体中航行时，随着航速的提高，雷诺数Re也随之提高，边界层则随之变薄，由壁面函数公式(6)可知，此时的y+值必然降低，而由文献[9]的相关研究，在选用SSTk-ε湍流模型时，y+应不大于1，故这里选取y+为0.5和1.表5为选用SSTk-ε湍流模型时，不同航速下阻力与纵倾角随y+值变化的计算值与试验值的误差.

表5 y+值对阻力和纵倾角的影响

由表5可知，在非滑行航行阶段(Fr▽<3)，计算结果整体的误差均较大，不如采用k-ε湍流模型时的结果；而在滑行航行阶段(Fr▽>3)，当y+=1时的结果误差较小低于8%，故在采用k-ε湍流模型模拟滑行航行阶段时，y+应接近1为最佳.

2.5 时间步长适应性研究

在数值模拟中，时间步长的选取对计算结果会造成一定影响，一个合适的时间步长能使非定常的数值模拟得到同时满足准确性、稳定性和计算成本的需求[10-11].

这里主要讨论滑行航行阶段采用k-ω湍流模型的时间步长适应性，这里选取时间步长为0.01,0.005,0.0025,0.001和0.000 5 s，时间离散采用二阶形式.表6为采用k-ε湍流模型时，滑行航行阶段阻力和纵倾角随时间步长变化的计算值与试验值的比较.

表6 时间步长对阻力和纵倾角的影响

由表6可知，随着时间步长的减小，计算结果误差并不一定越小，但是结果基本上呈现单调的收敛关系.同时，对于高速情况，时间步长不宜过大，否则会出现发散的情况，这也一定程度说明k-ε湍流模型稳定性较差的特点，此外，综合分析阻力与纵倾角角的计算值与试验值的误差为标准，再考虑到时间步长对计算收敛过程的影响，兼顾计算资源的经济性，选择最佳时间步长为0.005 s，时间离散为二阶形式.

2.6 试验验证

由上述网格独立性、湍流模型、y+值以及时间步长适应性的研究，对M型船从排水、过渡到滑行航行阶段的静水直航运动进行数值模拟，并将计算的阻力与纵倾角与试验结果对比见图4.在全航行阶段，采用上述数值模拟方法计算的阻力与纵倾角，与试验结果基本吻合，误差均在12%以内.

图4 数值模拟与试验验证

3 M型船特性研究

3.1 水气两相流场特性分析

由于滑行涵道特殊的喇叭形结构，M型船在航行过程中，中心排水体艏部产生的兴波，随空气一同进入涵道，而涵道出口一般仅仅略高于静水面，此外高速航行时，滑行面与液面的剧烈摩擦，产生大量水气，此时，涵道中形成复杂的水气混合物，并且压力逐渐上升，最终在艉部以喷溅的形式从出口溢出，形成集中的尾兴波.

图5为不同航速时滑行涵道中水气分布图.由图5可知，随着航速的提高，滑行涵道中水气混合物逐渐减少，同时对比不同位置处的水气分布情况，可见接近船尾的水气分布更密集更复杂，这是因为高速情况下，由于涵道中急剧升高的压强，水气混合物从出口压出.

图5 不同航速时滑行涵道中水气分布

图6为不同航速时船体表面水气分布图.由图6可知，在排水和过渡航行阶段，随着航速的提高，船体表面水气分布呈现有规则的水气分层现象，一方面是由于航态纵倾，另一方面是由于水气混合物此时主要集中于涵道中，并未出现剧烈的溢出现象；而在滑行航行阶段，船表水气分布显得更为紊乱，这是由于涵道中的水气混合物以喷溅的形式从出口溢出，同时滑行面与液面之间高速的相对摩擦运动，引起近壁面处出现局部的气化现象.

图6 不同航速时船体表面水气分布

图7为不同航速时的船后兴波图.由图7可知，在排水航行阶段，船后兴波的横波系和散波系与普通滑行艇类似；在过渡航行阶段，滑行涵道中的水气混合物开始少量溢出，由于涵道收缩的结构，使液面被压缩，流体以集束的形式从出口流出，由于两股束流之间的相互干扰，同时侧壁刚裙对中心排水体兴波的阻隔作用，船后兴波相比排水航行阶段要更为收缩；在滑行阶段，涵道中的水气混合物因高压从出口急剧溢出，使船后兴波进一步收缩，而不会出现明显的“空穴”和“鸡尾流”现象.

图7 不同航速时的船后兴波

3.2 减阻原理分析

3.2.1滑行升力

与普通滑行艇类似，特殊的滑行面形状在航行中提供一定的水动升力，使船体抬离水面，达到减阻的目的.图8为不同航行阶段船表压力.由图8可知，在不同航行阶段，随着航速的提高，船表压力逐渐提升.

图8 不同航行阶段船表压力

3.2.2表面效应升力

由于涵道的特殊形状，中心排水体产生的艏兴波随同空气，以螺旋状的运动方式进入涵道，并在涵道内被加压、分离，形成外旋的水气涡流，可视为纵向的二元涡，带自由面的二元涡产生表面效应升力.表面效应升力L表达式为

L∝kρω2r3l

(8)

式中：k为系数；ω为二元涡角速度；r为二元涡半径；l为二元涡的长度.

图9为不同航行阶段涵道中的涡流，由图9可知，随航速增加，涡流由船的艏部逐渐向后延伸，由式(8)可知，该升力也随之增加，说明航速越高，槽道涡表面效应升力越大.

图9 不同航行阶段涵道中的涡流

3.2.3水气的润湿作用

由图5～6涵道和船表水气分布情况可知，尤其在滑行阶段滑行面上产生的气化现象，能显著降低摩擦阻力系数，从而达到减阻效果.

3.2.4抑制船行波系

由于侧壁刚裙平直的外壁能降低船体压力梯度，减小形状阻力，而内壁可以屏蔽中心排水体产生的兴波.由图7可知，随着航速增加，船后兴波呈现收缩的趋势.图10为船后不同距离处，船行波波面高度分布图，由图10可知，随航速增加，船行波波面高度分布更集中.

图10 不同航速的船行波

4 结 论

1) 通过网格收敛性和独立性研究，建议M型船数值模拟时，船表网格尺寸宜取为0.012 5m，即0.625%L.

2) 在排水航行阶段(Fr▽<1)，宜采用k-ε湍流模型，且y+值应取200为宜；在过渡航行阶段,宜采用k-ε湍流模型,并且其前半段(13)，宜采用SSTk-ε湍流模型，且y+值应接近1为佳.

3) 选用SSTk-ε湍流模型时，通过时间步长适应性研究，推荐选用时间步长为0.005 s，时间离散为二阶.

4) 通过对M型船水气二相流的数值模拟研究，发现滑行涵道中的水气混合的现象，主要发生在排水和过渡航行阶段，而滑行阶段主要是在滑行面上发生了气化现象.

[1] 赵连恩,王庆.高性能船(舰)原理与技术[M].哈尔滨：哈尔滨工程大学出版社,2013.

[2] 迟云鹏,孟宪钦.高速槽道艇阻力及耐波性能试验研究[J].船舶工程,1995(3):27-31.

[3] 苏永昌,赵连恩.高性能槽道滑行艇的运动特性[J].中国造船,1996(1):11-16.

[4] 陈辉.M船型水气两相流场特性研究[D].武汉：武汉理工大学,2011.

[5] 唐建飞,黄武刚.M船型阻力模型试验研究[J].中国舰船研究,2014(5):49-52.

[6] 黄武刚.M型艇与槽道型艇的阻力和耐波性比较[J].船海工程,2015,44(1):56-59.

[7] GHASSABZADEH M, GHASSEMI H. Determining of the hydrodynamic forces on the multi-hull tunnel vessel in steady motion[J]. Journal of the Brazilian Society of Mechanical Sciences and Engineering,2014,36(4):697-708.

[8] 邹劲,姬朋辉,孙寒冰,等.网格因素对三体滑行艇阻力计算影响探究[J].船舶,2016,27(3):8-14.

[9] 章丽丽,孙寒冰,蒋一,等.三体滑行艇槽道的水气动力特性研究[J].哈尔滨工程大学学报,2017,38(1):31-36.

[10] 王福军.计算流体动力学分析：CFD软件原理与应用[J].热能动力工程,2005(5):468-468.

[11] 孙华伟.滑行面形状对滑行艇阻力与航态影响数值分析[D].哈尔滨：哈尔滨工程大学,2012.

Research of the RANS Simulation of Two-phase Flow Around M-hull Ship Based on Overset Grid

YUZeshuang1,2)MAOXiaofei1,2)

(SchoolofTransportation,WuhanUniversityofTechnology,Wuhan430063,China)1)(KeyLaboratoryofHighPerformanceShipTechnologyofMinistryofEducation,WuhanUniversityofTechnology,Wuhan430063,China)2)

M-hull ship is a new kind of planning craft with the resistance characteristics of the normal planning craft, the sea-keeping performance of the high-speed multi-hull ship and the excellent resistance performance of the micro-bubble drag reduction ship. In order to achieve the goal to dynamically capture the complex changing navigation state during different navigation phases, accurately predict the resistance performance and simulate the two-phase flow during non-navigation phase of displacement, the RANS simulation from the navigation phase of displacement to transition and planning in the static water is carried out based on the overset gird technology of Star-CCM+. The mesh independence study is performed according to the grid size closed to the ship surface and an optimal mesh scheme is determined combined with the computing resource. Comparing the compute values of the resistance and trim angle and the experiment values, the influence between turbulence model and they+value, and the time step adaptive is discussed during different navigation phases. The drag reduction principle during different navigation phase is also discussed according to the simulation of the two-phase flow around the M-hull ship.

M-hull ship; overset grid; two-phase flow; turbulence model; time step adaptive

U661.32

10.3963/j.issn.2095-3844.2017.05.030

2017-08-31

余泽爽(1992—)：男，硕士生，主要研究领域为船舶性能