平板微气泡减阻数值模拟及影响因素分析

傅慧萍

(上海交通大学海洋工程国家重点实验室，上海200240)

微气泡减阻的隐含机理仍然是研究热点［1］。Legner认为近壁流体可以被认为是空隙率可变的各向同性混合物，宏观量和湍流输运的变化导致摩擦阻力的下降［2］。Madavan等的结论类似，二者的分析都验证了空隙率的存在导致了粘性底层的增厚从而降低了摩擦阻力，并得出“当气泡停留在湍流边界层(TBL)的过渡层时最有效”的结论［3］。Kunz等采用欧拉双流体计算了BDR，得到了与减阻率密切相关的近壁气泡群的演化［4］。多相流在高精度计算方面的进展使气液相互作用的细观研究成为可能。Ferrante等采用欧拉－拉格朗日模型揭示了气泡与湍流在近壁区的相互作用细节［5］。Lu等采用直接数值模拟揭示了:当大的可变形气泡出现在过渡区，会与流向涡产生强烈的相互作用［6］。Mohanarangam等基于ANSYS CFX中的多尺寸组模型对二维平板BDR进行了数值模拟，并全面借鉴了试验模型的试验参数及丰富的试验数据进行了数学模型的校验［7］。丁力等采用混合物模型对二维简化船型进行了模拟，得到了喷气参数对船舶阻力的影响规律，在喷气量的相似换算方面取得了一些成果［8］。陈显文等同样采用混合物模型对某回转体及SUBOFF进行了模拟。由于采用的也是混合物模型，无法考虑微气泡运动过程中的变形、破裂和聚合，也忽略了气泡尺度的影响，而且数值计算结果没有试验数据进行支撑与校验［9］。王炳亮等分别采用Euler双流体模型(均一气泡直径)和混合物模型对三维平板及某散货船模型进行了模拟［10］。与前两者的区别是参照了Madavan等［11］的试验参数进行数值建模。但并没有采用后者的试验结果进行对比，而是采用了一系列的平板阻力经验公式进行阻力评估获得CFD结果的参考数值进行误差分析。我国目前所拥有的微气泡减阻技术离实用还有一段距离。根据现代测试手段获取的试验数据，建立更加符合物理本质的微气泡减阻模型;发展高效、高精度的高雷诺数网格生成技术与数值分析方法，用以分析微气泡减阻的流动特性和机理，是微气泡减阻数值模拟的发展趋势。

1 研究对象

本文的研究对象为平板，为与试验对比，依文献［11］进行建模。图1给出了平板底部通气的二维数值模拟示意图。计算域总长0.61 m，高0.057 m。沿长度方向，将平板分为3段:wall-1、wall-2、wall-3，长度分别为:0.178 m、0.178 m、0.254 m，其中wall-2为喷气入口。

图1 计算域Fig.1 Computational field

图2为三维建模中的底板部分，采用对称边界条件。阴影部分为所考察的平板，宽为0.051 m，计算域宽度为0.250 m。

图2 三维建模Fig.2 3D model

网格划分的边界层最小网格尺度为Δyp按下式计算:

式中:y+为当地雷诺数，L为特征长度，此处取计算域总长。由于本文的研究对象为5～500 μm的微气泡，故边界层网格要求较高，当地雷诺数要求为个位数的量级。计算中采用第1层网格高度是0.001 5 mm，图3所示为边界层网格划分。

图3 边界层网格Fig.3 Mesh in the boundary

边界条件的设置如图1所示，左端和wall-2分别为水流和喷气入口，均设置为速度入口;右端为流动出口，设置为outflow类型;wall-1、wall-3设置为无滑移壁面;计算域最大高度处设置为水流速度入口条件。当计算无通气平板的摩擦阻力时，wall-2设为壁面条件即可。

2 多相流模型

2.1 混合物模型

混合物模型用来模拟2种流体混合之后的两相流，它将各相设置为相互贯穿的连续体，从而只对混合物求解共同的控制方程。连续方程:

式中:Um为混合相速度，ρm为混合相密度。动量方程:

式中:μm为混合相黏度，p为压力，j为多相流的相数(j=1，2)，αj、ρj、Udr，j分别为第j相的体积分数、密度和漂移速度。

气相体积分数方程:

式中:Sg为气相的源项。相对速度(或滑移速度)，是指次要相(第2相)相对于主要相(第1相)的速度:

设任意相(j相)的质量分数定义为

则漂移速度与相对速度通过下式关联:

FLUENT里的混合物模型采用代数滑移公式，基本假设是为了能对相对速度进行代数表达，要求在一个较短的空间尺度上达到局部相间平衡。本文采用Manninen相对速度形式。

2.2 相群平衡模型

相群平衡模型是在欧拉双流体框架下开启，它将第二相(空气)进行尺寸分组，对每一组的体积分数进行求解。假设将气泡分成N个尺寸组，各组气泡之间存在聚并和破碎作用。那么第i组气泡的连续性方程为

式中:Si为第i组气泡由于聚并和破碎产生的源项，本文采用Luo模型;fi为第i组的体积分率:

假设ni和vi分别为第i组气泡的数密度和气泡体积，则有如下关系:

假设气泡是球形的，各组气泡体积比为

若拟在5～500 μm内对气泡直径进行分组(N =5)，则q=4.985。可以得到各尺寸分组如下:D0= 500 μm、D1=158 μm、D2=50 μm、D3=16 μm、D4= 5 μm。索特平均直径Dm是联系流场与微观气泡群的一个变量，即欧拉双流体模型和相群平衡模型通过它联系起来。该变量与各气泡分组的直径及其数密度有关，定义如下:

至于湍流的模拟，本文3种多相流求解方法均是对混合相进行湍流方程的求解，所采用的湍流模型为标准k-ε两方程模型。

3 计算结果与分析

3.1 重力的影响

为考察重力的影响，计算时先不考虑重力，待计算收敛后，再考虑重力至再度收敛。图4给出了相同水流速度和喷气速度下(Uw=9.6 m/s，Ua= 0.165 2 m/s)，气泡直径D分别采用5、50和500 μm得到的平板阻力系数收敛历程。图中曲线基本呈2段式:前一段不考虑重力，后一段考虑重力后，阻力明显呈上升趋势。这是由于考虑重力后，气泡在水中受到向上的浮力作用。对于底部通气方式，通气方向与浮力方向一致，浮力作用使气泡离开壁面，气体的润滑作用减小，阻力增加。从前后2段曲线的提升程度看，小气泡的阻力增加幅度小于大气泡，大气泡阻力系数大。

图4 重力的影响Fig.4 Effect of gravity on Cd

3.2 通气方式比较

由上节的研究可知重力对大气泡的影响显著，为此选择 500 μm泡径，在 Uw=9.6 m/s，Ua= 0.165 2 m/s时，将二维平板模型旋转180°，使底部通气变为顶部通气，探讨通气方式对BDR的影响。计算时仍然先不考虑重力，待第1阶段收敛后再考虑。图5给出了2种通气方式的比较。由图5可见:重力对2种通气方式的影响是完全相反的。对于底部通气，如前节所述，重力使阻力增大。对于顶部通气，通气方向与浮力方向相反，浮力是使气泡趋向壁面并停留在边界层内的，气体的润滑作用得以保持所以使阻力减小。对于船舶微气泡减阻，应该是适用于顶部通气方式，重力的影响是有利于减阻的。

图5 通气方式的影响Fig.5 Effect of injection method on Cd

3.3 减阻机理研究

3.3.1 气层厚度与减阻率

设Cd为通气情形下平板“wall-3”段的无量纲摩擦阻力系数;Cd0为不通气情形下的“wall-3”段无量纲摩阻系数;则减阻率DR由下式定义:

图6给出了Uw=9.6 m/s时，均一气泡直径D= 50 μm时的采用混合物模型得到的减阻率与气层厚度的关系。减阻率先是随喷气速度增大而增大，当喷气速度达到一定数值后，减阻率不增反降，即减阻率曲线存在一个峰值。气层厚度却不同，随喷气速度增加而单调增大。因此减阻率存在一个峰值，对应着一个最佳气层厚度或最佳喷气速度。这是因为喷气速度过小时，微气泡还不足以形成足够的气层，甚至类似于表面粗糙度，产生负增长;如果喷气速度过大，减阻效果也会下降，这是因为垂直向上的气流对水平方向的来流造成的扰动影响了流动的光滑性，导致了减阻效果的下降。

图6 减阻率与气层厚度的关系Fig.6 Relation between DR and gas layer thickness

3.3.2 体积分率以及气泡数密度

考虑到采用均一气泡直径的混合物模型得到的减阻率与试验值对比偏高，而PBM模型可以对气泡进行尺寸分组，更好地模拟了利用金属烧结板形成微气泡的气泡生成方式，本节将采用PBM模型对平板底部通气BDR进行数值模拟。利用PBM模型进行计算时，边界条件中通气口wall-2的各分组的体积分率均设为0.2，即假设通气口处，各种直径的气泡对于气相分数的贡献一致。气泡尺寸分组是PBM模型特有的，可以通过对体积分率以及气泡数密度的分析来观察通气口下游平板wall-3处各组气泡的分布及演化。本节给出了二维平板底部通气当Uw=9.6 m/s，Ua=0.055 1 m/s时，采用5～500 μm的尺寸分组(分10组)的PBM模拟结果。各组气泡的直径和平板wall-3处的体积分率如表1所示，对应的曲线见图7。由表1及图7可以看出:在平板wall-3处，直径最大和最小的气泡都有明显减少，而中间尺寸的气泡体积分率有所增大，各组分的体积分率不再相等。

表1 各组气泡的直径和wall-3处体积分率Table 1 fion wall-3 to bubble size Di

图7 wall-3处各组气泡对应的体积分率Fig.7 fion wall-3 to bubble size Di

去流向各个截面上不同尺度的气泡数密度及体积分率变化如图8所示。由图8可见，大部分尺度的气泡数密度和体积分率沿流向都呈减小趋势，只有大气泡组呈增大趋势，两者的变化趋势是一致的。即小气泡在流动方向上发生聚合变成大气泡，大气泡发生破裂，但由于聚合的影响大于破裂，总体上大尺度分组的数密度和体积分率增大。

图8 气泡数密度与体积分率沿流向的变化Fig.8 niand fialong x axis

3.4 三维效应

3.3 节采用PBM模型对二维平板底部通气BDR进行模拟，得到的减阻率与试验值对比仍然偏高。为此，本节开展基于PBM的三维数值模拟。图9给出了相同水流速度和喷气速度下(Uw=9.6 m/ s，Ua=0.055 1 m/s)，采用5～500 μm的尺寸分组(分10组)进行计算，得到的阻力系数收敛历程。由图9可见，三维较之于二维，阻力系数提升一个量级。多相流中α表示第二相即气相的体积分数:α= 1，表示纯气流;α=0，表示纯水流;α为0～1的任意值时，表示气水混合流。本节定义的气层厚度δα表示气相分数为α的等值面上的最大厚度。

图9 三维效应对阻力系数的影响Fig.9 Effect of 3-D on Cd

图10给出了当α分别为0.1、0.25、0.5时的等值面厚度分布。由图10可见，由于气泡的上浮运动，α等值面沿流向厚度逐渐增加，最大厚度出现在出口靠近中心线处;随着α的增大，等值面最大厚度减小，外缘沿流向向中心线收缩。

图10 三维BDR模拟Fig.10 3-D simulation of BDR

3.5 侧壁通气

侧壁通气BDR模拟必须是三维全流场模拟，由此也可进一步完善重力影响的研究。本节对平板侧壁通气BDR进行基于PBM的三维数值模拟。图11给出了相同水流速度和喷气速度下(Uw=9.6 m/s，Ua=0.055 1 m/s)，采用5～500 μm的尺寸分组(分10组)进行计算，得到的侧壁通气与底部通气的Cd收敛历程对比。由图11可见，侧壁通气经过了一个伪收敛平台再下降至最后的收敛值。可能因为松弛因子太低，收敛很慢。

图11 侧壁通气的阻力系数Fig.11 Cdwith plate vertical

图12给出了当α分别等于0.1、0.25、0.5时的等值面厚度分布。由图12可见，侧壁通气的气层厚度要大于底部通气的。在α=0.1的等值面上最大厚度为9.83 mm，为底部通气(3.86 mm)的2.55倍。但此最大厚度已偏离出所考察的平板(wall-3)。

另外，侧壁通气的减阻效果低于底部通气，这是因为:由于气泡的上浮运动导致大部分区域气相分数降低，甚至局部区域出现全沾湿状态，即未被气体润滑。

图12 侧壁通气BDR模拟Fig.12 Simulation of BDR with plate vertical

3.6 计算方法验证

基于PBM方法，采用50～500 μm的5分组气泡群，对Uw=9.6 m/s时的平板底部喷气进行数值模拟。表2和图13给出了数值结果与试验值［3］的比较。

表2 计算值与试验值的比较Table 2 Comparison between CFD and EFD

图13 数值模拟校验Fig.13 Validation of numerical method

数值模拟所得变化趋势与实验所得变化趋势基本相同，而计算所得与试验值相比偏高，这个误差并不影响本文对微气泡减阻的机理及尺度效应的定性研究。这可能是由于二维模拟所带来的系统误差所致。如果将二维模型变为三维模型，并进一步增大气泡直径尺寸分组，可以预见，阻力系数的预报精度将提高，但随之而来的是计算量增大和计算稳定性问题。

4 结论

1)重力对大气泡影响较大:考虑重力的影响后，气泡就受到了水中的浮力，浮力使气泡上升。越大的气泡所受浮力也越大。

2)重力对通气方式的影响:对于顶部通气方式，重力使气泡停留在壁面边界层，从而使减阻率相较于底部通气方式得到改善。

3)减阻率与气层厚度之间存在一定的关系:喷气速度增加，气层厚度也增加，但减阻率在达到一个峰值后开始下降。

4)三维效应使减阻率更接近试验值:故三维建模与计算是必要的。

5)侧壁通气的减阻效果不佳:重力效应使流场上下不对称，气泡沿流向向上漂移，考察对象局部区域全沾湿导致。

进一步的研究将拓展到三维船体底部通气建模与计算。

［1］CECCIO S L.Friction drag reduction of external flows with bubble and gas injection［J］.Annual Review of Fluid Mechanics，2010，42(1):183-203.

［2］LEGNER H H.A simple model for gas bubble drag reduction［J］.Physics of Fluids，1984，27(12):2788-2790.

［3］MADAVAN N K，MERKLE C L，DEUTSCH S.Numerical investigations into the mechanisms of micro-bubble drag reduction［J］.Journal of Fluids Engineering，1985，107(3): 370-377.

［4］KUNZ R F，GIBELING H J，MAXEY M R，et al.Validation of two-fluid Eulerian CFD modeling for microbubble drag reduction across a wide range of Reynolds numbers［J］.Journal of Fluids Engineering，2007，129(1):66-79.

［5］FERRANTE A，ELGHOBASHI S.Reynolds number effect on drag reduction in a microbubble-laden spatially developing turbulent boundary layer［J］.Journal of Fluid Mechanics，2005，543:93-106.

［6］LU Jiacai，FERNáNDEZ A，TRYGGVASON G.The effect of bubbles on the wall drag of a turbulent channel flow［J］.Physics of Fluids，2005，17(9):095102.

［7］MOHANARANGAM K，CHEUNG S C P，TU J Y，et al.Numerical simulation of micro-bubble drag reduction using population balance model［J］.Ocean Engineering，2009，36(11):863-872.

［8］丁力.微气泡减阻喷气参数换算关系研究［D］.武汉:武汉理工大学，2012:1-47.

DING Li.The conversion relation research about jet parameters of micro-bubbles drag reduction［D］.Wuhan:Wuhan University of Technology，2012:1-47.

［9］陈显文.回转体微气泡减阻和噪声的数值研究［D］.武汉:华中科技大学，2012:1-77.

CHEN Xianwen.Numerical simulation of axisymmetric body's drag reduction and noise by micro bubbles［D］.Wuhan: Huazhong University of Science and Technology，2012:1-77.

［10］王炳亮.船舶微气泡减阻数值模拟及机理研究［D］.哈尔滨:哈尔滨工程大学，2012:1-70.

WANG Bingliang.Numerical simulation and mechanism research on drag reduction of ship by microbubbles［D］.Harbin:Harbin Engineering University，2012:1-70.

［11］MADAVAN N K，DEUTSCH S，MERKLE C L.Reduction of turbulent skin friction by microbubbles［J］.Physics of Fluids，1984，27(2):356-363.