带有剪切应力输运性质的k-ε两方程湍流模型构造与应用

2019-08-21陈利丽

空气动力学学报 2019年3期

王伟,张扬,陈利丽

(1.航空工业第一飞机设计研究院,西安 710089;2.西安交通大学机械结构强度与振动国家重点实验室,西安 710049)

0 引言

在空气动力学所涉及的流动中，逆压梯度引起的分离流动是极为常见的流动现象之一，尤其是在航空器的气动设计领域，这一流动现象多带来气动性能损失、结构疲劳和操作性稳定性变差,会严重威胁航空器的飞行性能和安全。虽然直接数值模拟方法(DNS)和大涡数值模拟方法(LES)能够描述湍流场的结构和时间历程，但是对于工业界中广泛的气动力评估，仍然需要依赖于高效鲁棒的雷诺平均方法(RANS)。然而，普遍存在的非平衡效应是RANS方法所需要面对的最大求解障碍之一。RANS方法中对于湍流脉动的求解仍然基于雷诺分解假设，即关于雷诺应力的合理建模。早期的纯代数封闭模型，如Baldwin-Lomax(BL)模型[1]，由于缺乏对于湍流变量的历史效应反馈，在求解带有分离现象的流动时精度较差。但是，BL模型提供了一种雷诺应力的简化封闭思路——即完成湍流速度尺度、长度尺度和时间尺度中的至少两个量,即可完成在线性涡黏假设下的雷诺应力封闭。为了求解这三个变量，学者们发展了“N”方程模型，其中N代表求解基本变量所需的偏微分方程数目。例如，由Baldwin等发展的一方程模型(BB模型)[2]，以没有经过壁面衰减的涡黏系数为输运变量构造了一个输运方程，虽然表面上看似乎在效率与精度间达到了最好的平衡，并且引入了时间效应，然而经典的BB模型是由k-ε两方程模型推导而来，数项扩散相关项被省略，这必然引起湍流求解精度的下降。Menter等[3]的研究表明BB模型在剪切层边缘的模拟精度较差，Rahman等[4]认为BB模型的缺陷是由于破坏项的病态与对来流的极为敏感造成。Spalart等[5]基于经验和类涡黏性系数提出的Spalart-Allmaras模型(SA模型)是目前最受欢迎的一方程模型，其构造理念采用“堆积木”的形式完成，通过简单流动一步步进行标定，最终完成对于复杂流动的求解，壁面效应被融入了等价涡黏系数的比例系数中。SA模型在飞行器流动中的应用极为广泛。但是，Menter[6]认为SA模型在预测逆压梯度流动中精度欠佳，尤其是在回流问题中容易将回流速度过高估计。Coakley提出了涡黏系数的混合求解模型[7]，Menter[3,6]将这个理念发展为基于k-ω两方程的剪切应力输运湍流模型，以Bradshaw假设[8]为限制，改善了在逆压梯度和分离流动出现时，传统湍流模型将产生项做过大估计的缺陷，并且在该假设起效时，整个湍流模型实质上变成了剪切应力的输运方程。在此之后，有学者[9-11]将这种理念发展到一方程湍流模型，并在复杂流动的求解中得到了很好的应用。然而，在k-ε两方程模型中，还未有剪切应力概念的引入，本文将原始的Bradshaw假设采用直接数值模拟(DNS)数据重新标定，基于结构系综动力学[12]理论给出具有物理意义的标定方程，以此消除原始SST方程中的混合函数项。

另一方面，如图1所示，基于DNS数据的反求以及标准k-ε湍流模型的对比发现可知，近壁面附近的涡黏系数与DNS对应值差距较大。其原因在于通过摄动法分析后，采用k-ε两方程模型中的涡黏系数求解方法(Cμk2/ε)与涡黏性系数(μt)渐进解不一致，需要乘以一定的系数来满足固壁面的影响来得到正确的边界解。因此衍生了一种新的湍流模型子集即低雷诺数湍流模型(Low Reynolds Number,LRN)。该种模型的设计初衷并非针对特征雷诺数很低的流动，而是为了更正如图1所示的现象，修正当地雷诺数较低时无法产生正确的边界层渐进行为而提出。在近壁面区域通过逐步衰减系数来完成均衡态假设的近壁效应修正，该修正函数被称为壁面衰减函数(WDF)。因此低雷诺数湍流模型中的“雷诺数”指代的并非是模拟流动的特征雷诺数。LRN湍流模型中，Jones等[13]提出了第一种WDF函数并得以应用，在近十年来，几十种湍流模型相继被提出，随之带来的便是数种WDF函数以及额外的源项引入，但最终目的都是为了使湍流模型本身产生合理的近壁渐进行为。Patel[14]评估了9种LRN模型后，通过湍流边界层算例展示这些模型的精度，虽然所有模型都能够完美满足湍流对数率，然而，在湍动能的预测上表现的差距非常大，其根本原因来自于WDF函数的构造方法。常用来构造的基本元素是两种湍流雷诺数(Reτ和Rek)和无量纲的壁面距离y+，然而，y+引入的壁面摩擦速度在流动的分离点和再附点处为0，容易引起WDF函数的病态。另外，如图1所示，WDF的合理作用范围应该不高于对数层，但是，多数LRN模型在构造WDF函数时，过多关注WDF的衰减而忽略了WDF函数的峰值变化。Rodi等人[15]的研究清晰地表现出WDF的峰值是随着雷诺数的变化而变化，因此，本文另一个工作即为提出一个峰值可变的WDF函数形式。

新模型通过充分发展的湍流平板、NACA4412翼型、二维鼓包以及DLR-F6翼身组合体的计算来验证计算精度。

图1 基于DNS和标准k-ε模型求得的涡黏性系数μtFig.1 Eddy-viscosity coefficient got by DNS and standard k-ε turbulence model

1 湍流模型演化思路

1.1 k-ε两方程湍流模型

忽略掉浮力项的影响，标准不可压缩k-ε湍流模型的形式为：

(1)

(2)

其中:ρ为密度，Uj为速度矢量，k为湍动能，μ为层流黏性系数，μt为湍流黏性系数，Pk为湍动能产生项，ε为湍动能耗散率，t为时间变量，xj为空间张量，Cε1、Cε2为常数。

Pk=μtΩ2

(3)

1.2 WDF函数的构造

通常，对于二维边界层流动，将近壁面附近速度脉动量(即u′、v′和w′ )采用小扰动方式泰勒展开，可得：

(4)

由于受无滑移边界条件的约束，式(4)中的速度脉动右端均无常数项，同时，联合连续方程可得b1=0。通过式(4)，可以得到近壁附近湍动能和湍动能耗散率的阶数为：

(5)

式(5)表明，在近壁面附近保持Cμk2/ε的湍流涡黏性系数显然是不合理的，需要一定的降阶函数使μt保持正确的边界渐进性。此类函数即为WDF函数，WDF的精确表达式为：

(6)

Zhang等人[16]结合DNS数据，将三种常见湍流模型(Abidk-ε[17]、Yang-Shihk-ε[18]和Launder-Sharmak-ε[19])中的WDF函数通过式(4)求解出来。如图2，随着雷诺数的变化，三种模型的WDF函数并未随着一同改变，而是被机械地限制在了值为1附近。通过表1可以看出，三种模型的WDF函数在自变量趋于无穷时，峰值均为1，然而从DNS数据反求出的WDF峰值是随着雷诺数不同而变化的。

参考表1给出三种对比模型的WDF形式，从数学性质上可以看出，三个WDF的极值随着自变量的增加均为1。而图2的DNS数据反推表明，这样的极值是不合理的。

根据观察，Zhang等人[16]以因子Rek/Reτ来表征雷诺数的变化，通过构造两个不同的函数fμ1和fμ2，基于此形成新的WDF函数：

(5)

式(5)的构造理念在于——fμ1是将现有的WDF进行联合，以此达到近壁面附近更精确的衰减作用，通过fμ2控制峰值的增减。图3为不同雷诺数槽道中Rek/Reτ的表达。如图4所示，新构造的WDF函数的峰值与雷诺数变化具有一种自适应关联,其中fnc代表非均衡湍流识别函数。

(a)Change factor with the change of y+

(b)New WDF’s structure function fμ2图3 新构造的WDF带入DNS数据中的变化Fig.3 Changes of newly constructed WDF into DNS data

图4 新WDF函数在不同雷诺数下的变化Fig.4 Variation of new WDF function at different Reynolds numbers

1.3 改进的湍涡黏性系数计算方式

一般情况下，在k-ε两方程模型中，湍涡黏性系数的确定方式为：

μt=Cμfμk2/ε

(6)

式(6)是从均匀各向同性平衡湍流中推导而来，因此对于带有分离和逆压梯度的流场适应性不强。Menter[6]的研究表明，在逆压梯度的流动中，湍动能产生项可以数倍大于其破坏项，导致雷诺应力的过大估计，因此引入了Bradshaw假设进行平衡。然而，Brashaw将湍动能与雷诺应力间的比值设定为常值(0.31)，而Nagano[11]和文晓庆等[20]的研究均表明该常数是一个可变值。张扬[21]等根据已有的DNS数据，反推出Bradshaw常数，并以佘振苏[12]等人提出的结构系综动力学理论对新的常数进行标定。本文采用张扬等[21]的研究成果，将这个常数定义为：

(7)

基于式(7)，湍涡黏性系数计算公式为：

(8)

式(8)考虑了湍动能产生于耗散的平衡系数NP2D(式中fne代表非均衡湍流识别函数)。该系数的引入主要为了在产生项过大时将涡黏性系数的确定方式指向剪切应力输运模型，即将该系数与湍动能直接联系起来。

2 算例验证

2.1 湍流平板

采用湍流平板对湍流模型进行校验是发展新模型的第一步。本文采用Wieghardt等[22]发布的实验数据进行比对。计算的工况为：Ma=0.2，平板板长设定为1，因此其特征雷诺数Re=6×106，采用3套不同粒度的网格进行模拟，网格规模如表2所示。从图5给出了Rex=2.6×106处的速度型。可以看出，新模型在不同的网格下均能给出准确的湍流壁面率，证明新模型能够给出湍流边界层正确的速度结构。

表2 平板计算使用网格Table 2 Computational grids

2.2 DLR-F15三段翼型

二维多段翼型 DLR-F15是德国宇航院在 2005 年试验的一个三段翼型，该翼型是低噪声增升装置项目[23](Low Noise Exposing Integrated Design for Start and Approach ，LEISA)的研究对象。DLR 在荷兰 DNW-KKK 风洞对 DLR-F15 多段翼型进行了多个马赫数和雷诺数的测试，并且公开了部分数据。试验状态为：Ma=0.15，Re=2.094×106(参考弦长为c=0.6 m)，T=525R(R为兰氏温度单位)，α=6°。

(a)Coarse

(b)Medium

(c)Fine

图5 湍流平板的速度型
Fig.5 Velocity profiles of turbulent plates at two positions

从图6可以看出，新模型对于缝翼、主翼以及襟翼的峰值模拟精确，而SST模型过高地预测了前缘缝翼和主翼的峰值。

图6 DLF-F15多段翼压力分布Fig.6 Pressure distribution of DLR-F15

2.3 NACA4412翼型

NACA4412翼型[24]在特定的迎角下，后缘会出现一个稳定的分离泡，对于分离泡的起始点、分离区域的速度型以及再附点的预测一直都是一个具有挑战性的工作。计算的工况为：Ma=0.09，基于单位弦长的特征雷诺数Re=1.52×106，迎角α=13.79°。为了简化起见，采用一套y+小于1的网格进行计算，网格总量为57 921。远场边界约为翼型弦长的100倍。挑选6个站位进行速度型取样，分别为x/c=0.6753、x/c=0.7308、x/c=0.7863、x/c=0.8418、x/c=0.8973和x/c=0.9528。

从速度型上(图7)可以看出，新模型由于引入了对于分离流动这种非平衡效应的感知公式，计算所得速度型与试验值极为接近，而其他模型不同程度地低估了分离区的回流速度。

图7 NACA4412翼型不同站位的速度型Fig.7 Velocity type of NACA4412 airfoil at different locations

2.4 二维鼓包的计算

基于NASA官网发布的分离预测标模(二维鼓包)[25]，计算了其表面压力分布和表面摩擦力系数分布。

从压力分布上可见(图8)，Abid模型计算的压力最大值与新模型的计算结果近似，SST模型对于压力最大值的模拟偏低。对于负压峰值的模拟，新模型和Abid模型结果类似，SST模型模拟的负压峰值偏大。

表面摩擦力曲线(图9)显示新模型对于摩擦力峰值的预测好于其他两个模型，但是对于再附点的预测偏后。其他两个模型也未能很好地贴合。

2.5 DLR-F6翼身组合体

DPW会议考量的是从空客某构型中简化而来的一个机翼-机身-发动机短舱构型(DLR-F6)。风洞试验由法国 ONERA S2MA 风洞承担[26]。本文采用其中的翼身组合体(不含发动机短舱)模型进行计算，以考察湍流模型对于复杂构型的计算精度。计算状态和特征参数为：Ma=0.75，Re=3×106，T=520R，平均气动弦长c=0.1412 m，参考面积S=0.07265 m2，半展长B=0.5877 m，计算时定升力系数CL=0.5。图10为不同站位的压力分布计算结果，可以看出，新模型在三维构型上与其它三个模型结果基本相近。