0 引　言

通气减阻研究由来已久，这一现象首先由McCormick和Bhattacharyya在实验中发现，随后由Madavan、Deutsch和Merkle通过一系列实验证实^[1]。2021年，Sindagi等^[2]对散货船在不同航速（6～10 kn）下的缩比模型进行了实验研究，实验采用了2 mm和1 mm尺寸的喷射孔和不同的喷气速率，结果表明更大的孔可以产生更高减阻率。在数值分析方面，徐庚辉等^[3]基于CFD软件Fluent对两相流作用下的集装箱船船模开展微气泡数值模拟研究，分析了喷射两相流中的含气率、喷气孔位置、喷气速率对减阻性能的影响，发现流体中微气泡含量过高时会增大船体的摩擦阻力，流体含气率为10%～20%时减阻效果最好；赵晓杰等^{[4 − 5]}采用OpenFOAM的两相欧拉求解器对轴对称结构、二维平板的两相微气泡流动进行数值模拟研究，采用标准k-ε湍流模型，并通过求解界面输运方程来模拟气泡的聚并和破碎，结果表明边界层内的气流速度和空气体积分数对减阻效果起到重要作用，通气速度较大且气泡直径较小时减阻效率高。宋武超等^[6]则表明回转体的外形使其微气泡减阻特性受微气泡流形态的影响更大，与平板减阻特性有很大不同。并且在微气泡减阻有关实验研究中，关于微气泡大小和空气注入流量对减阻率的影响，公开发表文献中仍然存在矛盾^[7]。

综上，对于通气减阻的机理，数值模拟方法及数学模型本身，还存在着许多尚未揭示的尚无法准确建模的隐含因素。本文将多孔介质模型应用于通气减阻的模拟，试图考虑通气装置对减阻率的影响，同时考虑重力和自由面的影响，以期形成一种更加符合物理规律的通气减阻模拟方法，使其能真正模拟通气减阻，并对其实用化提供理论依据。

1 研究对象及多孔介质模型 1.1 研究对象

本文的研究对象为一个直径为80 mm的轴对称航行体^[8]：头部为半球，中间圆柱段长100 mm，尾部近似为长半轴等于60 mm的半椭球（见图1）。距离头部驻点15～20 mm处有2排孔径为0.4 mm呈放射状排列的小孔（每排36个，共计72个），小孔的中心轴线两两夹角10°与航行体轴线相交于2个横截面上。计算域如图2所示，为长宽高1200 mm×600 mm×500 mm的长方体域：水深400 mm，向上100 mm为空气层，航行体轴线位于240 mm水深。

假设来流速度均匀且沿航行体轴向，坐标系定义如下：坐标原点位于航行体头部驻点，x方向指向前，y方向铅垂向上，z方向按右手定则。由于流场关于xy坐标面对称，仅取一半流场（从首部看，取右侧）进行计算，边界条件设置如图3所示，位于头部的局部通气区域如图4所示。

1.2 多孔介质模型

所谓多孔介质就是固体物质内部和表面有许多孔隙，多孔介质模型（Porous Media Model）是将多孔结构简化为一个动量源，在建立几何模型时，可以不用建立复杂的几何结构，而主要通过在动量方程中添加动量源项来考虑其对流动阻力产生的影响。多孔介质分为各向同性和各向异性2种，前者指多孔介质各个方向的阻力相同，后者指各个方向阻力不同，有的方向流体容易通过，有的方向流体很难通过。式(1)以x方向（方向1）上的动量源项为例给出PMM模型表达式，其中第1项为粘性损失项，第2项为惯性损失项：

$ {S _1} = - {D_1}\mu u - {C_1}\frac{1}{2}\rho \left| u \right|u 。$

(1)

式中：μ为流体黏度；ρ为流体密度；u为方向1速度分量；D₁为方向1上的粘性阻尼系数；C₁为方向1上的惯性阻尼系数。

本文中的航行体为轴对称体，共72个小孔位于垂直于轴线的2个相距5 mm的横截面上，多孔区域如图1阴影及图4阴影所示。在此使用圆锥指定方法，方向1为锥轴方向，方向2为径向，方向3为圆周方向。通气主要沿径向由航行体内部通向头部表面，因此为显著的各向异性。当多孔介质为各向异性时，2个方向的阻尼系数若相差较大会引起收敛问题，因此一般两者取值相差不超过10³倍。

2 基本算例及流场模拟

本文的航行体速度计算范围从0.5～2.0 m/s，体积通气率为10 L/min（质量通气率为0.001021 kg/s），通气入口如图4所示为内径40 mm的圆柱面。由于通气速度较小，忽略惯性阻尼系数，仅考虑粘性阻尼系数，并将径向粘性阻尼系数D₂设为2×10¹⁰ m⁻²，另2个方向则设为2×10¹³ m⁻²。由于考虑自由面的影响，这里第一相为液体水，第二相为空气，第三相为由航行体内部通向头部表面的可压缩理想气体。图5以来流速度v_∞=1.0 m/s为例，给出流场数值模拟结果。

可见，基于PMM的三相流流场模拟可以获得丰富的流场信息：由航行体内部喷出的高压气体（图5(a)），经由多孔区域一路减压至航行体表面（图5(b)），在来流的作用下，向后流动（图5(c)）；通气流场受重力影响显著，表现在流场上下不对称−气体会上浮，向后流动的同时向上运动，逐渐接近自由面（α₂=0.5，图5(d)）；自由面在下方航行体的作用下则出现显著的横波和散波系（图5(e)）；由航行体尾部拖出双涡管，Q_criterion=125 s⁻²等值面（图5(f)）与α₃=0.1等值面（图5(c)）非常相似，后者在几何形态上略大于前者。

图6为来流速度分别为0.5 m/s和1.5 m/s时的自由面下通气等值面α₃=0.1的形态，结合图5(d)（v_∞=1.0 m/s）来看，随着来流速度的增加：通气等值面越来越贴近航行体表面，越来越接近来流方向向后发展；同时，航行体尾部出现上下2层尾流，但下层尾流不稳定，在高速时会发生断裂和脱落。

3 阻力与减阻率分析 3.1 PMM与传统方法的比较

当不考虑具体的通气装置（删除多孔区域porous zone，仅在航行体头部开出5 mm宽的环形开口以通气，该方法标识符为“nopz”），即使航行体表面出气口满足质量流量条件。本节分别采用传统nopz方法和PMM方法对4个来流速度下的三相流流场进行了数值模拟，并将通气状态标识为“Q₁”。非通气状态则标识为“Q₀”，PMM方法的非通气状态直接将通气率设为0得到；nopz方法的非通气状态，则通过令航行体全封闭来得到（即通气状态下nopz方法的环形开口设为“壁面”）。本节分别采用nopz方法和PMM对通气和非通气状态进行了4个来流速度下的数值模拟，计算工况如表1所示。图7给出了4个来流速度下2种方法的阻力系数计算结果（特征量取为$ \rho {v^2}S/2 $，其中特征面积S为轴对称体最大横截面积），表2给出nopz方法在4个来流速度下的通气和非通气阻力系数及其分量计算结果。表中，C_t为总阻力系数，包含2个分量：压差阻力系数C_p和摩擦阻力系数C_v。由图表可见：通气状态下来流速度为0.5和1.0 m/s时，nopz方法出现负阻力；而PMM方法4个来流速度下的阻力均为正。

表3以来流速度v_∞=0.5 m/s为例，给出了2种方法的通气阻力系数及其分量计算结果。可见，2种方法在航行体固壁各部分阻力系数及其分量的计算上，结果相差其实很小（不到1%），总阻力却存在正负号的差异，主要原因在于nopz方法忽略了通气装置的压差阻力（vent代表图4所示2个侧壁），从而导致符号相异这一巨大差异。

3.2 PMM减阻率分析

表4给出 PMM方法在4个来流速度下的通气和非通气阻力系数及其分量计算结果，由表可见PMM方法成功解决了不考虑通气装置时阻力为负的问题；可以成功地模拟较大来流速度范围内的通气减阻，即PMM方法可以定性正确地分析通气减阻。

由表4和图7可见PMM方法的4个来流速度仅有1个例外，在代表低速的0.5 m/s工况下，PMM方法并不减阻。图8(a)为通气和非通气状态下2条阻力系数求解的时间历程曲线，可见通气状态下阻力系数呈规则局部振荡，其下包络线几乎沿着非通气状态阻力曲线。这种振荡是由于航行体首尾压力周期性振荡引起的。为此在航行体头部和尾部驻点上各设置一个压力监测点，图8(b)和8(c)分别为2个监测点的压力时历，显而易见它们都具有规则局部振荡的特性。进一步对3个时间历程进行快速傅里叶变换后，得到三者的频率特性如图8(d)所示：首尾压力振荡与阻力系数的振荡完全同频，基频在15.65 Hz左右。结合图6(a)来看，低速下，气体受重力影响向上的运动分量较大，气泡脉动较大，引起头部和尾部压力振荡，且航行体表面包裹的气层面积相对较小，摩擦阻力的减小量不及头尾压差的增大，导致通气减阻效果不明显。

PMM之所以具有此优势，进一步分析可知此方法中，通气装置2个侧边壁在低速时产生阻力，高速时产生推力。下面以v_∞=1.0 m/s代表低速，v_∞=2.0 m/s代表高速为例，给出侧边壁上的压力分布如图9所示。为方便对比，图9将2个垂直于主流方向的侧边壁压力分布在同一张图上进行显示，右侧对应靠近首部的较小截面（图4所示通气侧壁1），左侧对应流向较大截面（图4所示通气侧壁2）。可见，v_∞=1.0 m/s时，前后侧壁上的压力分布沿径向分布比较相似，仅后侧壁上多出环状的一圈，所以产生正向阻力；v_∞=2.0 m/s时，前后侧壁上的压力分布沿径向差别比较大，相同半径上前侧壁的压力显著大于后侧壁，故产生负阻力。

3.3 非通气状态计算误差分析

本节研究非通气状态下，即PMM方法在通气率为0情况下的阻力计算结果与传统计算结果之间的相对计算误差。由图7可见，非通气状态4个来流速度下，分别采用nopz方法和PMM方法所得的阻力曲线呈现出平行的趋势；由表2和表4中的具体数值，可以定量地分析2种方法的相对误差，表5为2种方法所得阻力系数及其分量之间的相对误差百分比。由表可见各个速度下，PMM相较于nopz方法：压阻平均增大37.0%，摩阻平均降低6.7%，总阻力平均增大19.3%。可见，关于PMM的定量分析能力，在非通气状态计算结果与航行体全封闭阻力计算结果之间，存在确定性的系统误差，哪个更接近物理真值尚需实验的检验。换言之，如若PMM更接近实验值是最好，如若nopz方法更准确，则可通过确定的系统误差进行换算。

3.4 通气状态下阻力影响因素分析

本节固定来流速度v_∞=1.0 m/s，在之前采用径向阻尼系数D₂=2×10¹⁰ m⁻²的算例基础上，分别向上调节阻尼系数至2×10¹¹ m⁻²，向下调节阻尼系数至2×10⁷ m⁻²，以此来研究多孔介质阻尼系数对航行体总阻力的影响。

图10(a)为总阻力系数C_t随阻尼系数D₂的变化，图中曲线采用4次多项式进行拟合。可见C_t存在一个最低点，正对应着基准算例D₂=2×10¹⁰ m⁻²，大于或小于这个值阻力都有所增大，尤其是当D₂<2×10⁸或D₂>1×10¹¹后，通气阻力将高于非通气阻力（图中点划线），减阻效果消失。由图10(b)可见，总阻力最低点对应着通气出气口平均静压（$ p/(0.5\rho v_\infty ^2) $）最低点、平均速度（$ \left| v \right|/{v_\infty } $）最高点，此点上的通气相分数（$ {\alpha _3} $）也接近于1，几乎为纯粹通气相。也就是说，出气口压力低（头部压阻小），通气更顺畅（平均速度及气相分数高），通气减阻效果也越好。这里存在一个阻尼系数最佳值，取值过高或过低都将影响减阻率。本节还探讨了孔隙率对阻力计算的影响。这里的孔隙率表示流体所占据的体积分数，影响动量和能量方程，此前所有的计算均取值σ=0.3。本节另外计算了0.2和0.4两个数值，发现阻力所受影响甚小。由以上探讨可知，通气状态下的阻力定量分析在获得实验数据后可通过阻尼系数的调节来实现精度要求。后续的数值模拟研究将围绕此展开。

4 结　语

传统忽略通气装置的通气减阻模拟会出现负阻力，由此可知忽略通气装置的数值模拟存在缺陷，有必要对其进行改进。

为了考虑通气装置对通气减阻率的影响，本文将多孔介质模型应用于通气减阻的模拟，并同时考虑重力和自由面的影响。通过对多孔介质模型中的阻尼系数、孔隙率等参数进行合理设置，成功对自由面下航行体通气三相流流场进行了数值模拟，并得出了合理的阻力和减阻率数据。

PMM相较于传统方法，在非通气状态两者之间存在确定性的系统误差，总阻力偏大约19.3%；通气状态下的定量分析则可在未来获得实验数据后，通过阻尼系数的调节来满足精度要求。

基于本文提出的PMM方法，可进一步开展通气率、通气位置及航行体外形对减阻率影响，并开展水下高速航行体优化设计等研究。