0 引言

临近空间高超声速飞行器的研发在中国及世界范围内已受到越来越多的关注，但同时也面临着巨大的技术挑战和新的科学问题^[1-3]。周恒和张涵信在他们所写的《空气动力学的新问题》一文^[1]中分析了在情况下现有的空气动力学的不足，提出了若干需要考虑的新问题，其中之一就是如何考虑流场中可能出现的局部稀薄气体效应。他们分析了临近空间高超声速飞行器的边界层，指出：由于边界层中温度很高，再加上临近空间飞行器的飞行高度很高，因此边界层中的气体密度可能很低，即气体分子自由程较大。另一方面，边界层中的剪切很强，在这两个因素的影响下，边界层中法向相隔一个分子自由程的两层流体的速度差和当地的分子运动的平均速度相比不是一个小量，这会导致分子运动的速度分布函数显著偏离Maxwell分布，使得现有的各种基于实际分子自由运动速度的分布函数对Maxwell分布仅有小的偏离假设下提出的分析计算方法不再可靠，需要寻找新的方法。根据这一思想，陈杰和赵磊在文献[2]中提出了对应的参数Zh，用以刻画气体稀薄效应对剪应力的影响。

参数Zh定义为：相隔为一个分子自由程的两层流体的宏观速度差和分子热运动速度之比。这是一无量纲参数，即：

$ Z h=\lambda\left|\frac{\mathrm{d} u}{\mathrm{d} x}\right| / c_{m} $

(1)

其中:λ是分子自由程，|du/dx|是当地宏观速度u在其垂直方向的梯度，c_m是当地分子自由运动最可几热速度。文章通过DSMC数值模拟证实了随着Zh增加，分子速度分布函数对Maxwell分布的偏离的确会越来越大。同时获得了线性本构关系计算的剪应力随Zh增大出现的系统性误差。随后，又提出了通过在CFD中修正黏性系数以考虑局部稀薄效应的思路，该方法的可行性和对气动力的改进在超声速平板问题上得到了验证^[4-5]。而且重要的是，该方法无需改变现有的CFD方法基本框架，所增加的计算时间只不过30%左右，其计算成本远低于基于气体动理论的计算方法，有望应用于工程实际问题中。

在临近空间高超声速飞行器的设计中，不仅要考虑气动力，还要考虑气动热^{[3, 6-10]}。对于长时间在近空间中以高超声速机动飞行的飞行器，精确的热环境预测对低冗度热防护设计十分重要^[6]。临近空间高超声速飞行器气动加热受稀薄气体效应、高温真实气体效应、飞行器复杂外形带来的复杂流场结构等多种因素的影响^[7-10]。工程上往往针对驻点热流基于半经验参数构建桥函数等方式考虑稀薄非平衡计算^[11]。在数值计算方面也发展了多种考虑稀薄效应的宏观方程模型^[11-15]。

本文将在之前研究强剪切流中稀薄效应对应力计算工作的基础上^{[2, 4-5]}，进一步考虑气体稀薄效应对热流计算的影响，将研究以下问题：1)存在强温度梯度时，气体的稀薄效应如何分析和处理，分子速度分布函数是否还是Maxwell分布，如果不是，刻画偏离程度的无量纲参数是什么；2)如果该参数的值大了，会如何影响热流计算？傅里叶导热定律是否成立？是否能定量地找到等效导热系数和此参数的依赖关系。研究的路线和文献[2]相同，即以DSMC为手段寻找气体有稀薄效应下的热传导规律。3)初步验证将该规律纳入CFD计算中是否可以有效改善受稀薄效应的热流计算？

1 计算模型 1.1 一维导热问题

为单纯研究存在强温度梯度时气体稀薄效应对热流计算的影响，选取无宏观流动的导热问题进行分析。我们这里想要获得的是对于具有一定稀薄程度的气体在存在强温度梯度时，傅里叶定律所计算的热流误差，而不考虑壁面引起的效应。如图 1(a)所示，该研究对象为距壁面较远且存在强温度梯度的流体，即在厚度L的薄层内存在较大的温度差T_h-T_c。经典的两平板间导热问题如图 1(b)所示，相距L的两板分别保持不同的温度T_c、T_h，壁面边界条件通常采用Maxwell完全漫反射条件。对于图 1(b)所示问题，可基于两板间距定义传统Knudsen数Kn_L=λ/L。如果Kn_L较小，如Kn_L < 0.1, 两板间中心区域可认为不受固壁边界条件影响，等价于图 1(a)中所要研究的薄层。但随着Kn_L增大，壁面效应越明显，中心区域的热流计算将既受到强温度梯度的影响，也受到壁面影响。虽然可以通过增大两板间的距离减弱壁面效应，但是当两板距离L增大后，如在中心区域实现相同的温度梯度dT/dx，两板的壁面温度之差就要相应地增大，很难实现。如保持中心温度不变，一端温度甚至可能小于0°，显然不合理。如要求一端保持不低于0°，则另一端温度就非常高。即便这样也不容易实现较大Zh_T时的情况。因此本文将通过在虚线处施加合适的边界实现图 1(a)中所示问题的模拟。

1.2 DSMC方法

本文采用基于Bird提出的标准DSMC算法^[16]建立的自编程程序。作为第一步，为避免高温下可能引起的分子振动能被激发，甚至导致分子离解等复杂问题掩盖气体稀薄效应的本质问题，本文考虑单原子气体氩气作为模拟气体，采用硬球模型。氩气分子质量为6.63×10^-26 kg，直径为3.659×10^-10 m。

DSMC方法使用大量的模拟分子模拟真实的气体，用一个模拟分子代表一定数量的真实气体分子，该方法的本质是在时间步长Δt内，对分子的运动和分子间的碰撞进行解耦，要求时间步长小于平均碰撞时间，网格尺度小于平均自由程。因此，为保证计算有效导热系数的足够精确，采用了文献[17]中建议的网格大小和时间步长。为获得统计噪声较小的热流以及更准确的等效导热系数，模拟使用了大量的模拟分子，初始每个网格内分子数为1×100个，对导热稳定后的10⁶个时间步进行统计平均。

1.3 稀薄效应刻画参数

周恒和张涵信^[1]指出：相隔一个分子自由程的两层间的分子碰撞是产生剪切力和热传导的主要机制。在剪切流中，刻画稀薄效应的无量纲参数Zh为一个分子自由程上宏观速度的差与微观分子最可几速度之比。该参数有效表征了剪切流中的非平衡。根据同样的思想，这里将该参数类推至温度，流体的温度是分子平均热速度的统计。如果一个分子自由程上分子热运动速度的改变量与热运动平均速度本身相比不再是小量，那分子速度分布也不再会是Maxwell分布。

因此，定义Zh_T为相隔一个分子自由程的两层流体的特征分子热速度梯度和当地的平均热速度之比为Zh_T：

$ Z h_{T}=\lambda\left|\frac{\mathrm{d} \bar{c}}{\mathrm{d} x}\right| / \bar{c} $

(2)

如果进一步考虑$\bar{c}=\sqrt{8 R T / {\rm{ \mathsf{ π} }} m}$，则上式可写为：

$ Z h_{T}=\frac{1}{2} \lambda\left|\frac{{\rm d} T}{{\rm d} x}\right| / T $

(3)

需要指出的是Zh_T的形式虽然与文献[18]中的Kn_GLL-T只差一个常数，但是意义不同。这里是两个速度或者温度之比，而文献里Kn_GLL-T的意义为两个长度之比，即分子平均自由程与温度的标尺长度之比。Zh_T参数的物理意义将在2.1节中通过对速度分布函数的分析进一步详细解释。

1.4 边界条件

这一节探讨边界条件的实施。如1.1节所示，我们需要采用合适的边界条件模拟距壁面较远且存在强温度梯度的一层流体。在稀薄强温度梯度条件下，气体显然不可能处于平衡态，即在DSMC模拟的边界处不能通过平衡态Maxwell分布向计算区域加入粒子。为此，我们提出了新的边界条件以实现边界处所设定的温度以及粒子所处的非平衡特性。

该边界条件的主要思想是对于不受壁面影响的导热问题，模拟区域边界处与中心点具有相似的非平衡性，即分布函数形状以及偏离Maxwell分布的程度相同，因此通过中心点的粒子信息在边界处进行重构，使其与中心点具有相似非平衡态分布函数，但两组粒子的统计温度不同。在下文将此边界称为“非平衡克隆边界”。其具体实现方法分以下几个步骤：在每一次时间步推进后，1)统计边界网格中存在的粒子数，用N表示；2)从中间网格随机抽取N个分子复制到边界网格。对低温端，由于对应的温度低，所以一般N要大于中间网格的粒子数，需要另外随机抽取若干粒子，补足差数；3)统计抽取的这N个粒子重新组成的粒子团的平均速度。因为随机性，重新组成的粒子团平均速度可能不再是0，有很微小的偏差。因此需要将这个小的偏差速度减掉，以免造成边界处不正确的动量输入；4)统计N个粒子的温度T₁，根据壁面设定的温度T_w和这个统计温度比重整粒子的速度，即粒子速度乘以$\sqrt{T_{\mathrm{w}} / T_{1}}$。这样，壁面网格内的分子平均速度为0，温度为T_w，分子的速度分布函数与中间网格类似。

为验证此边界，我们分别模拟了图 1所示两种情况。对于图 1(a)所示模型, T_c=700 K，T_h=1200 K，L为中心点平均自由程的5倍，两边均用克隆边界。对于图 1(b)所示模型, T_c=400 K，T_h=2000 K，中心点附近温度为1200 K，为避免壁面影响，L为中心点平均自由程的10倍。图 2给出了两个计算模型相同温度1200 K处的概率密度分布函数(Probability Density Function, PDF)，并与1200 K的Maxwell分布对比。可以看到图 1(a)模型中克隆边界条件所实现的分布函数与图 1(b)模型中心点的温度函数重合，且相同程度地偏离了平衡态，验证了克隆边界可以重构不受壁面影响的非平衡态分布函数。峰值微弱的差别来自于温度的差，克隆边界处温度为1200 K，中心点附近所取的点温度为1197 K。

1.5 等效导热系数

本文最终要获得的是Zh_T对热流计算的影响，求得等效的导热系数，其定义为：

$ k_{e}=\frac{q_{\mathrm{DSMC}}}{\mathrm{d} T / \mathrm{d} x}=A_{k}\left(Z h_{T}\right) k(T) $

(4)

即DSMC获得的热流比上当地的温度梯度，可理解为气体在线性本构关系假设下所表现的等效导热系数，它是传统的导热系数乘以修正系数A_k(Zh_T)。

2 结果分析 2.1 速度分布函数

我们计算了不同Zh_T的参数下计算区域中心点的概率密度分布函数(PDF)。随着Zh_T的增加，分布函数逐渐偏离Maxwell分布。不同于剪切流中分布函数呈现单峰、双峰、三峰等多种形态^[2]，导热问题中不同Zh_T的参数下分布函数偏离的形态是一致的，如图 3所示，只是偏离程度不同。图 3给出了Zh_T=0.03和Zh_T=0.1时的概率密度分布函数。Zh_T=0.03时分布函数的峰值已经开始向右倾斜，Zh_T=0.1则明显偏离。

从分布函数与Maxwell分布对比来看，粒子运动的最可几速度不再是零而变为正。其原因可分析如下：设图 4中的A、B线分别为一个网格的两个边界，温度梯度方向为从A向B，即B处的温度高于A处的温度，用N_A、N_B分别表示在网格线上进出的粒子数目。由于气体宏观静止，因此在网格线处，单位时间从左往右和从右往左穿过的粒子数应相等，即图中所示的N_A⁺=N_A^-，N_B^-=N_B^-。设c为分子热运动平均速度，则可认为在一个时间步长dt内，位于B或A两侧一个宽为cdt的范围内的正向或反向运动的粒子在dt时间内均可分别向左或向右穿过网格线，对初始的Maxwell分布来说，正向和反向的粒子数应相等。而位于该范围的粒子数则和密度与c成正比。由于密度和温度成反比，而c则与$\sqrt {T} $成正比，因此在dt时间内穿过网格的粒子数应该和1/ $\sqrt {T} $成正比。由于B处的温度比A处的温度高，所以B处穿过网格的粒子数小于A处穿过网格的粒子数，也就是图中所标示的N_B^- < N_A⁺，即有净流入的正向粒子和净流出的反向粒子，使得分布函数偏于正的一侧。

根据以上分析，单位时间内一个网格的正向粒子净流入量和通过中心界面的正向粒子流量之比，即：

$ \lambda\left|\frac{\mathrm{d}(1 / \sqrt{T})}{\mathrm{d} x}\right| /(1 / \sqrt{T}) $

(5)

可以刻画分布函数偏离Maxwell分布的程度。而这一参数和公式(3)的定义一致。

我们还可从DSMC的计算过程中做统计以证实这一趋势。在计算中心点网格温度为1000 K，网格内模拟粒子数在100个左右，等导热稳定后我们对中心点网格低温侧和高温侧进出的粒子在1×10⁶时间步上进行统计平均。如图 3(a)所示的Zh_T=0.03的情况，中心点网格低温侧进出的粒子均为5.4088个，而高温侧进出的均为5.3893，即网格内总的进出粒子相等，但低温侧进入的比右侧多0.36%。图 3(b)所示的Zh_T=0.1的情况，低温侧进入的比高温侧多0.42%，当Zh_T=0.36时低温侧进入的比高温侧多0.52%。因此，Zh_T越大，网格中来自低温侧的粒子比来自高温侧的粒子越多。不过DSMC的结果是最终平衡态的结果，当数量更多的正向粒子进入网格后，会和反向运动的粒子发生碰撞，使得正反向的粒子数趋于相等(正向粒子数虽然较多，但平均速度较小，网格内粒子之间相互碰撞的总的效果是使得粒子数和大小趋于相等，和输入粒子的作用相反)。这是一个十分复杂的过程，很难用简单的物理过程说清楚。

2.2 壁面效应分析

本节主要分析壁面效应的影响，以及验证克隆边界是否能够有效消除壁面影响。对同一问题，如果使用完全漫反射条件，即所模拟的问题为受限于间隔距离很小的两平板间的导热问题，壁面附近气体的温度与给定的壁面温度不相等，即所谓的温度跳跃，壁面附近努森层中的温度分布也呈现非线性变化，如图 5黑线所示。如果采用克隆边界，温度分布则如黑点所示，可基本实现所施加的温度梯度，且成近似线性分布。在靠近壁面的一个点处有一个小的间断，但这并不影响我们在中心区域的分析。产生该小的间断的原因是，在计算中，网格尺度是当地自由程的1/5以下，按理一个网格两侧5个网格内的粒子都有可能不经过粒子间碰撞直接到达该网格对热量交换做出贡献。但边界网格内的那个网格在边界方向只受到临近一个网格的影响，而损失了一部分粒子对它的影响。比如低温端右侧的第二个网格，它左侧只受到第一个网格内进来的粒子的影响，而实际情况应该是有从更远处的网格中更低温度的粒子参与能量交换。失去了这些更低温度粒子的参与，自然要求相邻网格的更大温差来弥补。高温端的情况类似。

我们还进一步分析了壁面效应对于等效导热系数计算的影响。为此，在固定Zh_T的同时，改变两壁面间的间距，看大到什么程度后其等效导热系数不再改变(大小以间距L相当于多少个中心位置的自由分子程数λ₀计，即Kn_L=λ₀/L)。对初始Zh_T分别为0.04和0.08的情况采用漫反射边界和克隆边界做了两组算例。图 6给出了所得的导热系数k的修正值A_k随Zh_T的变化。图中空心圆和空心三角为完全漫反射边界计算的结果，并标出了相应的Kn_L值。红色实心圆和实心三角为克隆边界计算的结果。对比发现，漫反射边界计算所得的A_k随两板间距减小而急剧减小，而隆边界条件计算的A_k的值则几乎不受两板间距的影响。对我们所关心的情况，用克隆边界条件时，计算域最多只要2.5倍的中心位置处的分子自由程就够了，而对两端为固壁的情况，则即使Zh_T仅为0.04，计算域也要达到10个分子自由程以上才可以。而由于温度有不能小于零的限制，对两端固壁情况，有时不可能满足要求。

这部分结果说明：采用壁面漫反射条件在计算域不够大时不能获得无量纲参数Zh_T和修正系数A_k之间的规律性变化；而克隆边界可有效消除固壁影响，只需要较小的计算域就能消除壁面对热流计算产生的影响。对我们研究的情况，只要两个边界间距离大于2.5倍中心处的分子自由程，即可获得相对可靠的结果。

2.3 导热系数修正曲线A_k(Zh_T)

我们计算了不同中心点温度、不同温度梯度的50个工况，Zh_T在[0, 0.4]范围内变化，所得的导热系数修正曲线如图 7所示。图中不同颜色符号代表了不同的中心点温度，从1000 K到2000 K变化，我们可以看到，该曲线不依赖于当地的温度，只要Zh_T相同，所获得的导热系数修正值就相同。所有工况的导热系数修正值A_k随着Zh_T遵循统一的规律，呈单调递减变化，也就是等效导热系数比传统的导热系数偏小，傅里叶定律计算的热流值偏大。

3 圆柱绕流

为初步验证上述规律是否可以有效提高对驻点热流的预测精度，分别采用DSMC和考虑导热系数修正的CFD模拟了单原子气体氩气绕半圆柱的流动。CFD模型的计算框架及计算格式详见文献[4, 5]。模拟算例为文献[19]中Kn=0.05的工况, 圆柱半径为0.1524 m，壁面温度T_w=500 K，来流温度T_∞=200 K，来流气体的数密度为8.494×10¹⁹/m³，来流马赫数分别为Ma=10和Ma=25。

高超声速绕圆柱流动在驻点附近存在强的温度梯度，流体沿圆柱表面向后运动时又存在强的剪切。因此，该问题需要同时考虑强剪切和强温度梯度引起的稀薄气体效应的影响。在计算中，根据流场当地的参数Zh和Zh_T对导热系数进行了修正。2.3节中获得了Zh_T对导热系数k的修正值A_k(Zh_T)，Zh对导热系数k的修正值A_k(Zh)在文献[2]中获得。在此算例中，两者以加权平均的方式获得最终的导热系数的修正值，即：

$ A_{x}=\frac{Z h}{Z h+Z h_{T}} A_{k}(Z h)+\frac{Z h_{T}}{Z h+Z h_{T}} A_{k}\left(Z h_{T}\right) $

(6)

需要说明的是这里作为初步尝试采用了加权平均的方式考虑参数Zh和Zh_T对导热系数的影响，主要基于流动图像，认为在驻点附近应是Zh_T起主导，沿圆柱表面向后运动的边界层中应为Zh起主导。更加合理的方式将在今后的工作中探讨。

图 8和图 9分别给出了Ma=10和Ma=25时圆柱表面热流系数的计算结果，同时给出了完全漫反射边界条件计算的DSMC结果，以及使用滑移边界条件的N-S方程和不使用滑移边界条件的N-S方程结果。横坐标ϕ为从驻点开始沿圆柱表面顺时针旋转的角度，纵坐标为无量纲的热流系数C_H=q/0.5ρ_∞U_∞³。从图中可以看出，无论是否采用滑移边界，传统的N-S方程计算的热流均高于DSMC结果。而通过本文所建立的导热系数修正算法对不同马赫数的圆柱表面热流，尤其驻点附近热流预测有了很大的改善。

4 结论

本文以临近空间高超声速飞行器气动热计算为背景，考虑具有一定稀薄程度的气体在存在强温度梯度时热流的计算。主要获得以下结论：

1) 提出了表征该问题的稀薄效应刻画参数Zh_T；

2) 获得了依赖参数Zh_T的导热系数修正规律；

3) 通过导热系数修正可有效提高高超声速圆柱绕流壁面热流的预测精度。

值得说明的是，本文仅考虑了一种简单的情形，即单原子气体不存在速度梯度和压力梯度，仅存在强温度梯度的情况下对热流计算的影响。而实际飞行器问题则更为复杂，驻点附近可看作纯的导热问题，但是过了驻点后，法向还存在很大的速度梯度、压力梯度，并且伴随高温的振动能激发，电离离解等效应，需要在今后的研究中逐步分析，以解决真实飞行器流动中的问题。另外，飞行器的壁面条件应如何处理，也是要专门研究的难题。

致谢: 感谢周恒院士的指导和通过有启发性的讨论给予的帮助。