水面比降对船舶航行阻力影响的数值研究

冀 楠，杨 光，舒麟棹，钱志鹏，杨 春

(重庆交通大学 航运与船舶工程学院，重庆 400074)

0 引 言

船舶航行过程中经常会遇到有水面比降的情况，例如：明渠河道、受洪水涨落或潮汐影响的河道(长江南京至镇江河段[1])或者河底高程沿流向逐渐降低的河道水面均会产生比降；海洋中长波波面也可以认为是有比降水面；山区河流因河道高程和宽度多变, 水面比降现象更常见。船舶航行于带有比降的水面时会产生不同于一般水面的水动力特性，导致船舶航行姿态改变，进而影响航速甚至操纵性能的恶化，容易造成事故。因此，为了提高船舶航行在有比降水面时的安全性，有必要对水面比降对船舶的航行阻力与浮态影响进行研究。

对水面比降影响船舶航行问题的研究主要有基于经验公式的理论分析方法和模型试验法，研究范围基本集中在内河航道的上滩问题上。基于经验公式的理论分析方面，曹民雄等[2]利用急流滩模型的船模上滩试验资料，对比分析了几种内河船舶上滩阻力计算方法；周晓杰等[3]根据航道规范和理论计算，对流速-比降组合值对应的船舶上滩指标进行计算分析，得到了研究河段对应的规范设计尺度；董晓韬[4]采用船舶推力-航行阻力平衡的方法，计算多种功率、载重条件下船舶自航上滩的临界流速-比降组合，从而获取多种工况下的急滩通航水力指标；张鹏等[5]结合船舶在洪水急流滩自航上滩的实测资料，确定了三峡库区洪水急流滩河段船舶自航上滩的水力指标；李庆宁等[6]研究了三峡-葛洲坝航段航行的3种船型在泄洪期安全上水过滩的通航能力；童思陈等[7]以山区河流澜沧江为代表，提出了适合于澜沧江代表船舶的水流阻力计算式，并且结合坡降阻力计算, 建立了澜沧江急流滩的通航水力指标。模型试验方面，曹民雄等[8]进行了多种水力条件下的船模上滩试验，发现水面比降、水流速度及船舶水深吃水比是船舶上滩的主要影响因素；刘晓菲等[9]开展了乌江小幺滩河段航道整治物理模型试验，通过调整河段流速和比降，改善水流条件，以满足船舶安全航行要求。

截至目前，国内外学者对涉及船舶航行于比降水面的问题主要集中在评估水力特性和通航指标上，对航行阻力的计算采用的是经验公式或模型试验。笔者采用计算流体力学(CFD)方法，基于商业软件STAR-CCM+，以某型集装箱船为研究对象，通过改变河道断面纵向收缩角和入口流速获得不同的水面比降工况，分析了水面比降对船舶航行阻力与浮态的影响，为船舶航行于带有比降水面时的安全决策提供一定的参考。

1 船舶航行于比降水面阻力分析

船舶航行于比降水面时，由于水面坡度的影响，船舶的浮态及阻力都会区别于一般水面。如图1，根据文献[5]～文献[7]的分析，船舶处于比降水面临界状态时，会产生纵倾，船舶重力沿船体纵向的分力被定义为比降阻力。此时船舶航行的总阻力和比降阻力的经验计算公式为：

R=Rv+Rj

(1)

式中：Rv为水流阻力，N；Rj为比降阻力，N。

比降阻力Rj的计算公式为：

Rj=βWJ

(2)

式中：W为船舶排水量，N；J为船长范围内平均水面比降；β为考虑船舶上滩时水面比降局部增大的修正系数，一般取β=1.1～1.2。

图1 船舶航行于比降水面Fig. 1 The ship sailing on the water surface gradient

船舶航行时，由于水面存在一定比降，船首尾的水位差导致船舶倾斜。船倾斜后，浮心位置会发生偏移，致使浮心和重心不在同一铅垂线上，于是作用方向相反的重力和浮力会促使船舶纵倾角增大，增大后的纵倾角设为θ′，θ′理论上要大于比降对应的水面坡度角θ。因此，实际上船舶重力在船体纵向的分力应该为Wsinθ′，这才是比降阻力。由此可知，采用式(2)计算得到的比降阻力会小于实际值，从而导致利用经验公式计算得到的船舶航行于比降水面时的总阻力会偏小。

2 数值方法

2.1 控制方程

在空间固定坐标系下，数值模拟采用的是不可压缩流动RANS方程，控制方程为连续性方程和雷诺平均方程，如式(3)和式(4)：

(3)

(4)

2.2 湍流模型

考虑计算时间和计算精度，选用Realizablek-ε湍流模型，湍流动能k及湍流耗散率ε方程如式(5)和式(6)：

(5)

(6)

式中：k为湍流动能；ε为湍流耗散率；μt为湍流黏性系数；C1=1.44；C2=1.9;σk=1.0；σε=1.2。其中μt为：

(7)

式中：Cμ=0.09。

3 计算模型与验证

3.1 计算模型

以国际研讨会标准集装箱船KCS为研究对象，如图2。船模主要参数为：垂线间长Lpp=7.278 6 m，型宽B=1.019 m，型深D=0.341 8 m，设计吃水T=0.34 m，湿表面积S=9.438 m2，方形系数Cb=0.65，模型缩尺比λ=31.67。

KCS虽为海船船型，但其船型的长、宽、吃水与长江过闸船型集-18[10]的相应参数基本呈2∶1的尺度比。虽然KCS船型的方形系数0.65比集-18的方向系数0.53大，但根据船舶原理[11]中关于方形系数对船舶阻力影响的叙述可知，在低于1.93 m/s的航速范围内，方形系数差异造成的船型阻力变化仅为10%左右。

图2 KCS船几何模型Fig. 2 KCS ship geometric model

3.2 船舶浮态与阻力之间关系的验证

船舶在比降水面航行时会由于水面坡度而产生不同程度的尾倾，进而造成船舶阻力的变化。为了验证文中数值模型的合理性，将船舶在不同初始纵倾角下的阻力数值计算结果与文献[12]进行了对比。

表1 网格基本参数信息Table 1 Grid basic parameter information

表2 船模纵倾阻力数值计算与试验对比Table 2 Comparison between numerical calculation and experiment of trim resistance of ship model

由表2可以发现，三套网格的计算结果与试验的拟合度较好，最大误差仅为3.21%。网格1在纵倾角小于0.5°后，与试验结果存在一定偏离，特别是在纵倾角为0.25°和-0.5°时的误差较大；网格3在纵倾角0.25°时误差最大，大于0.25°时的计算结果较网格2偏差也较大。在综合考虑计算精度与时间成本后，选择网格2的属性设置参数作为数值计算的网格划分参数。

3.3 计算工况设置

计算水深设定为8 m，水深吃水比为h/T=23.5。通过改变河道断面纵向收缩角和入口流速V来获得不同的水面比降。河道断面纵向收缩角是边坡与x方向的夹角α，如图3。笔者设置夹角α分别为1°、2°、3°、5°，入口流速V分别为0.5、1.25、2 m/s，以α-V来表示对应工况。数值计算一共设置3组工况，第一组工况：1-0.5、2-0.5、3-0.5、5-0.5；第二组工况：1-1.25、2-1.25、3-1.25、5-1.25；第三组工况：1-2、2-2、3-2、5-2。

船体坐标系如图4，原点O位于船模重心，x轴指向船首，y轴指向左侧岸壁，z轴为垂直水面向上。船模首倾时纵倾角为正值，尾倾时纵倾角为负值。

图3 计算域平面Fig. 3 Plan of calculation domain

图4 船体坐标系Fig. 4 Ship hull coordinate system

3.4 边界条件及网格划分

假定船舶沿航道中线航行，且航道关于航道中线对称，因此，航道和船舶均可采用半模进行设置和计算。边界条件的设置如图5，计算域的进口、顶部及底部均设置为速度进口，出口设置为压力出口，侧边及船体表面设置为壁面，计算域中纵剖面设置为对称平面。计算域中，进口距船首2Lpp，出口距船尾3Lpp，出口的宽度固定为0.7Lpp，边坡长度固定为3.3Lpp。设定船模距岸壁的距离超过一倍船长，水深吃水比h/T=23.5，认为船舶不受岸壁和浅水效应的影响，只在带有比降水面的敞水水域航行。数值计算放开了船体升沉和纵倾两个自由度；求解时间离散采用一阶；扩散项采用二阶中心差分格式离散；对流项采用二阶迎风格式，以提高计算精度；压力-速度耦合方程采用SIMPLEC算法进行求解。

网格采用切割体网格，对自由液面及船体周围进行体加密，如图6。船体首尾型线曲度较大，需进行体加密，如图7。经过多次试算，最终确定无量纲化的壁面距离y+值的范围在30～60之间，说明船体边界层的网格设置是合理的。

图5 计算域边界条件Fig. 5 Boundary conditions of calculation domain

图6 计算域整体网格Fig. 6 Global grid of calculation domain

图7 船首尾网格Fig. 7 Bow and stern grid

4 结果分析

4.1 阻力分析

表3汇总了3组工况计算的阻力值，其中水流阻力作为对比项，是在提取到对应比降工况下船体附近平均流速数值后，另外计算的船舶在敞水中的阻力值。比降阻力的数值模拟值由总阻力与水流阻力的差值得到，比降阻力的经验公式值由式(2)求得，并且也给出了经验公式值与数值模拟值的差值百分比。从表3可知：每组工况中，随着比降的增加，船舶的总阻力逐渐增大，意味着航行难度增加；当入口流速一定，每组工况的水流阻力相差并不大，水流阻力主要由船速决定，或者说是固定船模数值计算中的入口水流流速决定。在每组工况中，船舶周围的流速与入口流速相比变化也不大，这说明式(1)将船舶在有比降水面上航行时的总阻力分解为水流阻力和比降阻力两个部分的思路是合理的。

表3给出了3组工况的比降阻力数值模拟值和经验公式值随水面比降变化的对比。由表3可以发现，采用经验公式计算的比降阻力明显偏小，这与第1节理论分析得出的结论一致。而且，在每组工况中，随着比降的增大，相对误差又在缩小，说明式(2)在大比降时计算准确度高，而并不适用于小比降时的情况。究其原因，应该是因为比降较小时，船舶重力浮力力矩造成的船舶纵倾角调整相对影响较大，而在大比降时，由于比降数值大，因此船舶重力浮力力矩造成的船舶纵倾角调整相对较小，比降阻力的数值趋向于直接用比降坡度角计算的结果。

图8是3组工况的摩擦阻力Rf和压阻力Rp随不同比降变化的曲线。

图8 不同比降下的摩擦阻力和压阻力曲线Fig. 8 The frictional resistance and pressure resistance curves under different gradient

从图8中可以看出，阻力的主要成分是压阻力Rp(粘压和兴波阻力之和)，且压阻力Rp会随着水面比降的增大而明显增大；摩擦阻力Rf随水面比降的变化较小，在高流速工况下，甚至有一定的减小。摩擦阻力Rf与船体附近流速以及船模的湿表面积Sm有关。表4给出了3组工况船模在不同比降下的湿表面积Sm。由表4可以发现：每组工况内，湿表面积Sm变化很小；工况之间，由于入口流速变大，船模湿表面积Sm有一定的增加，但增加数值仍然很小。由此可知，船舶在有比降水面航行时，船舶摩擦阻力Rf的变化主要由船体附近流速决定，而与湿表面积Sm的变化关系不大。

表4 不同比降下KCS船模湿表面积Table 4 The wet surface area of KCS ship model under different gradient

图9展示了第一组工况船体底部动压分布情况，从图9中可知：比降为0.1‰ 时，低压区覆盖船底首部到尾部的全部区域；比降为0.2‰ 时，船首部压力增大，船底低压区缩减到仅覆盖尾部区域；比降为0.4‰ 时，船底低压区进一步缩减到仅覆盖尾肩部；比降为0.6‰ 时，虽然流速几乎没变，但比降的增大仍然使得船底的压力数值整体增大。由此可见，比降的增大，不但会增大船体首尾的压力差，也会增大船体所受水压力的数值。

图9 船体底部压力随比降变化分布Fig. 9 Distribution of hull bottom pressure changing with gradient

4.2 自由液面波形和船舶浮态分析

图10为第三组工况的自由液面波形等值线。图11为3组工况距船中纵剖面1倍的船宽B处波高H沿船长X方向的分布。由图10、图11可知：每组工况下，随着比降的增大，波面的平均倾角会增大。3组工况间的对比规律为：随着船速的提高，船体兴波对水面的扰动作用逐渐显著，船首壅水导致首波峰的幅值逐渐增大，船尾后部的波浪起伏也逐渐明显。从能量的角度看，波浪起伏越大，需要消耗越多的能量，这一部分能量最终以阻力的形式体现，这也是船舶航行于比降水面时阻力增大的原因之一。

图10 自由液面波形等值线Fig. 10 Free surface waveform contour

图12是3组工况下船体纵倾角随水面比降变化的曲线，并且给出了水面坡度角作为对比，因为尾倾时船体纵倾角为负值，为了便于比较，图12中的水面坡度角取负值。由图12可发现，船舶在上滩航行时船身会产生尾倾现象。随着水面比降的增大，船身尾倾更严重，并且尾倾角度大于水面坡度角，这也进一步证明采用式(2)计算出的比降阻力值会偏小。

图13是不同流速下，船尾自由面随水面比降变化分布，数值仿真中采用VOF方法捕捉船体自由液面。从图13中可以明显看到，随着水面比降增大，自由液面高度在逐渐抬升，说明船体尾倾角度在增大。

图11 不同比降下的波高分布Fig. 11 Wave elevation distribution under different gradient

图12 纵倾随比降的变化曲线Fig. 12 Variation curves of trim changing with gradient

图13 船尾自由面随比降变化分布Fig. 13 Distribution of free surface of stern changing with gradient

4.3 流场特性分析

桨盘处伴流场是计算螺旋桨推力及主机功率的重要参数，因此需研究水面比降对桨盘处伴流场的影响[14]。单桨船的桨盘面关于中纵剖面左右对称，为方便对比，采用图14进行分析，图14中Wx为桨盘处标称伴流分数。由图14可知：图14(a)中，桨盘面中心处的高伴流区随着比降的增大呈现扩大的趋势，桨盘面下部的伴流随着比降的增大有向外扩展的趋势，即随着比降的增大，桨盘处的平均伴流增大；图14(b)中，桨盘面下部的伴流随着比降增大仍呈现向外扩展的趋势，但上部水流却随着比降的增大呈收缩趋势，说明高伴流区在向下移动，伴流的变化必将对螺旋桨的性能造成影响。

图14 第三组工况下标称伴流分数Wx随比降变化的分布Fig. 14 Distribution of nominal wake fraction changing with gradient in the third group of working conditions

图15～图17分别展示了第三组工况船底、船首、船尾涡量细节分布，涡量显示根据Q准则得到。从图15中可以看到船首舭涡、侧涡、肩涡和船尾舭涡、船尾涡等紊乱的涡流。对比4个工况下的船底涡量分布可以发现，随着水面比降的增大，船底涡量强度会逐渐减弱。水面比降为1.7‰ 时，船底涡量强度最高，船首、尾舭涡附着在船底；水面比降增大至7.4‰ 时，舭涡涡量减少，强度也变弱。从图16可以发现，水面比降为1.7‰ 时，船首侧涡和肩涡清晰可见，并且呈现融合状态；比降增大至3.3‰ 时，侧涡、肩涡分离，且侧涡强度已开始减弱；比降增大至4.6‰ 时，肩涡强度开始减弱；当水面比降为7.4‰ 时，船首侧涡分离成两小段，两侧翼状结构消失。从图17可以发现，水面比降增加时，船尾涡的强度也会减弱，水面比降增大到4.6‰ 时，船尾涡的尾翼状结构消失。

图15 船底涡系分布Fig. 15 Diseribution of vontex system at bottem of hull

图16 船首涡量Fig. 16 Bow vorticity

图17 船尾涡量Fig. 17 Stern vorticity

5 结 论

笔者以实验船模为研究对象，对船舶航行于有比降水面进行了数值模拟，分析了船舶的阻力、兴波、浮态和黏性流场的变化。得出了以下主要结论：

1)计算流体力学(CFD)方法可以用来模拟船舶遭遇有比降水面航行时的阻力变化和航行特性，评估船舶航行能力。

2)船舶航行在带有比降的水面时，水面比降越大，船舶所受到的总阻力也越大，增加了航行难度；水面比降对摩擦阻力的影响小，对压阻力的影响较大。

3)船舶在有比降水面航行会产生一定量的尾倾。随着水面比降的增大，船身尾倾会更严重，并且尾倾角度会大于比降对应水面坡度角，这将导致比降阻力的真实结果比经验公式计算出的数值大。

4)船舶航行在带有比降的水面时的流场特性表现为：船尾桨盘处的伴流分数会随着水面比降的增大而增大，影响螺旋桨的运转效率，增加航行难度；水面比降的变化会影响船体涡量强度和涡系的发展，涡量强度随着水面比降的增大会减弱。