0 引 言
可压缩壁湍流在航空航天工程中广泛存在。可压缩壁湍流的转捩、非定常分离、激波及其与边界层的相互作用等复杂流动已成为超声速、高超声速先进飞行器及其动力装置设计的关键问题[1]。与之相关的风洞，特别是超声速风洞，实验研究受到较大的限制，不同实验结果之间也有一定的差异。随着高性能计算机及高精度数值格式的发展，直接数值模拟已成为探索复杂流动机理，发现新的流动现象及流动规律的重要工具。直接数值模拟不仅可以给出各个流动时刻清晰的流动图画，而且可以计算包括脉动运动在内的所有湍流瞬时运动量在三维空间中的时空演化过程。因此，可以通过建立可靠的直接数值模拟数据库来很好地弥补上述实验研究的不足。
转捩主要可分为两种类型[2]：自然转捩和旁路转捩。前者为小扰动诱发的转捩，经过感受性阶段、不稳定波的线性增长阶段和非线性增长阶段，层流破碎转捩到湍流。后者为较大的局部干扰导致的转捩，边界层内的不稳定波直接进入非线性增长阶段，层流转捩为湍流。目前，对第二类转捩知之甚少，主要依靠经验处理，且风洞特别是超声速风洞实验研究比较有限[3]。
在直接数值模拟过程中，为了促进转捩的发生，通常需要添加一些人工扰动。常采用的扰动形式有：在壁面上增加随机吹吸扰动或具有随机相位的周期性吹吸扰动、在入口处增加周期性尾迹或格栅湍流形式的扰动。这两种扰动经常用来诱导旁路转捩的发生。对于自然转捩，通常根据线性稳定性理论或线性Navier-Stokes(N-S)方程组来确定所需的小扰动，例如在入口处增加二维或三维Tollmien-Schlichting(T-S)波形式的扰动等。在外流中，通常采用壁面上的吹吸扰动诱导旁路转捩的发生。例如，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, 13, 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方程组：
γ=1.4
Pr=0.72
Ts=124K
对式(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取固定值，它们和β的取值在后文详细给出。
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 |
下面给出固定相位周期性吹吸扰动作用下，在湍流充分发展区所得的平均及高阶统计量与相关参考文献、理论分析和实验数据的对比结果。在以下的计算过程中，取φ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)中的两个参考值分别为符合不可压缩平板边界层壁面律的U+VD=y+和不可压缩平板边界层对数律的U+VD=1/0.41 ln(y+)+5.1，其中U+VD为经过Van Driest变换后的速度。可以看出，本文直接数值模拟结果与理论分析符合甚好。由上可知：尽管本文采用固定相位周期性吹吸扰动，诱导转捩发生的扰动不再具有随机性，但所得湍流场与具有随机性的扰动诱导的湍流场具有相同的平均特征。
 |
| 图 6 边界层转捩及湍流充分发展区热力学变量和 速度平均值的统计分析Fig. 6 Density contours and the statistics of thermodynamical variables and velocity |
图 7给出了本文平板表面平均摩擦系数的分布 与层流和湍流边界层理论预测值及相关参考文献计算结果的对比。可以发现，在层流区本文的计算结果与层流边界层理论预测值符合甚好，但在湍流区本文的计算结果大于湍流边界层理论的预测值。这些差异与他人的直接数值模拟结果是一样的。与相关文献结果相比，本文的计算结果更接近理论预测值，原因可能与近壁区网格的分辨率相关。
 |
| 图 7 摩擦系数的分布Fig. 7 Averaged surface friction distribution |
摩擦系数曲线从层流值开始突然显著地增大，然后慢慢减小，逐渐靠近湍流边界层理论的预测值。转捩过程即发生在摩擦系数突然增大的区域。可以看出，三个直接数值模拟结果的平均转捩位置有明显的差异，本文的转捩位置最为靠后。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等的计算结果符合甚好，但这些变量脉动均方根值明显小于文献的计算结果，但二者的分布趋势是一致的。由此可以推测：与转捩的平均位置相似，不同的转捩诱导方式可以影响到热力学变量脉动均方根的大小。
图 9是雷诺正应力和剪切应力的统计结果和相关的实验数据和参考文献计算结果的对比。图 9(a)的参考值为文献[17]不可压缩平板边界层湍流DNS的结果。可以发现，流向雷诺正应力明显大于不可压缩流直接数值模拟的结果，法向和展向分量的差异相对较小。图 9(b)的参考值为文献[17]中雷诺剪切应力直接数值模拟结果和公认的不可压缩平板边界层湍流实验数据[17]。图 9(c)的参考值分别为公认的不可压缩平板边界层湍流实验数据[17]和文献[4]直接数值模拟的计算结果。图 9(b)、图 9(c)中的实验数据均没有考虑密度的影响。可以看出，本文的直接数值模拟的计算结果稍小于两个参考文献DNS的结果，但三个DNS结果的分布趋势是一致的，说明直接数值模拟采用的数值方法可以影响到雷诺应力的计算。与不可压缩平板边界层湍流实验值相比，考虑密度影响的雷诺剪切应力的峰值大于实验值，反之，则小于实验值。相对而言，考虑密度变化后的剪切应力分布与实验吻合较好。
 |
| 图 9 雷诺应力的分布Fig. 9 Distributions of Reynolds stress in normal direction |
采用高精度有限差分格式和基于Euler方程组的特征边界条件，通过直接数值模拟研究了周期性吹吸扰动的相位和频率对超声速平板边界层旁路转捩的影响，并建立了由固定相位周期性吹吸扰动诱导的、Ma∞=2.25、Re∞=2.5×107/m超声速平板边界层湍流的直接数值模拟数据库。数值计算结果表明：
(1) 对于周期性吹吸扰动，相位的随机性不是导致旁路转捩发生的必要条件；采用固定相位，尽管扰动不再具有随机性，但边界层内旁路转捩仍然会发生。
(2) 采用固定相位周期性吹吸扰动，随着吹吸频率的增大，转捩的平均位置向下游移动；但频率超过某一临界值时，边界层不再发生转捩。
(3) 无论采用随机相位或固定相位，湍流的平均特征是一致的；二阶统计量具有相同的变化趋势。相对于平均参数而言，平均转捩位置和二阶统计量对转捩的诱导因素、数值方法和边界条件较为敏感。
| [1] | Roman D, Allen J B, Liebeck R H. Aerodynamic design challenges of the blended-wing-body subsonic transport[R]. AIAA 2000-4335. |
| [2] | William J G. Innovative control effect ors (configuration 101) dynamic wind tunnel test report[R]. AFRL-VA-WP-TR-1998-3043, 1998 |
| [3] | Wong M D, Flores J. Application of overflow-MLP to the design of the 1303 UCAV[R]. AIAA 2006-2987. |
| [4] | Rizzi A, Tomac M, Nangia R. Engineering methods for SACCON configuration[R]. AIAA 2010-4398. |
| [5] | Liu Xiaojing, Wu Jianghao, Zhang Shuguang. Aerodynamic design and optimization of the blended wingbody aircraft for 250 passengers[J]. Acta Aerodynamica Sinica, 2011, 29(1): 78-84.(in Chinese) 刘晓静, 吴江浩, 张曙光. 250座级翼身融合布局客机气动设计与优化[J]. 空气动力学学报, 2011, 29(1): 78-84. |
| [6] | Gong Junfeng, Zhu Xiaoping, Tao Yujin. Split drag rudder hinge moment predict for flying wing aircraft[J]. Acta Aerodynamica Sinica, 2010, 28(4): 472-477. (in Chinese) 龚军锋, 祝小平, 陶于金. 飞翼飞机开裂式阻力方向舵铰链力矩的预测[J]. 空气动力学学报, 2010, 28(4): 472-477. |
| [7] | Zuo Linxuan, Wang Jinjun. Experimental study of the effect of AMT on aerodynamic performance of tail less flying wing aircraft[J]. Acta Aerodynamica Sinica, 2010, 28(2): 132-137. (in Chinese) 左林玄, 王晋军. 全动翼尖对无尾飞翼布局飞机气动特性影响的实验研究[J]. 空气动力学学报, 2010, 28(2): 132-137. |
| [8] | Kong Yinan, He Kaifeng, Wang Lixin, et al. Simulation of flyingwing aircraft maneuver with leading edge jet control[J]. Acta Aerodynamica Sinica, 2013, 31(2): 186-191. (in Chinese) 孔轶男, 何开锋, 王立新, 等. 增加喷流控制的飞翼飞机机动仿真[J]. 空气动力学学报, 2013, 31(2): 186-191. |
| [9] | Chu J, Luckring J M. Experimental surface pressure data obtained on 65 deg delta wing across Reynolds number and Mach number ranges[R]. NASA 1996-4645. |
| [10] | Willy Pritz. Numerical simulation of the peculiar subsonic flow-field about the VFE-2 delta wing with rounded leading edge[R]. AIAA 2008-393. |
| [11] | Simone Crippa, Arthur Rizzi. Steady, subsonic CFD analysis of the VFE-2 configuration and comparison to wind tunnel data[R]. AIAA 2008-397. |
| [12] | Polhamus E C. Application of the leading-edge-suction analogy of vortex lift to the drag due to lift of sharp-edge delta wings[R]. NASA TN-D-4739, 1968. |
| [13] | Hunt J C R, Wray A A, Moin P. Eddies, streams, and convergence zones in turbulent flows[R]. N89-24555/9/XAD, 1989. |
| [14] | Marcel Lesieur. Large-eddy simulations of turbulence[M]. Cambridge UniversityPress, 2005. |
| [15] | Kalkhoran I M, Smart M K. Aspects of shock wave-induced vortex breakdown[J]. Progress in Aerospace Sciences, 2000, 36(1): 63-95. |
| [16] | Délery J M. Aspects of vortex breakdown[J]. Progress in Aerospace Sciences, 1994, 30(1): 1-59. |
| [17] | Schiavetta L A, Boelens O J, Crippa S. Shock effects on delta wing vortex breakdown[J]. Journal of Aircraft, 2009, 46(3): 903-914. |
| [18] | Visbal M R. Computational and physical aspectsof vortex breakdown on delta wing[R]. AIAA 95-0585. |
| [19] | Robinson B A, Barnett R M, Agrawal S. Simple numerical criterion for vortex breakdown[J]. AIAA Journal, 1994, 32(1): 116-122. |