0 引　言

内河航运凭借运量大、成本低及低碳排放优势，已成为现代航运网络的关键组成部分，尤其在我国，其运输成本仅为铁路的1/3、公路的1/9，经济与环境效益显著^[1]。然而，受航道水深不足、断面狭窄及岸壁结构复杂等条件制约，船舶航行时阻力显著增加、航速受限，严重制约通航效率并加剧能源消耗。

针对限制性航道船舶的阻力特性，国内外学者已展开了广泛的研究。Elshcrbiny等^[2]结合实验、数值、分析及经验方法，研究了苏伊士运河平静水域的整体情况，揭示了运河深度和宽度对船舶航行性能的强耦合影响；Du等^[3]基于数值模拟，计算了2艘内河船舶在全封闭水道中的阻力，并通过拖曳水槽实验进行验证，研究了航道尺寸和船舶吃水对航速的影响；Hadi等^[4]利用计算流体力学（Computational Fluid Dynamics，CFD）方法和试验，重点分析了梯形截面与矩形截面航道宽度和水深对船舶阻力的影响；Terziev等^[5]通过CFD方法研究了航道水深突变对船舶水动力性能的影响，结果表明，当船舶进入更浅的水区时，航道底部会形成强边界层，且在临界速度附近，水深变化引起的阻力增加可达初始值的2倍；易信宇等^[6]以乌江航道为例，基于乌江标准船型开展了阻力与自航试验，分析了不同水深条件下的阻力与自航性能，同时比较了不同动力方案，最终推荐采用柴油-天然气双燃料模式；徐双喜等^[7]结合数值计算与模型试验，以漓江航道反铲式挖泥船为例，通过计算下沉量对航行阻力的修正，较准确地估算了浅水阻力；俞中奇等^[8]基于试验数据提出了适用于限制性Ⅲ级航道的阻力估算公式，吴俊等^[9]则对该公式进行了修正，研究了中间渠道船舶的航行阻力；冯伟^[10]探讨了内河I级航道船舶阻力与中间渠道断面系数的关系；肖应彪^[11]结合试验与数值模拟研究了乌江通航隧洞与船舶的关联性，并建立了隧道尺度与船舶参数之间的关联性模型。

从数值模拟到模型试验，或是将二者结合的研究方法，国内外学者已开展了大量研究。然而，大部分研究主要基于模型尺度，通过二因次法或三因次法将阻力外推至实船时，可能会导致实船阻力出现偏差。为提高计算精度，本文采用CFD方法对实尺度下的集装箱船绕流场进行模拟，以期获得更加准确的阻力结果。同时，基于数值结果分析实尺度条件下不同航道断面形状对航行阻力、周围流场特性、自由液面兴波及船体压力分布的影响。此外，考虑到研究对象主要适用于断面尺寸较小的运河或水系，本文所研究的航速均处于低亚临界航速范围，以更贴近实际航行条件。

1 数值方法

本文采用雷诺平均纳维-斯托克斯（Reynolds-Averaged Navier-Stokes，RANS）方法求解粘性不可压缩流场，其控制方程表达如下：

$ \frac{\partial u_i}{\partial x_i}=0，$

(1)

${ \dfrac{\partial u_i}{\partial t} + \dfrac{\mathrm{\partial}}{\partial x_j}(\rho u_iu_j) = - \dfrac{1}{\rho}\dfrac{\partial p}{\partial x_j} + \dfrac{1}{\rho}\dfrac{\partial\mathrm{ }}{\partial x_j}\left(\mu\dfrac{\partial u_i}{\partial x_j} - \rho\overline{\mu_i^{'}\mu_j^{'}}\right) + S_i}。$

(2)

式中：${u_i}$和$ {u_j} $为速度分量的时均值；${x_i}$和${x_j}$为空间坐标，$u_i^{'}$和$u_j^{'}$为速度分量的脉动值，$ \rho \overline {\mu _i^{'}\mu _j^{'}} $为雷诺应力，$i,j = 1,2,3$；$p$为时均压力；$\rho $为流体密度；$\mu $为动力黏性系数；${S_i}$为源项，需要引入湍流方程来封闭RANS方程。

本文采用可实现的湍动能-耗散率（Realizable $k - \varepsilon $）湍流模型，其数学表达式为：

$k$方程：

$ \frac{\partial(\rho k)}{\partial t} + \frac{\partial(\rho ku_j)}{\partial x_j} = \frac{\partial}{\partial x_j}\left[\left(\mu + \frac{\mu_t}{\sigma_k}\right)\frac{\partial k}{\partial x_j}\right] + G_k + G_b - \rho\varepsilon - Y_M。$

(3)

$\varepsilon $方程：

$ \begin{split} & \frac{\partial(\rho\varepsilon)}{\partial t}+\frac{\partial(\rho\varepsilon u_j)}{\partial x_j}=\frac{\partial\mathrm{ }}{\mathrm{ }\partial x_j}\left[\left(\mu+\frac{\mu_t}{\sigma_{\varepsilon}}\right)\frac{\partial\varepsilon}{\partial x_i}\right]+ \\ & \quad\quad\rho C_1\overline{S}\varepsilon-\frac{\rho C_2\varepsilon^2}{k+\sqrt{v\varepsilon}}+\frac{C_{\varepsilon1}C_{\varepsilon3}G_b\varepsilon}{k}。\end{split} $

(4)

式中：$k$为湍动能；${\mu _t}$为湍流黏度；${G_b}$为浮力引起的湍动能衍生项；${Y_M}$为可压缩湍流扩散产生的波动；$\overline S $为平均应变率；${\sigma _k}、{\sigma _\varepsilon }、{C_1}、{C_2}、{C_{\varepsilon 1}}$和${C_{\varepsilon 3}}$均为常数。

近壁面流动采用STAR-CCM+软件中的all-y⁺壁面处理方法，根据y+值的大小自动选择处理方式，当y⁺<1时直接求解粘性底层流场，当y⁺>30时使用壁面函数求解边界层流场，当1<y⁺<30时使用混合函数求解边界层流动，这样可减小近壁面区域流动模拟对网格的依赖程度。船舶航行时产生的自由面兴波用流体体积（Volume of Fluid，VOF）法捕捉^[12]，同时，采用高分辨率界面捕捉技术（High-Resolution Interface Capturing，HRIC）来提高自由面的模拟精度，用以保证体积分数求解的二阶精度^[13]。

2 研究对象及工况

本文以1艘96 TEU内河集装箱船为研究对象，其三维效果图如图1所示，船型的主尺度参数见表1。

分别计算实尺度下，该船位于航道1与航道2不同工况下的航行阻力，并分析其周围流场、兴波和船体压力的变化特征。航道1与航道2简化后的断面形状如图2和图3所示，其中H为航道水深，B_c为航道水面半宽，y_s为船体一侧靠近岸壁的距离。本文设置了4种水深和3种航速条件，总共计算了24个算例，具体的航行计算工况详见表2，其中船长傅汝德数$Fr = U/\sqrt {g{L_{{\text{pp}}}}} $，水深傅汝德数$F{r_h} = U/\sqrt {gH} $，雷诺数$Re = U{L_{{\text{pp}}}}/\nu $。

3 计算域及网格划分

本文采用的直角坐标系O-xyz固定在船基平面和该船体的中纵剖面处（X_船尾柱、Y_船中线面、Z_船基平面），x轴指向船首，y轴指向左舷，z轴向上。在数值模拟中，采用船静止，水流移动的方法进行仿真计算；同时，由于本文研究的航道为对称航道，因此采用在y=0～Y_max处的半流域进行仿真计算。以H/T表示水深与吃水之比，y_s/B表示集装箱船右侧距岸壁距离与船宽之比。当H为5 m时，航道水深最大，H/T为1.78（小于2），集装箱船处于浅水状态，所以航道底部采用有航速无滑移的边界条件；同时，当集装箱船行驶在航道2时，y_s/B为2.10（大于2），岸壁采用滑移边界条件，而对于航道1而言，y_s/B为1.67（小于2），故岸壁采用有航速无滑移的边界条件^[14]。其计算域空间大小和边界条件如图4和表3所示。对入口边界与出口边界应用强制消波和阻尼消波2种消波方式，在航道边界前后600 m范围内进行消波，以消除较远边界处波浪反射的影响。

针对船体进行实尺度的仿真计算，网格基本尺寸设置为1 m。在边界层网格加密方面，总厚度设定为0.1 m，划分为10层棱柱层，增长率为1.3。对自由液面、开尔文波区域以及湍流区域分别进行加密，总网格数为240万左右，船体周围网格如图5和图6所示。根据国际船模拖曳水池会议（International Towing Tank Conference，ITTC）建议，在船舶阻力和粘性绕流的RANS计算中，时间步长推荐为小于0.01L_pp/U的值^[15]，故时间步长取0.2 s。

4 网格及时间步长收敛性分析

本文选择航道1在水深4.5 m，航速9 km/h的工况进行网格收敛性和时间步长分析。网格收敛性研究中，采用5套加密的网格进行数值计算，不同网格尺度的计算结果如表4所示，用h_i/h₁表示网格加密比，其中下标“1”代表最密的网格，h_i+1/h_i=$\sqrt[4]{2} $。时间步长取0.2 s。图7给出了阻力随网格加密的收敛趋势。

将5套网格计算得到的阻力用最小二乘曲线拟合的方法得到数值计算精度p，网格计算误差${\delta _{RE}}$以及不确定度${U_G}$^[16]。如表5所示，可以看到阻力的数值收敛精度较大，且不同网格尺寸计算结果的网格计算误差与不确定度相近，说明若数值振荡比较显著，数值不确定度也较大，用最小二乘曲线拟合可以较好地计入振荡数据对网格收敛的影响。

时间步长收敛性研究中，时间步长分别设定为0.4、0.28、0.2、0.14、0.1 s，使用验证的网格3作为本次时间步长收敛性验证的计算网格。不同时间步长的计算结果如表6所示，用${{{t_i}} / {{t_1}}}$表示时间步长加密比，其中下标‘1’代表时间步长最小，$ {{{t_{i + 1}}} / {{t_i} = \sqrt 2 }} $。图8给出了阻力随时间步长加密的收敛趋势。

将5种不同时间步长计算得到的阻力用最小二乘曲线拟合的方法得到数值计算精度p，网格计算误差${\delta _{RE}}$以及不确定度${U_G}$^[16]。如表7所示，虽然阻力的数值收敛精度较小，但不同时间步长计算结果的网格计算误差与不确定度相近，并且随时间步长减小，阻力差异并不明显。因此，综合考虑计算效率和计算精度，在后续计算中网格数量设为240万左右，时间步长设为0.2 s是合理的。

5 数值计算结果 5.1 计算结果验证

为了验证仿真结果的有效性，本文采用兹万科夫水流阻力经验公式对航行阻力进行初步估算^[17]，如式（5）所示，此法系根据苏联伏尔加河非自航驳船试验结果，并分析了其它近似公式而得出，适用于自航和非自航的内河船舶水流阻力的粗略估算。

$\left\{ \begin{split} & R_v=g(f_1A_sV_s^{1.83}+\xi_1\delta A_mV_s^{1.7+4F_r})，\\ &\xi_1=\dfrac{17.7\delta^{2.5}}{\left[\dfrac{L_W}{6B}\right]^3+2}。\end{split}\right. $

(5)

式中：${f_1}$为机动船摩阻系数，钢质船${f_1}$=0.17；${\xi _1}$为剩余阻力系数；${A_s}$为船舶湿表面积, m²；${V_s}$为船舶对水的速度，其值为船舶表面流速和对岸航速之和；${A_m}$为船舶浸水部分舯剖面积, m²；${F_r}$为傅汝德数；${L_W}$船舶水线长，m；$g$为重力加速度，取$ g= 9.81\; \mathrm{m/s}^2 $基于式(5)对6 ～12 km/h区间内的部分航速估算了对应的阻力值，并与数值模拟得到结果进行了对比，结果如表8所示。

由表8可知，在航道1中，数值计算所得的阻力与经验公式计算结果的相对误差介于1.0%～25.2%，而在航道2中，该误差为2.2%～17.9%。这种误差主要受到航道水深、宽度和船型等因素的影响，而兹万科夫水流阻力经验公式未考虑这些因素，因此在水深较小或航速较大的极限工况下，可能存在一定的适用性偏差，从而导致误差较大。结合网格及步长的收敛性分析以及与经验公式的对比，本文的计算结果具有一定的可靠性。

5.2 阻力特性分析

由图9可知，阻力随着水深的增加呈下降趋势。在水深4 m时，对于航道1，船体在航速12 km/h时所受阻力为6 km/h时的5.22倍；当水深增至4.5 m，该比值降至4.87；进一步增加至5 m时，比值进一步降低至4.70。这表明，不同航速之间的阻力差异随水深增大而减小，浅水效应的影响逐渐减弱。此外，航道1的航行阻力普遍高于航道2，且该差异在浅水深时更为显著，随着水深增加整体呈现减小的趋势。在航速为6 km/h和9 km/h时，由于航速较低，其水深弗劳德数较小，由水深和航道宽度导致的阻力差异并不显著；而当航速提升至12 km/h时，两航道间的航行阻力差异普遍超过10%。

当航速为6 km/h时，航道2在水深4.5 ~5 m之间的航行阻力变化不明显，而航道1的航行阻力则有所降低。造成这一现象的原因可能与低航速下绕流速度较小有关，使得水深变化对阻力的影响不如高速时显著。此外，航道1的宽度小于航道2，在相同的水深变化下，航道1受到更强的岸壁效应影响，这可能使得在低航速下其船体绕流速度变化更大。相比之下，航道2更宽，岸壁效应较弱，水深变化对其阻力的影响趋于减小。

5.3 船体绕流分析

由图10可知，随着水深的增加，船体周围的绕流速度逐渐减小。由于航道1的宽度小于航道2，相同水深下，航道1的岸壁效应更加显著，因此船体周围的绕流速度普遍高于航道2，并且在水深变化相同时，航道1中绕流速度的变化幅度相对较大。

由图11可知，当航速为6 km/h时，水深由4.5 m增加至5 m，其周围绕流速度的变化幅度较小。相比之下，航道1的绕流速度变化大于航道2，这主要受航道宽度的影响。较窄的航道在相同水深变化条件下，其周围绕流速度的变化更为显著。此外，当航速为6 km/h时，处于超低亚临界航速范围，其最大水深弗劳德数仅为0.267，而在水深4.5 m时，水深弗劳德数更低，使得整体绕流速度的变化幅度进一步减小。因此，在水深由4.5 m增加至5 m过程中，航道2的阻力变化非常小，而航道1由于受岸壁影响，阻力变化则相对更为明显。

5.4 自由液面分析

由图12和图13可知，当水深为4 m时，航道1中船侧水面波动最大波谷幅值接近1 m，船体中部水面的波动范围约为–0.2～–0.6 m，周围水面波动较为剧烈。随着水深的增加，船侧水面波动的波峰与波谷幅值逐渐减小，波动周期缩短，波动减缓，从而降低了船体航行时的兴波阻力。此外，在相同水深条件下，航道1中船侧水面波动幅值普遍大于航道2，其周期更长，波幅也更大。

由图14和图15可知，当航速为6 km/h时，水深由4.5 m增加至5 m，船侧水面波动幅值整体降低，波动特性趋于平缓。在水深4.5 m时，船体中部水面的波动范围为–0.02～–0.1 m，而当水深增加至5 m后，波幅范围为–0.03～–0.07 m，且在船体中部至尾部区域，船体兴波逐渐衰减至接近0。这一现象表明，在航速6 km/h的条件下，对于航道2，当水深超过5 m时，船体周围的兴波特性趋于稳定，水深变化对自由液面的扰动影响显著减弱，因而兴波阻力的变化幅度相对较小。

5.5 船体压力分析

由图16可知，随着水深的增加，船体底部的高压区域逐渐扩大，并由船体中部向船首和船尾扩散，其中船体中部的压力最大。在相同水深条件下，航道2中船底高压区域的面积大于航道1，且随着水深增加，航道1的高压区域扩展幅度大于航道2。

一般而言，由于航道1的宽度小于航道2，所以对于同一水深变化时，其航道2横截面积的变化要大于航道1，所以由船底流速变化导致的压力变化也应该大于航道1。然而，由于航道宽度的影响，航道1在4 m水深时船底与船侧速度均大于航道2，并且在水深变化为0.5 m时，其船底高压区域面积增大幅度大于航道2，所以使得航道1的船底压力变化大于航道2。这表明窄航道浅水效应对于水深变化的敏感度更高，其流场和船底压力变化更加显著。

由图17可知，在航道2中，水深从4.0 m增加至4.2 m时，船底压力变化小于4.2～4.5 m的变化。这表明，水深小于4.2 m时，航道2中船体周围的绕流速度变化幅度减小，从而导致航行阻力的增长速率趋缓。

6 结　语

本文采用CFD数值方法，对实尺度下的集装箱船在2种不同断面航道中的水动力性能进行模拟和预报。由于缺乏相关试验数据，首先进行了网格尺度及时间步长收敛性分析，并确定了合理的网格尺度和时间步长，同时将计算结果与经验公式进行对比验证。在此基础上，计算了不同水深条件下两种航道中集装箱船航行时的阻力，分析了实尺度条件下航道断面形状对航行阻力、周围流场、自由液面兴波及船体压力的影响，并进行了对比研究。基于上述分析，得出以下结论：

1）在相同水深条件下，较窄的航道1相比较宽的航道2，船舶航行阻力较高，且阻力差异在水深较浅时尤为明显。随着水深减少，航道阻力普遍上升，但航道1的阻力变化速率较快，表明其浅水效应更加显著。

2）相同水深变化下，航道1船体周围绕流速度变化幅度更大，自由液面兴波幅值更高，船底高压区扩展更快。相比之下，航道2的浅水效应较弱，绕流场和自由液面波动相对平缓，尤其在航速6 km/h且水深超过5 m时，阻力与流场特性趋于稳定，浅水效应明显减弱。

本文仅针对低亚临界航速范围内的情况进行分析，所得结论仅适用于断面尺寸较小的运河或水系，对于更高航速条件不具备普适性。由于计算资源的限制，本文未考虑船体垂荡和纵摇的自由度。因此，未来的研究将考虑浅水效应引起的船体下蹲对航行阻力的影响，以更接近实际情况。