激波传播与云空化脱落过程脉动冲击研究

邱 宁 朱 涵 周文杰 潘中永 袁寿其 刘 祥

(1.江苏大学流体机械工程技术研究中心， 镇江 212013； 2.江苏大学能源与动力工程学院， 镇江 212013)

0 引言

空化会产生冲击波，引起噪声、空蚀和振动[1]。空化流动具有复杂的非定常流动特性和高度的可压缩性，涉及流体力学、声学和材料学等多门学科[2-4]。目前，对于云空化的产生、脱落和溃灭等机理的研究在逐步深入[5]。

云空化被认为是附着在材料表面的不稳定空穴断裂脱落的结果[6-8]。文献[9]最先在实验中观察到了云空化脱落的周期性行为，同时发现了回射流的存在[9]。此后，云空化的形成和脱落机理受到各国学者的关注，大量学者使用理论分析、数值模拟和实验观测的方法进行了相关探索，主要是以水翼为载体进行了研究[10-12]。研究表明，云空化的脱落存在两种主要机制：回射流和激波。一些学者对激波进行了研究[13-16]，发现在σ/(2α)大于4时(σ表示空化数；α表示水翼攻角)，回射流是空化失稳的原因，当σ/(2α)较小时，激波主导了附着空穴的断裂；云状空化的溃灭会产生向周围传播的压力波；激波的传播造成了空化流中的材料损伤和噪声。

云空化演变的数值模拟，主要是基于Navier-Stokes(N-S)方程。由于激波的产生与不可压缩气体有关，因此，需要在模拟中考虑汽相的压缩性[17]。文献[18]首次利用可压缩流动方程，探究了空泡回弹过程中的激波形成过程。文献[19]提出了一种基于隐式亚网格尺度的大涡模拟方法，很好地捕捉到了激波的产生。文献[20]对有限质量传递模型进行修改，成功模拟了水翼流动的激波，并对其进行了动力学分析。文献[21]引入液相和汽相的状态方程进行数值模拟，捕捉到了空穴溃灭引发激波的现象。目前为止，对于NACA0015水翼多数研究只是对云空化的发展形态进行捕捉，较少使用压力传感器捕捉水翼表面压力变化和激波传播。

本文使用实验与模拟计算相结合的方法对NACA0015水翼的云空化流动和压力脉动进行研究。实验使用高速摄像机和压力传感器进行同步测量获得云空化的流场信息；利用二次开发技术，对湍流模型进行汽相的压缩性修正以捕捉其流场信息。通过实验和模拟相结合，分析相关流动机理。

1 实验装置和测试系统

实验在浙江大学航空航天学院空化水洞中进行[22]。实验段长度为1 000 mm，截面为200 mm×200 mm的正方形。实验段结构如图1所示。实验采用的NACA0015水翼弦向长度为100 mm，展向长度为200 mm。在弦长中心处设置旋转轴，将水翼固定于实验段并可绕轴调整所需要的攻角。实验时，可通过流道的透明窗观察空化形态。使用高速摄像机和压力传感器进行同步测量。实验段结构和测试系统如图1所示。

2 数值模拟方法

2.1 连续性方程和动量方程

空化流动存在汽相和液相的相互转换，在模拟中使用基于牛顿流体的N-S方程。该方程在笛卡尔坐标系下的表述为

(1)

(2)

(3)

ρm=ρlαl+ρvαv

(4)

μm=μlαl+μvαv

(5)

式中i、j——坐标方向

u——速度p——压力

ρl——液体密度ρv——蒸汽密度

αv——蒸汽体积分数

αl——液体体积分数

μl——液体层流粘度

μv——蒸汽层流粘度

μT——湍流粘度

m+——凝结速率m-——蒸发速率

ρm——汽-液混合相密度

μm——汽-液混合相层流粘度

t——时间

xi、xj——在i、j方向上的坐标

δij——切应力

2.2 空化模型

本文使用Zwart空化模型[23-24]进行计算。通过建立汽液两相的输运关系来描述空化相变过程。蒸发速率和凝结速率定义为

(6)

(7)

式中RB——气泡半径

αnuc——气核体积分数

pv——饱和蒸汽压力

Cprod——当地静压大于饱和蒸汽压力时的蒸汽凝结系数

Cdest——当地静压低于饱和蒸汽压力时的蒸汽蒸发系数

模型中各系数取值为：RB=1×10-6m，Cprod=0.01，Cdest=50，αnuc=0.01%。

引入湍动能对饱和蒸汽压力进行修正，修正后的饱和蒸汽压力计算公式为

ptur=0.39ρmk

(8)

pv=psat+0.5ptur

(9)

式中k——湍动能

ptur——局部湍流压力波动值

psat——未引入湍动能时的饱和蒸汽压力

2.3 湍流模型

SST(剪切应力传输)湍流模型是一种基于标准k-ε模型和标准k-ω模型的混合模型。该模型能够较为准确地模拟近壁面流动特征和流体流动过程中的逆压梯度。但该模型存在对湍流粘度的过度预测，同时无法准确模拟分离区域内湍流的输运特性。为解决以上问题，对SST模型进行合理修正。使用FBM(滤波器)模型[25]对远离壁面的大尺度空化旋涡区域进行滤波修正，捕捉大尺度涡流的流动特性，对近壁面使用DCM(密度修正)模型[26]进行混合密度修正，以考虑汽相和液相混合的压缩性。湍流粘度定义为

(10)

其中

fhybrid=ξ(ρm/ρl)fFBM+
(1-ξ(ρm/ρl))fDCM

(11)

(12)

(13)

(14)

式中，λ为滤波尺度，参数C1、C2和Cμ的推荐值分别为：4、0.2和0.09[27]。使用桥接函数ξ(ρm/ρl)将DCM模型和FBM模型结合在一起，对不同的区域使用不同的方式进行求解以捕捉不同尺度的流动特性。

2.4 数值设置与网格验证

为了保证数值计算的精度并提高计算效率，使用ANSYS ICEM软件对该计算模型进行六面体结构化网格划分。水翼周围采用O-block型拓扑结构，并在水翼近壁面区域进行网格加密，从而较好地控制水翼表面边界层流动处网格结构，网格细节如图4所示。由于SST模型对网格具有较高的要求，在处理近壁面网格时要使y+(无量纲化的壁面距离)分布尽可能小于1，壁面法向网格尺度变化比和沿流向的网格尺度变化比不大于1.2[30]。

考虑网格对计算的影响，采用初始的SST模型对水翼计算域进行了网格无关性验证。表1给出了NACA0015水翼流场的网格信息。

表1 NACA0015水翼流场网格信息Tab.1 NACA0015 hydrofoil’s mesh information

以一组网格总数为894 780的网格为基础，对水翼周围O-block型拓扑结构进行加密，在O-block结构内的x、y、z方向上每次各增加10个节点，共设置5组网格。对不同网格数下的流场进行水的单相定常流动计算，并对升阻力系数进行对比。升力系数和阻力系数的定义为

(15)

(16)

其中

A=CS

式中Fl——水翼升力

Fd——水翼阻力

A——水翼等效面积

S——水翼展长

不同网格数下翼型的升阻力系数变化趋势如图5所示，随着网格数的增大，升力系数Cl和阻力系数Cd的差异越来越小，呈现出明显的网格无关性趋势。考虑计算精度和资源，最终选用网格方案Ⅴ(网格总数为4 355 420)进行后续计算。该方案下水翼的y+在0.046 9～7.236之间，其中大部分在1以下，表明该网格方案下的近壁面流动情况符合湍流模型的计算要求。水翼壁面的y+分布如图6所示。

3 结果与讨论

3.1 云空化形态的演变

在计算过程中，提取了整个流体域中的空泡体积随时间的变化，如图7所示。随着时间的发展，空泡体积展现出规则的周期性变化。空泡体积的最小值稳定在1.5×103mm3左右，并不为零，说明在流体域中，始终有空泡的存在。这与实验中观察到的现象一致：空化云在水翼尾部溃灭的同时，头部有附着空穴产生，两部分空化形态同时存在，在时间上没有断层。对空泡体积进行频域分析，得到该工况下空化周期性变化的主频为27.78 Hz。实验测得的空化周期主频为29.98 Hz，较为接近，表明模拟方法具有较高的可信度。

图8为一个周期内高速摄影和模拟计算的空化形态俯视图，每幅图像间隔为1/8周期。模拟计算的空化云图取空泡体积分数为10%的等值面。在t0～t2时刻，水翼前缘有附着空穴形成和发展。同时可以观察到由于水翼两侧壁面的存在，影响了水翼展向端点附近的剪切力和湍动能，导致水翼展向两个端点附近的附着空穴生长较快。激波在吸力面上的传播造成了附着空穴的断裂脱落。在t3时刻，水翼的后缘出现了较为剧烈的波动，此时激波在该区域传递，附着型空化开始被激波带来的水流抬升离开壁面。在t4时刻，通过高速摄影拍摄到的图像可以明显观察到：在Ⅰ区域，附着空穴形态较为稳定，此时激波还未传播到该区域。Ⅱ区域附着空穴形态出现轻微波动，表明此时激波传递到此处开始改变该区域的流动状态。在t5时刻，激波将附着空穴截断，水翼前缘空穴回缩，此时，水翼吸力面大部分区域被空化云覆盖。在t6～t7时刻，脱落的空化云向下游移动并聚集。

图9为实验和模拟结果的对比图。其中图9a为经过灰度化处理后的高速摄像图片，图9b为捕捉到的仿真结果。对应时间的激波传播趋势基本一致。激波传递的前缘与附着空穴存在明确界限，这一发现与文献[31]的实验结果一致。在t=0 s时刻，吸力面的大部分区域被空穴覆盖，水翼远场溃灭的大尺度云空化产生的激波传递到水翼尾缘，高强度冲击压力在接触水翼吸力面后向前缘传递，使得空穴尾部与壁面分离，从图9b中速度矢量可以清晰看到这一现象。在t=0.003 2 s激波传播至水翼吸力面中部，附着型空穴已经在激波作用下发生断裂，从速度矢量可以看出，云空化周围的水流呈顺时针旋转趋势，近壁面处水流继续沿壁面向水翼上游移动。在t=0.004 8 s和t=0.005 6 s时刻，激波带来的逆向水流继续向前缘移动，导致水翼前缘的附着型空穴回缩，长度逐渐变短。随着云空化向下游移动，下游吸力面重新被空穴覆盖，成为低压区。由于此时水翼中游无空化覆盖，导致此处压力较高。当空化发展至t=0.007 2 s时刻，激波传递至水翼前缘，此时与来流相遇，在水翼前缘形成一个小范围高压区。

3.2 云空泡溃灭的动力学分析

大尺度空化云溃灭到最小体积后由于不可压缩气体的存在而回弹，产生激波向外传播[18]。为了更直观地分析激波的产生过程，对空化云溃灭过程的压力分布情况进行研究，以水翼展向中线为参考，设置压力面，探究大尺度云空化溃灭产生的压力变化，为更好描述，将该压力面命名为A。

图10为云空化的发展趋势和A面及水翼表面的压力分布。在T1时刻，激波刚击穿附着在水翼表面的空穴，使其断裂，在水翼吸力面上形成分散的涡状空化，此时水翼上表面的空化区域距壁面较近，水翼上表面附近大部分区域压力较低，空化未覆盖区域压力有所升高。在主流推动下，云状空化向水翼尾缘移动。在T1+9.6 ms时刻，脱落的空化云在水翼尾缘聚集成大尺度的马蹄形空穴结构，该区域压力较低，由于在此过程中伴随着部分空泡的溃灭，周围的压力高于该区域。同时由于吸力面中段不再有空泡覆盖，可以观察到该区域压力回升。在T1+13.2 ms时刻，大尺度空化云向下游输运过程中，移动较慢的云空化结构首先发生断裂溃灭，如图10中a1区域所示。而后，在T1+18 ms时刻，云空化中的旋涡空化结构发生断裂破坏，如图中的b1区域所示。在空穴完全溃灭的瞬间，即图中的T1+19.2 ms时刻，观察到以溃灭点为中心的高压区域，此时激波产生，溃灭点压力达650 kPa。压力波在抵达翼型表面后，从水翼后缘向前缘迅速传播。

空化的发展和水翼周围的流场结构有着密切的关系。为了更好地分析流场结构，引入Q判据来描述[32]。

Q是涡度张量和形变率张量的差值。当Q>0时，意味着旋转的趋势大于轴向形变的趋势，此时可认为该区域的流动是由涡流主导的。

图11为一个典型时刻T3空泡体积分数为10%的等值面上的Q分布情况。从图中可以看出，在该时刻，水翼前缘被附着空穴覆盖，吸力面中下游存在上个周期脱落的云空化。在附着空穴前缘，存在一个细长的区域，具有较大的Q值，Q值沿水翼弦向向后逐渐降低。该区域具有较大的Q值可能是由于来流在接触到水翼前缘时方向突然改变，使得此处流动极不稳定。在大尺度空化云结构中，涡量最大的部分出现在管状涡空化结构中，该结构具有明显的旋转特征，旋转效应远大于形变效应，导致该区域流动出现强烈的不稳定性。

图12为一个周期内的涡分布情况和对应的空泡体积分数分布。在T1～T5时刻，上一周期脱落的空化云在吸力面跟随主流向下运动，同时在水翼前缘有附着空穴生长，在这个过程中，水翼前缘Q等值面的分布较为连续，随着时间的推移，该区域覆盖的吸力面范围也呈增大趋势，与空穴的生长过程基本一致。这说明在附着空穴的生长过程中，空化区域边缘的旋涡流动趋势较为连续和稳定；在空化云移动和发展区域，由于空泡在此区域不断溃灭，该部分的流动具有剧烈波动。在T6～T8时刻，激波开始作用于水翼表面，空穴闭合区的结构变得不再稳定，开始出现断裂，随着时间的推移，断裂区域的范围开始向前缘逐渐扩大。

图12给出了A面上的Q值分布情况，从图中可以看出，Q值的分布规律与空穴长度基本对应，且Q值较高的区域分布在水汽交界面附近。在T1～T5时刻Q值较高的区域主要位于上一周期脱落的云空化外围，由于部分空泡溃灭，导致此处的流动不够稳定，同时随着附着空穴的生长，空穴的相分界线附近也有较大的Q值出现，而近壁面处Q值较低，说明此时水翼表面上的流动相对稳定。在T6～T8时刻，激波传递到水翼尾缘，推动近壁面流体向水翼前缘运动。尾缘区域出现相反方向的涡结构，Q值最大的区域和Q值最小的区域交叉存在，旋转效应和形变效应相互制约。随着激波的传播，带有正向涡量的流体(Q>0)向水翼前缘移动，Q值较大的区域从水翼尾缘逐步向上游扩大，同时在水翼近壁面处也出现了较大的Q值，在空化发生的核心区域，激波将附着空穴切断。

3.3 压力脉动分析

受到逆压梯度的影响，在空穴生长至一定长度后，会在空穴尾部产生回射流，贴近壁面向上游移动。图13给出了空穴发展过程中回射流的分布，对应于图14中的a、b和c时刻。可以看出，在空穴生长过程中，回射流一直存在。随着附着空穴的发展，回射流覆盖区域逐渐扩大，但回射流前缘的位置并未出现继续向水翼头部推进和切断附着空穴的趋势。同时在图14对应时刻，相应位置的压力稳定保持在较低水平，表明该工况下的空化演化并不是由回射流主导的。

水翼吸力面空化结构的演变与表面的压力脉动密切相关。图14为其中5个周期的监测点#1、#2和#3的压力脉动信号，从图中可以看出，数值模拟方法很好地捕捉到由激波造成的压力脉动趋势，当空化增长到覆盖监测点时，压力最低；当空化不再覆盖监测点处壁面时，压力开始上升。由于监测点设置在不同的弦向位置，所测结果也对应于不同的空化发展阶段。实验得到的监测点#1、#2、#3处的压力脉动幅值约为61.74、69.14、37.06 kPa。计算得到的监测点#1、#2、#3处的压力脉动幅值约为152.65、122.27、111.93 kPa。数值计算和实验结果得到的幅值存在差异的原因可能在于：模拟中监测的为点的静压，实验中的压力传感器记录的是以所测点中心为圆心的面平均静压。由于测量面的面积远大于空泡投影到水翼表面的面积，导致无法精确捕捉到空泡的溃灭压力，使得脉动偏小。

图中相邻两条绿线代表一个完整的空化周期，绿色箭头表示附着空穴的发展趋势。空化在水翼头部产生，沿弦长向尾缘生长。由于空穴的生长需要时间，导致靠后的监测点更晚被附着空穴覆盖。当空穴覆盖监测点时，监测点测量到的压力最小，对应于该温度下的饱和压力。空化覆盖到吸力面的大部分区域时，逆压梯度的作用会导致空穴闭合区的流动失稳和空泡卷起。上一个周期脱落的大尺度空化云溃灭后形成的激波传递至水翼表面时，开始在尾缘出现向前缘的流动，将空穴与壁面分离，传递趋势如图中蓝色箭头所示。首先传递至监测点#3，导致该点处压力上升，而后依次向前传递。激波传递的速度可通过相邻监测点开始出现大幅度压升的时间差来计算(如图中蓝线和黑线的交点)。NACA0015水翼在攻角12°、空化数1.4下，实验得到激波在监测点#2和#1之间传播的平均速度为11.31 m/s，模拟得到激波在监测点#2和#1之间传播的平均传递速度约为11.53 m/s，二者在数值上较为接近。

4 结论

(1)基于FBM-DCM方法对SST湍流模型进行可压缩修正后，很好地捕捉到了激波的产生过程。大尺度空泡云在主流作用下向下游运输，在水翼尾缘远场溃灭至最小体积，释放出高强度压力脉冲。

(2)Q值较大的区域主要位于水汽交界面附近。附着空穴覆盖的水翼表面流动较为稳定，云空化区域空化流动的波动较大。

(3)空泡覆盖区域，压力处于较低水平，激波在吸力面上从后缘向前缘传播，传递过程中会导致当地压力的回升。该工况下，实验得到激波在监测点#1和#2之间速度为11.31 m/s，模拟得到的激波在监测点#1和#2之间速度为11.53 m/s，较为接近。