0 引　言

在全球应对气候变化的背景下，我国提出在2030年前实现“碳达峰”和2060年前实现“碳中和”的奋斗目标^[1]。为推进双碳目标的实现，航运业正面临着向低碳绿色环保化转型和发展的挑战，并对船舶节能减排技术提出了新需求。空气润滑减阻技术是一种通过在船体表面注入气体形成覆盖气层从而降低船舶阻力的新型节能减排技术。开展空气减阻性能研究对船舶节能减排技术的提高以及早日实现“双碳”的战略目标具有重要意义。

1876年，傅汝德最早提出利用气层降低船舶阻力的想法。由于当时技术和理论条件的限制，这一想法未被实现。Hartmut等^[1]开展了气泡减阻机理研究工作，研究发现在物体表面湍流边界层中注入微小气泡后，流体密度降低和湍流修正相互作用，导致减阻现象的产生；Madavan等^[2-3]提出一种混合长度模型，研究表明微气泡在湍流边界层中的作用主要是改变局部流体的有效黏度和密度；Conroy等^[4]提出湍流流体中的能量传递是非高斯分布的，并尝试使用一种递归积分方法来量化气液之间非高斯能量转移；Evseev等^[5]发现边界层内近壁面区域的空隙率是减阻的关键因素；Ferrante等^[6]研究同样表明，在气液两相流中，气泡可以降低湍流强度。

国内，董文才等^[7-8]认为平板喷气减阻的主要原因在于喷气时在平板和水之间所形成的气液混合流改变了平板湍流边界层的流动结构，抑制了湍流流动。陈少峰等^[9]通过梳理EEDI相关要求和计算准则，提出了一种基于EEDI指数的气层减阻技术评价方法，并通过整理国内首艘气层减阻改造船的实测数据进行了EEDI的计算，经过换算，该船装配气层减阻后EEDI减小率达到5.9%，减排效果显著。

本文以某大型滚装船为研究对象，设计一种新型空气减阻装置，采用数值模拟法开展空气减阻机理研究，通过对比通气前后船体表面物理量的动态变化过程揭示空气减阻机理，并基于数值模拟的方法，对比了空气减阻前后船体总阻力的变化情况。

1 基本控制方程

在计算流体力学中，需要遵循三大守恒方程，包括连续性方程、动量方程以及能量方程，由于在船舶阻力性能研究领域以及船舶空气减阻领域一般不存在物质之间的热量或能量交换现象，故在此研究方向一般不考虑能量守恒问题。这三大守恒方程均为对流扩散方程，其中包括非稳态项、对流项、扩散项以及源项，其基本方程形式如下：

连续性方程

$\frac{{\partial \rho }}{{\partial t}} + \frac{{\partial \left( {\rho u} \right)}}{{\partial x}} + \frac{{\partial \left( {\rho v} \right)}}{{\partial y}} + \frac{{\partial \left( {\rho w} \right)}}{{\partial z}} = 0。$

(1)

再利用哈密尔顿算子式可将上式改写为矢量形式：

$\frac{{\partial \rho }}{{\partial t}} + \nabla \cdot \left( {\rho {\mathbf{u}}} \right) = 0。$

(2)

动量方程：

$\begin{gathered} \rho \frac{\rm d{u}}{\rm d{t}} = - \frac{{\partial p}}{{\partial x}} + \frac{{\partial {\tau _{xx}}}}{{\partial x}} + \frac{{\partial {\tau _{yx}}}}{{\partial y}} + \frac{{\partial {\tau _{zx}}}}{{\partial z}} + {S_{{M_x}}} ，\\ \rho \frac{\rm d{v}}{\rm d{t}} = - \frac{{\partial p}}{{\partial y}} + \frac{{\partial {\tau _{xy}}}}{{\partial x}} + \frac{{\partial {\tau _{yy}}}}{{\partial y}} + \frac{{\partial {\tau _{zy}}}}{{\partial z}} + {S_{{M_y}}} ，\\ \rho \frac{\rm d{w}}{\rm d{t}} = - \frac{{\partial p}}{{\partial z}} + \frac{{\partial {\tau _{xz}}}}{{\partial x}} + \frac{{\partial {\tau _{yz}}}}{{\partial y}} + \frac{{\partial {\tau _{zz}}}}{{\partial z}} + {S_{{M_z}}}。\\ \end{gathered}$

(3)

式中：p为压力（表面力）；τ_ij为应力分量（表面力）；S_Mi为体积力分量。

2 船型参数

本文选用7800 m车道货物滚装船作为研究船舶阻力性能的目标船型，采用缩尺模型进行数值模拟研究，缩尺比为1∶32。该船型的主尺度参数如表1所示。7800 m车道货物滚装船采用隐形球艏船首、方尾船尾形式和平底船底形式，滚装船的几何模型如图1所示。

3 空气减阻方案设计

设计一种空气释放装置如图2所示。空气释放装置平齐装配于船底板上，关于船体中纵剖面对称分布于船体两侧，装置整体结构由气室和导流板构成，气体经由导流板引流后，附着于船体外底板处，形成覆盖气膜，实现空气减阻的效果。在实际工程应用中可通过船体内部的通气管道将高压气体输送至气室中，并经由出气壁面释放，释放后的气体经由导流板引流后，附着于船体外底板处，形成覆盖气膜。

数值模拟计算设定中，空气释放装置区域的边界条件需设定为：释放气体的壁面为质量流出面（Mass Flow Inlet，通气速率为3.92×10⁻⁴ kg/s），气室内部壁板和导流板与船体外表面保持一致，设定为壁面，详细的通气装置计算边界设定见图2。本文的阻力性能分析均基于设计航速（20.8 kn）进行，换算为模型速度为1.8914 m/s，喷气出口通气速率均为3.92×10⁻⁴ kg/s。

为了更加清晰地监测喷气过程中船身周围物理量的变化特性，以空气释放装置中纵剖面为中轴线，在减阻装置出气口附近以及船底平板区域均设置了监测点，详细的监测点布置如图3所示。

4 网格收敛性分析

为验证计算网格的收敛性，以设计航速工况为例，采用3套不同尺寸的网格方案（Mesh1～3）对滚装船及其细化区域进行网格划分，3套网格方案中的网格尺寸按网格细化率${r_G}$为$\sqrt 2$进行划分，3套网格方案中的网格尺寸可按细、中、粗尺寸进行区分。3套网格划分方案的网格划分如图4所示。

3套网格方案中的网格总数和总阻力计算结果统计如表2所示。

根据表2中3套网格尺寸方案的计算结果进一步验证网格收敛性。具体验证计算过程如下，用ε表示各总阻力计算结果之间的比较误差，比较误差计算式为：

$\begin{gathered} {\varepsilon _{21}} = {F_2} - {F_1} = 0.103 ，\\ {\varepsilon _{32}} = {F_3} - {F_2} = 0.1576 。\\ \end{gathered}$

(4)

式中：F₁～F₃分别为3套网格方案中滚装船的总阻力计算结果。网格收敛率R_G为：

${R_G} = {\varepsilon _{21}}/{\varepsilon _{32}} = 0.6535。$

(5)

由上式计算所得到的网格收敛率R_G处于区间（0，1）之间，满足网格单调收敛条件，所以3套网格方案为单调收敛，满足网格收敛性。3套网格方案的准确度阶数精度P_G为：

${P_G} = {\text{In(}}{\varepsilon _{32}}/{\varepsilon _{21}}{\text{)/In}}({r_G}) = 1.227\ 2 。$

(6)

修正因子C_G为:

${C_G} = ({r_G}^{{p_{\text{G}}}}{\text{ - 1)/}}({r_G}^2{\text{ - 1)}} = 0.530\ 1 。$

(7)

网格误差δ_REG和改进网格误差δ_G分别为：

$\begin{aligned} &{\delta _{REG}} = {\varepsilon _{21}}/\left( {{r_G}^{{P_G}} - 1} \right) = 0.194\ 3，\\ &{\delta _G} = {C_G} \times {\delta _{REG}} = 0.103 。\end{aligned}$

(8)

数值不确定度U_G和修正后数值不确定度U_GC分别为：

$\begin{aligned} &{U_G} = \left( {2 \times \left| {1 - {C_G}} \right| + 1} \right) \times {\delta _{REG}} = 0.376\ 9 ，\\ &{U_{GC}} = \left| {1 - {C_G}} \right| \times {\delta _{REG}} = 0.091\ 3。\end{aligned}$

(9)

3套网格方案中的网格收敛率R_G满足单调收敛条件，改进的网格误差δ_G和修正后数值不确定度U_GC均小于数值不确定度U_G，本文所采用的网格尺寸方案均满足数值误差和不确定度要求。网格方案2的计算结果与网格方案1的计算结果较为接近，综合考虑数值模拟所需的时间和模拟精度等因素，选取网格方案2进行阻力数值研究。

5 空气减阻机理研究

考虑数值计算稳定性，在流场初步稳定10 s之后船底气层减阻装置开始进行喷气。图5～图6为各监测点在喷气前后的密度和湍流粘度变化曲线图。

图5表明，开启空气润滑减阻系统后，气口中心P₀处的密度急剧下降，由纯液相密度997.561 kg/m³降至纯气相密度1.184 kg/m³，PY₁点处的密度变化与P₀一致，PY₂点处的密度并无显著变化。沿船长方向，PX₁～PX₃点处的密度分别降为：408.902 kg/m³、474.265 kg/m³和574.632 kg/m³。释放空气后，在船体表面会形成包含气体和液体的气液混合流，距离气口越近，混合流的密度越小，随着气体的扩散，逐渐在船体表面形成覆盖气层。气层的存在导致液体与船体表面的直接接触减少，造成船体浸湿面积的降低。

图6表明，开启空气润滑减阻系统后，各监测点处的湍流粘度变化与密度变化趋势基本一致，各监测点处的湍流粘度在喷气后呈现不同程度降低，气液混合流的存在不仅减少了与船体表面接触的浸湿面积，并且造成了流体局部湍流粘度的降低。

为了更直观地体现船体物理量的变化，本文提取稳定状态计算结果进行分析。图7为空气释放装置内外部的流线分布情况，可以看出，气体经导流板引流后通过气口并随着水流流动方向扩散，并且经由气口窄缝的加速作用，位于气口处的气体释放速度最大。

图8分别为船底区域的气相分布图、湍流粘度分布图、动力粘度分布图和剪切力分布图。

可知，湍流粘度和动力粘度的分布与气相分布高度相似，释放空气后，气层所能覆盖到区域的湍流粘度和动力粘度均出现了降低，气层覆盖区域处的剪切力明显低于未能覆盖的区域，气液混合流的存在导致了船体表面剪切力的减少。

为了使减阻机理的分析更加清晰，本文再次从理论分析的角度，揭示阻力降低的原因。流体粘性导致的剪应力计算式为：

$\tau {\text{ = }}\mu \frac{{\partial \bar u}}{{\partial y}} - \rho \overline {u'v'}。$

(10)

式中：$\mu$为流体动力粘性系数；$u'$和$v'$分别为流向和法向上的脉动速度分量；$- \rho \overline {u'v'}$又被称为雷诺应力项。

当气层覆盖的船体区域时，其湍流粘度和动力粘度均会显著降低，导致式(10)中第1项$\mu (\partial \bar u/\partial y)$计算值的降低。释放空气后，船体表面接触介质由纯液体变为气液混合流，降低了流体的密度，密度的降低导致式(10)中第2项雷诺应力的降低，最终导致船体所受剪切力的降低。

经过数值和理论分析，空气润滑减阻基本原理如下：喷气装置在船体释放空气后，船体表面形成覆盖气层和气液混合流，流体密度降低、湍流粘度和动力粘度降低、雷诺应力减小，综合效果导致了船体摩擦阻力降低，实现了润滑减阻的目的。

6 空气减阻性能研究

为了更好地体现空气润滑减阻的性能，本文对比了设计航速下，空气减阻装置开启前后船体阻力的变化特性。图9为开启空气润滑减阻系统喷气后滚装船的阻力及其分量变化时历曲线。

可以看出，释放空气后（即10 s后），摩擦阻力下降明显，而后维持于一定值；剩余阻力和总阻力的变化趋势相似，释放空气后逐渐降低，最后小幅度波动。本文所设置的通气速率为3.92×10⁻⁴ kg/s，通气速率相对较小，仍能看出相较于摩擦阻力的显著降低，剩余阻力的减小幅度较小。研究表明，空气润滑减阻主要是通过降低船舶航行过程中的摩擦阻力来实现减阻效果。

本文计算总阻力、摩擦阻力和剩余阻力数据统计如表3所示。可以看出，在该通气速率下，滚装船获得总减阻率9.53%的空气减阻效果，减阻效果明显。

7 结　语

1）在船体底部释放空气后，在船体表面会形成包含气体和液体的气液混合流，距离喷气口越近，混合流的密度越小。气层的存在将导致液体与船体表面的直接接触减少，船体浸湿面积降低。气液混合流的存在减少了与船体表面接触的浸湿面积，降低了流体局部的湍流粘度。

2）空气减阻技术基本原理为：喷气装置在船底部释放空气后，船体表面形成覆盖气层和气液混合流，流体密度降低、湍流粘度和动力粘度降低、雷诺应力减小，导致船体摩擦阻力减小，实现润滑减阻的目标。

3）船底通气后，船体所受摩擦阻力降低明显，在设计航速、通气速率相对较小的工况下，空气减阻性能达到了9.53%，减阻效果明显。基于上述研究推测，如果增大通气量，降低空气逃逸速率，空气减阻效果应该会更加明显。