赵晓龙，张君安，董 皓，刘 波

(西安工业大学 机电工程学院，西安 710021)

近年来，随着航空航天事业的发展，超声速飞行器的战略价值受到了各科技大国的高度重视，成为各国近期的研究热点[1].超声速与低速最大的不同之一便是激波.在超声速流场中，气流的黏性效应强，激波与物面附近的边界层流动会发生相互干扰，形成近壁区的复杂流场，导致干扰区内激波的位置、形状、激波前后的流动参量以及边界层特性均发生变化，甚至有可能支配整个流场[2-5].由于超声速的地面实验存在难度大，成本高，周期长和飞行条件不易控制等困难，因此，工程中常常依靠数值方法来弥补这些不足.

研究人员针对超声速计算问题做了一些研究.文献[6]研究了超声速大攻角旋转飞行器气动特性的数值模拟方法，以Woodward有限基本解方法为基础，计算了旋转对弹体和弹翼涡流生成、发展的影响.文献[7]对跨声速和超声速流中激波/边界层干扰进行了数值模拟，针对湍流模型中的Bradshaw 常数进行了修正，并对跨声速和超声速流中激波/边界层干扰进行了数值模拟研究，计算得到的壁面压力分布、分离区长度和速度剖面均与实验值吻合较好.文献[8]使用雷诺平均方法研究了平板超声速流中的激波和边界层问题.文献[9]研究了超声速二维底部压力的简化解析计算方法，应用质量和能量守恒原理，推导出了一个超声速二维底部压力的简化解析算法.文献[10]针对超声速气流绕局部透气凹角流动问题进行了数值模拟，用具有压力分裂形式的简化纳维-斯托克斯方程为控制方程，数值模拟了超声速来流条件下的激波-边界层干扰被动控制，数值计算时采用了多重扫描法对控制方程差分离散.文献[11]在迎风型格式和矢通量分裂技术的基础之上，对捕捉激波方法进行一种新的尝试，该方法抑制非物理波动，将其转换成守恒形式，得到了基本上无振荡的激波捕捉格式.

综上所述，目前对于超声速问题计算主要从数学模型、控制方程及差分格式上进行研究.由于超声速问题本身的复杂性，本文采用简化的纳维-斯托克斯方程作为控制方程难以真实地模拟出气体流动的实际状态.本文采用直接数值模拟方法，计算连续方程、含有黏度项的纳维-斯托克斯方程和能量方程，不对应力分量和黏度进行简化，同时联立完全气体状态等其他基本方程，使用高精度的差分格式进行离散，保证计算数据的准确性，同时研究分析了不同雷诺数对尖前缘平板超声速绕流状态的影响.

1 计算模型及控制方程

本文的物理模型如图1所示，一个带有尖角的二维平板放置在马赫数Ma大于1的超声速气流中.其中板长为L，初始来流的压强p、密度ρ、速度V和温度T均为标准值，超声速气流流过平板上表面.

图1 物理模型

建立包含了黏性效应和应力分量的连续方程、动量方程和能量方程的向量形式为

(1)

其中

式中：U，E，F分别为所构建向量形式的三个分量；t为时间；x,y为位置坐标；ρ为密度；u为x方向的速度；v为y方向的速度；V为速度；Et为单位体积流体动能和内能之和；e为内能；p为流场内压力；τxx,τxy,τyx,τyy为不同方向的应力；qx,qy为热流分量；k为热传动率；λ为第二黏性系数；μ为动力黏度；T为温度.

完全气体假定的状态方程为

由此可见，“301条款”是一项美国制定的单方面发动调查的国内法规则，美国总统凭此规定可采取单方行动使美国从他国对其不公平贸易行为中解脱出来。

p=ρRT

(2)

其中R为理想气体常数.

进一步假定为定压比热的完全气体，则有

e=cvT

(3)

(4)

式中：|V|为速度的模；u为x方向的速度；v为y方向的速度；cv为定容热容.

假定气体是定压比热的完全气体，则动力黏度为

(5)

其中μ0，T0分别为标准黏度和标准温度.

(6)

式中：Pr为普朗特常数；cp为定压比热.

由式(1)～(6)联立求解，采用2阶精度的数值差分，可以求得U.其差分格式为

(7)

时间步长的计算公式为

Δt=min[K(Δt0)i,j]

(8)

2 数值计算分析

为使计算方便，本文选取L=0.01 mm进行直接计算，则该模型的计算坐标如图2所示，x为平板表面坐标位置；y为边界层位置.计算网格既能光滑无振荡地捕捉激波，又具有相当高精度以保证能精确描述边界层内不稳定因素的微小变化.本文选用的计算参数见表1.

图2 计算坐标

计算区域进行初始化，同时分别计算该区域内的边界条件、动力黏度、热传导系数和应力；使用时间推进的计算方法，单独求解每次推进的时间步长；设定收敛准则.本文数值计算的收敛准则为：在每一时间步后，若所有网格点上的密度较前一时间步长上的密度值的变化不超过1×10-6kg·m-3，则认为收敛，若不收敛，则返回第二步再次计算推进时间步长.比较入口质量流量和出口质量流量两者偏差，解耦输出所需参数.

表1 计算参数Tab.1 Calculation parameters

根据控制方程，计算马赫数等于2，4，6，8时坐标内气体速度、压力、温度和密度的分布状态，可得气体状态如图3所示.

由图3可知，在计算坐标区域内，气体的速度、压力、温度和密度产生了明显的非线性变化，存在明显的速度、压力、温度和密度分解线，这个分界线就是通过数值计算捕捉到的激波.在激波的分界线的上方气体的各个状态基本不发生变化，约等于流入平板时的气体状态；在激波分界线下方，随着越靠近平板的位置，来流的气体状态变化越强烈，导致这种变化的原因就是因为气体黏度的存在，气流的速度越高，黏度对气体状态的影响越大.

图3 不同马赫数时计算坐标内的气体状态

图4为马赫数分别等于2，4，6和8时平板表面的压力分布，由图4可知，当超声速的气流流过平板表面时，平板表面的压力先激增后急剧下降到某一值，然后趋于稳定；且马赫数越大，沿平板表面的压力分布越大.当马赫数等于8时，平板表面的最大压力与标准压力的比值为5.65倍；马赫数等于2时，平板表面的最大压力与标准压力的比值为2.05倍.随着马赫数的增大，在沿平板表面位置上，压力急剧下降后，又略微升高，然后缓慢下降.当马赫数等于8时，在平板表面位置等于0.13×10-5m时，平板表面的最大压力与标准压力比值为3.1倍；当平板表面位置等于0.2×10-5m时，平板表面的最大压力与标准压力比值增大至3.3倍.

图5为不同马赫数时边界层的流体状态.随着边界层坐标数值增大，流场内的气体与标准压力和标准温度的比值先增大后减小；马赫数持续增大，直至等于外界来流气体的马赫数.当马赫数等于2时，计算区域的边界层为0.12×10-5m，流场内最大压力与标准压力比值为2倍，最高温度与标准温度比值为1.2倍，在边界层坐标等于0.8×10-5m时马赫数达到最大；马赫数等于4时，计算区域的边界层为0.8×10-6m，最大压力与标准压力比值为2.9倍，最高温度与标准温度比值为1.5倍，在边界层坐标等于4.3×10-6m时马赫数达到最大；马赫数等于6时，计算区域的边界层为6.83×10-6m，最大压力与标准压力比值为4.1倍，最高温度与标准温度比值为2.4倍，在边界层坐标等于3.3×10-6m时马赫数达到最大；马赫数等于8时，计算区域的边界层为5.7×10-6m，最大压力与标准压力比值为5.6倍，最高温度与标准温度比值为3.3倍，在边界层坐标等于2.9×10-6m时马赫数达到最大.

图4 平板表面压力分布

图5 不同马赫数时边界层的流体状态

3 结 论

1) 应用直接数值模拟求解含有黏度项的纳维-斯托克斯方程，反应了所求模型的实际流体状态，直接数值模拟所求模型的精度优于其他方法.

2) 气体以超声速的速度流入平板上时，在平板上产生激波.随着马赫数的增大，产生激波的范围越来越小，靠近平板位置的速度小于流入平板的速度，但压力整体大于远离平板位置的压力.

3) 流场内的气体的最大压力和最大温度出现在平板表面的边界层中，随着马赫数增大，出现的位置越来越靠近平板表面，随着边界层坐标数值增大，马赫数持续增大，直至等于外界来流气体的马赫数.