基于νt 尺度方程的雷诺应力模型初步研究*

2022-08-28陈彦君王圣业符翔刘伟

物理学报 2022年16期

陈彦君王圣业符翔刘伟

(国防科技大学空天科学学院,长沙 410073)

雷诺应力模型一直是湍流模式理论研究的前沿和难点,而提高数值鲁棒性是其广泛开展工程应用的关键.借鉴经典的 k-kL湍流模型,本文构造了一种新的 νt 尺度方程,并将其用于耦合SSG/LRR 模式从而形成SSG/LRR-νt 雷诺应力模型.通过零压力梯度湍流平板边界层、翼型尾迹流、超声速方腔流和NACA0012 翼型45°迎角分离流动4 个标准算例对新模型进行了验证与确认.同时,为了测试模型的数值鲁棒性,采用高精度数值格式对模型方程进行了离散求解,并与SA 涡粘模型和SSG/LRR-ω 雷诺应力模型进行对比.结果表明: νt 尺度方程在黏性壁面边界严格为零,相比传统的 ω 尺度,具有更好的数值鲁棒性,从而可实现新模型与高精度数值格式的匹配并获得更好的网格收敛效率;新模型具备雷诺应力模型的传统优势,可对拐角流动进行很好的模拟;具备尺度自适应能力,对于非定常分离流动的模拟存在一定的潜力.

1 引言

湍流是经典物理学中著名的难题,数百年来人们一直致力于对它的研究.其中,湍流模式理论是人们研究湍流运动规律得到的主要成果之一,并且依托飞速发展的计算机技术,成为目前解决实际湍流问题的主要手段.无论是雷诺平均模拟(RANS),大涡数值模拟(LES)还是混合RANS/LES 模拟,核心均是对雷诺应力张量(或亚格子应力张量)进行封闭建模.这包括基于Boussinesq 假设的线性涡粘模式,直接建立应力输运方程的雷诺应力模式(或亚格子应力模式)以及其他非线性模式等[1].

本文主要研究RANS 方法中的雷诺应力模式(RSM),它由周培源先生[2,3]首次建立起一般方程框架.RSM相比目前广泛应用的线性涡黏模式,在诸多方面具备优势,如曲壁面流动、旋转/涡流动、拐角流动等[1,4,5].但也存在明显不足,包括计算花费和数值稳定性两个主要方面.前者属于固有问题,是需要离散求解更多的偏微分方程所导致的.但当前随着计算机硬件水平的巨大提升,其超出涡黏模式的计算代价已经可以接受,且远未到达LES 的需求水平.而对于后者,包括美国宇航局(NASA)[6]、德国宇航局(DLR)[7]等均指出需要持续开展研究投入.

无论何种形式的RSM,方程中往往仍包含一个未知量,即各向同性耗散率ε.在物理上实质是需要提供一个标量的湍流长度尺度或时间尺度.多年实践表明,尺度提供方程对RSM 的数值稳定性有非常大的影响.例如,早期的RSM 往往耦合ε尺度方程求解[8].但与在k-ε模型中类似,ε方程存在两个问题[9]: 一是缺乏自然边界条件;二是在壁面附近需要处理高阶关联,易造成数值刚性问题.这些问题使得早期的RSM 应用效果很不理想.Wilcox[1]发展的ω方程是在尺度建模中里程碑式的工作,它避免了ε方程在近壁附近的数值刚性问题.其后Menter[10]又进一步发展了ω尺度方程,通过设计过渡函数,使得ω尺度方程在远场也具良好表现(过渡到类似ε模型,减弱来流湍流度敏感性).2005年,Eisfeld和Brodersen[7]将Menter 的ω尺度方程用于RSM,提出了SSG/LRR-ω模型.该模型经过10 多年的研究和改进,在很多工程问题中取得了良好的模拟效果[11-14].

尽管如此,ω方程本身在理论上仍然存在一些问题.Togiti和Eisfeld[15]指出ω在壁面附近存在奇性且缺乏自然边界条件.因此将g方程[9](g1/在壁面边界为0)用于RSM(SSG/LRR-g模型),并表明可降低对近壁区网格分辨率的依赖.其后,Eisfeld和Togiti等[14]又将Ilinca和Pelletier[16]发展的 l og(ω)尺度变换用于RSM(SSG/LRR-log(ω)模型),同样获得了明显的鲁棒性提升.Abdol-Hamid[17]也提出了将kL尺度方程用于RSM 的思路,但并未深入研究.舒博文等[18]对SSG/LRR-g模型在多个典型航空问题中进行了应用研究,指出了该模型预测分离问题的优秀能力.本文作者借鉴 Rotta[19],Menter和Egorov[20]以及 Abdol-Hamid等[21]发展的k-kL模型,推导了新的νt尺度方程(νt),并用于RSM.该尺度方程从变量形式上看,类似Spalart-Allmaras (SA)模型[22],在壁面边界为0 具有更好的数值潜力.同时相比g方程和l og(ω)方程,具备尺度自适应能力,还可用于非定常分离问题的模拟.

2 湍流模型及数值方法

2.1 νt 尺度方程

方程右端依次为生成项、耗散项、耗散修正项和扩散项.“—”代表一般Reynolds 平均量;“∧”代表Favre 平均量.dw为网格到壁面的距离.

LvK为von Karman 长度尺度,通过下式得到:

同时,为防止出现非物理现象,需要对LvK进行限制: 0.1Lt＜LvK＜1.3κdwfp.其中,限制函数fp为

耗散修正项中的经验函数fε由下式得到:

其余经验系数通过丰富的湍流算例进行标定.本文分别取κ0.41,Cp10.775,Cp21.47,Cε和σν2/3.最后需要指出,νt在壁面边界为0;在远场或入口处根据实际的湍流黏性比赋值,缺省值为 0.009ν.

2.2 SSG/LRR-νt 雷诺应力模型

对NS 方程进行Reynolds 平均后会出现雷诺应力项.传统线性涡黏模型,如SA 模型,基于Boussinesq 假设对雷诺应力项封闭;而雷诺应力模型则是直接建立雷诺应力输运方程:

雷诺应力模型中,最关键的是再分配项.本文采用Eisfeld和Brodersen[7]发展的SSG/LRR 混合模型,即通过过渡函数F实现在近壁区使用Launder-Reece-Rodi 模型[8],在远离壁面过渡到Speziale-Sarkar-Gatski 模型[23].

表1 再分配项系数Table 1.Coefficients in redistribution item.

过渡函数根据Menter[10]在SST 模型中提出的F2得到:同时,扩散项系数也通过该函数得到:σR0.5F+0.44(1-F)/(3Cµ).

2.3 高精度数值方法

本文采用团队自研的高精度CFD 软件[4,24,25].该软件主要基于单元中心型有限差分格式,包括二阶精度MUSCL 格式、五阶和七阶WCNS 系列高阶精度格式[26]等.需要强调的是,本文对NS 方程和湍流模型方程求解采用松耦合模式,但针对两者的空间离散精度是一致的.换句话说,湍流模型将同样采用高阶精度离散,以此验证新模型的数值鲁棒性.另一方面,本文采用邓小刚等[27]提出的对称守恒网格导数算法(SCMM)计算网格导数,以降低网格变换引入的数值误差.

3 计算结果

3.1 零压力梯度湍流平板边界层

零压力梯度湍流平板边界层是湍流模型研究中最基础的验证算例.本节参考NASA 湍流模型资源网站[28]中推荐的流动条件: 入口马赫数Mainlet0.2,每米雷诺数Rem5×106.通过控制第一层网格距壁面的高度,生成了5 套网格.网格轮廓见图1(a),其中平板长度为2 m.

图1 平板边界层网格收敛性分析 (a)平板网格示意图;(b)阻力系数结果Fig.1.Convergence analysis on plate boundary layer meshs: (a)Sketch of plate mesh;(b)drag coefficient results.

图1(b)给出了SSG/LRR-νt模型在5 套网格上得到的阻力系数结果.该模型分别用二阶MUSCL 格式和七阶WCNS 格式进行离散求解.同时添加了二阶和七阶精度下的SA 模型和二阶精度下的SSG/LRR-ω模型进行了对比.值得一提的是,SSG/LRR-ω模型由于ω在壁面附近存在奇性,无法搭配低耗散的七阶格式进行稳定的计算.而νt尺度在壁面边界严格为零,使得整个模型具有更好的鲁棒性.再搭配高精度数值格式,可实现类似SA 模型的优秀网格收敛特性.

图2 给出了SSG/LRR-νt模型结合七阶WCNS格式在0.37 网格上的计算结果.在x0.97处速度型与经典的Coles 公式符合;整个平板上摩阻分布与经验公式符合.证明了SSG/LRR-νt模型能够对最基础的湍流平板边界层进行有效模拟.

图2 SSG/LRR-νt 模型在=0.37 网格上的结果 (a)x=0.97 处速度型;(b)摩擦阻力分布Fig.2.Results of SSG/LRR-νt model on grid of=0.37 : (a)u-velocity profile at x=0.97;(b)friction drag coefficient along the plate.

图3 表明,相同网格时,七阶精度格式的计算花费约为二阶格式的1.8 倍,但随着网格尺度的减小,高精度格式的计算误差下降更快.所以,在达到相同误差水平的条件下,七阶格式的花费比二阶格式小,即效率更高.另一方面,相同网格时,SSG/LRR-νt模型的计算花费约是SA 模型的1.6 倍.其花费的增加比率小于方程数的增加(前者共12 个偏微分方程;后者6 个),说明SSG/LRR-νt模型的求解效率是良好的.

图3 x=0.97 处摩擦阻力误差与网格尺度和计算时间的关系 (a)误差与网格尺度;(b)误差与计算时间Fig.3.Relationship between friction drag error at x=0.97 and grid scale as well as CPU time: (a)Error vs.grid scale;(b)error vs.CPU time.

3.2 翼型尾迹流

翼型尾迹流同样来自NASA 湍流模型资源网站,用来考核湍流模型在自由剪切流动中对雷诺应力分量的模拟准确度.流动条件为: 来流马赫数Ma∞0.088,基于弦长的雷诺数Rec1.2×106.网格采用该网站提供的粗网格(2 81×49),如图4 所示.

图4 翼型尾迹流网格及切面Fig.4.Airfoil wake mesh and slices.

图5 给出了图4 所示切面上雷诺切应力分布.与平板算例类似,采用七阶WCNS 离散,仅在粗网格上即可实现较好的模拟效果.证明了高精度离散对于湍流量的预测同样有意义.SSG/LRR-νt模型由于更好的数值鲁棒性,能够成功与高精度格式结合进行求解.而SSG/LRR-ω模型却无法稳定求解.

图5 翼型尾迹流雷诺切应力分布Fig.5.Reynolds shear stress distribution on airfoil wake case.

3.3 超声速方腔流

方腔是典型的考察拐角流动模拟能力的算例.对于传统基于Boussinesq 假设的湍流模型,无法分辨雷诺正应力各向异性,导致很难得到正确的拐角涡结构.该算例的条件设置和网格同样参考NASA湍流模型资源网站.入口马赫数Mainlet3.9,基于方腔宽度的雷诺数ReD5.08×105,方腔的宽或高D25.4 mm.由于方腔为对称结构,可取四分之一部分进行计算.网格和边界条件设置如图6(a)所示.

图6 方腔流动示意图 (a)计算网格及边界条件;(b)SA 模型在 x/D=50 面上的速度矢量图;(c)SSG/LRR-ω 模型在x/D=50 面上的速度矢量图;(d)SSG/LRR-νt 模型在 x/D=50 面上的速度矢量图Fig.6.Sketch of square duct flow: (a)Mesh and boundary conditions;(b)velocity vector distributions on x/D=50 by SA;(c)velocity vector distributions on x/D=50 by SSG/LRR-ω;(d)velocity vector distributions on x/D=50 by SSG/LRR-νt .

图6(b)—(d)给出了三种湍流模型在x/D50剖面上的速度矢量图,颜色通过横向速度大小标识.在拐角处,会形成一对反向旋转且对称的角涡.对于SA 这类的线性涡黏模型而言,是无法有效捕捉到的.而对于雷诺应力模型,则天然具备优势.本文作者发展的SSG/LRR-νt模型同样继承了该特点.

图7 给出了x/D40和50 两个剖面上,流向速度沿对角线的分布.网格采用的 2 41×41×41 的粗网格.SSG/LRR-νt模型由于可以进行高精度离散,因而获得了较好的结果.而SSG/LRR-ω模型则无法采用七阶WCNS 格式稳定求解.

图7 流向速度沿对角线的分布 (a)x/D=40;(b)x/D=50Fig.7.Distribution of u-velocity along diagonal: (a)x/D=40;(b)x/D=50 .

3.4 NACA0012 翼型45°迎角分离流动

最后考察非定常分离湍流问题.算例选择45°迎角的NACA0012 翼型,来流马赫数Ma∞0.5,基于弦长的雷诺数Rec1.3×106.该算例常被用来考核像分离涡模拟(DES)等混合RANS/LES模型[29],而对于单纯的RANS 模型而言具有挑战性.图8(a)给出了计算网格,参考Strelets[29]工作,网格节点数为 1 93×103×31,展向长度为1c.计算采用非定常双时间步法,每个外迭代步推进0.01T(Tc/U∞).计算 5 0T后开始平均统计,一直到 1 50T结束.图8(b)对比了SSG/LRR-νt和SA模型得到的平均压力系数分布.可以看到: SA 模型明显高估了背风面(上表面)的吸力;而SSG/LRR-νt模型给出了与实验符合的结果.

图8 NACA0012 翼型网格及表面压力分布对比Fig.8.NACA0012 airfoil mesh and comparison of surface pressure distributions.

图9 分别给出了两种模型在100T时刻得到的展向涡量云图.传统的RANS 模型会在分离区高估涡黏性,导致小尺度涡结构被耗散掉.而本文提出的SSG/LRR-νt模型,νt尺度方程参考了k-kL模型具有尺度自适应能力((1)式生成项中包含湍流长度尺度和von Karman 长度尺度的比值).因而在WCNS 格式的配合下,分辨出了丰富的涡结构.