0.   引言

随着绿色船舶规范的实施以及内河标准船型建设的推进，内河船舶的能效指标越来越受到业界的关注。交通运输部更是要求^[1]：对于新建、改建的内河运输船舶，应当满足内河运输船舶标准船型指标体系中的强制性指标要求。强制性指标^[2]包括船舶主尺度系列标准、CO₂排放指标和燃料消耗指标，其中CO₂排放指标和燃料消耗指标均与船舶能效指数相关联。船舶能效指数评估包括能效指数前期评估和建造后评估2类^[3-4]。用于能效设计指数前期评估的功率曲线一般基于水池试验。在船型设计初期，可以通过水池试验评估船舶快速性能，以获得优良的船型。但进行船模试验需要制作船体模型，费用高、周期长，且不易对船型进行多方案调整。因此，在船舶设计初期，如何精确预报船舶阻力，获得流场信息，以掌握船舶快速性性能指标并方便后续阶段的船舶优化改进是船舶设计和能效指标验证的关键。

随着计算机技术的飞速发展和计算数学理论的不断完善，计算流体动力学(Computational Fluid Dynamics，CFD)得到了蓬勃发展，成为研究船舶航行性能的重要手段之一。相比于模型试验方法，采用CFD方法不仅费用低、周期短，而且还可以准确获得船舶周围流场信息。目前，国内一些单位采用基于RANS方程加自由面的处理方法对船舶快速性能进行的预报在工程上基本能满足其精度要求^[5-8]，但以上这些研究大多是基于海船船型。受吃水及航道的限制，内河船舶多采用双艉或者三艉设计，船型宽大而型深较小，主要船型系数与海船相比差异较大，而目前有关双艉船型阻力预报及流场捕捉方面的研究还较少。

FINE/MARINE是船舶与海洋工程水动力专业使用的粘性流计算软件^[9]，主要包括网格生成器HEXPRESS、后处理工具CFView和不可压粘性流场求解器ISIS-CFD。其特色在于，HEXPRESS为全六面体非结构网格生成工具，采用体到面的网格生成方式可以快速生成高质量的网格，并且其对船舶海洋工程问题的模拟可以通过界面设置来实现，无需进行二次开发，方便快捷。

本文将基于计算流体动力学理论，提出一种双艉船型航行阻力预报及流场捕捉的工程方法。采用FINE/MARINE软件，对某双艉客船的航行阻力进行计算，探讨船舶第1层网格节点高度、船舶航速以及附体对计算结果的影响。通过不同航速下的阻力预报，并将数值结果与试验数据进行比较，以验证所提方法的正确性。

1.   数学模型

本文对水面船舶自由面扰流问题的数值模拟是将自由面流动作为两相流来处理，自由面为水和空气的交界面，采用流体体积法(Volume of Fluid，VOF)求解。水面船舶自由面扰流问题数学模型的控制方程包括：连续性方程、体积分数方程、动量方程和湍流模型方程。

1.1   控制方程

不可压缩粘性流体的连续性方程和动量方程分别为：

式中：t为时间；ρ为密度；v为控制体；A为控制体面积；n为法矢量；U_d为控制体表面法矢量方向上的速度；U为速度矢量；p为压力；U_i为x_i坐标轴方向上的平均速度分量；τ_ij为粘性应力张量；g_i为重力矢量；I_i和I_j为方向向量。

1.2   湍流模型

本文采用SST k-ω湍流模型，其湍动能k和湍流耗散率ω的输运方程如下：

式中：μ为分子粘度；U_j为x_j坐标轴方向上的平均速度分量；μ_t为湍流涡粘度；t_ij为湍流雷诺应力张量；S_ij为平均应变率张量；F_I为辅助混合函数；P_ω为ω的导出项；β为湍流应力常数；β^*，σ_k和σ_ω²均为为湍流模型常数。

1.3   自由面捕捉算法

自由液面捕捉法是将空气和水作为单一流体同时计算。该单一流体的性能(质量体积和粘性因数)在空间的变化取决于构成函数c_i。在自由液面计算中，c_i在空气中的取值为0，在水中的取值为1，通过求解下式确定：

自由液面捕捉算法具有更好的灵活性和适应性，可以较好地模拟破碎波等复杂的自由液面。

1.4   离散格式

控制方程的离散项采用隐式有限体积法，直接求解三维粘性不可压缩多相流的雷诺平均方程，其具有二阶空间和时间精度。动量方程离散采用GDS格式，时间离散采用时间步进算法。自由液面捕捉采用BRICS可压缩型离散格式，能减小自由液面附近构成函数的数值扩散。

2.   研究对象

研究对象为航行于三峡库区的某双艉客船，该客船为双机、双桨、双舵，船舶线型为双艉线型、V型船艏。该船由于船型特殊，缺乏相关的资料，为准确预报其船舶阻力，在武汉理工大学的拖曳水池进行了阻力测试。本文对船模进行了CFD模拟计算，主要船型参数如表 1所示，试验船模及布置图如图 1所示。

表  1  双艉客船船型主要参数

Table  1.  Main parameters of the twin-skeg passenger ship

	实船	船模
垂线间长Lpp/m	80.800	5.387
型宽B/m	16.000	1.067
吃水T/m	2.500	0.167
湿表面积Sw/m2	2 339.793	5.466
方形系数C b	0.593	0.593
缩尺比λ	15.0	15.0

下载: 导出CSV 

| 显示表格

图  1  船模及试验布置图

Figure  1.  Model ship and test arrangement

下载: 全尺寸图片 幻灯片

3.   数值计算

3.1   计算区域及边界条件

考虑到船型及流动的对称性，进行数值模拟时只计算一半区域即可。计算区域如图 2所示，其中，计算边界定义如下：

图  2  数值计算区域及网格划分

Figure  2.  Numerical calculation domain and mesh scheme

下载: 全尺寸图片 幻灯片

1) 前端——模型首部前约1.0L_PP处；

2) 后端——模型尾部后约3.0L_PP处；

3) 侧边界——模型侧方约1.5L_PP处；

4) 上边界——水线以上约0.5L_PP处；

5) 下边界——水线以下约1.5L_PP处；

6) 对称面——模型中纵剖面的延展面。

计算域的上、下边界取为压力边界，前、后截面及远离船体的一侧取为远流场边界，靠近船体的一侧取为对称边界，船体甲板面为滑移物面，其他表面为不可滑移物面。因在求解过程中存在着时间偏导项，故船体从静止到某一航速，给定加速时间，然后船舶按照设计航速航行直至收敛。

3.2   网格划分

计算网格质量影响数值计算的正确性和精确性。本文采用软件自带的前处理器HEXPRESS生成全六面体非结构化网格。通过将网格粗化、细化、吸附以及在自由面附近进行细化形成贴体网格。其中，第1层边界层网格间距及边界层层数根据无量纲特征数Y⁺确定。分别对双艉客船带舵和不带舵这2种状态下的阻力进行计算，其中带舵的计算网格总数为82万，不带舵的网格总数为76万。网格质量通过正交性来保证，即绝大部分网格的正交性为90°。其中，不带舵网格的最小正交性≥22.5°，对于带舵的网格，由于舵尾缘尺寸较小，故最小正交性≥8°，网格质量满足数值计算要求。船体表面网格划分如图 3所示。

图  3  船体表面网格划分

Figure  3.  Mesh scheme of the hull

下载: 全尺寸图片 幻灯片

4.   计算结果分析

为方便对结果的处理，引入了无量纲系数：傅汝德数Fr、雷诺数Re、船舶总阻力系数C_T和摩擦阻力系数C_{F_ITTC}，分别定义如下：

式中：V_m为船模速度；g为重力加速度；ρ_w和μ_w分别为水的密度与粘度；A_w为船舶静水湿表面积；R为船模总阻力。

4.1   对带舵和不带舵的影响分析

双艉客船的水池试验考虑带舵的影响。在进行CFD数值计算时，为便于比较分析舵对船艉流动的影响，分别进行了裸船体和带舵船体2种状态下的多航速计算，其中Fr=0.099~0.217。总阻力系数的数值计算结果及其与试验数据的比较如表 2所示。表中：C_{T_EXP}为水池试验结果，考虑了带舵的情况；C_{T_CAL1}为裸船体阻力计算结果；C_{T_CAL2}为带舵船体阻力计算结果。定义误差：

表  2  双艉客船不同航速下阻力计算结果及其与试验数据的对比

Table  2.  Comparison of numerical results and experimental data for total resistance of the twin-skeg passenger ship at different speeds

Fr	Vm（/ m·s-1）	Re	CT_EXP	CT_CAL1	ERR1 /%	CT_CAL2	ERR2 /%
0.099	0.717	4.075×106	4.113×10-3	4.099×10-3	-0.334	4.112×10-3	-0.024
0.118	0.861	4.891×106	4.090×10-3	3.982×10-3	-2.641	4.000×10-3	-2.200
0.138	1.004	5.706×106	3.990×10-3	3.900×10-3	-2.256	3.958×10-3	-0.802
0.158	1.148	6.521×106	3.891×10-3	3.879×10-3	-0.308	3.949×10-3	1.491
0.178	1.291	7.336×106	3.970×10-3	3.892×10-3	-1.965	3.965×10-3	-0.126
0.217	1.578	8.966×106	4.018×10-3	3.948×10-3	-1.724	4.019×10-3	0.025

下载: 导出CSV 

| 显示表格

ERR₁=(C_{T_CAL1}-C_{T_EXP})/C_{T_EXP}

ERR₂=(C_{T_CAL2}-C_{T_EXP})/C_{T_EXP}

由表 2可以看出，不带舵和带舵2种状态下的阻力数值计算结果与试验结果均吻合良好，误差在3%以内。在不考虑舵影响的计算中，也即裸船体时，虽然网格质量和计算效率都很高，但计算结果与试验结果相比普遍偏小。相比于裸船体的计算，考虑带舵计算时网格数量会有所增加，网格划分难度也较大，但其计算结果与试验结果更为吻合。图 4和图 5所示分别为裸船体和带舵船体在Fr=0.178下航行时的船艉流场对比。从图中可以看出：一方面，舵作为附体会增加船舶阻力；另一方面，也会改变船艉的流线分布。为更真实地反映水池试验状态，在后续的算例中，将均采用带舵船体进行计算分析。

图  4  Fr=0.178时裸船体艉部流线图

Figure  4.  Flow visualization of the ship stern without rudder at Fr=0.178

下载: 全尺寸图片 幻灯片

图  5  Fr=0.178时带舵船体艉部流线图

Figure  5.  Flow visualization of the ship stern with rudder at Fr=0.178

下载: 全尺寸图片 幻灯片

4.2   对Y⁺的影响分析

在进行船舶CFD绕流问题模拟时，通常将船体表面视为无滑移壁面。根据流动特性的不同^[10]，壁面边界层可以分为3个区域：粘性底层、过渡层和惯性子层。边界层的位置通常用无量纲特征数Y⁺来表示。无量纲分析和试验工作表明，在边界层特定的Y⁺范围内，具有特定的流动特点。对大多数高雷诺数流动来说，一般在粘性底层内使用壁面函数，以节约计算资源。在惯性子层区域，由于湍流应力近乎常数，故主流方向的流速分量呈对数分布。因此，为捕捉壁面附近流动的物理特性，可在物理特征变化较大的地方分布大量网格点，尤其是第1层网格点应落在此区域，且必须布置足够的网格数以保证分辨所有的流动梯度。这就要求正确布置边界层第1层网格点的位置以及边界层层数。对于无量纲特征数Y⁺，可采用下式进行估算：

式中：y_wall为第1层网格距离；V_ref为来流的参考速度；L_ref为参考长度；υ为粘性系数。

针对本文所研究的双艉船舶，以V_m=1.721 m/s为例，对应的Fr=0.237，Re=9.781×10⁶，该航速下对应的水池试验结果为C_{T_EXP}= 4.247×10^-3。研究此航速下船舶第1层网格节点高度以及网格边界层层数对计算结果的影响，其中船舶初始网格设置及网格细化准则均与前述带舵航行时船舶阻力计算的设置相同。分别设置船舶边界层网格层数为6~10层，网格变化率为1.2，逐步改变船舶第1层网格节点高度和边界层层数，以查看边界层网格方案对计算结果的影响。计算结果如表 3所示。

表  3  不同网格方案计算结果

Table  3.  Simulation results of different mesh scheme

边界层网格数	ywall	CT_CAL	ERR/%	Y+
10	0.010 5	4.151×10-3	-2.247	48
10	0.001 65	4.152×10-3	-2.228	48
10	0.001 065	4.157×10-3	-2.110	47
6	0.001 065	4.130×10-3	-2.736	49
10	0.000 767 9	4.255×10-3	0.194	45
8	0.000 767 9	4.206×10-3	-0.943	47

下载: 导出CSV 

| 显示表格

由表 3可以看出，随着边界层第1层网格点位置的减小，阻力计算结果有所增加，而Y⁺则是逐步减小。当第1层网格节点距离相同时，不同的边界层网格数对计算结果的影响也不一样，较多的边界层网格数能够更好地捕捉边界层内速度梯度的变化，计算结果将更加精确。对于船体近壁面流动计算，推荐的Y⁺取值范围为20~60，边界层网格数约为8~10层。

4.3   不同航速下阻力的数值计算分析

双艉客船的设计航速范围较大，在进行船舶阻力模拟计算时，为进一步验证文中数值模拟方法的可靠性与可行性，对较高航速下的船舶总阻力也进行了计算，例如Fr取为0.237，0.257和0.276时。在高航速的数值计算中，考虑了最优Y⁺取值及边界层数量。表 4给出了带舵双艉客船在全航速范围内(Fr=0.099~0.276) 不同航速下的总阻力数值计算结果，并与试验值进行了比较分析。

表  4  双艉客船阻力计算结果

Table  4.  Resistance calculated results of the twin-skeg passenger ship

Fr	Vm/（m·s-1）	Re	CT_EXP	CT_CAL	ERR/%
0.099	0.717	4.075×106	4.113×10-3	4.112×10-3	-0.024
0.118	0.861	4.891×106	4.090×10-3	4.000×10-3	-2.200
0.138	1.004	5.706×106	3.990×10-3	3.958×10-3	-0.802
0.158	1.148	6.521×106	3.891×10-3	3.949×10-3	1.491
0.178	1.291	7.336×106	3.970×10-3	3.965×10-3	-0.126
0.197	1.434	8.151×106	4.027×10-3	3.975×10-3	-1.291
0.217	1.578	8.966×106	4.018×10-3	4.019×10-3	0.025
0.237	1.721	9.781×106	4.247×10-3	4.255×10-3	0.188
0.257	1.865	10.596×106	4.397×10-3	4.298×10-3	-2.252
0.276	2.008	11.411×106	4.631×10-3	4.529×10-3	-2.203

下载: 导出CSV 

| 显示表格

由表 4可以看出：在所计算的速度范围内，尽管航速变化较大，但和水池试验结果相比，数值计算结果的误差均在3%以内，说明本文中数值模拟精度较高，能满足工程需要。

不同航速下，船舶各阻力占总阻力的百分比不一样。一般对水面排水型舰船而言，低速航行时，兴波较小，摩擦和粘压阻力占主要成分。而当航速增加时，船体兴波成分会逐渐增加，船舶兴波阻力随之增大。图 6所示为不同航速下船艏兴波高度的变化。从中可以看出，随着船舶航速的增加，船艏兴波高度也相应增大，因此在双艉客船高速航行时，应考虑船艏消波的必要性。

图  6  不同航速下船艏兴波高度的变化

Figure  6.  Variation of wave height of the bow at different speeds

下载: 全尺寸图片 幻灯片

图 7所示为不同航速下的自由液面波形图。从中可以看出，不同航速下船舶自由液面的波形差别较大。这同时也表明文中所提计算方法是正确的，能准确反映物理现象，可为后续船型的改进提供参考。

图  7  不同航速下船舶自由面波形

Figure  7.  Free surface wave of the ship at different speeds

下载: 全尺寸图片 幻灯片

5.   结语

本文基于计算流体动力学理论，提出了一种双艉船型航行总阻力及流场信息捕捉的工程方法。采用粘性流计算软件FINE/MARINE，对某双艉客船的阻力特性进行预报分析，探讨了船体壁面第1层网格点的位置及边界层层数、船舶航速以及附体对总阻力计算结果的影响。最后，对某双艉客船带舵船体在全航速范围内的总阻力进行了数值计算，并与试验数据进行了比较，结果显示，误差均在3%以内，且阻力是随着航速的增加而增大。研究表明，本文所提方法计算效率较高，易于实现，预报精度较高，能满足工程需要，具有较强的工程实用性。