0 引言

湍流的多样性使得描述一个湍流流动入口是非常困难的。对于湍流的数值模拟，通常需要在入口给定复杂的湍流脉动生成方法来模拟真实的湍流场。常用的湍流入口边界生成方法有自然转捩方法、合成湍流方法和流场参数回收方法等等^[1-3]。

自然转捩方法最为简单^[4]，但计算量大，计算域包含了层流区、转捩区和湍流区。流场参数回收方法是采用计算本身来产生入口所需的湍流脉动信息，它其中又包括利用时间发展湍流场方法(流向周期性方法)^[5-7]和回收/调节方法^[8-9]，该方法适合的研究对象是不关心转捩区的充分发展湍流。在工程上，经常使用在壁面上引入粗糙度、凸起或壁面变形等使流动和几何参数发生改变的方法，对转捩产生影响^[10-11]。

由于压缩性的影响，高马赫数下边界层流动更为稳定，因此在高超声速下生成湍流场就更为困难。一些常用的生成湍流场的方法在高马赫数条件下不能达到预期效果。比如回收/调节方法，该方法把下游某处速度场、热力学参量和各种湍流统计量按照一定的比例缩放后引入入口，但这样引入的前提是Morkovin假说成立^[9]，其实该方法在M＞5时就不是很有效了。在实际数值计算中也发现，在低马赫数情况下(M=3)使用该方法能产生湍流场，而在高马赫数情况下(M=6)就不能激发出湍流场。

本文主要针对高马赫数流动，对自然转捩、波纹壁面结构和利用时间发展流场的参数回收这几种方法进行了数值计算，并对结果进行了分析和评价。

1 控制方程和数值方法

考虑一般曲线坐标系(ξ, η, ζ)，在此坐标系下三维Navier-Stokes方程形式为：

$ \frac{{\partial \hat U}}{{\partial t}} + \frac{{\partial \hat E}}{{\partial \xi }} + \frac{{\partial \hat F}}{{\partial \eta }} + \frac{{\partial \hat G}}{{\partial \zeta }} = \frac{{\partial {{\hat E}_v}}}{{\partial \xi }} + \frac{{\partial {{\hat F}_v}}}{{\partial \eta }} + \frac{{\partial {{\hat G}_v}}}{{\partial \zeta }} $

(1)

采用完全气体假设，方程的具体表达式参见文献[12]。方程(1)的数值求解采用有限差分方法，通量采用Sterge-Warming矢通量分裂形式，通量离散波纹壁计算时使用五阶WENO格式^[13]，其他算例使用五阶迎风格式；黏性项采用六阶中心差分格式；时间方向采用具有TVD性质的三步三阶Runge-Kutta方法^[14]推进求解。

边界处理：壁面采用无滑移条件和等温壁条件；出口采用外推方法；展向采用周期边界条件。

2 几种可压缩湍流入口条件生成方法

气体参数选取26km高空处值，来流马赫数为6，壁温为600K。计算时高空来流速度、密度和温度值作为相应无量纲参考量，长度用L=1m作为无量纲参考量。

2.1 自然转捩方法

在计算域入口存在边界层问题时，将计算域入口取在上游层流段，根据Orr-Sommerfeld方程得到不稳定特征模态，叠加至层流剖面上，激励扰动经过线性、非线性增长转捩至成为完全发展的湍流。

在高马赫数情况下，最不稳定的是第二模态二维扰动^[15]。为促使扰动快速进入非线性增长阶段，在计算域入口处加入一个第二模态二维扰动和一对第一模态三维扰动，扰动的形式为：

$ u' = A\hat u\left( y \right){{\rm{e}}^{{\rm{i}}\left( {\alpha x + \beta z - \omega t} \right)}} $

(2)

其中：A为扰动幅值，β为展向波数，ω为扰动的频率，α为在给定参数下Orr-Sommerfeld方程的特征值，$ \hat u\left( y \right) $为相应的特征函数。具体参数见表 1，α=α_r+iα_i是复数，α_i＜0表示扰动是增长的。

图 1中黑色实线给出了自然转捩方法得到的物面摩擦系数分布，虚线是层流和湍流壁面摩擦系数的理论值^[16]。在x=0.55m左右曲线开始抬升，转捩开始，经过一段距离后达到湍流状态。图 2给出了计算到统计定常后速度梯度张量二阶不变量Q2=5000的三维等值面，可以看到，计算域入口附近流动是以二维结构为主，随着向下游的推进，开始有三维结构出现，Λ涡形成，在流向x＞0.9m到达完全湍流区。图 3给出了流向为x=0.9m处的平均速度剖面沿壁面法向的分布，实线为平均速度，三角形为Van Driest变换后的平均速度，虚线分别为壁面律(u⁺=y⁺)和对数律(u⁺=2.5ln(y⁺)+5.8)，平均速度剖面与对数率和壁面率符合得很好。在计算域入口处添加上述扰动，数值模拟获得了平板边界层从层流失稳转捩到湍流的整个发展过程。

2.2 壁面上引入波纹结构

在结构表面上引入粗糙单元或改变壁面局部结构，都可使局部流场状态发生改变，促使流向涡的产生，使流动剪切作用增强，起到促进转捩的效果。

在流向位置为x₁=0.45m至x₂=0.5m之间布置波纹壁面，波纹壁的分布函数为：

$ y = \left\{ \begin{array}{l} 0, \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;x < {x_1}\\ h\sin \left( {{k_x}\left( {x - {x_1}} \right)} \right)\sin \left( {{k_z}z} \right), \;\;\;\;\;\;{x_1} \le x \le {x_2}\\ 0, \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;x > {x_2} \end{array} \right. $

(3)

其中：高度h=2.45mm(约为0.7倍当地边界层厚度)，k_x=5×2π/(x₂－x₁)，k_z=2×2π/l_z，l_z为展向长度。

图 4给出了Q2=5000的三维等值面，在下游由于波纹壁面的干扰作用，使得流动剪切作用增强，流场变得不稳定，层流快速转捩到湍流，从图 4中可以看出x≈0.6m左右开始流动结构就比较饱满了。图 1中实心方形表示的是波纹壁方法得到的物面摩擦系数曲线，由于波纹壁的存在，能引发壁面摩擦系数曲线迅速抬升，在经过波纹壁面之后，曲线会缓慢下降，经过一段距离后达到完全湍流状态。图 5给出了x=0.68m处的平均速度剖面，可以看到平均速度剖面同对数率和壁面率符合得较好。

2.3 利用时间发展湍流场进行参数回收方法

采用时间发展湍流场作为空间发展流场的入口条件，预先模拟的时间发展湍流场流向采用周期性边界条件，在进行回收时只需用预先模拟的湍流场的某一时刻的流动信息即可，预存的数据量小。

首先，计算一个随时间发展的湍流场。与2.1节计算自然转捩时引入扰动的方法类似，初始时刻在整个计算域中引入一些扰动，扰动形式为：

$ u' = A\hat u\left( y \right){{\rm{e}}^{{\rm{i}}\left( {\alpha x + \beta z} \right)}} $

(4)

具体参数见表 2，此处ω_i代表扰动的增长率。

图 6给出了壁面摩擦系数曲线随时间的分布。可以看到t=0.2s左右壁面摩擦系数曲线开始抬升，也就意味着转捩的开始，到t=0.5s左右达到峰值进入湍流区，到达湍流区后系数曲线的变化趋于平缓且有所下降。图 7为0.8s时速度梯度张量二阶不变量Q2=5000的三维等值面，此时已是充分发展湍流。取0.8s时刻的瞬时流场数据作为空间发展湍流流场的入口条件，假定此瞬时湍流流场数据为U_T(x, y, z)，则空间模式DNS入口条件可表示为：

$ {U_{{\rm{inlet}}}}\left( {t, y, z} \right) = {U_T}\left( {x, y, z} \right) $

(5)

这里通过时空变换关系，x=ct，将某一时刻随时间发展的湍流场整个计算域的流场信息转换为空间发展的湍流流场入口处的时间序列。黄章峰等^[5]认为，壁湍流是具有大尺度相干结构的，而相干结构对湍流的性质起着决定性的作用，而相干结构具有某些波的特点，其平均传播速度约为0.92，在实际计算中也发现设定为此对流速度能得到较好的效果，因此取传播速度c为0.9。

图 1中的空心圆形为回收方法得到的壁面摩擦系数沿流向的分布，可看到由于入口条件的人为给定，流动需要经过一小段调整，之后就维持在完全湍流状态。图 8为Q2=5000的三维等值面，计算域入口开始就是湍流状态，而且在向下游的推进中能维持住湍流流动。图 9给出了流向为0.55m处平均速度剖面，同对数率和壁面率符合得较好。

3 可压缩湍流入口条件生成方法评价

对第2节中计算出的湍流场进行归纳总结，从生成湍流场品质和能否快速生成湍流场两方面来对这几类入口条件给定方法进行评价。首先，图 10分别给出了y⁺=40时3种方法生成的湍流场速度分布云图，自然转捩方法生成的湍流场最为均匀饱满；引入波纹壁面结构促发转捩生成的湍流场不如自然转捩方法得到的湍流场均匀、饱满；同样，使用随时间发展湍流场进行参数回收方法时得到的湍流场也不如自然转捩方法得到的湍流场饱满。

另从图 1中的壁面摩擦系数曲线对比图，可以直观地看到引入波纹壁面结构方法相比自然转捩方法能更为快速的生成湍流场，不需要经历很长的流向距离(从流动入口位置0.4m到完全湍流处0.7m左右)。而采用时间发展流场进行回收作为入口条件的方法，从计算域入口开始就是湍流场，只要经过一小段的调整(入口0.4m到0.45m左右)就能获得一个比较真实的充分发展湍流，此种方法也是一种快速生成湍流场的方法。

采用时间发展流场进行回收方法可以推广到其它外形和不同计算工况。董明等^[6]证实此方法是对不同马赫数、雷诺数、壁面条件和典型边界层类型均可行的有效的空间模式直接数值模拟入口条件生成方法，在实际数值计算中也发现，使用文献中对平均密度和温度以及脉动量进行相应修正的方法，可导出不同计算条件下的入口条件。在这里用2.3节计算得到的随时间发展的湍流场作为入口，可以计算出来流马赫数7、高空27km、壁面温度为700K时的湍流场。图 11为Q2=5000的三维等值面，可以看到同样从计算域入口开始就是湍流状态。图 12为壁面摩擦系数沿流向的分布，可看到由于入口条件的影响，流动需要经过一段比之前相同计算工况条件下稍长一些的调整距离之后就维持在完全湍流状态。可见，采用随时间发展湍流场进行参数回收方法，可以推广到更高的马赫数。

4 结论

数值模拟研究了自然转捩、波纹壁诱发转捩、利用时间发展边界层湍流这三类促发湍流场生成的入口条件给定方法。研究表明，在高马赫流动下这三种方法都可以实现湍流场的生成。其中自然转捩方法处理起来最为简单，得到的湍流场也最为均匀饱满，但是计算域长，网格量大，不能快速生成湍流场；波纹壁面的引入，能很快诱导湍流生成，是一种能快速生成湍流的方法，但生成的湍流场不如自然转捩方法得到的湍流场均匀饱满；利用时间发展流场进行参数回收这种湍流生成方法，从计算域入口就是湍流场，也是一种能有效生成湍流场的方法，该方法还有一个好处就是可由同一预先模拟的时间发展湍流场导出不同马赫数、雷诺数、壁面条件所需要的湍流入口条件，这对于要计算很多参数下的湍流情况是很有用的。