不同湍流模型在水翼前缘空化流场数值计算中的应用研究

2021-11-19郑汉文赖桂桦曾永顺姚志峰

东北电力大学学报 2021年6期

郑汉文，刘婧，赖桂桦，曾永顺，姚志峰，3

(1.中国农业大学水利与土木工程学院，北京 100083；2.中国航发商用航空发动机有限责任公司，上海 200241；3.中国农业大学北京市供水管网系统安全与节能工程技术研究中心，北京 100083)

伴随能源需求量的增加，大量扬程高、尺寸大、转速高的水力机械被投入能源生产，机组运行的稳定性问题日益突出.在水力机械工作的过程中，当液压降低至汽化压力时，叶片前缘可能出现附着型空化区域[1]，且该区域内空泡也可能脱离前缘向下游延伸，这种现象称为“前缘空化”[2].前缘空化是水力机械中主要空化类型，空化可侵蚀结构表面，引起结构高幅值振动[3].前缘空化一般发生在层流边界层分离后，该区域压力分布与空泡长度有密切关系[4].前缘空化末端汽液混合，空泡与水流发生强烈相互作用[2]，表现出高度非定常特征[5].

大量学者通过实验研究了前缘空化的发展形态、流场特性以及与水翼结构相互作用[6-8].Kawanami[6]等研究了绕水翼云状空泡的脱落机制，认为空泡断裂与脱落由空泡尾部反向射流引起.Wu[7-8]等探究了弹性水翼空化，发现空化发生会显著增加弹性水翼的形变，同时弹性水翼动态特性也会改变空化发展的进程.

由于实验在测量空化流场方面存在局限，数值模拟逐渐成为空化研究的重要手段[9].大量学者对数值模拟算法进行了研究，探讨了空化数值模拟中湍流模型与空化模型的改进方法.Geng[10]等通过改变空化模型中凝结系数和蒸发系数参数，对比分析水翼空化流场特性后得出当凝结系数为0.08至0.1、蒸发系数为300至500时实验与数值模拟较为吻合.

目前，求解湍流的方法分为直接求解和间接求解两大类.直接求解方法(DNS)可以得到湍流流场的精确信息，但对计算资源的要求非常高，在实际工程中往往不采用直接求解方法.间接求解方法可分为三大类：第一类是系综平均法，主要代表是雷诺平均法(RANS)，RANS可以计算高雷诺数的复杂运动，但计算的结果一般都是时均值，不能反应湍流的细节信息；第二类是空间筛滤法，主要代表为大涡模拟法(LES)，LES基于湍动能传输机制，直接计算大尺度涡流的运动，小尺度涡流的运动对大尺度涡流的影响通过建模体现出来，可以得到比RANS方法更多的动态信息；第三类为混合方法，主要代表为分离涡方法(DES)，它兼有RANS和LES的特点[11].

为了更加明确空化形态定量分析方法，在采用相同空化模型的条件下，本文对比了RANS、LES和DES三种湍流模化方法在空泡形态预测上的差异.通过实验数据，明确了空化模型中蒸发系数和凝结系数参数，定量分析了空泡形态与实验结果的差异，为工程上预测前缘空化流场提供指导.

1 湍流计算方法

1.1 雷诺平均(RANS)法

雷诺平均法(RANS)广泛使用在工程中，其原理是采用湍流统计理论，将非稳态N-S方程作时间平均，求解工程中需要的时均量[12].该方法将流场中的变量分解为平均值和脉动值，并在此分解的基础上求解均化后的运输方程[13].雷诺平均后的N-S方程[14]：

(1)

(2)

(3)

(4)

公式中：ut为湍动粘度；k为湍流动能；δij为判定符号，当i=j时，δij=1，当i≠j时，δij=0.

根据湍动粘度ut求解方程的数量，可将涡粘模型分为零方程模型、一方程模型、两方程模型.目前在工程应用中两方程模型使用最为广泛，主要有标准k-ε模型、RNGk-ε模型、k-ω模型和SSTk-ω模型等[13].

在k-ω两方程模型中，湍动粘度ut的数值与湍动能k和比耗散率ω相关，其关系如下式所示

(5)

湍动能k方程和比耗散率ω方程分别为

(6)

(7)

公式中：β′=0.09，α=5/9，β=0.075，σk=2，σw=2为常数；密度ρ、速度矢量ui可以通过求解N-S方程得到；pkb及pωb为湍动能浮力生成项；pk为湍动能生成项[16]，其计算公式

(8)

基于k-ω的SST湍流模型将湍流切应力的运输考虑在内，并对逆压梯度下流动分离的起始位置和数量有精确的预测.该湍流模型通过公式(9)定义涡粘系数，以得到湍流切应力为

(9)

公式中：νt=μt/ρ；S为应变率的不变量度；F2为混合函数，用以将该式的运用限制在边界层内，其表达式为

(10)

(11)

湍流模型的正确使用依赖混合函数，其表达式取决于计算点离最近的固体表面的距离以及流场流动参数，式中y即为距固体表面最近距离.

1.2 大涡模拟(LES)方法

LES湍流计算方法的基本原理是通过滤波函数将瞬时湍流运动分分解为大尺度和小尺度两种类型.对大尺度运动采用直接求解的方法进行计算，对过滤的小尺度运动建立亚格子尺度(Subgrid-Scale，简称SGS)应力模型进行计算，以体现小尺度运动对大尺度运动的影响.LES模型未进行时均化处理，只对空间进行滤波，计算精度上相对于RANS模型具有一定优势，同时LES模型较直接求解可节省大量计算资源.

(12)

(13)

公式中：D为流体域；G为滤波函数；V为控制体；x′为实际流体域中的坐标值；x为滤波后大尺度空间中的坐标值.

不可压缩流体滤波后的连续性方程及动量方程如公式(14)和公式(15)所示；其中τij为亚格子尺度应力，用以体现小尺度涡的运动对大尺度涡的影响，其表达公式为

(14)

(15)

(16)

在LES湍流模型求解过程中，大尺度涡流被直接求解，小尺度涡流则采用亚格子模型.CFX中采用涡流粘度方法，通过下式将亚格子尺度应力与大尺度应变率张量相关联.

(17)

(18)

(19)

公式中：Ls为亚格子尺度混合长度，计算公式为Ls=CSΔ；CS为Smagorinsky常数，该参数的取值对模型有较大影响，当取值接近0.1时能够较好模拟大多数流动；Δ为当地网格尺度.

本文所采用的LES WALE(Wall-adapted local eddy viscosity)模型，相对于Smagorinsky模型在层流剪切层的求解更为准确，其涡粘系数μt采用公式(20)求解.

(20)

(21)

1.3 分离涡模拟(DES)方法

求解高雷诺数的边界层流动，在使用LES湍流计算方法成本过高，使用RANS湍流计算方法不能满足精度要求的情况下，工程中常采用DES湍流计算方法.DES湍流计算方法将RANS和LES湍流计算方法的特点结合在同一个混合方程中，即在流动边界层内未发生分离或流动分离较为轻微时使用RANS求解，而在流动分离剧烈的区域使用LES求解.在SST-DES模型中，当RANS模型计算得到的湍流长度Lt大于该区域内网格尺度时，计算由SST模型转换为LES模型[15].此时，湍动能耗散项方程的长度尺度变为该区域的网格尺度Δ

ε=β*kω=k3/2/Lt→k3/2/(CDESΔ)for(CDESΔ

(22)

Δ=max(Δi)，

(23)

(24)

k方程变化为

(25)

公式中：CDES=0.61.

2 计算方案

2.1 计算模型

以二维NACA 0009钝型尾部形状水翼作为本文计算对象，研究水翼前缘空化流场.翼型弦长L为100 mm，尾部厚度h为3.22 mm，如图1所示.

图1 二维NACA 0009钝型尾部形状水翼模型

流体计算域，如图2所示.水翼前缘距测试段进口2.5L；水翼尾部距测试段出口4L；水翼安放位置距测试段顶部0.7L，攻角为2.5°.在距水翼尾部中点延长线上距水翼尾部10 mm处设置监测点，并定义x方向为流体流动方向，y方向为垂直于流动的方向，得到如图1的坐标系.

图2 二维计算域

2.2 网格划分

空化计算中在流速20 m/s、攻角2.5°、空化数σ=0.81的条件下不同湍流模型的网格划分方案和边界条件，如表1所示.为消除计算域网格划分方式对计算结果的影响，本文以水翼所受升力级及压力为关键参数，采用基于理查德森外推加速法的网格收敛性指数[17]对上述网格划分方案进行可靠性分析.结果表明，RANS型、DES湍流模型、LES湍流模型的收敛性指数分别在13%、1%、23%以内，均满足计算要求.

表1 不同湍流计算方法网格数及边界条件(流速20 m/s、攻角2.5°、空化数0.81)

2.3 计算设置及边界条件

采用ZGB空化模型进行空化数值计算，选择Water at 25 ℃及Water Vapour at 25 ℃为测试段中水及水蒸气介质，其中，水的密度ρ=997 kg/m3，水蒸气密度为ρ=0.023 08 kg/m3.

RANS湍流计算方法瞬态项采用二阶欧拉后差分格式，对流项及湍流模型相关量均采用高阶精度求解；DES湍流计算方法瞬态项采用二阶欧拉后差分格式，对流项采用中心差分格式，湍流模型相关量采用高阶精度求解；LES湍流计算方法瞬态项采用二阶欧拉后差分格式，对流项采用高阶精度求解.

计算域进口边界条件为速度进口，来流速度20 m/s；出口边界条件为压力出口，压力值与Dupont[18]实验中测试段空化数σ=0.81保持一致；计算域中上下壁面及水翼表面采用忽略粗糙度的无滑移壁面.空化计算前先进行流场的定常计算，收敛后，在定常计算的基础上引入空化模型进行空化流场的非定常计算.非定常计算设置每个时间步内最大迭代次数为50次，收敛残差标准为1×10-4.

3.4 空化模型凝结系数、蒸发系数调整

结合实验数据进行大量计算后得出，对于RANS湍流计算方法，当凝结系数为0.15、蒸发系数为140时，水翼表面压力系数与实验值较为吻合；对于DES湍流计算方法，当凝结系数为0.3，蒸发系数为140时，水翼表面压力系数与实验值较为吻合；对于LES湍流计算方法，当凝结系数为0.3，蒸发系数为30时，水翼表面压力系数与实验值较为吻合.后文将采用这组参数进行空化模拟.

3 计算结果与分析

3.1 前缘空化形态定量分析

(26)

(27)

Dupont[18]实验中水翼NACA0009在2°攻角、20 m/s流速、空化数σ=0.81时的空化区域，如图3所示.其中空泡长度lcavity=0.029 4 m，空泡高度hcavity=0.002 6 m.下文将以此数据作为衡量对象分析不同湍流模型对前缘空化形态的影响.

图3 水翼前缘空化区域长度及高度定义

分别取三个湍流模型的一个运动周期T，截取该周期内6个时刻对应的空化区域图像算出空泡长度和空泡高度，并通过公式(26)、公式(27)求得该周期内的lcavity和hcavity，结果如表2所示.计算结果显示，与RANS和DES湍流模型相比，LES湍流模型计算所得空泡长度lcavity和空泡高度hcavity与实验值更为吻合；对于RANS和DES湍流模型，空泡长度的最大相对误差为16.8%，空泡高度的最大相对误差为26.7%.

表2 不同湍流计算方法预测空泡长度及高度(流速20 m/s、攻角2.5°、空化数0.81)

在流速20 m/s、攻角2.5°、空化数0.81条件下，一个周期内3个时刻RANS湍流模型、DES湍流模型和LES湍流模型的前缘空化流场，如图4～图6所示.对比分析三种湍流模型的空化区域和绕流流场发现：(1)RANS湍流计算方法模拟的空泡形态不随时间改变，在空泡尾部x/L在0.22～0.3范围内存在明显的回流区域；(2)DES湍流计算方法模拟的空泡在其尾部存在脱落区域，空泡脱落后不会再向流动方向移动，并且x/L在0.24～0.3这个区域内，空泡的长度和高度会随时间有一个较小的改变；(3)LES湍流计算方法可以模拟空泡从初生经过生长到脱落的整个过程，在这个过程中空泡的长度和高度会随着时间有一个较大的变化，x/L在0.08～0.31范围内均存在回流现象并且回流区域的位置会随空泡形态的改变而改变.

图4 RANS湍流计算方法空泡形态及绕流流场

图5 DES湍流计算方法空泡形态及绕流流场

图6 LES湍流计算方法空泡形态及绕流流场

3.2 近壁区速度分布

为探究不同湍流模型对前缘空化近壁区速度的影响，在水翼相对弦长x/L=0.1、0.4、0.7、0.99位置处作垂直于x轴的采样线，每条采样线上均布100个采样点，测取3个流动周期内每个采样点对应的x方向速度和y方向速度，并将所测速度值求平均得到每点的平均速度Caverage.将速度无量纲化作为横坐标，C=Caverage/Cinlet，其中Cinlet为计算域进口处速度；将y坐标与水翼弦长L之比作为纵坐标.将三种湍流模型对应的无量纲速度和实验值进行比较[18]，如图7所示.图7(a)给出了采样线段处的空泡的平均高度.

图7 不同湍流模型水翼边界层速度(流速20 m/s、攻角2.5°、空化数0.81)

对比x方向速度的实验值和计算值发现：在x/L=0.4和x/L=0.7时三种湍流模型给出的预测值均与实验值相差较大；在y/L>0.1时，三种湍流计算方法给出的预测值均与实验值吻合的较好，但在y/L<0.1区间内，预测值均显著高于实验值.

对比y方向上速度的实验值和计算值发现：在x/L=0.4和x/L=0.7时，实验值中存在沿y轴正向的速度分量，即存在速度大于0的区域，但从图中可以看出RANS、DES、LES湍流模型均未准确的预测到这一现象；在y/L>0.1时，三种湍流计算方法给出的预测值与实验值吻合得较好.

当采样点与壁面间的法向距离超过0.1，三种湍流模型计算所得x方向速度和y方向速度均与实验值较为吻合，即不同湍流模型在远离壁面的区域对流场的模拟较为准确.由于本文空化计算的假设是均相流，计算结果反映的是均相速度，与真实值存在差异，因此在近壁区，RANS、DES和LES湍流模型计算的边界层速度与实验结果存在较大差异，不能准确反映近壁区速度分布.

三种湍流模型下同一流动周期内水翼表面时均速度分布云图，如图8所示.显示了水翼前缘空化流场中速度的分布规律.三种湍流模型对应前缘空化流场中的速度分布规律基本相同：空化区域均存在于水翼前缘且该区域的流速显著高于该流场中的其它区域；水翼的前端和尾部均具有速度较低的的区域：前端的低速区成因是来流的撞击，尾部的低速区成因是漩涡脱落；水翼的尾部存在因漩涡脱落形成的速度为负的区域.

图8 不同湍流模型水翼表面时均速度分布云图(流速20 m/s、攻角2.5°、空化数0.81)

三种湍流模型对应空化流场速度分布存在差异：在水翼的上表面，RANS和DES湍流计算方法速度为负的区域集中，而LES湍流计算方法速度为负的区域分散；RANS和DES湍流计算方法回流区域固定，LES湍流计算方法回流区域随时间变化，说明LES湍流计算方法的出的空泡脱落位置随时间变化，可以反映出空化的非定常特性.