0 引言

气层减阻技术可大幅降低船舶阻力，是实现我国航运业节能减排新目标的重要途径。相对于模型试验，目前气泡船的数值研究相对较少，以下学者开展了有意义的研究：Jin-Keun Choi等^[1]联合采用面元法与求解UN-RANS（Unsteady Reynolds Averaged Navier-Stokes equation）的方法分别计算了喷气下瘦长型船的兴波、船底气层形态、粘性阻力；D.Kim等^[2]采用直接数值模拟（DNS）的方法计算了断阶后气层的非定常流动形态，并与试验结果进行对比，数值结果较为精确的与试验结果吻合；Hoang Cong Liem等^[3]结合边界层积分方程及经验公式，计算了喷气对船舶阻力的影响，该方法可用于计算实尺度气泡船的减阻效果；李云波等^[4-6]基于势流理论，采用边界元方法计算了船底凹槽内的气层形态；王家楣等^[7]基于PHOENICS软件中的IPSA（相间滑移模型）计算了二维气泡船的喷气减阻效果，计算结果与试验结果在趋势上一致；蔡红玲等^[8-10]开展了高速气泡船流场数值模拟，在Fluent平台上采用气液两相流模型进行计算，所得结果在规律上合理，但是由于在船体表明仅采用非结构网格，因而不能较好地模拟边界层流动。

上述研究均为湍流模型、网格、壁面函数等进行系列探讨，所得结果尚需提高。本论文针对气液混合流，系统研究湍流模型、网格、壁面函数对计算结果的影响，最终形成可较好模拟气液混合流的数值方案，对肥大型气泡船气层形态及减阻率预报具有指导意义。

1 Mixture两相流模型

Mixture模型的控制方程由连续性方程和动量方程构成，其具体表达形式如下：

连续性方程为：

$ \frac{\partial }{{\partial t}}({\rho _m}) + \nabla \cdot ({\rho _m}{\vec v_m}) = 0 \text{，} $

(1)

其中$ {\vec v_m} $为质量平均速度，即

$ {\vec v_m} = \frac{{\sum\limits_{k = 1}^n {{\alpha _k}{\rho _k}{{\vec v}_k}} }}{{{\rho _m}}} \text{。} $

(2)

式中： $ {\alpha _k} $为第k相的体积分数；$ {\rho _m} $为混合流密度；$ {\rho _m} = {\alpha _l}{\rho _l} + {\alpha _g}{\rho _g} $；$ {\alpha _l} $为液相体积分数；$ {\alpha _g} $为气相体积分析。

混合相的动量方程由各分相的动量方程相加得到，可表示为：

$ \begin{aligned} \frac{\partial }{{\partial t}}({\rho _m}{{\vec v}_m}) \!+\! \nabla \! \cdot \!({\rho _m}{{\vec v}_m}{{\vec v}_m}) \!=\! -\nabla p \!+\! \nabla \cdot [{\mu _m}(\nabla {{\vec v}_m} \!+\! \nabla \vec v_m^T)] + \\ \quad \quad \quad \quad \quad {\rho _m}\vec g + \vec F + \nabla \cdot (\sum\limits_{k = 1}^n {{\alpha _k}{\rho _k}} {{\vec v}_{dr,k}}{{\vec v}_{dr,k}}) \text{。} \\ \end{aligned} $

(3)

式中：n为混合流中相的数量；$ \vec F $为体积力；$ {\mu _m} $为混合相的粘度系数；$ {\vec v_{dr,k}} $为第k相的漂移速度，即

$ {\mu _m} = \sum\limits_{k = 1}^n {{\alpha _k}{\mu _k}} \text{，} $

(4)

$ {v_{dr,k}} = {\vec v_k}-{\vec v_m} \text{。} $

(5)

式中： $ {\vec v_k} $为第k相速度；$ {\vec v_m} $为各相的平均速度。

该模型同时考虑了相间的滑移速度，即第p相相对于第q相的滑移速度$ {\vec v_{pq}} = {\vec v_p}-{\vec v_q} $，根据Manninen等的研究结果，滑移速度可表示为：

$ {\vec v_{pq}} = \frac{{{\tau _p}}}{{{f_{drag}}}} \cdot \frac{{({\rho _p}-{\rho _m})}}{{{\rho _p}}}\vec a \text{，} $

(6)

$ {\tau _p} = \frac{{{\rho _p}d_p^2}}{{18{\mu _q}}} \text{。} $

(7)

式中： ${\tau _p}$为颗粒松弛时间；d为辅相颗粒（如气泡、水滴等）的直径；$ \vec a $为颗粒加速度。

$ \vec a = \vec g-({\vec v_m} \cdot \nabla ){\vec v_m}-\frac{{\partial {{\vec v}_m}}}{{\partial t}} \text{。} $

(8)

式中：$ {f_{drag}} $为曳力函数。采用Schiller和Naumann的研究结果，可表示为：

$ {f_{drag}} = \left\{ \begin{gathered} 1 + 0.15{{Re} ^{0.687}} \text{，} \quad {Re} \leqslant 1000 \text{，} \\ 0.0183{Re} \text{，} \quad \quad \;\;\;{Re} > 1000 \text{。} \\ \end{gathered} \right.\quad \; $

(9)

2 湍流模型及壁面处理方法

采用光滑平板为对象研究气液混合流数值模型的构建方法。平板的外形及尺寸如图 1所示。平板总长1 200 mm，宽380 mm，厚度10 mm，在距离头部410 mm处安装喷气板。喷气板宽120 mm，长70 mm。

计算流域的网格布局如图 2所示。流域总长为6倍板长，总宽为4倍板宽。流域入口距离平板头部1倍板长，设置为速度入口；出口距离平板尾部4倍板长，为静压力出口；喷气入口设置为质量流量入口。由于气体不可压缩，则根据质量流量可换算得到相应的体积流量。

众多研究表明^[11-15]：只有当气泡分布在湍流边界层内时才具备较好的减阻效果，因此采用Mixture模型计算混合流边界层时，壁面附近网格的处理方法较为关键。

目前，在船舶粘性流场计算中，壁面的处理方法主要有标准壁面函数法、增强壁面函数法及无壁面函数的直接计算方法。

根据湍流模型、壁面函数的自身特征以及数值计算经验，壁面附近网格的处理办法一般基于如下原则：

1）标准壁面函数法。y⁺应大于30～60，最好接近30，湍流边界层内应布置一定的网格数量。标准壁面函数不适用于层流底层，因此第1层网格应布置在对数律层内。

2）增强壁面函数法。该方法可用于计算层流底层的流动，当用于计算层流底层时，y⁺最好取值为1左右，但如果第1层网格布置在了层流层内，略高的y⁺值也可行。计算过程中，层流底层内应布置至少10个网格节点。

3）对于S-A模型（Spalart-Allmaras），若采用增强壁面函数，则y⁺=1；若采用标准壁面函数，则y⁺＞30。

4）对于低雷诺数k-ω模型，网格要求与增强壁面函数相同；采用高雷诺数k-ε模型时，网格要求与标准壁面函数相同。

5）对于LES模型，壁面附近的网格要求与增强壁面函数的要求相同。

本文计算了8种壁面函数及湍流模型的组合方案，如表 1所示。表中y⁺值所对应的来流速度为1.241 m/s。

3 气层形态分析 3.1 气层宏观形态

采用上述湍流模型及壁面函数计算了气流量Q=0.5 m³/h，速度V=1.241 m/s时平板底部的气层形态，所得结果如图 3～图 10所示。图 11为相同气流量及来流速度下气层的模型试验观测结果。

从气层的分布上看：计算方案1与计算方案3所得结果与试验结果在变化规律上相差甚远，说明使用S-A模型与LES模型模拟气液混合流时具有较大局限性，在本文建立的网格方案下，尚需进行优化。计算方案2、方案4～方案8在气层的分布规律上与试验结果较一致，但各方案之间存在差别。

从气泡浓度与气层流态上看：计算方案2、方案3与方案5所得气体在壁面上的分布浓度大于0.95，平板表面上几乎形成了分层流；计算方案7与计算方案8所得气体在壁面上的浓度为0.5～0.8，主要为气体与液体的混合流动。从图 11的试验图片可以看出：气体主要气泡群的形式存在，未形成分层流。说明计算方案7或计算方案8所得结果与试验结果较为接近。

从图 7～图 9（计算方案5～方案7）还可以看出：在湍流模型及壁面函数不变的情况下，壁面法向的第1层网格厚度（y⁺值）对平板表面气泡浓度的计算结果有较大影响。当气流量及来流速度相同时，第1层网格厚度增加，壁面上的气泡浓度降低。

3.2 气层厚度及宽度

采用3.1节的计算方案8计算不同流动工况下平板底部的气层，并将计算所得气层厚度与实验结果进行对比（试验中气层厚度的测量见文献[6]），如图 12～图 14所示。图中横坐标x表示与平板头部的距离，其中x=415 mm处为喷气入口；纵坐标表示沿流动方向不同位置处的气层厚度。

从图 12～图 14可看出，采用Mixture两相流模型与RNG k-ε湍流模型相结合的数值方法，可有效模拟气层厚度沿流动方向的变化规律。

图 15和图 16给出了不同流动工况下气层宽度数值计算结果与实验结果的对比。图中x=0 mm处为喷气入口；纵坐标表示在流动方向上的气层宽度。

从图 15和图 16可看出：采用Mixture模型与RNG k-ε湍流模型相结合的数值方法可有效模拟气层沿平板宽度方向的扩散边界。

4 结语

1）采用RNG k-ε湍流模型、标准壁面函数、第1层网格1 mm、y +为31～35的计算方案可较好地模拟平板底部气液混合层的宏观形态。

2）采用Mixture模型与RNG k-ε湍流模型相结合的方法可有效模拟气层沿平板宽度方向的扩散边界和气层厚度沿流动方向的变化规律，说明搞方法对模拟气液混合流气层形态有效。