低速船舶微气泡减阻数值研究

赵晓杰，宗 智，王加夏，洪智超，4，胡俊明

（1.江苏科技大学船舶与海洋工程学院，江苏镇江 212000；2.大连理工大学船舶工程学院，辽宁大连 116024；3.工业装备结构分析国家重点实验室，辽宁大连 116024；4.江苏省船舶与海洋工程装备技术创新中心，江苏南通 226100）

0 引 言

航行中的船舶会受到阻力。为了实现减阻，研究者们提出了许多减阻技术[1]，如为了减少兴波阻力，采用球鼻艏或优化船体型线；为了降低摩擦阻力，在船舶表面设置超疏水涂层，或在船底喷入高分子聚合物、微气泡等。其中，微气泡减阻技术受到人们的广泛关注，因为它的减阻效率较高并且无污染。对于低速船舶，摩擦阻力占主要部分，微气泡减阻技术更具应用前景。

几十年来，国内外学者对微气泡减阻进行了许多试验研究，研究对象主要有平板、回转体、船体等。McCormick 等[2]首先通过试验发现了微气泡减阻现象，然后研究学者又相继对平板和回转体的微气泡减阻进行了试验研究[3]。对于船舶微气泡减阻，减阻情况更加复杂，因为减阻效果会受船型、喷气形式、复杂工况等多个因素影响。Sayyaadi等[4]对一双体船进行了试验研究，发现喷气量对减阻效果有显著影响，并存在一个临界喷气量；王家楣等[5]采用多孔硅板生成微气泡，对不同吃水的船模微气泡减阻进行了试验研究，发现船舶在浅吃水下减阻更好。

前人也对微气泡减阻做了许多数值研究。Kunz 等[6]利用两相欧拉模型研究了平板的微气泡减阻，通过和试验结果对比验证了模型的准确性；傅慧萍[7]利用不同多相流模型对平板微气泡减阻进行了数值模拟并对影响因素进行了分析；对于船舶微气泡减阻，郑晓伟等[8]利用相间滑移算法研究了二维船的微气泡减阻，所得结果与试验定性符合；吴浩与欧勇鹏[9]将船底简化为平板并利用Mixture模型研究喷气减阻，较好地模拟出了气液混合流。可以看出，目前船舶微气泡减阻模拟都对数值模型进行了一定简化，几何模型简单，并且没有充分考虑气液间及气泡间相互作用。

因此，为了进一步研究船舶微气泡减阻规律，本文基于OpenFOAM 中两相欧拉模型对低速散货船的微气泡减阻展开数值研究，考虑五种相间作用力以及气泡的聚合和破碎，并考虑气泡对液体湍流的影响。区别于二维平板，船舶具有不同阻力成分，本文研究气泡对船舶摩擦阻力、粘压阻力和总阻力的减阻率影响，同时研究喷气量、气泡直径、航速及吃水等因素对减阻率的影响，分析船舶附近气体体积分数、湍流粘度和气泡直径的分布等，为微气泡减阻技术在船舶上的应用提供理论指导。

1 数值模型

1.1 两相欧拉模型

在OpenFOAM 的两相欧拉模型中，对气液两相（本文指水和空气）分别进行建模，建立的质量和动量守恒方程[10]分别为

式中：i为相指数，这里i=l或g分别表示液体和气体；αi表示i相体积分数；Ui表示i相流体速度；ρi表示i相流体密度；p为压力；g为重力加速度；∑Fi为i相的相间作用力，用于考虑相间动量传递；νi表示i相运动粘度系数；νt,i表示i相湍流运动粘度系数。

1.2 相间作用力

为了充分考虑气液间动量传递作用，本文考虑五种相间作用力，公式如下：

式中：Fd为相间阻力，Fl为升力，Fvm为虚质量力，Ftd为湍流离散力，Fwl为壁面润滑力。

相间阻力是气泡在液体中运动时受到的阻力，与相对速度方向相反，计算公式为

式中：αg为气体体积分数，dg为气泡直径，ρl为液体密度，Ul和Ug分别表示液体速度和气体速度，Cd为相间阻力系数。本文采用Ishii和Zuber[11]提出的阻力模型，该模型考虑三种气泡形态，对于球形气泡，阻力系数计算公式为

对于椭球型气泡，阻力系数计算公式为

式中，Re是气泡雷诺数，Eo表示无量纲Eötvös数，f(αg)表达式为

对于帽子形气泡，阻力系数计算公式为

不同气泡形状区域划分基于以下公式：

升力是由于液体速度梯度引起的剪切力，它是除阻力外最重要的作用力，可以影响气体垂向分布，计算公式为

式中：Cl是升力系数。采用Tomiyama升力模型[12]，可以考虑不同直径的气泡，计算公式为

式中，Cvm为虚拟质量力系数，这里假设气泡是球形，根据势流理论取Cvm=0.5。

湍流耗散力为流体湍流脉动作用对气泡产生的作用力，采用Bertodano模型[13]，计算公式为

壁面润滑力是为了防止气泡接触壁面而产生的使气泡远离壁面的作用力，计算公式[14]为

式中，壁面润滑力系数Cwl的计算公式为

式中：yw为气泡离壁面的垂直距离；Cw1和Cw2是校准系数，根据文献[14]推荐，取值分别为Cw1=-0.01 和Cw2=0.05，表示壁面润滑力在距离壁面5倍直径距离时变为0。

1.3 气泡群体平衡模型

为了考虑气泡直径变化，采用基于界面面积密度方程的气泡群体平衡模型，控制方程[15]为

式中：a表示单位体积气泡界面面积密度；ψ表示气泡形状因子，忽略气泡变形，ψ取值为1/36π；Sj表示气泡聚并和破碎的源项，可以分成三部分：

其中，SRC和SWE分别表示随机碰撞和尾迹夹带作用，它们会引起气泡聚并，表达式分别为

式中：ng为单位体积气泡数密度，ut与ur分别为湍流速度和气液相对速度。根据Ishii等[15]实验研究，对于随机碰撞，常数取值CRC=0.04和C=3。对于尾迹夹带，常数取值CWE=0.04和αmax=0.75。STI为湍流作用，它导致气泡发生破碎，表达式为

式中，CTI为湍流作用参数，Wecr为临界韦伯数，取值分别为CTI=0.06和Wecr=6.0。

1.4 湍流模型

对于湍流模型，本文采用Lahey[16]提出的改进k-ε模型。为了考虑气泡影响，在控制方程中加入源项，并增加额外湍流粘度来计算气泡导致的湍流增加。湍流模型的控制方程为

式中，湍流粘度系数νt的计算公式为

源项Φkl与Φεl的计算公式为

式中，常数Cp取值为0.25，其它参数为Cμ=0.09，Cε,1=1.44，Cε,2=1.92，σk=1.0，σε=1.3。

2 计算设置

2.1 几何模型

本文研究的是95 000 DWT 散货船，缩尺比为1:50，实船垂线间长为231 m,型宽为38 m，型深为20.7 m，设计吃水为12.5 m，排水体积为91 893.0 m3，湿表面积为12 543.0 m2，设计航速为14.5 kn，船体型线见图1。在船模17站处设置长方形喷气板，长400 mm，宽90 mm，如图2所示。

图1 散货船横剖线图Fig.1 Sectional curves of the bulk carrier

图2 船舶喷气板示意图Fig.2 Sketch of the jet plate on the ship model

2.2 计算域和网格划分

计算区域如图3 所示，不考虑自由面影响，采用叠模模型，将自由面当作对称滑移壁面。为了减少计算量，采用一半船体进行计算。速度入口距离船体1倍船长，出口边界距离船体3倍船长，底面和侧面距离船体1倍船长。

图3 计算区域Fig.3 Computational domain

网格划分利用切割体网格，网格分布如图4所示。最大网格尺寸为231 mm，在船体附近进行网格加密，最小网格尺寸约为7.2 mm，船体表面设置三层边界层，第一层网格厚度满足y+>30。

图4 网格分布Fig.4 Mesh distribution

2.3 边界条件和计算方法

对于边界条件，船体对液体速度采用无滑移壁面条件，对于气泡速度采用滑移边界。入口为速度入口，出口为压力出口。对于喷气面，本文假设采用微米多孔板产生微气泡，为了匹配实际物理条件，假定孔隙率为0.3，与Kunz等[6]数值研究中条件一致。对于其他边界，设置为对称边界条件。

控制方程采用有限体积法进行离散。对于空间离散，梯度项和拉普拉斯项采用高斯线性格式（Gauss linear），对流项采用限制高斯线性格式（Gauss limited linear），在相输运控制方程中，采用Gauss vanLeer离散格式。空间离散为二阶精度。对于时间离散，采用一阶精度的隐式欧拉格式。

控制方程求解采用隐式速度压力耦合算法（PISO），在OpenFOAM[17]中利用PIMPLE 循环进行设置，外迭代次数设为1，压力矫正次数为2。时间步采用自适应时间步，最大库朗数设为0.9，从而保证计算收敛。

2.4 计算工况设置

为了研究航速、喷气量、气泡直径和吃水等因素对船舶微气泡减阻的影响，设置了不同计算工况。考虑三种航速，分别为0.946 m/s、1.055 m/s 和1.164 m/s，其中1.055 m/s 为设计航速。考虑三种吃水，分别为0.138 m、0.25 m 和0.29 m，其中0.25 m 为模型船的设计吃水。考虑三种气泡直径，分别为200 μm、500 μm和800 μm。计算工况设置总结见表1。

表1 计算工况总结Tab.1 Summary of calculation conditions

采用无量纲喷气量Cv作为参数，表达式为Cv=Qg/(Qg+Ql)，式中Qg是气体喷气流量，Ql是边界层内液体流量。Ql计算公式为Ql=U( )δ-δ∗b，式中U为船舶航速，b为船宽，δ为边界层厚度，δ∗为边界层排挤厚度，表达式为δ∗=0.125δ，根据边界层理论[18]，δ计算公式为δ=0.37x/Re0.2。无量纲喷气量Cv范围为0～0.05，间隔为0.01。

本文研究三种阻力成分的减阻率，计算公式为ΔRi=(Ri-R′i)/Ri，式中Ri表示未喷气时阻力，R′i表示喷气后阻力，下标i取f、p和t时分别表示摩擦阻力、粘压阻力和总阻力。

3 计算结果与分析

3.1 计算模型验证

为了验证本文模型，采用两种算例：一是计算不喷气下船舶摩擦阻力，并和ITTC 推荐的摩擦阻力计算公式进行对比；二是计算喷气情况下平板微气泡减阻，并和前人试验数据进行对比。其中ITTC推荐的摩擦阻力系数公式为

表2 表示船舶在三种速度下摩擦阻力和通过ITTC 推荐公式计算得到的摩擦阻力结果对比，可以看出计算误差在2%之内，说明本文的数值模型可以准确计算船舶阻力。

表2 船舶计算摩擦阻力和ITTC公式计算摩擦阻力对比Tab.2 Comparison of calculated ship resistance and the values from ITTC formula

平板微气泡减阻试验验证采用Merkle 等[3]的实验，通过多孔介质板产生微气泡实现减阻，采用实验工况，考虑流速为9.6 m/s和14.2 m/s，本文气泡群体平衡模型需要先给定平均气泡直径，根据实验观测[3]以及Mohanarangam 等[19]在相同工况计算所得气泡直径分布，两种流速下平均气泡直径分别取600 μm和400 μm。图5和图6分别表示流速为9.6 m/s和14.2 m/s时的数值结果，横坐标为喷气量，纵坐标为平板摩擦阻力减阻率，可以看出数值结果和试验值吻合良好，说明本文微气泡减阻数值模型是可行的和准确的。

图5 速度为9.6 m/s时数值结果Fig.5 Numerical results at a speed of 9.6 m/s

图6 速度为14.2 m/s时数值结果Fig.6 Numerical results at a speed of 14.2 m/s

3.2 不同气泡直径下船舶微气泡减阻率

气泡直径是影响微气泡减阻率的重要因素。图7（a）～（c）表示在不同气泡直径下船舶摩擦阻力、粘压阻力和总阻力的微气泡减阻率，对应计算工况1～3。从图中可以看出喷气情况下，摩擦阻力、粘压阻力和总阻力都实现了减阻，并且随着喷气量的增大，减阻率也都随之增大，说明喷气量是影响船舶微气泡减阻率的决定因素。另外，摩擦阻力和总阻力的减阻率增大速度一开始很快，但随后会变得缓慢，直到达到最大值后几乎不再改变，说明对于摩擦阻力和总阻力，存在一个最佳临界喷气量，这和Sayyaadi等[4]试验研究所得的结论相一致。对于粘压阻力，减阻率增加存在波动，这是因为微气泡对流场存在扰动、局部压力出现不稳定波动造成的，但总体上减阻率单调增加。

图7 不同气泡直径下不同阻力成分微气泡减阻率Fig.7 Microbubble drag reduction rates of different drag components for different bubble sizes

对于不同气泡直径，从图7（a）～（c）可以看出，较小的气泡可以获得较大的摩擦阻力、粘压阻力和总阻力减阻率，说明小气泡可以更有效地减阻，这和Song等[20]三维回转体的微气泡减阻试验结论相一致。由于船舶是低速散货船，摩擦阻力占比较大（约80%），而粘压阻力占比较小（约20%），尽管粘压阻力减阻率存在不稳定波动，但是摩擦阻力减阻率和总阻力减阻率的变化趋势还是彼此接近。

3.3 不同气泡直径下船舶附近气体分布

图8（a）～（c）给出了无量纲喷气量为0.05时不同气泡直径下船舶附近气体体积分数分布。可以看出在最小气泡（d=200 μm）时，船底气体分布均匀，计算得到表面气体体积分数均值为0.327，此时减阻率也较大；在最大气泡（d=800 μm）时，靠近喷气板处气体覆盖较多，随着远离喷气板气体逐渐减少，导致船底气体分布不均匀，计算得到表面气体体积分数均值为0.308，此时减阻率也较小，说明近壁面气体体积分数分布是影响船舶微气泡减阻的关键因素。

图8 不同气泡直径下船舶附近气体体积分数分布Fig.8 Distribution of air volume fraction around the ship for different bubble sizes

不同气泡直径导致船体附近气体体积分数分布不同，原因是不同气泡直径的气液相互作用力不同。由公式（4）可以看出，在相同气体体积分数时，气泡直径越小，气液间相间阻力越大，使得小气泡从喷气板喷出后难以远离壁面，更容易沿着液体流动方向延伸至船尾，从而造成气体分布均匀，实现较高减阻率；而对于大气泡，气液相间阻力较小，使得大气泡难以沿着液体流动方向向船尾延伸，从而导致气体在喷气口附近聚集，船尾附近气体体积分数较低，导致减阻率也较低。

3.4 不同气泡直径下湍流运动粘度分布

为了进一步分析不同气泡直径下船舶微气泡减阻率不同的原因，图9（a）～（d）给出了不喷气时及喷气时（Cv=0.05）不同气泡直径的船舶附近湍流运动粘度分布。可以看出，喷气后船底液体湍流运动粘度都会变大，说明气泡对液体湍流具有增强作用，这和Mohanarangam 等研究得到的结论相一致。根据公式（23），可以得到在相同气体体积分数下，气泡直径越大湍流运动粘度也越大的结论。图9也表明，800 μm、500 μm 和200 μm 直径气泡产生的湍流运动粘度依次减小，因此这也是小气泡可以更有效减阻的原因之一。

图9 不同气泡直径下船舶附近湍流运动粘度分布Fig.9 Distribution of turbulent kinematic viscosity around the ship for different bubble sizes

3.5 船舶附近气泡直径变化

图10（a）～（c）为无量纲喷气量为0.05 时不同初始气泡直径下船舶附近气泡直径分布。从图中可以看出，在船底沿着流动方向气泡直径都有所增大，这和Merkle 等[3]平板微气泡减阻实验的气泡变化趋势相同，原因是气泡沿流向会发生随机碰撞，再加上相邻气泡尾迹夹带的作用，气泡会发生聚并，导致气泡直径变大。在船舶尾部附近，由于存在较大的几何形状突变，船尾湍流作用剧烈，就会容易导致气泡破碎，气泡直径又会变小。

图10 Cv=0.05时不同初始气泡直径下的船舶附近气泡直径分布Fig.10 Distribution of bubble size around the ship for different initial bubble sizes when Cv=0.05

另外，从图中可以看出不同初始气泡直径下的气泡直径变化情况不同。对于较小初始直径气泡如200 μm直径气泡，在喷气板不远处就开始发生聚并，最大气泡直径可达约260 μm，气泡直径增加约30%，大气泡在船底分布也多；而对于较大初始直径气泡如800 μm 直径气泡，在距离喷气板较远处才开始发生聚并，最大气泡直径约为910 μm，气泡直径增加约13%，增大的气泡仅分布在船尾附近。可见较小气泡的聚并作用更强，原因在于不同气泡直径导致的气泡聚并和破碎作用不同，根据公式（18）～（19），气泡随机碰撞和尾迹夹带作用与气泡数密度有关，相同喷气量下小气泡的气泡数密度更大，所以气泡聚并作用更强。又根据3.4 节可得，小气泡的湍流增强作用较小，所以湍流作用导致的气泡破碎作用也较弱，较小气泡的聚并作用更加剧烈。

3.6 航速与吃水对船舶微气泡减阻的影响

图11 为不同船舶航速下不同阻力成分的微气泡减阻率，对应工况1、4 和5。可以看出，不同航速下摩擦阻力和总阻力减阻率都随喷气量的增加而增加，到最后增加速度会变慢，减阻率都会接近一个恒定值；粘压阻力减阻率也随喷气量的增加而增加，但也存在波动，这是因为气泡导致流场局部压力不稳定造成的。此外，航速较高时船舶减阻率较高，说明此时船舶附近气体体积分数较大。这是因为不同航速下气液相互作用力不同，根据公式（4）可知流速较高时，气泡和水的相对速度较大，此时气液相间阻力较大，气泡喷出后远离壁面趋势较小，气泡在较强的气液相互作用下更容易延伸至船尾，造成船底气体体积分数较大且气体分布均匀，因此船底获得更大气体覆盖面积，从而获得较大减阻率。

图11 不同航速下不同阻力成分微气泡减阻率Fig.11 Microbubble drag reduction rates of different drag components for different ship speeds

图12 为不同船舶吃水下不同阻力成分的微气泡减阻率，对应工况1、6 和7。可以看出，对于摩擦阻力、粘压阻力和总阻力，小喷气量时减阻率差别较小，随着喷气量增加，减阻率差别逐渐变大。另外，小吃水时减阻率较大，这和王家楣等[5]试验的结论相一致，这是因为船舶减阻主要原因是船底摩擦阻力的减少，小吃水时船底摩擦阻力在整个船体总阻力中占比较大，更容易获得较大的减阻率。

图12 不同船舶吃水下不同阻力成分微气泡减阻率Fig.12 Microbubble drag reduction rates of different drag components for different ship drafts

4 结 论

本文基于OpenFOAM 中两相欧拉模型对三维低速散货船的微气泡减阻进行了数值研究，考虑五种相间作用力、气泡聚合和破碎以及气泡对液体湍流的影响，研究了气泡对不同阻力成分减阻率的影响，研究了喷气量、气泡直径、航速及吃水等因素的影响，分析了船舶附近气体体积分数、湍流粘度和气泡直径的分布，得出的主要结论如下：

（1）气泡可以同时减少低速船舶摩擦阻力、粘压阻力和总阻力，并且随着喷气量增大，减阻率也都随之增大，说明喷气量是影响船舶微气泡减阻率的决定因素。

（2）船舶附近气体体积分数分布是影响船舶微气泡减阻的关键因素。气泡较小时船舶附近的平均气体体积分数较大，且气体分布更加均匀，减阻率较高。另外，较小气泡产生的湍流运动粘度也较小，这也是小气泡可以更有效减阻的原因之一，因此在工程应用中应尽量保证气泡直径足够小。

（3）气泡会沿着液体流动方向发生聚并，导致气泡直径增大。较小气泡的聚并作用更剧烈，原因是相同喷气量下它的气泡数密度更大，随机碰撞和尾迹夹带作用更强，另外较小气泡的湍流增强作用较小，导致气泡破碎作用也较弱。

（4）较高航速和小吃水更有利于船舶微气泡减阻。较高航速时，气液相间阻力较大，船底气体体积分数较大且气体分布均匀，从而获得较大减阻率。小吃水时，船底摩擦阻力在整个船体总阻力中占比较大，更容易获得较大的减阻。