王新光，陈 琦，万 钊，张爱婧，陈坚强

(1. 空气动力学国家重点实验室，绵阳 621010；2. 中国空气动力研究与发展中心计算所，绵阳 621000)

0 引 言

近年来，国家对于研发具有自主知识产权的自主品牌计算流体力学软件，给予了前所未有的重视。众所周知，湍流流动广泛存在于航空、航天、航海以及人们的日常生活中，湍流模拟方法是计算流体力学界长期以来面临的研究难题和研究热点。商业软件中通常使用RANS方法模拟湍流，常常需要在黏性底层中布置一定数目的网格点(通常y+≈1)，从而使得壁面附近的网格十分细密，不仅会极大地增加迭代收敛的计算步数，还会带来较为严重的数值刚性问题，从而导致迭代计算的稳定性下降。在复杂外形的数值计算中，网格量骤增，导致90%的计算时间都消耗在壁面附近1%计算域内。因此，目前主流商业软件在工程湍流模拟应用中默认采用壁面函数方法对湍流边界层进行模拟[1-2]，壁面函数方法可以大幅放宽壁面第一层网格的尺度(通常30

基于对数定律的标准壁面函数已经被大量不可压缩数值实验确认，其中包括DNS模拟零压力梯度边界层和槽道流等结果。随着近年DNS模拟能力的增强，涌现了大量可压缩湍流边界层DNS结果,其中文献[4]使用DNS模拟了可压缩低雷诺数槽道湍流，当马赫数较小时标准壁面函数和DNS结果较为吻合，但是随着马赫数的增加基于不可压流动的标准壁面函数和DNS之间的偏差逐渐增大。对于Ma2.48的槽道湍流，在黏性底层边缘附近y+=10，相较于壁面参数，密度减小了50%，温度增高了60%，表明低雷诺数超声速槽道湍流存在强可压缩性。

尽管壁面函数的研究工作已经持续开展了50多年，但其中大多数研究成果的应用工作并不乐观，如壁面函数与用于全局流动问题模拟的湍流模型不相容的问题[5-6]，容易引起数值模拟结果对壁面附近第一层网格节点位置的强烈依赖性，导致分离流动的数值模拟结果精度很差，因而提出了网格和流动自适应的壁面函数方法，保证在壁面函数有效的情况下对近壁流动物理特性进行恰当的分辨，特别是针对非平衡效应显著的流动区域、逆压梯度导致的流动分离和再附点附近区域，以获得比较准确的气动力/热特性。此外，文献[7]也针对壁面函数在RANS模拟中的应用问题开展了相关研究工作。

近年来可压缩壁面函数的代表性工作，文献[8]使用Crocco-Busemann方程，忽略压力梯度的影响，将速度和温度函数耦合起来，发展了适用于可压缩湍流边界层的壁面函数。国内壁面函数研究起步较晚，其中文献[9-10]将文献[8]中考虑流动可压缩性和热传导的壁面函数耦合到了k-ω两方程模型中；文献[11]和文献[12]分别通过数值实验和风洞实验对壁面函数的系数进行了修正研究；文献[13]对文献[8]中的壁面函数进行了修正研究。由于壁面律在存在逆压梯度的复杂流动中是否适用无法确定[14-15]，文献[8]忽略对流项和压力梯度影响的壁面函数对于分离点和再附点附近流动的模拟会存在困难。

除了标准壁面函数，由于近年来解析壁面函数[16]在壁面处不涉及太多假设，保留了边界层方程中的对流项和压力梯度项，因此在不可压缩流体中已经得到应用。由于解析壁面函数在粗网格上的预测精度可以接近低雷诺数模型的结果，而计算时间比低雷诺数低一到两个数量级。同时由于解析壁面函数在实际编程，计算中的鲁棒性和可操作性等优点，都使得其工程应用得到广泛关注。而将解析壁面函数用于可压缩流体中时，由于需要考虑流体可压缩性，其方法研究对于先进壁面函数在复杂超声速和高超声速流动中的应用也具有实际意义。本文基于OpenFoamV 5.0软件平台将目前应用于不可压缩流动的解析壁面函数，通过可压缩湍流边界层控制方程的简化，考虑密度变化、对流项变化和黏性耗散项对解析壁面函数的影响，探索其在可压缩流动中的应用，并通过二维超声速激波边界层干扰算例进行验证。

1 物理模型

本文解算器使用Favre平均NS方程，无黏通量由Roe-Pike格式计算，黏性通量采用二阶中心差分格式，时间推进采用一阶欧拉格式。根据完全气体假设，状态方程为p=ρRT，空气黏性系数使用Sutherland方程计算。壁面函数方法基于高雷诺数RANS模型，即标准k-ε模型。对于存在壁面的流动问题，局部雷诺数Rt=k2/ν/ε很小，黏性效应不能忽略，因此文献[17]采用基于Rt的阻尼函数来计算近壁面湍流尺度，并在ε方程中添加Yap修正[18]，来降低在分离流动中ε方程得到的近壁面湍流，即低雷诺Launder-Sharmak-ε模型。

2 壁面函数

2.1 标准壁面函数(Standard Wall Function)

在局部平衡湍流边界层内，壁面定律为：

(1)

其中:U+和y+定义为：

(2)

(3)

(4)

(5)

其中:常数cl=2.55，下标n表示壁面第一次网格面,如图1所示。近壁面网格点P的壁面剪切应力为：

(6)

2.2 解析壁面函数(Analytical Wall Function)

解析壁面函数[16]在壁面附近通过可压缩湍流边界层假设简化输运方程，考虑对流项和压力梯度的影响，壁面附近x-y平面上简化的动量和能量方程为：

(7)

(8)

其中:P表示压力;T表示温度;μ和μt分别表示层流黏性系数和湍流黏性系数;Pr和Prt分别表示层流普朗特数和湍流普朗特数。湍流黏性系数为：

μt=0 当y

(9)

(10)

(11)

其中:C表示动量简化式(7)中的压力梯度和对流项。

(12)

(13)

(14)

(15)

(16)

湍动能方程中的在壁面附近生成项为：

(17)

对于壁面附近的耗散项，在全湍流区为:

(18)

在黏性底层内耗散项为：

(19)

为了消除在黏性底层交界面处的间断(图1(a))，重新选取一个位置(使用下标d表示)，令耗散项在壁面网格内连续分布(图1(b))，需满足条件：

图1 耗散项在壁面网格内分布示意图Fig.1 Schech map of dissipation term in the wall cell

(20)

因此，在壁面网格内的平均耗散项为：

(21)

对于壁面温度的解析表达式可以通过类似于壁面速度的两次积分得到，对于绝热壁有：

(22)

其中，Cth下标1表示黏性底层的值，2表示全湍流区的值：

等温壁为：

(23)

原始解析壁面函数详情可见文献[16]。

2.3 可压缩修正(Modification of AWF)

可压缩流动的能量方程通常可表示为：

(24)

式中:E表示总能,σ表示流体应力张量,q表示热流。

为了增强壁面函数的鲁棒性，本文发展的解析壁面函数和主网格的能量方程统一使用式(24)，对于可压缩湍流边界黏性耗散项不能忽略[19]，因此在解析壁面函数能量方程中需包含粘性耗散项。

湍流壁面函数方法使用粗网格进行数值计算，第一层网格壁面距离通常满足20

(25)

(26)

式(26)中右端第二项表示黏性耗散项，使用式(25)积分得到的解析速度计算，方程中的对流项和压力梯度分别表示为：

上式D、C、Dth中下标1表示黏性底层的值，2表示全湍流区的值,在黏性底层内和全湍流区对简化动量式(25)和能量式(26)积分，可得到解析速度和温度表达式，具体形式如下：

(27)

(28)

其中：

系数N表示：

黏性底层温度解析表达式为：

(29)

其中：

全湍流区温度解析表达式为：

(30)

其中：

通过壁面第一层网格内的解析速度和温度分布可以得到壁面剪切应力和壁面热流。同时考虑到在使用粗网格求解时，为了更加准确描述主网格控制方程的黏性耗散项，即主方程中壁面网格内的平均粘性耗散项为：

在程序中通过解析速度表达式(27)、(28)数值求解。

3 数值结果

本文使用斜射激波-湍流边界层干扰来验证发展的可压缩修正解析壁面函数(MAWF)，三种来流马赫数实验来流条件列于表1中[20-22]。来流边界为充分发展的湍流边界层，其边界层厚度δ0和动量厚度θ0如表1所示，具有1.5%的自由来流湍流强度和黏性系数比率μt/μ=10。

表1 斜激波边界层干扰来流条件Table 1 Inflow conditions for the impinging shock interactions

计算中使用结构网格，对于低雷诺数湍流模型，壁面第一层网格较密，第一层网格的y+≈1,对于高雷诺数湍流模型，第一层网格y+大约为30。为保证计算结果与网格无关，在开展结果对比之前，分别对表1中算例开展了网格无关性研究，计算结果表明所有算例均取得网格无关性结果，详情可参考文献[23]，其中Ma=2.9算例标准壁面函数法粗网格结果如图2所示,其中xI表示无黏性情况下斜激波入射位置坐标。

图2 Ma=2.9斜激波边界层干扰的网格收敛性研究Fig.2 Grid independency study for the Ma=2.9 impinging shock interaction

3.1 Ma=2.9斜激波边界层干扰[20]

计算网格入口边界的下部使用充分发展的湍流边界层，入口边界的上部使用无粘斜激波关系式得到的斜激波后参数，在壁面处，使用等温无滑壁边界条件，壁温为Tw=271K。在顶部和出口边界处，使用Neumann边界条件。本算例密网格为200×90，粗网格为150×60。

图3比较了10°和13°斜激波的壁面压力、摩擦系数和壁面热流，图中湍流模拟方法分别为标准壁面函数(SWF)、原始解析壁面函数(AWF)、可压缩修正解析壁面函数(MAWF)和低雷诺数Launder-Sharmak-ε模型(LS)。对于β=10°，所有模型均返回相似的壁面压力，和实验数据[20]吻合。从摩擦系数结果来看，在无粘激波反射点附近存在一个小的分离区。LS模型的结果与实验结果非常接近，壁面函数不能复现分离区，其中SWF倾向于低估分离区上下游的摩擦系数。壁面热流不同模型的结果差异较大，其中MAWF结果更接近于密网格LS模型。对于β=13°，所有模型预测的壁面压力和摩擦系数结果与实验值较为吻合，MAWF预测的壁面热流和密网格LS模型较为接近，而SWF和AWF均不能准确预测壁面热流，在分离区前后甚至出现热流符号相反的现象。

图3 Ma=2.9算例不同入射角壁面压力(上)、摩擦系数(中)和壁面热流(下)Fig.3 Wall pressure(top)，skin-friction(middle) and wall heat-flux(bottom) distribution for the Ma=2.9 case with different impinging angles

3.2 Ma=5.0斜激波边界层干扰[21]

本算例网格类似于Ma=3.0算例，其入口边界使用充分发展的可压缩湍流边界层，其中动量厚度如表1所示，在壁面处使用等温无滑壁边界条件，壁温为Tw=300 K。本算例密网格为240×80，粗网格为120×45。在相同设置的情况下，不同湍流模拟方法消耗的计算机时间对比如图4所示，通过对比表明使用粗网格的壁面函数将大幅减少计算时间，仅为密网格低雷诺数模型计算时间的5%左右，一是由于网格量大幅减少，二是由于壁面距离的增加，在相同CFL数时计算时间步长的增加，收敛速度加快。

图4 Ma=5.0算例不同模型计算消耗时间对比Fig.4 The CPU time comparison using different turbulence models for the Ma=5.0 case

图5比较了不同入射角的壁面压力、摩擦系数和壁面热流，从图中可知四种数值模拟方法得到的分离区较试验结果[21]较小。当激波边界层干扰加强时，AWF得到的壁面摩擦系数和Stanton数在分离起始位置出现非物理振荡，主要原因是AWF使用壁面第一个网格点信息计算对流项，如式(11)，由于分离区附近流向速度变化较大，对流项在整个壁面网格内取常数，将导致AWF预测能力降低。而MAWF认为对流项在壁面处为零，以无量纲壁面距离的平方项递增，从物理上更符合边界层内流动规律，从而提高预测精度。对于三种不同角度斜激波，MAWF相较于AWF，壁面热流的预测能力大幅上升，主要得益于能量方程简化中保留了黏性耗散项，使得粗网格预测的壁面热流精度大幅改善，接近密网格LS的结果。

3.3 Ma=7.2斜激波边界层干扰[22]

本算例由于激波生成器拐角处的膨胀波会对下表面产生影响，因此在数值模拟的过程中加入了激波产生器的固壁边界，激波产生器夹角β=7.5°和15°，且由于本算例风洞试验模型是轴对称体，因此z方向按照试验外形旋转2°，其余边界设置和Ma=2.9和Ma=5.0算例相同，壁面温度为300K。本算例密网格为420×105，粗网格为280×40。

图6给出了四种不同湍流模型的数值结果与试验值[22]的对比，其中MAWF大幅度提升了壁面热流的预测精度，在分离区前后得到和密网格LS模型相同的结果，且避免了AWF在分离区起始位置附近的非物理振荡。和实验值相比，数值结果高估了分离区的大小和分离区内的壁面压力、摩擦系数和壁面热流，这主要是RANS模型的局限性导致，文献[24]总结了不同RANS模型预测的β=15°壁面热流结果，所有模型均出现高估分离区热流的情况。

图6 Ma=7.2算例不同入射角的壁面压力(上)、摩擦系数(中)和壁面热流(下)Fig.6 Wall pressure(top)， skin-friction(centre) and Stanton number (bottom) distribution for the Ma=7.2 case with different impinging angles

4 结 论

本文通过可压缩湍流边界层特点，考虑边界层内密度的变化，将基于不可压缩流动的解析壁面函数推广到可压缩流动原始解析壁面函数(AWF)，并考虑边界层内对流项变化和粘性耗散项的影响，构造了适用于可压缩修正的解析壁面函数(MAWF)，通过三个不同来流马赫数的斜激波边界层干扰算例进行测试，并与LS和SWF进行了比较。主要结论如下：

1)整体来看四种数值模拟方法预测的壁面压力差异不大，且与实验吻合良好。

2)整体来看四种数值模拟方法预测的壁面摩擦系数差异不大，标准壁面函数预测的数值偏小，可压缩解析壁面函数由于考虑了对流项在边界层内分布，消除了原始解析壁面函数在分离区起始位置的非物理振荡。

3)壁面热流的数值结果差异较大，整体来看本文构造的可压缩解析壁面函数得到的结果和密网格Launder-Sharmak-ε模型结果最为接近，而原始解析壁面函数和标准壁面函数，不能准确预测壁面热流。

4)壁面函数模型使用粗网格，相较于低雷诺数模型可以大幅减少计算时间，对于Ma=5算例，计算时间约为Launder-Sharmak-ε模型的5.3%，而相较与标准壁面函数仅增加了1%左右的时间。

总的来说，四种方法都显示出良好的壁面压力和摩擦系数预测能力。然而，标准壁面函数和原始解析壁面函数严重低估壁面热流，即使在分离区前后其结果也和密网格Launder-Sharmak-ε模型结果相差甚远。本文构造的解析壁面函数充分考虑可压缩湍流边界层流动的特征，弥补了原始解析壁面的缺点，从而大幅提升了壁面热流的预测精度，且不需要在近壁面区域内的精细网格。