钝头体中的广义雷诺比拟关系1)

陈星星 陈 皓 范晶晶 温玉芬 张 正 马友林

(中国运载火箭技术研究院，北京 100076)

引言

雷诺比拟描述了流体在壁面上的剪切和传热之间的关联关系，最早由雷诺在1874 年提出[1]，在各类流体传热传质问题中受到广泛关注和研究.在边界层理论研究的早期，Blasius[2]在不可压缩的半无穷长平板绕流中找到了边界层方程的第一个理论解，从中可以计算得到平板壁面的雷诺比拟关系，结果表明摩阻系数和热流系数的比值，即雷诺比拟系数，为常数1/2.

不可压缩平板绕流中的雷诺比拟关系具有十分简单的形式，直观地揭示了边界层中动量和能量输运之间的关系.更进一步的研究表明，即使考虑更加复杂的效应，如可压缩流动[3-4]、尖前缘的稀薄气体效应[5-6]、燃烧[7-8]等，平板绕流中的雷诺比拟系数仍然为常数.

对于工程实践中更常见的湍流边界层流动，雷诺比拟也受到广泛的关注和应用.文献[9-14]利用雷诺比拟关系构建解析理论预示平板和尖锥表面的热流分布，这些理论方法在工程中得到了广泛的应用.Chi 等[15]利用动量方程的湍流模型和雷诺比拟关系构造了能量方程的湍流模型.作为目前湍流流动中少数的理论结果之一，雷诺比拟关系还被用来与各类数值计算、试验测量结果进行对比验证[16-18].

近年来，陈星星等[19]进一步对弯曲壁面上的雷诺比拟关系开展了研究，提出了广义雷诺比拟以描述雷诺比拟系数与壁面外形之间的关系.传统的雷诺比拟关系适用于平板、尖锥绕流等边界层流向压力梯度为零的流动，而对于存在流向压力梯度的流动，则热流和摩阻的比拟关系一般不是常数[20].以二维的圆柱绕流和轴对称的圆球绕流为例，其热流峰值处于驻点位置，沿壁面向下游单调下降，而摩阻则在驻点处为零，沿壁面向下游先增加后减小，二者的比值随壁面位置变化.

陈星星等通过对钝头体壁面附近层流边界层方程的理论分析，推导出在高超声速条件下，壁面摩阻和热流的雷诺比拟系数随与壁面当地倾斜角(壁面当地切线与来流方向夹角的余角)呈正比关系，其比例系数与壁温相关而与马赫数、雷诺数无关.理论分析结果与DSMC 数值仿真计算结果吻合[19].进一步地，陈星星等[21]还将钝头体中的广义雷诺比拟关系推广到的包含化学反应流动的情形.

广义雷诺比拟关系理论与DSMC 数值模拟计算结果吻合较好.但是受DSMC 方法的限制[22-24]，主要适用于稀薄气体或近连续区流动的计算，雷诺数较低，而对于工程中常见的较高雷诺数流动，乃至湍流流动，在当前计算机硬件水平下，采用DSMC 方法很难实现.

另一方面，高超声速钝头体绕流问题是高速气体动力学中的经典问题，在各类航天飞行器中广泛存在，长期以来受到深入研究.文献[25-26]分别研究了钝头体壁面和驻点热流问题，给出了理论预示方法.国内王智慧等[27-29]提出了驻点热流受稀薄气体效应影响的理论判据.在数值求解N-S 方程方面，文献[30-32]研究了不同离散格式、不同数值网格对钝头体热流计算结果的影响.Schwartzentruber 等[33]对比了采用DSMC 方法和数值求解N-S 方法计算得到的圆柱表面摩阻和热流大小.

为了对钝头体中的广义雷诺比拟关系开展进一步研究，本文采用数值求解N-S 方程的方法计算了不同雷诺数、不同外形条件下的钝头体绕流问题，并将壁面上的雷诺比拟系数分布与理论预示结果进行对比，为广义雷诺比拟关系研究提供数据支撑.

1 圆柱壁面的广义雷诺比拟关系

1.1 理论预示

分别定义摩阻系数Cf=和热流系数Ch=，其中τw为壁面剪切力，qw为壁面热流，为来流密度，为来流速度.

对于如图1 所示的高超声速绕钝头体流动，壁面摩阻和热流系数的比值为

其中，Pr为气体普朗特常数，ue为沿流向的边界层外缘速度，分别为壁面上气流无量纲速度和温度的法向梯度.式中摩阻和热流系数的比值Cf/Ch即为雷诺比拟系数.

根据可压缩流动的Bernoulli 方程得到的钝头体边界层外缘速度分布公式[34]

式中，γ 为气体常数，pw为壁面压力，p0为驻点压力.在来流马赫数的条件下，上式可以近似为关于壁面当地倾斜角θ 的线性分布[35]:，θ 的定义见图1.此时可以得到线性的广义雷诺比拟关系式[19]

其中，Cr为一个与壁面温度相关的比例系数，根据前期研究结果,对于二维圆柱外形

对于轴对称圆球外形

式中，Tw为壁面温度，T0为来流气体总温.

图1 钝头体绕流示意图Fig.1 Hypersonic flow around a blunt-nosed body

1.2 N-S 方程数值计算

为了对更一般条件下的广义雷诺比拟问题开展研究，本文采用数值求解N-S 方程的方法模拟圆柱绕流问题，计算壁面的摩阻分布和热流分布.采用不同的数值方法进行计算，并与理论预示以及文献[32]中公开发表的数值计算结果进行对比.计算工具为CFD++软件(TVD 格式)[36].物理模型为以空气为介质的理想气体定常层流流动.

为了方便对比，计算参数与文献[32]保持一致：圆柱半径8 cm，来流条件等效为海拔70 km 高度处大气条件：来流静压p=4.85 Pa,来流密度=7.48×10−5kg/m3，来流马赫数为3,6,12，其中不同马赫数对应的壁面温度Tw依次分别为300 K,500 K,1000 K.

计算使用对称半模型结构化网格，如图2 所示，在圆柱θ=90◦的位置和计算出口边界之间使用一段平直壁面作为过渡.网格数量为261(沿壁面)×125(壁面法向).沿壁面法向的第一层网格高度为1.5×10−5m.

图2 计算域与网格示意图Fig.2 Computation regime and mesh

不同马赫数条件下圆柱壁面的雷诺比拟系数计算结果见图3.根据Schwartzentruber 等[32]的计算结果，相比于DSMC 方法，采用N-S 方程计算得到的驻点热流偏大3.5%～10.1%，表明在70 km 高度处，8 cm 半径的圆柱绕流存在局部稀薄气体效应，处于近连续流动区.对于摩阻和热流的雷诺比拟系数Cf/Ch，采用DSMC 方法和N-S 方程计算得到的结果基本一致.

本文采用TVD 格式数值求解N-S 方程的结果也同时显示在图中.可见3 种数值方法计算得到的雷诺比拟关系较为一致.在高超声速(=6,12)条件下，不同数值方法，包括DMSC 方法以及不同数值格式求解N-S 方程计算得到的雷诺比拟关系均与理论预示的结果吻合较好.

图3 不同马赫数条件下的雷诺比拟关系Fig.3 General Reynolds analogy relation in different Mach number

对不同来源、不同方法的计算数据对比分析发现，对于超声速圆柱绕流中的雷诺比拟关系，采用基于粒子碰撞和统计抽样的DSMC 方法和基于连续介质假设求解N-S 控制方程方法计算的结果基本一致.在高超声速条件下，数值计算结果与理论预示结果吻合较好，而在较低马赫数(=3)条件下，数值计算结果与理论预示结果则存在一定区别，表明当前的高超声速广义雷诺比拟关系的理论结果还需要进一步改进才可以适应一般超声速流动的情形.

为了对计算结果进行确认，检查了本文采用的TVD 格式计算结果的收敛性以及不同网格数目对计算结果的影响.不同马赫数条件下驻点热流以及下游(θ=90◦)位置的热流计算结果收敛历程见图4.图中纵坐标为各时间步计算得到的热流值与最后一步计算的热流值之间的相对误差：（ .可见，计算过程中热流值稳步收敛，最终收敛结果的波动幅度小于10−7.使用不同密度网格计算得到的雷诺比拟系数结果见图5.三套网格对应的壁面附近第一层网格高度由疏到密各自为2.2×10−5m,1.5×10−5m 和1.0×10−5m.可见，在不同的马赫数条件下，网格数目对雷诺比拟系数的影响较小.

图4 热流系数计算相对误差随迭代步数的变化Fig.4 Heat Fluxes vary with iteration steps

图5 不同网格数目下的雷诺比拟系数计算结果Fig.5 Reynolds analogy coefficients under different amount of mesh

1.3 雷诺数对广义雷诺比拟关系的影响

在文献[32]提供的来流条件下，数值求解NS 方程以及DSMC 方法的雷诺比拟计算结果与理论预示结果均吻合较好.但是计算条件的来流密度和来流雷诺数较低.参考王智慧等[27,29]提出的驻点稀薄气体流动判据Wr=，对于海拔高度H=70 km，参考长度8 cm 的计算条件，马赫数=3，6，12 时，雷诺数分别为398，796，1592；Wr分别为0.013，0.018，0.026.对于计算的3 种状态，Wr参数的大小表明流动处于近连续区.而在常见的工程问题中，来流雷诺数往往大若干量级.为了在更一般的流动问题中考察广义雷诺比拟关系的适应性，在上述计算条件的基础上，计算了更大雷诺数条件下的圆柱绕流问题，计算条件见表1.计算选择高超声速的=6 和12 状态，在维持来流温度以及壁面温度条件下，通过提高来流压力条件增大雷诺数.

表1 计算条件Table 1 Simulation cases

图7 中显示了不同雷诺数条件下，驻点下游热流和摩阻系数的分布曲线.图中采用驻点热流系数和最大摩阻系数对热流系数和摩阻系数作归一化处理.由图可见，不同来流压力条件下，沿圆柱壁面的热流和摩阻分布规律基本一致，而在驻点下游较远的位置(θ >60◦)，不同来流压力条件下，归一化的热流和摩阻曲线出现了一定程度的差别.

图6 驻点热流系数计算结果Fig.6 Heat flux coefficients on stagnation point

图7 归一化的热流和摩阻分布Fig.7 Normalized heat fluxes and skin frictions

归一化的热流和摩阻分布曲线预示不同雷诺数条件下的雷诺比拟关系规律近似一致.由图8 可见，不同来流压力p，也就是不同雷诺数条件下，广义雷诺比拟关系具有如下规律：

(1)在θ ≤60◦时，不同雷诺数条件下的Cf/Ch曲线分布一致，数值计算的雷诺比拟关系与理论吻合较好.在θ=60◦时，数值计算与理论预示结果的相对误差为2.1%～6.5%；

(2)在θ>60◦时，Cf/Ch随雷诺数增大而增大，且偏离理论预示的线性关系.在θ=90◦时，数值计算与理论预示结果的相对误差最大为31.6%.

图8 不同来流条件下的广义雷诺比拟关系Fig 8 General Reynolds analogy relation with p=4.83～4830

针对不同雷诺数条件的计算结果表明，对于圆柱外形，即使是较大雷诺数条件，理论推导的线性广义雷诺比拟关系仍然可以较好地预示Cf/Ch分布.但是在比较靠近下游的区域，数值计算得到的广义雷诺比拟关系则与理论预示存在差异，且二者之间的差异随雷诺数增大而加大.

对于圆柱下游区域数值计算的雷诺比拟关系与线性理论预示结果的区别，可能是由于在靠近下游的部位，边界层的自相似程度发生变化，导致线性的广义雷诺比拟关系理论基础与实际情况有所区别.另外，圆柱壁面的热流分布随θ 快速降低，在靠近下游的区域，数值计算得到的热流相对误差可能随之加大，导致雷诺比拟关系的计算结果与实际情况出现差别.

2 幂次体壁面的广义雷诺比拟关系

理论上，广义雷诺比拟关系适用于一般的钝头体外形.除圆柱以外，本文另外选取一种幂次外形钝头体，采用数值求解N-S 的方法计算了不同来流条件下的壁面雷诺比拟关系.幂次体曲线方程为y=0.16x0.5,x=0～0.23 m.计算的来流条件与表1 一致.

幂次体外形与圆柱外形计算得到的流场马赫数云图如图9 所示.幂次体壁面附近的流动图像仍然是由脱体弓形激波主导，但由于外形相对更加尖锐，激波脱体距离比圆柱外形明显减小.

图9 圆柱和幂次体绕流马赫数云图Fig.9 Mach number counters on flows around a cylinder and power-law body

图10 幂次体壁面的雷诺比拟关系Fig.10 Reynolds analogy relations on a power-law body

不同来流马赫数和压力条件下计算得到的壁面雷诺比拟系数见图10.结果表明，在不同的雷诺数条件下幂次体壁面的雷诺比拟系数基本一致，且相比于圆柱外形，比拟系数更接近线性分布.另一方面，在幂次体壁面，雷诺比拟系数Cf/Ch沿θ 分布的斜率小于圆柱外形上分布.在=6,Tw=500 K 的条件下，式(2)预示的圆柱壁面的雷诺比拟关系为=1.942，而幂次体外形的计算结果则近似为≈1.48；=12,Tw=1000 K 条件下，圆柱壁面的结果为=1.522，而幂次体外形的计算结果近似为≈1.24.在θ=80◦时，数值计算与理论预示结果的相对误差最大为10.2%.

利用CFD 方法计算钝体壁面压力分布，根据式(1)分别计算得到的圆柱体和幂次体ue速度分布见图11.结果表明不同外形边界层外缘的速度分布存在一定差异，可能是导致不同钝头外形下雷诺比拟系数分布规律有所区别的原因.具体的影响机理还有待后续进一步的理论研究.

图11 边界层外缘速度分布Fig.11 Velocities along boundary layer edge

3 结论

本文采用数值求解N-S 方程的方法研究了钝头体壁面的广义雷诺比拟关系.采用不同的数值方法计算了不同来流条件和钝头外形下的壁面的热流、摩阻以及雷诺比拟系数的分布.数值计算结果与理论预示以及文献中公布的数据吻合.

通过改变来流压力条件，计算了不同雷诺数条件下圆柱壁面的热流和摩阻分布，数值计算得到的驻点热流与经典理论预示结果一致.不同来流条件下的归一化热流分布和摩阻分布在距离驻点下游较远的位置随雷诺数增加存在一定差异.受其影响，在不同的雷诺数条件下，壁面雷诺比拟关系在较下游的位置也存在一定区别，较高雷诺数下计算得到的Cf/Ch随θ 分布曲线偏离理论预示的线性分布规律.

本文同时计算了抛物线形钝头体壁面的雷诺比拟关系.在相同的来流条件下，相比于圆柱外形，幂次体壁面的雷诺比拟关系与线性规律吻合较好，且在不同雷诺数条件下的区别较小.受到钝头体外形对边界层外缘以及边界层内速度和温度的影响，不同外形壁面的雷诺比拟系数存在一定区别.

本文的数值计算和对比研究结果表明，理论预示的广义雷诺比拟关系存在于较大范围的雷诺数区间以及不同类型的钝头体外形中.在实际工程应用中，如果针对具体的流动条件和钝头外形进行修正，则可以进一步提高预示精度.