激波/湍流边界层干扰分离泡直接数值模拟

2022-09-05童福林董思卫段俊亦李新亮

航空学报 2022年7期

童福林，董思卫，段俊亦，李新亮

1. 中国空气动力研究与发展中心空气动力学国家重点实验室，绵阳 621000 2. 中国科学院力学研究所高温气体动力学国家重点实验室，北京 100190 3. 中国空气动力研究与发展中心计算空气动力研究所，绵阳 621000 4. 中国科学院大学工程科学学院，北京 100049

数十年来，激波与湍流边界层的相互作用一直是工程界和学术界都十分关注的热点问题和基础性问题，至今仍未得到彻底解决。在强激波的诱导作用下，干扰区内会出现大范围的流动分离现象，并形成复杂的波系结构。大量研究表明，对于无分离或小分离边界层流动，现有湍流模型往往有较好的预测精度，数值结果与试验测量数据吻合较好。而对于存在大范围分离流的激波/湍流边界层干扰问题，工程上常用的湍流模型往往不能很好地重现干扰区内某些重要的流动现象，例如多数模型能够较为准确的预测出再压缩过程后的峰值热流和压力，但采用不同湍流模型往往给出了不同的分离区长度和分离过程中的压升结果。因此，进一步深入研究干扰区内分离泡的典型特征非常有助于加深对该问题的理解认识，同时也可为湍流模型的改进提供重要的理论依据。

国内外学者针对压缩拐角和平板入射激波这两类经典构型开展了大量的风洞试验和高精度数值模拟研究，特别是在初始分离准则和分离区长度的准确预测方面。Roshko等通过试验测量数据发现，大尺寸分离区出现的临界角随雷诺数增长具有增大的趋势，而Settles等的试验结果进一步证实了，对于完全湍流情况，该临界角与雷诺数无关。Zheltovodov等在总结归纳大量试验数据和数值模拟结果的基础上，给出了绝热壁条件下考虑了雷诺数、马赫数、拐角角度和边界层厚度影响的分离区长度预测公式。随后，Jaunet等通过试验还发现壁面加热对分离判据的影响较小。Zhu等通过直接数值模拟(Direct Numerical Simulation，DNS)研究了壁温对分离泡大小的影响规律，发现随着壁温增加，分离泡变大。他们认为壁温的影响主要体现在对近壁雷诺数的改变上，并由此建立了分离泡大小与壁温的0.85 次方成正比的半理论模型。童福林等对入射激波与超声速湍流边界层干扰进行了一系列的DNS研究，分析了强膨胀作用下分离区长度的变化规律，给出了考虑膨胀效应的分离区长度归一化公式。

近些年，分离泡空间结构特征及其非定常运动特性越来越受到重视。Grilli等采用动态模态分解方法对其大涡模拟数据进行模态分析，利用4个低频模态成功重构了压缩拐角内分离泡的膨胀和收缩。Priebe等数值研究了24°压缩拐角分离激波低频振荡现象与分离泡的关联性。低通滤波结果表明，低频振荡运动与下游分离泡的舒张和收缩运动密切相关。然而，以往研究大多是针对展向平均后的二维流场开展的，并没有考虑分离泡展向三维结构的影响。实际上，现有风洞试验和数值模拟结果都证实，即便是在压缩拐角和平板入射激波这类准二维构型下，分离区内流动仍存在显著的三维特征，如分离流线发散和汇聚、回流斑块和类“Owl eye”结构等。

本文采用DNS方法对入射激波/平板湍流边界层干扰区内分离泡进行数值研究。通过定量比较3个不同展向站位分离泡的典型特征，细致探讨展向三维结构对分离泡非定常运动特性的影响规律，条件统计分析分离泡内瞬态分离微团的几何形态。此外，采用本征正交分解(Proper Orthogonal Decomposition,POD)方法，比较分析各展向站位脉动速度场的差异，同时对分离泡的非定常运动过程进行低维重构。为了便于比较和验证结果，计算参数的选取与Dupont等的试验和Fang等的DNS相近。

1 计算设置

控制方程为三维曲线坐标系(,,)守恒型可压缩Navier-Stokes方程组:

(1)

式中：守恒变量，无黏通量、和，黏性通量、和的具体表达式见文献[18]。采用自由来流参数以及单位特征长度对方程进行无量纲化，黏性系数的计算采用Sutherland公式。DNS计算时，求解器为课题组自主开发的高精有限差分软件OpenCFD-SC，该程序已在压缩拐角、平板入射激波、超声速膨胀角等复杂流动中得到广泛应用，可以保证计算结果的准确性和可靠性。黏性项的计算采用WENO_SYMBO_LMT格式以及Steger-Warming流通量分裂方法，无黏项的计算采用八阶中心差分格式，时间推进采用三阶Runge-Kutta方法，详细的数值方法介绍见文献[18, 22]。

图1给出计算模型示意图。模型流向长度=137.6 mm，法向高度为=12.7 mm，展向宽度为=4.4 mm，来流方向为从左往右。如图1所示，计算域入口处层流边界层在壁面吹吸扰动(Blowing and suction)作用下转捩生成非定常湍流边界层，随后在计算域上边界数值生成一道入射斜激波，进而与湍流边界层产生相互作用。这里入射斜激波波角取为33.2，上边界层激波入射点取为=82.2 mm，为入射激波在壁面上的名义入射点。来流马赫数为=2.25，来流静温为=169.44 K，基于单位米的来流雷诺数为=2.5×10。图1中:为上游湍流边界层参考点;为层流边界层厚度。

图1 计算模型示意图Fig.1 Computational model

图2 计算网格示意图Fig.2 Computational grid

具体边界条件如下：首先，在计算域入口取为相同来流条件下的层流解，层流边界层厚度为=0.6 mm。计算域出口采用超声速出口条件及缓冲区。计算域上边界采用简单无反射边界条件，同时在激波入射点前后分别设置为自由来流参数和按照Rankine-Hugoniot关系式给出波后参数。计算域下边界采用无滑移等温壁条件，壁温为=321.9 K，吹吸扰动带的起始位置分别=7.62 mm和=20.32 mm，多频扰动波的分布函数和扰动形式与Fang等的DNS完全一致。由于本文转捩区相对较短，为了匹配上游参考点处湍流边界层参数(见表1中边界层厚度雷诺数、位移厚度雷诺数*和动量厚度雷诺数)，扰动幅值和扰动基频取为=0.2和=0.628/。在本文分析讨论中，为来流速度，为参考点处的边界层厚度;、、分别对应流向、法向、展向上的速度。

表1 参考点xref处湍流边界层参数

2 结果验证

在流场达到统计定常后(约两个无量纲时间/)，对三维瞬态流场开始进行统计取样，共获得600个流场样本，统计时间为/=750(为采用时间)。

图3给出了上游层流边界层的转捩过程，这里采用物面法向距离渲染的瞬态密度梯度等值面云图显示。在壁面吹吸扰动作用下，层流边界层间歇性显著增强，同时边界层外层出现了大尺度湍流凸块(Bulges)。图4给出了参考点处的平均速度剖面和雷诺应力分布情况。本文中平均指的是时间和展向平均。可以看到，计算结果与Bookey等的试验数据以及Fang等的DNS结果吻合良好。从图4(b)中还可以看到，雷诺正应力峰值出现在近壁区=13处，这与Pirozzoli等的研究结果也较为一致。

图5还给出了干扰区内无量纲化平均物面压力=(-)(-)的分布情况，其中和分别为物面压力和来流压力;为入射激波波后压力值。与Fang等相同，图中流向坐标=(-)，其中为=(+)/2对应的流向位置。在强激波的干扰作用下，物面压力急剧升高，随后在干扰区下游逐步逼近于无黏解。计算得到的压力分布与Dupont等的试验结果和Fang等的DNS结果均符合较好，进一步证实了本文计算结果的可靠性。

图3 上游边界层转捩过程瞬态密度云图Fig.3 Instantaneous density gradient in upstream transitional boundary layer

图4 参考点xref处湍流统计特征Fig.4 Turbulence statistics at reference station xref

图5 物面压力分布Fig.5 Wall pressure distribution

3 平均结构

图6为干扰区内时间平均分离泡云图，这里采用速度等值面=0进行显示。图中:和代表平均分离点，和代表平均再附点。图7 为物面平均摩阻分布及不同展向站位截面内分离泡高度的比较情况。这里沿展向分别取了3个站位：=0.7，=1.1和=1.75，其中站位对应为展向中截面。各站位截面内分离泡高度通过速度等值线=0的法向位置来确定。

图6 时间平均分离泡高度云图Fig.6 Contour of mean separation bubble height

图7 平均摩阻及不同展向截面内分离泡高度分布Fig.7 Distribution of mean skin-friction and separation bubble height in xOy plane at different spanwise locations

从总体构型来看，分离泡呈现复杂的三维结构特征，其流向尺度明显大于其他两个方向，特别是在法向，其峰值高度仅约为0.12,说明分离泡整体形态以扁平型为主。同时，以摩阻曲线的过零点来确定分离区流向范围，可以看到，分离泡在流向上可以划分为两个区域：次分离泡(Second Bubble, SB)，位于-区间；主分离泡(First Bubble, FB)，位于-区间。可以看到，尽管次分离泡的展向范围略大于主分离泡，但其流向长度和法向高度则明显小于主分离泡，以站位为例，前者分别约为后者的42%和3%。从图7中还可以看到，主分离泡法向高度在展向也存在较大的差异，站位和的峰值高度约为站位的25%和19%，这说明主分离泡沿展向呈现中间高两边低的山峰型结构。

4 非定常特性

图8给出了站位截面内不同时刻的瞬态流向速度云图，图中黑色箭头为速度矢量方向，白色实线为时间平均分离泡的高度。可见，在湍流的强间歇作用下，瞬态分离泡分布极为不规则，在下游再附区内仍存在一定的出现概率。值得特别关注的是，不同时刻下分离泡面积变化非常剧烈。以流向速度<0区域的面积和为表征(图中深蓝色区域)，在/=100时，瞬态分离泡远大于平均分离泡，其峰值高度约为0.3；而在/=300时，其面积又急剧减小，此时分离泡以小尺度多块结构随机分布为主。这表明分离泡具有较强的非定常特性，对应为分离泡的膨胀和收缩过程，这与Priebe等在压缩拐角流动中的研究结论较为类似。

为了定量描述分离泡的膨胀/收缩过程，这里对3个展向站位的瞬态分离泡进行了高频采样，总采样时间跨度为/=442, 约为6个无量纲时间()。图9分别给出了各站位分离泡面积的时间序列，其中黑色虚线表示时间平均值。采样时，还进一步将瞬态分离泡按照图6中的定义区分为主分离泡和次分离泡分别进行统计。

图8 不同时刻展向中截面瞬态流向速度云图Fig.8 Contours of instantaneous streamwise velocity in spanwise mid-plane

图9 分离泡面积时间序列Fig.9 Time series of separation bubble areas

图中蓝色曲线对应为次分离泡，红色曲线为主分离泡，黑色曲线则代表未区分时得到的结果。

从图9中可以清楚看到，相较于主分离泡，次分离泡的瞬时面积非常小，其整体贡献可以忽略不计。这主要是由于尽管次分离泡在流向上具有一定的范围，但是其法向高度非常低，因而导致其占比非常小。而从主分离泡来看，各展向站位处主分离泡的瞬时面积差异较大，但各时间序列与未区分时的统计结果均吻合较好，这也证实了主分离泡在干扰区流动分离现象中占主导作用。因此，在本文后续的分析中主要以主分离泡为研究对象。

图10给出了各站位分离泡面积预乘功率谱的比较情况。图中:PSD为功率谱密度，为频率，为Strouhal 数，定义为=/。以往研究表明，干扰区分离激波的低频振荡运动时间尺度约为10/～100/。从本文计算结果来看，各站位处分离泡的特征频率也基本符合这一规律。可以看到，站位和的预乘功率谱较为接近，主要以双峰结构为主，峰值频率分别约为=0.02和=0.04，而站位处则以单峰特征为主，峰值频率出现在=0.04处，这很可能是由于其分离区尺度相对较小的缘故。结果表明，干扰区分离泡的非定常特性表征为低频膨胀/收缩过程，不同展向站位之间的峰值频率略有差异。

图10 分离泡面积预乘功率谱Fig.10 Pre-multiplied power spectral density for separation bubble areas

为了进一步给出各站位之间膨胀/收缩过程的相关程度，图11还分别给出了站位分离泡面积与站位和分离泡面积的互相关系数。从图11(a)中可以看到，站位与站位存在正的强相关性，两者之间时间延迟约为τ=-11，而图11(b)的结果表明，站位与站位之间则为负相关，延迟时间约为=-3/。显然，站位分离泡的膨胀收缩过程明显滞后于站位和。

图11 分离泡面积互相关系数Fig.11 Cross-correlation coefficient between separation bubble areas

5 条件统计分析

Waindim等采用经验模态分解(Empirical Mode Decomposition, EMD)对分离激波瞬时位置进行了条件分析，研究了分离泡收缩(激波往上游运动)和分离泡膨胀(激波往下游运动)的两者情况下的流场结构特征。与Waindim等的不同之处在于，本节采用EMD方法直接对分离泡面积的时间序列进行分解，条件统计分析低频膨胀/收缩过程中分离泡瞬态几何特征。

首先，通过EMD方法，将原始分离泡脉动信号自适应分解为一系列具有不同特征时间尺度的本征模态函数(Intrinsic Mode Function, IMF)。随后，选取特定的本征模态函数对分离泡的低频膨胀/收缩过程进行了重构，并据此分别对膨胀和收缩两个特定过程中的分离泡瞬态结构进行提取和统计。EMD方法的求解过程如下：

1) 获得分离泡面积脉动原始时间序列()的所有局部极大值()和极小值()。

2) 采用三次样条差值方法，依据局部极值的位置，分别构建原始信号的上包络线()和下包络落线()。

3) 从原始型号()中减去上下包络线的均值()=05(()+())，判断剩余部分()=()-()是否满足本征模态函数的条件。如不满足则将其代替()后再重复步骤1)～步骤3)，直至()满足本征模态函数的条件，并将其记为分离泡面积脉动原始时间序列的第一个IMF分量(IMF)，依次类推，最后将()可分解为若干个本征模态函数和残余项。更为详细的经验模态分解方法介绍可参见文献[26]。

EMD方法具有良好的自适应频率分辨率，无需人为给定滤波阀值。图12分别给出了采用EMD方法对站位处分离泡面积脉动进行分解后得到9个本征模态函数(IMF～IMF)。可以清楚看到，不同的IMF都有着特定的频段，且每个IMF包含的频率成分和频宽都是各不相同的，但整体来看，IMF～IMF对应为从高频到低频的排列，同时脉动幅值呈现增大的趋势。

图12 站位Z2分离泡面积本征模态函数Fig.12 IMFs of separation bubble areas at station Z2

图13分别给出了IMF～IMF对应的预乘功率谱。可见，IMF～IMF的峰值频率由≈4逐渐减小到≈0.01，与之前图10中的研究结果较为符合，这也证实了EMD方法得到的本征模态函数是准确可靠的。从定量比较来看，IMF～IMF的主频基本维持在0.1<<10范围内，这里将其定义为高频模态，而IMF～IMF的能量主要集中在<0.1区间，因而定义为低频模态。图14还给出了低频模态和高频模态重构信号的预乘功率谱，低频重构信号为()=IMF()+IMF()+IMF()+IMF()+IMF()，高频重构信号为()=()-()。如图14所示，低频重构信号的预乘功率谱在<0.1区间内基本与原始信号完全重合，而高频重构信号在>0.3的范围内则与原始信号吻合较好，这说明基于IMF～IMF的低频重构信号能够准确表征分离泡的低频非定常运动过程。

图13 本征模态函数预乘功率谱Fig.13 Pre-multiplied power spectral density for IMFs

图14 分离泡面积脉动低频和高频成分预乘功率谱Fig.14 Pre-multiplied power spectral density for low-and high-frequency of separation bubble area fluctuations

在低频重构出分离泡脉动信号后，采用式(12) 分别确定分离泡的膨胀(Dilation)和收缩(Contraction)运动：

(2)

图15给出了站位处分离泡脉动的低频重构信号以及提取得到的膨胀(红色曲线)和收缩(蓝色曲线)运动。可见，瞬态膨胀和收缩运动的时间尺度差异显著。例如，在/=100，此时膨胀运动对应的时间跨度约14，而随后的收缩运动则只持续了4，这说明此时分离泡的非定常运动以长时膨胀为主。而在/=420，收缩运动的时间尺度明显大于前后时刻的膨胀运动，因而长时收缩运动占主导。总体来看，由于分离泡脉动处于统计定常态，因而膨胀和收缩运动在总时间上的占比基本一致。

图15 站位Z2分离泡面积脉动条件取样分析(红：膨胀；蓝：收缩)Fig.15 Conditional analysis of separation bubble area fluctuations at station Z2(red: dilation; blue: contraction)

图16 分离微团瞬态云图Fig.16 Contour of instantaneous separation micro-clusters

图17 分离微团高度/长度比概率密度分布函数Fig.17 Probability density function of aspect-ratio for separation micro-clusters

图18 分离微团面积和高度的散点图(黑色：压缩；红色：膨胀)Fig.18 Scatter plots of separation micro-clusters area and height(black: contraction; red: dilation)

6 本征正交分解

为了进一步揭示分离泡非定常运动中的典型相干结构，采用POD方法对瞬态流向速度场进行了分析。在笔者前期的研究中，通过对展向平均后的非定常流场进行低阶近似，得到了分离泡非定常演化过程中能量占优的特征模态。但由于POD分析针对的是展向平均场，其结果并不能精确反映出分离泡三维结构差异的影响规律。从本文之前的分析来看，此时分离泡沿展向存在变化剧烈的三维结构(见图6)，因此这里将分别针对3个展向站位(～)的流向-法向剖面内非定常流向速度场开展POD分析。样本总数为4 425，采样区间为-3<(-)<06和0<<04。采样时间取为0.1/，对应的可分辨范围为0.009<<5。流向-法向剖面内非定常流向速度场(,,)的具体分解过程如下：

(3)

图19给出了不同站位归一化后的POD模态能量和累积能量的比较情况。这里第个模态的归一化能量定义为：=∑，其中为第个模态的特征值。从图19(a)中可以看到，模态能量随着模态阶数的增加而急剧降低，>100时，模态能量下降了约两个数量级，其能量衰减率基本符合律，这与Mustafa等的研究结果是一致的。不同站位能量分布曲线也基本重合，低阶模态能量略有差异。以第1阶模态为例(能量贡献最大，下文简称主能量模态)，站位处主能量模型的能量占比约为20.8%，略高于站位和站位处的15.6%和16.4%。从图19(b)中累积能量分布还可以看到，站位～的前10阶模态能量占比分别约为42.1%、44.9%和43.9%，这说明相较于高阶模态，低阶模态的贡献明显占优。

图19 模态能量分布Fig.19 Mode energy distributions

为了考察POD模态的非定常特性，图20分别给出了不同站位模态系数的预乘功率谱。可见，各站位的结果基本类似，随着模态阶数增加，占优频率从低频区往高频区快速移动。对于<10的低阶模态，其峰值频率主要集中在0.01<<0.1范围内，这与图10中分离泡脉动的特征频率较为接近；而对于>100的高阶模态，其能谱则以>1的高频特征为主，低频区内则没有明显的脉动能量。图21还给出了各模态时间系数()与分离泡面积()的相关系数。这里相关系数定义为

图20 模态时间系数预乘谱云图Fig.20 Contours of pre-multiplied spectra for time coefficient of POD modes

(4)

由式(4)可以看到，各站位主能量模态与分离泡面积脉动强相关，站位处为正相关=0.76，站位和处为负相关=-0.89和=-0.91。对于<10时，尽管相关系数急剧衰减，但两者之间仍存在着弱相关。然而，高阶模态与分离泡面积脉动的相关性系数则基本维持在零附近。研究结果表明，分离泡的低频膨胀/收缩运动与低阶模态(<10)密切相关，特别是主能量模态。

图21 模态时间系数与分离泡面积相关系数Fig.21 Correlation between time coefficient and separation bubble area

图22给出模态1和模态100的空间分布情况。图中黑色虚线为时间平均声速线，粉色实线为分离泡外边界。从定性比较来看，尽管分离泡在展向存在复杂的三维结构，其流向和法向尺度差异显著，但各站位主能量模态的空间结构基本一致，分离泡展向三维结构的影响可以忽略不计。如图22所示，站位～主能量模态均对应为大尺度含能结构，模态能量集中在声速线下方和分离泡上方的狭长区域，而在上游湍流边界层和分离泡内并没有出现明显的大尺度含能结构。研究表明，主能量模态与分离泡上方剪切层密切相关，特别是其脚部区域，这与前人在平板入射激波中的研究结论是较为类似。造成这一现象的主要原因很可能是分离泡上方的剪切层仍以准二维结构特征为主，见图22(a)～图22(c)中声速线的位置。高阶模态的空间则与主能量模态完全不同，对应为含能较低的小尺度正负交替结构，表征了分离泡的高频脉动过程。

图22 POD模态空间分布Fig.22 Spatial distribution of POD modes

本文研究结果也进一步支持了激波湍流边界层干扰低频振荡现象的下游机制。对于上游机制，大量学者认为低频振荡现象来源于上游湍流边界层速度型、压力脉动、拟序结构等因素，而下游机制则认为其诱因来源于下游干扰区内，如分离泡、剪切层等下游因素。从计算结果来看，低阶模态与分离泡的低频膨胀/收缩运动存在较强关联(见图21)，同时低阶模态含能结构集中出现在分离泡上方的剪切层，而非上游湍流边界层内。为了定量考察低阶模态和分离泡低频膨胀和收缩的内在关联，这里采用前10个POD低阶模态对各站位的非定常脉动速度场进行了重构，具体过程为