限制性航道中船行波传播特性的数值研究

王孟飞，邓 斌,2,3，蒋昌波,2,3，伍志元,2,3

(1.长沙理工大学 水利工程学院，湖南 长沙 410114； 2.洞庭湖水环境治理与生态修复湖南省重点实验室，湖南 长沙 410114； 3.水沙科学与水灾害防治湖南省重点实验室，湖南 长沙 410114)

1 研究背景

限制性航道是指因水域狭窄、断面系数小而对船舶航行有明显限制作用的航道。限制性航道中船行波具有波高大、流速急等特点[1]，易破坏航道环境[2]和航道边坡[3]或对其他涉水建筑物造成干扰或损害[4]。目前，我国内河通航标准中关于限制性航道尺度的确定，仅给出双航线下的最小取值(限制性航道断面系数大于6)，航道尺寸的确定仅依据航迹带宽度附加一定安全距离，没有充分考虑航道尺寸与船舶航行条件的影响。为防止船行波对航道及周围船舶产生影响，往往仅规定航道内船舶的限制航速[5]。近年来，随着“一带一路”的深入推进，我国内河航运的任务日益繁重，船舶大型化、快速化的发展特点，使得现有航运条件日益紧张，加大对限制性航道的开发与恢复通航越来越受关注。因此，深入开展限制性航道通航关键技术研究极为必要，特别是限制性水域条件中船行波的相关研究。

目前，国内外很多学者通过现场调查[3,6]、物理实验[7-8]及数值模拟[9-10]等方法对船行波的传播特性[11, 12]、船舶的稳定性[13]以及船行波对周围护岸[3]、泥沙[2]及船舶本身[4]的影响进行了大量的研究。船行波有首波系、尾波系、散波和横波等不同分类形式，通常是随弗汝徳数Fr的不同呈现不同的运动形式[11]。部分学者针对开阔水域条件建立了可描述船行波波高与弗汝德数、离船舷距离、船型等因素的相关关系的数学模型[14-15]。但限制性水域条件下船行波的传播特性不同于开阔水域[16]。因此，越来越多的学者研究限制性水域条件内船行波的相关特性。如：Zhu等[17]解释了浅水条件下单船体与双船体船的首波系与尾波系的不同传播规律。Ji等[1]通过求解Navier-Stokes方程和采用不同的湍流方法分析了不同船型和航速下船行波波高的演化规律，并讨论了船舶尾流的水动力特性。骆婉珍等[18]基于CFD技术应用重叠网格的方法模拟分析狭窄、浅水航道内船行波波形的叠加、反射等现象，但该类方法计算量大，难以求解大船速下船行波的传播演化[2]。

近年来，随着Boussinesq方程不断的改进，如消除了浅水的限制又考虑了波浪的破碎，减小了数值计算量[19]，越来越多的学者运用Boussinesq方程模拟波浪的传播变形等。如：Dam等[16]通过建立二维Boussinesq模型，模拟细长船行波及其传播过程，并得到不同船速下波高的变化规律。Shi等[20]基于完全非线性的Boussinesq模型FUNWAVE-TVD模拟了船行波引起的紊动破碎，并针对两种不同的破碎准则(黏性破碎准则和波高破碎准则)进行对比分析，结果表明两种破碎标准均能准确的模拟船行波在不同Fr下的传播过程。卓明泉等[21]基于COULWAVE开源程序代码实现浅水航道船行波的数值模拟，再现了浅水条件下船首孤立波现象。然而，上述研究大都基于Boussinesq模型下船行波传播特性的研究，较少对船舶吃水以及跨临界航速范围内船行波的特性进行详细描述。

综上所述，尽管现有研究对限制性水域内船行波的传播规律有了一定的认识，但基于限制性水域内船舶航行诱导的船首孤立波波高以及横向流速的研究相对较少。因此，为弥补现有研究的不足，本研究运用完全非线性的Boussinesq模型FUNWAVE-TVD[22-23]，以压力源项近似表示船舶引起自由面的波动，采用Shi等[20]提出的破碎准则，重点讨论限制性水域内船行波的相关特性。

2 数值模型建立

基于完全非线性Boussinesq方程建立了限制性航道内船舶航行的数值模型，具体控制方程如下。

质量守恒方程：

(1)

(2)

动量守恒方程为：

(3)

涡流黏度如下所示：

v=Bδ2(h+η)ηt

(4)

式中：δ为混合长度系数，本研究取δ=1；B为避免数值振荡的常数[19]，本研究取B=1。

船舶模型的实现由压力源确定，即在完全非线性Boussinesq模型中加入压力源项。其中，船舶中心点为(x*,y*)的压力分布可由公式(5)～(7)给出：

(5)

(6)

(7)

数值模型基于三阶Runge-Kutta方法的自适应时间步长和MUSCL-TVD数值解法；干湿边界网格最小深度为0.001 m，CFL数设置为0.25，通量计算采用HLL(Harten Lax & van Leer)构造方法，数值消波技术采用耗散型海绵层消波方法。

3 数值模型工况设置与验证

3.1 数值模型工况设置

参考我国《内河通航标准》(GB 50139-2014)规范的要求，模型计算区域采用Ⅵ级单线航道尺寸，船模采用Ⅵ级航道标准驳船(船长L=32 m，宽B=7.0 m，满载设计吃水p为1.0 m，船舶形状系数α和β均为0.5)，并在计算区域左侧设置消波区域。考虑到船行波的传播在航道内通常呈对称分布，计算在航道一侧布置15个测点，数值模型计算平面布置图如图1所示。

图1 数值模型计算平面布置图(单位：m)

考虑Ⅵ级航道内标准驳船的空载吃水p≈0.4～0.5 m，满载吃水p≈0.9～1.0 m，以及标准航行船速，共设置5种计算工况(A1～A5)，拟定的模型计算工况如表1所示。

表1 模型计算工况设置

3.2 数值模型验证

为验证数值模型计算限制性航道内船行波的可靠性，分别从数值实例和网格收敛性验证两个方面进行详细说明。二维验证算例取自参考文献[7]，实验水槽宽D为3.5 m，水深h为0.114 m；船舶长度L为1.6 m，宽度B为0.4 m，吃水p为0.1 m，船舶的形状系数α和β均为0.5；初始时刻船舶位于x=6 m、y=1.75 m的位置，其他参数设置均与参考文献[7]保持一致。

图1给出了3种工况下测点(x=33.125 m，y=0.72 m)波高历时曲线的对比，结果显示，数值模型计算波高与文献[7]中船行波的变化趋势相吻合。采用skill数验证模型精度要求：

(8)

图2为船行波波高计算值的实验验证结果。计算结果表明各工况下skill数均大于0.90，均方根误差均小于0.038，说明模型能较好地模拟船行波的传播。此外，图2(a)～2(c)黑色方框内为船首孤立波，这是典型浅水航道内船行波的重要特征。

3.3 网格收敛性验证

为探究网格大小对计算结果的影响，确保所有航速下CFL数均小于0.25。选取满载吃水工况A4、Fr=1.2下x=80 m、y=1 m处4种不同大小网格的船行波波高历时曲线进行收敛性验证，验证结果见图3。从图3可以看出，各波高历时曲线基本吻合，其中，以dx×dy=0.2 m2网格大小为参考，网格dx×dy=0.19、0.21、0.22 m2的波高历时曲线的均方根误差分别为0.004，0.004，0.008 m；同一位置处不同网格下船行波最大波高相对误差小于3%。可见，模型在风格大小dx×dy=0.2 m2下已达到计算精度要求。因此，模型计算采用dx×dy=0.2 m2的网格尺度。

图2 不同Fr数船行波波高模拟值的实验验证结果

图3 4种不同大小网格的收敛性验证

4 结果分析与讨论

4.1 波面的演化及最大波高

由于船行波以Fr=1为界呈现不同的传播特性[7]，其中，在跨临界航速(Fr=1)下船行波中横波和散波叠加易形成较大波高。因此，选取工况A3、Fr=1分析船舶运动过程中船行波波面的历时演化过程，结果见图4。由图4可看出，在船舶航行的起始阶段(t=5 s)，船舶前方形成弯曲的孤立波波列面。随着船舶的运动，弯曲的波列面开始变直，船舶前方孤立波数量逐渐增多。与此同时，船舶后方的尾波呈现出开尔文尾流模式，船尾波传播受侧壁反射、叠加等影响，在限制性航道两侧及航道中间各点处(如：x=100 m，t=15 s)易形成较大波高，从而对近岸航道边坡等设施产生不利影响。

为深入分析船首孤立波的传播特征，以A3、Fr=0.75下船行波波高历时曲线为例，定义船行波中船首孤立波和船尾波的最大波高，如图5(a)所示。由于限制性航道内船行波的传播易形成复杂的横向流动，对航道两岸设施及船舶航行造成严重的影响。因此，为分析近岸横向流动的分布规律，定义最大近岸横向流速如图5(b)所示。

图6为工况A3、Fr=0.75下距离近岸1 m处不同位置(x=80, 100, 120, 140, 160 m)的船行波波高历时曲线图。从图6可以看出，在浅水条件下，船首孤立波波高在80 m

图5 最大波高及近岸横向流速的定义

图6 距离近岸1 m处不同位置(x=80, 100, 120, 140, 160 m)的船行波波高历时曲线(A3、Fr=0.75)

图7为工况A3、Fr=0.75下距离近岸不同位置(y=1, 2, 3 m)的船行波波高历时曲线图。可见船尾波波高在最大值出现之后易受航道宽度及边壁影响，且越靠近岸壁船尾波波高越大。另外，随着船首孤立波不断向前传播，船首孤立波逐渐形成新的次孤立波，造成航道水位的上升。

4.2 船行波最大波高分布

基于上述最大船首孤立波波高的定义，图8(a)给出了不同Fr数条件下最大船首孤立波波高随相对吃水的分布规律。可以看出，随着相对吃水的增加，最大船首孤立波波高呈线性规律增加，斜率大致为0.4～0.7。图8(b)给出了不同吃水条件下最大船首孤立波波高随Fr数的变化规律。可以看出，随着Fr的增加，最大船首孤立波波高呈幂函数形式增加，且幂次方约为3.5。

基于上述分析，通过多元非线性回归方程得到了p、Fr与最大船首孤立波波高的经验关系式：

(9)

式中：ηmax为最大船首孤立波波高，m；h为水深，m；Fr为弗汝德数；p为船舶吃水，m。

方程式(9)各系数标准差见表2，多元非线性回归方程拟合优度R2=0.967。为了更加直观地展现多元回归方程的拟合精度，图8(c)给出了数值计算与多元回归方程式(9)的结果对比，可见方程式(9)可以较好地预测不同吃水条件下最大船首孤立波的波高。

表2 船首孤立波回归方程的系数

图9为不同吃水条件下最大船尾波波高随Fr的变化规律。由图9可以看出，当p≤0.7 m时，最大船尾波波高以Fr=0.9为界呈现不同的分布规律。当Fr<0.9时，最大船尾波波高逐渐增加；相反，当Fr>0.9时，最大船尾波波高逐渐减小，这是由于当Fr>0.9时船尾波破碎导致波高减小，这与David等[24]的研究大体一致，但David等并未讨论船舶相对吃水与最大船尾波波高之间的影响关系。然而，当p>0.7 m时，最大船尾波波高随Fr增加逐渐增加，与上述p≤0.7 m下最大船尾波波高呈现不同的分布规律，这是由于当船舶吃水较大时，同一波高下最大船尾波波长相对较长，船尾波不易破碎。

4.3 船行波最大近岸横向速度

影响限制性航道设计的主要因素除最大船行波波高外，还受航道内横向速度大小的影响。图10为限制性航道内不同吃水条件下最大近岸横向流速随Fr的变化规律。从图10中可以看出，受船行波破碎影响，最大近岸横向流速仅在p≤0.7 m、Fr≤0.9下整体呈线性关系，具体为随Fr，p的增加分别呈线性减少和线性增加的趋势。

图7 距离近岸不同位置(y=1, 2, 3 m)的船行波波高历时曲线(A3、Fr=0.75)

图8 最大船首孤立波波高变化规律

图9 不同吃水条件下最大船尾波波高随Fr的变化规律 图10 不同吃水条件下最大近岸横向流速随Fr的变化规律

5 结 论

基于完全非线性的Boussinesq方程建立了限制性航道内船舶运动的数值模型，重点研究了船首孤立波、船尾波以及近岸横向流速在不同p、Fr下的分布规律。与传统开阔水域及深水条件下船舶航行的数值模型相比，该模型能更合理的描述限制性航道中船舶航行过程中波高、流速等分布规律。具体结论如下：

(1) 基于我国《内河通航标准》(GB 50139-2014)规范的要求，采用FUNWAVE-TVD模型建立了精确描述船行波传播及演化的数值模型。研究表明该模型能很好地模拟限制性航道条件下船行波的传播特征。

(2) 本文研究条件下，最大船首孤立波波高与船舶吃水之间呈线性关系，斜率约为0.4～0.7之间；最大船首孤立波波高与Fr之间呈幂函数关系，幂次方约为3.5；此外，基于上述规律与分析建立了限制性水域条件下最大船首孤立波波高与船舶吃水、弗汝德数的幂函数经验公式。

(3) 船尾波波高及近岸横向流速受波浪破碎影响，以p=0.7 m和Fr=0.9为界分别呈现不同的分布规律。其中，Fr主要通过船行波波高影响船行波的破碎，而p主要通过船行波波长影响船行波的破碎。同样，受波浪破碎影响最大近岸横向流速仅在p≤0.7 m、Fr≤0.9内整体呈线性分布规律。

本研究基于完全非线性的Boussinesq方程、压力源项以及黏性破碎准则模拟了限制性航道内船首波系以及船尾波系的传播演化过程，为进一步开展限制性航道中船舶安全通航提供了理论基础，但本研究仅考虑一定船型、航道断面尺寸、静水条件下船行波的演变规律，更复杂的航道条件、船舶运动以及水动力特性之间的内在关系还有待进一步研究。