高超声速飞行器热环境计算方法、装置和计算机设备

1.本发明涉及高超声速空气动力学、计算流体力学技术领域，具体是一种高超声速飞行器热环境计算方法、装置和计算机设备。


背景技术：

2.高超声速飞行器具有速度快、突防能力强、射程远，快速响应、精确打击等特点，是掌握战略制高点的有效工具，代表了航空航天领域未来的发展方向，在多个领域中均发挥着重要的战略作用，是当今竞相发展的前沿技术。然而，高超声速飞行器的研制仍然面临诸多挑战，其中高超声速飞行器热安全问题是制约高超声速飞行器技术发展的重要难题之一。提升新型临近空间高超声速飞行器的性能需要发展更高效、更高精度的气动热环境预测方法。
3.随着物理模型、数值算法和计算机软硬件水平的不断发展，计算流体力学（cfd）被广泛地用于解决各种航空航天飞行器设计中涉及的大量气动力热问题。数值算法作为cfd的核心，是影响计算精度的重要因素之一，一直以来都是研究的重点。然而，由于高超声速流动具有高马赫数、大能量的特征，高超声速流场中常常呈现出强激波和黏性干扰等现象，同时伴随有高强度的激波/激波干扰、激波/边界层干扰等复杂流场结构。这就要求数值格式在精确分辨边界层等黏性流动的同时，又能够实现对激波等强间断的稳定捕捉，有效消除间断附近的非物理振荡。
4.目前对于气动热环境预测，大多采用的是传统的激波捕捉格式，激波捕捉格式分为低耗散格式和高耗散格式。低耗散格式如roe格式、hllc格式，由于激波稳定性不足，导致产生非物理解和数值痈现象，使得飞行器表面热流计算结果紊乱。高耗散格式，如hll格式、vanleer格式，耗散性性大，可以抑制激波不稳定，但是对边界层流动的速度型具有抹平作用，导致计算出的飞行器表面热流严重偏小。
5.在高超声速流动的数值模拟中，激波捕捉的质量与稳定性直接影响到气动特性特别是气动热的预测精度与可靠性。目前，应用于高超声速流动数值模拟的激波捕捉方法通常会遭遇激波不稳定问题，如著名的“carbuncle”与激波后振荡等异常现象。不仅影响了数值计算效率，而且严重降低气动热预测的可靠性。值得关注的是，数值试验表明无论是一阶的迎风通量格式还是二阶高分辨率格式、高阶格式，在模拟多维强激波时均呈现出明显的激波不稳定现象。并且，数值格式的激波稳定性对计算网格具有明显的依赖性。这些因素严重限制了激波捕捉格式特别是高分辨率、高精度格式在复杂外形高超声速气动力热预测中的应用。
6.目前克服激波不稳定现象的途径主要是混合方法，即在激波处使用耗散大的格式，在光滑流动区域使用耗散低的格式，混合格式本质上是利用耗散大的格式在激波处弥补多维数值黏性的不足，以实现抑制激波不稳定的目的。然而，混合类格式的求解精度受该转换函数的影响，尤其是在激波与黏性流动相干扰的复杂流动情况下。此外，由于此类方法引入的数值黏性通常基于“准一维”的黎曼通量，因此提升激波稳定性的效果仍然对计算网
格具有很强依赖性。此外，即便是处理器hlle、处理器ausm 处理器等鲁棒性高的格式在不同网格形式下均表现出不同程度的不稳定问题。实际上，目前常用的“准一维”黎曼通量的数值黏性与网格特性有密切联系，是导致激波捕捉稳定性具有很强网格依赖性的根源。因此，众多改进黎曼通量在某些网格下具有良好的激波稳定性，但是在网格发生变化时，同样会出现激波失稳问题。


技术实现要素：

7.针对上述现有技术中临近空间飞行器热环境预测问题，本发明提供一种高超声速飞行器热环境计算方法、装置和计算机设备，解决激波捕捉格式所存在的激波不稳定问题，同时精确分辨边界层黏性流动，进而提高高超声速气动热的计算精度。
8.为实现上述目的，本发明提供一种高超声速飞行器热环境计算方法，包括如下步骤：步骤1，读取飞行器模型，生成流场计算网格，并计算每个网格单元的面权重系数；步骤2，离散流场的三维流动控制方程，得到网格单元的离散方程，并基于当前的流场变量计算离散方程中的粘性通量与无粘通量；步骤3，基于当前的流场变量计算每个网格单元的人工粘性系数，并基于面权重系数与人工粘性系数计算每个网格界面的人工粘性项；步骤4，将人工粘性项添加至对应网格单元的离散方程后求解离散方程，完成流场的更新并得到新的流场变量；步骤5，根据时间推进重复步骤2-4对离散方程进行迭代求解，直至流场收敛；步骤6，基于收敛后的流场变量计算高超声速飞行器表面的热流。
9.在其中一个实施例，步骤2具体包括：获取积分形式的可压缩三维流动控制方程，为：(1)式中，为守恒变量矢量，为网格单元，为网格面矢量，、为界面通量矢量；、、的表达式分别为：（2）
式中，为密度，为总能；为速度矢量，为方向速度，为方向速度，为方向速度；为静压，为控制体表面外法线单位向量，为法向量的分量，为法向量的y分量，为法向量的z分量；为速度矢量在法向量上的投影，即；h为总焓，、、分别为、、三个方向的正应力，、、、、、为各个方向的切应力，、、为能量项；采用基于非结构网格数据结构的格心型有限体积法对式(1)进行离散，因此对于网格单元，可以得到离散方程为：（3）式中，为当前单元的编号，为与单元相邻的面的编号；为网格单元的体积，为网格单元的守恒变量，为网格的相邻面数目，为网格面的无粘通量，为网格面的粘性通量，为网格面的面积；直接通过式(2)计算网格单元的粘性通量，并采用amhllc近似黎曼求解器求解无粘通量。
10.在其中一个实施例，步骤3中，所述基于当前的流场变量计算每个网格界面的人工粘性系数，具体为：（4）式中，为网格界面的人工粘性系数，、为该网格界面的尺度，、、分别是当地密度、流动速度和声速，为速度散度；网格界面的尺度，，其中，为该网格界面的左右相邻网格单元中的最大面对角线长度，为该网格界面的左右相邻网格单元中的最小边长。
11.在其中一个实施例，步骤3中，所述基于面权重系数与人工粘性系数计算每个网格单元的人工粘性项，具体为：
（5）（6）式中，f
av
为人工粘性项，为网格单元的面权重系数，为速度u在y方向的导数，为速度v在x方向的导数，为速度u在z方向的导数，为速度w在x方向的导数，为速度v在z方向的导数，为速度w在y方向的导数。
12.在其中一个实施例，步骤4中，所述将人工粘性项添加至对应网格单元的离散方程即将人工粘性项f
av
添加式(3)中，为：（7）求解式(7)即能得到新的流场变量，完成流场的更新。
13.在其中一个实施例，步骤1中，所述网格单元的面权重系数为：式中，为网格单元的面权重系数，为该网格单元的面法向矢量与该网格单元距离最近的壁面法向矢量的夹角。
14.在其中一个实施例，步骤6中，所述计算高超声速飞行器表面的热流，具体为：（8）式中，为高超声速飞行器表面的热流，为气体比热比，r为气体常数，为物理粘性系数，为普朗特数，为温度在壁面法向的梯度，t为温度，n为壁面法向量。
15.实现上述目的，本发明还提供一种高超声速飞行器热环境计算装置，采用上述的方法进行高超声速飞行器热环境计算。
16.实现上述目的，本发明还提供一种计算机设备，包括：至少一个处理器；以及与所述处理器连接的至少一个存储器、总线；其中，所述处理器、存储器通过所述总线完成相互间的通信；所述处理器用于调用所述存储器中的程序指令，以执行上述方法的部分或全部步骤。
17.本发明提供的一种高超声速飞行器热环境计算方法、装置和计算机设备，突破准一维通量格式对于网格质量的依赖性，对于任何类型的网格单元，如四面体、金字塔、三棱柱、六面体单元，以及质量较差的网格，如大长宽比和大扭曲度网格，都能很好的抑制激波不稳定，同时不损失边界层速度型的分辨精度。对于高超声速飞行器常见的激波-激波干扰、激波-边界层干扰等复杂流动，本发明提供的方法同样可以精确分辨，不会发生常规人工粘性方法的精度损失的问题。采用本发明的方法，可以有效解决高超声速飞行器热环境预测时的激波不稳定问题和边界层计算不精确问题，从而获得高超声速飞行器表面的精确热流分布。
附图说明
18.为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图示出的结构获得其他的附图。
19.图1为本发明实施例中高超声速飞行器热环境计算方法的流程图；图2为本发明实施例中网格单元夹角计算示意图；图3为本发明实施例中空天飞机外形及网格划分示意图；图4为本发明实施例中传统激波捕捉方法计算的空天飞机头部热流示意图，其中：(a)是传统激波捕捉方法模拟的结果示意图，(b)是本发明方法模拟得到的结果示意图；图5为本发明实施例中计算得到的空天飞机头部热流示意图，其中：(a)是传统激波捕捉方法计算的结果示意图，(b)是本发明方法计算得到的结果示意图；图6为本发明实施例中超声速飞行器热环境计算装置的模块图；
图7为本发明实施例中计算机设备的架构图。
20.本发明目的的实现、功能特点及优点将结合实施例，参照附图做进一步说明。
具体实施方式
21.下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
22.需要说明，本发明实施例中所有方向性指示（诸如上、下、左、右、前、后
……
）仅用于解释在某一特定姿态（如附图所示）下各部件之间的相对位置关系、运动情况等，如果该特定姿态发生改变时，则该方向性指示也相应地随之改变。
23.另外，在本发明中如涉及“第一”、“第二”等的描述仅用于描述目的，而不能理解为指示或暗示其相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”的特征可以明示或者隐含地包括至少一个该特征。在本发明的描述中，“多个”的含义是至少两个，例如两个，三个等，除非另有明确具体的限定。
24.在本发明中，除非另有明确的规定和限定，术语“连接”、“固定”等应做广义理解，例如，“固定”可以是固定连接，也可以是可拆卸连接，或成一体；可以是机械连接，也可以是电连接，还可以是物理连接或无线通信连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通或两个元件的相互作用关系，除非另有明确的限定。对于本领域的普通技术人员而言，可以根据具体情况理解上述术语在本发明中的具体含义。
25.另外，本发明各个实施例之间的技术方案可以相互结合，但是必须是以本领域普通技术人员能够实现为基础，当技术方案的结合出现相互矛盾或无法实现时应当认为这种技术方案的结合不存在，也不在本发明要求的保护范围之内。
26.如图1所示为本实施例公开的一种高超声速飞行器热环境计算方法，具体包括如下步骤1-步骤6。
27.步骤1，读取飞行器模型，生成流场计算网格，并计算每个网格单元的面权重系数，具体地：对于流场计算网格上的每一个网格单元，先计算其面法向矢量与该网格单元距离最近的壁面法向矢量的夹角，夹角计算示意图如图2所示。再将夹角的正弦值作为该网格单元的面权重系数，即面权重系数。
28.步骤2，离散流场的三维流动控制方程，得到网格单元的离散方程，并基于当前的流场变量计算离散方程中的粘性通量与无粘通量，其具体实施过程为：首先，获取积分形式的可压缩三维流动控制方程，为：（1）式(1)中，u为守恒变量矢量，为网格单元，为网格面矢量，、为界面
通量矢量。其中，、、的表达式分别为：（2）式(2)中，为密度，e为总能；为速度矢量，为方向速度，为方向速度，为方向速度；为静压，为控制体表面外法线单位向量，为法向量的分量，为法向量的分量，为法向量的分量；v为速度矢量在法向量上的投影，即；h为总焓，、、为各个方向的粘性正应力，、、、、、为各个方向的粘性切应力，、、为能量项；然后，采用基于非结构网格数据结构的格心型有限体积法对式(1)进行离散，因此对于网格单元，可以得到离散方程为：（3）式中，为当前单元的编号，为与单元相邻的面的编号；为网格单元的体积，为网格单元的守恒变量，为网格i的相邻面数目，为网格面的无粘通量，为网格面的粘性通量，为网格面的面积；在对粘性通量与无粘通量的计算过程中，粘性通量可以直接通过式(2)计算得到，无粘通量则可以采用amhllc近似黎曼求解器求解得到：
式中，和表示间断右侧和左侧的数值通量，和指黎曼解的两个波速，分别对应间断的左侧和右侧，s*为中间状态波速。为中间状态左右两侧的数值通量，其计算公式为：式中，为面左右两侧网格的状态量，为中间间断左右两侧的状态量，为间断左右两侧的速度矢量与法向矢量点积，为修正压力。
29.步骤3，基于当前的流场变量计算每个网格单元的人工粘性系数，并基于面权重系数与人工粘性系数计算每个网格界面的人工粘性项。
30.本实施例中将一种改进的人工粘性方法应用到可压缩粘性流动计算中，从而实现高超声速热环境的精确计算。人工粘性方法通过在激波位置添加ns方程的粘性来达到抑制无粘流动激波不稳定的目的。基于当前的流场变量计算每个网格界面的人工粘性系数，具体为：（4）式中，为网格界面的人工粘性系数，、为该网格界面的尺度，、、分别是当地密度、流动速度和声速，通过左右网格单元的平均或者线性插值得到，为速度散度，速度散度可通过高斯定理求出。
31.一般在存在激波的地方，流动的压缩性表现明显，因此当网格界面上的时，该网格界面无需添加人工粘性，即。在传统的人工粘性方法中，通常是直接将网格界面的尺度设置为，其中，为该网格界面的左右相邻网格单元的最大面对角线长度，其通过引入最大特征尺度，人工粘性方法可以在任意网格单元上抑制无粘流动的激波不稳定现象。
32.然而，对于可压缩粘性流动，近壁区域流动梯度大，需要在法向加密网格，尤其是为了计算气动热，第一层网格高度达到量级，网格长宽比非常大，表现出明显的各项异性。壁面法向极小的网格尺度与边界层的流动特性导致速度散度的量级增大，人工粘性添加的判据很容易满足，造成真实物理粘性的精度损失，同时流场收敛也很困难。因此本实施例中对尺度进行修正：即保持尺度不变，为；对于尺度，令，其中，为该网格单元的左右相邻网格单元中的最小边长。这一修正能够使得人工粘性在抑制激波不稳定的同时，不损失边界层流动的精度。
33.在已知同一网格单元的面权重系数与人工粘性系数后，即能计算得到该网格单元的人工粘性项，为：（5）式(5)类似于ns方程的粘性项，与之不同的是，本发明去除了粘性正应力项，只保留、、、、、这几项人工剪切粘性应力，其表达式为： （6）式(5)-(6)中，为人工粘性项，为网格单元的面权重系数，为速度u在y方向的导数，为速度v在x方向的导数，为速度u在z方向的导数，为速度w
在x方向的导数，为速度v在z方向的导数，为速度w在y方向的导数。其中，面权重系数对于气动热的精确计算至关重要，对于边界层内流动梯度大的网格面，即与固壁边界平行的网格面，=0，此时式(6)中的人工粘性应力项为0，克服了人工粘性方法对于边界层流动问题，尤其是激波-边界层干扰问题的精度不足的缺点。
34.步骤4，将人工粘性项添加至对应网格单元的离散方程后求解离散方程，完成流场的更新并得到新的流场变量，具体地：将人工粘性项f
av
添加式(3)中，为：（7）求解式(7)即能得到新的流场变量，完成流场的更新。
35.步骤5，基于新的流场变量再次计算粘性通量、无粘通量与人工粘性项后，再次代入式(7)求解，即根据时间推进重复步骤2-4对离散方程进行迭代求解，直至流场收敛。
36.步骤6，基于收敛后的流场变量计算高超声速飞行器表面的热流，具体地：基于流场的压力、密度、温度速度等原始变量以及变量的梯度计算高超声速飞行器表面的热流，为： ）（8）式中，为高超声速飞行器表面的热流，为气体比热比，r为气体常数，为物理粘性系数，为普朗特数，为温度在壁面法向的梯度，t为温度，n为壁面法向量。
37.下面结合具体的示例对本发明中的高超声速飞行器热环境计算作出进一步的说明。
38.以空天飞机算例为例，对所提出的高超声速飞行器热环境计算方法进行验证。空天飞机由机头、机身和机翼三部分组成，机头为修型钝头圆锥，机身由半圆柱（上半部）与半方柱（下半部）组成，机翼具有68
°
后掠角，以及双翼尖上反的特点。模型全长290mm，横向最大尺寸184.8mm，最大高度58mm。实验条件为来流马赫数，单位雷诺数，总压，总温，其计算网格如附图3所示，网格总量为356万，壁面第一层网格法向高度为m。
39.空天飞机对称面的压力云图如图4所示，上边的图4(a)是传统激波捕捉方法模拟的结果，飞行器头部弓形激波出现了非物理的数值不稳定现象，导致流场紊乱；下边的图4
(b)是本发明方法模拟得到的结果，头部弓形激波可以清晰、光滑地捕捉。
40.空天飞机头部的壁面热流分布云图如图5所示，上边的图5(a)是传统激波捕捉方法计算的结果，热流分布出现了明显的紊乱，驻点热流值与实验值偏差较大；下边的图5(b)是本发明方法计算得到的结果，热流分布规整，驻点热流值与实验值吻合良好。
41.本实施例还公开了一种高超声速飞行器热环境计算装置，该装置在计算高超声速飞行器热环境是采用上述高超声速飞行器热环境计算方法所执行的各个步骤。参考图6，该高超声速飞行器热环境计算装置包括：网格生成模块，用于读取飞行器模型与最初的流场变量，并生成流场计算网格；面权重计算模块，与网格生成模块相连，用于计算每个网格单元的面权重系数；粘性通量计算模块，与网格生成模块、流场更新模块相连，用于根据当前的流场变量计算离散方程中的粘性通量；无粘通量计算模块，与网格生成模块、流场更新模块相连，用于根据当前的流场变量计算离散方程中的无粘通量；粘性系数计算模块，与网格生成模块、流场更新模块相连，用于根据当前的流场变量计算每个网格面的人工粘性系数；人工粘性项计算模块，与流场粘性系数计算模块、面权重计算模块相连，用于根据面权重系数与人工粘性系数计算每个网格面的人工粘性项；流场更新模块，与粘性通量计算模块、无粘通量计算模块、人工粘性项计算模块，用于将粘性通量、无粘通量与人工粘性项代入离散方程求解，对流场进行更新；热环境计算模块，用于在流场收敛后根据流场变量计算高超声速飞行器表面的热流。
42.参考图7，本实施例还提供一种计算机设备，包括：发送器、接收器、存储器以及处理器。其中，发送器用于发送指令和数据，接收器用于接收指令和数据，存储器用于存储计算机执行指令，处理器用于执行存储器存储的计算机执行指令，以实现上述实施例中高超声速飞行器热环境计算方法所执行的各个步骤。其具体实施过程与前述高超声速飞行器热环境计算方法相同。
43.需要注意的是，上述存储器既可以是独立的，也可以跟处理器集成在一起。当存储器独立设置时，该终端设备设备还包括总线，用于连接存储器和处理器。
44.以上所述仅为本发明的优选实施例，并非因此限制本发明的专利范围，凡是在本发明的发明构思下，利用本发明说明书及附图内容所作的等效结构变换，或直接/间接运用在其他相关的技术领域均包括在本发明的专利保护范围内。