涡轮机械中的空化流动会导致性能损失和叶片载荷的变化。这些后果可以归因于准稳态空化效应,主要取决于空化结构的时间平均形态。此外,在空化条件下的涡轮机械还会受到瞬态现象的影响,如可压缩性效应、流量波动、噪声、振动、侵蚀等,其中空化流动的非稳态行为和两相结构起到了影响作用。因此,研究和建模空化的非稳态行为对于估计涡轮机械中的液压非稳态性及其相关效应至关重要。
在这种背景下,已经开发了一些考虑空化两相结构和非稳态行为的数值模型。其中大多数基于单一流体方法:忽略了液相和气相之间的相对运动,将液气混合物视为具有可变密度的均质介质。混合物的密度直接与局部空隙率相关,并通过不同的方法进行处理,例如通过状态方程处理,使用附加方程将空隙率与气泡簇的动态演化联系起来,或通过液气之间的质量传递定律处理多物种方法。最后一种模型可以用于全两流体方法,其中包括相间相对运动。
主要的数值难点在于压力场与空隙率之间的强耦合,以及两相介质的强可压缩性与纯液体流动的准不可压缩行为的共存。此外,湍流对非稳态两相可压缩流动的影响尚未被很好地了解。
本文展示了使用二维代码进行的非稳态空化流动模拟。该代码基于Delannoy和Kueny为无粘流体开发的模型。Reboud和Delannoy、Reboud等人,以及Coutier-Delgosha等人进行了几项物理修改,以扩大应用范围并改进物理建模。
在该数值代码中,空化的两相特性通过引入一种强烈将流体密度与压力变化联系起来的重力方程来处理。从数值角度来看,这种物理方法与基于SIMPLE算法的压力校正方案相关,并进行了适当的修改以考虑空化过程。该模型已经在许多案例中得到了验证,例如文丘里管、翼型或叶片级联。数值结果与实验结果显示出良好的一致性,特别是代码能够可靠地模拟非稳态空化流动的循环自振荡行为。
这种循环空化行为所受控制的复杂非稳态机制受到所采用的湍流模型的强烈影响。本文的目的是研究不同湍流模型(标准的k-ε RNG模型;修正的k-ε RNG模型;考虑或不考虑可压缩性效应的k-ω模型)对空化流动数值模拟的影响。事实上,由于模型中采用了正压状态方程,汽化和凝结过程主要对应于高度可压缩的流动区域。因此,在湍流模型中考虑气液混合物的可压缩性至关重要。
本文描述了在数值模拟中应用的湍流模型,并展示了在二维文丘里型截面中获得的不同数值结果的比较,该截面的实验行为已由Stutz和Rebouch研究。在这种几何结构中,流动特征表现为不稳定的空化行为,伴随着几乎周期性的蒸汽云脱落。
除了标准的k-ε模型外,还应用了另外三种湍流模型,以研究它们对两相湍流流动的影响。主要应用了Wilcox提出的包含可压缩性考虑的湍流模型,并获得了令人鼓舞的结果。文中展示了与空隙率和腔内速度的实验测量值的比较。
本文讨论了这些湍流模型预测复杂两相流动的能力,并研究了可压缩性效应的作用。
本研究采用了一种基于先前数值和物理研究的单一流体模型。在计算域内,流体密度ρ根据正压状态方程ρ(P)随局部静压的变化而变化(见图1)。当某个单元中的压力高于蒸汽压附近的值(P>Pv+ΔPv/2)时,流体被视为纯液体,整个单元被液体占据,密度ρl通过Tait方程计算。如果压力低于蒸汽压附近的值(P<Pv-ΔPv/2),单元内充满蒸汽,密度ρv由理想气体定律(等温假设)确定。在纯蒸汽状态和纯液体状态之间,单元内为液体/蒸汽混合物,该混合物被视为单一流体,具有可变的密度ρ。密度ρ与空隙率α直接相关,α=(ρ-ρl)/(ρv-ρl),表示混合物中蒸汽所占的局部比例。
图1 水的正压状态方程ρ(P),水温20℃
为了模拟混合物状态,正压状态方程在蒸汽压附近(ΔPv/2的范围内)提供了一个平滑的链接。该方程的主要特征是其最大斜率1/A²min,其中A²min=∂P/∂ρ。Amin因此可以解释为混合物中的最小声速。该值在先前的研究中进行了校准。结果表明,该最佳值与水动力条件无关,对于冷水,大约为2m/s,对应的蒸汽压Pv=0.023bar,对应的ΔPv约为0.06bar(图1中的虚线)。本文后续计算中均采用了这一数值。
由汽化和凝结过程产生的质量通量通过正压状态方程隐式处理,无需额外的假设。对于动量通量,该模型假定局部液体和蒸汽的速度相同:在混合区中,蒸汽结构被认为完全随主流携带。这一假设通常用于模拟片状空化流动,在这种情况下,界面被认为处于动态平衡状态。因此,相间的动量传递与质量传递密切相关。
为了求解与上述正压状态方程相关的时间相关雷诺平均纳维-斯托克斯方程,该数值代码在二维结构化曲线正交网格上应用了SIMPLE算法,并对其进行了修改以考虑空化过程。时间离散化采用隐式方法,采用了Zhu提出的HLPA非振荡二阶对流格式。该数值模型在Coutier-Delgosha等中有详细描述:方法进行了完整的验证,并广泛研究了数值参数的影响。本文主要描述了应用的不同湍流模型。
边界条件。在代码中,速度场在计算域入口处被设定为已知值,而在出口处设定静压。沿固体边界,湍流模型与壁面定律相关联。
初始瞬态条件。为了启动非稳态计算,应用以下数值程序:首先,进行一个稳态步骤,并设置足够高的出口压力,以避免整个计算域内出现任何蒸汽。然后,在每一个新的时间步中逐渐降低该压力,直到达到对应于所需空化数σ的值。在压力下降过程中会出现蒸汽。之后,整个计算过程中空化数保持恒定。
数值模拟是在一个文丘里型截面上进行的,该截面的收敛和发散角分别约为18度和8度(如图2(a)所示)。文丘里底部从喉部向下游的形状模拟了具有倒角前缘几何形状的诱导器叶片吸力面,弦长为Lref = 224 mm。
图2 (a)文丘里管外形段的曲线正交网格(160×50个网格)
(b)喉部区域网格的分布情况
根据实验观察,在这种几何结构中,流动表现出非稳态的空化行为,伴随着准周期性波动。每个循环包括以下连续步骤:附着的片状空化从文丘里喉部开始增长。在空化闭合处产生回流射流,并沿着文丘里底部向空化上游端流动。回流射流与空化表面的相互作用导致空化破裂。生成的蒸汽云随后被主流携带,直到其塌陷。
在空化数σ约为2.4(基于时间平均的上游压力)和入口速度Vref = 7.2 m/s的情况下,实验观察到的蒸汽脱落频率约为50 Hz,空化长度为45-65 mm。
标准计算网格由160×350个正交单元组成(如图2(a)所示)。在喉部之后的主流方向上对网格进行了特殊拉伸,以便有效地模拟两相流区域:在这个方向上使用了大约50个网格点来模拟后续获得的45 mm长的平均空化(如图2(b)所示)。在另一个方向上,靠近墙壁处也进行了拉伸,以便在第一个网格点处获得边界层的无量纲参数y1,范围在30到100之间,并使用标准的壁面定律。为了提高空化区域的精度,文丘里截面底部的网格比上部更精细:在以下章节中获得的空化厚度包含大约30个网格单元。
(a)标准k-ε RNG模型
在数值代码中应用的第一个模型是标准的k-ε RNG模型,并与壁面定律相结合。在此模型中,雷诺方程中应用的有效黏度定义为μ = μt + μl,其中μt = ρCμk²/ε 是湍流黏度,Cμ = 0.085(Yakhot 等)。
该模型最初用于完全不可压缩的流体,对于高度可压缩的两相混合物,未在此应用特别的修正。因此,流体的可压缩性仅通过湍流方程中的平均密度ρ变化来考虑。
使用该模型时,实验中观察到的不稳定空化行为无法得到正确的模拟:在经历了初始的空化长度瞬态波动后,数值计算导致了空化片的准稳态行为,其整体上趋于稳定(图3)。
图3空腔长度的时间演变。横坐标表示时间,纵坐标表示空化隧道中的X位置。颜色代表密度值:纯液体为白色,蒸汽为红色至深蓝色。在给定的时间和位置,颜色表示空化隧道相应横截面中的最小密度。计算条件:σ ≈ 2.4;Vref = 7.2 m/s;网格分辨率为160×50,时间步长为Δt = 0.005 Tref (Tref = Lref/Vref)。
由此产生的空化长度与实验观测值相比要小得多(在这种情况下,误差超过50%)。此外,与Stutz和Reboud通过双光学探针获得的实验数据进行比较,数值计算中获得的空化片大部分区域的平均空隙率被高估。计算结果显示,空化片上游部分的时间平均空隙率很高(超过90%),但在尾流中骤然下降至0%,而测量的空隙率从未超过25%,且从空化片的上游端到尾流逐渐减小。
与实际情况的不符似乎与空化片后部湍流黏度的过高预测有关。云状空化过程的循环行为与空化闭合处的回流射流发展密切相关。事实上,湍流流动模拟中的主要问题在于沿固体壁面的逆流过早地被移除:回流射流过早停止,未能导致空化片的破裂。
值得注意的是,[10]中报告的数值测试证实,模型中应用的网格尺寸、空间格式和时间离散化并不会改变结果。即使使用更精细的网格(264×390)和一阶或二阶精度的时间离散化格式进行计算,仍然会导致空化片的完全稳定化。
(b)修正的k-ε RNG模型
为了改进湍流建模,并尝试更好地模拟回流射流行为和蒸汽云脱落,我们通过简单地降低混合物的湍流黏度,特别是在低空隙率区域,来修正标准的k-ε RNG模型:
其中,
事实上,根据实验结果[15],回流射流主要由液体组成(α趋近于0),因此减少混合物的湍流黏度会对模拟产生显著影响。现在可以预测出非稳态的回流射流行为,并且蒸汽云的脱落得到了良好的模拟(见图4)。
图4 混合物粘度的修正(n=10)
图5展示了在这一非稳态计算中观察到的瞬态演化。图5(a)展示了在某一时刻文丘里型管道各个截面内的最小密度值。与图3相比,它提供了有关蒸汽云脱落过程的信息:脱落的空化部分清晰可见,并且可以评估波动频率。此外,还提供了每个截面内的最大空隙率。曲线5(b)显示了瞬态上游压力的演化。
图5 文丘里型管道内非稳态空化流动的瞬态演化。(a) 空腔长度的时间演化(横坐标),长度以纵坐标为刻度。附着空穴和云状空穴的瞬时密度分布在T=11Tref时绘制在左侧(速度矢量仅在两个方向上以2个单元绘制1个单元)。(b) 上游压力的时间演变。
片空化的实验自振荡行为得到了正确的模拟(图5)。在这种情况下,得到的是一个波动的空化,其平均长度约为45 mm(即 Lcav / Lref = 0.2)。平均波动周期约为0.59 Tref,相应的脱落频率为55 Hz。长度和频率均与测量值(分别为45 mm和50 Hz)非常吻合。基于参考速度和平均空化长度的减频约为0.33,略微高于经典的斯特劳哈尔数 St = 0.3(实验测量值为0.28)。瞬态演化几乎是周期性的,但有一些随机扰动影响了振荡的规律性。这一行为与实验观察一致,实验指出空化周期相当规律,频率围绕一个中心值波动[14]。
进行了多项测试,以研究数值参数对非稳态空化结果的影响。包括网格大小、时间步长、时间离散化方案的阶数、ρv/ρl比率以及正压状态方程最大斜率Amin的影响,详细内容见[10]。本文中的表1报告了关于网格分辨率和时间步长影响的结果。网格大小的影响似乎很小,只要网格足够精细:使用两种最精细网格时,空化振荡频率几乎保持不变。相比之下,时间步长的影响无法完全消除:使用Δt = 0.005 Tref时,周期频率的系统性略微被高估。然而,必须考虑到实验测量的不确定性与此同级别。
表1使用k-ε修正模型对文丘里型截面进行的数值试验 σ=2.4,Amin=2m/s,ρv/ρl=0.01
(c) k-ω模型:可压缩性效应
如前所述,本文提出的数值模拟采用了单流体物理模型来描述空化现象。根据采用的正压状态方程,在蒸汽/液体混合区域内,声速A = √(∂P/∂ρ)非常低,流体局部高度可压缩。因此,流动在汽化或凝结区域内经常达到超过5的马赫数。实际上,Birch和Eggers表明,马赫数的增加在某些配置中,如混合层,可以显著改变湍流结构。为了分析和评估可压缩性效应对空化流动湍流建模的影响,我们在数值代码中实现了Wilcox提出的k-ω湍流模型。与之前的模型相比,其主要优势在于Wilcox在湍流方程中提出的修正措施,用于模拟密度波动的效应。模型的主要特征如下所述,关于湍流模型的所有详细信息可参见[16]。
基本模型方程
方程类似于k-ε模型中求解的控制方程:所有修正均源自使用湍流耗散率ω代替湍流耗散率ε。湍流黏度表达式和k方程被修改,ε方程被ω方程所取代。
可压缩性效应
通过分析混合层增长率随马赫数的变化,Wilcox观察到,仅通过考虑混合长度的平均密度演变无法很好地模拟这一现象,必须施加显式的附加修正。他得出结论认为,雷诺应力传输方程与可压缩性效应直接相关,并在湍流方程中引入了马赫数的影响。Wilcox提出的修正基于Sarkar等和Zeman的先前研究,他们旨在在k-ε湍流模型中考虑可压缩性效应。他将这些发展应用于k-ω模型,提出了以下方程,将参数βω和βω*定义为湍流马赫数的函数:Mt² = 2k/A²(其中A是局部声速)。
其中指数i表示之前给出的不可压缩模型的值:
结果
图6展示了使用k-ω湍流模型进行的计算。在第一个案例中(图6(a)),模型中未考虑可压缩性的影响。数值结果不理想:与标准的k-ε RNG模型一样,未能正确模拟非稳态空化行为,模拟得到的空化长度比实验结果小50%以上。另一方面,考虑可压缩性效应获得的数值结果与实验结果非常吻合。它们与之前修正的k-ε RNG模型中呈现的结果非常相似:如图6(b)所示,蒸汽云脱落过程得到了正确的模拟。平均空化长度仍然正确,脱落频率几乎保持不变。
图6 使用k-ω模型计算非稳态空化行为。(a)不考虑可压缩效应;(b)考虑可压缩效应。
与使用双光学探针获得的实验数据进行了比较。这种侵入式技术可以测量空化片内部两相结构的局部空隙率和速度。图7展示了四个剖面的时间平均值和标准偏差值。这些比较表明实验结果与数值结果之间总体上有良好的一致性。时间平均空隙率的分布预测得很好,空化片下的回流射流结构在定性上是正确的。波动水平也与测量结果非常接近。x = 0.065 m处的负速度区域对应于空化片的后部,交替受到回流射流的发展和蒸汽云脱落的影响。由于在该配置下平均空化长度略微被低估,因此这两种模型都未能模拟这些负速度。第二个差异可以在0.03 m处观察到,涉及靠近墙壁的空隙率分布;然而,由于重要的波动水平,实验对回流射流中的速度和空隙率的估计在该区域中具有很高的不确定性。因此,实验结果与数值结果沿墙的比较可能并不具有代表性。
在一个文丘里型截面中,应用了四种湍流模型来模拟非稳态空化流动。使用标准的k-ε RNG模型和未考虑可压缩性效应的k-ω模型所获得的结果与实验观测结果不符:这些模型未能再现实验中观察到的蒸汽云脱落行为。此外,[10]中报告的数值测试和本文中提出的结果表明,数值模型引起的整体耗散水平似乎不是其无法模拟非稳态空化行为的原因。
另一方面,修正后的k-ε RNG模型和包含可压缩性效应的k-ω模型能够非常可靠地模拟文丘里的非稳态空化行为。使用修正的k-ε RNG模型获得的令人满意的结果已经通过在其他几何结构中的模拟得到确认,如水翼、翼型级联以及另一个导致更稳定空化片的文丘里型截面。在所有情况下,均正确获得了实验中的整体行为,并且在非稳态配置中准确预测了振荡频率。
本文中使用k-ω模型获得的结果突出了可压缩性效应在湍流模型中的重要作用。实际上,处理可压缩性效应的修正考虑了密度波动:平均方程中出现了一个附加项,增加了湍流耗散。最终效果是在可压缩区域(即蒸汽/液体混合区)内减少了湍流黏度。
为改进k-ε RNG模型而提出的修正也基于在混合区域中降低湍流黏度,这些区域的特点是极低的声速和较大的马赫数。实际上,使用修正后的k-ε RNG模型获得的结果与通过可压缩k-ω模型计算的结果类似。
根据数值计算,流体的可压缩性对湍流结构具有显著影响,必须在模拟非稳态空化流动时予以考虑。