一种考虑滑移修正的辐射平衡温度计算方法

摘要

本发明公开了一种考虑滑移修正的辐射平衡温度计算方法，主要用于在基于Navier‑Stokes控制方程组数值模拟高超声速滑移流动过程中，滑移温度边界和壁面辐射平衡温度边界的耦合处理过程。该方法在传统辐射平衡能量方程的基础上，引入壁面气体稀薄滑移温度修正，形成基于滑移温度的辐射平衡方程，通过隐式迭代求解滑移温度辐射平衡方程得到气体滑移温度，再通过表面温度滑移关系得到壁面辐射平衡温度分布。该方法较好地兼容了表面辐射能量平衡和表面温度稀薄滑移效应两种物理机制，计算稳定性好，收敛迅速。

权利要求

1.一种考虑滑移修正的辐射平衡温度计算方法，其特征在于在传统能量辐射平衡方程的基础上，引入壁面气体稀薄滑移温度修正，得到基于滑移温度的能量辐射平衡方程，通过隐式迭代求解基于滑移温度的能量辐射平衡方程得到气体滑移温度，再通过气体温度滑移关系得到壁面辐射平衡温度分布，包括以下步骤：S1：在高超声速滑移流动数值模拟过程中，获取壁面网格微元上壁面温度初值 、气体滑移温度初值 和对应壁面法向次层微元的气体温度为 ；

S2：网格微元上温度滑移方程的余量为： ，表面能量辐射平衡方程的余量为： ，当网格微元满足  时，需要进行迭代计算，在迭代次数 时，；

当网格微元上满足 ，不需要进行迭代计算，即气体滑移温度，壁面辐射平衡温度 ，

S3:对于需要进行迭代计算的流场网格微元，在传统能量辐射平衡方程的基础上，引入壁面气体稀薄滑移温度修正，得到基于滑移温度的能量辐射平衡方程，隐式迭代计算滑移温度：

S4：根据网格微元上迭代残差与迭代计算精度要求 的大小比值，判断第 步迭代是否收敛；

S5：反复迭代S3和S4，直至计算收敛，得到气体滑移温度 ，再由壁面气体温度滑移关系，计算得到壁面辐射平衡温度 的分布；

其中：为滑移修正函数， 为迭代计算精度要求，为气体热传导系数， 为与壁面温度无直接关系的热流项，是表面辐射系数，是斯蒂芬‑玻耳兹曼常数,  为第一层壁面网格的法向间距， 为迭代次数， 和 分别为迭代的第 和 步的滑移温度。

2.根据权利要求1所述的一种考虑滑移修正的辐射平衡温度计算方法，其特征在于S1所述的高超声速滑移流动数值模拟，其流动控制方程为考虑滑移效应的Navier‑Stokes方程组。

3.根据权利要求1所述的一种考虑滑移修正的辐射平衡温度计算方法，其特征在于： S3中的 为与壁面温度无直接关系的热流项，由实际的表面物理机制决定，包括表面组分扩散热流、壁面质量引射能量通量、壁面机械流失能量通量中的单项或多项组合。

4.根据权利要求1所述的一种考虑滑移修正的辐射平衡温度计算方法，其特征在于在S4中，当迭代残差满足 时，该微元迭代计算收敛；当时，该微元迭代计算未收敛，进行下一步迭代，即 ，返回步骤S3。

5.根据权利要求1所述的一种考虑滑移修正的辐射平衡温度计算方法，其特征在于：所述滑移修正函数 的表达式由数值模拟采用的滑移模型决定，忽略热蠕动效应，将滑移模型中一阶温度滑移条件通过 的形式表达，所述 为气体温度壁面法向梯度。

说明书

一种考虑滑移修正的辐射平衡温度计算方法

技术领域

[0001] 本发明涉及数值模拟计算领域，具体涉及高超声速稀薄滑移流动壁面辐射平衡温度的计算方法。

背景技术

[0002] 壁面辐射平衡温度条件和表面气体温度滑移条件，是高超声速流动模拟中常用的两种温度边界条件。针对这两种边界条件，单独的研究较为常见：在使用壁面辐射平衡温度
条件时，常忽略表面气体的稀薄滑移效应，即认为表面气体温度与壁面温度相同；而在考虑
表面气体温度滑移效应时，则一般采用等温壁面条件，即壁面温度为某一固定的温度条件。
二者的综合运用，较为少见。

[0003] 随着飞行速度与高度的提高，在高空高焓条件下，飞行器壁面温度往往接近于辐射平衡温度，同时飞行器表面气体存在温度滑移现象（一般认为飞行高度60km以上出现）。
精细模拟高空高焓高超声速流动时，需要同时考虑这两种物理机制。

[0004] 一般来说，壁面辐射平衡温度可由迭代求解壁面温度的能量辐射平衡方程得到，而表面气体滑移温度则由壁面温度和表面气体温度梯度，通过气体滑移模型给出，因此二
者的计算顺序通常为先由能量辐射平衡方程计算得到壁面温度，再由壁面温度通过滑移模
型计算得到气体滑移温度。

[0005] 这种方法对于高超声速稀薄滑移流动与壁面辐射平衡温度的解耦计算是适用的，但在二者耦合的数值模拟过程中，很容易出现非物理的振荡、不收敛甚至发散。这主要是由
于在流动控制N‑S方程每一步求解过程中，流场中非边界点微元的值是由N‑S方程计算得到
的，因此在壁面温度边界条件处理时，壁面次层网格微元的温度是一个已知的确定值。此
时，在壁面温度边界处理的能量辐射平衡方程求解过程中，如果由于某种原因（如数值误差
等）导致壁面温度出现非物理的波动，这一波动将通过滑移关系传导至气体滑移温度，从而
使滑移温度出现波动；又由于“壁面次层微元的温度是一个的确定值且第一层壁面法向距
离极小”，这将导致壁面热流的剧烈波动；而热流的剧烈波动有反过来影响能量辐射平衡方
程的求解，从而使非物理波动迅速反馈放大，甚至导致计算失效、发散。

[0006] 因此，仍有必要发展较为高效、稳定的考虑滑移修正的壁面辐射平衡温度计算方法。

发明内容

[0007] 本发明的目的是先通过隐式迭代求解滑移温度辐射平衡方程得到气体滑移温度，避免传统方法的“误差迭代反馈放大”的问题，再通过表面温度滑移关系，反向计算得到壁
面辐射平衡温度分布，兼容表面辐射能量平衡和表面温度稀薄滑移效应两种物理机制。

[0008] 为了实现上述目的，本发明采用如下技术方案：

[0009] 步骤一：从数值模拟过程中，获取壁面网格微元温度初值、气体滑移温度初值、壁面次层微元温度及其相关量；

[0010] 步骤二：判断需要进行迭代计算的流场区域；

[0011] 步骤三：对于需要进行迭代的流场区域，构建基于滑移温度的辐射平衡方程，通过隐式迭代计算滑移温度；

[0012] 步骤四：判断迭代是否收敛；

[0013] 步骤五：重复步骤三、四，直至迭代收敛，得到气体滑移温度，再由壁面气体温度滑移关系，计算得到壁面温度。

[0014] 本发明主要用于高超声速滑移流动数值模拟过程中壁面温度边界的处理，这里流动数值模拟控制方程为考虑滑移效应的Navier‑Stokes方程组。

[0015] 综上所述，由于采用了上述技术方案，本发明的有益效果是：

[0016] 本发明较好地兼容了表面辐射能量平衡和表面温度稀薄滑移效应两种物理机制, 实现了滑移温度边界和壁面辐射平衡温度边界耦合处理。

[0017] 本发明在传统辐射平衡能量方程的基础上，引入壁面气体稀薄滑移温度修正，形成基于滑移温度的辐射平衡方程，由于该方程综合表征了滑移温度、表面热流、表面辐射之
间的能量平衡关系，因此避免了传统方法的“误差迭代反馈放大”的问题。

[0018] 本发明对辐射能量进行隐式处理，形成滑移温度辐射平衡方程隐式迭代计算方法，避免了“随温度变化，辐射能量变化剧烈”带来的迭代稳定性差的问题。

附图说明

[0019] 本发明将通过例子并参照附图的方式说明，其中：

[0020] 图1为本方案的计算流程示意图；

[0021] 图2为实例一返回器外形示意图；

[0022] 图3为实例一采用本发明计算得到的热流与文献比较；

[0023] 图4为实例一采用本发明计算得到的壁面辐射平衡温度和气体滑移温度。

具体实施方式

[0024] 本说明书中公开的所有特征，或公开的所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以以任何方式组合。

[0025] 本说明书（包括任何附加权利要求、摘要和附图）中公开的任一特征，除非特别叙述，均可被其他等效或具有类似目的的替代特征加以替换。即，除非特别叙述，每个特征只
是一系列等效或类似特征中的一个例子而已。

[0026] 如图1 所示，是本实施例的计算流程，以“考虑Maxwell经典滑移模型的完全气体稀薄滑移流动辐射平衡边界处理”为实例，说明其具体的实现过程：

[0027] 步骤S1：在高超声速滑移流动数值模拟过程中，获取壁面网格微元上壁面温度初值 、气体滑移温度初值  和对应壁面法向次层微元的气体温度为 ；

[0028] 这里主要针对控制方程为考虑滑移效应的Navier‑Stokes方程组的高超声速流动数值模拟。在控制方程每一个时间步的推进过程开始时，壁面网格微元上的壁面温度初值
和气体滑移温度初值 ，均为已知条件(由上一步的推进给出)；经控制方程离散计算
之后，就能得到壁面法向次层微元的气体温度为 。在这一过程，还需得到符合“表面辐射
能量平衡”和“表面温度稀薄滑移效应”两种物理机制的壁面温度边界，即壁面温度 和气
体滑移温度  分布，以用于流动控制方程组下一步推进。因此，本实施例的计算，主要用于：
由  、  、  等已知条件，通过迭代计算，得到同时满足上述两种物理机制的气体滑移
温度 和壁面温度  。

[0029] 步骤S2：根据网格微元上温度滑移方程与表面能量辐射平衡方程的余量和与迭代计算精度要求的大小比值，判定是否需要进行迭代计算；

[0030] 当网格微元上满足 时，则微元上的 、 、 等已经达到了同时满足温度滑移关系与表面能量辐射平衡方程的计算精度要求，即同时满足了“表面
辐射能量平衡”和“表面温度稀薄滑移效应”两种物理机制，因此不需要进行迭代计算，气体
滑移温度  ，壁面温度  。反之，当网格微元上满足
时，则需要进行迭代计算，设定 为迭代计算的初值，即在迭代次数  时，  。
这里R1为温度滑移方程的余量， R2为表面能量辐射平衡方程的余量。

[0031] 之所以在步骤S2中判断是否需要迭代，是因为在控制方程推进过程中，网格微元上的 、  等均由流动控制方程上一个时间步推进的迭代计算得到，因此可能直接就满
足了温度滑移方程与表面能量辐射平衡方程，无需进入迭代计算过程，可以减少计算量。随
着控制方程的不断推进，流场逐渐趋于稳定，这种“无需计算”的区域，还可能增大，从而使
计算量进一步减小。

[0032] 其中温度滑移方程的余量 ，在本实例中，可由Maxwell经典滑移模型得到。忽略热蠕动效应，Maxwell经典滑移模型中滑移温度条件的一阶形式为：

[0033]

[0034] 其中 为气体滑移温度，  为气体温度壁面法向梯度，  为壁面温度，  为气体比热比， Pr是普朗特数，  为热协调系数。  为当地气体分子自由程，其形式一般为：

[0035]

[0036] 这里  分别为气体粘性系数、气体密度和气体普适常数。

[0037] 由上述两个公式可得到温度滑移方程:

[0038]

[0039] 这里为滑移修正函数，根据Maxwell经典滑移模型，其形式为：

[0040]

[0041] 通常情况下，高超声速流场壁面压力  满足  条件，因此壁面压力主要受流场结构影响，受壁面气体温度影响相对较小。在本实例中，由完全气体状态方程，

[0042]

[0043] 可知在壁面处，气体密度受壁面气体温度影响相对较大。因此，为了更好地表征壁面气体滑移温度对M的影响，这里将上式代入，得到M关于p、Ts等变量的函数形式，

[0044]

[0045] 将温度滑移方程进行一阶离散，并代入 ，得到其方程余量

[0046]

[0047] 其中 为第一层壁面网格的法向间距。可见，当 趋近于0时，  近似满足温度滑移方程。

[0048] 表面能量辐射平衡方程余量 ，可由表面微元辐射能量平衡方程得到。在壁面处构建无限薄的贴体微元，根据能量守恒关系，得到辐射能量平衡方程

[0049]

[0050] 其中  为气体热传导系数，  是表面辐射系数，  是斯蒂芬‑玻耳兹曼常数， 为壁面向外的辐射能量。 为与壁面温度无直接关系的热流项，由实际的表面物理机
制决定，包括但不限于表面组分扩散热流、壁面质量引射能量通量、壁面机械流失能量通量
中的单项或多项，消除上式中的  ，则辐射能量平衡方程可写为不含微分的形式：

[0051]

[0052] 将 、  代入上式，得到其方程余量

[0053]

[0054] 可见，当  趋近于0时，  、  等近似满足辐射平衡能量方程。

[0055] 步骤S3:对于需要进行迭代计算的流场网格微元，在传统能量辐射平衡方程的基础上，引入壁面气体稀薄滑移温度修正，得到基于滑移温度的能量辐射平衡方程，隐式迭代
计算滑移温度：

[0056]

[0057] 为迭代步数，  和  分别为迭代的第 和  步的滑移温度。

[0058] 这里对迭代关系式的构建方法进行介绍。由消去辐射能量平衡方程中的壁面温度 ，得到关于滑移温度的辐射平衡方程

[0059]

[0060] 将上式的右端项的辐射能量记为

[0061]

[0062] 其一阶导数形式为

[0063]

[0064] 滑移温度的辐射平衡方程可写为：

[0065]

[0066] 对上式中传导热流和辐射能量项进行迭代隐式处理，得到隐式方程形式：

[0067]

[0068] 将辐射能量项  在  处进行Taylor展开

[0069]

[0070] 这里 ，如果忽略高阶无穷小量  ，则上式改写为：

[0071]

[0072] 将上式代入合并，可得：

[0073]

[0074] 将上式代入隐式方程形式，并采用空间一阶离散得

[0075]

[0076] 将上式进一步整理，得到隐式迭代计算滑移温度的公式。

[0077] 由于隐式迭代计算滑移温度的公式的迭代综合表征了滑移温度、表面热流和表面辐射之间的能量平衡关系，因此避免了传统方法的“误差迭代反馈放大”的问题。同时，隐式
迭代计算滑移温度的公式将辐射能量进行隐式处理，形成滑移温度辐射平衡方程隐式迭代
计算方法，因此避免了“随温度变化，辐射能量呈‘温度的4次方’剧烈变化”带来的迭代稳定
性差的问题。

[0078] 步骤S4：根据网格微元上迭代残差与迭代计算精度要求的大小比值，判断第m+1步迭代是否收敛；

[0079] 迭代残差满足  时，该微元迭代计算收敛；

[0080] 当  时，该微元迭代计算未收敛，进行下一步迭代，即，返回步骤S3。

[0081] 步骤S5：反复迭代S3和S4，直至计算收敛，得到气体滑移温度，再由壁面气体温度滑移关系 ，计算得到壁面辐射平衡温度的分布。

[0082] 这里壁面气体温度滑移关系，由温度滑移方程空间一阶离散得到。

[0083] 实施例一：返回器非平衡稀薄滑移流动与表面辐射平衡耦合数值模拟。返回器外形如图2所示，计算模拟状态为飞行高度85km，来流速度7600.0m/s。高温气体模型为5组分
空气化学反应Park模型和单温度模型，滑移模型为经典的Maxwell滑移模型；飞行器表面完
全催化条件(FCW)和完全非催化条件（NCW）。

[0084] 图3给出了实例一采用本发明的计算结果与文献（Votta R, Schettino A, Ranuzzi G, Borrelli S, Hypersonic low‑density aerothermo‑dynamics of orion‑
like exploration vehicle[J]. Journal of Spacecraft and Rockets, 2009, 4(46):
781 787）的比较。可以看出，其热流结果与DSMC结果更接近。
~

[0085] 图4给出了实例一采用本发明计算得到的壁面辐射平衡温度和壁面气体滑移温度分布。可以看出，本发明计算得到的温度分布光滑平整，没有出现非物理的数值振荡。这说
明本发明能稳定地实现“非平衡稀薄滑移流动”与“表面辐射平衡”两种物理机制耦合的数
值模拟。由图还可以看出，文献采用等壁面温度（1184K）的设定，只能在一定程度上近似模
拟壁面热平衡时的状态。

[0086] 本发明并不局限于前述的具体实施方式。本发明扩展到任何在本说明书中披露的新特征或任何新的组合，以及披露的任一新的方法或过程的步骤或任何新的组合。