基于传质源项的水翼空蚀风险数值预报

郑恩慧，曹彦涛，程怀玉，季 斌，彭晓星

（1.中国船舶科学研究中心船舶振动噪声重点实验室，江苏无锡 214082；2.武汉大学水利水电学院，武汉 430072）

0 引 言

随着现代舰船吨位和航速的不断提升，空化现象也愈发常见。空化往往会显著增大舰船推进器的振动、噪声，引起性能下降，严重时甚至会直接造成推进器叶片表面材料的剥蚀，引起空蚀，严重危害舰船的安全。在推进器的设计过程中，需要在保证推进效率的同时，尽可能减小或避免空蚀的影响。因此需要一个准确客观的空蚀风险预报方法，以获取空蚀可能出现的区域和强度，从而为设计方案的优化提供依据。

传统推进器设计中主要依靠模型涂层试验来判断空蚀风险的潜在位置及强度，即通过在模型表面喷涂软面涂层开展设计工况下的空化试验，利用涂层剥离区的状况判断空蚀风险。Cao 等[1]利用涂层试验得到了两种工况下NACA 0009 三维扭曲水翼的空蚀风险区域；Li 等[2]分别对8°攻角NACA 0015 水翼及6.5°攻角NACA 0018-45 水翼进行设计工况下的涂层试验。但采用试验方法研究空蚀风险问题通常面临着试验周期长、费用高、数据有限等问题。近年来随着CFD 的不断发展，空化流动数值模拟技术日趋成熟，其在推进器空化性能的设计评估中得到了广泛的应用[3-7]。在可靠的空化流场数值模拟结果基础上，对空蚀风险进行数值预报也成为了可能。目前，研究者们通常以空化流场压力、蒸汽体积分数等特征参数为依托，通过一些表征参数构建空蚀风险预报的方法。这些方法近年来获得了较大发展，在空化流场空蚀风险的评估中得到了一定的应用。Dular[8]将蒸汽体积分数的标准差与空蚀风险联系起来；Nohmi[9]通过观察数值模拟中的汽泡行为，简单定义了固壁压力及蒸汽体积分数的函数来衡量空蚀风险；Ochiai[10]提出了离散泡法，其空蚀风险参数是通过R-P 方程计算微观汽泡内部的压力来估计的。

准确可靠的数值模拟结果是开展空蚀风险数值预报的基础。对于空化流动的数值模拟而言，当前基于均质平衡流模型框架的质量输运方程空化模型得到了广泛应用[11-15]。结合改进型RANS、DES或LES等湍流模拟方法[16-20]，可以较好地模拟非定常空化流动的主要特征。

而在空蚀风险表征参数的构建方面，当前存在两种主要研究思路：一是精细化捕捉空化泡群溃灭特征，利用欧拉-拉格朗日方法精细化捕捉空化溃灭的瞬态载荷特征，但该类方法所需计算资源巨大，当前尚未达到工程实用化的程度[21-25]；另一种思路是抛开复杂的流动细节，基于宏观流场的参数构建能够反映空蚀风险特性的表征方法。其中，以势能假设为基础的空蚀风险表征参数方法已经获得较大发展。Patella[26]提出了空泡溃灭后冲击波的能量传递机制；Li[2]提出了基于压力时间导数的空蚀风险表征参数；Leclercq[27]利用局部势能变化率结合立体角模型计算了壁面吸收到的辐射能量；Melissar⁃is[28]假定空泡结构的初始势能在溃灭前首先转化为空泡界面处的动能，随着空泡体积的减小，动能逐渐向空泡中心聚焦，一旦达到溃灭条件，全部能量以压力波的形式向外释放，压力波以投影的方式辐射到壁面上，从而获得壁面上的空蚀风险区域；Mohammad[29]则直接计算空泡溃灭过程中周围液体的动能，避免了势能计算过程中驱动压强的不确定性问题，并且同时考虑了冲击波与微射流两种空蚀作用机制。

上述研究加深了人们对于空蚀风险预报理论与方法的认识，为本文的研究提供了非常有意义的参考。本文基于单泡溃灭模型，分析影响空泡溃灭的因素，在当地压强大于饱和蒸汽压及压强变化率大于某个阈值的前提下，利用传质源项计算一段时间内由于空泡溃灭引起的壁面上的能量累积，通过为壁面累积能量划定阈值确定空蚀风险区域的范围，通过壁面累积能量的大小判定空蚀风险的相对严重程度。通过该方法预报NACA 0009三维扭曲水翼的空蚀风险结果，并与试验结果进行比较，验证该方法的准确性。

1 水翼空化演化过程数值模拟

1.1 控制方程与大涡模拟方法

本文空化流动计算采用的两相流模型为mixture 模型。控制方程为连续性方程和动量方程，以张量形式表述如下：

式中，ui是混合物在i方向的速度，p是汽液混合物的压力。混合物的粘性μ与密度ρ视为蒸汽体积分数αv的函数：

其中，下标v和l分别代表蒸汽相及液体相。

为了能较好地模拟各种湍流尺度及更好地捕捉涡旋结构，本文采用LES 方法模拟湍流。大涡模拟的控制方程为

式中，应力项τij用亚格子模型求解。通过WALE 模型模拟亚格子尺度湍流特性。该模型在此类空化流动中已经得到了较为广泛的应用，取得了令人满意的预报结果[16,18]

1.2 Schnerr-Sauer空化模型

本文采用的空化模型为Schnerr-Sauer(S-S)空化模型[15]，汽、液两相间的质量传输由蒸汽体积分数输运方程来控制：

式中：αv表示蒸汽体积分数；ρl为液体密度；ρv为蒸汽密度；ṁ为汽液两相间的净质量输运源项，又分为蒸汽生成项ṁ+及凝结项ṁ-，表征如下：

其中：pd为泡外驱动压力；pv为饱和蒸汽压；n为汽泡数密度，本文取n=1×1013。

1.3 计算域及网格划分

本文采用的水翼模型为三维扭曲水翼，其截面翼型为NACA 0009。弦长为112.5 mm，展长为225 mm，左右对称，水翼中部攻角为11°，两端攻角为0°，水翼模型的外形如图1 所示。图2 给出了数值模拟中使用的计算域及相应的边界条件的设置，C为水翼弦长，边界条件采用速度进口与压力出口，进口速度为14 m/s，通过调节出口压力值使空化数等于1.2，空化数定义为

图1 水翼模型Fig.1 Hydrofoil model

图2 计算域及边界条件Fig.2 Computational domain and boundary conditions

需要注意的是，为了节省计算资源，本文在数值模拟中仅对该水翼的一半进行了模拟计算，中间切面选取对称面边界条件，其余壁面均选取无滑移壁面条件。

图3 给出了计算域内的网格划分情况。为保证网格质量，全域均采用结构化网格。为精细地捕捉水翼吸力面的空泡结构演变特性，在水翼吸力面附近进行网格细化，在水翼壁面附近添加边界层，使得壁面函数y+在1 左右。为兼顾计算效率与计算精度，参考Long 等[30]的工作，本文的总网格数取1100万，非定常计算时间步长选取为Δt=1×10-5s。

图3 网格划分情况Fig.3 Mesh generation

1.4 计算结果与对比分析

以数值模拟手段对空蚀风险预报的前提是该模拟获得的空化流场形态应当与实际一致。为了验证空蚀风险预报结果的准确性，本文在涂层试验过程中使用高速摄影拍摄水翼吸力面的空化形态。试验工况与数值模拟一致。

图4 为一个脱落周期（T0）内水翼吸力面试验空化形态与数值模拟蒸汽体积分数αv=0.1 的等值面图像的对比。结果显示，试验和数值结果较一致地展现了片空化的产生、发展、脱落及溃灭的准周期过程。该工况下空化形态主要特点为：水翼两侧流动状态相对稳定，片空化的脱落、再生过程集中在水翼中央；片空化充分发展时被向上游发展的回射流从导边处切断形成脱落云，见图4（a）～（b）；随后脱落云中间部分向上卷起，两端附着在水翼表面，逐渐形成U型涡，见图4（c）～（d）；脱落云与U 型涡随主流向下游运动过程中，伴随着大量空泡溃灭，直到U 型涡离开随边时，空泡几乎全部溃灭，见图4（e）～（f）。

图4 水翼吸力面试验中的空化形态与数值模拟αv=0.1 空化形态Fig.4 Cavitation morphology captured by test and iso-surface of αv=0.1 captured by numerical simulation on hydrofoil suction surface

2 空蚀风险涂层试验

2.1 试验设备及试验方法

本次试验主要设备包括循环式高速空泡水筒、高速摄影相机等。图5为高速空泡水筒的示意图。高速空泡水筒试验段的长度为1600 mm，宽和高均为225 mm，最高水速为25 m/s，试验段压力可调整为5 kPa至500 kPa，边壁为透明有机玻璃，以便高速摄影拍摄试验情况。高速摄影采用3200 fps的拍摄频率，对水翼的吸力面进行拍摄。

图5 高速空泡水筒示意图Fig.5 Schematic diagram of high speed cavitation tank

水翼模型参数与数值模拟模型完全相同，铝合金材质，在水翼表面喷涂软面材料，便于更快地观测空蚀风险情况，试验前水翼吸力面如图6所示。本次水翼安装采用0°攻角，吸力面朝下，在试验段的安装位置如图7所示。

图6 试验前水翼吸力面Fig.6 Hydrofoil suction surface before test

图7 水翼安装位置Fig.7 Installation position of hydrofoil

试验水速为14 m/s，试验段空化数为1.2。涂层试验始终在该工况下进行，试验开始后，空化形态很快趋于稳定。试验过程中，通过高速摄影判断水翼空蚀风险的大致情况，当发现高速摄影中出现明显的空蚀风险区域时，停止运行空泡水筒，再在无空化情况下确认空蚀风险区域是否形成。本次试验进行1 h后，发现水翼吸力面出现了大量的空蚀风险区域，随即停止试验，将水翼从高速水洞中取出，观察并记录水翼表面空蚀风险情况。

2.2 涂层试验空蚀风险结果及分析

空蚀风险试验结果如图8 所示。其中，黑色区域为涂层材料受到剥蚀的区域。图8（a）中标出了空蚀风险区域，图9则给出了相应区域的高速摄影典型空化形态。

图8 涂层试验空蚀风险结果Fig.8 Results of cavitation erosion risk

图9 试验空蚀风险结果与空化形态的对应Fig.9 Correspondence between cavitation erosion risk results and cavitation morphology

涂层试验结果显示，空蚀风险区域主要分为两部分：主空蚀风险区域及条带状空蚀风险区域。其中主空蚀风险区域大约处于水翼弦向2/11弦长到8/11弦长之间，其空蚀风险最为严重，对应于片空化脱落区域。与主空蚀风险区域相连接的两个条带状的空蚀风险区域一直延伸到随边，存在明显的点状空蚀风险区域，对应于U 型涡两腿经过的区域。此外，U 型涡脱离开随边的位置处，空蚀风险明显加重，对应于条带状空蚀风险区域与随边相连的位置，如图9（b）中红色箭头指向的位置。

3 空蚀风险预报方法

3.1 空蚀机理回顾

空泡向水翼下游移动进入压力恢复区域时，由于泡外压力不断增大，空泡开始溃灭。人们普遍认为，空蚀是蒸汽空泡溃灭引起的。对于初始半径为a、初始泡壁速度为u0=0 的蒸汽空泡，不考虑表面张力、粘性、不可凝气体在内的溃灭动力学方程，即R-P方程为

式中，r0为汽泡半径，pd(t)为泡外压强，pv为饱和蒸汽压，ρl为液体密度。式（10）的解，即溃灭过程泡壁速度为

空泡收缩至溃灭的时间为

显然，当空泡半径r0→0，泡壁面流体沿径向速度u0→∞。实际流动中，由于粘性泡内存在不可凝结气体，使得空泡不能溃灭至半径为0，泡壁面速度也不能达到无穷大。尽管如此，空泡溃灭至最小半径时，其泡壁速度也将是很大的数值，这实际上就是造成空蚀的微冲击波及微射流的来源。

当然，只有非常靠近壁面的空泡溃灭才对空蚀有贡献，否则微冲击波与微射流在尚未到达壁面时将衰减耗尽。此外，无论是球形绕流体近壁面的顺流而下，还是球形泡的收缩都会产生趋壁效应[31]。定性而言，趋壁效应使得近壁面侧压力进一步降低，原本中心对称的冲击波在靠近壁面方向的强度会大于其他方向，而微射流的方向也会指向壁面，使得向壁面方向传递的能量更加集中。

所以，决定空蚀风险的主要是近壁面的流场，本文提出的空蚀风险参数，只需在贴壁流场中去判别与寻找。Phillip 和Lauterborn[32]也通过试验证明了最大的冲击载荷是由与壁面直接接触的空泡溃灭造成的。

3.2 蒸汽泡溃灭的能量转化过程

根据前述，本文两相流的数值模拟方法与三维扭曲水翼空化流试验结果吻合较好。但该两相流数值模拟无法直接显示蒸汽泡溃灭引起的径向流场速度u0的时空分布，更无法表征微冲击波与微射流。为此需要寻找其它表征溃灭的物理参数。

3.2.1 蒸汽泡势能

将半径为r0的单个蒸汽球泡与周围半径为R（R→∞）的球状液体空间组成一个系统，如图10所示。由于系统蒸汽泡内部气压为饱和蒸汽压Pv，而系统外则为大于饱和蒸汽压的当地压力Pd(t)，因此，在t时刻系统内外存在压差Pd(t)-Pv。

图10 单空泡溃灭模型Fig.10 Single bubble collapse model

收缩过程中，蒸汽球泡泡径为r0时，若泡径变化dr0，pv对液体作负功为

泡外充满水的充分大半径为R的同心球体因连续性方程同时需缩小dR，有

pd(t)对液体作正功为

故泡收缩时，外力对液体实际作功为

式中，Vv为蒸汽泡的体积。

因此，当蒸汽泡收缩dVv时，外力对泡外液体作功dW；而蒸汽泡损失的势能为

蒸汽泡溃灭释放出全部势能

3.2.2 溃灭过程中蒸汽泡的势能转化为流场的动能

初始泡壁速度为0的空泡，当半径由a缩减至任意值r*≤a时，空泡的势能减小为

式中，u*为泡壁r*处的速度。

当空泡半径由a缩减至r*时，流场动能增加为

则ΔEp=-ΔEk。可知蒸汽泡减小的势能全部转化为周围流场的动能。所以，寻求蒸汽泡溃灭所引起的流场动能的时空分布，可以等价转换为寻求流场蒸汽泡势能减少的时空分布。

3.3 空蚀风险表征

本文采用的S-S 空化模型中选用单位体积液体气泡数密度为n=1×1013，空蚀风险是群泡的物理行为，就是对单位体积中存在的1×1013个蒸汽泡的时空特性的描述。

3.3.1 空蚀风险表征参数的选取

（1）环境压强

式（10）可改写为

（2）环境压强变化率

（3）壁面累积能量

由于目前数值模拟尚不能给出空泡群溃灭过程准确定量的描述，且空蚀风险预报无需定量确定壁面用于剥蚀能量的大小，只需要一个相对值。但有一点可以肯定的是，近壁面流场释放的空泡溃灭的能量越大越多，壁面累积的用于空蚀的能量便越多，空蚀就越容易发生，空蚀程度就越严重。

式（18）给出了单泡溃灭时势能释放的表达式，则对于群泡，在单位体积中释放的势能变化率为

式中，αv为单位体积的蒸汽体积分数，ep为单位体积泡群释放的势能。

因为空蚀是微冲击波与微射流长期作用在壁面上的结果，参数一、二的满足可视为单次冲击具备了引起空蚀风险的潜力，但需要多次类似的冲击叠加才可能导致空蚀。考虑到空化流的非定常特性，故选取一充分长的时间T，将壁面累积的空泡溃灭过程中释放的势能值超过阈值C2作为参数三。

3.3.2 空蚀风险表征方法

（3）壁面空蚀风险表征参数三为Ic.e＞C2。

根据前人的研究，本文在计算壁面累积能量Ic.e时采用了时间平均压力场pˉd(t)作为驱动空泡溃灭的环境压力场[28]。

需要说明的是，式（26）中蒸汽体积分数变化率项若使用物质导数的展开式进行计算，其时间偏导数项及对流项均需采用差分格式，若无法保障数值计算结果的时间步长足够小及网格分辨率足够低，将会对计算精度产生较大影响。本文采用传质源项对蒸汽体积分数变化率进行替代，式（6）将蒸汽体积变化与相变过程中的传质源项联系起来，传质源项的表达式见式（7），可发现其全部变量数据取自当前时间步，且无需求解对流项，因此解决了差分格式带来的误差问题，相对而言，精度更高。根据式（6）～（7），溃灭过程中的蒸汽体积分数变化率可表示为

在满足当地压强大于饱和蒸汽压，压强随时间变化率大于阈值C1的条件后，计算壁面累积能量Ic.e，壁面累积能量大于C2可视为该区域存在空蚀风险。

3.4 空蚀风险预报结果

3.4.1 阈值C1、C2的选定

阈值C1、C2作为压强随时间变化率及壁面累积能量的下限，共同决定了空蚀风险区域的范围。在此基础上，壁面累积能量Ic.e的相对大小则可表征空蚀风险的严重程度。

分别选取阈值C1=2×105、1×106、5×106及C2=1×1010、2×1010、3×1010，将阈值C1、C2两两配对计算水翼表面的空蚀风险。由此得到9个空蚀风险预报结果，并将其与试验空蚀风险区域进行比较，如图11所示。可以看出，随着阈值C1的增加，水翼中前部的空蚀风险区域逐渐缩小，中后部空蚀风险越来越明显；随着阈值C2的增加，空蚀风险区域的外围轮廓，尤其是水翼导边附近区域随之缩减。通过比较，与试验结果吻合最好的为图11（e）。因此选定阈值C1=1×106，C2=2×1010。

图11 9组阈值计算的空蚀风险预报结果与试验结果比较Fig.11 Prediction results of cavitation erosion risk calculated by 9 groups of threshold values compared with the experimental results

3.4.2 空蚀风险预报结果

图12 为采用基于传质源项的空蚀风险预报方法得到的空蚀风险预报结果与试验结果的对比。可以看出，空蚀风险预报结果在片空化脱落的区域预报出了空蚀风险，同时在U型涡经过的位置也预报出了两个条带状的空蚀风险区域，与试验结果显示的空蚀风险区域基本吻合。此外，预报结果显示片空化脱落区域的空蚀风险程度更为严重，与试验结果亦较为吻合。尽管预报结果与试验结果在大部分区域吻合较好，但预报结果显示在水翼两端即红框范围内存在空蚀风险区域，与试验结果不同。数值模拟和试验显示在水翼两端的空化形态为片空化的周期性波动，与水翼中央出现的云空化不同，片空化周期性变化引起的空蚀风险预报方法需进一步研究。

图12 空蚀风险预报结果与试验结果的对比Fig.12 Comparison between cavitation erosion risk prediction results and painting test results

4 结 论

本文以NACA 0009 三维扭曲水翼为数值模拟及涂层试验对象，将数值模拟空化形态与试验高速摄影结果进行了对比。基于单个蒸汽球泡的溃灭模型，分析了溃灭过程中能量转化过程及相质量的转化过程，提出了空蚀风险的三个表征参数，并合理选定了表征参数的阈值。基于该方法，本文对绕NACA 0009 三维扭曲水翼表面的空蚀风险进行了预报，并与涂层试验结果进行了对比分析。主要结论如下：

（1）采用大涡模拟方法可较好地捕捉绕NACA 0009 三维扭曲水翼空化流动的初生、发展、脱落、溃灭等过程，与试验结果吻合较好；

（2）提出空蚀风险预报的三个表征参数：当地压强、当地压强变化率及壁面累积能量Ic.e。在当地压强大于饱和蒸汽压及当地压强变化率大于阈值C1的前提下，通过壁面累积能量Ic.e大于阈值C2可以划定空蚀风险的区域，通过Ic.e的大小可以判定区域内空蚀风险的相对严重程度。其中，Ic.e的求解采用了基于传质源项的表达，避免了数值误差。

（3）通过与涂层试验比较，本文建立的空蚀风险数值预报方法获得的空蚀风险结果与试验空蚀风险结果在大部分区域吻合良好。两者显示的空蚀风险区域主要分为两部分，一是位于片空化脱落的区域，二是位于U型涡两腿经过的区域，其中片空化脱落区域空蚀风险更为严重。

虽然本文空蚀风险预报仅对三维扭曲水翼进行了验证，而且阈值C1、C2的确定带有经验性质，但因为三个空蚀风险参数均建立于空蚀机理之上，所以具有一定的普适性。下一步将研究如何将阈值C1、C2与影响空蚀的物理量联系起来，从理论上确定阈值C1、C2。