0 引　言

船舶专业化、大型化是未来水路运输发展的趋向，可以更有效的节约运输成本、减少污染排放。但由于天然航道条件限制，船舶的大型化发展受到制约，寻求新方法提高大型船舶的运力、扩展船舶通航隧道以及开建人工水道将会成为未来水路发展的研究方向。

隧道作为一种特殊的受限通道，受到地质条件和经济条件的制约，因船舶通航隧道的复杂性、修建难度大等特点，且在世界范围尚处于起步阶段，通航隧道具有不同于现有通航规范中限制性航道的特点^[1]，当下关于船舶隧道的学术研究及实践理论较为缺乏，而目前对通航隧道建设和管理还没有相关规范和标准，关于通航隧洞航道尺度与通行船舶参数的关联性研究成果也极少。李焱等^[5]采用水工物理模型和船模试验方法研究了构皮滩水电站单线通航隧道内船舶航行的船舶阻力、下沉量以及水面波动特性。汤建宏等^[6]以溪洛渡枢纽的通航隧道为例，结合现有通航规范中船舶通航临界航速的计算公式，提出最经济的隧道断面尺寸系数。郭洪雨^[7]依托富春江通航隧道特点，从结构方面开展了隧道结构验算以及施工方案的分析研究。肖应彪等^[8]针对乌江通航隧洞航道尺度与通行船舶参数的关联性开展了试验研究，探究小断面系数下隧洞航道尺度与通行船舶参数关联性和揭示乌江通航隧洞与千吨级船舶的关联性。李果等^[10]开展长江人工水道与万吨级船舶相关性试验，探讨缩尺比和傅汝德数对与万吨级船舶之间的关联性，提出长江人工水道万吨级船舶操纵安全技术要求建议。

长期以来，船舶领域的流动问题主要是依靠模型试验或试验与理论相结合的方法来解决。但试验研究往往周期长、费用高，并且计算精度不高。计算流体力学是一门用数值计算方法直接求解流动主控方程（Eluer Navier-Stokes 方程）以发现各种流动现象的规律的学科，应用数值计算的方法不仅可以节约设计成本、缩短设计周期、提高经济性，还可以从微观视角观察流动信息，对研究尺度效应、船舶线型优化设计等具有重要的发展和积极的应用，已经成为船舶线型设计及优化的新兴手段，在国内外得到广泛应用^[12]。

1 数值理论 1.1 控制方程

船舶过洞数值模拟主要研究分析船行波波动特性、船舶下沉量特性、船舶航行阻力等船舶航行特征，受到流体粘性影响较大，因此采用粘性流理论方法，根据质量、动量和能量守恒三大物理守恒定律建立控制方程。

引入雷诺平均法，构建雷诺平均方程（RANS方程）对流体的运动进行描述，针对不可压缩的粘性流体，其连续性方程与动量方程分别为：

$ \frac{\partial \rho }{\partial t}+\frac{\partial }{\partial {x}_{i}}\left(\rho {u}_{i}\right)=0 ，$

(1)

$ \frac{\partial\left(\rho u_i\right)}{\partial t}+\frac{\partial}{\partial x_j}\left(\rho u_iu_j\right)=-\frac{\partial p}{\partial x_i}+\frac{\partial\tau_{ij}}{\partial x_j}+F_i。$

(2)

式中：$ i,j=\mathrm{1,2},3 $；$ \rho $为流体密度；$ p $为流体微元上的压力；$ {u}_{i} $、$ {u}_{j} $为速度分量；$ {\tau }_{ij} $为粘性应力张量，介质为牛顿流体时，其定义为

$ {\tau }_{ij}=\mu \left(\frac{\partial {u}_{i}}{\partial {x}_{j}}+\frac{\partial {u}_{j}}{\partial {x}_{i}}\right)-\frac{2}{3}\mu \frac{\partial {u}_{k}}{\partial {x}_{k}}{\delta }_{ij} 。$

(3)

式中：$ \mu $为流体动力粘性系数；$ {\delta }_{ij} $为克罗内克尔符号；当$ i=j $时，$ {\delta }_{ij} $=1，否则$ {\delta }_{ij} $=0。

1.2 湍流模型

RANS方法采取近似方法对各尺度湍动信息进行时均化处理，得到经平均化处理后的流场数据，其对网格的要求相对较低，计算量相对较低，且数值结果满足工程计算的要求，所以该方法在实际工程领域获得了广泛的应用。因此船舶过洞数值模拟采用RANS方法。

RANS方法是通过采用雷诺平均法将各物理量表征为时均值与脉动值的叠加，即为RANS方程组：

$ \frac{\partial \rho }{\partial t}+\frac{\partial }{\partial {x}_{i}}\left(\rho {\bar {u}}_{i}\right)=0 ，$

(4)

$ \frac{\delta \left(\rho {\bar {u}}_{i}\right)}{\delta t} + \frac{\partial }{\partial {x}_{j}} \left(\rho {\bar {u}}_{i}{\bar {u}}_{j}\right) = -\frac{\partial \bar {p}}{\partial {x}_{i}} + \frac{\partial \overline{{\tau }_{ij}}}{\partial {x}_{j}} + \frac{\partial }{\partial {x}_{j}} \left( -\rho \overline{{{u}_{i}}^{{'}}{{u}_{j}}^{{'}}}\right) 。$

(5)

式中：$ i,j=\mathrm{1,2},3 $；$ \rho $为流体密度；$ -\rho \overline{{{u}_{i}}^{{'}}{{u}_{j}}^{{'}}} $为雷诺应力项，由于其为未知量，导致方程无法封闭，因此需进一步引入湍流模型以封闭求解RANS方程。

由于速度脉动产生了雷诺应力导致原始方程不再封闭，这种情况下必须添加新的方程才能使得控制方程封闭可解，此时便需要引入湍流模型。本研究采用Realizable $ k-\varepsilon $模型，该模型适用范围比较广泛，可以模拟一些复杂流动，如边界层流动、旋流、分离流等。Realizable $ k-\varepsilon $模型中关于$ k $和$\varepsilon $的输运方程分别为：

$ \frac{\partial \left(\rho k\right)}{\partial t}+\frac{\partial \left(\rho k{u}_{i}\right)}{\partial {x}_{i}}=\frac{\partial }{\partial {x}_{j}}\left[\left(\mu +\frac{{\mu }_{t}}{{\sigma }_{k}}\right)\frac{\partial k}{\partial {x}_{j}}\right]+{G}_{k}-\rho \varepsilon ，$

(6)

$ \begin{split} & \frac{\partial\left(\rho\varepsilon\right)}{\partial t}+\frac{\partial\left(\rho ku_i\right)}{\partial x_i}=\frac{\partial}{\partial x_j}\left[\left(\mu+\frac{\mu_t}{\sigma_k}\right)\frac{\partial\varepsilon}{\partial x_j}\right]+\rho C_{1\varepsilon}E\varepsilon- \\ & l\rho C_{2\varepsilon}\frac{\varepsilon^2}{k}。\end{split} $

(7)

式中：$ {\sigma }_{k} = 1.0 $；$ {\sigma }_{\varepsilon} = 1.2 $；$ {C}_{2} = 1.9 $；$ {C}_{1} = \mathit{\rm max}\left(0.43,\dfrac{\eta }{\eta + 5}\right) $；$ \eta=\left(2E_{ij}\cdot E_{ij}\right)^{\frac{1}{2}}\dfrac{k}{\varepsilon} $；$ E_{ij}=\dfrac{1}{2}\left(\dfrac{\partial u_i}{\partial x_j}+\dfrac{\partial u_j}{\partial x_i}\right) $；$ \mu_t=\dfrac{\rho C_{\mu}k^2}{\varepsilon} $；$ {C}_{\mu }={1}/\left({{A}_{0}+{A}_{S}\dfrac{k{U}^{*}}{\varepsilon}}\right) $；$ {U}^{*}=\sqrt{{S}_{ij}{S}_{ij}+{\widetilde{\mathrm{\Omega }}}_{ij}{\widetilde{\mathrm{\Omega }}}_{ij}} $；$ {\widetilde{\mathrm{\Omega }}}_{ij}= {\mathrm{\Omega }}_{ij}-2{\varepsilon}_{ijk}{\omega }_{k} $；$ {\mathrm{\Omega }}_{ij}= {\bar{\mathrm{\Omega }}}_{ij}-{\varepsilon}_{ijk}{\omega }_{k} $。

2 计算模型及网格 2.1 模型设置

设计代表船型为130 m×19.2 m×5.5 m（船长×船宽×吃水）。按照1∶25的缩尺比创建船模，船模具体形状如图1所示，船型主尺度参数见表1。

2.2 计算域及边界条件设置

本文采用重叠网格的方法需要有2套网格，因此需要对整个计算域进行划分，分为内、外两域，包裹船体四周的区域为内域，其余的背景域为外域。

为了保证来流均匀，外域入口设置在距离船体1.5 L处，其边界条件设为速度入口，并根据研究工况设置好通航隧洞入口处的流动速度；将外域的出口设置在距离船体3L处，其边界条件设为压力出口。将压力出口选在离船舶足够远的下游，是为了使其对船舶周围的流场无干扰和影响，又能保证尾流充分发展。

外域的顶部、底部、岸壁乃至船体表面的边界条件均设为无滑移壁面条件，在此边界上，壁面不能被穿透，且流体不能滑动。外域的各边界面与船体间的距离与确定的模拟工况保持一致。

内域各面中除贴近岸壁一侧的表面设为无滑移壁面条件、与船中纵剖面重合的一面设为对称面外，其余4个面的边界条件均设为重叠网格。除此之外，为避免船行波对计算域入口以及出口边界条件的影响，采用源项阻尼消波的方式，分别在计算域入口和出口处沿船长方向设置一段消波区。设置完成后的船模边界条件如图2所示。

2.3 网格划分

本文对流场中压力梯度较大、流场细节较复杂的区域都进行了加密处理。

船体航行时，船体扰动自由液面会形成由横波和散波波系组成的Kelvin 波，故在垂向上对波高方向进行加密；船体首尾部分压力梯度较大，对阻力影响很敏感，且船尾后的漩涡、分离流等需要细密的网格来捕捉。因而本研究对船首、尾网格进行了加密；近壁面湍流流动从内向外依次是粘性底层、过度层、对数律层。为了区分这些湍流结构距离壁面的距离，引入无量纲长度参数$ {y}^{+} $对近壁区进行处理，边界层$ {y}^{+} $控制在30～300之间，使之符合壁面网格尺寸要求；同时针对船舶航行、六自由度运动等动态问题可以提供高效便捷的处理，本文在进行通行船舶二自由度运动的求解时使用重叠网格来实现网格的耦合。网格划分后的模型如图3所示。

3 试验结果对比

为保证船舶数值模拟计算结果的有效性，对构建的船舶数值模拟模型开展网格收敛性分析，并与船模阻力、下沉量试验结果进行对比验证。

在本文中网格选用中等粗细网格，时间步长取t=0.02 s进行计算，通过迭代残差以及船体受力收敛的情况来判断迭代收敛性。求解结束时残差与初始时的残差相比均下降到3～4个数量级，同时各残差已下降至10⁻⁴数量级或更低。在每一个迭代过程均采用最后的几个振荡周期作为考察对象，将此区间的均值作为最终用于分析的计算结果。

网格收敛性分析时，已初步表明所用数值方法的阻力计算误差较小。除此以外，还针对在拖曳水池进行阻力试验，试验对象与数值模拟保持一致，采用的缩尺比为1∶25。考虑到船模试验为约束模试验（即试验过程中船模固定，浮态不发生变化），为保持一致性数值计算也采用约束模型。

3.1 阻力对比结果

船舶在通航隧洞中航行时，其阻力性能主要受隧洞航道尺度和船舶自身航速的影响。随着航速的提高，船体受到的阻力会增加。将计算结果与试验结果换算到同一模型，计算结果与试验结果对比情况如表2所示。

通过验证，船模阻力数值模拟计算结果与试验结果变化趋势一致，但因船舶实际航行试验时受外界风速以及水流波动等因素影响，阻力结果与仿真模拟结果总体计算误差控制在1.53%～3.41%。

3.2 下沉量对比结果

船舶的下沉量与自身的航速及航道的尺度密切相关。设定与阻力对比试验相同参数，进行换算，换算后数值计算对应的船舶下沉量与船模试验结果对比结果如表3 所示。

在验证工况下，船模下沉量数值模拟计算结果与试验结果变化趋势同样基本一致，计算结果先小于试验结果，随着船舶航速的提升后大于试验结果，主要由于CFD计算中船舶航行于理想状态，下沉量为船舶平均下沉量，在船模试验中，因船舶实际航行试验时受外界风速以及水流波动以及传感器位置等因素影响，下沉量结果与仿真模拟结果总体计算误差控制在2.27%～6.82%，符合CCS发布的《基于计算流体力学（CFD）的内河船舶评估指南》，认为该数值计算方法在阻力求解方面的有效性得到验证，可以用于后续多工况下的研究。

4 数值模拟计算及结果分析

根据通航隧洞最高水位和汛期限制水位要求，选取的航道水深为11 m和41 m，作为典型浅水航行水深，对应模型的水深为0.44 m和1.64 m。通行船舶的航速会直接影响船舶在通航隧洞中的航行性能与姿态，同时通行船舶航速的变化会导致隧洞内流场的变化，选取实船安全航速0.8 ～5.0 m/s范围内6个航速开展数值模拟计算分析，最终确定研究的实船航速为0.8、1.5、2、2.5、3、3.9 m/s，相应数值模拟的船舶航速为0.16、0.3、0.4、0.5、0.6 、0.78 m/s。

综上所述，本研究选取船舶尺度、航道宽度、船舶航速作为影响参数，并以此展开后续关联性数值模拟研究。

4.1 阻力模拟计算

根据通航隧洞实际最高水位和汛期限制水位要求，针对船舶设计代表船型，采取控制变量法，分别从不同的航道宽度与航道水深条件下进行探究，选取实船安全航速0.8～5.0 m/s范围内6个航速开展数值模拟计算分析，求得阻力的具体数值结果如表5 所示，根据表5数据绘制4种工况下的阻力图，见图6。

采用控制变量法，分别研究航道水深和航道宽度不同时船舶阻力变化。

由图7可知，船舶阻力随着船舶航速的提高而增大，但趋势因航道宽度和水深而异。采用控制变量法，分别研究航道水深和航道宽度不同时船舶阻力变化。

当航道宽度一定时，水深作为变量，随着水深的增加船舶阻力在减小，随着航速的增加，在2种航道宽度下均表现为水深越小，船舶阻力增幅越大。在航速增加到0.5 m/s后，水深为0.44 m的船舶阻力出现骤增现象，船舶行驶可能会出现危险。

水深一定时，航道宽度作为变量，当水深为0.44 m时，低速航行时2种航道宽度下的船舶阻力几乎一致，在航速增加到0.5 m/s后，不同航宽下的船舶阻力均出现了骤增情况，但航道宽度较低时的船舶阻力随着航速提高的增幅明显大于航道宽度较大的增幅；当水深增加，为1.64 m时，不同航道宽度下的阻力变化几乎一致。

4.2 下沉量模拟计算

船舶的下沉量与自身的航速及航道的尺度密切相关。数值计算的各工况参数与对应的船舶下沉量及无量纲结果如表6 所示。根据表6数据绘制4个工况下沉量图，见图9。

同样采用控制变量法，分别研究航道水深和航道宽度不同时船舶下沉量变化。

船舶航行于隧洞内时，下沉量随航速的提高而改变，但增幅因航道水深和航道宽度的不同而异。

由图10 可知，在航道宽度不变，水深作为变量时，2种航道宽度下下沉量变化基本一致。水深为1.64 m的船舶下沉量全程随航速的增加而下降，并且在2～3 mm附近小幅度波动；水深为0.44 m的船舶下沉量明显远大于水深为1.64 m的船舶下沉量，同时随着航速的增加，水深为0.44 m的船舶下沉量增加越为明显，在航速达到0.5 m/s后出现骤增现象。

在水深不变的条件下，航道宽度作为变量，水深为0.44 m时，航道宽度越大增幅越为明显，在船舶航速达到 0.5 m/s以后2种航道宽度下的船舶下沉量都显著激增。水深增加后，航道宽度较小的船舶下沉量明显大于航道宽度较大的船舶，初始时随着航速的增加下沉量明显减少，在航速达到0.5 m/s后船舶下沉量开始缓慢增加；航道宽度较大的船舶随着航速的增加下沉量缓慢增加，在航速达到0.3 m/s后下沉量减少，变化趋势与航道宽度较小的船舶基本一致。

5 结　语

本文从 CFD 的基础理论出发，结合CFD 软件STAR-CCM+进行了船舶浅水区域阻力与下沉量的计算与仿真，详述了建模和求解过程，为后续浅水狭窄区域船舶航行提供参考。

根据自航过洞船舶数值模拟计算结果可知：

1）采用 CFD 方法对三峡隧洞中通航船舶水动力特性进行分析，得出了不同航速下阻力与下沉量结果，试验表明CFD 结果与试验值之间误差不超过 4%。所以，采用 CFD 方法得出的数值模拟结果可以满足工程应用要求。

2）由控制变量法，结合CFD计算结果可知，随着航速的增加，水深对于船舶航行影响大于航道宽度。水深越浅，同一航速下船舶阻力和下沉量越大，当航速达到一定程度后，船舶阻力和下沉量都会出现骤增现象；水深增加后，航宽对船舶阻力的变化几乎无影响；航宽越小，船舶初始下沉量越大，随着航速的增加，航宽对船舶下沉量影响也逐渐减小。可根据此方法建议行驶在不同的情况下安全行驶航速。

根据CFD方法，可以准确预测船舶在狭窄浅水区域航行情况，预测航行规律，给出安全性建议。