网格对高超声速钝头体表面热流数值模拟结果的影响

马崇立，刘景源

南昌航空大学 飞行器工程学院，南昌 330063

发展具有一定高超声速巡航能力、可重复使用的新型天地往返运输系统及在大气层内飞行的高超声速飞行器已成为航天航空强国技术储备中需要研究的关键课题，其中高超声速复杂外形飞行器表面热流精确预估是其中的一个核心问题。

对高超声速飞行器表面热流精确预估有试验、理论分析及数值模拟等3 种方法。但是无论高超声速飞行器的地面试验还是飞行试验，均存在成本高、周期长等问题。并且地面试验条件的限制无法模拟真实飞行环境；而飞行试验虽然可以模拟飞行环境，但有投入巨大、可重复性弱等缺点。由于高超声速流动方程的非线性性质以及高超声速飞行器流场的复杂性，理论分析无法给出可用于指导工程设计的结果。因此，数值模拟技术成为高超声速飞行器表面热流精确估计的重要手段。

高超声速飞行器头部激波强、流动减速剧烈，气动加热严重。为了有效防热，高超声速飞行器采用钝头型。因此，数值模拟精确预估高超声速钝头体飞行器表面热流即成为高超声速复杂外形表面热流必须先行解决的问题。

数值模拟预估高超声速钝头体的精度受模拟所用网格、选取的离散方法、边界条件及数值格式等的影响［1-2］。而钝体表面法向第1 层网格尺度被认为是重要的影响因素［3］。文献［4-5］给出了基于壁面参数的壁面法向网格雷诺数，认为壁面法向网格雷诺数对气动热的预测具有重要的影响，并经数值分析指出，为了精确预估壁面热流，壁面法向网格雷诺数必须小于一定的数值。Klopfer 等［6］提出了基于来流参数的法向网格雷诺数，数值模拟结果表明，为了精确预估高超声速钝体表面热流，基于来流参数的法向网格雷诺数亦必须小于一定的数值。潘沙等［7］研究了基于来流参数的法向网格雷诺数对高超声速钝头体热流的影响，并给出了法向网格划分的建议。张翔等［8］应用数值模拟方法研究了高超声速二维半圆柱绕流，并分析了法向网格雷诺数与激波位置网格尺度对壁面热流的影响。闫文辉等［9］应用数值模拟方法分析了二维半圆柱超声速绕流的摩阻及热流，结果表明，为准确计算摩阻及热流法向第1 层网格间距、网格长宽比与扩张度等需要取合适的值。Qu［10］系统地研究了几种二阶迎风数值格式在计算高超声速壁面热流时，对限制器及网格雷诺数的依赖，结果表明，SLAU数值格式计算高超声速热流具有高精度及强鲁棒性。Zhao 等［11］基于数值模拟方法研究了不同网格雷诺数对高超声速化学非平衡绕流流动的飞行器表面及驻点热流的影响，结果表明，随网格雷诺数减小，所计算的飞行器表面热流及驻点热流均增大并最终计算的热流具有网格无关性。Qu［12］研 究 了Roe、AUSM+及AUSMPW+这3种数值格式计算高超声速真实飞行器表面热流的模拟精度，结果表明，为精确计算飞行器表面热流，基于来流参数的网格雷诺数不应大于20。

由于缺少理论上的支撑，上述方法确定的壁面法向网格尺度具有一定的经验性及不确定性。程晓丽等［1］则提出将分子平均自由程作为壁面法向网格尺度，数值模拟结果表明，可将该尺度作为法向网格划分的最优下限。张智超等［2］则进一步发展了程晓丽等提出的方法，改善了基于分子平均自由程准则需要通过粗网格试算获得驻点位 置 壁 面 密 度 的 缺 点。Yang 与Liu［13］亦 基 于 分子运动论，引入平均自由程给出了一种驻点附近法向网格划分方法，并应用所给出的方法数值模拟了高超声速半圆柱绕流热流流场，结果表明，为准确计算高超声速表面热流，壁面法向第1 层网格的尺度应为基于壁面参数的分子平均自由程的大小，并且所提出的方法不但适用于介质为空气的高超声速绕流，对氮气亦适用。Ren 等［14］根据一维流动理论，给出了基于来流及壁面参数的法向网格雷诺数及基于壁面参数的分子平均自由程三者在计算高超声速连续流及过渡区时的关系，并指出为精确计算壁面热流，壁面法向第1 层网格尺度必须严格满足一定的网格准则。Huang 与Hu［15］基于Kemp-Riddell 经验公式，提出了计算高超声速飞行器表面热流的壁面法向第1 层网格划分的新准则，并基于新准则计算了球头高超声速绕流表面热流，结果表明，所提出的新方法，比基于网格雷诺数的方法给出的驻点热流更精确。张亮等［16］则给出了一种适用于高超声速高冷壁下表面热流法向网格尺度上限确定方法。

对高超声速钝头体表面热流进行精确预估的绝大多数的研究，均关注表面法向网格的划分方法，而流向网格尺度对气动热的影响研究有所欠缺。文献［17-18］虽然基于差分格式误差匹配原则，提出了网格划分依据，但该方法仅基于Navier-Stokes 方程组的流向动量方程，忽略了对高超声速流动的能量方程的分析，而能量方程对高超声速热流预估具有重要影响。另外，该方法需要事先给定基于自由来流参数的法向网格雷诺数，而法向网格雷诺数的选择具有一定的不确定性的缺陷。

本文基于法向网格划分方法、高超声速能量方程及动量方程修正方程各误差项的误差匹配，提出了基于壁面参数的法向网格雷诺数及基于钝头体特征长度的来流雷诺数作为网格划分的两个依据，从而给出了一种对高超声速钝头体表面绕流的较为完整的网格划分方法，然后应用所提出的方法对二维半圆柱、球头等高超声速钝头体绕流进行了数值分析及讨论，并指出本文方法不足之处及进一步的研究方向。

1 数值格式的误差匹配

应用数值模拟方法对高超声速钝头体表面热流进行预估，必然与流场控制方程的精确解存在一定的误差，而合理控制数值模拟的误差对精确预估飞行器表面热流具有重要意义。由于高超声速飞行器表面存在边界层，所以在流场数值计算时边界层内计算网格的长宽比较大，而计算误差与网格的长宽比密切相关。保证控制方程的差分格式的修正方程各误差项为同一量级［18］，是预估飞行器表面热流的前提。

1.1 能量方程误差匹配

二维定常无量纲高超声速能量方程［19］的差分格式的修正方程为

式中：为密度；为比定压热容；、分别为流向及法向速度分量；为温度；γ为比热比；Ma∞为来流马赫数；为压强；Pr=0.71 为普朗特数；为热传导系数；、、、分别为流向及法向与温度相关的误差项；为能量方程中压力离散误差项；为能量方程耗散函数离散误差项；Re为基于飞行器特征长度的自由来流雷诺数。为耗散函数，表达式为

为方便起见，以下分析及论述忽略上述方程符号上标“-”，而用无上标量表示无量纲量。

不失一般性，将无量纲的热传导系数k当作近似常数处理，式（1）中流向与法向温度相关项的误差项分别为

式中：a1、a2、b1、b2、c1、c2、d1、d2为差分格式的系数。

式（1） 中RΦ=Rux+Rvy+Rvx+Ruy+Ruxvy+Rvxuy，等号右端各项的具体形式如下：

式中：g1、g2、h1、h2、r1、r2、s1、s2、e1、f1为差分格式的系数。

用δ表示边界层内流动的特征尺度，根据边界层理论，可得无量纲参数量级关系为［18，20］

式中：上标p、q、m、n为各物理量求导阶次。

1.1.1 能量方程温度项误差匹配

本节首先对能量方程式（1）中温度相关项进行分析。在数值计算中，为了能够准确捕捉到壁面附近的流动参数，需要将边界层内沿法向方向进行加密，即边界层内法向网格尺度Δy要小于δ；所以Δy δ应为小量，从而与温度相关的误差项可简化为

为了保证数值计算的精度，并使法向与流向2 个方向所选择的网格尺度最优，各误差项应保持为同一量级［18］，即式（7）～式（10）左侧误差项的比值量级为1。对以上各误差项进行匹配，可得如下结论：

1）RTx1与RTy1匹配，则

2）RTx2与RTy1匹配，则

3）RTx1与RTy2匹配，则

4）RTx1与RTy2匹配，则

在数值计算中，若流向与法向差分格式采取相同的阶数，即a=b=2c=2d；故网格长宽比量级应取式（7）～式（10）匹配结果中较小量以保证边界层内网格尺度满足数值计算需求，即

1.1.2 能量方程耗散项误差匹配

耗散项与能量方程的温度项相比，对飞行表面热流的影响亦重要，因此为了保证数值计算的精度，能量方程耗散项中流向与法向的误差项也应为同一量级。

根据边界层理论以及Δy δ应为小量，结合式（5），耗散项各误差项可记为

同1.1.1节，式（12）各误差项亦应保持为同一量级。对以上各误差项进行匹配，可得如下结论：

1）Rux与Ruy匹配，则

2）Rvx与Rvy匹配，则

3）Rux与Rvy匹配，则

4）Rvx与Ruy匹配，则

若耗散项流向与法向差分格式采用相同的阶数，即g=h=r=s，为保证数值计算精度，网格长宽比量级应取上述匹配结果中较小量以保证边界层内网格尺度满足数值计算需求即匹配式（15），其长宽比量级关系式同式（11）。

不失一般性，对耗散函数各导数项的离散均采用相同的差分精度，因此式（12）最后两项量级相同，并且与式（12）其他4 项进行误差比较时，所得结论包含在式（13）～式（16）内，限于篇幅，此处不再赘述。

1.2 动量方程误差匹配

对动量方程的分析与文献［18］类似，二维定常无量纲高超声速边界层流动的动量方程差分格式的修正方程为

式 中：、分 别 为 流 向 与 法 向 对 流 项 误 差项；、分别为流 向及法向 黏性项误 差项。

同1.1 节，以下分析及论述忽略上述方程符号上标“-”，用无上标量表示无量纲量，则式（17）中误差项分别为

式中：a*1、b*1、c*1、c*2、d*1、d*2为动量方程差分格式的系数。

根据边界层理论与式（5），各误差项可简化为

为保证数值计算精度，对流项与黏性项误差项比值量级为1。对以上各误差项进行匹配可得：

1）Rux1与Ruy1匹配，则

2）Rux2与Ruy1匹配，则

3）Rux1与Ruy2匹配，则

4）Rux2与Ruy2匹配，则

若流向与法向差分格式采用相同阶数，即a*=b*=2c*=2d*，网格长宽比量级应取式（20）～式（23）匹配结果中较小量，即

1.3 能量方程与动量方程误差匹配结果

通过对能量方程与动量方程进行误差分析可得，数值计算中为保证数值计算精度，网格长宽比量级上应取较小量以保证边界层内网格尺度满足数值耗散较小的要求。分别对温度项与耗散项误差匹配中长宽比量级较小项，即式（7）、式（9）、式（15）、式（20）及式（22）等的分析表明：①当流向与法向导数各项采用相同差分精度时，网格长宽比应与飞行器特征长度的自由来流雷诺数平方根同量级；②当流向各导数采用高阶差分精度时，如法向第1 层网格间距确定，为了满足误差匹配，则流向网格尺度应增大即放宽对流向网格尺度的要求；③当法向采用高阶精度时，流向网格尺度应降低。

当对流项、温度二阶空间导数及耗散函数项等均取相同差分精度时，网格长宽比应与基于钝头体特征长度的来流雷诺数的平方根同量级，即式（11）或式（24）。当能量及动量方程的对流项、温度二阶空间导数及耗散函数项等取不同差分精度时，由于网格长宽比量级应取上述匹配结果中较小量，网格长宽比的结论亦为式（11）或式（24）。

对于三维问题，由于2 个切向方向的截断误差表达式类似，所以只需要将上述匹配结果中的流向网格尺度Δx用切向网格尺度hτ=max(Δx， Δz)代替即可把本节的结论推广到三维数值模拟中。

另 外，应 用 式（5）的 量 级 估 计，式（1）及式（17）与压强相关项的误差项进行误差匹配的结果表明，流向与法向网格尺度比为来流雷诺数。即所得结果不对本节结论产生影响，具体分析此处不再赘述。

2 高超声速钝头体绕流流场网格划分

第1 节的分析虽然给出了流场流向与法向网格划分时，应考虑流向与法向网格尺度的比，但数值模拟中还不能确定流向及法向等方向的网格尺度的具体数值。若能确定流向或者壁面法向的网格尺度，则高超声速绕流流场的网格划分尺度即可确定。

由壁面网格雷诺数的定义及壁面分子平均自由程的公式［20］可知

式中：aw为基于壁面参数的声速。

由式（25）及式（26）可知

由程晓丽等［1］给出的壁面法向第1 层网格为壁面上分子平均自由程，因此基于壁面参数的网格雷诺数为1.5。

根据文献［21］对流域划分，则为了使Navier-Stokes（N-S）方程组成立，努森数Kn＜0.1。根据文献［22］，上述流域划分方法可推广到数值模拟的壁面附近的网格划分中，则对于数值计算：

由式（27）及式（28）可知

但考虑到数值计算中并非直接求解N-S 方程组，而是求解其修正方程组，因此实际计算中，应结合数值格式的精度确定基于壁面参数的网格雷诺数选取的具体数值。

因此，本文提出基于壁面法向网格雷诺数Recell，w，与基于绕流流场特征长度的自由来流雷诺数相结合确定绕流流场网格法向及流向网格尺度，用于数值模拟前网格划分的参考。具体为，首先由来流及壁面参数，根据指定的Recell，w具体数值，确定法向网格尺度，然后再根据流向与法向网格尺度比确定流向网格尺度。周向网格尺度应用1.3 节的三维推广即可。

3 算例验证

为验证第2 节所提出的网格划分方法，本节应用完全气体模型，对高超声速半圆柱、球头等进行数值模拟，对钝锥及升力体等外形引用文献给出的经验网格划分，对本文提出的网格划分方法进行验证。数值模拟均采用三维完全气体NS 控制方程组；对流项离散采用效率较高且相对稳定的van Leer 提出的迎风数值格式，应用MUSCL 方法使格式精度为3 阶，并采用文献［23］给出的熵修正方法减小该格式的数值黏性耗散，黏性项采用二阶中心差分格式，时间导数项采用隐式推进，以提高计算效率。

3.1 二维半圆柱绕流

3.1.1Ma∞=8 半圆柱绕流

首先选择文献［24］的圆柱绕流实验作为验证算例。圆柱直径D=0.076 m，自由来流马赫数Ma∞=8，静压P∞=855 Pa，静温T∞=125 K，壁面温度Twall=294 K，基于圆柱直径的来流雷诺数为Re=3.76×105，限于篇幅，圆柱绕流几何模型可见文献［7， 13］。根据第2 节给出的流场网格划分方法，可得周向与法向第1 层网格尺度比为ΔxΔy～O(613)；再基于壁面参数的网格雷诺数Recell，w=15，代入自由来流参数，可给出Δnw=7.16×10-7m，由第2 节给出的方法可得周向网格尺度为Δx=4.38×10-4m，对应的周向网格节点数约为270。

根据文献已有的数值模拟对网格划分方法的建议，可知用第2 节流场网格划分方法，给出的半圆柱周向网格点数过多。究其原因，第2 节给出的方法，首先需要确定法向网格尺度，才能根据周向与法向网格尺度比，再确定周向网格尺度。若确定的法向网格尺度过于保守（网格尺度过小），必然导致给出的周向网格尺度也偏于保守，从而致使周向网格点数过多。

为验证上述的论述及分析，进行以下数值模拟。第1 步，首先给定周向足够多的网格点数，然后改变不同的法向网格尺度，应用数值模拟与实验结果进行对比，用数值方法找到此算例的最大法向网格尺度，确保此法向网格尺度的网格无关性及与实验结果的一致性。第2 步，再改变周向网格尺度（网格点数），验证本文提出的周向网格尺度与法向网格尺度的比是否适用。

图1 为不同法向网格尺度下，半圆柱表面斯坦顿数St周向分布图。由图可见，在确保周向网格尺度足够小时，当法向网格尺度小于8×10-6m时，圆柱的斯坦顿数几乎不随法向网格变化而变化，根据文献［1］的方法得到的Δnw=1.05×10-7m 过于保守，即应用分子平均自由程确定法向网格尺度偏保守。而本文基于壁面参数的网格雷诺数Recell，w=15 给出的法向网格尺度虽然有所改进，但还是过小。

图1 Ma∞=8 时法向网格尺度对斯坦顿数的影响（半圆柱）Fig.1 Influence of normal grid scale over wall on Stanton number at Ma∞=8 （semi-cylindrical）

根据数值试验找到的法向网格尺度为Δnw=8×10-6m，然后再根据第2 节给出的周向与法向网格尺度比可得，Δx=4.9×10-3m，即半圆柱周向节点数为25。为验证本文给出的周向网格是否合理，分别选用壁面法向第1 层网格尺度为Δnw=8×10-6， 4×10-6m，通过改变周向节点数的斯坦顿数分布结果如图2 所示。由图可见，当法向网格取值合理时，用本文的方法给出的周向网格的取值亦合适。并且从图2（b）可见，当法向网格取值过于保守时，与实验结果的对比可知，周向划分25 个网格点，对计算结果几乎无影响。但根据法向网格与周向网格的比可知，此时要求周向50 个网格点。因此，若确定的法向网格尺度过于保守，必然导致给出的周向网格尺度也偏于保守，从而致使周向网格点数过多。

图2 Ma∞=8 时周向网格尺度对斯坦顿数的影响Fig.2 Influence of circumferential grid scale on Stanton number at Ma∞=8

3.1.2Ma∞=16.34 半圆柱绕流

为进一步验证第2 节的网格划分方法，以及3.1.1 节的初步结论。本节选择文献［24］给出的实验结果进行验证。圆柱直径D=0.076 m，自由来流的马赫数、静温、静压分别为Ma∞=16.34、T∞=52 K、P∞=82.95 Pa；壁 面 温 度Twall=294.4 K，基于圆柱直径的来流雷诺数为Re=2.97×105。 根 据 第2 节 确 定 的Δnw=1.77×10-6m、周向网格尺度Δx=9.61×10-4m，对应的周向节点数约为125。与3.1.1 节类似，当给出的法向网格尺度过小时，周向网格尺度亦过小。

基于3.1.1 节给出的方法，应用数值模拟方法，根据图3 可知数值计算确定合理的法向网格尺度Δnw=2.0×10-5m，可得周向网格尺度上限为Δx=1×10-2m，对应的周向网格节点数为11。由图4 亦可得，当法向网格取值合理时，用本文的方法给出的流向网格的取值亦合适。

图3 Ma∞=16.34 时法向网格尺度对斯坦顿数的影响Fig.3 Influence of normal grid scale over wall on Stanton number at Ma∞=16.34

图4 Ma∞=16.34 时周向网格尺度对斯坦顿数的影响Fig.4 Influence of circumferential grid scale on Stanton number at Ma∞=16.34

3.2 三维球头绕流

3.2.1Ma∞=8 球头绕流

为验证第2 节的网格划分方法以及3.1.1 节的初步结论，本节选用文献［17］给出的三维球头算例进行数值验证。球头直径D=15 mm，自由来流参数为Ma∞=8，T0=892 K，P0=7.8 MPa，Twall=294 K，基于球头直径的特征雷诺数为Re=1.7×105。根据第2 节确定的Δnw=7.5×10-7m，对 应 的 切 向 网 格 尺 度 为hτ=3.09×10-4m，对应的节点数约为38。与3.1.1 节类似，即当给出的法向网格尺度过小时，周向网格尺度亦过小。

采用与3.1 节相同方法，应用数值模拟方法确定法向网格尺度上限。根据图5 可知，数值计算确定合理的法向网格尺度Δnw=4.3×10-6m，对应的切向网格节点数为9。由图6 可得，当法向网格尺度取值合理时，用本文的方法给出的切向网格尺度的取值亦合适。

图5 Ma∞=8 时法向网格尺度对斯坦顿数的影响（球头）Fig.5 Influence of normal grid scale over wall on Stanton number at Ma∞=8 （sphere）

图6 Ma∞=8 时切向网格尺度对斯坦顿数的影响Fig.6 Influence of circumferential grid scale on Stanton number at Ma∞=8

图7 给出了球头表面斯坦顿数分布云图，由图可见，数值方法较好地模拟了球头表面的斯坦顿数分布。

图7 Ma∞=8 球头壁面斯坦顿数分布Fig.7 Stanton number distribution over sphere at Ma∞=8

3.2.2Ma∞=13 球头绕流

为验证第2 节的网格划分方法在来流马赫数较高的球头算例适用性，本节选用文献［25］中欧洲风洞F4 试验的标准球头进行了数值计算，其直径D=0.35 m，试验来流参数为U∞=4 320 m/s，T∞=265 K，P∞=52.81 Pa，Twall=294.4 K，可得基于球头直径的特征雷诺数为Re=6.26×104。由第2 节的方法则法向网格尺度为Δnw=4.23×10-6m，切向网格与法向网格尺度比约为250，可得切向网格尺度为hτ=1.06×10-3m。

由前述可知，该尺度偏于保守，采用同前两节的方法，采用数值模拟方法确定网格尺度上限。由图8 可知，数值计算确定的Δnw=3.9×10-5m，可得切网格尺度hτ=9.75×10-3m，所对应的切向网格节点数为28。由图9 可得，当法向网格尺度取值合理时，用本文的方法给出的切向网格尺度的取值亦合适。

图8 Ma∞=13 时法向网格尺度对斯坦顿数的影响Fig.8 Influence of normal grid scale over wall on Stanton number at Ma∞=13

图9 Ma∞=13 时切向网格尺度对斯坦顿数的影响Fig.9 Influence of circumferential grid scale on Stanton number at Ma∞=13

图10 给出了球头表面斯坦顿数分布云图。由图可见，数值方法较好地模拟了球头表面的斯坦顿数分布。

图10 Ma∞=13 球头壁面斯坦顿数分布Fig.10 Stanton number distribution over sphere at Ma∞=13

3.3 文献算例验证

为进一步验证高超声速流动中误差匹配原则的适用性，本节引用已发表文献中相关数值计算算例，验证复杂外形高超声速数值计算中误差匹配原则的适用性。

3.3.1Ma∞=10 钝锥绕流

本节通过文献［18］中钝头圆锥外形验证第2 节中给出的网格划分方法的适用性，其流场的数值模拟通过隐式TVD 格式和LU 分解算法进行求解三维N-S 方程，钝锥由直径D=0.055 88 m 的球头和半锥角为15°的锥体组成，总长为0.447 04 m，几何外形可见文献［18］。其来流参数为Ma∞=10.6，T∞=47.3 K，P∞=109.9 Pa，Twall=294.4 K，攻角为0°。可得基于球头直径的特征雷诺数为Re=2.2×105，由第2 节的方法，可得法向网格尺度为3.18×10-6m，切向网格与法向网格尺度比约为460，切向网格尺度为hτ=1.46×10-3m。文献［18］给出的经网格无关性验证的壁面法向第1层网格尺度为Δnw=4.2×10-6m，切向网格尺度hτ＜1.397×10-3m，与应用第2 节方法给出的切向网格尺度与法向网格尺度比相近。

3.3.2Ma∞=15 升力体绕流

本节通过文献［26］中升力体外形验证第2 节给出的网格划分方法的适用性，采用有限体积法对控制方程进行离散，数值格式为van Leer，应用3 阶精度的MUSCL 方法使数值精度为3 阶，时间推进采用隐式LU-SGS 方法。其升力体头部直径D=140 mm，模型总长为3.3 m，几何外形见文献［26］。来流参数为Ma∞=15，T∞=250.35 K，P∞=287.4 Pa，Twall=300 K；基于升力体头部直径的自由来流Re=1.68×105。由第2节的方法，可得法向网格尺度为6.2×10-7m，切向网格与法向网格尺度比约为410、切向网格尺度为hτ=2.54×10-4m。文献［26］给出的经网格无关性验证的壁面法向第1 层网格尺度为Δnw=8.4×10-6m，切向网格尺度hτ＜4.1×10-3，与应用第2节方法给出的切向网格尺度与法向网格尺度比一致。

4 讨论及分析

4.1 表面热流收敛性

高超声速数值模拟中，表面热流密度的收敛性较慢。本文采用的表面热流收敛要求是，首先计算至残差不再下降，并稳定在一定的数值范围后，再取2 个相邻5 000 步的热流数值不再变化作为收敛判断准则。3.1 节中Ma∞=8，16.34 圆柱绕流与3.2 节中Ma∞=8，13 球头绕流算例残差历程如图11 所示。

图11 数值计算残差历程Fig.11 Residual history for different numerical simulations

以上方法确保了本文高超声速热流计算的收敛性。

4.2 数值格式对表面热流的影响

表1 给出了本文及文献［17-18， 26］计算钝头体表面热流的网格划分尺度。由表1 可见，切向网格尺度对高超声速钝头体表面热流的影响不如法向网格影响显著，并且只要钝头体壁面法向网格尺度合适，根据第2 节给出的切向网格划分尺度亦合理。

表1 不同算例网格划分尺度汇总Table 1 Collections of grid scales for different numerical cases

表1 中文献［17］及文献［26］虽然采用了改进了熵修正函数的二阶TVD 格式，但是数值模拟给出的壁面网格雷诺数相比其他算例还是过小。而文献［6］由于采用未改进二阶TVD 格式，为了能准确计算表面热流，所要求的法向网格尺度为壁面分子自由程的大小［1］（此时基于壁面参数的网格雷诺数为1.5）。其原因在于，上述2 个文献采用的二阶精度TVD 格式在热流计算上由于格式的高耗散型而存在困难［27］，并且文献［28］发现即使物面网格很小，原TVD 格式计算驻点热流仍较实验和文献给出的结果偏低。而文献［17， 26］由于采用了焓守恒修正的van Leer 格式数值离散，并用MUSCL 插值到3 阶精度，从而有效抑制了数值格式的黏性效应，所以对壁面网格雷诺数放宽了要求。文献［11］研究结果表明，当采用4 种差分格式，并用MUSCL 插值到二阶精度时，Roe 格式的壁面雷诺数可放宽到约50（此处根据该文献基于来流参数的网格雷诺数，换算为基于壁面参数的网格雷诺数），而SLAM格式由于有效抑制了数值格式的黏性，壁面雷诺数甚至可放宽到约90。文献［10］的研究结果则表明，不同迎风型数值格式对网格雷诺数的要求不同，应选择数值黏性小，则可放宽对网格雷诺数的要求。

4.3 壁温及驻点热流密度的影响

由于高超声速飞行器表面热流影响参量较多，以下研究圆柱及球头不同表面温度，以进一步验证本文结论的适应性。图12 给出了不同外形、不同壁面温度时，法向网格雷诺数对斯坦顿数的影响。由图可见，壁面温度不同时，网格雷诺数上限范围几乎不变，即根据本文提出的方法确定的网格尺度对壁面温度不敏感。

图12 不同壁温下法向网格尺度对表面斯坦顿数的影响Fig.12 Influence of normal grid scale on Stanton number under different skin temperature conditions

表2 给出了本文及参考文献［18， 26］中算例流场参数及驻点热流密度等数值，其中Tt为总温。由表2 可见，本文得出的结论对不同来流总温与表面温度比TtTw、变化范围较大（约50 倍）的驻点热流密度qs均具有一定的适应性。

表2 不同算例流场参数汇总Table 2 Collections of flow field parameters for different numerical cases

综上，对采用数值模拟进行高超声速飞行器表面热流预估时，若所选用的数值格式未进行数值黏性影响研究，那么在数值模拟的具体实施中需进行网格无关性验证。对于数值黏性小的格式，可取的壁面网格雷诺数分别为160、80、40 进行网格无关性验证；若采用数值黏性较大的格式，则需要采用40、20、10 的壁面网格雷诺数进行验证。而流向网格的划分，则可采用3 套网格中间的壁面网格雷诺数给定的壁面法向第1 层网格尺度，应用流向和法向网格尺度比的方法，确定流向网格尺度，并采用等间距的流向网格划分。若所选择的数值格式黏性较小，壁面法向网格雷诺数可增大到40 以上。

5 结 论

应用理论分析及数值模拟，提出了一种基于壁面网格雷诺数及基于钝头体特征长度的来流雷诺数的网格划分方法，并给出了壁面网格雷诺数的取值范围。应用所提出的网格划分方法对不同来流条件的高超声速半圆柱及球头表面热流进行了数值模拟，结论如下：

1）法向方向网格尺度对高超声速钝头体表面热流的影响大于切向网格尺度的影响。

2）基于能量方程的修正方程各误差项匹配分析表明，为使各误差项为同一量级从而减小数值误差，流向与法向第1 层网格尺度比应为基于钝体特征长度的来流雷诺数的平方根。

3）基于壁面网格雷诺数确定的法向网格尺度的取值受数值格式影响，数值黏性小的格式可取的壁面网格雷诺数大。对于工程上采用的能有效抑制数值黏性的格式，壁面网格雷诺数可供参考的取值范围为20～80。

应该注意，本文给出的方法受所选择的数值格式的耗散性影响较大，在应用本文给出的网格划分方法前，应选择经评估的、数值耗散性小的数值格式。对数值黏性大的格式则不能进一步放宽网格尺度的要求。

另外，虽然本文给出了一种计算高超声速钝体表面热流的网格划分方法，但高超声速真实飞行器表面热流的准确计算，受数值方法、网格质量、边界条件等多种因素的影响，因此真实高超声速飞行器表面热流的精确预估还需要深入研究，以满足工程应用的需求。