考虑挡板对空化脱落的抑制效应数值分析

2020-12-16万德成

海洋工程 2020年6期

李勇，刘成，万德成

(上海交通大学船海计算水动力学研究中心 (CMHL) 船舶海洋与建筑工程学院海洋工程国家重点实验室，上海 200240)

空化流动是一种复杂的气液两相流动，它涉及到流体力学中诸多的流动问题，包括相变、质量传输、可压缩性和非定常性。当液体中的局部压力低于当前温度下的饱和蒸汽压力时，空化现象便发生了，所以空化常常发生在速度和压力梯度大的水利设备中，例如泵和螺旋桨[1]。空化流动是非定常的，因为空化现象伴随着空泡周期性的产生与溃灭，这会给水利设备造成不利的压力震荡与疲劳破坏。不管是试验手段还是理论手段，都很难直接研究空化现象的机理机制，现在对空化现象的研究都是从结构简单的水翼着手，逐步往结构复杂的螺旋桨或者泵发展。

根据空化不同的形态和形式，空化可以被分为片空化，云空化，超空化和涡空化等。片空化本身就是非定常的，在生成并发展到一定的程度后就会坍塌破碎成一束束细碎的空泡，并从附着在水利设备上的片空化脱离开来，形成云空化[2]。大量的观察和研究被用来分析片空化发展成云空化的机理机制，最先被发现的片空化脱落机制是回射流机制。Knapp等[3]发现在空化流动中存在与主流方向相反的回射流。Furness等[4]对一个二维的缩放喷管进行理论和试验研究，发现在片空化发展的过程中，一股与来流方向相反的回射流，从片空化末端的底部，紧靠着管壁向上游移动，并在喷管喉部将片空化切断。Kawanami等[5]通过高速摄像技术和压力检测技术对云空化进行一系列的试验观察，发现片空化的脱落和云空化的产生是由于片空化的底部存在着由尾部向导边移动的回射流。同年，Kawanami等为了验证回射流在空化流动中的作用，在水翼上表面的不同位置布置了长短不同的挡板，来观察挡板的存在对回射流以及片空化脱落效果的影响，研究发现，片空化底下的挡板能阻止回射流向上游移动，从而一定程度上缓解片空化的脱落。以上研究都定性地说明回射流的存在会使片空化脱落并产生云空化。Callenaere等[6]通过观察回射流相对于片空化厚度的大小，定量地分析了回射流的相对厚度对片空化稳定性的影响。

近年来，片空化断裂和脱落的另外一种机制被发现，也就是大尺度空化瞬间溃灭的激波机制。Arndt等[7]通过综合试验和数值模拟的手段，对一个二维NACA0015水翼的空化发展过程进行研究，研究表明，根据σ/(2α)取值的大小(式中的σ为空化数，α为水翼的攻角)，可以把片空化过渡到云空化的机制分为两类，当σ/(2α)的值比较大时，回射流机制是片空化发展成云空化的主要原因；当σ/(2α)的值比较小时，激波机制是片空化脱落生成云空化的主要原因，同时Arndt等[7]指出σ/(2α)=4是取值大小的过渡阶段，σ/(2α)≈4时两种控制机制都会对片空化的脱落产生影响。也就是说，水翼在相同的攻角下，空化数σ越小，激波机制越明显；空化数σ越大，回射流机制则成为空化脱落的主要原因。Ganesh等[8]通过高速摄像技术和X射线来观察片空化向云空化的过渡过程，通过系统地改变空化数σ的大小，在楔形顶点的流动分离区观察空化的水动力性能，发现在空化数较小时，激波机制成为了片空化脱落的主要因素。

对空化的研究一般采用试验与数值模拟相结合的方式。在试验中激光多谱乐测速(LDV)技术和粒子图像测速(PIV)技术是观察空化最常用的方法。Kubota等[9]通过高速摄像机和LDV技术观察到云空化一般处在最大涡强中间，并且发现云空化是由很多细小的空泡构成的，而片空泡则要纯净得多，几乎是块状的单质空泡。Foeth等[10-11]则通过时间分辨粒子图像测速(TR-PIV)技术观察回射流导致片空化的脱落和云空化的产生。现在采用数值模拟对空化现象进行研究的方法也流行起来，新型湍流模型和空化模型的更新和选取，让数值模拟的结果与试验结果更好的吻合。Hsiao等[12]提出了一种基于欧拉-拉格朗日耦合方法的多尺度空化模型，通过小尺度模型捕捉微观空泡结构，大尺度模型追踪大空泡团的动力学特性，从而清晰地描述小空化核生长成为宏观空泡的过程。

针对回射流对片空化脱落的作用，研究阻断回射流对片空化发展的影响。为此，分别对平板和水翼的吸力面添加挡板，用来阻挡逆流而上的回射流，从而研究不同布置不同长度的挡板对回射流的影响，以及比较分析添加挡板后的水动力性能。数值模拟采用开源CFD平台OpenFOAM的interPhaseChangeFoam求解器，选取Schnerr-Sauer空泡模型和修正后的SST k-ω湍流模型对带挡板的水翼进行数值求解，并将模拟结果与试验结果相互比较。

1 数学模型和研究方法

1.1 控制方程

在给出控制方程之前，先定义一个体积分数α，在网格中代表液体的相分数。在控制方程中，把气液两相看作是单质混合模型，共享速度和压力，而速度场需要相分数α来确定。考虑到气液两相系统需要在同一个网格里面表示，网格内全是液体时，则α=1；网格内全是气体时α=0；当相分数α介于0和1之间时，网格内是气液混合物，同时也看作是水气界面。由于在控制方程中只能用到一个体积分数，因此把液相的体积分数αl写作α，而气相的体积分数αv由αv=1-α得出。求解空化现象的控制方程由三部分组成：液体的相方程，连续性方程和动量方程。

(1)

(2)

(3)

U=αlUl+αvUv

(4)

p_rgh=p-ρgh

(5)

1.2 空化模型

采用的空化模型是基于输运方程模型(TEM)的Schnerr-Sauer空化模型[13]，该模型将水气混合物看做是包含大量球形蒸汽泡的混合物。根据Rayleigh-Plesset方程得到的不同源项用来模拟空化的汽化与液化过程，当局部压力低于饱和蒸汽压力pv时，发生剧烈汽化；当局部压力高于饱和蒸汽压力pv时，空泡液化为液态。空化模型可以由式(6)表示：

(6)

其中，Cv和Cc是经验常数，一般取1；ρl，ρv和ρ分别是液体、气体和气液混合物的密度；p是局部的压强，pv是液体在当前温度下的饱和蒸汽压力；Rb是气核半径，可由式(7)求出，式中的n是液相中气泡的数目。

(7)

1.3 修正的湍流模型

Larsen和Fuhrman[14]指出，常用的标准二方程闭合模型(k-ω和k-ε)都是无条件不稳定的，基于标准k-ω湍流模型变形来的SST k-ω湍流模型也有着同样的问题，即无条件不稳定，且有着渐近的不稳定增长率，使得部分区域的黏度总是大得不符合物理现象。为了避免过度预测湍流闭合区域的湍动能和涡黏度，采用一种新的稳定化的SST k-ω湍流模型来闭合控制方程。在新的SST k-ω湍流模型中，Larsen和Fuhrman推荐采用一种简单的办法来将湍流模型稳定化处理：通过对湍流黏性系数进行简单的修正，从而达到稳定湍流模型的目的。具体的思路就是将涡黏度νT由式(8)重新定义为式(9)。

(8)

(9)

在修正后的涡黏度公式中，分母中的max函数，除了原来就存在的两个判断大小的依据，新添加了第三个判断依据，在接近势流的流动条件下才会被激活。式中的a1，λ2和β*都是经验常数，其值分别取0.31，0.44和0.09，p0为湍流产生量，k为湍流动能密度，ω为湍流耗散率，p0为湍流产生量，pΩ为修正后的湍流产生量，F2是与k和ω有关的正切函数。根据解析分析，并通过这样的修正，湍动能和涡黏度不再随着计算的进行增长，而是渐进指数衰减[13]。

1.4 计算模型

为了研究挡板对回射流的影响以及回射流与片空化脱落之间的关系，采用了两个物理模型来进行分析。第一个是根据Zhang等[15]的试验设计一个平板水翼，其形状如图1所示，是一块尺寸为150 mm×200 mm×12 mm(长×宽×高)的平板水翼，其中150 mm是平板水翼的弦长，攻角固定为0°，一条全展长的挡板设置在平板水翼的上表面，挡板高2 mm，宽2 mm，距离导边0.37倍弦长。

在这个工况下，来流速度为10 m/s，饱和蒸汽压力取3 540 Pa。出口压力根据空化数求得，空化数的定义由式(10)给出，其中P取出口压强，Pv是饱和蒸汽压。采取的空化数和文献[14]中采用的一致，取值0.7。在整个计算域中，进口距离平板水翼的导边为2倍弦长，出口距离平板水翼的随边为4倍弦长，计算域的高度为2倍弦长，展长选取0.3倍弦长。具体的边界条件由图1给出。

图1 平板水翼的尺寸及其计算域Fig. 1 Sizes of the flat hydrofoil and its computation domain

(10)

由OpenFOAM自带的blockMesh工具来画整个计算域的背景网格，然后用snappyHexMesh工具，对平板水翼附近的一个矩形的加密区进行网格细化。分别采用10×105，15×105，18×105和20×105的网格进行计算，18×105网格的计算结果和20×105网格的计算结果已经比较相近，为了保证水翼附近的网格足够细，最终选取了20×105的网格来进行数值计算，平板水翼中剖面处的网格分布如图2所示。

图2 中剖面处的网格分布Fig. 2 Mesh distribution in mid-plane

第二个物理模型采用的是Kawanami[5]所做的一系列试验中用到的对称型椭圆鼻水翼(elliptic nose foil)，椭圆鼻水翼的形状和理论坐标由图3给出。椭圆鼻水翼弦长150 mm，并在一个高和宽为600 mm×150 mm的螺旋桨空泡水筒中进行试验。

图3 椭圆鼻水翼的理论坐标系统Fig. 3 E.N. foil and coordinate system

在第二种工况下，计算域的尺寸与进行试验的空泡水筒尺寸一致。计算域的入口距离椭圆鼻水翼的导边为2倍弦长，计算域的出口距离椭圆鼻水翼的导边6倍弦长，计算域的宽度和高度与空泡水筒的宽度和高度一致，也就是1倍弦长和4倍弦长。计算域的布局和边界条件由图4给出。

图4 椭圆鼻水翼的计算域和边界条件Fig. 4 Computational domain and boundary conditions of the E.N. hydrofoil

为了研究不同长度不同位置的挡板对回射流的影响，同时考虑到试验中的最大空化长度为0.5倍弦长，分别选择在椭圆鼻水翼吸力面距离导边0.37倍弦长处和0.60倍弦长处设置全展长的挡板。同时为了研究不同长度挡板的抑制效果，一个只有十分之一展长的挡板同样布置在0.37倍弦长处。挡板布置示意如图5所示。椭圆鼻水翼的计算工况由表1给出，表中的k和ω分别为湍动能和比耗散率。不同于平板水翼的第一个工况，椭圆鼻水翼的计算工况为20°C下的饱和蒸汽压力2 338.8 Pa，水翼的攻角为6°，来流速度为5 m/s。

图5 挡板长度与位置分布Fig. 5 Arrangement of the obstacles used in simulations

网格的生成方法与工况一相同，都是通过OpenFOAM自带的blockMesh工具生成背景网格，然后用snappyHexMesh生成加密区。计算采用的网格有33×105，35×105，37×105和40×105这4套，为了使得计算准确同时提高计算效率，最终选取了一套37×105的网格，添加挡板与不添加挡板的网格量都近似在37×105。图6给出了椭圆鼻水翼在中剖面处的网格，其中，图的左右下角，分别是椭圆鼻水翼导边和随边的网格分布。图7为在0.37倍弦长处，十分之一倍展长的挡板局部网格分布。

图6 中剖面以及水翼两边的网格分布Fig. 6 Mesh distribution in mid-plane and near the leading edge and the trailing edge

图7 布置(d)中的挡板附近网格分布Fig. 7 Mesh distribution around the obstacle in configuration (d)

2 数值计算及结果分析

2.1 升阻力系数

首先分析平板水翼的物理模型。图8给出了平板水翼在空化数为0.7时的升力系数历时曲线。从图中可以看出，挡板的存在会改变升力系数曲线的相位，尤其是在空化发展到0.06 s之后，添加挡板的水翼与不添加挡板的水翼呈现出完全相反的相位。

图8 空化数为0.7时的升力系数历时曲线Fig. 8 Time history of lift coefficient at cavitation number 0.7

值得注意的是，在给平板水翼0.37倍弦长处添加挡板的前后，数值模拟计算得到的平均阻力系数几乎没变，不加挡板前的平均阻力系数为0.076 6，添加挡板后的阻力系数为0.079 67，变化不大。而计算得到的升力系数却发生了很大的变化，尽管添加挡板前后，升力系数的值都比较小，但对平板水翼添加挡板使得水翼的平均升力系数由0.003 565上升到了0.010 733，升力系数为无挡板时的3倍。可以得出，在给水翼添加适当的挡板后，升力系数的增幅会远大于阻力系数的增幅，也就是说，挡板的存在，可以以相对较小的阻力代价，换取巨大的升力提高，挡板的合理使用，在工程应用中有着巨大的潜能。在后文对平板水翼表面压力分布和空泡形态的分析中，解释了挡板的存在会使得升力系数得到大幅度提高的原因。

2.2 压力分布

图9和图10分别为平板水翼上表面的平均压力分布和瞬态压力分布。从图中可以看出，挡板的存在改变了平板水翼上表面的压力分布，无论是从瞬态上观察还是直接分析平均压力分布，都能发现挡板的存在降低了挡板后面区域的压力。

图9 平板水翼上表面的平均压力分布Fig. 9 Average pressure distribution on the upper surface

图9的平均压力分布更加能说明问题，在图9中椭圆标记出的区域，挡板的存在会在挡板后的一片区域制造低压区，在椭圆外，添加挡板与不添加挡板的压力区别不大。添加挡板会在平板水翼上表面造成低压区，使得上下表面的压差增大，这很好地解释了为什么挡板的存在会使平板水翼的升力系数大幅度的提高。后文通过对上表面空化形态的观察会发现挡板的存在会抑制片空化向后发展，也因为空化在挡板后面难以发展，压力得以降低，低压区得以形成。图11为平板水翼在0.25 s时的瞬态压力分布云图。

图10 平板水翼上表面的瞬态压力分布Fig. 10 Transient pressure distribution on the upper surface

图11 0.25 s瞬时压力分布Fig. 11 Contour of pressure at 0.25 s

2.3 回射流分析

图12展示了平板水翼空化界面与回射流的相互作用。从图中可以看出，在片空化发展到一定程度时，片空化底部的回射流开始逆着来流方向向上游移动，挡板的存在会阻碍回射流继续向上发展，从而有效抑制片空化的脱落。虽然回射流的速度会比主流速度要小，但回射流速度和主流的速度是同一量级的。

图12 平板水翼周围的空化形态和速度矢量Fig. 12 Cavitation shape and velocity vectors around the hydrofoil

椭圆鼻水翼的抑制效应则更加明显。如前文中提到的空化脱落的两种机制，根据σ/(2α)取值的大小，当σ/(2α)的值比较大时，回射流机制是片空化发展成云空化的主要原因；当σ/(2α)的值比较小时，激波机制是片空化脱落生成云空化的主要原因，两种机制的分水岭在于σ/(2α)=4，基于文中的数值模拟条件计算得出该工况下σ/(2α)=5.11，略高于这个分水岭，可认为回射流为该工况下空化脱落的主要机制。图13给出了椭圆鼻水翼在0.37倍弦长处的挡板对回射流的阻挡现象。在片空化底部生成的回射流受挡板的阻碍后改变了流向，不再逆流而上。与挡板相撞后的回射流不再有继续向上的冲劲，会顺着主流，随脱落的空泡一起翻卷流向下游，并在挡板的后方形成一个涡。此后，片空化不再发展到挡板之后，挡板的存在抑制了空化区的形成以及空化的脱落。

图13 布置(b)工况下挡板对回射流的阻挡过程Fig. 13 Block phenomenon of placed obstacles on the re-entrant jet in configuration (b)

2.4 不同布置挡板的效果

当没有产生空化现象时，在一个光滑的水翼表面安装一小块挡板无疑会增加水翼的阻力，从而导致水翼的升阻比减小，效率降低。而当空化现象发生后，情况就出现了反转。在水翼表面布置合适的挡板可以以较小的阻力代价来获得巨额的升力提升，从而提高升阻比，提高水翼的效率。这是因为在水翼吸力面低压区的局部压力下降到饱和蒸汽压后，压力就因为空化的产生而无法继续下降。而添加合适的挡板可以改变这种情况，挡板的存在可以改变吸力面的压力分布，减小空化区，并在挡板后的一小块区域形成新的低压区，从而有效地提高升力，提高升阻比，提高水翼的效率。图14是图5给出的不同挡板布置下，椭圆鼻水翼在上表面的瞬态压力分布图。

从图14可以看出，不同布置的挡板会不同程度地改变椭圆鼻水翼吸力面的压力分布，且趋势都是类似的，都会在挡板的前端让压力突增，而在挡板后的一小块区域内形成低压区。为了比较不同布置挡板的升阻力变化，将计算得到的升阻力系数在表2中给出。

表2 不同挡板布置下的平均升阻力系数Tab. 2 The calculated average lift and drag coefficients under different configurations

布置a就是不加挡板的工况，布置b相对于它来说，升阻力系数都有提高。遗憾的是，因为椭圆鼻水翼有一定的攻角，升力系数的增幅并没有远高于阻力系数的增幅，但可以肯定的是在空化区设置挡板的确有效地提高了升力系数。在试验中，椭圆鼻水翼的最大空化长度为0.5倍弦长，布置c的挡板设置在了空化充分发展的区域以外，这种情况就如同没有空化产生时在水翼表面添加了多余的附件，不仅阻力系数大幅增加，升力系数也因为挡板未对空化区域产生影响而减小，这种布置的挡板在实践应用中应该尽可能的避免。布置d和布置b一样都将挡板设置在了0.37倍弦长处，但其长度只有展长的十分之一，因为挡板设置在空化区域内，有效地提高了升力系数，同时又不会过多的增加阻力系数，根据升阻比来看，这种布置的挡板是最有希望提高水翼效率的一种布置方式。