高焓湍流边界层壁面摩阻产生机制分析1)

刘朋欣 孙 东 李 辰 郭启龙 袁先旭

(中国空气动力研究与发展中心空气动力学国家重点实验室,四川绵阳,621000)

(中国空气动力研究与发展中心计算空气动力研究所,四川绵阳,621000)

引言

高超声速飞行器在以极高马赫数飞行时,飞行器周围会产生较强的激波并压缩来流空气,使其温度急剧升高.温度的升高会激发复杂的热力学过程,比如振动能激发和气体分子离解,使飞行器表面的气体处于非平衡状态[1].当飞行器在低空稠密大气高速飞行时,飞行器表面流动会发生转捩并形成湍流.湍流会使得壁面摩阻升高数倍[2],损害飞行器飞行性能.湍流与非平衡效应耦合在一起,使得表面摩阻的预测更加困难,高温非平衡效应对壁面摩阻产生机制的影响也尚不清楚.

目前同时考虑高温化学非平衡效应和湍流流动的研究尚处于起步阶段,相关报道还十分有限,且均是采用直接数值模拟方法(direct numerical simulation,DNS)进行分析研究.Duan 和 Martin[3-4]使用5组分空气反应模型进行了DNS,结果表明高焓化学非平衡条件下基于Morkovin 假设的大部分相关尺度关系仍然成立,并评估了湍流对化学反应速率的影响.Kim[5]使用单一组分的O2或N2,并考虑热化学非平衡过程,发现半当地尺度可以较好地消除热化学非平衡效应的影响.刘朋欣等[6-7]发现在边界层的对数区,各流动参数的拓展自相似理论的相对标度指数基本符合标度规律.吴正园等[8]的研究表明湍流边界层中的高温气体效应对壁面平均压强和脉动压强有显著的增加效果.Renzo 和Urzay[9-10]研究了含化学非平衡的边界层转捩过程,并基于模拟结果分析了壁面摩阻、热流的流向演化.Passiatore 等[11]对比了两种不同气体模型(冻结流动和化学非平衡流动)下的湍流边界层,发现吸热反应会使得温度脉动和密度脉动低于冻结流模型,但化学非平衡效应对速度的一阶和二阶统计量影响较小,并简单地采用RD (Renard &Deck)方法对摩阻进行了分解研究.Volpiani[12]开发了一套模拟高超声速激波/化学非平衡湍流边界层相互的直接数值模拟程序,发现化学非平衡效应使得激波干扰区的流动特性发生了显著变化.

上述研究结果表明,高温非平衡效应会显著影响流场热力学特性分布.温度分布的改变又会直接影响黏性系数的大小,继而对摩阻分布产生影响.因此有必要研究高温非平衡效应下的壁面摩阻产生机制.目前高温非平衡湍流边界层的研究大都关注于湍流边界层宏观统计特性规律.而对工程中较为关心的摩阻特性缺乏深入分析.

摩阻分解技术将壁面摩阻的产生与流动机制结合起来,是分析摩阻产生机理的一种有效手段.目前常见的有FIK 分解[13]和RD 分解[14].Fukagata 等[13]通过对平均动量方程的连续三次积分,建立了摩阻系数与雷诺应力之间的关系,将摩阻分解为层流部分的贡献和湍流部分的贡献.后续的改进工作提升了公式的易用性[15-16],并将该公式推广至任意复杂外形[17].Renard 和Deck 指出[14]FIK 公式没有明确的物理意义,并提出了另外一种摩阻分解方法,即RD 分解该方法在假设远离壁面附近流体保持静止的绝对参考坐标系下,通过对流向动量方程做变换和积分,将摩阻的产生分为分子黏性耗散作用和湍动能生成作用,具有更明确的物理意义.此方法在旋涡结构输运过程[18]、大尺度湍流运动[19]和破碎[20]对摩阻的影响研究得到了应用.Li 等[21-22]将其推广至可压缩湍流中,并研究了雷诺数效应对摩阻产生的影响.

本文利用RD 摩阻分解技术,分析了高温化学非平衡湍流边界层摩阻生成机制,给出了对摩阻生成起主要作用的流动过程;通过设置低焓完全气体模型对比算例,从摩阻空间分布特点、热力学和流动特性方面分析了高焓效应对主要摩阻产生项的影响.

1 数值方法与算例设置

1.1 控制方程

低焓算例工质为空气,采用完全气体模型,控制方程为可压缩Naiver-Stokes (N-S);高焓算例由于需要考虑空气的化学非平衡效应,采用含化学反应的N-S 方程作为控制方程.后者的形式如下

其中 ρ和p是密度和压力;uj为分速度;fk是第k种组分的质量分数;E为单位体积的总内能;hk是比焓,通过温度的多项式拟合得到;Sk是第k种组分的生成速率; τij为黏性应力;K和Dk分别为热传导系数和组分质量扩散系数.在完全气体计算中,通常取Pr=0.72,比热比 γ=1.4,黏性系数采用Sutherland 公式计算.而对于多组分气体,则首先通过Lennard-Jones 分子黏性系数计算公式求得各组分的黏性系数,并通过Wilke 公式计算混合物黏性系数; γ 也不再为常数;通过Pr数和S c数与黏性系数的关系式来计算K和Dk.多组分输运系数(黏性系数、热传导系数和组分扩散系数)具体计算公式可参考文献[23].

1.2 数值方法

计算采用课题组发展的Inhouse 程序.该计算程序已经成功开展大量的可压缩湍流[24-25]、激波/边界层干扰[26-27]、旋转爆震精细流场结构[28]、高温非平衡湍流边界层[6-7]等复杂流动的数值模拟,计算精度和鲁棒性均得到了充分验证.在本文计算中,无黏通量采用了7 阶有限差分格式WENO-Z[29],黏性项采用4 阶中心差分格式,时间推进采用具有TVD 性质的3 阶Runge-Kutta 方法.高焓算例计算不考虑温度非平衡效应,且只考虑空气组分的离解过程,化学反应源项求解采用Gupta 等[30]提出的5 组分(N2,O2,N,O,NO) 6 基元反应模型.

1.3 算例设置

高焓算例的计算状态示意图见图1.为了激发空气的化学非平衡态,选择30 km 高空以马赫数20 飞行的楔形体,其头部产生一道斜激波.斜激波后的温度较高且已经足够激发空气发生化学反应,选择此状态作为湍流边界层外缘流动状态.高焓算例来流组分质量分数设置为f(O2)=0.267 5,f(N2)=0.732 5,并使用空气化学反应模型模拟化学非平衡过程,记为TH 算例.作为对比,设置同等边界层参数下的低焓算例,以评估高温非平衡效应对湍流边界层统计特性和摩阻的影响,并记为TL 算例.两算例的来流状态及壁温设置见表1.

表1 来流状态及壁温设置Table 1 Inflow condition and wall temperature

图1 计算模型状态示意图Fig.1 Schematic of computational model

计算初场采用相同来流状态下RANS[31]计算结果的剖面生成.入口的脉动通过数字滤波合成湍流方法[32]产生;出口采用特征边界条件;上边界固定为初始值;壁面采用黏性等温条件;展向采用周期性边界条件.

计算域设置为20 δ × 4 δ × 2 δ (流向 × 法向 × 展向,δ 为边界层名义厚度),相应的网格为901 × 161 ×301.流向和展向网格均匀,法向网格在壁面附近加密.文献[6]进行了同等工况下的网格无关性验证,本文所采用的网格可以满足直接数值模拟的要求.湍流充分发展段的网格分辨率及边界层厚度和雷诺数见表2.此时的边界层厚度 δ 约为5 mm.保持两个算例的Reδ2一致,且可以发现Reθ和Reτ也相差不大.分析所采用流向位置如图2(b)中黑色实线所示.且如无特殊说明,下文中分析默认基于此流向截面处的流场.在下文的表述中约定,代表变量q的Reynolds 平均,表示Favre 平均.两种平均方式所对应的脉动量分别为:.

图2 瞬时流场结构图Fig.2 Instantaneous flow structures

表2 边界层厚度、雷诺数和网格尺度Table 2 Thickness,Reynold number and grid resolution

2 结果分析与讨论

2.1 湍流边界层瞬时与统计结果

图2 给出了瞬时流场旋涡结构和密度梯度分布图.从图中可以看出,在入口添加的扰动,经过一段距离的发展,形成了充分发展的湍流,流动结构变得紊乱无序.从整体上看,高焓算例和低焓算例的流动结构没有明显的变化,说明化学非平衡过程对流动运动过程影响较小,并不会显著改变流场结构特性.

图3 和图4 分别给出了平均流向速度和雷诺应力分布.Van Driest 变换的定义为:.从图3 中可以看出,此时速度边界层仍然存在明显的线性区、缓冲区和对数区,只不过截距有所增大,与文献中高焓算例[3]结果符合较好.除了在缓冲区,高焓算例的平均流向速度略高于低焓算例之外,焓值并不会明显改变湍流边界层速度剖面的整体分布特性.从图4 中可以看出,经过密度加权的雷诺应力分布具有相似的分布趋势,与文献完全气体结果[33]及高焓算例结果[3]符合较好.高焓效应对雷诺应力分布的影响不大,并没有对湍流速度脉动产生明显影响,这与图2 中的瞬时流场结构对比结果是一致的.

图3 Van Driest 变换后的速度分布Fig.3 Profile of Van Driest transformed velocity

图4 雷诺应力分布Fig.4 Reynolds stress distribution

2.2 湍流边界层摩阻空间分布

本节首先对比高低焓值下摩阻的空间分布,分析高温非平衡效应对摩阻分布的整体影响.图5 中给出了不同焓值工况下的局部瞬时摩阻系数在壁面上的分布,从图中可以看出瞬时摩阻的分布是相似的,都呈现出局部高摩阻斑点,并且存在一些摩阻较弱的条带结构.另外,高焓情况下的条带结构尺度有减小趋势.

图5 瞬时壁面摩阻系数分布Fig.5 Instantaneous distribution of skin friction coefficients

图6 给出了摩阻脉动的概率分布函数(probability distribution function,PDF),使用其均方根进行归一化,并与高斯分布函数进行对比.结果表明,高焓非平衡流动和低焓流动的PDF 曲线均呈现出倾斜趋势;两者分布基本一致且均偏离高斯分布的中心位置,峰值位置约在-0.5 处;高焓算例的峰值略低于低焓算例,但两者的峰值概率密度略高于高斯分布;计算结果与文献[34]完全气体模型(Reθ=2300)较为一致.这说明摩阻脉动的PDF 受化学非平衡效应的影响有限.PDF 曲线向负方向一侧倾斜,说明当地低速条带发生的概率较大;＞0 对应于高速条带流动事件,且摩阻脉动的幅值达到4 左右,表明高速条带产生的脉动更强.

图6 摩阻脉动PDF 分布Fig.6 PDF profile of skin friction fluctuations

图5 的结果表明,摩阻脉动的高低值条带在尺寸上存在差异.进一步可以通过摩阻脉动的空间两点相关系数来定量分析高温非平衡效应对摩阻脉动条带尺度的影响.壁面上相关系数定义为

其中代表摩阻脉动,Δx,Δz分别代表流场和展向间距,(x0,z0) 在所选流向站位的展向中间位置.式(3)是指摩阻脉动的空间两点相关系数值等于0.2 的流向尺度和展向尺度.

图7 给出了相关系数的分布云图.不同焓值下的相关系数均呈现出椭圆形分布,且流向尺度远大于展向尺度.高焓流动的分布范围在流向和展向两个方向上均小于低焓流动.这与图5 中瞬时摩阻脉动的条带尺寸一致.沿两个方向进行长度尺度积分[34],可以得到定量的尺寸,具体见式(3).低焓流动的流向尺度为1.41δ,展向尺度为0.12δ;高焓流动的流向尺度为1.20δ,展向尺度为0.11δ.可见高焓算例的流向尺度减小了14.9%,展向尺度减小了8.3%.

图7 壁面摩阻脉动两点相关系数空间分布

进一步,采用壁面黏性尺度δv进行对比.低焓流动的流向尺度为1130.6δv,展向尺度为96.0δv;高焓流动的流向尺度为993.5δv,展向尺度为95.5δv.可见高焓算例的流向尺度减小了12.1%,展向尺寸基本一致.

高焓条件下摩阻条带结构尺度减小的原因可能来自于高温非平衡效应.热力学量的脉动引起局部流向速度梯度增大,如图3 中所示在缓冲区较高,条带结构的能量更高.相较低焓流动,输运相同能量,高焓边界层条带结构更小.具体原因还需要进一步验证.

2.3 湍流边界层摩阻分解

研究摩阻的产生机理以及与之相关的近壁区动力学过程,对于基础研究和工程实际都有重要意义,可以为减阻方法的设计提供指导.对于可压缩壁湍流,假设展向均匀,可将壁面摩阻写成如下分解积分关系式[22]

其中 ρ∞和U∞分别为来流的密度和速度,本算例中分别等于 ρe和Ue.Cf,V代表分子黏性耗散项;Cf,T代表湍动能生成项;Cf,G1表示平均流向动能的变化项;Cf,G2表示流向黏性正应力不均匀性项,Cf,G3表示流向雷诺正应力不均匀性项,两者合起来称为流向不均匀项.其中.

在对湍流边界层摩阻空间分布分析的基础上,本节继续采用摩阻分解式(4)来分析边界层中摩阻生成机制.表3 比较了采用摩阻分解技术得到的摩阻平均值(Cf,decomp)、直接采用时间平均得到的摩阻值(Cf,avg)以及两者之间的相对误差.相对误差的定义为error=(Cf,decomp-Cf,avg)/Cf,avg.结果显示,本文实现的摩阻分解具有较高的精度,相对误差均在1.5%以下.

表3 摩阻分解公式相对误差Table 3 Relative error of decomposition formula

表4 给出了式(4)中不同流动和输运过程对摩阻产生的贡献比例.可以看出,分子黏性耗散项(Cf,V)和湍动能生成项(Cf,T) 所占比例最大,基本上在40%以上;其次是平均流向动能的变化项(Cf,G1),比例在10%;流向黏性正应力不均匀性项(Cf,G2)和流向雷诺正应力不均匀性项(Cf,G3)占比很小,基本上可以忽略.对比不同焓值算例可以发现,高焓值情况下,Cf,V和Cf,T的占比略有下降,减少的这部分贡献基本上由Cf,G1补上.摩阻分解分析表明,湍流边界层中对摩阻生成起主要作用的流动过程是Cf,V和Cf,T.下面将进一步分析两者产生的原因以及高低焓值边界层摩阻产生差异的原因.

表4 摩阻分解公式中各项贡献Table 4 Contributions of different terms to skin friction

2.4 分子黏性耗散项对壁面摩阻的贡献分析

分子黏性耗散过程是壁面摩阻产生的主要流动过程之一.图8 给出了分子黏性耗散项预乘积分函数沿法向的分布.两算例的分布趋势较为一致,约在y≈0.015δ(y+≈10)处达到峰值.且分子黏性耗散项的主要作用区域在y≤0.1δ(y+≤30)的近壁区,而在外区基本没有贡献.这与文献中基于槽道、平板湍流[21-22]的分析结果一致.高焓会使得y＜0.015δ 以内的近壁区分子黏性耗散略大于低焓流动;而在y＞0.015δ的区域中略小于低焓流动.

图8 黏性耗散项预乘积分函数Fig.8 Pre-multiplied integrand of viscous dissipation terms

进一步可以将分子黏性耗散项分解为平均部分Cf,V,m和脉动部分Cf,V,f[21]

表4 中给出了两部分的占比.可以发现,Cf,V,m占据了分子黏性耗散项的绝大部分,占比高达95%以上.这就说明焓值对分子黏性耗散项的影响主要通过Cf,V,m项体现.

将Cf,V,m的积分函数分为两部分:和,分别体现了流动特性和热力学特性的作用.图9给出了两者沿壁面法向的分布,并给出了温度的分布,均采用无穷远处参数作为无量纲参考值.从流动特性可以看出,高焓使得近壁区(y＜0.015δ)的∂u/∂y增大(∂v/∂x较小,影响可以忽略);而对外区的流动特性影响很小.从热力学特性可以看出,高焓对整个边界层内的黏性系数都有较大影响,基本上均低于低焓算例的黏性系数.黏性系数主要受到温度的影响,两者的分布具有明显的一致性.当高焓时,边界层内发生空气化学非平衡过程,离解反应会吸收一部分热量,使得高焓边界层内的黏性系数低于低焓边界层.综合流动特性和热力学特性的分布,可以发现,在y＜0.015δ 时,流动特性起到主导作用,导致高焓流动黏性耗散项贡献略大于低焓流动;而在y＞0.015δ 时,热力学特性起到主要作用,导致低焓流动黏性耗散项贡献略大于高焓流动.

图9 Cf,V,m积分函数分解及温度分布Fig.9 Integrand parts of Cf,V,mand the temperature distribution

2.5 湍动能生成项对壁面摩阻的贡献分析

湍动能生成项是壁面摩阻产生的另一个重要部分.图10 给出了湍动能生成项预乘积分函数沿法向的分布.两算例的分布趋势较为一致;但低焓流动时,峰值位置更加远离壁面.结合图4 中雷诺应力分布基本一致的情况来看,峰值位置的外移应该来自于的分布不同;结合图9 来看,高焓流动在预乘积分函数附近的要大于低焓流动.两算例在边界层外区均出现了二次峰值,这与较高雷诺数下的槽道流动[21]一致.二次峰值的产生是由于较高雷诺数下,边界层外区存在大尺度湍流运动导致.大尺度湍流运动会引起外区产生大量湍动能,从而通过摩阻的湍动能生成项影响摩阻;且外区大尺度湍流流动对内区小尺度湍流运动具有调制作用,促发湍动能在内外区的重新分配[35].

图10 湍动能生成项预乘积分函数Fig.10 Pre-multiplied integrand of TKE production terms

为了进一步分析湍动能生成项中起主导作用的流动事件,本节采用象限分析技术[36-37]将边界层中的主导流动事件与摩阻分解相结合以分析不同流动事件对摩阻中湍动能生成项的影响.

在边界层中,可以按照脉动的正负将流向速度脉动和法向速度脉动分解到4 个象限中.第1 象限:u′′＞0,v′′＞0表示外向运动(outward motion);第2 象限:u′′＜0,v′′＞0表示上抛运动;第3 象限:u′′＜0,v′′＜0表示内向运动(inward motion);第4 象限:u′′＞0,v′′＜0表示下扫运动.图11 给出了Q1～Q4 事件在边界层内发生的概率.不同焓值下不同事件的概率分布曲线是相似的.Q2 和Q4 发生的概率显著大于Q1 和Q3.

图11 Q1～Q4 事件发生概率Fig.11 Probability distribution of Q1～Q4 events

图12 中给出了不同事件对应的雷诺应力与总雷诺切应力的比值.可以看出,不同焓值的比值曲线很好的重合在一起,说明该比值对焓值变化不敏感;Q2 和Q4 事件是雷诺切应力的主要贡献项,并且在y≈0.03δ以内,Q4 事件的比值要大于Q2 事件;在y≈0.03δ以外,Q2 事件的比值大于Q4 事件,且这一趋势会在整个边界层范围内得到保持.

表5 定量给出了不同事件的湍动能生成项贡献的摩阻与整体摩阻的比值,表中以Cf,T,Qi表示Qi事件对湍动能生成项的贡献.可以看出,Q2 和Q4 事件起到绝对主导作用,且Q2 贡献要略大于Q4 事件,这与图12 中所指出的Q2 对雷诺应力的贡献大于Q4 一致.且Q1 和Q3 对壁面摩阻的贡献为负,两者的值基本相同,且都远小于Q2 和Q4.

表5 Q1～Q4 湍流生成项对壁面摩阻的贡献Table 5 Contributions of TKE production terms produced by Q1～Q4 to the time-averaged value

图12 Q1～Q4 事件对湍动能生成项的贡献Fig.12 Contributions of Q1～Q4 events to TKE production terms

图13 给出了不同事件的被积函数沿法向的分布,以考察主要起作用的边界层区域范围.可以看出,不同焓值下,被积函数的曲线分布趋势都是相似的,高焓效应仅会稍微改变幅值;Q4 事件的峰值位置要比Q2 事件更为靠内;Q1 和Q3 在整个边界层范围内的值都是负值.

图13 Q1～Q4 湍流生成项预乘积分函数Fig.13 Pre-multiplied integrand of TKE production terms produced by Q1～Q4

3 结论

本文通过采用直接数值模拟方法,研究了高焓湍流边界层壁面摩阻产生机制,并通过设置同等边界层参数下的低焓算例,对比分析了高温化学非平衡效应对摩阻产生的影响,可以得到以下结论.

(1)不同焓值下的摩阻脉动条带尺寸存在差异.高焓流动的条带尺寸在流向和展向两个方向上均小于低焓流动.

(2)湍流边界层中对摩阻生成起主要作用的流动过程是分子黏性耗散项Cf,V和湍动能生成项Cf,T.

(3) 分子黏性耗散项Cf,V的主要作用区域在y≤0.1δ的近壁区.焓值的影响主要通过Cf,V的平均部分Cf,V,m项体现.且在y＜0.015δ 时,流动特性(速度梯度)起到主导作用,导致高焓流动Cf,V,m略大于低焓流动;而在y＞0.015δ 时,热力学特性(黏性系数)起到主要作用,导致低焓流动Cf,V,m略大于高焓流动.

(4)边界层中的上抛和下扫运动是影响摩阻分解湍动能生成项Cf,T的主导事件,外向运动和内向运动对Cf,T的影响很小,且起到负贡献.