超声速平板边界层旁路转捩直接数值模拟

朱海涛，单 鹏

(北京航空航天大学能源与动力工程学院，北京 100191)

超声速平板边界层旁路转捩直接数值模拟

朱海涛，单 鹏*

(北京航空航天大学能源与动力工程学院，北京 100191)

采用具有七阶精度的有限差分格式和基于Euler方程组的特征边界条件，通过直接数值模拟研究了由固定相位周期性吹吸扰动诱导的、Ma∞=2.25、Re∞=2.5×107/m(635 000/inch)空间发展平板边界层的旁路转捩过程，并建立了相应的直接数值模拟数据库。数值结果表明:扰动相位的随机性并不是转捩发生的必要条件，固定相位周期性吹吸扰动即可诱发边界层的转捩;当采用固定相位时，随着吹吸频率的增加，转捩位置向下游移动;当扰动频率足够大时，周期性吹吸无法诱导出边界层的转捩。随后，在所得数据库中提取了速度及热力学变量的平均值，与已有的数据库和理论分析吻合良好;速度脉动均方根与公认的不可压缩平板边界层实验基本吻合;热力学变量脉动均方根和雷诺应力与已有数据库的统计结果及实验数据具有相同的分布特征。对比发现，不同转捩的诱导因素对平均转捩位置和热力学变量均方根、雷诺应力等高阶统计量的影响较大。

平板;超声速;直接数值模拟;湍流边界层;雷诺应力;转捩

0 引 言

可压缩壁湍流在航空航天工程中广泛存在。可压缩壁湍流的转捩、非定常分离、激波及其与边界层的相互作用等复杂流动已成为超声速、高超声速先进飞行器及其动力装置设计的关键问题［1］。与之相关的风洞，特别是超声速风洞，实验研究受到较大的限制，不同实验结果之间也有一定的差异。随着高性能计算机及高精度数值格式的发展，直接数值模拟已成为探索复杂流动机理，发现新的流动现象及流动规律的重要工具。直接数值模拟不仅可以给出各个流动时刻清晰的流动图画，而且可以计算包括脉动运动在内的所有湍流瞬时运动量在三维空间中的时空演化过程。因此，可以通过建立可靠的直接数值模拟数据库来很好地弥补上述实验研究的不足。

转捩主要可分为两种类型［2］:自然转捩和旁路转捩。前者为小扰动诱发的转捩，经过感受性阶段、不稳定波的线性增长阶段和非线性增长阶段，层流破碎转捩到湍流。后者为较大的局部干扰导致的转捩，边界层内的不稳定波直接进入非线性增长阶段，层流转捩为湍流。目前，对第二类转捩知之甚少，主要依靠经验处理，且风洞特别是超声速风洞实验研究比较有限［3］。

在直接数值模拟过程中，为了促进转捩的发生，通常需要添加一些人工扰动。常采用的扰动形式有:在壁面上增加随机吹吸扰动或具有随机相位的周期性吹吸扰动、在入口处增加周期性尾迹或格栅湍流形式的扰动。这两种扰动经常用来诱导旁路转捩的发生。对于自然转捩，通常根据线性稳定性理论或线性Navier-Stokes(N-S)方程组来确定所需的小扰动，例如在入口处增加二维或三维Tollmien-Schlichting(TS)波形式的扰动等。在外流中，通常采用壁面上的吹吸扰动诱导旁路转捩的发生。例如，Rai［4］等采用高精度有限差分格式实现了Ma∞=2.25的平板边界层流动的直接数值模拟，发现速度脉动均方根的分布与不可压缩平板边界层湍流实验结果很吻合，从而验证了Morkovin假设的正确性;Rai采用的扰动即为在平板壁面上增加具有随机相位的周期性吹吸扰动。采用相同的扰动，Gatski等［5］和Pirozzoli等［6］分别对超声速平板边界层湍流的高阶统计量和湍动能的生产和耗散过程进行了细致的研究。在入口处增加尾迹或格栅湍流形式的扰动在内流中比较常用。例如Rodi和Wissink等［7-8］采用直接数值模拟的方法研究了在周期性尾迹作用下低压涡轮叶栅中分离流、分离泡、转捩等复杂的流动现象，并研究了这些流动现象对涡轮叶片表面热传导效应的影响，数值结果和实验结果吻合良好。Zaki等［9］采用基于交错网格的有限体积法，通过求解不可压缩N-S方程组，对自由湍流作用下和周期性尾迹作用下V103叶栅流场的非定常分离及旁路转捩过程进行了直接数值模拟，数值结果反应的平均流场、分离及转捩等特征与实验结果是一致的。Egorov等［10］采用根据线性N-S方程组得到小扰动波，研究了超声速平板边界的感受性问题。国内陈健［11］等完成了周期性吹吸扰动作用下光滑圆管流动转捩过程的直接数值模拟。傅德熏等［1，12-14］采用具有随机相位的周期性吹吸扰动完成马赫数Ma∞=2.25和Ma∞=6的平板边界层湍流的直接数值模拟，并研究了Ma∞=2.25平板湍流边界层中的声效应，以及高马赫数下湍动能的生成和耗散机制。李新亮等［15］采用随机吹吸扰动对Ma∞=6的钝锥绕流进行了直接数值模拟，辅以线性稳定性理论分析了钝锥表面局部区域转捩位置相对靠后的原因。

尽管已经得到了若干直接数值模拟的数据库，但是仍然需要新的、采用不同数值方法、边界条件和转捩诱导因素的直接数值模拟结果对已有的数据库进行验证和检验，以期对所建立的DNS数据库和转捩的发生有较为全面的认识。本文采用具有七阶精度的有限差分格式和基于Euler方程组的特征边界条件［16］，以固定相位周期性吹吸为转捩的诱导因素，对Ma∞=2.25、Re∞=2.5×107/m的空间发展平板边界层湍流进行了直接数值模拟，并研究了吹吸扰动的相位和频率对转捩发生的影响。

1 数值方法及问题描述

考虑如下三维非定常N-S方程组:

补充完全气体状态方程、温度与总能量、热传导系数和Prandtl数之间的关系式及粘性系数满足的Sutherland公式，

上述公式中用到的常数值如下:

至此，完成方程组(1)的封闭过程。在下文的计算过程中，采用了在一个方向上进行加密的非均匀网格。此时，需要在方程组(1)中的对应方向上无粘和粘性空间导数项乘以该方向上的网格变换系数。此外，粘性应力张量和热通量中该方向上的空间导数项也要做类似的变换，此处不再详述。本文计算过程中参考长度L0=0.0254 m(1 inch)，基于L0的雷诺数Re∞= 635 000;自由来流的马赫数Ma∞=2.25，温度T∞= 169.44 K。上述方程组根据自由来流的密度ρ"、速度V∞、温度T∞及上述参考长度L0进行无量纲化，并以方程组(1)进行直接数值模拟。

对式(1)中的无粘项，首先采用Steger-Warming方法进行流通矢量分裂，然后对正负流通矢量均采用守恒型的七阶迎风差分格式计算其空间导数。粘性项采用具有八阶精度的中心差分格式离散。最后通过具有TVD性质的三步三阶Runge-Kutta公式进行时间积分。

本文的直接数值模拟分两步进行，首先进行二维层流计算，然后以所得的稳态二维层流解为初场，通过在壁面增加固定相位的周期性吹吸扰动进行三维边界层转捩的直接数值模拟。二维和三维计算采用的数值方法是相同的。

二维计算区域如图1所示，平板长度取为Lx=0.281 m，入口取在平板前缘前方0.0762 m的地方，上边界与平板的距离为Ly=0.127 m。入口边界条件为自由来流边界;平板壁面为等温壁面条件，Tw=322 K;从入口到平板前缘采用对称边界条件(图中虚线);平板上方边界采用基于Euler方程组的特征边界条件;出口边界当流向马赫数小于1时采用特征边界条件，反之，由插值获取出口边界上的流场信息;二维计算的初场取为自由来流。在进行一定时间的流场推进后即可获得稳态的二维层流解。

图1 二维计算区域Fig.1 Two-dimensional calculation area

在三维计算中，平板长度取为Lx=0.14224 m，法向高度取为 Ly=0.00508 m，展向长度为 Lz= 0.00889 m。壁面条件与二维计算一致;展向采用周期性边界条件;壁面上方边界采用如下基于Euler方程组的特征边界条件其中下标b、e、"分别表示边界处、在边界处进行外插和自由来流对应物理量的值;n为边界处指向计算域内部的单位法向量;c为当地声速;γ=1.4。出口边界条件与二维计算一致:若流向马赫数小于1.0，则由式(3)所示的特征边界条件确定，其中，pb为沿法向最后一个流向马赫数小于1.0的网格点处的压力值;若流向马赫数大于等于1.0则直接进行外插确定出口边界上的流场信息。

入口边界和初场的设置如下所述:获得二维稳态结果后，取二维平板Lx=0.1778m处的数值解(从平板前缘算起，此处已远离平板前缘的激波)沿展向复制作为三维计算的入口边界条件，展向速度为0，然后把入口边界处的物理量沿流向复制，作为三维直接数值模拟的初场。入口边界条件在三维瞬态计算过程中保持不变。

计算网格在流向和展向为均匀网格，分别取Nx=1681和Nz=320个网格点。法向网格沿壁面按双曲正切函数进行加密，网格数为Ny=71。本文三维计算的网格数及无量纲化的空间步长和相关文献的对比如表1所示。相较于其他文献，本文在流向和展向上的空间步长较大，但由于本文采用高精度的计算格式，在一定程度上可以降低对空间步长的要求［2］。

为了加快层流到湍流的转捩，计算过程中在平板流向Lx=0.0127～0.0254 m的区域添加了如式(4)所示的固定相位周期性吹吸扰动:vbs为平板壁面法向的吹吸速度;xa、xb为扰动区域沿流向上的起始和终止坐标;φl、φm分别为扰动的空间和时间相位，β为扰动的基频。在本文的计算中，A=0.04，xa=0.5，xb=1.0，lmax=10，mmax=5。φl、φm取固定值，它们和β的取值在后文详细给出。

表1 本文及相关文献的网格参数Table 1 Non-dimensional mesh parameters of present and references

2 直接数值模拟结果

2.1 扰动相位和频率的影响

对于周期性的吹吸扰动，通常时间和空间相位φl、φm均取［0，1］内的随机数。在本文的直接数值模拟过程中，采用固定相位的吹吸扰动，即诱导转捩发生的扰动不再具有随机性。为了较为快速地检验不同的时间和空间相位组合能否导致转捩的发生，本小节的计算中采用了较为稀疏的网格:流向上取Nx= 1121，展向上取Nz=160，法向上仍取Ny=71，吹吸频率β=2.5。图2给出了i)φl=0，φm=0;ii)φl=1/4，φm=0;iii)φl=0，φm=1/2三种相位组合情形下，平板表面摩擦系数的时均分布图。在ii)和iii)情形下，时间和空间扰动中幅值最大的扰动波的相位分别相差π/2、π。可以看出三种相位组合均导致了转捩的发生。图3给出了z=0.0055截面上转捩过程中和充分发展湍流场中展向涡量的分布云图。可以看出经过边界层平均流和周期性扰动的相互作用，诱导出了不稳定的剪切层;在x=2.2附近，该剪切层失稳导致了层流的破碎;破碎的边界层不断向下游发展，逐渐形成充分发展的湍流场，如图3所示。因此，对于周期性吹吸扰动，相位的随机性并不是转捩发生的必要条件，湍流场所具有的随机性是其内在的本质特征，与扰动的随机性并无直接的联系。

图2 扰动相位对平板表面时均摩擦系数的影响Fig.2 Time averaged surface friction with different time and space disturbance phases

图3 不同时刻z=0.0055截面上展向涡量的分布云图Fig.3 Distribution of spanwise vorticity on the plane of (x，y)at z=0.0055，t=4.2 and t=8.0

图4给出了情形i)下，不同吹吸频率对平板表面平均摩擦系数的影响。可以发现，当β=3.5时，摩擦系数不断地减小，与层流对应的摩擦系数值很相近，表明边界层内没有转捩的发生。β=2.5、1.5时，摩擦系数在平均转捩位置处突然变大，随后保持在一个显著地大于层流摩擦系数的水平。显然β=1.5对应的转捩位置比较靠前。因此，当吹吸频率增大时，转捩位置向平板的下游移动，但当频率继续增大时，固定相位的周期性吹吸扰动无法诱导出转捩的发生，这与高慧等人［13］的结论是一致的。

图4 不同吹吸频率下平板表面平均摩擦系数的分布Fig.4 Averaged friction with different frequencies of blow and suction disturbances

2.2 统计结果分析

下面给出固定相位周期性吹吸扰动作用下，在湍流充分发展区所得的平均及高阶统计量与相关参考文献、理论分析和实验数据的对比结果。在以下的计算过程中，取φl=0，φm=0，β=2.5。

由相关文献的计算可知，在流向坐标0.12192 m处，边界层已成为充分发展的湍流边界层。图5给出了y+=5.1323，流向坐标为0.12192 m处密度、速度和温度脉动的二阶自相关量的统计平均值沿展向的分布。可以看出，各个二阶自相关量的统计平均值均趋向于0。因此，本文计算域展向长度是足够的。

图5 二阶自相关量统计平均值沿展向的分布Fig.5 Correlation tensor distribution along spanwise

图6(a)为t=8.0时，流场z=0.0055处x-y截面上密度分布云图;上方为层流边界层失稳及破碎时的密度分布;下方为充分发展的湍流区中密度的分布云图。可见，尽管所添加的扰动不再具有随机性，但边界层中仍然诱导出了非线性增长的不稳定波，使层流逐渐失稳、破碎，最终发生转捩，发展为充分发展的湍流边界层。在充分发展的湍流区中可以发现低密度流体的上抛和高密度流体的下卷，并由此形成锐利的湍流边界。图6(b)和图6(c)分别是在流向坐标为0.12192m处热力学变量的统计平均值沿法向的分布和经过Van Driest变换［2］后的平均速度剖面。图6 (b)中的参考值为Gatski等［5］人的直接数值模拟结果。可以发现，二者基本吻合。图6(c)中的两个参考值分别为符合不可压缩平板边界层壁面律的= y+和不可压缩平板边界层对数律的=1/0.41 ln(y+)+5.1，其中为经过Van Driest变换后的速度。可以看出，本文直接数值模拟结果与理论分析符合甚好。由上可知:尽管本文采用固定相位周期性吹吸扰动，诱导转捩发生的扰动不再具有随机性，但所得湍流场与具有随机性的扰动诱导的湍流场具有相同的平均特征。

图6 边界层转捩及湍流充分发展区热力学变量和速度平均值的统计分析Fig.6 Density contours and the statistics of thermodynamical variables and velocity

图7 摩擦系数的分布Fig.7 Averaged surface friction distribution

图7给出了本文平板表面平均摩擦系数的分布与层流和湍流边界层理论预测值及相关参考文献计算结果的对比。可以发现，在层流区本文的计算结果与层流边界层理论预测值符合甚好，但在湍流区本文的计算结果大于湍流边界层理论的预测值。这些差异与他人的直接数值模拟结果是一样的。与相关文献结果相比，本文的计算结果更接近理论预测值，原因可能与近壁区网格的分辨率相关。

摩擦系数曲线从层流值开始突然显著地增大，然后慢慢减小，逐渐靠近湍流边界层理论的预测值。转捩过程即发生在摩擦系数突然增大的区域。可以看出，三个直接数值模拟结果的平均转捩位置有明显的差异，本文的转捩位置最为靠后。Rai［4］和Gao［13］的结果均是基于具有随机相位周期性扰动的，由此可知:不同的转捩诱导方式对转捩平均位置具有一定的影响。

图8是在流向坐标为0.12192 m处湍流充分发展区速度和热力学变量脉动均方根(Root Mean Square，RMS)沿法向的分布曲线。图8(a)中的参考数据为不可压缩平板边界层的实验结果［17］。可以发现，尽管自由流的Ma∞=2.25，但在平板边界层内，速度脉动均方根沿法向的分布与不可压缩平板边界层的分布特征十分接近。说明当来流马赫数为2.25时，流场的可压缩性对速度脉动的影响也是比较小的。此时，可压缩效应对湍流脉动动力学特征的直接影响较弱，其影响主要反映在对流场平均流动特征的改变。因此，通过可以把平均密度变化影响考虑进去的Van Driest变换，流场的平均速度分布与不可压缩平板边界层的壁面律和对数律速度剖面是相吻合的，Morkovin假设也是成立的。图8(b)中的参考值为Gatski等［5］人直接数值模拟结果。可以发现，尽管本文热力学变量的统计平均值与Gatski等的计算结果符合甚好，但这些变量脉动均方根值明显小于文献的计算结果，但二者的分布趋势是一致的。由此可以推测:与转捩的平均位置相似，不同的转捩诱导方式可以影响到热力学变量脉动均方根的大小。

图8 速度及热力学变量脉动均方根的分布Fig.8 Distributions of the RMS of the fluctuations of velocity and thermodynamics variables

图9是雷诺正应力和剪切应力的统计结果和相关的实验数据和参考文献计算结果的对比。图9(a)的参考值为文献［17］不可压缩平板边界层湍流DNS的结果。可以发现，流向雷诺正应力明显大于不可压缩流直接数值模拟的结果，法向和展向分量的差异相对较小。图9(b)的参考值为文献［17］中雷诺剪切应

图9 雷诺应力的分布Fig.9 Distributions of Reynolds stress in normal direction

力直接数值模拟结果和公认的不可压缩平板边界层湍流实验数据［17］。图9(c)的参考值分别为公认的不可压缩平板边界层湍流实验数据［17］和文献［4］直接数值模拟的计算结果。图9(b)、图9(c)中的实验数据均没有考虑密度的影响。可以看出，本文的直接数值模拟的计算结果稍小于两个参考文献DNS的结果，但三个DNS结果的分布趋势是一致的，说明直接数值模拟采用的数值方法可以影响到雷诺应力的计算。与不可压缩平板边界层湍流实验值相比，考虑密度影响的雷诺剪切应力的峰值大于实验值，反之，则小于实验值。相对而言，考虑密度变化后的剪切应力分布与实验吻合较好。

3 结 论

采用高精度有限差分格式和基于Euler方程组的特征边界条件，通过直接数值模拟研究了周期性吹吸扰动的相位和频率对超声速平板边界层旁路转捩的影响，并建立了由固定相位周期性吹吸扰动诱导的、Ma∞=2.25、Re∞=2.5×107/m超声速平板边界层湍流的直接数值模拟数据库。数值计算结果表明:

(1)对于周期性吹吸扰动，相位的随机性不是导致旁路转捩发生的必要条件;采用固定相位，尽管扰动不再具有随机性，但边界层内旁路转捩仍然会发生。

(2)采用固定相位周期性吹吸扰动，随着吹吸频率的增大，转捩的平均位置向下游移动;但频率超过某一临界值时，边界层不再发生转捩。

(3)无论采用随机相位或固定相位，湍流的平均特征是一致的;二阶统计量具有相同的变化趋势。相对于平均参数而言，平均转捩位置和二阶统计量对转捩的诱导因素、数值方法和边界条件较为敏感。

［1］ Li Xinliang，Fu Dexun，Ma Yanwen.Direct numerical simulation of compressible turbulence flows［J］.Acta Mechanica Sinica，2010，26(6):795-806.

［2］ Fu Dexun，Ma Yanwen，Li Xinliang，et al.Direct numerical simulation of compressible turbulent flow［M］.Beijing:Science Press，2010.(in Chinese)

傅德熏，马延文，李新亮，等.可压缩湍流直接数值模拟［M］.北京:科学出版社，2010.

［3］ Egorov I V，Fedorov A V，Soudakov V G.Direct numerical simulation of unstable disturbances in supersonic boundary layer［R］.AIAA 2004-588.

［4］ Rai M M，Gatski T B，Erlebacher G.Direct numerical simulation of spatially evolving compressible turbulent boundary layers［R］.AIAA 95-0583.

［5］ Gatski T B，Erlebacher G.Numerical simulation of a spatially evolving supersonic turbulent boundary layer［R］.NASA TM 2002-211934.

［6］ Pirozzoli S，Grasso F.Direct numerical simulation and analysis of a spatially evolving supersonic turbulent boundary layer at M=2.25［J］.Physics of Fluids，2004，16(3):530-545.

［7］ Rodi W.DNS and LES of some engineering flows［J］.Fluids Dynamics Research，2006，38(2-3):145-173.

［8］ Wissink J G.DNS of separation，low Reynolds number flow in a turbine cascade with incoming wakes［J］.Engineering turbulence modeling and experiments，2002，(5):731-740.

［9］ Zaki T A，Wissink J G，Durbin P A，et al.Direct computations ofboundary layers distorted by migrating wakes in a linear compressor cascade［J］.Flow Turbulence Combust，2009，83:307-322.

［10］Egorov I V，Fedorov A V，Nechaev A V.Receptivity of supersonic boundary layer on a blunt plate to acoustic disturbances［R］.AIAA 2004-249.

［11］ Chen Jian，Cui Guixiang，Xu Chunxiao，et al.Direct numerical simulation of spatial evolutional transition in pipe flow［J］.Chinese Journal of Theoretical and Applied Mechanics，2003，35(1):6-13.(in Chinese)

陈健，崔桂香，许春晓，等.空间发展圆管转捩的直接数值模拟［J］.力学学报，2003，35(1):6-13.

［12］Li Xinliang，Fu Dexun，Ma Yanwen，et al.Acoustic calculation for supersonic turbulent boundary layer flow［J］.Chinese Physics Letters，2009，26(9):094701.

［13］Gao Hui，Fu Dexun，Ma Yanwen，et al.Direct numerical simulation of supersonic turbulent boundary layer flow［J］.Chinese Physics Letters，2005，22(7):1709-1702.

［14］Li Xinliang，Fu Dexun，Ma Yanwen.Direct numerical simulation of spatially evolving supersonic turbulent boundary layer at Ma=6［J］.Chinese Physics Letters，2006，23(6):1519-1522.

［15］Li Xinliang，Fu Dexun，Ma Yanwen.Direct numerical simulation of hypersonic boundary layer transition over a blunt cone with a small angle of attack［J］.Physics of Fluids，2010，22(2):1-18.

［16］Hirsch C.Numerical computation of internal and external flows，Vol 2:computational methods for inviscid and viscous flows［M］.Chichester:John Wiley＆Son，1990.

［17］Spalart P R.Direct numerical simulation of a turbulent boundary layer up to Reθ=1410［J］.Journal of Fluid Mechanics，1988，187: 61-98.

Direct numerical simulation of by-pass transition of the supersonic flat plate turbulent boundary layer

Zhu Haitao，Shan Peng*
(School of Jet Propulsion，Beihang University，Beijing 100191，China)

Using the characteristic boundary conditions and high-order finite difference schemes，direct numerical simulation(DNS)was performed to investigate the effect of phase and frequency of the periodic blow and suction disturbance which was used to trigger the by-pass transition of the flat plate boundary layer with free stream Ma∞=2.25，Re∞=2.5×107/m(635 000/inch).The numerical results show that the randomness of the disturbance phase is not the necessary condition to trigger the transition，and with the deterministic disturbance phase，the average location of transition moves backward when the disturbance frequency increases.However，if the frequency is bigger than a critical value，the transition in the boundary layer will not happen.A DNS database was established under the condition of deterministic disturbance phase.The Van Direst velocity agrees well with the wall and logarithm law of incompressible flat plate boundary layer，and the turbulence intensities coincide with the incompressible flat plate experimental datum.The averaged thermodynamic variables are very close to that of reference results.These verify the present DNS results.Then，the statistics of thermodynamic variables fluctuation and Reynolds stress are extracted from the DNS database.Compared with the relative references'DNS results，the disturbance，boundary conditions and numerical methods have observably effect on the second order statistics and the averaged transition position.

flat;supersonic;DNS;turbulent boundary layer;Reynolds stress;transition

0258-1825(2015)03-0345-08

O357.4;V211.3

A

10.7638/kqdlxxb-2013.0083

2013-08-07;

2013-11-07

朱海涛(1982-)，男，山西省襄汾县人，博士研究生，研究方向:计算流体力学，大规模科学和工程计算.E-mail:zht83@126.com

单鹏*(1955-)，男，北京人，教授，研究方向:航空发动机，叶轮机械和计算空气动力学.E-mail:pshan@buaa.edu.cn

朱海涛，单鹏.超声速平板边界层旁路转捩直接数值模拟［J］.空气动力学学报，2015，33(3):345-352.

10.7638/kqdlxxb-2013.0083 Zhu H T，Shan P.Direct numerical simulation of by-pass transition of the supersonic flat plate turbulent boundary layer［J］.Acta Aerodynamica Sinica，2015，33(3):345-352.