激波/湍流边界层干扰物面剪切应力统计特性

童福林，周桂宇，周浩，张培红,*，李新亮2,

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

激波与湍流边界层的相互作用是现代高速飞行器气动设计中不可或缺的基础问题，至今仍未被充分理解。激波/湍流干扰区内会出现流动分离/再附、局部强压力脉动和热流峰值，这会对飞行器气动性能、防热层以及结构疲劳等方面产生显著影响。因此，进一步深入研究干扰区内复杂流动现象有助于加深对该问题的理解认识，为工程应用提供重要的理论参考依据。Dolling[1]和Gaitonde[2]从热流预测、激波的非定常运动特性以及流动控制等方面对该问题进行了详细的综述。

依据激波的产生方式不同，可以将激波与边界层干扰划分为[3]：压缩拐角、入射激波干扰、双锥、后掠压缩拐角、单楔、双楔以及内流道问题。压缩拐角和入射激波干扰是其中两类最具代表性的流动构型。自20世纪40年代以来，国内外大量学者对这两类构型进行了系统的风洞试验和数值模拟研究。

在风洞试验方面，Settles和Fitzpatrick[4]研究了不同激波强度下压缩拐角干扰区内物面压力、摩阻及平均速度的演化规律。Ardonceau[5]和Smits等[6]分析了激波干扰对湍流脉动的增强机制。结果表明，激波干扰湍流剪切应力的影响要明显强于雷诺应力的其他分量，干扰区下游湍流边界层内质量通量脉动强度显著增强。分离激波的低频振荡现象及其物理机制一直以来都是风洞试验的研究热点。Andreopoulos和Muck[7]发现上游湍流边界层的猝发现象与低频振荡运动密切相关。但随后，Erengil和Dolling[8]的研究结果表明，上游压力脉动才是激波非定常运动的主要来源。目前，对于该问题的实验研究尚未达成共识。大量研究表明，激波运动机制与上游边界层速度型剖面[9]、大尺度高低速条带结构[10]等因素密切相关。

在数值模拟方面，随着计算速度和数值格式的飞速发展，直接数值模拟(Direct Numerical Simulation，DNS)方法已逐渐成为激波/湍流边界层干扰复杂流动机理方面的重要研究手段。相较于风洞试验，DNS可以直接给出干扰区内试验难以测量的流场信息。Adams[11]首次采用DNS方法研究了压缩性对干扰区下游湍流结构的影响机制，发现边界层猝发频率与激波运动频率较为接近。Ringuette等[12]数值研究了雷诺数对干扰区物面压力脉动频谱特性、分离区长度以及湍流脉动的影响规律。Priebe等[13-14]着重探讨了激波运动的低频振荡现象。低通滤波后的DNS瞬时流场表明，分离激波的非定常运动与分离泡的膨胀/收缩存在较强关联。此外，李新亮等[15]研究了激波低频振荡机制以及干扰区湍动能的输运机制。结果证实了激波的低频振荡与上游边界层拟序结构无关。近年来，童福林等[16-18]开展了大量的激波/湍流边界层干扰直接数值模拟研究，探讨了激波强度、壁面温度、马赫数等因素对干扰区内复杂流动现象的影响规律。

总体来看，国内外在激波与湍流边界层相互作用问题上取得了长足的进步，并在一些复杂问题流动机理方面达成共识。但在干扰区内物面剪切应力统计特性方面，相关试验及DNS的研究报道较为少见。Murthy和Rose[19]对马赫数Ma=2.9下的入射激波湍流边界层干扰问题进行了实验研究，获得了干扰区内物面剪切的平均量和脉动量。进一步深入开展干扰区内物面剪切应力统计特性的演化规律研究，有助于为改进现有湍流模型和亚格子模型提供理论支撑。

本文采用直接数值模拟方法对入射激波/平板湍流边界层相互作用问题进行数值研究。着重探讨分离激波低频振荡运动对物面剪切应力功率谱密度的影响机制，研究分离泡内流向及展向剪切应力分量的概率密度分布规律。采用本征正交分解方法，分析比较了流向剪切应力脉动与物面压力脉动的差异。为了便于比较和验证结果，计算参数的选取与Bookey等[20]的实验结果和Priebe等[13]的DNS结果相近。

1 计算参数

直接数值模拟的控制方程为三维可压缩无量纲Navier-Stokes方程组：

(1)

式中：Q为守恒变量；F、G和H为3个方向上的无黏通量；Fv、Gv和Hv为3个方向对应的黏性通量，具体表达式参见文献[21]。方程的无量纲化采用无穷远处来流参数以及单位特征长度。计算时，为了抑制激波间断区的数值振荡同时保证对湍流边界层内不同尺度流动结构的高分辨率，采用Martin等[22]优化构造的WENO_SYMBO_LMT格式以及Steger-Warming流通量分裂方法计算无黏项。同时，采用八阶中心差分格式对黏性项进行离散，时间推进采用三阶Runge-Kutta方法计算。需要特别指出的是，本文DNS采用的高精度差分求解器OpenCFD-SC软件已在多个激波/湍流边界层干扰问题[15-18]中得到了成功的验证和确认，可以保证DNS结果的准确和可靠。

图1 计算模型示意图Fig.1 Illustration of computation model

如图1所示，计算模型为入射激波与平板湍流边界层的相互作用问题。气流方向为从左往右，来流马赫数为2.9，基于单位长度的来流雷诺数为5 581.4 mm-1，来流静温为108.1 K，壁面温度为307 K。计算域流向跨度为-363 mm

网格点数为3 200×200×140(流向×法向×展向)，计算网格采用代数解析方法生成，流向网格在激波与湍流边界层的干扰区内均匀分布(如图1所示)，法向网格在近壁区采用了双曲正切函数的加密处理，展向网格均匀分布。这里以x=-60 mm处壁面量为参考量，干扰区内网格尺度分别为Δx+=4.5、 Δy+=0.5、 Δz+=5.0,与Priebe等[13]的DNS结果较为接近。如无特别说明，本文中上游湍流边界层的统计变量均取自x=-60 mm处(位于充分发展湍流边界层内)。表1分别给出了上游湍流边界层的马赫数Ma、边界层厚度δ、位移厚度δ*、动量厚度θ和物面摩阻系数Cf。

表1 上游湍流边界层参数Table 1 Parameters of incoming turbulent boundary layer

2 结果验证

本节通过与以往数值模拟结果和风洞试验数据的对比分析，进一步验证计算结果的准确性，其中包括上游湍流边界层的统计特性、干扰区内平均压力和摩阻分布以及物面压力脉动的功率谱等。

图2 湍流边界层统计特性Fig.2 Statistical characteristics of turbulent boundary layer

图3分别给出了激波/湍流边界层干扰区内物面压力pw/p∞和摩阻系数沿流向的分布情况。为了便于比较说明，这里将压力和摩阻分布的流向坐标均进行了平移和无量纲处理，其中x*为流向坐标平移后的值，xsep为平均分离点流向坐标。计算得到的压力和摩阻分布与Priebe等[13]的数值结果均基本重合。与Bookey等[20]的试验数据比较来看，干扰区下游的物面压力值要明显高于试验值。在直接数值模拟时，展向取为周期性边界条件，而实际风洞试验时展向为真实固壁，因此洞壁干扰[13]是造成该差异的主要因素之一。

另一方面，图4给出了计算得到的上游湍流边界层和分离平均起始点物面压力脉动的预乘谱(Pre-Multiplied Power Spectral Density,fPSD)，其中f为频率，PSD为功率谱密度。可见，在无干扰区内，压力脉动的无量纲峰值频率出现在1.0U∞/δ附近,U∞为来流速度。由于干扰区内分离激波的低频振荡运动，分离点物面压力脉动的低频能量在(0.004～0.01)U∞/δ的范围内急剧增强。计算值与以往激波湍流边界层干扰直接数值模拟得到的低频峰值频率范围(0.002～0.006)U∞/δ[13]较为接近，这也证实了本文DNS计算采用的数值方法和网格分辨率能够准确捕捉到干扰区内的分离激波低频振荡现象。

图3 物面平均压力及摩阻系数分布Fig.3 Distribution of wall average pressure and skin friction coefficient

图4 物面压力脉动预乘谱Fig.4 Pre-multiplied power spectral density of wall pressure fluctuations

3 流场结构

图5分别给出了入射激波与平板湍流边界层干扰区内的无量纲瞬态密度梯度场和时间平均密度场。图中红色和蓝色曲线为Ma=1和u=0的瞬态等值线，这里u为流向无量纲速度。可以看到，在入射激波和分离激波的相互作用下，干扰区内存在着强逆压梯度，边界层内出现了大范围的流动分离，同时在下游边界层外缘还存在较强的压缩波系。此外，从两者的定性比较来看，分离泡内存在强烈的间歇性和非定常特征。

研究结果表明，物面剪切流向分量要比其展向分量大了约一个量级，这表明干扰区内物面剪切以流向剪切为主。图6分别给出了物面流向剪切应力τx的瞬态和时均分布云图。为了便于下文比较说明，沿流向选取了5个典型特征位置，其中E1～E3分别位于上游无干扰边界层内、分离泡和下游再附区内，S和R分别位于平均分离和再附点。从整体分布规律来看，上游无干扰边界层内物面剪切以条带结构特征为主，这与边界层近壁区的高低速条带紧密相关。随后，该条带结构在分离区内被破坏并消失，分离区内流向剪切表征为强间歇特性。在干扰区下游，流向剪切的量值呈逐渐增大趋势，但其展向分布规律与上游存在明显差异。从时均结果也可以看到，干扰区下游的流向摩阻沿展向表现为强烈的非均匀性。

图5 瞬态密度梯度和时均密度流场Fig.5 Instantaneous density gradient and mean density flow fields

图6 物面流向剪切应力瞬态和时均分布云图Fig.6 Instantaneous and mean contours of streamwise component of wall shear stress

4 物面剪切应力统计特性

4.1 预乘谱分析

分析物面剪切应力信号的功率谱密度有助于理解干扰区内各剪切分量脉动特征的演化规律。此外，分离激波的低频振荡运动对物面剪切脉动的影响同样值得关注。

图8给出了干扰区内各特征位置物面剪切应力预乘谱分布的比较情况。可以看到，对于流向分量，相较于上游E1处，分离区内和再附区下游边界层内峰值频率略有升高，低频区脉动能量也有一定的增强，但从整体分布规律来看，脉动能量仍然以高频特征为主。干扰区内展向分量的演化规律与流向分量基本类似，只是在干扰区下游S和E3处的峰值频率略有降低，但高频脉动仍占主导。研究结果也进一步表明，分离激波的低频振荡运动不会对干扰区内物面剪切应力各分量的脉动产生实质影响。

图7 湍流边界层内物面剪切应力信号预乘谱Fig.7 Pre-multiplied power spectral density of wall shear stress in incoming turbulent boundary layer

图8 干扰区内物面剪切应力信号预乘谱Fig.8 Pre-multiplied power spectral density of wall shear stress in interaction region

如图4所示，分离激波的大尺度低频振荡运动对物面压力脉动影响显著，尤其是对平均分离点附近的低频能量。而从图8的结果来看，激波的非定常运动对流向和展向分量的影响则要弱得多。尽管在低频部分，有一个数量级的升高，但脉动仍以高频能量为主。从定性分析来看，认为造成这种差异的原因很可能有以下两方面。首先，分离激波非定常运动的影响作用主要体现在流场中零阶分量，如压力、质量通量[28]等，而物面剪切应力是流向或展向速度的法向梯度，为流场参数的一阶分量。可以观察到，尽管干扰区内剪切脉动的低频能量有一定的升高，但其增长速率要明显低于压力脉动。另一方面，从图5中还可以看到，分离激波在边界层外层逐步弱化为弱压缩波系，因而其非定常运动对近壁区内流动参数的影响程度会急剧减弱。从图4和图8的比较来看，边界层内激波强度的弱化对物面剪切的影响程度要远强于物面压力脉动，具体更为详细的定量作用机制有待下一步深入研究。

4.2 概率密度分布

为了研究干扰区内物面剪切应力演化的统计特性，图9给出了E1～E3处流向剪切分量及其脉动的概率密度函数(PDF)。在本节中剪切应力分量采用当地的时空平均值τx,av和τz,av进行归一化处理，分别表示为

(2)

从图9可以看到，在上游干扰区E1和下游再附区E3，流向分量的PDF曲线都近似地呈现对数正态分布。尽管流动在这两个特征位置均以附着流为主，但仍存在一定的概率出现负的流向剪切。以往零压力梯度平板的试验结果表明[29-30]，边界层内存在较小概率的回流现象(Backflow)，上述研究成果也进一步证实了该流动现象的存在。另外，在分离泡内E2处，从统计意义来看，流动以负剪切为主要特征，但其函数分布在正剪切范围内仍存在较高的可能性，这主要由于分离泡的强间歇性的缘故。对于流向剪切的脉动量，由于采用当地均方根进行了归一化处理，可以清楚看到，干扰区上游和下游的函数分布近似重合，且与Carlos等[26]的不可压平板数值结果吻合较好，两者只是在脉动量变化剧烈的区域存在较明显的差别，这表明在这个区域内，局部应力脉动的统计特征是相似的。此外，在分离泡内的脉动量变化范围要明显大于前两个区域。

图10给出了干扰区内展向剪切及其脉动值的概率密度分布曲线。函数与流向分量的分布规律则完全不同。如图所示，干扰区内的展向剪切均近似呈现正态分布。偏斜因子的计算表明，3个特征位置处的偏斜因子绝对值小于0.1，这说明分布函数以对称特征为主。但在变化较剧烈的区域，分离泡内的发生概率要明显高于其他两个区域。与此同时，干扰区内归一化后的展向脉动值均与Carlos等[26]的统计结果吻合。综上所述，在激波湍流边界层干扰区内，研究表明，相较于上游充分发展湍流边界层，分离泡内的流向剪切统计特性变化剧烈，而展向剪切统计特性的变化则可忽略不计。

图9 干扰区内流向剪切应力概率密度函数Fig.9 Probability density functions of streamwise component of wall shear stress in interaction region

为了更好地考察流向剪切与展向剪切之间的相互关系，定义两者的夹角ψτ(t)为

ψτ(t)=arctan(τz(t)/τx(t))

(3)

图10 干扰区内展向剪切应力概率密度函数Fig.10 Probability density functions of spanwise component of wall shear stress in interaction region

图11给出了干扰区内夹角的概率密度函数分布。总体来看，函数呈对称分布，在变化剧烈的区域要明显高于高斯正态分布。对于上游湍流边界层，夹角主要集中出现在[-45°,45°]范围内，这与Jeon等[31]的研究结果是一致的。在分离泡内E2处，可以看到，小夹角事件概率略有降低。值得注意的是，E2处夹角大于25°的可能性则要明显高于其他两个区域，这主要是由于分离泡内流向剪切急剧降低，而展向剪切变化较小，两者较为接近，导致大夹角事件的发生概率也相应升高。在干扰区下游再附区E3处，流向剪切逐渐恢复到初始值，因而大角度事件的发生概率又随之降低。

图11 物面剪切应力夹角概率密度函数Fig.11 Probability density functions of angle between wall shear stress components

图12 物面剪切应力矢量夹角和幅值联合概率密度函数Fig.12 Joint probability density functions of angle and magnitude of wall shear stress vector

图13 联合概率密度随物面应力矢量幅值的变化Fig.13 Variation of joint PDFs with magnitude of wall shear stress vector

4.3 本征正交分解

为了进一步深入分析物面剪切脉动的能量结构，采用本征正交分解(POD)方法探究了非定常物面剪切脉动场的典型相干结构。通过POD方法可以对复杂高维度的流场进行低阶近似，提取出非定常演化历程中能量占优的特征模态。假设非定常物面剪切场为T(x,z,t)，POD分析可以确定一族正交基函数φj(x,z)，j=1,2,…，具体分解过程如下[32]：

T(x,z,t)=T(x,z,t)

(4)

POD分析针对400个流向/展向平面内瞬态物面流向剪切脉动场进行操作。依据模态特征值对模态能量Ej进行排序，归一化的模态能量定义为

(5)

图14给出了物面流向剪切场POD模态能量的分布情况。可以看到，高能量模态主要集中在前20个模态，约占总能量的50%。

图15分别给出了物面流向剪切脉动场的能量占优主模态的空间结构。图中云图为归一化后的正交基向量，虚线表征了分离激波非定常运动的流向范围，通过计算压力脉动信号的间歇因子[16]得到。为了便于比较差异，图中还给出了物面压力脉动场对应的能量模态结构，这里SS1、SS2、SS3和SS40分别代表流向剪切脉动场的第1、2、3和40个POD模态，而P1、P2、P3和P40分别为压力脉动场对应的POD模态。如图15(a) 和图15(b)所示，剪切场的第1个模态空间结构与压力场较为相似，沿展向近似呈现二维分布，约占总能量的16%。从分布规律来看，能量结构主要集中在平均分离点S及间歇区内，这也说明了分离激波沿流向的非定常低频振荡运动主导了物面压力和流向剪切脉动场。可以清楚看到，两者不同之处在于，流向剪切场在平均再附点下游还有一定的能量结构。

图14 归一化POD模态能量分布及累积能量Fig.14 Distribution of normalized energy of POD modes and cumulative energy

图15 物面压力和流向剪切应力POD模态空间分布Fig.15 Spatial distribution of POD modes of wall pressure and streamwise shear stress

如图15(d)和图15(f)所示，物面剪切场第2和第3模态的空间结构则与第1模态完全不同，此时，空间结构主要以再附点下游沿展向正负交替大尺度结构为主，平均分离点附近结构强度则要弱得多。以往的动态模态分解(DMD)研究表明[33]，再附点下游的摩阻分布与Görtler-like涡结构密切相关。本文的POD结果也进一步证实了该研究结果。此外，尽管压力场的第2和第3模态结构也集中在再附点附近，但其空间分布沿展向仍呈二维结构，在流向表征为正负交替结构，这很可能与下游分离泡沿流向的膨胀/收缩运动有关。图15(g)和图15(h)给出了第40个模态的空间分布规律，该模态为低能量模态，仅占总量的0.5%。可见，剪切脉动场与压力脉动场的空间分布差别较小，均以小尺度结构特征为主。

5 结 论

本文采用直接数值模拟方法研究了来流马赫数2.9、12°激波角的入射激波与平板湍流边界层相互作用问题，详细地分析了干扰区内物面剪切应力场的典型统计特征，如预乘谱、概率密度分布和相干结构等，得到以下结论：

1) 数值模拟准确捕捉到了分离激波的非定常运动。与物面压力脉动不同的是，分离激波低频振荡运动对物面剪切应力预乘谱没有实质影响，干扰区内流向及展向切应力的脉动能量仍以高频特征为主。

2) 干扰区内分离泡对物面剪切应力各分量统计特性的影响机制差异明显。流向剪切概率密度函数变化剧烈，分离泡内不再满足对数正态分布规律，而展向剪切概率密度函数的变化则较小，近似于正态分布。

3) 主能量模态的空间结构表明，分离点附近物面剪切脉动与分离激波的低频振荡运动密切相关，而下游再附区内则由大尺度Görtler-like流向涡结构占主导。

致 谢

感谢国家超级计算天津中心、国家超级计算长沙中心、中国科学院网络中心超级计算中心以及山西吕梁超算中心提供计算机时。