文丘里管内空化流动大涡模拟的验证与确认

梁蕴致，龙云，龙新平*，程怀玉

1 武汉大学 动力与机械学院，湖北 武汉 430072

2 武汉第二船舶设计研究所，湖北 武汉 430205

3 武汉大学 水利水电学院，湖北 武汉 430072

0 引 言

随着计算机技术的发展，数值模拟的应用领域越发广泛。为解决模拟结果的准确性问题，部分研究者从误差的来源出发，提出了一套数值模拟结果的验证与确认（verification and validation，V&V）方法。验证用于说明数值方法求解结果与对应方程准确解的偏离程度，确认则用于说明数值解与真实流动情况间的偏离程度[1]。Schaefer等[2]认为数值计算过程中存在几何简化等模型形式、湍流模型系数等参数、离散方法等求解方法、边界条件等输入量、后处理插值等结果处理这5 个方面的不确定性因素，从而使得数值模拟结果与实际流动情况产生偏离。

陈鑫等[3]认为目前几何外形、湍流模型、边界条件等带来的不确定度研究已较为成熟，但针对计算网格的不确定度研究仍然有待完善。目前，对于网格不确定度的研究主要基于估算离散误差的Richardson 外推方法[4-5]，对多套系统加密网格的数值模拟结果进行处理，从而得出采用最密一套网格数值模拟结果的不确定度。另外，现有V&V 方法大多基于RANS 模型的稳态数值模拟，例如Roache[6-7]提出的网格收敛性指数（GCI）方法及其改进版本[8-9]；Wilson 等[10]提出的修正因子法（CF）；Tao 等[11]提出的安全系数法（FSI）。ASME等组织基于上述研究提出了各自的V&V 方法的实施规程[12-13]。张涵信等[14]另辟蹊径，考虑了网格和加权平均误差进行两次加权平均，得出真值，从而得到不确定度。

随着计算机性能的不断提高，在水动力学领域，非稳态、多相流数值模拟研究越来越普遍，大涡模拟（large-eddy simulation，LES）也被越来越多地用于数值模拟研究中，但相关的V&V 方法不够全面，仍有待进一步研究。对于LES，类比稳态RANS 模拟结果的V&V 方法，Vreman 等[15]指出，除网格不确定度外，亚格子模型不确定度也是LES 不确定度的重要组成部分。Klein[16]将这两个准确度阶数假定为2，给出了一种LES 不确定度的三方程计算方法。Xing[17]认为，LES 不确定度还应包括网格与模型耦合带来的不确定度，提出了七方程和五方程的LES 不确定度计算方法。对于空化这一典型的包含相变的多相流动，Long 等[18]及程怀玉等[19]基于RANS 模型对5 套系统加密的网格在水翼空化流动中的不确定度进行了计算，发现非定常空化流动增加了计算误差，从而增大了V&V 方法的实施难度。在Xing[17]对于无空化LES V&V 实施方法的研究基础上，Long 等[20]综合大量计算结果，针对水翼空化流动，将2 个准确度阶数固定，从而将五方程模型简化为三方程模型，并对可变螺距（CP）桨叶进行了五方程和三方程的数值模拟准确度计算[21]。

综上，V&V 方法被广泛应用于计算流体力学的各领域，是检验数值模拟准确度的有效手段，但现有V&V 方法主要基于RANS 模型的稳态外部流动模拟，LES V&V 实施方法的研究较少，尤其是对于有限空间内部空化流动LES V&V 实施方法仍有待研究。为此，本文将对3 种空化数工况下的文丘里管内流动进行数值模拟，基于五方程模型，对有限空间内部空化流运动LES 的验证与确认展开研究。

1 数值模型与边界条件

1.1 几何模型与计算域

本文计算所用文丘里管截面图与尺寸如图1[22-23]所示，喉管直径Dth=10 mm，喉管部分长度Lth=Dth，进口直径Din=5Dth，同时进、出口尺寸相等，该文丘里管的收缩段收缩角为45°，扩散段扩散角为12°。

图1 文丘里管截面图Fig. 1 Sketch view of Venturi

本文计算域如图2(a)所示，为保证流动的充分发展和计算的稳定性，在文丘里管的进口及出口位置均增设了5 倍进出口直径的直管段。本文数值模拟中，以喉管进口截面中心为坐标起点，共设置4 个监测点，分别位于x/Dth=0.5，x/Dth=5，x/Dth=9，x/Dth=13 处的壁面位置，如图2(b)所示。

图2 计算域及监测点位置Fig. 2 Computational domain and positions of monitoring points

1.2 网格划分

基于Richardson 外推方法发展的一系列V&V方法需要在同样的工况下对多套系统加密的网格展开数值模拟计算[24]，网格的系统加密即基于同一套拓扑结构，采用恒定的比例增加各个方向网格数量。为了满足五方程大涡模拟V&V 方法的实施要求，本文数值模拟共使用了5 套网格加细比r为1.2 的系统加密结构化网格，这5 套网格以最密一套网格为G1，最稀疏一套网格为G5，其网格单元数如表1 所示。

表1 网格单元数Table 1 The number of grid cells

同时，为满足大涡模拟的计算要求，网格在喉管和壁面位置对网格进行了加密。图3(a)为G1 喉管位置的壁面网格，从收缩管到喉管网格平滑过渡，尺寸不断减小，从喉管到扩散管，网格尺寸逐渐扩大。图3(b) 为G1 喉管中部截面网格，靠近壁面位置网格高度减小。

图3 文丘里管局部网格划分Fig. 3 Part of grid division of Venturi

1.3 控制方程

文丘里管内空化流动为典型的两相湍流流动，本文采用Mixture 模型将气液两相假定为均质流动，无相间滑动，采用体积分数的变化模拟相关流动参数。湍流模型采用LES 方法，并用Zwart-Gerber-Belamri (ZGB)空化模型模拟气液相变过程。 LES 滤波处理后的N-S 方程为：

式中：p为压强；ui为i方向的速度，其中，在三维笛卡尔坐标系中，下标i 为x轴正方向，j为y轴正方向，k为z轴正方向； δij为Kronecker-delta 函数；ρm和µm分别为均质混合流体的密度和动力黏度，其中ρ1为水密度，µl为水动力黏度，αv为蒸汽体积分数，ρv为蒸汽密度；式(2)带有波浪型上划线的项表示滤波处理项，而τij为亚格子应力，代表被过滤的小尺度脉动的作用。基于Boussinesq 假设，亚格子应力项如式(5)所示，其中µSGS为 亚格子湍流黏度，S˜ij为应变率张量。本文采用的亚格子应力模型为Dynamic Smagorinsky-Lilly 模型。

空化流动模拟方面，大量算例证明，基于Rayleigh-Plesset 方程得出的ZGB 空化模型[25]有较好的普适性，且模拟精度和收敛性均较好，复杂流动中也表现出了较好的收敛性，故本文也选用ZGB 模型模拟空化过程。ZGB 模型的质量输运方程为：

式中，me和mc分别为蒸发源项和凝结源项，表达式为：

式中：pv为饱和蒸汽压；RB为气泡直径；Fvap和Fcond为经验系数；αnuc为成核点的体积分数。

1.4 边界条件

本文仿真分别计算了空化数为2.086，0.765及0.263 这3 种工况。空化数定义为：

式中：pout为出口压力；uth为喉管中间截面平均流速。当空化数σ 为2.086 时，试验中未观察到空化现象，而空化数为0.765 及0.263 时，试验中均发现有较为剧烈的空化现象，即本文数值模拟同时涵盖了空化和无空化流动。

本次计算涉及的流体为水和水蒸气，均采用25 ℃时的物性参数，饱和蒸汽压为3 540 Pa，水的密度为997 kg/m3，黏度为8.899 ×10−4kg/(m∙s)，水蒸汽的密度为0.023 08 kg/m3。

表2 为本文数值模拟的边界条件，其来源为试验测量数据[22-23]。数值模拟时，进口设置为速度进口，出口设置为压力出口，壁面设置为无滑移壁面。对应于本文数值模拟采用的网格加细比，当采用更细密的一套网格时，时间步长除以1.2，即G5 网格采用的时间步长为1.44×10−5s，G4 网格采用的时间步长统一为1.2×10−5s，依此类推。

表2 计算边界条件Table 2 Boundary conditions

本文计算使用CFX 软件，并采用192 线程服务器进行并行计算，内存128 G，计算共耗时1 440 h。

为排除时间步长的选取对非稳态数值模拟计算结果的影响，选取G4 网格的无空化工况，另外选取了1.44×10−5s 以及1×10−5s 两个时间步长，进行了时间步长无关性验证。相比于时间步长为1.2×10−5s，两时间步长下进出口压力比m的变化均在0.5%以下，4 个监测点的压力系数Cp变化也都在3%以下，可以认为，本次数值模拟的各套网格满足时间步长无关性。

压力比m是文丘里管主要性能指标之一，为出口压力pout与进口压力pin之比，定义如式（10）。对于进出口等径的文丘里管，该参数表现了文丘里管中的压能损失。压力系数表达式如式(11)。

式中，pmon为瞬态计算结果趋于稳定后最后2 个周期内的监测点压力平均值。

2 计算结果与试验对比

2.1 进出口压力比对比

图4 为3 种空化数下，文丘里管压力比的试验值及5 套网格的数值模拟结果[22-23]，而G1～G5分别为从最密一套网格到最稀疏一套网格共计5 套网格的模拟值。

图4 5 套网格模拟压力比与试验值的对比Fig. 4 Comparison of pressure ratios between experimental results and simulations in five grids schemes

图4(a)中，空化数σ=2.086 时，随着网格加密压力比上升，在网格G2 处达到峰值。图4(b)中，空化数σ=0.765 时，模拟结果不断降低，逐渐向试验值靠拢。图4(c) 中，空化数σ=0.263时，压力比变化幅度最小。所有网格计算结果与试验值间的误差均在5%以内，表明本文数值模拟的3 种空化数下，压力比这一全局量的模拟结果与试验值均吻合较好。

2.2 监测点压力系数对比

文丘里管几何模型较为简单，但其独特的收喉扩结构使喉管和扩散管部分的流场较为复杂，对于这一典型部分内部流动结构的研究很有必要。图5 分别为σ=2.086，0.765，0.263 时，G1～G5网格下压力系数Cp的模拟值与试验结果。

图5 文丘里管壁面监测点压力系数Fig. 5 Pressure coefficient of monitoring points on the wall of Venturi

图5(a)和图5(b)中，随着网格数的增加，数值模拟结果更加接近试验值，最密一套网格G1 计算结果与试验值间的误差均很小，与试验结果吻合较好。而图5(c)中，x/Dth=0.5，x/Dth=5 监测点压力系数误差很小，而监测点x/Dth=9，x/Dth=13 处误差相对较大。空化较弱或无空化时，各监测点压力系数模拟值与试验值均较为吻合，而空化较剧烈时，喉部和扩散管前部监测点仍较为吻合，仅扩散管中后部2 个监测点压力系数偏差较大。

上述对比表明，本文所计算的3 种工况下，进出口压力比这一全局量的试验值与模拟值均较为吻合；对比各监测点压力系数这一局部量时，仅空化数σ=0.263 时，扩散管中后部2 个监测点的压力系数误差较大，其余各监测点和各工况下的模拟结果与试验值均较为吻合。即当无空化或空化较弱时，本文所采用数值模拟方法能够较为准确地反映实际流动；空化较剧烈时，部分局部量偏差较大，这与求解方法和空化模型都有一定关系，有较大的改善空间，但全局量仍较为符合。

3 V&V 实施方法与结果

3.1 验 证

LES 的误差可以分为数值误差δSN及模型误差δSM两个部分[17]：

目前，广泛认为由直接数值模拟（DNS）方法得到精度较高的离散解可以作为数值基准值，但并非所有问题都能够采用DNS 求解，故在大部分问题中， ϕc为需要求解的未知量。根据二元泰勒级数展开法及广义Richardson 外推方法，同时参考已有规程[12-13]对RANS 模型的误差计算方法，展开两误差项并略去高阶量，则LES 方法的计算误差可以表示为

式中：cN和cM分别为数值误差系数和模型误差系数；pN和pM分别为数值误差准确度阶数和模型误差准确度阶数；Δn为滤波尺寸，下标n表示对应数值模拟所采用的网格；h∗n为综合考虑当地网格尺寸和瞬态计算时间步长得到的时空尺度，表达式如下：

式中：hn为当地网格尺寸；Δtn为时间步长。

对于系统加密的几套网格，当采用的加细比为r时，则计算时相邻网格有如下关系：

如要求解式(13)，理论上应有至少5 个等式组成方程组。本文以G1 网格的当地尺寸为网格尺寸h，采用该套网格时，数值模拟的时间步长为Δt，由5 套系统加密网格得出的数值解为 ϕn，可构成的方程组如下：

方程组为复杂非线性方程，借助Matlab 最优化工具求解[26]，采用信任域方法和Jacobian 矩阵迭代法，得出压力比m和各监测点压力系数Cp的准确度阶数如下。图6 为压力比m的准确度阶数。如图所示，当σ=2.086，0.263 时，pN，pM比较接近，但σ=0.765 时，pN，pM明显低于其他2 个工况。图7 为Cp的准确度阶数。如图所示，同一工况下，不同监测点处的准确度阶数有较大差别，且同一监测点处，不同工况下的准确度阶数也不相同。

图6 压力比准确度阶数Fig. 6 The observed order of accuracy of pressure ratio

图7 压力系数准确度阶数Fig. 7 The observed order of accuracy of pressure coefficient

文献[21]中，pN，pM在各个工况下相差较小，可以固定这2 个准确度阶数并将五方程求解简化为三方程；但对于文丘里管这种速度、压力变化剧烈的内部流动，准确度阶数会随监测点位置、流动工况的改变而改变，且变幅较为明显。故该方法存在不足，需要寻找准确度阶数随监测点位置和流动工况变化的规律，从而加以改进。

LES 验证需要进一步计算验证不确定度，计算方法为

式中：ULES为LES 的验证不确定度，表明 LES 结果与精确解之间的偏离程度；FSLES为不确定度安全系数，根据Roache[6]的推荐，该系数可以取为1.5 或者3。考虑到文丘里管内空化流动为三维非定常相变过程，计算中会带来较大的不确定性，故当σ=0.765，0.263 时，取安全系数为3，而σ=2.086 时，取安全系数为1.5。

进出口压力比m这一全局量的不确定度ULES均较小（图8），σ=0.765 时，ULES相对较大；对应于图4 中，σ=0.765 时模拟结果与试验值之差也最高，这说明验证不确定度ULES的大小能够反映压力比m计算误差的大小。压力系数Cp的验证不确定度如图9 所示，在σ=2.086，0.765 工况下，4 个监测点位置的验证不确定度ULES均较小；当σ=0.263 时，x/Dth=0.5 及x/Dth=5 这2 个监测点压力系数的ULES较小，而x/Dth=9 和x/Dth=13 这2 个监测点处ULES明显增加。对应图5(c)，这2 个监测点处的计算值与试验值的误差也偏大，且误差最大的x/Dth=9 处，ULES也最大，即验证不确定度ULES的大小也能够很好地反映压力系数Cp的误差。

图8 压力比验证不确定度Fig. 8 The verification uncertainty of pressure ratio

图9 压力系数的验证不确定度Fig. 9 The verification uncertainty of pressure coefficient

3.2 确 认

验证结束后，还需要对LES 结果进行确认，即计算并对比分析确认不确定度和对比误差。Long 等[21]指出，相较于RANS 模型的验证过程，LES 验证过程已经将模型误差纳入了考虑范围，因此 LES 的确认不确定度可以表示为

式中：Uv为确认不确定度；UD为试验的不确定度。本文中UD为试验所用压力测量仪表的不确定度，压力漂移值为压力表量程的0.5%。

当对比误差（试验值D与LES 模拟值之差）小于确认不确定度时，则认为在该不确定度Uv的水平下，LES 的确认得以实现。对比误差为

压力比m的确认不确定度及对比误差E如图10所示。在3 种空化数的工况下，对比误差E均远小于确认不确定度Uv，即3 种空化数工况下的进出口压力比m均实现了对应不确定度下的确认。

图10 压力比的确认不确定度及对比误差Fig. 10 The validation uncertainty and comparison error of pressure ratio

各监测点处压力系数的确认不确定度与对比误差如图11 所示。在空化数σ=2.086，0.765 时，各监测点压力系数的对比误差均小于对应的确认不确定度；而当σ=0.263 时，x/Dth=9 和x/Dth=13 这2 个监测点处的对比误差大于确认不确定度。即σ=0.263 时，x/Dth=9 和x/Dth=13 这2 个监测点处的压力系数没有实现对应不确定度水平下的确认，其余各工况和各监测点的压力系数均实现了对应不确定度水平下的确认。

图11 压力系数的确认不确定度及对比误差Fig. 11 The validation uncertainty and comparison error of pressure coefficient

在空化数σ=0.263 时，文丘里管处于剧烈的空化状态，其蒸汽体积分数云图、试验高速摄像和截面网格如图12 所示。图12 (a)和图12(b)中，空化云延伸范围较长，x/Dth=9 监测点位于空化云的末端，在此位置空化云脱落溃灭，同时，部分脱落的空化云还会运动到x/Dth=13 位置处再溃灭，空化云的脱落溃灭在壁面位置造成了不稳定回流，进而导致这2 个位置处的不确定度明显上升；此外，监测点x/Dth=9 和x/Dth=13 位于扩散管的中后部，说明当空化数σ=0.263 时，本文的数值方法仍有较大的改善空间，需要对网格及求解方法做进一步调整。

图12 σ=0.263 时文丘里管空化云形态及部分G1 网格Fig. 12 Cavitation clouds and part of G1 grid scheme in Venturi at σ=0.263

4 结 论

本文采用LES 耦合ZGB 空化模型，对3 种空化数工况下的文丘里管内瞬态空化流动进行了数值模拟，运用Matlab 最优化工具求解五方程V&V 方法中的数值基准和准确度阶数，进而计算出验证不确定度ULES、确认不确定度Uv和对比误差E，进行了文丘里管内部空化及无空化流动大涡模拟的验证与确认，得出如下结论：

1) 最密的一套网格数值模拟结果与试验结果的直接对比表明：3 种空化数工况下的压力比这一全局量的模拟值与试验数值均较为吻合。无空化或空化较弱时，监测点压力系数与试验结果均吻合较好；而当σ=0.263，空化较为剧烈时，x/Dth=9 和x/Dth=13 处的压力系数误差偏大，其余2 个监测点数据吻合较好，数值模拟方法有较大的提升空间。

2) 验证结果表明，压力比这一全局量的验证不确定度均较小；当空化较弱或无空化时，监测点压力系数验证不确定度均较小，而空化较为剧烈时，扩散管中后部监测点压力系数的验证不确定度较大。总体上验证不确定度较大时，对比误差也较大，说明验证不确定度能够反映数值模拟的准确程度。

3) 确认结果表明，在大部分工况和监测点处，均实现了对应验证不确定度水平下压力比的确认；对比监测点压力系数时，仅在空化数σ=0.263，监测点x/Dth=9 和x/Dth=13 处没有实现确认，说明空化数适中或无空化时，使用细化的网格，LES 方法能较好地模拟非定常空化流动。但空化数较小和空化极为剧烈时，需要针对该工况改善网格，并调整求解方法，以精确捕捉空化的演化过程。

4) 目前，实施有限空间内部空化流动的LES验证与确认时，压力比等全局量的验证与确认较易实现，但空化较为剧烈时，完全实现各监测点压力系数的验证与确认还需要继续完善数值模拟求解方法，未来还将进一步修改空化模型，使之符合更大范围的空化工况。此外，为实现五方程V&V 方法，一个工况需要起码5 套系统加密网格的计算结果，计算资源消耗较大，而三方程方法还需寻找准确度阶数随监测点位置和流动工况变化的规律，以做进一步改进，拓展V&V 方法的应用范围。