基于当地流动特征的Schnerr-Sauer空化模型改进研究

胡俊， 侯夏伊， 于勇

(北京理工大学 宇航学院，北京 100081)

空化是指液体中的压力降低到饱和蒸汽压时，产生充满气体空穴[1]，是一种因流体因素作用而在波体内部或在波体与固体界面上发生的液体与其蒸气的相变过程与现象[2-3]. 空化涉及非定常、可压缩、相间传质和湍流脉动等许多复杂的流动现象，并且会引起水下翼型、水中机械的材料损伤，产生空泡噪声[4-8]. 研究空化现象具有很现实的工程意义，发展空泡流的数值模拟技术是目前空化研究领域中的热点.

研究空化问题的数值模拟方法可分为状态方程法(barotropic model)和输运方程法(transport equation model,TEM). 其中，输运方程法通常不考虑流场的可压缩性，假设气液两相速度相等，需要添加额外的质量输运模型即空化模型来描述气液两相间的相变规律. 目前常用的空化模型均涉及经验参数，如气泡的数密度、非冷凝气体质量分数和气泡半径等，这些参数的取值很大程度上影响了数值模拟的准确性. Asnaghi等[9]基于大涡模拟(LES)的方法，引入当地瞬时剪切应变率修正了考虑非冷凝气体的Schnerr-Sauer空化模型，对二维NACA0009 MOD水翼的定常空化流场计算. 结果表明该修正方法可以有效消除考虑非冷凝气体的Schnerr-Sauer空化模型中经验参数的取值对数值模拟结果的影响.

基于这一思路，在RANS框架下模拟空化流时，也可以引入平均剪切应变率修正空化模型. 本文针对考虑非冷凝气体的Schnerr-Sauer空化模型，将修正后的空化模型通过UDF嵌入Fluent18.2商业软件，采用Mixture多相流模型和Realizablek-ε湍流模型对攻角为2.5°，Re=2×106的二维NACA0009 MOD水翼的定常空化流场进行数值模拟. 通过与已有的实验结果对比，分析了经验参数的取值对采用两种空化模型的计算结果的影响，验证了引入当地流动特征修正空化模型的方法在二维定常空化流计算中的可行性.

1 数学模型和数值计算方法

1.1 控制方程

Mixtrue模型是一种简化的多相流模型，将流体中各相视为相互混合的单一流体，通过求解混合物的连续方程，动量方程及气相的体积分数方程来模拟气液两相流体的运动. 其中，混合相连续方程及动量方程为

(1)

(2)

式中：ρm、ρl和ρv为混合相密度、液相密度和气相密度；μm为混合相动力黏度；Ui、Uj为混合相速度. 其中，混合相密度及动力黏度为

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

(3)

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

(4)

式中：αv为气相体积分数；μl和μv分别为液相动力黏度及气相动力黏度.

除混合相连续方程和动量方程外，还需要额外的气相体积分数方程来描述气液两相间的质量传递

(5)

1.2 空化模型

1.2.1气泡动力学基本方程

Rayleigh-Plesset单气泡动力学模型假设液体区域内均为球形气泡，气泡间无相互作用:

(6)

式中：RB为气泡半径；PB为气泡表面压力；P为环境压力；S为气液分界面上的表面张力系数；ϑl为液相运动黏度；t为时间. 忽略单气泡动力学方程中的黏性项，表面张力项及高阶导数项，可得到单个气泡半径变化率为

(7)

1.2.2Schnerr-Sauer空化模型

Schnerr-Sauer模型基于Rayleigh-Plesset单气泡动力学方程，将式(1)与式(5)联立，得到气相质量变化率与气相体积分数变化率之间的关系为

(8)

其中气相体积分数αv与气泡数密度n及气泡半径RB的关系为

(9)

由式(9)可得到气相体积分数变化率与气泡半径变化率之间的关系

(10)

Schnerr-Sauer空化模型的最终表达式为

(11)

式中：Cevap=1.0和Ccond=-1.0分别为蒸发项及凝结项系数；Pv为饱和蒸气压. 由于Schnerr-Sauer空化模型没有考虑到实际流体中存在的非冷凝气体对空化初生及发展的影响，且需要确定模型中参数气泡数密度n的值，在模拟空化流场时存在一定的局限性.

1.2.3考虑了非冷凝气体的Schnerr-Sauer空化模型

经研究表明，实际流体中存在的非冷凝气体对空化的初生与发展起到了重要作用，该部分气体在空化发生的过程中充当空化核，使液体抗拉强度下降，当液体压力下降到一定值时，相变首先从空化核处发生. 由于原始的Schnerr-Sauer空化模型没有考虑非冷凝气体对空化过程的影响，在模拟空化流时存在一定缺陷. 为了提高Schnerr-Sauer空化模型模拟空化流的能力，本文将流场视为气体、液体及非冷凝气体三相混合物，推导出一种考虑了非冷凝气体的Schnerr-Sauer空化模型.

Asnaghi等假设空化流场中的单个气泡由非冷凝气体发展而来，故流场中气泡数密度n与非冷凝气体数密度nnuc相等，并将非冷凝气体的体积fnuc及体积分数αnuc定义为[8]：

(12)

(13)

参考原始的Schnerr-Sauer空化模型对气相体积分数的定义即式(9)，将修正后的Schnerr-Sauer空化模型的气相体积分数与非冷凝气体体积分数之和定义为[8]

由式(14)可得到单个气泡半径的表达式

(15)

由于非冷凝气体对蒸发过程的影响远大于空泡溃灭过程，本文仅在Schnerr-Sauer空化模型的蒸发项中引入非冷凝气体，最终得到考虑非冷凝气体的Schnerr-Sauer空化模型：

(16)

式中：dnuc、fnuc、αnuc和nnuc分别为非冷凝气体直径、体积、体积分数以及数密度.

1.2.4引入剪切应变率修正模型系数

t∞=c/U∞

(18)

式中：Ce-mod为修正后的蒸发系数；t∞为主流时间尺度；c为翼型弦长；U∞为来流速度.

2 计算网格及边界条件

计算采用了和试验相同的二维NACA0009 MOD型水翼及流动条件，翼型攻角为2.5°，翼型弦长c=100 mm，对应雷诺数Re=2×106. 图1给出了计算域及其边界条件，计算区域入口距翼型前缘2.0c，出口距翼型尾缘4.0c. 如图2所示，整个计算域被划分为9个块，总网格数为34 720. 翼型附近的区域采用C型结构化网格划分，可以较好地匹配翼型的形状，并在翼型周围的近壁区域进行网格加密，近壁面y+值为 30～200之间，满足壁面函数要求，计算残差为10-4.

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

图2 计算网格Fig.2 Computational meshes

数值计算边界条件采用速度入口及压力出口，上下边界均为对称条件，翼型表面采用绝热、无滑移固壁条件. 为了保证与实验条件一致，入口速度为U∞=20 m/s，计算了4种空化数下的流场.

3 结果与讨论

3.1 非冷凝气体数密度n对考虑非冷凝气体的Schnerr-Sauer空化模型的影响

为了探究非冷凝气体数密度n的取值对考虑非冷凝气体的Schnerr-Sauer空化模型的影响，在空化数σ分别为0.75、0.80、0.85、0.90时，取dnuc为定值10-5，计算了n取106、108、1010和1012时的空化流场.

图3所示的计算结果表明：随着空化数σ由0.75增大到0.90，空化流场的流动状态更加稳定，计算出的空穴长度也有所减小. 同时，参数n的取值对空穴长度的影响随空化数的减小而减小. 4种空化数下计算出的空穴长度均随n的增加而增加，其中，n=1012时，空穴长度达到最大，且压力系数的分布最接近于实验值. 但由于引入当地流动特征修正前的空化模型的蒸发项系数为1，即使n取较大值，即n=1012时，在4种空化数下计算出的空穴长度都明显小于实验值. 当σ=0.75、n=1012时，采用该空化模型能够捕捉到空穴闭合区域的压力峰值及逆压梯度；当σ=0.80、n=1012时，这种逆压梯度较小且与实验值相差较大；当σ=0.85、0.90时，即使非冷凝气体数密度n取较大值，采用该空化模型也难以捕捉到空穴闭合区域的逆压梯度及压力峰值.

图3 采用系数修正前模型计算出的翼面压力分布Fig.3 Pressure distribution of hydrofoil calculated by using the original model

3.2 非冷凝气体直径dnuc对考虑非冷凝气体的Schnerr-Sauer空化模型的影响

为了研究非冷凝气体直径dnuc取不同值时，考虑非冷凝气体的Schnerr-Sauer空化模型的表现，在空化数σ=0.75，n=108和1010时，改变非冷凝气体直径dnuc，使dnuc=10-3、10-4、10-5、10-6和10-7m，对二维水翼的定常空化流场进行数值模拟，如图4所示.

图4 采用系数修正前模型计算出的翼面压力分布Fig.4 Pressure distribution of hydrofoil calculated by using the original model

图4表明dnuc的取值对数值模拟结果的影响较小. 当n=108，dnuc取10-3～10-5时，以及n=1010，dnuc取10-3～10-6时，计算出的翼型表面空穴长度及压力系数分布接近. 但即使dnuc取值较大，即dnuc=10-3时，压力系数分布与实验值差距依然较大，翼型表面的空穴长度明显小于实验值，且难以捕捉到空穴闭合区域的压力峰值及逆压梯度. 当n=108，dnuc=10-6、10-7以及n=1010，dnuc=10-7时，由于非冷凝气体直径dnuc取值过小，导致非冷凝气体体积分数过小，空穴长度也会明显缩短. 相比非冷凝气体数密度n，考虑非冷凝气体的Schnerr-Sauer空化模型对非冷凝气体直径dnuc的取值并不敏感.

3.3 引入平均剪切应变率修正后的考虑非冷凝气体的Schnerr-Sauer空化模型

由于考虑非冷凝气体的Schnerr-Sauer空化模型中的参数非冷凝气体数密度n的取值对数值模拟结果的影响较大，非冷凝气体直径dnuc的取值对数值模拟的结果影响较小，本文仅探究了当参数n取不同值时，引入剪切应变率修正后的考虑非冷凝气体的Schnerr-Sauer空化模型的表现. 在空化数σ分别为0.75、0.80、0.85和0.90时，取dnuc为定值10-5，计算了n=106、108、1010和1012时的空化流场，如图5所示.

当σ分别为0.75、0.80和0.85，非冷凝气体数密度n取106、108、1010和1012时，采用系数修正后的空化模型清晰地捕捉到了空穴闭合区域的压力峰值及逆压梯度；当σ=0.90，n=106时,难以捕捉到这种现象. 当n=106时，由于非冷凝气体密度n的取值过小，采用修正蒸发项系数后的空化模型计算出的压力分布虽比采用修正前的模型计算出的结果更接近实验值，但空穴长度仍明显小于实验值.

当非冷凝气体数密度n取108、1010和1012时，引入当地流动特征，采用平均剪切应变率修正考虑非冷凝气体的Schnerr-Sauer空化模型的蒸发项系数可以有效地消除参数非冷凝气体数密度n对空化流数值模拟的显著影响，准确预测了空穴长度及翼面压力分布，提高了该空化模型的适用性与准确性.

图5 采用系数修正后模型计算出的翼面压力分布Fig.5 Pressure distribution of hydrofoil calculated by using the modified model

4 结束语

在考虑非冷凝气体的Schnerr-Sauer空化模型中，经验参数非冷凝气体数密度n的取值对数值模拟结果影响较大，非冷凝气体直径dnuc的取值对数值模拟结果的影响较小.

引入当地流动特征，采用剪切应变率修正考虑非冷凝气体的Schnerr-Sauer空化模型后，参数非冷凝气体数密度n的选取对数值模拟的结果不再有显著影响，同时，该修正方法提高了考虑非冷凝气体的Schnerr-Sauer空化模型的准确性，准确地预测了空化造成压力梯度突变的位置.