吕全亮，刘雯玉，徐 犇，刘鹏鑫

(1.江苏省镇江船厂(集团)有限公司，江苏 镇江 212000；2.江苏科技大学 船舶与海洋工程学院，江苏 镇江 212003；3.江苏省船舶设计研究所有限公司，江苏 镇江 212003)

0 引言

随着计算机和信息技术的发展、湍流模型的突破、非线性波浪数字模拟技术与动网格技术的运用，基于计算流体动力学(CFD)的船舶阻力性能研究技术愈加成熟。

许炀垲等[1]利用CFD软件，针对3种不同湍流模型，对KCS船的阻力性能进行了预报，通过与船模试验数据进行比较，得到了最适合速度范围内的湍流模型。倪崇本[2]基于粘流理论，将三因次船体阻力预报方法与CFD数值模拟手段进行了有效的结合，提出了基于CFD数值模拟计算的一种快速预报实船阻力的方法。王金宝等[3]对低速肥大型船舶进行了阻力性能研究，并对尾流场进行了数值模拟，将其方法应用于其他实船，结果较为精确。苏绍娟[4]等采用数值模拟进行船舶阻力性能预报，发现在恰当的模型与高质量的网格情况下，CFD数值模拟能很好且快速地模拟船舶航行状态。王健等[5]利用Fluent软件分别对高速单体船和高速三体船进行了数值模拟，预报了各阻力并分析了各流场的变化，为多体船的研究做出了贡献。魏可可等[6]以5415船模为研究对象，用STAR-CCM+对其船模阻力、升沉及横倾角进行了计算，并与试验值进行了对比，得到了总阻力系数、剩余阻力系数随弗劳德数Fr增大而增大的结论。查若思等[7]基于RANS方程的黏性数值模拟方法，利用OpenFOAM对多种船型在不同航速下静水阻力和兴波情况进行了模拟，分析了各船型快速性。JOSIP[8]利用流体体积表面捕获技术的RANS方程，预测了受损船舶在静水中前行时周围的流场，并将结果与Brodarski研究所的模型实验进行对比，预测了受损船舶的总阻力性能。

以往对于船舶阻力性能多限于实船，且对于船模数据与实船数据的转换多依靠于二因次法。本文首先基于粘性流体理论，采用随体动网格技术模拟船舶航行状态；其次通过三因次法进行船模数据换算，与试验数据对比验证了数值算法及换算方法的可靠性，完成了不同载况下客渡船静水中阻力性能的数值预报研究，采用理论计算给出了不同载况下车客渡船的有效功率，对车客渡船的推进器选型具有指导意义。

1 数值计算方法

1.1 基本控制方程

流体流动要受物理守恒定律的支配，最基本的守恒定律包括：质量守恒定律、动量守恒定律及能量守恒定律，在流体力学中具体体现为连续性方程、动量方程和能量方程。水通常被当做不可压缩流体，且本文也不考虑能量的交换，因此，数值计算仅需要满足连续方程和动量方程。

(1)连续性方程

(1)

(2)

(2)动量守恒方程

(3)

式中：p为流体微元体上的压力；τxx、τxy、τxz、τyx、τyy、τyz、τzx、τzy、τzz为因分子粘性作用在微元体表面上的粘性应力τ的分量；Fx、Fy、Fz为微元体上的体力，若只受重力，且z轴竖直向上，则Fx=0,Fy=0,Fz=-ρg。

对于不可压缩流体，该方程可进一步简化为：

(4)

式中：μ为k的粘性系数。

1.2 湍流模型

本文采用混合了计算近壁区域粘性流动的k-ω模型和用来计算远场自由流动的k-ε模型的SSTk-ω模型。在这个模型中，湍流动能k和特定耗散率ω的输运方程为：

(5)

Yω+Dω+Sω

(6)

与标准k-ω模型相同，其中Γk、Γω分别为k和ω的扩散系数，且表达式相同。其中，μt的定义为：

这里，普朗特数σk与σω不再为常数，其定义为：

上列两式中F1与F2的定义为：

(7)

Gk和Gω为k和ω方程的产生项，其定义为：

其中，

α∞=F1α∞,1+(1-F1)α∞,2；

Yk、Yω为k和ω方程的耗散项，其定义为：

Yk=ρβ*kω；Yω=ρβω2；β1=F1βi,1+(1-F1)βi,2

Dω为横向扩散项，是将k-ε模型转换成了k-ω模型形式产生的，定义为：

1.3 离散方法

计算流体动力学的核心内容就是离散方法，常用的离散方法为有限差分法、有限元法和有限体积法。本文采用有限体积法进行计算，其控制方程使得任何一组控制体积都能满足特征变量的积分守恒，即使在网格较粗的条件下，也能准确表现出积分守恒。

1.4 数值方法验证

在进行正式计算之前，需要对以上数值方法进行验证。选取KCS标准船型进行缩尺计算，缩尺比为31.6。其几何模型见图1，主要参数见表1。

将缩尺后的KCS船体模型导入HEXPRESS软件进行网格划分。建立计算域长宽高分别为4Lpp、3Lpp、2.5Lpp，船模距前方入口1Lpp，距出口3Lpp，距上下边界分别0.5Lpp和2Lpp，左右边界各1.5Lpp。计算域内船体采用标准壁面函数，上下边界设置为指定压强(液体静压)，入口边界、出口边界及两侧边界设定为远场边界。为更精确地得到流场数据，需要对船体球鼻艏、水线面、尾流场进行局部加密。计算域划分与船体网格分别见图2、图3。

图1 KCS船体模型

表1 KCS船主要参数

图2 计算域网格划分

图3 船体表面网格

得到缩尺船模数据之后需要将船模阻力换算成实船阻力，计算方法分别为二因次法与三因次法。

二因次法具体步骤如下：

(1)由CFD模拟求得船模的总阻力系数Ctm

(8)

式中：Rtm为船模阻力；ρ为淡水密度；Sm为船模湿表面积；Vm为船模速度。

(2)利用桑海公式计算船模和实船摩擦阻力系数Cf

(9)

式中：Re为雷诺数。

(3)求出船模的剩余阻力系数Crm

Crm=Ctm-Cfm

(10)

式中：Ctm、Cfm分别为船模阻力系数、船模摩擦阻力系数。

(4)计算实船总阻力系数Cts

Cts=Crs+Cfs+Ca+ΔCf

(11)

式中：Crs为根据弗劳德数假设Crs=Crm求得的实船剩余阻力系数；Cfs为实船摩擦阻力系数；Ca为空气阻力系数，压载状态下取0.142×10-3；ΔCf为粗糙度补贴系数，压载状态下取值为0.046×10-3。

(5)计算实船阻力公式

(12)

式中：Rts为实船阻力；Ss为实船湿表面面积；Vs为实船速度。

由于三因次法认为粘压阻力系数Cpv与摩擦阻力系数Cf之比为常数ξ，即

(13)

实船总阻力系数为：

Cts=(1+k)Cfs+Cws+ΔCf

(14)

式中：Cws为实船兴波阻力系数。

船模总阻力系数为：

Ctm=(1+k)Cfm+Cwm

(15)

式中：Cwm为船模兴波阻力系数；k为形状系数，(1+k)为形状因子，仅与船体形状有关。

又因：

Cwm=Cws

故实船总阻力系数为：

Cts=(1+k)(Cfs-Cfm)+Ctm+ΔCf

(16)

式中(1+k)形状因子可通过普鲁哈斯卡的假设求得。普鲁哈斯卡认为弗劳德数在0.1～0.2之间时，可近似认为兴波阻力Cwm与弗劳德数Fr的4次方成正比关系。故有：

(17)

通过线性拟合可求得(1+k)值。

摩擦阻力系数可用1957ITTC公式计算：

(18)

(19)

式中：ks为船体表面粗糙度，通常取1.50×10-4m。

因此，实船的总阻力为：

(20)

由以上公式将船模换算成实船数据后再与试验值进行对比，结果见表2。表中结果证明，该数值方法与船模数据换算方法能较为准确地预报船舶阻力性能。

表2 KCS船模阻力计算结果

2 计算模型

本文目标船型为60 m镇扬车客渡船模型，选用Solidworks软件完成三维建模，缩尺后导入Hexpress软件里进行网格划分，缩尺比为1∶16。在船体附近及自由液面进行特殊加密，其中：上部分为空气域，下部分为水域。计算船舶阻力时，建立计算域的长、宽、高分别为15、11.25、9.375 m。船模距前方入口1LOA(LOA为总长)，距出口3LOA，距上下边界分别0.5LOA和2LOA，左右边界各1.5LOA。计算域内船体采用标准壁面函数，上下边界设置为指定压强(液体静压)，入口边界、出口边界及两侧边界设定为远场边界。图4为船体模型三视图，图5为表面网格划分示意图，实船与船模主要参数见表3。

图4 船体三视图

图5 计算域及网格剖面

表3 车客渡船主要参数

本次计算5种吃水状态、实船航速Vs=2～12 kn共30种工况，见表4。湍流模型选用k-omega(SST-Menter)模型，k与ω的值与雷诺数有关，采用中心差分(AVLSMART)格式离散；自由面采用混合自由面捕捉与重构相结合(BRICS)格式离散；采用压力速度耦合算法求解。

3 结果及分析

3.1 车客渡船阻力及有效功率的计算与分析

本文在预报了船舶分别在A航区满载到港、满载出港、空载到港，以及B航区满载出港、满载到港的情况下，以实船航速2、4、6、8、12 kn和设计航速10 kn共30种工况下的运动状态，船模模拟速度Vm分别为0.257、0.514、0.772、1.029、1.286、1.543 m/s。当流场稳定，阻力变化稳定后，取稳定时间内数据的平均值。A航区、B航区所得船模阻力数值计算结果分别见表5、表6。

表4 计算工况

表5 A航区不同状态船模阻力系数表

本文所研究的车客渡船设计航速为10 kn，弗劳德数为0.210，属于低速船。而对于低速船来说，粘性阻力所占的比例往往比较大(可达70%以上)；兴波阻力所占的比例则比较小。此船的数据也正符合这一规律，并随着速度的增加，兴波阻力所占比例逐渐变大，总阻力值增长速度也增大。

表6 B航区不同状态船模阻力系数表

对该船采用二因次法换算其实船数据，并与三因次法换算所得数据进行比较，绘出阻力随航速变化关系曲线，见图6。

图6 三因次法和二因次法的计算结果对比

从图6中可以看出，随着速度的增加，阻力在逐渐增大。实船阻力换算方法中二因次法小于三因次法。首先，二因次法把船体阻力分成了互不相关的两个独立部分：一部分仅与重力或弗劳德数有关；另一部分仅与粘性或雷诺数有关，而忽略了两者的相互影响，事实上这两种影响是存在的。其次，弗劳德数将兴波阻力和粘压阻力这两种不同性质的阻力合并为剩余阻力，在理论上是不恰当的。最后，船体形状是相对复杂的三因次物体，其周围流动情况与平板相比显然有一定差别，因而用相当平板的摩擦阻力来代替船体摩擦阻力，必然是有误差的。因此，三因次法计算阻力更加符合实际情况。本文采用三因次法给出船舶不同工况下的阻力，见图7。船舶附体阻力取船舶阻力的10%，得到总阻力后乘以对应航速得到船舶的有效功率，见图8。

图7 车客渡船阻力特性曲线

从图7、图8中可以看出：

(1)船舶阻力随着航速和排水量的增加而增加，呈二次型增加趋势。

(2)空载到港阻力最小，设计航速下的阻力约为67.6 kN，有效功率约为347 kW。

(3)A航区航行中满载出港阻力最大，设计航速下的阻力约为85 kN，有效功率约为450 kW。

图8 车客渡船有效功率特性曲线

3.2 车客渡船流场分析

图9给出了不同航速下车客渡船自由面波形图。随着航速增加船体兴波愈加显著，与常规单体船开尔文波形不同。由于车客渡船下潜体呈梭型，具有分流作用，在船体尾流场重新汇聚，尾波形较为分散且不断向外扩散。

车客渡船船体表面压力分布云图见图10。从图中看出，船底梭型潜体艏部为压力驻点形成高压区，对流体形成分流作用，降低了高压区的面积，流体绕过梭型最宽区域后，梭型体后半部分由于流体分离形成涡流区域，出现一个低压区，所以整体梭型结构相对于方形结构具有明显的导流减阻的效果。

图9 车客渡船自由面波形图

图10 车客渡船船体表面压力分布云图

4 结论

本文基于粘性流体理论，完成了梭型车客渡船不同工况下的阻力数值预报，采用三因次法将船模阻力换算至实船阻力，并对车客渡船尾流场进行分析，得到以下结论：

(1)本文所建立数值预报方法可较为准确地预报船舶阻力及流场分布规律。

(2)三因次法船模阻力换算至实船阻力考虑因素更加贴合实际情况，准确性较高。

(3)车客渡下潜体采用梭型结构相对于常规方形底，具有分流减阻效果。