高阶矩湍流模型研究进展及挑战

王圣业 符 翔 杨小亮 郑皓榜 邓小刚∗

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

2 军事科学院,北京100190

1 引 言

湍流是流体在空间和时间上都具有急剧不规则和高度随机性脉动的一种流动状态,在自然界和工程技术各领域中极为常见(郑晓静和王国华2020).湍流问题作为流体力学中有名的难题,经过包括许多伟大科学家在内的学者长达100 多年的努力,已形成独特的研究体系(是勋刚1992).然而到今天,距离彻底弄清和完全精准地预测湍流仍有很长的路.流体力学中,对湍流问题的研究和预测包含三大手段: 理论、实验和数值计算.随着20 世纪70 年代计算机技术的飞速发展,数值计算(计算流体力学,CFD) 在研究湍流运动规律,解决实际湍流问题中发挥着越来越重要的作用(Slotnick et al.2014).

湍流的数值模拟方法一般包括三大类: 雷诺平均模拟(RANS) 方法、直接数值模拟(DNS) 方法和大涡模拟(LES) 方法.DNS 方法不添加模型而试图直接求解所有尺度的湍流.其网格分辨率及最大允许时间步长必须足够小,从而捕捉最小湍流尺度(Kolmogorov 尺度).而RANS 方法则仅关心最大尺度湍流,通过建立模型去模化所有湍流谱的特征.然而由于湍流的复杂性,很难寻找到一种尺度去近似模拟所有的流动情况(Wilcox 2006).LES 试图仅对最小尺度湍流进行模化,来消除DNS 求解宏观湍流问题时不必要的花费,因此相比RANS 方法可以建立更加简洁和普适的模型.理论上,为了精准可靠地预测湍流,最好是直接将所有湍流波动求解直到最小尺度(即DNS).然而,由于湍流最小和最大尺度的比与湍流雷诺数有关(λ/L=),使得所研究湍流问题的雷诺数会受到计算机硬件水平的限制.因此DNS 目前仅用于低雷诺数的湍流机理研究,在湍流结构生成和演化等方面取得了巨大的成就(傅德薰等2010).

在工程技术领域,RANS 方法扮演了中流砥柱的作用.RANS 方法的核心是对雷诺应力进行封闭建模,其思想可追溯到一百多年前Boussinesq (1877) 提出著名的涡黏性假设.其后,诸多流体力学先驱发展了一系列半经验理论(陈懋章2002).但由于只考虑了一阶湍流统计量,这些涡黏模型也被称为一阶矩模型.过去的三十多年间,涡黏模型在预测稳态流动问题方面获得了很大的成功.典型代表为在航空工程领域得到广泛应用的Sparlart-Allmarsa(SA)一方程模型(Sparlart 1994,2007) 和剪切应力输运(SST) 两方程模型(Menter 1994,2003).然而当流场中出现明显的流动分离时,传统涡黏模型往往会表现出不足,进而无法准确预测分离流动下飞行器的气动特性.比如,在分离涡处会高估涡黏性,无法分辨不同尺度的涡结构(周铸等2017); 在部件结合处无法估计湍流正应力各向异性,形成虚假分离(Rumsey et al.2016) 等.物理上更为完善的封闭建模方法是直接对湍流应力项建立微分输运方程.该方法最早可追溯到1940 年由我国著名力学大师周培源所建立的雷诺应力输运微分方程组(周培源1940).由于雷诺应力属于二阶湍流统计量,而在建立输运方程组的过程中会出现三阶(更高阶) 统计量.如果在二阶范围内对更高阶统计量建模,则称为二阶矩模型;如果对三阶统计量继续建立微分输运方程,则为三阶矩模型;依此类推.目前,学术界普遍将二阶矩及其以上的湍流封闭模型称为高阶矩模型(Hanjali´c & Launder 2011).

在RANS 方法框架下,雷诺应力(输运)模型(RSM)主要考虑以二阶矩封闭为主.Rotta (1951)进一步发展了周培源的工作,提出了完整的RSM 模型.其后,Donaldson 和Rosenabum (1968)提出了模型不变性概念,Lumley (1978) 发展了封闭性近似等,但都处于初级研究阶段.到了20 世纪90年代,计算机水平的飞速发展和一阶矩模型的不足使得人们又重新开始重视二阶矩模型.Launder(1989) 对过去二阶矩模型的发展进行了总结,并认为在随后几年中将一些模型纳入商业软件是相当可行的.Speziale 等(1994)认为RSM 方法给出了用于复杂湍流流动计算的LES 方法和RANS 方法之间的桥梁.他进一步讨论了雷诺应力模型的可实现性问题,并提出了简化的设计准则.其后,诸多学者基于上述准则,对雷诺应力输运方程中的未知项(如湍流耗散项、压力应变关联项和扩散项) 的建模进行了研究.而在这其中,最为核心并且困难的是压力应变关联项的建模.一方面由于压力应变关联项的量级和生成项相当,其在很多实际流动问题中十分重要,并且也是区别于涡黏模型的关键(Wilcox 2006); 另一方面,由于压力应变关联项很难通过实验测量得到,如何构造一个合理的模型需要高超的创造力.近年来,包括美国航空航天局(NASA)、德国宇航院(DLR)等机构纷纷提出重点发展二阶矩雷诺应力模型的计划.例如NASA 兰利研究中心在2012 年指出了几个关键研究热点(Malik 2012): “改进RANS 结果计算机制,使二阶矩封闭方法可实际应用于刚性方程问题; 在雷诺应力方程方面,改进长度尺度方程、压力应变模型、壁面处理以及压缩性影响、动态激波相互作用效应.”而在2014 年又进一步指出(Slotnick et al.2014): “对RSM 模型研究工作的投资包括模型的不同变种在CFD 软件中的应用、验证和确认.” DLR 更是走在前列,发展出的JHh-v2 雷诺应力模型 (Jakirli´c & Hanjali´c (2002) 提出的 JHh 模型改进版) 基于先进的近壁耗散修正,对近壁湍流尤其是压力梯度诱导分离等效应具有很好的模拟效果;通过欧洲FLOMANIA 计划发展的SSG/LRR-ω模型(Eisfeld&Brodersen 2005)借鉴了Menter 的混合函数思想,具有良好的鲁棒性,在一些航空工程问题中取得了较好的模拟效果(Togiti et al.2014,Rumsey 2015,Lee-Rausch et al.2016,董义道等2016,Wang et al.2018).我国学者也在雷诺应力模型发展方面提供了诸多开创性工作,例如符松教授提出的Fu-Launder-Tselepidakis 模型等(Fu 1987,1998,符松1994); 聂胜阳等将 SSG/LRR-ω模型与γ-Reθ转捩模型结合 (Nie 2017),扩展了转捩模型范围; 作者等提出了基于新型涡黏尺度的雷诺应力模型(EV-RSM),在鲁棒性方面获得了明显提升,以此成功实现了高阶矩模型的高精度离散及稳定求解(Wang et al.2020).

LES 被普遍用于非稳态流动的精细模拟,并且随着计算机水平的发展而逐渐进入工程领域.诸多学者及机构提出LES 方法作为下一代CFD 软件核心的预测或计划(Slotnick et al.2014).LES方法的核心是对亚滤波(或亚格子) 尺度应力进行封闭建模,其思想首先由气象学家Smagorinsky(1963) 提出.由于仅对最小尺度湍流进行模化,LES 框架下的湍流模型相比RANS 框架下更加普适.再加上LES 巨大的计算资源要求,使得模型被构造得更为简洁.目前,在实际流动问题中广泛被应用的亚滤波模型均是一阶矩模型,通过Boussinesq 的涡黏性假设与滤波NS 方程建立联系.然而,在近年来的一些精细流动研究中,很多学者指出了传统涡黏性模型存在的不足(Georgiadis et al.2010).与 RANS 方法类似,LES 方法中也存在高阶矩湍流模型.早在 1972 年便由气象学家Deardoff (1973,1974) 提出并在气象流动中进行了尝试,但他的工作大大地走在了时代的前面.近年来计算机硬件水平及湍流结构识别技术的进一步发展,为二阶矩亚滤波模型的研究提供可能(Stoellinger et al.2015).

湍流模拟方法正经历从RANS 到LES 的关键性转变,然而该过程也并非能一蹴而就(Slotnick et al.2014).工程问题中的特征雷诺数普遍较高,即使附着边界层内的最大尺度也会变得很小,采用LES 模型对网格尺度的要求也并不比DNS 减弱太多.如果在边界层附近采用各向异性的模型,如RANS 模型,而在远离壁面的区域采用LES 模型,则能大大降低对计算机硬件的要求.这便是混合RANS/LES 方法.Sparlart (1997) 基于SA 湍流模型提出了分离涡模拟方法(DES),是在实际工程中准确预测非定常湍流的里程碑工作.这之后,各种形式的混合RANS/LES 方法如雨后春笋般迅速发展(Frohlich & von Terzi 2008),同时在一定程度上被应用于真实工程分离流动模拟当中,取得了令人鼓舞的成果.然而,受限于线性涡黏假设,基于传统RANS 方法的混合形式出现了诸多不足(Greschner et al.2008).正如NASA 报告(Slotnicket al.2014) 指出的: 首先,模拟任何从边界层内部发展出来的分离流动,都需要RANS 方法的不断改进; 其次,为提高方法的鲁棒性,在边界层内,RANS 和 LES 需要自动进行光滑过渡.因此,改进混合 RANS/LES 方法,尤其是 “灰区” 问题,关键仍是RANS 模型.而高阶矩模型在对过渡区弱分离流的捕捉方面具有比涡黏模型更大的优势.如Chaouat (2011) 将亚滤波尺度应力模型与部分积分输运模型(PITM) 结合; Probst (2011)将DES 思路引入εh-RSM 模型,在翼型尾迹小分离模拟中展示了相对涡黏性模型DES 的优势; 作者等(Wang et al.2017) 将IDDES 引入SSG/LRR-ω模型,在翼型和三角翼失速模拟中取得了较好的效果.

图1

综上,高阶矩模型是湍流模式理论发展中举足轻重的部分.尤其是进入新世纪后,随着计算机硬件水平和数值求解技术的巨大提升,高阶矩模型又迎来了发展的小高潮.多个国外科研机构提出过复兴高阶矩模型的计划,尤其是美国NASA 在CFD 的2030 愿景中(Slotnick et al.2014) 将高阶矩模型作为RANS 部分研究和资助的主要方向,并将2020 年设置为关键里程碑节点(如图1).基于此,本文对高阶矩模型的发展情况进行了总结,包括高阶矩模型各项的建模过程、尺度提供方程的演化历程、数值求解技术的研究需求、其相对于传统涡黏模型的优势以及部分CFD 软件的集成情况等.抛砖引玉,以求在促进我国的湍流模式理论研究中做出一点贡献.

2 湍流量输运方程

2.1 二阶矩模型一般框架

如引言所述,周培源先生首次建立了一般湍流的雷诺应力所满足的微分输运方程组.但在该方程组中又出现了未知的三阶统计量(周培源1940).为了使方程封闭,周培源基于自由湍流条件引入了一些假设,得到了二阶矩模型的雏形(Chou 1945).其后,Rotta (1951)发展了周培源的思路,提出了更一般的二阶矩模型.至今,RANS 方法下的二阶矩模型基本都是基于他们的工作发展而来.本文设雷诺应力项(“ˆ”代表Favre 平均),六分量方程具有如下通用形式(Cecora et al.2015)

其中右端分别为生成项、压力与应变关联项(再分配项)、耗散项、扩散项及质量项.其中,生成项可精确得到

压力与应变关联项、耗散项和扩散项在二阶封闭范围内,均是雷诺应力、湍动能、耗散率、平均密度、平均速度等的函数.而质量项通常被忽略.

与RANS 方法中的二阶矩模型类似,LES 方法中的二阶矩模型(或称为亚格子应力模型)早在1973 年代就被Deardorff 提出,并被用于气象流动.本文设亚格子应力为则一般框架表示为(Chaouat et al.2017)

等式右端仍然包含生成项、压力与应变关联项(再分配项)、耗散项、扩散项等.其中生成项也是通过亚格子应力精确获得.而其余右端项需要建模处理.

由于 LES 方法下的二阶矩模型与 RANS 方法下有着相同的结构,如果引入过渡函数F使RANS 方法作用于边界层附近,而在远离边界层区采用LES 方法,则可构造混合RANS/LES 框架下的二阶矩模型

作者曾基于Sanchez-Rocha 和Menon (2009) 在可压缩流动中得到的混合RANS/LES 框架,通过选取SSG/LRR-ω模型作为RANS 部分,Deardorff 模型作为LES 部分,Menter 的F2函数作为过渡函数对塔形凸起物绕流进行了模拟.图2 表明混合RANS/LES 框架下的二阶矩模型具备准确分辨涡流动的潜力.

2.2 压力与应变关联项建模

上述二阶矩模型中除了生成项可以在二阶湍流量范围内精确表示,其他右端项均需建立模型处理.而其中最关键的一项是压力与应变关联项的建模.首先,该项与生成项阶数相同,在大多数复杂湍流流动中扮演着关键角色,也是雷诺应力模型区别于涡黏模型的核心.其次,由于该项涉及压力脉动,而试验中很难对其观测,因此建立一个合理的封闭需要很大程度的创造性(Wilcox 2006).

图2

Rotta (1951) 通过对压力脉动泊松方程的推导,将压力与应变关联项分解成慢项和快项

慢项主要决定湍流是否会趋于各向同性,通过定义一个无量纲来度量雷诺应力各向异性

并认为慢项直接与该度量成正比,即Aij=−C1εaij.Lumley 和Khajeh-Nouri (1974) 进一步提出了非线性的慢项表达式,可视为在线性项后添加了aij的张量积项

对于快项,在 1970 年代前学者普遍认为其远小于慢项,可以忽略 (Donaldson & Rosenbaum 1968,Daly & Harlow 1970).但其中也有不同的声音,如Crow (1968) 和Reynolds (1970) 提供的湍流例子表明快速压力应变的影响远远大于慢压力应变.Launder,Reece 和Rodi (1975) 完全基于运动学的考虑,设计了一种巧妙的封闭逼近,即著名的LRR 模型(慢项仍基于Rotta 的线性模型)

其中

该封闭最显著的特征之一是只存在一个封闭系数,即0.40.6.同时,LRR 还提出了一种简化版本

为便于区分,式 (8) 通常称为 LRR-QI 模型,式 (10) 称为 LRR-IP 模型.在计算壁面附近流动时,LRR 模型无法考虑压力– 反射效应,需要进行修正.包括 Gibson 和 Launder (1978) 以及Craft 和Launder (1992) 的修正形式.

另一方面,Lumley (1978) 从可实现性的角度出发,提出了Lumley 压力与应变关联模型.其后,Speziale,Sarkar 和Gatski (1991) 提出了类似但更简单的非线性模型,即SSG 模型

从形式上看,SSG 模型包括LRR 模型中的相关项,即,C2,为零.尤其是,SSG 模型可以不需要对压力– 反射效应进行修正,因此在壁面附近比LRR 表现更好.同一时期,还有Shih-Mansour-Chen 模型(Shih et al.1987),Fu-Launder-Tselepidakis 模型(Fu et al.1987)等,重点也均是对压力与应变关联项的非线性封闭.

Wilcox (2006) 通过研究指出,LRR 模型的压力– 反射问题与耗散尺度ε有关.通过将ε尺度替换为ω尺度,可避免在壁面处添加有关修正,即 Stress-ω模型.由于基于ω的 LRR 模型与基于ε的 SSG 模型均无需压力 – 反射修正,因此数学形式十分相似.Eisfeld 等 (2005) 将 Menter 在BSL(或 SST) 模型中的过渡函数引入到了压力与应变关联项的建模中,得到了 SSG/LRR 混合模型.该模型与式(11) 相同,只是系数Ci和通过过渡函数F1获得

图3 展示了SSG/LRR-ω模型的构造思路,其兼具SSG 模型和LRR 模型的优势,在航空复杂流动等实际应用中获得了良好的评价(Togiti et al.2014,Rumsey 2015,Lee-Rausch et al.2016,董义道等2016,Wang et al.2018).

图3

LES 方法中的二阶矩模型,关键也是对于压力与应变关联项的封闭建模.Deardoff(1973,1974)提出的亚格子应力模型与式(5) 类似,该项也被分为慢项和快项两部分

第一项系数cm与雷诺应力模型中Rotta 系数C1作用相同.第二项是针对均匀应变湍流在初始各向同性状态时获得的.相对于RANS 框架下的压力与应变关联项,LES 框架下的封闭要简陋很多.这也是由于计算机硬件以及试验观测的限制造成的.当前,随着DNS 的发展给压力与应变关联项的研究和标定提供了新的途径.如Jakirli´c 和Hanjali´c (2002)提出了改进近壁湍动能及应力耗散模拟的JHh-v2 雷诺应力模型,并利用完全发展的槽道流动(Rem= 5600) 进行了标定.Stoellinger 等(2015) 基于Hanjali´c等的工作得到了一种亚格子应力模型,并在低雷诺数槽道流动中与DNS 结果进行了对比分析.

2.3 扩散项建模

式(1) 中扩散项包含两部分: 湍流扩散和分子扩散.分子扩散与NS 方程一致,通常利用牛顿黏性应力公式得到(陈懋章2002).而湍流扩散含有压力脉动相关量和三阶速度脉动关联量,需要进行建模.在湍流扩散建模中最常用的方法是假设一个梯度输运过程.Daly 和Harlow (1970) 认为,雷诺应力的湍流扩散速度应与雷诺应力的梯度成正比,于是有

其中k2/ε为量纲决定系数,由 Mellor 和 Herring (1973) 采用.这种形式在数学上很简单,同时数值稳定也很好,尤其在复杂外形应用中获得了成功.但是,该形式不符合旋转不变性.Donaldson(1972) 提出了对称形式,其后Launder 等(1975) 对此进行了完善,提出了更一般的形式

2.4 耗散项建模

因为耗散发生在最小的湍流尺度,大多数二阶矩模型使用了Kolmogorov(1941)的局部各向同性假设,即

其中ε被称为各向同性耗散率,与标准两方程模型(k-ε或k-ω模型)类似,通过再建立尺度方程得到(参见本文第3 节).基于该形式,可自然利用k=−τij/2 关系得到标准的湍动能输运模型.

湍流耗散的准确建模对于整个雷诺应力模型至关重要.在至今很长的一段时期,关于雷诺应力模型的改进主要集中在湍流耗散上.一方面由于耗散在现实中是各向异性的,特别是在壁面附近.Hanjali´c 和 Launder (1976) 提出了低雷诺数耗散形式

其中bij为无量纲雷诺应力各向异性张量,即

同时,fs是一个低雷诺数阻尼函数 (或称为过渡函数),根据经验选择它随湍流雷诺数变化fs=(1+ReT/10)−1,ReT≡/(ε¯ν).其后,Kebede 等(1985)、Gilbert 和Kleiser (1991)、Launder 和Tselepidakis (1993) 等均对该形式进行了标定和发展.尤其是,Jakirli´c 和 Hanjali´c (2002) 在其 1993年版本(Hanjali´c & Jakirli´c 1993)的基础上提出了近壁耗散修正雷诺应力模型(JHh-v2 模型),将各向异性耗散率分解成均匀部分和非均匀部分

其中前者为均匀部分,仍由式(17) 得到; 后者是对应应力分量的黏性扩散,根据应力输运方程的解计算得出.该模型在翼型小分离流动中取得了良好的计算效果(Probst et al.2011).

在LES 框架下,耗散项与滤波/网格尺度相关.在Deardorf f提出的亚格子应力模型中使用了各向同性耗散模型

该式也被用于k方程亚格子模型(Sanchez-Rocha & Menon 2009).Stoellinger 等(2015) 基于Hanjali´c等工作得到的亚格子应力模型,同样采用了各向异性的耗散模型

该模型在壁面处耗散不为0,在低雷诺数槽道流动中获得了良好的效果,如图4.

引言中提到,受限于RANS 方法的局限,高阶矩RANS 模型在模拟非定常流动中也需要针对耗散项进行修正.得益于雷诺应力模型和亚格子应力模型在各向同性耗散项上的相似,发展高阶矩模型下的混合RANS/LES 方法可借鉴涡黏模型的做法.如Probst(2011)将DES 思路引入εh-RSM模型,Wang 等 (2017) 将 IDDES 方法引入 SSG/LRR-ω模型

图4

图5

2.5 更高阶矩模型一般框架

在对二阶矩模型开展研究和应用的同时,人们也一直在探索通过更高阶的展开,建立更为准确的湍流方程.然而受限于计算花费和数值刚性等问题,二阶矩模型尚未能在工程中被广泛采用,而更高阶矩模型目前更仅限于学术研究.图5 是Jeyapaul 等(2014)统计的更高阶矩模型所需的动量输运方程数.可以看到,随着维度的增加,微分方程数的增长是巨大的.

本文中对更高阶矩模型的讨论仅限于三阶和四阶,由于对更高阶项的建模还很不成熟,下面仅给出一般框架.对照式(1),三阶矩模型的一般框架为

右端项依次为生成项、压力应变关联项、耗散项和扩散项,其中后三项需要在不高于三阶湍流量的范围内进行建模处理.

可以证明,该假设对高斯速度分布是精确的.而参考二阶矩模型,Kawamura 等(1995) 提出湍流扩散速度应与三阶矩应力的梯度成正比,并依据DNS 数据给出了比例系数.除此之外,Grossman 和Narayan (1993),Cheng 等(2005),Gryanik 等(2005) 均提出了相关模型.

其次是耗散项.高阶矩湍流模型中的大多数耗散模型均基于各向同性假设,同时对其进行近壁修正.高阶耗散模型的一个特点是大于二阶矩的耗散壁值为零.这形成了一个简单的耗散边界条件.然而,近壁渐近线为0 需要通过阻尼或其他方式建模(Jeyapaul 等2014).一个比较准确的耗散模型由Dekeyser 和Launder (1985) 提出,其采用了广义梯度扩散假设.然而由于三阶耗散并非恒为负,该假设在简单槽道流动中存在问题.后针对该问题发展了改进版.Kawamura 等(1995) 考虑了近壁情况,提出了更为复杂的耗散模型.

最后是压力应变关联项.在三阶矩模型下,它仍然是最为困难的部分.Dekeyser 和 Launder(1985) 曾参照二阶矩的思路提出过相应的模型,其同样包括快项和慢项.其后,Kurbatskii 和Poroseva (1997),Kawamura 等(1995) 等基于该思路进行了发展.

四阶矩模型一般框架为

其中右端的压力应变关联项、耗散项和扩散项需要在不高于四阶湍流量的范围内进行建模处理.尤其是在四阶矩输运方程中会出现五阶相关项,可以假设其是二阶动量矩和三阶动量矩的函数(Jovanovi´c et al.1993).该假设认为概率分布函数为Gram-Charlier 级数,即截断近似关于高斯分布.表1 汇总了RANS 框架下二至四阶封闭模型的各项研究情况.随着阶数的提高,封闭建模的难度往往也更大.

3 尺度方程

经过上节涉及的建模过程,RANS 框架下的雷诺应力方程中往往仍包含一个未知量,即各向同性耗散率ε.在物理上实质是需要一个标量的湍流长度尺度或时间尺度来封闭整个模型.而在LES 框架下,这一尺度是通过滤波或网格提供的.因此,本小节的讨论均是在 RANS 框架下开展的.考虑到耗散率的计算方式以及相应的附加输运方程所求解变量存在不同,将这一类用来封闭模型系统的方程统称为湍流尺度方程.

在建立耗散率输运方程时需要考虑所模拟流动的湍流雷诺数,即Ret=k2/νε.由于湍流耗散发生在最小尺度上,对于高雷诺数流动通常采用Kolmogorov 的当地各向同性假设(Kolmogorov 1941) 来模拟耗散率张量,即εij=2/3εδij.Hanjali´c 和 Launder (1976) 基于ε的精确输运方程 (Daly1970) 建立了标准耗散率模型方程,即

表1 高阶矩RANS 模型各建模项研究汇总

其中方程右侧依次为扩散项、生成项和破坏项.该式与k-ε模型中的耗散率方程相比存在两点不同,即: (1) 生成项中的P直接由雷诺应力计算得到; (2) 扩散项包含一个各向异性扩散系数(Pope 2000).该模型的不足之处在于其生成项建模方式,即当湍流处于不平衡或者平均流动的时间尺度量级小于等于湍流的特征时间尺度时,小尺度湍流可能没有足够时间调整以适应大尺度上的变化(Al-Sharif 2011).标准耗散率模型方程在 LRR 模型 (Launder et al.1975) 和 SSG 模型(Speziale et al.1991) 中得到了应用.

与对k-ε模型的改进工作相似地,研究者们也对雷诺应力模型中的耗散率模型方程提出改进建议,如Launder (1996) 建议对方程系数进行修改以提高总体模拟效果,Hanjali´c 等(1997) 建议在方程中增加新的源项以提高对再附分离流的预测能力,等等.这其中针对低雷诺数效应所开展的改进工作具有代表性.比如,流动在接近壁面时,湍流雷诺数逐渐降低至0,耗散率张量的各向异性也越来越明显,因此有必要对湍流模型进行相应的修正以适应流动特征的变化.考虑低雷诺数效应时,由耗散率的精确输运方程(Daly 1970) 出发对模型进行修正,其具体形式如下

图6

式中右侧的前三项分别对应于耗散率的生成、破坏和扩散,模拟方式与高雷诺数模型相同; 在高雷诺数模型中被忽略的黏性扩散保留其精确形式;为了加入第四项的影响,可以将参数Cε1和Cε2相应地设置为Ret的函数; 第五项的一种常见建模方式为

作为代表,在 Hajali´c-Jakirli´c低雷诺数模型 (Hanjali´c 和 Launder 1976) 和 Shima 低雷诺数模型 (Shima 1998) 中采用了针对低雷诺数效应改进的ε输运方程来封闭整个模型系统.

在对k-ε模型的研究(Speziale et al.1992) 中发现,采用湍动能耗散率作为特征尺度存在两个问题:一是对于耗散率缺乏自然边界条件;二是在壁面附近需要处理高阶关联.这两个问题对于模型的稳定求解和壁面附近的模拟精度产生了不利的影响.针对这两个问题,研究者们采用了不同的尺度来构造尺度方程,比如均匀耗散率Chien 1982,Kalitzin et al.1996)、比耗散率ω(Menter 1994,Wilcox 1988)、湍流耗散时间τ(Kalitzin et al.1996,Kok & Spekreijse 2000) 等.当然,并不是所有在模型中出现的尺度都会被应用到雷诺应力模型中.图6 给出了应用于雷诺应力模型的尺度变量发展脉络.

为了解决湍动能耗散率缺乏自然边界条件的问题,部分研究者(Craft & Launder 1996,Craft1998) 采用了均匀耗散率来作为湍流特征尺度,其形式为

在两方程k −ω模型大获成功的基础上,比耗散率ω也被引入到雷诺应力模型的封闭处理中,其形式为

典型代表是 Wilcox 的 Stress-ω模型 (Wilcox 2006) 和 SSG/LRR-ω模型 (Eisfeld & Brodersen 2005).后者参考了SST 模型的设计思路(Menter,1994): 在近壁区采用ω方程,在远离壁面的区域采用ε方程,两者通过混合函数连接(参考图3),有效地解决了方程对来流湍流度敏感的问题以及方程在壁面附近数值刚性较大的问题.

尽管采用ω方程的SSTk-ω模型和 SSG/LRR-ω模型已经取得了较为广泛的应用,但ω方程本身在理论上仍然存在一些问题.Togiti 和Eisfield (2015) 就指出: 因为ω在壁面附近存在奇异行为且缺乏自然边界条件,基于ω方程的雷诺应力模型的精度在某种程度上取决于近壁网格分辨率.为了解决这两个问题,Togiti 和Eisfield 将两方程模型框架中的尺度g(Speziale et al.1992) 引入SSG/LRR 模型,其具体形式为

数值模拟结果表明,采用g方程的雷诺应力模型可以得到与ω方程相当的结果,并且其对于近壁网格分辨率的敏感性降低.

本文作者(Wang et al.2020)基于Menter 提出的k-kL模型思路,发展了涡黏尺度νt并用于构建新型雷诺应力模型.其具体形式为

该尺度不仅有自然的壁面边界条件(壁面处为0),同时具备利用Boussinesq 假设和Bradshaw 近似独立提供湍流应力的能力.以此能够实现传统涡黏模型和雷诺应力模型的结合,在计算稳定性和效率方面具备一定的优势.

利用雷诺应力模型模拟转捩问题,其实质也可归为湍流尺度问题.目前针对涡黏模型的转捩模型以Langtry 和Menter (2009) 提出的γ-Reθ模型为代表,通过在已有的湍流模型中附加间歇因子输运方程,达到控制湍流生成和耗散的目的.Nie 等(2017)成功将γ-Reθ模型扩展至SSG/LRR-ω模型,给雷诺应力模型求解转捩问题提供了新思路.

4 高阶矩模型数值求解

高阶矩湍流模型相比线性涡黏模型在数值求解时也更为困难.当采用线性涡流黏性模型时,线性雷诺应力张量与(平均流) 黏性应力张量成比例,有助于NS 方程解算器的稳定; 相反地,当使用高阶矩模型时,非物理的雷诺应力张量很容易引起失稳效应.因此,数值求解高阶矩模型方程时要保证获得的雷诺应力张量合理(即可实现性),这对相关物理限制、数值离散方法、边界条件和网格因素等提出了更高的挑战.图7 给出了影响高阶矩模型数值求解的主要因素,下面分别进行讨论.

图7

4.1 可实现性限制

Schumann (1977) 首先提出了可实现性这一概念.他要求雷诺应力模型得到的能量分量应为正,且非对角分量满足Schwarz 不等式.即

其后,Lumley (1978,1983) 指出在雷诺应力模型设计中应系统地使用可实现性限制.然而,随着实际应用逐渐增多,人们发现即使在模型的设计过程中满足可实现性,在数值求解时也可能产生违反物理的雷诺应力结果.

为解决该问题,Chassaing (2003) 提出了一个显式的强制可实现方法: 在每个计算网格点上检查可实现性约束,并且在违反的情况下,雷诺应力模型的局部解被强制设为零.而后,Chaouat(2006) 改进了这一措施,对正应力分量强制实施同性湍流,切应力分量仍强制为零.如下

该限制措施由于简单高效在复杂问题中获得了广泛的应用.如作者(Wang et al.2018) 基于该可实现性限制,采用三阶精度WCNS 格式对SSG/LRR-ω模型进行了离散求解,在跨声速涡轮叶片、翼身组合构型等真实外形求解中仍能保证良好的数值稳定性.

Mor-Yossef(2014)将最初针对两方程模型开发的无条件正收敛隐式格式(UPC)拓展至雷诺应力模型.UPC 方法具有任意时间步正应力保正性的特点,但其他可实现性约束仍有可能会违反.因此,针对式33(b) 定义了和满足的临近点平均值并施加了如下限制条件

4.2 数值离散方法

对雷诺应力模型的数值离散包括空间方面和时间方面.空间离散方面最受关注的是对流项的处理.目前,基于近似黎曼解的数值方法应用最为广泛,但在NS 方程中表现较好的低耗散通量格式可能会促进数值不稳定(Brun et al.2000,Gerolymos & Vallet 2007).Naser 等(2014) 指出,数值不稳定性应归因于对接触间断的错误处理.对于雷诺应力模型的对流通量由更具耗散性的质量通量来近似更为合适.而NS 方程仍可采用低耗散通量.另一方面,为提升黎曼求解的精度,插值过程需要更多的模板点.而在不光滑区需要在插值模板中增加限制器(包括加权等方式).由于雷诺应力各分量在无量纲后相比平均流动量具有量级小、变化大的特点,在插值过程中要格外注意经验参数的选择(Wang et al.2018).

时间求解方面,传统的显式格式虽然鲁棒,但在真实流动问题中很难突破时间步长的限制.Chassaing (2003) 提出了一种高效、鲁棒的隐式格式.该格式依赖于双时间步方法,用于定常流动时双时间步是人工时间增量(LDT).同时结合可实现性限制,实现了二阶精度雷诺应力模型的稳定求解.在隐式时间过程中,如何处理源项对收敛性和稳定性有较大的影响.正如Chaouat(2011)提出的建议: 将时间推进格式分裂,其中控制方程的对流项、扩散项在时间上显式处理,而雷诺应力方程的源项隐式处理.同时,他又针对性的设计了隐式格式来保持湍流耗散和法向雷诺应力的正性.Mor-Yossef (2014) 发展的UPC 格式同样是基于对源项的隐式处理,并在推广到雷诺应力模型时获得了很好的收敛性表现.

4.3 边界条件

边界条件是影响模型稳定性的重要因素,尤其是壁面边界条件.对于雷诺应力各分量,在壁面处应保证为0.但尺度方程由于变量形式不同,其边界条件处理方式和数值特性有明显差异.如第3 节所述,早期的RSM 通常基于耗散率尺度ε,但由于数值刚性相关问题,基于ε的RSM 很难用于复杂外形模拟.目前在真实工程模拟中,尤其航空领域,最常用的尺度提供变量是ω.特别是Eisfeld等(2005)引用了Menter 的混合思想提出的SSG/LRR-ω雷诺应力模型,在许多实际的航空流动中,如三维跨音速流动中的激波– 边界层相互作用、高升力下的自由涡和拐角流动中都得到了证明.然而,耗散率ω在壁面边界为无穷大值.Wilcox 及Menter 等采用的近似条件ωwall=60µ/(βy21)(非自然边界条件),在一定程度上会影响数值稳定性和灵敏度.这点在高精度离散时尤为明显.针对该问题,Bassi 等(2011) 提出了ω对数形式转换,有效避免了非物理状态; Hartmann 等(2011) 提出了新的ω边界条件,并在k方程中添加人工黏性.他们均实现了k-ω模型的高精度离散求解.而针对雷诺应力模型,作者(Wang et al.2020) 发展的νt尺度具备自然的边界条件(壁面处为0),基于经典的LU-SGS 方法和可实现性限制成功地实现了雷诺应力模型5-7 阶的高精度离散和稳定求解,其中源项均采用了Chaouat 建议的隐式处理.图8 展示了不同离散精度下翼型尾迹模拟结果对比.对RANS 方程和RSM 方程采用一致7 阶精度离散求解,在该粗糙网格上就能与实验吻合较好.

图8

4.4 网格因素

如同数值方法,雷诺应力模型对网格方面的要求也更为严苛.虽然目前尚未发现雷诺应力模型存在针对网格的特殊要求,但下面几方面是要格外注意的.首先,基于结构网格的有限差分格式是最早试图求解雷诺应力方程的方法(Wilcox 2006).但近年来,随着离散格式精度的不断提高,网格导数的几何守恒性对离散准确性,尤其是稳定性的影响也更加显著.针对NS 方程,Visbal 和Gaitonde (2002) 发现通过采用守恒形式的网格导数可大幅减小该误差.而多数CFD 软件中,雷诺应力模型往往采用与NS 方程相似的求解过程.因此减小雷诺应力模型在网格转换过程中的误差也是十分必要的.本文作者(Wang et al.2018) 将Deng 等(2011,2013) 提出的对称守恒网格导数算法(SCMM) 扩展至雷诺应力模型,并获得了良好的效果.

非结构网格相比结构网格更适用于工程问题,但由于高阶矩模型方程数较多,无论基于有限体积方法还是有限元方法都需要巨大的存储和计算花费.Mor-Yossef (2016) 成功地将其基于结构网格上(有限差分法) 提出UPC 格式扩展至非结构网格(有限体积法),但并未像基于有限差分法时开展三阶精度的尝试.主要原因是有限体积方法需要更复杂的模板点.采用DG 为代表的有限元方法可以克服模板向外扩展的困难,但又会引入新的问题,即壁面距离的准确计算.一种方式是通过求解程函方程(Eikonal-equation) 得到任意点到最近壁面边界的距离.Liu 等(2010) 提出了应用与内部流动问题的壁面距离算法.Schoenawa 和Hartmann (2014) 将该算法扩展至外流问题,并成功用于SST 模型的高精度求解.但在面对复杂外形时,还需添加人工黏性来保证稳定.而在非结构网格上,高于三阶的雷诺应力模型离散还少有研究.

图9

5 高阶矩模型对比线性涡黏模型

Boussinesq 假设是将湍流脉动运动类比分子热运动产生黏性的过程.对于剪切应变占主导的流动具有良好的适用性,这也是工程中最关心的边界层等为代表的流动类型.但随着CFD 技术求解问题的日渐复杂,基于 Boussinesq 假设的线性涡黏模型在众多流动问题中暴露出了不足.Wilcox(2006)总结了一些线性涡黏模型失败的并且值得注意的应用类型,包括: 平均应变率突然变化的流动、曲壁面上的流动、旋转流动(或涡流动)、三维流动等.

5.1 平均应变率突然变化的流动

考虑Tucher 和Reynolds (1968) 实施的实验.该实验中平均应变率维持在一个有限距离中,然后超过临界点后突然去除.湍流由均匀拉伸而变成各向异性,并在停止拉伸点(临界点) 的下游逐渐趋向各向同性.Wilcox 和Rubesin (1980) 采用他们提出的k-ω2模型对该流动进行了模拟,结果如图9 所示.当平均应变率突然在x ≈2.3 m 去除后,线性涡黏模型预测了一个各向同性的瞬时,即所有法向的雷诺应力变得平等.进而,湍流以有限速率接近各向同性.

更值得注意的是图9 中K值的变化,定义如下

如果x ≈2.3 m 的下游没有消除平均应变率,线性涡黏模型预测的K值应为连续的.但前面提到,由于Boussinesq 假设是对平均运动黏性的比拟,因此当平均应变率去除后,K值会产生突变.同时由于Boussinesq 假设对正应力的估计是各向同性的,因此K值在突变后将为0.后面,Wilcox(2006) 利用其发展的ω基雷诺应力模型对该问题又进行了模拟,所得结果很好地克服了线性涡黏模型的缺陷.

5.2 曲壁面上的流动

这里考虑一个经典的弯曲管道流动,由Smits 等(1979) 进行的实验.实验使用了一个高度为0.127 m 的恒截面积方形风管,中部具有30◦弯曲.CFD 模拟的主要焦点是评估凸壁曲率的湍流模型(图10a 中为下壁),被NASA 湍流模型网站(Rumsey 2020) 选用并收录了多个模型的计算结果.

本文引用该网站上SA,SA-RC 和 SSG/LRR-ω模型的结果进行对比分析.观察壁面摩阻,SA模型在凸壁面附近的预测明显高于实验值.真实的物理过程中,凸壁面附近的湍流水平往往低于直壁面.而SA 模型缺失了捕捉该现象的能力.同时对于凸壁面附近的雷诺应力,SA 模型也给出了较大的偏差.对于基于k的两方程模型,这种情况也依然存在.通过引入判断函数在曲面边界附近修正涡黏模型方程(称为旋转/曲线修正,如SA-RC (Shur et al.2000)) 是一种常用的工程方法,然而它们是不能用于任意流动的特殊修正方法(Wilcox 2006).冯·卡门曾提出关于曲壁面上流动的经典稳定性论点,认为k方程应该是的方程更合适.而一些 RC 修正也正是基于概论点构建的,如Wilcox 和Chambers (1977) 修正等.从这个角度讲,雷诺应力模型自然具备预测流线曲率影响的能力.图10 中SSG/LRR-ω模型的结果也印证了这点.对于凹壁面(图10a 中为上壁),由于存在Gortler 涡导致所有二维模拟都较差.但雷诺应力模型往往仍好于线性涡黏模型.

5.3 旋转/涡流动

考虑一个航空工业中最典型的翼梢涡流动,实验在 NASA 兰利研究中心的低速风洞中开展(Chow 1993,1997).实验模型是一个半展长机翼(翼型为NACA0012),翼根固定于风洞壁面,安装角 10◦(即攻角 10◦).如图11(a)所示,RANS 模拟的主要难点在于翼梢涡的保持,即过大的物理黏性会导致翼梢涡在向下游的过程中被逐渐耗散.图11(b)中给出了DLR 采用TAU 软件开展的CFD 模拟 (Cecora et al.2015),RANS 模型包括原始版 SA,SST,SSG/LRR-ω以及 JHh-v2 模型,后两者是雷诺应力模型.横坐标表示以翼梢后缘为原点的x方向坐标,纵坐标表示x对应截面上涡核中心流向速度.可以看到,涡黏模型由于过高估计涡核处的涡黏系数,导致翼梢涡迅速衰减.而两个雷诺应力模型无此现象.

在非定常涡流动中,基于高阶矩模型的混合RANS/LES 方法依然具备潜力.图12 给出了SSTURANS,SST-IDDES,SSG/LRR-URANS 和SSG/LRR-IDDES 方法对24.6◦迎角钝前缘三角翼的模拟情况(Wang et al.2017).实验在NASA Langley 中心的风洞实施(Chu & Luckring 1996),来流马赫数为0.85,基于三角翼气动弦长的雷诺数达6 百万.实验迎角增加至24.6◦时,主涡由于激波诱导作用而发生破裂,破裂位置位于翼根弦长60%∼80%间.SST-URANS 得到的展向压力分布,始终存在两个峰值(主涡和二次涡),说明其推迟了涡破裂的发生.SSG/LRR-URANS 在80%占位处峰值消失,说明捕捉到了涡破裂.但对破裂后非定常流动的模拟较差.而添加IDDES 修正可弥补这点,即SSG/LRR-IDDES 模型.而对比传统的SST-IDDES,SSG/LRR-IDDES 的结果更接近实验值.

图10

图11

5.4 三维流动

仍然考虑一个航空工业中的典型流动,翼身结合处流动(典型的三维T 型区,如图13a).该问题来源于第五届AIAA 阻力预测会议(DPW-5),计算模型为波音公司设计的通用(宽体客机)研究模型CRM (Vassberg et al.2008).该飞机模型的“设计巡航速度” 为Ma= 0.85,在美国NASA 的国家跨声速风洞以及埃姆斯研究中心风洞中均进行了详细的风洞实验.实验中在小迎角下(4◦附近),翼身结合处并未发现分离涡.而在DPW-5 会议提交的CFD 模拟中,涡黏模型普遍给出了分离涡,并且导致了“升力突降” 的现象(Levy et al.2013).随后,该问题被广泛关注,包括NASA 专门组织了拐角流动系列实验(Rumsey 2016)等等.图13 为德国宇航院Togiti 等(2014)对该问题开展的研究,采用SA 和SST 等涡黏模型在翼身连接处出现分离涡.主要原因是Boussinesq 假设无法考虑雷诺正应力的各向异性,即在三维流动中存在潜在的缺陷.而雷诺应力模型自然无此问题.参考雷诺应力模型,Spalart 等(2000)提出了对Boussinesq 假设的修正,即二次本构关系(QCR),很好地解决了该问题.但如同RC 修正,这是一种不通用的方式.

本节通过数个典型流动案例说明了高阶矩湍流模型相比于传统涡黏模型存在的几方面优势.然而,优势的获得往往伴随着一定的代价.最主要的是计算花费和鲁棒性要求的提升.对于前者,多位学者(Mor-Yossef 2014,C´ecora et al.2015,Wang et al.2020)针对不同的算例开展过评估.虽然在相同网格上,高阶矩模型需要高出涡黏模模型30%∼60%的计算花费.但考虑到相同网格时高阶矩模型获得更准确结果的潜力,尤其针对涡流动.因此在达到同等准度条件下,其对网格量的需求往往不会高于涡黏模型.而对于鲁棒性,则仍需要在模式理论、数值算法等多个方面继续开展研究.这也是目前影响高阶矩模型广泛应用最大的挑战.

图12

6 CFD 软件中的高阶矩湍流模型

CFD 软件作为湍流模型应用于实际流动的载体,对于一些模型的发展起到了巨大的推动作用.例如由Menter(1994)提出的SST 湍流模型,以及其衍生的γ-Reθ转捩模型(Menter et al.2006)、尺度自适应 (SAS) 模型(Menteret al.2003) 等,借助ANSYS 公司旗下的CFX,Fluent 等 CFD 软件在世界范围内得到了广泛的应用.而上述模型之所以被纳入该软件,除了Menter 本人作为ANSYS公司技术顾问外,模型本身也必须具备鲁棒性强、适应范围广等优秀的性能.在CFD 软件发展的早期,集成的湍流模型基本上为低阶矩模型.直到1990 年代后,随着计算机硬件水平的提升和二阶矩模式理论在一定程度上的成熟,其中一些表现较好的模型才开始纳入商业CFD 软件(Launder 1989).

图13

在ANSYS Fluent 中(以2019 R3 版本为例),二阶矩湍流模型被分为两大类: 基于ε的模型和基于ω的模型.前者默认采用Fu 等(1987)在LRR-IP 模型基础上改进的线性压力–应变模型(也称 Fu-Launder-Leschziner 模型),还包括二次压力 – 应变模型 (即 SSG 模型).后者则包括 Stress-ω模型和Stress-BSL 模型.这两种模型均基于LRR-QI 压力– 应变模型,区别在于尺度方程ω分别采用的Wilcox (2006) 版本和Menter (1994) 版本.在软件中对模型的可选项还包括边界条件、壁面反射效应和近壁处理等内容.同时针对特殊问题,可添加SAS 修正(分离涡问题)、低雷诺数数修正或剪切流修正等.

ANSYS CFX 与Fluent 相似,雷诺应力模型也被分为了两大类,即基于ε的模型和基于ω的模型.前者包括了LRR-IP 模型、LRR-QI 模型和SSG 模型,后者则包括了Stress-ω模型和Stress-BSL模型.值得注意的是,在该软件中对所有模型均采用了先进壁面函数.

NUMECA FINE/Marine 中也集成了雷诺应力模型,即Rij-ω模型.该模型是在 Shima (1998)的低雷诺数模型基础上将ε方程替换为ω方程(Deng & Visonneau,1997) 并采用了新的近壁模型而得到的.而压力与应变关联项是在LRR-QI 模型基础上的改进版.

除了商业CFD 软件,专用CFD 软件也为一些二阶矩模型在该领域的发展和推广起到了积极作用.例如由美国NASA 主导开发的CFL3D,FUN3D 软件、德国DLR 开发的TAU 软件均集成了SSG/LRR-ω模型.通过系统的验证、确认及应用,使该模型成为针对航空航天气动问题中高阶矩模型的标杆(Eisfeld et al.2016).表2 给出了部分CFD 软件中高阶矩模型集成情况.

表2 部分CFD 软件中高阶矩模型情况

7 高阶矩湍流模型挑战与展望

高阶矩模型是湍流模式理论发展中举足轻重的部分,对其物理机制和数值特性的研究不仅是认识湍流的基础科学问题,更是解决复杂湍流的工程需求.从上述第2∼6 节的进展情况看,高阶矩模型未来面临的最大挑战是如何确保在不低于涡黏模型准确度的同时,实现鲁棒、高效的求解,进而得到更广泛的应用.这需要在湍流建模理论、尺度提供方程、数值离散方法、CFD 软件集成及验证与确认等方面持续开展工作,如图14.

(1)应力输运方程中各源项(尤其压力与应变关联项)的物理机制及建模理论仍然将是高阶矩模型研究的核心.例如,1980 年代末开始兴起的重整化速度群(Renormalization group)理论.Yakhot&Orszag(1986)首先将其用来分析湍流并成功地建立了基于该理论的k-ε模型(RNGk-ε模型).该理论基于统计学方法,通过将自然空间中湍流方程转化为傅里叶空间(谱空间) 进行分析,得到独立于经验标定的湍流方程参数.其后,Rubinstein & Barton(1990) 将重整化速度群理论用于发展二阶矩模型,取得了令人鼓舞的效果,表明了该理论在推导高阶矩模型中具有很好的潜力.系统性逆尺度(Systematic upscaling) 是一种新的多尺度计算方法,用于精确推导在越来越大尺度上给定的物理系统的控制方程(Brandt 2005).该方法结合了过去出现的两种多层计算模式的优势,即应用数学中的多重网格和理论物理中的重整化速度群.美国NASA 在CFD 的2030 愿景报告(Slotnick et al.2014) 中将其列为有希望推动湍流模式理论变革的方法之一.机器学习(Machine learning) 是通过构建算法从数据中建立模型.近年来,得益于人工智能发展的高潮,机器学习也被成功用于湍流建模,并迅速成为本领域的热门(Lav et al.2019,张伟伟等2019).Tracey 等(2015) 通过神经网络模型构建了SA 方程中新的源项,而模型训练所采用的数据可通过实验结果或高分辨数据反演得到.相似的思路也可用于k方程模型(Zhang 2015) 甚至高阶矩模型.虽然压力与应变关联作用很难通过实验直接测定,但基于DNS 数据对压力与应变关联项进行机器学习建模仍有相当大的潜力.

图14

(2) 尺度提供方程一直是RANS 理论中的关键.从图6 的发展趋势来看,不同的尺度方程通过互补、融合可以实现更优秀的特性,如更符合物理的边界条件、更小的数值刚性制约以及更好的分离涡求解能力等等.Wang 等(2020) 推导的νt尺度方程,通过引入Menter 的k-kL模型近壁耗散修正,提升了雷诺应力模型的近壁渐进特性.Eisfeld 和Rumsey (2020) 将Yap (1987)在ε方程采用的长度尺度修正引入到了ω方程,使新的SSG/LRR-ω模型能更好地模拟后体分离问题.

(3) 数值离散方法是实现高阶矩模型广泛应用的关键.从软件集成的情况来看,二阶精度下的雷诺应力方程离散已比较成熟,无论针对结构网格采用有限差分离散还是非结构网格采用有限体积离散.但上升到高精度离散时,稳定性和收敛性则远不如低阶矩模型.而高阶矩模型的优势,又恰好是针对旋转、分离这类复杂的流动,需要高精度格式的低耗散特性作支撑.因此,未来的关键和难点是如何实现高阶矩模型的高阶精度离散及稳定求解.

(4) CFD 软件能够为湍流模型的应用与推广提供平台.但从表2 总结的情况看,国外开发的商业软件极少会采用我国原创的湍流模式.未来对于湍流模式的推广,一种是基于开源CFD 软件平台(如Openfoam 等),但主要适合于学术研究和交流; 一种是走自主发展的道路.目前,诸多科研院所自主开发了CFD 软件或代码,同时集成了各自发展的湍流模式.但由于软件规模与国外商业CFD 软件差距悬殊,影响力仍十分有限.只有从国家层面加以推动整合,才能形成有竞争力的软件产品,为我国原创的湍流模式提供良好的应用平台.

(5) 无论是开展新型高阶矩模型研究,还是将已有模型纳入CFD 软件作为产品模块,都缺少不了系统的验证与确认.例如,SSG/LRR-ω模型通过参与美国NASA 组织的湍流模型资源(TMR)网站(Rumsey 2020)、第五次阻力预测会议(DPW5) (Togiti et al.2014)、第二次高升力预测会议(HiLiftPW-2)(Thompson & Hassan 2015) 等的验证与确认,得到了广泛的好评和迅速的推广.我国在近些年也加大了这方面的投入,如空气动力研究与发展中心组织的第一届航空CFD 可信度研讨会 (AeCW-1) (王运涛等2019),已经有高阶矩模型参与对比(聂胜阳等2019).且标准模型和实验数据均为自主设计和测试.因此,未来为支撑原创湍流模式和自主CFD 软件的发展,验证与确认工作必须走在前面.

8 结论

本文从CFD 由RANS 到LES 的发展趋势谈起,对高阶矩湍流模型进行了综述.重点总结了高阶矩模型各项的建模过程、尺度提供方程的演化历程和数值求解技术的研究需求.同时,通过典型湍流问题展示了其相对于传统涡黏模型的优势,并且给出了部分CFD 软件对高阶矩模型的集成情况.最后对高阶矩湍流模型未来的研究与发展方向进行了展望.本文得到的主要结论如下:

(1) 进入新世纪后,随着计算机硬件水平和数值求解技术的发展,高阶矩模型又再次迎来了发展的高潮.多个国家的科研机构均提出过复兴高阶矩模型的计划.然而,在工业领域,高阶矩模型仍未被广泛推广,因此需要进一步加强应用性研究.

(2) 应力输运方程中各源项(尤其压力与应变关联项) 的建模理论是高阶矩模型研究的核心.近年来随着DNS 的不断突破,高分辨实验观测技术的发展给方程建模提供了更准确、丰富的信息.而诸如重整化速度群、系统性逆尺度、机器学习等理论的发展,有望为高阶矩建模提供新思路.

(3) 尺度项/方程对高阶矩模型的实际表现有着综合性的影响.根据历史上尺度方程的演化规律可以得出其研究重点包括更符合物理的边界条件、更小的数值刚性制约、更好的分离涡求解能力以及对转捩模拟的能力等.

(4)当前专门针对高阶矩模型的数值方法研究还很少.从学术的角度,高阶矩模型还很难做到像涡黏模型一样,与NS 方程保持一致的任意阶精度求解.这将导致在湍流量的计算上产生过多误差,甚至有时比涡黏模型的效果差.从工程的角度,高阶矩模型仍然缺乏鲁棒性.这是制约高阶矩模型广泛应用的关键问题.

(5) 高阶矩模型相比涡黏模型的优势仍然不够清晰,这需要持续的开展验证与确认研究.而随着湍流模拟方法经历由RANS 到LES 的转变,LES 框架下Boussinesq 假设是否更普适性以及建立输运方程模型是否有优势等问题也需要开展探索.

(6) 高阶矩模型的推广需要优秀的CFD 软件平台,而国外开发的商业软件极少会采用我国原创的湍流模式.对于我国而言,只有坚持走自主创新之路,形成有竞争力的软件产品,才能真正为原创的湍流模式提供良好的应用平台.

致 谢国家自然科学基金(12002379)、湖南省自然科学基金(2020JJ5648)、国防科技大学科研计划(ZK20-43) 和国家专项工程(GJXM92579) 资助项目.