李大尉，赵孟石，曾卿丰，郑智颖，裴 禹，姚顺宇，姚立明

(1.黑龙江省科学院 高技术研究院，黑龙江 哈尔滨 150020；2.哈尔滨工业大学 能源科学与工程学院，黑龙江 哈尔滨 150001)

0 引言

当物体以足够高的速度在水下运动时，根据伯努利方程物体周围局部压力下降，在此作用下水中气核增大，破坏水的连续性从而发生空化。空化产生的空泡附着在物体上，且受物体运动速度增大的影响而得到充分发展后，空泡融合并扩大成为一个整体，最终形成尺寸超过物体特征尺寸的超空泡。超空泡的存在可减少水下航行体与水的接触面积，从而显著减小航行阻力，极大地提高水下航行速度[1]，减阻率最高可达90%左右[2]，因此可以通过形成超空泡包裹水下航行体的方式，给鱼雷或潜射弹体进行减阻[3]。此外，超空泡表面高效的相变特性可以应用在海水淡化领域，对超空泡内产生的蒸汽进行收集和冷凝即可得到淡水[4-5]。空泡溃灭过程中产生的高温高压还可以应用于溶液的消毒杀菌[6]。

实验中自然超空泡的产生一般通过高速射弹和高速来流冲击两种方法，但前者产生的超空泡不够稳定且难以测量，后者则需要庞大的水洞试验系统[7]。而通过旋转机械来产生超空泡有系统体积小、空泡稳定以及容易控制和利用等独特优势[7]。目前对于采用旋转空化器产生旋转超空泡的研究相对较少，在前期研究[8-9]中已针对四叶片楔形叶片旋转空化器开展了数值仿真研究，并研究了减速板对于旋转空化器中的作用。本文设计了一种双叶片空化器，对其三种不同叶型在不同转速下的旋转空化流动进行了三维定常数值模拟计算，对比分析了其水动力学特性，总结了转速和叶型的影响规律，所得结果对后续楔形叶片旋转空化器的叶型设计具有指导意义。

1 几何模型和网格划分

本文针对双叶片旋转空化器建立了三种不同叶型，其示意图如图1所示，两个叶片呈中心对称。其中，叶型1为原型楔形叶片，其横截面为顶角45°的等腰三角形，如图1(a)所示，已有研究表明顶角45°为获得超空泡的最佳角度[3-4]。叶型2为改型叶片，其沿平行于原型楔形叶片进口边的方向，在进口边上下两侧对称地加工了两个沟槽，从而形成一个主进口边和两侧的两个副进口边，如图1(b)所示。叶型3也为改型叶片，在保持原型楔形叶片叶尖端形态不变的基础上，将叶片其他部分的进口边加工为向叶片内部凹陷的沟槽，从而形成两个平行的进口边，如图1(c)所示。图2以叶型1(即原型楔形叶片)为例给出了数值模拟中建立的计算模型，空化器直径为90 mm，转轴直径和高度分别为28 mm和20 mm。选取高25 mm和直径100 mm的圆柱体区域作为计算域。计算域的上、下表面均设置为压力值为101 325 Pa的压力出口边界条件，圆柱侧面设置为固定(转速为0)无滑移壁面边界条件，其他边界设置为转速相同的无滑移壁面边界条件，转速值根据不同工况而改变。

图1 旋转空化器叶型示意图

图2 计算模型示意图(以叶型1为例)

图3以叶型1为例给出了计算域网格划分的示意图，计算域的划分使用非结构化网格，对空化器附近的流域使用了更高密度的网格进行划分，以满足数值计算中y+值的要求。为了保证计算的准确性和经济性，针对原型楔形叶片旋转空化器在转速ω=6 000 r/min下的工况进行了网格无关性验证。图4给出了数值模拟得到的超空泡体积Vc和表面积Ac随网格数的变化，从图中可以看出，当网格数超过800万后，Vc和Ac不再随网格数变化，因此在后续数值计算中针对叶型1选取了数量为9 690 752的网格。采用相同网格尺寸和划分方法对叶型2和3的计算域进行网格划分，其相应的网格数量分别为10 650 748和8 246 471。需要注意的是，在计算超空泡的体积和表面积时空泡轮廓由气相体积分数αv=0.1的等值面确定。

图3 计算域网格划分示意图(以叶型1为例)

图4 网格无关性验证

2 数值方法

数值计算过程中，采用PRESTO!格式离散压力方程，采用SIMPLEC算法求解速度与压力的耦合，湍流模型采用Realizablek-ε模型，并采用Scalable壁面函数来处理近壁面流动。采用基于无滑移假设的Mixture多相流模型来描述空化两相流，该模型假设气液界面处流体速度连续，气、液两相之间不存在速度滑移，视为混合相进行处理，相应控制方程为：

混合相连续性方程

(1)

混合相动量方程

(2)

气相体积分数方程

(3)

式中p——压力;

ui和uj——混合相速度;

μt——湍流有效粘度;

αv——气相(水蒸气)体积分数；

ρv和ρm——气相和混合相的密度;

μv和μm——气相和混合相的动力粘度;

Rc和Re——气相的冷凝率和产生率。

混合相的密度ρm和动力粘度μm由下式定义

ρm=αvρv+(1-αv)ρl

(4)

μm=αvμv+(1-αv)μl

(5)

式中μl和ρl——液相(水)的动力粘度和密度。

气相体积分数方程中的气相产生率Re和冷凝率Rc需要空化模型来描述，本文中采用了Schnerr-Sauer模型：

当p≤pv时

(6)

当p≥pv时

(7)

式中 pv——饱和蒸汽压；

R——气泡半径，其与αv之间存在如下关系

(8)

式中nb——气泡的数量密度，计算中其值为nb=1×1013m-3。

数值计算中水和水蒸气的物性根据IAPWS(国际水和蒸汽性质协会)数据库进行设置。

3 计算结果及分析

为了探究不同叶型旋转空化器的水动力学特性，对不同转速(3 500 r/min、4 000 r/min、5 000 r/min、6 000 r/min、8 000 r/min、10 000 r/min、12 000 r/min)工况进行了数值模拟计算，得到了空泡形态随叶型和转速的变化。

图5给出了不同转速下三种不同叶型空化器形成的空泡形态的对比。从图中可以看出，相同转速下三种叶型空化器叶片后方形成的空泡形态相似，空泡尺寸随转速的升高而增大。较高转速下叶尖处空泡尾部光滑，而半径较小处空泡尾部则向空泡内部凹陷，这是由于在空泡的尾部形成了回射流[4]，且凹陷区域随着转速的升高而逐渐缩小。叶型1产生的空泡形成于叶片出口边，叶型2和3产生的空泡还可以形成于两个副进口边(叶型2)和上下两个进口边(叶型3)。转速低于4 000 r/min时，空泡无法完整地包裹叶片出口边，叶尖部分出口边无空泡形成，而叶根处产生的回射流引起了较大面积的空泡脱落(图5(a)、(e)和(i))。转速低于5 000 r/min时，叶型2副进口边和出口边形成的两部分空泡分割开来，如图5(e)和(f)所示；转速高于6 000 r/min时，两部分空泡则能够连接成一个整体，如图5(g)和(h)所示。而对于叶型3而言，转速低于12 000 r/min时，空泡虽然能在上下两个进口边处形成，但无法与叶片出口边处形成的空泡连通，如图5(j)和(k)所示；只有转速达到12 000 r/min时，两部分空泡才连通，如图5(l)所示。此外，叶尖端面处也有空泡的产生，转速较高时叶尖端面处的空泡与出口边处的超空泡会连接起来。随着转速提高(高于8 000 r/min)，叶型2产生的空泡持续发展并与计算域边界(包括圆柱面固壁边界和上下压力出口边界)接触，空泡边界未被捕捉到而产生空缺，因此空泡实际表面积比计算得到的表面积更大，如图5(g)和(h)所示。值得注意的是，较高转速下叶型2产生的空泡与上下压力出口边界的接触表明，需增大计算域的高度以满足较高转速下对空泡形态的捕捉；而空泡与圆柱面固壁边界的接触表明，在当前旋转空化器装置的设计尺寸下，空泡将与四周固体壁面接触，进而在一定程度上影响空化器的工作性能，这也为后续旋转空化器装置的设计提供指导。

图5 不同转速下不同叶型旋转空化器形成的空泡形态

图6针对三种不同叶型给出了衡量旋转空化器性能的三个特性参数(空泡表面积Ac和体积Vc以及空化器阻力矩Td)随转速的变化规律。从图中可知，空泡体积Vc随转速的升高而增大。当转速低于5 000 r/min时，叶型1和叶型2形成的空泡体积几乎相等，当转速为5 000～10 000 r/min时，叶型2产生的超空泡体积明显大于叶型1产生的超空泡体积，而叶型3形成的空泡体积则始终小于相同转速下叶型1和2形成的空泡体积。此外，转速低于12 000 r/min时叶型2产生的空泡表面积明显大于另外两个叶型，叶型1和叶型3形成的空泡表面积的差异不大；而转速为12 000 r/min时，叶型2产生的空泡表面积略小于叶型1和3产生的空泡表面积，这是由于叶型2产生的空泡与计算域边界接触导致计算得到空泡的表面积比实际表面积要小，但不能说明真实情况下转速为12 000 r/min时叶型2产生的空泡表面积较小。空泡尺寸的增大也使得叶型2的阻力矩增大，但是在叶型1产生的空泡体积要大于叶型3产生的空泡体积的情况下，叶型3的阻力矩却始终大于叶型1的阻力矩。综上，从空泡尺寸的角度来说，采用叶型2的旋转空化器的性能优于采用另外两种叶型的旋转空化器。

为了进一步分析三种叶型空化器阻力矩的构成，图7给出了不同叶型旋转空化器在不同转速下由压力和粘性摩擦力产生的阻力矩。从图中可以看出，对于不同叶型的空化器，压力产生的阻力矩Tp起主要作用，而粘性摩擦力产生的阻力矩Tv对总阻力矩的贡献较小。三种叶型的Tp由大到小的顺序为：叶型2>叶型3>叶型1，因此三种叶型之间阻力矩Td的变化呈现相同的趋势(图6)。但三种叶型之间Tv的变化趋势却与Tp的变化趋势完全相反，这是由于叶型2和3的叶片表面有更大的面积被空泡覆盖(图5)。此外，随着转速的增加，三种叶型阻力矩之间的差异呈现变大的趋势。

图6 不同转速下不同叶型的特性参数

图7 不同叶型旋转空化器受到的由压力和粘性摩擦力产生的阻力矩随转速的变化

一般来说，空化器所受阻力随空泡尺寸的增大而增大，但图6的结果表明，叶型3形成的空泡尺寸小于叶型1形成的空泡，但采用叶型3的旋转空化器受到的阻力矩反而更大，为了分析其成因，图8给出了转速为6 000 r/min时三种不同叶型旋转空化器叶片表面的压力分布。从图中可以看出，相比于叶型1，叶型3在较大半径处的进口边区域的压力更大，叶型1、2和3表面的最大压力分别为400 147 Pa、441 767 Pa和470 497 Pa。进一步分析图9所示的不同叶型空化器叶片横截面压力分布可知，对于叶型1，流体在冲击进口边后分成两股，分别沿进口边上下两侧流向下游，因而在进口边处存在一个高压区；而对于叶型3，其叶尖部分的形状与叶型1相同，因此该处压力分布也与叶型1的压力分布相似，但在除叶尖外的其他叶片部分，流体冲击上下两个进口边之间的沟槽，因而在沟槽内形成高压区域，且压力大于叶型1进口边处的压力，从而造成叶型3受到的阻力矩大于叶型1受到的阻力矩。此外，叶型2的主进口边和两个副进口边之间的沟槽中同样由于受到流体冲击而存在高压区，且高压区的强度和大小均大于叶型1进口边处的压力分布。

图8 转速为6 000 r/min时不同叶型旋转空化器叶片表面压力分布

图9 转速为6 000 r/min时不同叶型旋转空化器叶片横截面压力分布(黑线表示空泡轮廓)

4 结论

为了探究叶型改进对楔形叶片旋转空化器水动力学特性的影响，本文针对三种不同叶型旋转空化器在不同转速下的旋转空化流动开展了数值模拟计算，对其产生的空泡形态特性及其阻力特性进行了比较和分析，主要得到了以下结论：

(1)叶型1(原型叶片)产生的空泡形成于叶片出口边，而叶型2和3(改进叶型)产生的空泡除了形成于出口边外，还可以形成于两个副进口边(叶型2)和上下两个进口边(叶型3)；对于叶型2，转速较低时副进口边和出口边形成的两部分空泡分割开来，转速较高时两部分空泡连接成一个整体；对于叶型3，在上下两个进口边处形成的空泡无法与出口边处形成的空泡连接成一个整体。

(2)三种叶型产生的空泡尺寸均随转速的升高而增大；相同转速下，叶型2产生的空泡尺寸最大，叶型1产生的空泡尺寸次之，叶型3产生的空泡尺寸最小；空泡尺寸的增大使得叶型2的阻力矩增大，但叶型3的阻力矩始终大于叶型1的阻力矩；压差阻力矩占总阻力矩的大部分，而粘性摩擦阻力矩的贡献较小，受叶片表面被空泡覆盖面积不同的影响，三种叶型之间粘性摩擦阻力矩的变化趋势与压差阻力矩的变化趋势完全相反。

(3)通过分析压力分布发现，流体冲击叶型3上下两个进口边之间的沟槽，在沟槽内形成高压区域，从而造成叶型3受到的阻力矩大于叶型1受到的阻力矩，因此采用叶型3的叶片改型产生了一定的负面效果，采用叶型2的叶片改型更为合理。