亚格子模型对三维扭曲水翼空化现象的影响
时间：2023-02-13 16:50:09 来源：千叶帆 本文已影响人
何朋朋，李子如，张孝旺，贺伟*
1 武汉理工大学 高性能舰船技术教育部重点实验室，湖北 武汉 430063
2 武汉理工大学 船海与能源动力工程学院，湖北 武汉 430063
在船舶领域中，螺旋桨空化不仅会恶化推进性能，而且还会伴随产生强烈的振动、噪声及剥蚀等问题。由于螺旋桨为三维扭曲形式，在不同半径位置处的空化类型各异，其中针对非定常云状空泡初生与溃灭机制的探讨，对于螺旋桨空蚀及空化噪声的合理预测至关重要。为了更好地观察及研究非定常云状空化的演变过程，许多学者对展向攻角渐变的三维扭曲水翼进行了数值与试验研究。
例如，文献[1–3]分别通过LDV（laser doppler velocimetry）技术、高速摄影技术和PIV（particle image velocimetry）技术在空化水洞中对三维水翼的空泡形态、攻角及空化数、空泡初生及分离现象开展了试验研究，从结果观察到了马蹄形云空泡，发现了空化流动特性主要与水翼攻角及空化数有关。文献[4]以Delft Twist11N 三维扭曲水翼为对象，引入了侧射流（side-entrant jet）的概念，基于试验结果分析了侧射流及主回射流（re-entrant jet）在空泡脱落中的作用。文献[5–8] 分别结合经过湍流黏度修正的SSTk-ω 和RNGk-ε 等湍流模型和Schnerr-Sauer（S-S），Singhal 及Zwart 空化模型，对Twist11N 水翼的空化现象进行了数值模拟，并使用经过黏度修正的雷诺平均N-S （RANS）方法模拟了大尺度空泡脱落现象，但空泡脱落频率与试验结果的误差较大，且云空泡形态明显小于试验结果，其主要原因是RANS 方法对所有尺度的湍动特征参数均进行了雷诺平均，对空泡精细结构和非定常瞬态特性的描述略显不足。文献[9]针对三维水翼梢涡空泡的初生问题，分别使用SSTk-ω 湍流模型、分离涡（detached-eddy simulation，DES）和大涡模拟（large-eddy simulation, LES）模型对水翼的梢涡流场进行了计算研究，结果表明，基于LES 方法的计算结果与试验值吻合较好。文献[10] 基于解析精度较高LES 方法，对NACA 66 水翼的非定常空化湍流进行数值研究，结果证实了压力脉动与空泡脱落过程有关，且空泡体积加速度是造成空泡水翼周围压力脉动的主要原因。文献[11] 使用两套不同密度的网格研究了网格分辨率对三维扭曲水翼周围云空泡流动数值模拟的影响，结果表明，对于空化流动数值模拟，细网格数值模拟结果的精度比粗网格的有显著提高。文献[12] 从拉格朗日观点对三维扭曲水翼的空化流动进行分析，结果更直观显示出回射流和侧射流的轨迹线。文献[13]对NACA 66水翼的通气空化流动进行数值研究，通过与试验结果对比验证了数值方法的有效性。总之，虽然基于LES 方法的三维扭曲水翼空化数值模拟已展开了较多的研究，但目前对于网格划分及亚格子模型的选择大多还依赖研究人员的经验，缺乏客观依据。
综上，针对使用LES 方法的三维扭曲水翼空化数值模拟中网格密度及亚格子模型的适应性问题，本文将基于LES 的黏性多相流数值方法对围绕Twist11N 扭翼型水翼的空化现象进行研究，重点探讨不同网格密度及亚格子模型对此类水翼空化流特性的预报能力和适应性，分析该水翼云状空化现象的非定常瞬态特性，期望研究结果可以在基于LES 方法的三维水翼及螺旋桨空化预报中为网格划分及网格子模型的选取提供参考，并有助于剖析船舶螺旋桨中的非定常空化流场。
1.1 控制方程
本文将包含汽液两相的空化流视为密度可变的单一流体，运用基于均质平衡流理论的均质混合模型来描述汽液两相。汽液混合空化流场的连续性方程和动量方程如下：
式中，下标m，v，l 分别指汽液混合相、汽相、液相。
1.2 大涡模拟方法
空化流动是一种复杂的湍流，本文使用具有较高解析精度的LES 方法对湍流进行处理。LES方法利用滤波器对Navier-Stokes（N-S）方程进行滤波，直接求解滤波后的N-S 方程，获得可解尺度物理量，而滤出的亚格子尺度物理量则通过相应的亚格子模型进行求解。对于不可压流动，经滤波处理的连续性方程及动量方程如下：
1.3 空化模型
采用基于质量输运方程的S-S 空化模型来描述汽液两相间的转换，该模型是以简化的Rayleigh-Plesset 气泡动力学方程为基础而建立的。
式中，nb为不气泡数密度，一般取值为1013。
本文以Twist11N 扭曲型水翼为研究对象，其展向不同位置处的截面形状均为NACA 0009 翼型，但经过了变攻角设计，从两端向中部的展向截面攻角变化范围为–2°～9°，具体详见文献[4]。图1 所示为Twist 11N 水翼不同角度的视图。其中：坐标原点位于水翼端部截面的弦长中点，XC表示入流方向，YC表示翼展方向，ZC表示垂直向上；
LE 与TE 分别表示水翼的导边和随边；
弦长C=0.15 m，展长l=0.30 m。
图1 Twist11N 水翼几何图Fig. 1 Geometrical views of Twist11N hydrofoil
考虑到Twist11N 水翼几何模型具有对称性的特点，计算时仅取水翼的一半。图2 所示为计算域示意图。其中，侧壁处翼截面弦长的中点与上下壁面的距离分别为C，翼展方向为C，入口至导边为2C，出口至导边为5C。
图2 一半Twist11N 水翼的计算域示意图Fig. 2 Schematics of computational domain for half-span ofTwist11N hydrofoil
图3 Twist11N 水翼周围O 型分块网格的划分形式Fig. 3 Grid scheme of the O-shaped block around Twist11N hydrofoil
表1 Twist11N 水翼网格划分参数Table 1 Grid parameters of Twist11N hydrofoil
数值研究使用ANSYS FLUENT 19.0 软件。Twist11N 水翼在空化数 σ=1.07 时，非定常空泡流数值计算的边界条件设置如下：入口采用速度入口条件，入口速度V∞=6.97 m/s；
出口条件为压力出口，压力值由空化数公式σ=(pout−pv)/0.5ρlV∞2计算得到，出口压力pout=29.0 kPa；
计算域顶部、底部、侧壁设置为自由滑移壁面，水翼表面为无滑移壁面，翼展中剖面为对称面边界条件。汽液两相的密度、动力黏度系数及汽化压力等参数参照试验温度T=24℃取值。参照文献[17]中关于时间步长选取方法（∆t=Tref/200,Tref=C/V∞），计算得到的时间步长为∆t=1.076×10−4s，最终选择时间步长为∆t=1.0×10−4s。权衡计算精度和效率，每个时间步中内迭代步数设置为100。
定义水翼的升力系数CL=FL/(0.5ρlV∞2A)、阻力系数CD=FD/(0.5ρlV∞2A)以及表征翼型上下表面压力分布的压力系数Cp=(p−pout)/(0.5ρlV∞2)，其中，FL和FD分别为水翼所受的升力和阻力，A为翼型投影面积。鉴于LES 方法计算得到的升阻力系数和压力系数具有随机振荡特性，故定义了多个周期内的时均升力系数C¯L、时均阻力系数C¯D和时均压力系数C¯p进行对比分析。定义数值计算结果与试验结果的误差为：
在进行空化数值模拟前，采用第22 届ITTC推荐的收敛比Rk对基于LES 方法的Twist11N 水翼湿流场数值模拟的网格收敛性进行判断，使用3 套不同密度网格计算得到的时均阻力系数为单调收敛，而时均升力系数为振荡收敛，满足网格收敛性要求。
3.1 网格密度对空化数值模拟的影响
基于WALE 亚格子模型的LES 方法，使用G1，G2 和G3 网格对Twist11N 水翼在σ=1.07时的空化现象展开数值模拟，分别探讨不同网格密度对空泡脱落频率、云状空泡形态和水翼升阻力系数脉动特性等的影响。
使用3 套网格计算得到空泡体积时历曲线，经傅里叶变换得到空泡脱落频率f，如图4 所示。结果显示，3 套网格计算的空泡脱落频率分别为30.14，31.78，31.24 Hz，与试验结果[4]（f=32.55 Hz）的误差分别为−7.40%，−2.37%，−4.02%。由此可见， G1 粗网格计算得到的空泡脱落频率与试验值的误差较大，而G2 与G3 计算得到的空泡脱落频率误差明显降低，且在现有计算设置下，网格的进一步加密并没有使空泡脱落频率的误差减小，反而更加偏离试验值。
图4 不同网格密度下计算的空泡脱落频率Fig. 4 Shedding frequency of cavity for various grid schemes
图5 给出的是一个周期内不同网格密度下预报的空化演变过程，采用汽相体积分数α=0.1等值面呈现的典型时刻（1～9）云状空泡形态，并与试验结果[4]进行对比。在所示俯视图中，来流方向为由上至下。
图5 不同网格模拟的空化演变过程与试验结果对比Fig. 5 Comparison of cavitation evolution by various grid schemes and experimental results
由图可见，与试验结果相比，3 套网格基本都能模拟出非定常云状空化演变过程的典型特征。具体过程包括：
1） 主脱落：附着片空泡发展到一定范围后，发生了破裂并脱落泄出大尺度空泡云团。
2） 马蹄形云空泡形成：破裂后的附着片空泡呈现凹型结构并继续发展，而脱落泄出的空泡云团则发生卷曲，形成马蹄形结构。
3） 二次脱落：马蹄形云空泡在主流体的作用下向下游移动，而附着片空泡在“耳垂”位置发生了小尺度结构的脱落。
4） 空泡溃灭：二次脱落的小尺度空泡云团发生溃灭，附着片空泡由凹型发展成为凸型结构；
主脱落的大尺度马蹄形云空泡在下游高压区发生溃灭，同时，发展起来的附着片空泡再次发生破裂，新的周期开始。
由图5 的试验结果还可以观察到，受回射流机制的影响，在时刻1 和时刻9 即将发生主脱落的片空泡表面较不规整，其它时刻的片空泡表面较为“光滑”，3 套网格大致都模拟出了上述这种现象，且随着网格的加密，片空泡上表面的结构形态逐渐细化，更加接近试验中相应区域的空泡形态。然而，数值模拟结果与试验结果相比，仍存在一些形态上的差异。例如，在试验中观测到的片空泡呈现出条状或泡状形态，但数值模拟结果并未捕捉到这些特征，主要是因为计算得到的空泡形态是由连续相分数界定得到的。
因网格密度不同，数值模拟空化演变过程中捕捉到的空泡形态存在一定差异，主要体现在：
1） 空泡长度Lc（如图5(b)中时刻8 所示）：从图6 给出的不同网格密度下数值计算得到的Lc/C值在一个完整周期内随时间的变化情况可见，Lc生长速度随时间的变化存在减慢的趋势，这是因为片空泡发展到一定长度后，在回射流和主来流作用下翼展中部片空泡被抬升翻卷，阻碍了片空泡继续往下游方向发展。而且，由G1 网格预报的Lc/C值明显低于G2 和G3 网格相应时刻的预报值，且后二者数值结果与试验结果[4]吻合较好，其中，与G2 网格数值计算结果相比，采用G3 网格在多数时刻点的Lc/C预报值与试验结果的偏差程度更大，这从图5 中的空泡形态上也可以明显看出。
图6 不同网格密度下计算的Lc/C 值对比Fig. 6 Comparison of cavity length-chord ratio for various grid schemes
2） 二次脱落：从空泡形态演变过程可见，G1网格数值计算结果在“耳垂”位置的空泡现象较简单，G2 和G3 网格计算结果在“耳垂”位置发生多次二次脱落的现象更明显；
从图5(c) 中的时刻1和时刻8 箭头所指位置可见，在网格加密后可以捕捉到更精细的空泡现象。
3） 马蹄形云空泡尺度：对于主脱落泄出而形成的马蹄形云空泡，从图5 中3 套网格模拟结果可见，在时刻3，脱落泄出的云空泡上部都有一个V 型的空泡结构，形态也都较为相似；
在时刻4，可以明显观察到马蹄形云空泡的雏形，其中G3网格数值结果中时刻4 的马蹄形云空泡倾斜程度不同于G1，G2 网格的结果。
总体上，在时刻1～时刻4 之间，3 套网格模拟的云空泡现象较为一致。在时刻4 之后，云空泡演变过程中伴随有收缩和溃灭现象，空泡溃灭形成的小尺度空泡与网格尺度量级更接近，网格尺度不同造成了包括流场中可捕捉到的漩涡尺度等流体湍动特征在预报精细程度上存在差异。可见，在时刻4 之后，马蹄形云空泡的形态差异较明显，在其展向尺度稳定后，G1 网格模拟的展向尺度小于试验结果，而G2，G3 网格模拟的结果都与试验结果相近。
图7 和图8 分别给出了3 套网格数值模拟计算得到的一个周期内升阻力系数时历曲线。图中，以无量纲时间t/T为横坐标，T为相应网格密度下的周期。由图可见，G1 网格数值计算的升阻力脉动特征相对较少，3 套网格的升阻力瞬态脉动位置较为一致，且主要集中在t/T=[0.2,0.3]及[0.8,1.0]这两个区间范围；
3 套网格数值计算的升阻力脉动差异主要集中在后半周期，与马蹄形云空泡的溃灭现象密切相关。此外，对比3 套网格所模拟的马蹄形云空泡的演变过程，还可发现G2 和G3 网格所模拟的马蹄形云空泡溃灭及回弹的现象比G1 网格的更明显（如图5 中时刻5～时刻8），这也解释了G2 和G3 网格所预报的升阻力脉动更剧烈这一现象的原因。
图7 不同网格密度下一个周期内升力系数时历曲线Fig. 7 Time histories of the lift coefficients for various grid schemes in a shedding cycle
图8 不同网格密度下一个周期内阻力系数时历曲线Fig. 8 Time histories of the drag coefficients for various grid schemes in a shedding cycle
表2 不同网格密度下时均升阻力系数计算结果的对比Table 2 Numerical results of the time-averaged lift and dragcoefficient for various grid schemes
图9 和图10 分别给出了空泡流中不同网格密度下，展向YC=0.12 m 和YC=0.15 m 位置处水翼表面的时均压力系数−C¯p与两组试验（即Delft和EPFL 这两个不同空泡水筒所做的试验）数据的对比，其中Delft 试验水翼的尺寸是EPFL的两倍，二者以相同雷诺数进行试验。由于部分压力监测点的试验结果带有一定的不确定性[20]，在对比时，既与相同尺寸的Delft的试验结果[18]作比较，也与EPFL的结果[18]作比较。由两个图可知，在水翼展向两个截面处，不同网格密度下压力面的时均压力系数计算结果基本一致，而吸力面压力的结果略有差异。主要体现在：
图9 YC=0.12 m 展向位置处水翼表面 的分布Fig. 9 The time-averaged pressure distribution around hydrofoil surface at YC=0.12 m
图10 YC=0.15 m 展向位置处水翼表面 的分布Fig. 10 The time-averaged pressure distribution around hydrofoil surface at YC=0.15 m
1）在0 2）在0.4 ≤x/C≤0.5区域，作为片空泡末端，该区域压力沿弦长方向分布变化较大，对于YC=0.12 m 展向位置，G3 网格所预报的压力转变趋势与试验结果更接近，而对于YC=0.15 m 展向位置，G2 网格较好地模拟了x/C=0.4 处的压力转变趋势。 3）在0.5 基于中密度网格G2 的LES 方法，采用WALE，WMLES 和SL 这3 种亚格子模型对Twist11N 水翼σ=1.07时的空化现象展开数值模拟，分别探讨不同亚格子模型对空泡脱落频率、云状空泡形态以及水翼升阻力系数脉动特性等的影响。 图11 给出了3 种亚格子模型计算所得的空泡脱落频率f。WALE，WMLES 和SL 这3 种亚格子模型计算得到的空泡脱落频率分别为31.78，29.19 和34.87 Hz，与试验结果[4]（f=32.55 Hz）的误差分别为−2.37%，−10.32% 和7.13%，其中WALE亚格子模型预报的脱落频率误差最小。 图11 不同亚格子模型计算的空泡脱落频率Fig. 11 Shedding frequency of cavity for various subgrid-scale models 图12 给出了一个周期内不同亚格子模型所预报的空化演变过程，同样采用汽相体积分数α=0.1等值面呈现典型时刻云状空泡形态，并与试验结果进行对比。 图12 不同亚格子模型模拟的空化演变过程与试验结果对比Fig. 12 Comparison of cavitation evolution by various subgridscale models and experimental results 由图可见，与试验结果相比，3 种亚格子模型基本都能够模拟出非定常云状空化演变过程的典型特征，但模拟得到的空泡形态差别较大，主要包括： 1） 片空泡闭合曲线轮廓形状：对比片空泡的弦向范围（如图12 中所示时刻9），可以发现WALE亚格子模型计算的结果与试验结果较吻合，片空泡范围最大时，弦长方向可延伸至弦长的50%左右，相应地，WMLES 和SL 亚格子模型计算的结果则偏小，只延伸至弦长的40%左右，但三者导边附近的片空泡展向延伸的长度相当。此外，由图13 给出的一个周期内不同亚格子模型计算的Lc/C值随时间的变化情况可知，WMLES 与SL 亚格子模型计算的Lc/C值相当，但明显低于试验结果[4]，而WALE 亚格子模型计算得到的Lc/C值与试验结果吻合较好。 图13 不同亚格子模型计算的Lc/C 值对比Fig. 13 Comparison of cavity length-chord ratios for various subgrid-scale models 2） 二次脱落：观察图12 中时刻3 至时刻7 的试验结果，可见片空泡“耳垂”位置受到回射流扰动而呈现出不稳定性的区域逐渐扩大，但WMLES与SL 亚格子模型计算结果在该区域上的空泡结构与试验相比要稀疏得多，特别是SL 亚格子模型，而WALE 亚格子模型计算结果在该区域观察到了更丰富的小尺度空泡结构。 图 14 YC=0.12 m 展向位置处水翼表面的分布Fig. 14 The time-averaged pressure distribution around hydrofoil surface at YC=0.12 m 3） 马蹄形云空泡：从俯视角度看，WALE 和WMLES 亚格子模型得到的马蹄形云空泡的尺度范围与试验结果吻合较好，而SL 亚格子模型计算得到的马蹄形云空泡明显小于试验结果。对于翼展中部马蹄形云空泡的发展而言，WMLES 亚格子模型在时刻1 能够模拟出上一个周期主脱落形成的马蹄形云空泡即将在翼面上通过随边的现象，而另两个亚格子模型得到的马蹄形云空泡则提前发生了溃灭。 表3 不同亚格子模型时均升阻力系数计算结果的对比Table 3 Numerical results of the time-averaged lift and drag coefficients for various subgrid-scale models 图15 YC=0.15 m 展向位置处水翼表面的分布Fig. 15 The time-averaged pressure distribution around hydrofoil surface at YC=0.15 m 本文基于LES 方法和S-S 空化模型，通过3 套网格密度及3 种亚格子模型，研究了Delft Twist11N 三维扭曲水翼的空化演变过程、空化脱落频率及升阻力系数等，得到如下结论： 1） 在最大壁面y+<1 的前提下，当弦长方向平均网格无量纲间距满足 ∆x+≤94 条件时，所用网格密度下的数值模拟既能够捕捉到较多的细小空泡脱落和马蹄形云空泡的初生及溃灭等非定常空化演变现象， 又能够获得具有一定精度的空泡脱落频率、时均升阻力系数和时均压力分布。 2） 3 种亚格子模型均能够模拟非定常云状空化演变过程中的典型特征。相较于WMLES 与SL 亚格子模型，WALE 亚格子模型捕捉到的片空泡及云空泡的演变与试验结果更吻合，且预报的空泡脱落频率、时均升阻力系数和时均压力系数精度更高。 因此，推荐选择基于WALE 亚格子模型的LES方法进行非定常云状空化现象的数值模拟。3.2 亚格子模型对空化数值模拟的影响
