分子黏性对解析湍流壁面函数的影响

王新光, 陈琦, 万钊, 高晓成, 燕振国

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

0 引言

高超声速湍流边界层模拟时,如需精确模拟壁面摩阻和热流,要求壁面网格足够密,通常无量纲壁面距离y+≈1,导致壁面附近网格骤增,使得计算收敛速度变慢,同时还会影响数值计算稳定性。因此,主流商业软件例如Fluent和CFD++,对工程外形进行湍流模拟时通常使用壁面函数,大幅放宽壁面网格尺度,以提高计算效率。

壁面函数最初基于Prandtl混合长度理论和若干假设,得到不可压缩湍流边界层内的相似解[1],在湍流计算中通过引入壁面函数达到提高粗网格气动力预测精度的目的[2]。在标准壁面函数的基础上,文献[3]通过对摩擦速度构造局部迭代实现壁面函数在超声速流动中的应用,其中在复杂外形气动力预测中可有效提高粗网格的预测精度。此外,通过添加压力梯度[4]、旋转修正[5]、可压缩性和壁面传热效应[5]等影响,进一步扩展了壁面函数在可压缩湍流边界层中的应用。文献[6]通过使用Crocco-Busemann温度方程,添加可压缩性对壁面函数的影响,并将速度和温度壁面函数耦合求解,发展了适用于可压缩湍流边界层的壁面函数。国内学者[7-8]将这种壁面函数耦合到了k-ω(k表示湍动能,ω为比耗散率)两方程模型中,对于无量纲壁面距离y+<100范围,取得满意的结果。文献[9]比较了两方程k-ωSST、线性Launder-Sharmak-ε(ε表示湍动能耗散率)模型(LS)k-ε和非线性k-ε模型在超声速和高超声速流动中的应用,比较发现上述k-ω和k-ε模型的结果是类似的。文献[10-11]基于文献[6]中的壁面函数,使用数值实验和风洞实验数据开展了修正研究,通过对速度、温度的壁面函数修正来提高其预测精度。目前大多壁面函数研究的适用性仅针对湍流边界层,对于存在分离和再附等逆压梯度的复杂流动,适用性仍无法确定[12-14]。

近年来发展的解析壁面函数[15]对于存在分离的流动表现良好[16-18]。文献[19]考虑壁面网格内对流项变化和能量方程中黏性耗散项的影响,发展了适用于可压缩流动解析壁面函数(MAWF),通过二维超声速激波边界层干扰算例进行验证,数值结果表明发展的可压缩修正解析壁面函数消除了原始解析壁面函数的非物理振荡,且大幅提升了壁面函数壁面热流的预测精度,接近密网格低雷诺数模型结果,且可节约大量计算时间,文中马赫数为5的算例其计算时间仅为密网格的5%。

可压缩湍流边界层直接数值模拟(DNS)结果[20]和标准壁面函数相比,当来流马赫数较小时标准壁面函数依然适用,与DNS结果吻合,但随着马赫数的逐渐增大,标准壁面函数和DNS之间的偏差逐渐增大,表明低雷诺数超声速槽道湍流存在显著的可压缩性。对于Ma为2.48的槽道流动,在黏性底层边缘附近y+≈10,相较于壁面参数,密度减小50%,温度增高60%。通常黏性底层内温度的快速变化必然导致流体黏性系数的变化,这一点对于壁面函数尤为重要。

本文基于发展的可压缩解析壁面函数,通过构造不同的黏性系数方程,考虑密度、黏性系数等可压缩流动参数对无量纲壁面距离的影响,构造了两种不同黏性系数分布,发展了两种适应于高超声速流动的解析壁面函数,以达到提高解析壁面函数粗网格壁面热流预测精度的研究目的,并通过高超声速激波边界层干扰算例,比较不同黏性系数方程对数值模拟结果的影响,为解析壁面函数的工程化应用起到抛砖引玉的作用。

1 可压缩流动无量纲壁面距离

解析壁面函数[15,19]理论中,无量纲壁面距离的定义对于解析壁面函数至关重要。原始的解析壁面函数中无量纲壁面距离的定义为

(1)

式中:ρw表示壁面密度;y为壁面距离;kP为湍动能,下标P表示壁面第1层网格点;μw表示壁面黏性系数。

随着马赫数的增加,黏性底层内温度、密度等参数发生很大变化,文献[20]中通过对可压缩湍流边界层的DNS结果分析,认为可压缩流动中无量纲壁面距离采用当地值相较于壁面值更合适。但是壁面函数使用粗网格,壁面第1层网格内通常包含黏性底层和全湍流层,因此一种方式是使用主网格计算得到的当地P点的值定义无量纲壁面距离:

(2)

式中:ρP表示壁面第1层网格点密度;μP表示壁面第1层网格点黏性系数。

数值实验结果表明,这种定义方式在计算激波边界层干扰算例时稳定性较差,为避免类似的问题出现,本文采用黏性底层边缘处的值定义,即

(3)

式中:ρv表示黏性底层位置的边缘密度;μv表示黏性底层位置的黏性系数。

黏性底层处的密度和黏性系数使用上一时间步的解析温度Tv计算,即后文式(9)在黏性底层边缘处的值计算。使用这种方式定义的无量纲壁面距离避免了数值不稳定性,同时考虑了壁面网格内流动参数变化的影响。

2 分子黏性系数建模

随着来流马赫数的增加,可压缩湍流边界层内温度梯度增加(见图1),图1中δ为边界层厚度,T为温度,Tw为壁面温度,Ma为来流马赫数,T∞为来流温度,湍流边界层来流条件可参考文献[21-23]。而发展的可压缩壁面函数[19],仅考虑壁面网格内湍流黏性系数的变化,可能导致解析壁面函数对于高超声速流动的预测精度降低。

图1 可压缩平板湍流边界层温度分布

对于完全气体,黏性系数仅是温度的函数,但直接将温度解析表达式和黏性系数通过Sutherland公式连接后,很难得到解析的速度和温度表达式。本文通过解析壁面函数得到黏性底层的温度Tv,通过Sutherland公式得到黏性底层位置的黏性系数μv,并与壁面黏性系数关联,构造黏性底层内的黏性系数表达式。其中一种简单的思路是与壁面处μw线性连接,如图2(a)所示。图2中,N为壁面第2个点,P为壁面第1个点,S为壁面中心,μN为N点的黏性系数,yv为黏性底层的距离。

图2 分子黏性系数分布示意图

则黏性底层内的黏性系数为

(4)

当使用式(4)代入解析壁面函数中时,会使得解析温度的表达式更加复杂,并使得解析壁面函数的稳定性变差[24]。考虑到上述问题,文献[24]对不同类型的黏性系数表达式进行了探索研究,其中两种方式的黏性系数表达式在数学上接近于线性表达式(见式(1)),但不存在线性表达式的稳定性问题,如图2(b)和图2(c)所示。

图2(b)表示双曲型分布:

(5)

图2(c)表示抛物型分布:

(6)

在全湍流区使用μv,如图2所示。在实际程序中仅需保存上一时间步的解析温度Tv。

图3给出了用式(5)和式(6)构造的解析壁面函数模拟Ma=8.18的激波边界层算例(来流条件见表1,其中β为激波产生器角度,θ0为动量边界层厚度,p∞为来流压力)。干扰区前后壁面网格内温度和分子黏性系数分布,并与密网格LS结果进行了对比,其中粗网格第1层网格高度为1.2 mm。由图3可知,激波边界层干扰区后黏性底层法向距离变小,且黏性底层内温度变化梯度更大,而全湍流区的温度变化较小,对应的分子黏性系数,在黏性底层边缘大约较壁面值增加了30%,其中采用式(5)的双曲型分布时,分子黏性系数呈上凸形式,采用式(6)的抛物型分布时,呈下凹形式,且更接近于密网格低雷诺数模型结果,而壁面网格内两种分布形式解析温度分布差异很小,几乎可以忽略,温度分布结果接近于密网格结果,而文献[19]中发展的MAWF其温度分布和密网格有明显差异。

表1 Ma=8.18时2D斜激波边界层干扰来流条件

图3 Ma=8.18算例粗网格解析温度(左)和分子黏性系数(右)分布与密网格低雷诺数模型结果对比

3 双曲型分布可压缩解析壁面函数

解析壁面函数通过对边界层内简化的动量方程和能量方程进行积分,分别得到边界层内速度和温度表达式[17-19],当黏性系数采用双曲型分布(见式(5))时,通过对简化的动量方程进行两次积分,可得到速度的解析表达式如下:

(7)

(8)

式中:常数α=0.229 5;下标1表示黏性底层(即y

vp为第1层网格点y方向的速度,p表示压力;系数A1、A2、B2分别为

注意到本文发展的解析壁面函数公式中bμ=0时,将不考虑黏性系数的分布,使用黏性底层的μv值来代替壁面值μw,其公式与可压缩解析壁面函数[18]相同。

对简化的能量方程进行类似的积分,可得到解析温度表达式为

(9)

(10)

式中:Pr为普朗特数;cp为定压比热,cp=1 004.06 J/(kg·K);αt为常数0.193 6;Dth1、Dth2为能量方程的对流项系数,

注意到解析温度表达式(式(9)和式(10))和解析速度表达式(式(7)和式(8))耦合在一起,因此在程序中使用数值积分计算。这里仅给出双曲型分布的最终表达式,抛物型的最终解析表达式可参考文献[25]。

4 数值结果与分析

本文使用不同来流马赫数的激波/湍流边界层干扰算例来验证考虑可压缩流体边界层参数变化的解析壁面函数、MAWF,与本文发展的双曲型(hyper-MAWF)或抛物型(para-MAWF)分布解析壁面函数进行对比。

数值计算中无黏通量使用Roe-Pike格式,黏性通量采用2阶中心差分格式。粗网格第1层网格y+大约为30,湍流模型使用标准k-ε模型,壁面模型采用上文描述的3种壁面函数。密网格第1层网格的y+≈1,湍流模型采用添加Yap修正[26]的LSk-ε模型。来流入口给定充分发展的湍流边界层,具有1.5%的湍流度和黏性系数比率μT/μ=10。本文算例网格无关性研究表明计算结果与网格无关,详情可参考文献[25],其中Ma=8.18算例密网格和粗网格网格无关性结果如图4所示,其中Q表示热流,Qw表示平板壁面热流。

图4 Ma=8.18算例网格无关性研究

4.1 Ma=8.18斜激波边界层干扰

斜激波边界层干扰作为超声速飞行器唇口处的简化模型,斜激波以不同角度入射平板,会产生激波反射、分离和再附等复杂流动现象。Ma=8.18斜激波边界层干扰[23],在激波边界层干扰区域前局部来流条件如表1所示。由于激波产生器产生的斜激波和拐角的膨胀波对壁面均有干扰,在计算中上边界给定激波产生器的壁面边界模拟入射激波和激波产生器拐角的膨胀波,计算网格如图5所示。

图5 Ma=8.18算例网格示意图

图6给出了密网格480×210使用LS模型与粗网格160×60壁面函数方法湍流边界层内的速度和温度分布,可见数值解和实验值基本吻合,其中密网格结果和实验值吻合较好,在湍流边界层内无量纲温度呈现先减小后增大的趋势,粗网格在边界层内网格点较少,仅有10个点,大致可以刻画出边界层内的温度变化。

图6 Ma=8.18算例边界层厚度0.94 mm位置边界层内速度和温度分布

图7给出了β=10°时算例不同位置近壁面解析温度与密网格LS模型的对比,其中x分别取值25 cm和55 cm时,截面位置位于激波边界层干扰区前后,为全湍流边界层。从图7中可知,MAWF温度梯度明显较大,表明MAWF计算得到的壁面热流值相对较大,在图9壁面热流分布中可对照观察到。当x分别取值35 cm和45 cm时,即激波边界层干扰区内,此时近壁面温度梯度较大,本文设计的两种壁面函数即hyper-MAWF和para-MAWF,温度梯度相对其他两种方法较小,对应图9壁面热流在分离区内较小。

图7 Ma=8.18,β=10°算例解析温度分布

图8给出了β=10°算例不同位置近壁面解析速度分布,其中在激波干扰前后位置x=25 cm和x=55 cm处,完全湍流区处速度分布有较明显的差异,其中MAWF解析速度梯度明显较其他数值结果大,虽然图6中边界层4种数值方法计算得到的速度分布大致与实验相同,但由于解析速度在壁面处的梯度决定了摩阻的数值,即MAWF计算得到更大的摩阻,而本文设计的两种壁面函数其解析速度分布几乎重合在一起,整体速度梯度较LS和MAWF较小。

图8 Ma=8.18,β=10°算例解析速度分布

图9比较了Ma=8.18算例不同激波产生器壁面压力、摩擦阻力Cf和热流结果(其中Q∞表示来流热流),相较于MAWF,考虑流体性质的解析壁面函数得到的壁面压力在分离区略有下降,壁面热流在分离区显著降低,更接近于实验值,两种不同分布形式黏性系数对应的壁面函数也出现差异,其中使用双曲型分布得到的结果更低一些,相较于抛物型分布低约4%,而相较于MAWF降低了约35%。密网格LS结果出现高估壁面热流现象,与文献[27]中密网格低雷诺数模型结论一致。壁面摩阻和壁面热流有相似的分布,除MAWF外,其他数值计算出现一个较小的分离区(摩阻为负),且在整个干扰区内本文设计的两种壁面函数计算得到的最大摩阻比LS结果小约30%,有较大差异,但本文实验没有对应的摩阻试验结果,暂时还无法确定壁面摩阻数值计算的精度。

图9 Ma=8.18算例激波产生器5°壁面压力(上)、摩阻(中)和热流(下)分布

4.2 Ma=9.22压缩拐角

压缩拐角流动常见于高超声速飞行器控制舵面、翼面和飞行器表面之间,本文选取Ma=9.22压缩拐角[28],计算网格较为简单,入口给定充分发展的湍流边界层,来流参数如表2所示,其中Re表示雷诺数。其中低雷诺数模型使用的密网格为240×150,壁面函数使用的粗网格为180×80。其中密网格压缩拐角角度β分别取值15°和32°时的马赫数云图如图10所示,从中可知随着马赫数的增大,压缩拐角产生的斜激波逐渐增强,与来流湍流边界层之间的干扰加剧。

表2 Ma=9.22压缩拐角来流条件

图10 Ma=9.22算例马赫数云图

图11比较了压缩拐角分别为15°、32°、34°和38°时不同模型的壁面热流分布,其中横坐标S表示距离拐角处的位置(cm),低雷诺数LSk-ε模型高估了壁面热流,当压缩拐角38°时最大热流几乎是实验值的3倍。文献[29]使用LSk-ε模型和Wilcoxk-ω模型,文献[30]使用Rodik-ε模型模拟了34°压缩拐角。这3种低雷诺数模型均严重高估了壁面热流值。MAWF也返回较高热流值,在压缩拐角上接近于LS模型结果,均高于实验值,而本文发展的两种壁面函数均给出合理的壁面热流值,更接近于试验结果,两种黏性系数数值结果差异较小,在5%之内,且双曲型分布更接近于实验值。

图11 Ma=9.22压缩拐角壁面热流

5 结论

本文基于高超声速湍流边界层特征,使用解析温度计算黏性底层边缘密度和黏性系数,重新定义了无量纲壁面距离,并构造了双曲型和抛物型两种黏性系数方程,将黏性底层内快速的速度梯度变化考虑到解析壁面函数中,构造了抛物型解析壁面函数para-MAWF和双曲型解析壁面函数hyper-MAWF,通过高超声速边界层干扰算例对构造的壁面函数进行了测试,并与低雷诺数LS和MAWF数值结果进行了对比。得出主要结论如下:

1)对于高超声速算例,本文构造的解析壁面函数预测的壁面热流结果更接近试验结果,密网格LS模型和MAWF数值结果在激波边界层干扰区均出现高估热流现象。

2)本文构造的两种不同黏性系数表达式,数值模拟结果差异较小,在5%之内,其中双曲型分布预测的壁面热流更接近于实验值。

3)本文构造的解析壁面函数预测的速度、温度和黏性系数分布更接近于密网格结果,而MAWF未考虑黏性系数分布,因此黏性系数和密网格之间差异较大。

整体来看,本文构造的解析壁面函数相较于MAWF壁面热流的预测结果更加精确,在干扰区甚至显示出优于密网格LS模型的预测结果,且壁面函数采用粗网格进行数值模拟,可显著节约计算时间[19]。由于两种壁面函数数值计算结果差异较小,且形式上双曲型更为简单,对于高超声速算例更推荐使用双曲型分布的解析壁面函数。