NNW-HyFLOW高超声速流动模拟软件框架设计

李鹏，陈坚强，丁明松，何先耀，赵钟，董维中,*

1. 中国空气动力研究与发展中心 空气动力学国家重点实验室，绵阳 621000

2. 中国空气动力研究与发展中心 计算空气动力研究所，绵阳 621000

在高超声速(马赫数Ma≥5)飞行条件下，飞行器绕流流场中将出现热化学非平衡效应(高温气体内能模态激发与松弛，以及组分间的复杂化学反应过程等[1])，对高超声速飞行器的气动力特性、气动热环境、流场结构、光辐射特性、等离子体环境、电磁散射特性和电磁通信等均有重要影响[2]。受运行成本和技术条件的制约，地面高焓测试设备可获取的热化学非平衡流场信息较为有限，难以实现真实飞行条件下对飞行器高温气体效应的系统性研究。因此，常采用数值模拟方法对高超声速飞行器的高温气体效应进行预测和评估。

自20世纪八九十年代以来，美国和欧盟等在高温气体非平衡效应研究方面做了大量的理论和应用研究，相关软件的发展从未中断，具备对高超声速飞行器高温非平衡流动及其相关的气动力、气动热环境以及等离子体分布等进行数值模拟和分析的强大功能，例如美国NASA兰利研究中心的LAURA[3]、美国Aerosoft公司的GASP[4]和GUST[5]、美国NASA兰利研究中心HAPB小组开发的 VULCAN[6]、美国NPARC联盟的WIND[7]以及法国ONERA的MSD[8]等。这些软件功能十分完备，应用面广，但在国内应用受限，难以获取技术支持。例如，GASP和GUST软件可以供北美的公司和大学免费使用，但是国内用户却很难得到它的PDF形式的使用手册。LAURA和VULCAN软件同样很难获取，国内多是通过公开发表的学术文献了解它们的相关动态。因此，国内许多大学和研究单位也在发展自己的高温非平衡流数值模拟软件，如北京航空航天大学的ACANS[9-10]、山东大学的PHAROS[1-2]、中国空气动力研究与发展中心的AEROPH_Flow[11-12]和HyperFLOW[13-14]等。这些In-house软件，在物理化学模型完整度、适应复杂流动状态、多种物理机制耦合以及多种网格适应性等方面，或多或少存在困难，很难得到有效推广；其代码低层次重复较多，可扩展性差，不利于软件的进一步开发。

为了建设功能完善、适用面广、可扩展性强、自由共享的国产高超声速流动数值模拟软件，在国家数值风洞(National Numerical Windtunnel，NNW)工程[15]的支持下，中国空气动力研究与发展中心计算空气动力研究所联合国内相关单位，基于风雷(PHengLEI)开源软件[16-17]基础框架开发以高温气体热化学非平衡效应及其气动力、气动热和气动物理特性模拟分析为主要功能，具有自主知识产权的新型高性能工业应用CFD软件，即国家数值风洞高超声速流动模拟软件(NNW-HyFLOW)。

区别于国内传统In-house软件，NNW-HyFLOW定位于国内共享和自主创新，其研制将面向更广的用户群体，包括但不限于高校、研究机构、相关工业部门和公司企业，满足高超声速流动和高温气体热化学非平衡效应及其相关的气动力、气动热、气体辐射以及等离子体环境等特性数值模拟和综合分析的需求。该软件将具备如下特点：① 集成国内外在化学反应模型、热力学模型、壁面催化、烧蚀、稀薄滑移、非平衡湍流模型、高效数值算法等方面的理论创新和先进方法；② 实现网格架构的多样性、功能的完整性、代码的可扩展性和执行高效性的协调统一。为了兼顾商业计算和科学研究，建立适用于结构、非结构、结构/非结构混合、重叠网格等多种网格类型的计算框架，满足用户对不同网格类型的偏好和选择；为了扩展复杂物理现象综合评估的功能需求，设计了多种物理效应清晰完备的接口定义，满足非平衡解算器与多种物理化学模型和气动物理等多学科软件耦合的计算需求；为了适应不同物理化学模型和气体属性等对象的多样性和高效计算需求，大量采用了面向对象类和高效通用的多维数组等数据结构设计，以减少代码的低效重复，增加程序模块的可扩展性，提升软件执行的计算效率。软件研发进度安排为：2017—2023年完成软件主体框架构建和常规高温气体非平衡流动模拟模块开发，形成软件1.0版本；2023—2027年将集成国内相关优势单位在高温非平衡效应研究方面的理论创新和先进方法，形成功能完善、适用面广、可扩展性强、高效的高超声速流动及其相关物理效应的综合模拟应用软件。目前已完成软件的主体框架和数据结构设计，形成软件测试版本。

本文将主要介绍NNW-HyFLOW的总体设计框架、特点和软件研发的相关进展。首先从框架设计、数据结构、耦合算法、并行计算以及接口设计等方面分别陈述软件的主要设计框架及特点。其次，对求解器解算模型及其理论方法进行简要介绍，包括模型方程、物理化学模型以及核心数值算法等。最后，选取典型的热化学非平衡流算例开展模拟测试，检验当前软件测试版本的正确性，以及对复杂工程外形热化学非平衡流场数值模拟的适用性和可靠性。

1 软件设计

NNW-HyFLOW软件基于风雷开源软件(风雷软件或PHengLEI)的基础框架进行设计开发，结合“黑匣子”原理，实现软件各模块的封装，以适应多种网格类型、多种物理机制、多种物理化学模型以及多种流动状态的模拟需求，因而具备适用性广、功能完整、扩展性好等优点。以下从框架设计、数据结构、耦合算法、并行计算以及接口设计等方面，详细介绍软件设计框架及特点。

1.1 框架设计

作为计算流体力学(CFD)软件，计算框架是NNW-HyFLOW软件设计的核心。如图1所示，NNW-HyFLOW的计算框架由前处理、求解器和后处理等3个系统模块构成。其中，前处理负责前台数据接收和转换，包括不同网格类型的格式转换、并行计算要求的网格分区转换以及化学模型参数读取与转换等。基于前置网格处理模块可将外部商业软件生成的网格格式，如Gridgen的grd、Pointwise的cgns、Fluent的cas以及GMSH的msh等，转换为风雷开源软件定义的混合格式，可满足不同类型网格的耦合计算需求。另外，前处理模块基于文献[17]所述的网格分区方法实现均衡的网格分区划分，能满足大规模并行计算的需要。与风雷开源软件的最大不同，NNW-HyFLOW前处理增加了化学模型参数读取与转换功能模块，可以满足对不同类型化学反应模型(包括自定义化学模型)的计算需求。

NNW-HyFLOW求解器包含黏性模型、化学模型、气体模型、壁面条件、数值格式以及控制条件等功能模块。为避免底层代码低效的重复利用，方便程序的进一步拓展与维护，其求解器继承了风雷开源软件面向对象的编程思想、多层次的数据结构及其在并行底层、边界处理和结果显示方面的优势，采用分层模块设计思想进行开发。如图2所示，各功能模块按照内核层、算法层、功能层和应用层进行分层构建，其中非平衡求解器基于已有的Navier-Stokes(N-S)求解器进行扩充和完善，在数值算法层面增加和封装多种物理化学模型及函数，在边界条件中增加多种非平衡流动特定的壁面条件，如辐射平衡温度、壁面催化、稀薄滑移、壁面材料烧蚀等条件类型，可以满足对多种物理机制、多种物理效应计算模拟的需要。为了兼容原有的完全气体CFD求解器，在数值算法层面、功能层面以及应用层面扩展相应模块，通过类的重载、继承和派生及其热化学流动空间、模板的定义，实现数据结构拓展和功能性扩充。

图2 功能模块分层设计示意图Fig.2 Hierarchical design of major functional modules

NNW-HyFLOW后处理主要用于气动力、热环境以及等离子体环境等非平衡流场特性的输出和显示。基于风雷开源软件的后置处理模块，NNW-HyFLOW增加了非平衡流动相关的结果处理函数，包括但不限于非平衡流场特性、气动力特性、气动热环境、等离子体分布特性等参数的提取和输出功能。同时增加了与气体辐射场、电磁场、固体温度场等模块耦合计算的接口模块，能够满足多种物理机制耦合效应模拟的需要。

NNW-HyFLOW对标LAURA软件进行开发，其在框架设计上几乎覆盖了LAURA软件的全部功能点，以LAURA[3]版本为例，两者具备的主要功能模块如表1所示。

表1 NNW-HyFLOW与LAURA的功能对比

1.2 数据结构

从类结构设计和数组两个方面分别介绍NNW-HyFLOW软件的底层数据结构。

NNW-HyFLOW采用面向对象编程思想设计和封装各功能模块，利用C++类的重载、继承和派生等多态特性实现拓展应用范围的目的。例如针对不同求解器采用了如图3所示的类设计，其中CFDSolver是所有基于CFD方法求解器的父类。为了方便程序的编写和维护，N-S方程和RANS方程采用解耦计算策略，因此由父类CFDSolver派生出基于N-S方程解算的N-S求解器(NSSolver)和基于RANS方程解算的湍流求解器(TurbSolver)。为了满足结构/非结构网格耦合计算的需要，再通过基类NSSolver和TurbSolver分别派生出基于结构网格解算的求解器(NSSolverStruct与TurbSolverStr)和基于非结构网格解算的求解器(NSSolverUnStruct与TurbSolverUnstr)。

图3 求解器类结构设计示意图Fig.3 Design diagram of solver class structure

NNW-HyFLOW采用完全耦合的策略将N-S 方程与组分输运/振动能量/电子能量方程联合解算，因此在CFDSolver基类下不再独立为非平衡求解器派生设计新类，非平衡求解功能仍基于NSSolver类进行开发，相关代码在NSSolverStruct和NSSolverUnStruct中扩充和添加。另外，在Solver基类下派生了气动物理求解器类(AeroPHSolver)、有限元求解器类(FESolver)等与气体辐射场、固体温度场等计算相关的其他求解器类型，可作为非平衡流场与气动物理、结构力学等解耦计算的预留接口。在CFDSolver基类下派生了磁流体求解器类(MHDSolver)，可以实现磁流体相关计算需求。

为适应完全气体、平衡气体以及化学非平衡和热化学非平衡气体等多种计算气体模型的一体化计算需求，NNW-HyFLOW采用了图4所示的类结构设计。在气体模型基类(Gas)下派生出完全气体(PerfectGas)和真实气体(RealGas)两个主要子类，其中完全气体又衍生出量热完全气体(CaloricallyPerfectGas)和热完全气体(ThermalPerfectGas)，原有的完全气体计算功能被量热完全气体类涵盖，而热完全气体类则作为接口用于实现平衡气体计算功能。另外，真实气体类主要用于热力学能量激发以及化学反应气体模型的计算，因此根据模拟大气环境和化学反应类型的差异，真实气体类下又派生了地球反应气体(EarthReactionGas)、火星反应气体(MarsReactionGas)以及燃烧反应气体(CombustionGas)等子类，分别用于模拟地球大气环境、火星大气环境以及各种燃气反应体系。

图4 气体模型类结构设计示意图Fig.4 Design diagram of gas model class structure

上述类结构设计为NNW-HyFLOW软件功能扩展提供了灵活方便的接口，充分满足了多种计算气体模型、多种物理化学效应一体化计算的需要。NNW-HyFLOW与LAURA软件在程序设计上的最大不同在于编程规范。采用C++面向对象编程思想使得NNW-HyFLOW在功能扩展方面更加灵活方便。LAURA软件的最新版本使用Fortran 95语言规范，同样也具备了面向对象的编程思想，能够设计派生数据类型以适应复杂数据结构构建的需要，但程序设计的主体仍然是面向过程，因此在底层代码复用和功能拓展方面，它并不如NNW-HyFLOW易用和方便。

另外在变量存储方面，NNW-HyFLOW仅对原有的数组空间进行拓展，对其维度不做变更，可保持原有数组数据结构的高效性。采用图5所示的数组拓展方式设计底层数据结构，兼顾了完全气体、平衡气体、单温度化学非平衡气体、两温度热化学非平衡气体、三温度热化学非平衡气体以及多振动温度模型等多种计算气体模型的一体化计算需要。相比于LAURA软件，计算功能更加全面，适用范围广。

图5 数组扩展方式示意图Fig.5 Diagram of array expansion method

1.3 耦合算法

基于图3所示的类定义和设计，NNW-Hy-FLOW利用Solver/CFDSolver等父类对象的虚函数定义及多态特性实现不同类型求解器在同一控制流程下的合作计算，满足结构/非结构网格耦合解算，以及非平衡求解器与其他求解器的耦合计算等要求。如图6所示，在NNW-HyFLOW框架下不同求解器耦合计算的控制流程包含3层循环：

图6 耦合计算控制流程示意图Fig.6 Diagram of coupled computing control process

1) 第1层循环：遍历每个迭代步，执行CFD计算的一次时间推进。

2) 第2层循环：遍历不同求解器类型，执行一类求解器的所有计算操作。不同求解器通过内部接口定义实现相关数据共享。

3) 第3层循环：遍历每个分区网格，执行对应求解器类型指定的计算流程。该层循环结束时，执行并行通信，实现不同分区网格交接面的数据交换。

在求解器内进行并行通信的设计，使得开发人员只需关注某一类求解器的并行通信数据，却不用考虑其他求解器的通信要求，从而给不同求解器耦合计算程序的开发提供便利。

由于LAURA软件仅应用于结构网格，因此它并没有上述耦合算法设计框架。NNW-HyFLOW的耦合算法设计不仅满足了不同类型网格一体化计算的需要，同时给求解器层面的耦合计算提供了方便易用的接口，进一步增强了软件的可扩展性。

1.4 并行计算

NNW-HyFLOW的并行算法设计如图7所示，将求解器列表捆绑到分区上，每个分区对应一个进程，且每个分区里包含多块网格，其并行计算过程同样见图6，当一种求解器类型的计算执行完成后，每个分区之间执行一次并行通信，同时在分区内部网格间执行一次交接面数据交换。当求解器列表遍历完成后，也要进行求解器之间的数据通信，以实现不同求解条件下流场或特殊边界的数据传递和交换，满足多物理场耦合计算的需要。

图7 并行分区示意图Fig.7 Diagram of parallel regions

在数据通信层面，NNW-HyFLOW直接采用了风雷开源软件提供的并行通信接口，其核心步骤包括数据压缩、MPI通信和数据解压[17]，如图8所示。为适应不同类型网格耦合并行计算的需要，NNW-HyFLOW同样采用了风雷开源软件设计的“数据容器”[16](DataContainer)结构，专门用于MPI层面的各种数据传递和交换，满足差异化数据的统一存取和并行通信。这种并行通信接口设计大大简化了并行计算的设计难度，提高了软件开发的效率。

图8 并行计算通信接口示意图Fig.8 Diagram of parallel communication interface

目前，NNW-HyFLOW已经完成了结构网格CFD求解器的并行计算代码实现，其根据图5所示的数组退化或拓展方式对网格交接面的密度、速度、压强、温度以及组分质量分数等相关数据进行封装，然后调用并行通信接口函数进行推送即完成并行计算代码实现。基于风雷开源软件提供的并行通信接口，NNW-HyFLOW软件的并行设计具有简单易用、扩展性强等优势，这是LAURA软件不具备的能力。

1.5 接口设计

与气动物理、流固耦合等NNW多学科软件的耦合计算主要采用完全解耦策略。因此，NNW-HyFLOW软件在后处理模块中预留了相关输出函数接口，将根据其他软件的格式要求输出流场或边界面的相关参数。同时，在应用层面同样预留了相关求解器接口，当其他软件的求解器有基于风雷开源平台进行程序开发的需要时，可以快速开展代码移植和求解器耦合计算。快速代码移植和便利的耦合计算程序设计是NNW-HyFLOW区别于LAURA软件的又一项优势。

NNW-HyFLOW软件的最大功能特色是热化学非平衡效应模拟。针对复杂的化学模型计算，NNW-HyFLOW设计了专门的函数接口用于处理不同类型的化学反应模型。类似“数据容器”的设计方式[16]，NNW-HyFLOW采用“化学模型容器”对内置化学模型或其他类型的化学反应模型进行标准化统一管理。如图9所示，“化学模型容器”包含3类主要信息：基本信息、化学反应式、反应速率模型。基本信息包括组分总数、反应式总数，以及各组分的化学名称。化学反应式部分以3个二维数组的方式分别列举和存储各个组分在每个化学反应式的前、后反应中的化学计量系数，以及作为第三碰撞体的系数。反应速率模型部分将列举各反应式前、后反应式速率的拟合参数形式，这些参数同样采用二维数组的方式存储。

图9 化学模型容器及接口设计示意图Fig.9 Design diagram of chemical model container and its interfaces

当前，NNW-HyFLOW软件通过“化学模型容器”内置和提供了Gupta、Dunn-Kang和Park等3类常用化学反应模型，其接口函数能够根据用户输入的气体组分列表自动匹配相关反应式，同时构建用于当前计算的最小集合化学模型，实现化学反应模型的自动化设置。这一设计提高了软件的计算效率，同时也大大增强了软件的易用性和灵活性。

另外，用户利用“化学模型容器”也可自定义化学反应模型，其输出接口函数可以自动生成参考化学反应模型的模板文件，用户可利用界面程序对相关参数进行修改，或直接打开模板文件修改有关参数，计算时再利用输入接口函数导入“化学模型容器”即可用于后续流场计算。LAURA软件同样具备化学反应模型自定义功能，但在易用性、灵活性和通用性方面，NNW-HyFLOW的自定义操作更具优势。

2 理论方法

基于前述软件框架的描述，以下仅从控制方程、物理化学模型以及核心数值方法等方面，对NNW-HyFLOW软件的解算模型及其理论方法进行简单介绍。

2.1 控制方程

NNW-HyFLOW软件主要基于Navier-Stokes控制方程组，通过对完全气体层流控制方程的继承、派生与扩充，实现二维/三维热化学非平衡流控制方程体系[1-2,5,11-12]的构建。其最核心的控制方程形式为

(1)

式中：Q为守恒量向量；F、Fv分别为对流项通量和黏性项通量；S表示源项。

高温气体化学反应和热力学激发及其非平衡过程的表征，主要通过非平衡源项S，与流动方程中各气体组分质量、动量和能量的时间演化与空间输运过程相联系。通过数据结构和函数的重载，实现多种气体模型和控制方程的统一。对于完全气体，Q=[ρ,ρu,ρv,ρw,ρE]T，S=0；对于单温度化学非平衡气体，Q=[ρ,ρu,ρv,ρw,ρE,ρi]T,S=[0, 0, 0, 0, 0,Θi]T；对于两温度热化学非平衡气体，Q=[ρ,ρu,ρv,ρw,ρE,ρi,ρev]T,S=[0, 0, 0, 0, 0,Θi,Θv]T；对于三温度热化学非平衡气体,Q=[ρ,ρu,ρv,ρw,ρE,ρi,ρev,ρee]T,S=[0, 0, 0, 0, 0,Θi,Θv,Θe]T；对于多振动温度热化学非平衡模型[12],Q=[ρ,ρu,ρv,ρw,ρE,ρi,ρevi,ρee]T,S=[0, 0, 0, 0, 0,Θi,Θvi,Θe]T。这里ρ为气体密度，u、v、w分别为直角坐标系x、y、z方向速度；E、ev、ee分别为气体的内能、振动能和电子能，对于不同气体模型，其物理实质、涵盖范围及表达形式存在差别；ρi和evi分别为第i个气体组分的密度和振动能；Θi、Θv、Θe分别为化学反应生成源项、振动非平衡源项和电子非平衡源项，Θvi为第i个组分的振动非平衡源项。

控制方程中组分ρi的质量连续性方程可以采用多种模式[18]：

1)i=1, 2,…,n。当气体组分不包含离子时,n=ns-1(ns为总的气体组分数)，结合ρ=∑ρi计算得到ρns；当气体组分包含离子时,n=ns-2，进一步结合电中性原则计算得到ρns-1。这样处理的优点在于可较好地继承完全气体方程形式及其软件模块，缺点在于相关非平衡偏导数形式较为复杂。

2)i=1, 2, … ,ns，不计算质量连续性方程，由ρ=∑ρi计算得到。该方法对于完全气体模块的继承性稍差。

3)i=1, 2, …,ns，同时计算ρ的连续性方程。其优点在于较好的继承完全气体软件模块，同时非平衡偏导数形式相对简单；缺点在于计算量较大，同时方程超定，需要额外进行质量归一化和电中性修正。

由于NNW-HyFLOW软件将兼顾多种物理效应耦合的模拟需求，因此除核心的非平衡控制方程组外，还需耦合考虑湍流控制方程、气体辐射方程、电磁麦克斯韦方程以及固体温度场热传导方程等控制方程。为满足多物理效应模型的多样化需求，这些控制方程采用各自独立的构架和耦合接口设计。

2.2 物理化学模型

在NNW-HyFLOW软件框架上，完全气体相关的热力学特性参数在图4所示的PerfectGas类中进行包装和实现，而热化学非平衡计算相关的物理、化学模型及其相关参数则在RealGas类中进行封装和实现。以下分别对热力学特性参数、化学反应生成源项及化学反应模型、振动-电子能量方程源项的计算方法进行简要介绍。

完全气体条件下的热力学状态参数采用理想气体状态方程计算。对于热化学非平衡气体，气体压强采用Dalton分压定理[1]计算，各组分的热力学状态参数仍由理想气体状态方程计算。对于比热、内能以及焓等其他热力学特性参数，热力学温度模型采用分子理论建模和计算，组分i的平-转动能etr,i、振动-电子能eve,i、绝对焓hi等特性参数的计算方法参考文献[1-2,12]。为简化计算过程，单温度模型通常采用CHEMKIN软件[19-20]给出的拟合公式进行计算。

为计算黏性系数，NNW-HyFLOW集成了文献[11,20]提出的两种拟合方法，其中Blotter方法作为默认选项。Eucken关系式[2]用于计算组分热传导系数，混合气体的总黏性系数和热传导系数均采用Wilke混合律计算。在假设Schmidt数为常数和考虑湍流对输运系数影响的基础上，组分质量扩散系数采用文献[2]给出的方法计算。

化学模型计算包含在源项的计算过程中，组分输运方程生成源项Θi的计算形式为

(2)

NNW-HyFLOW软件目前可提供的化学模型包括Gupta模型[19]、Dunn-Kang模型[11,21]和Park模型[22]，应用时可根据计算需要选取相应的化学模型和组分个数构建反应式组合及其反应速率系数。生成速率Rr项中的正、逆反应速率系数kr,f和kr,b均采用Arrhenius经验公式[18]进行计算，Park模型的逆反应速率系数已转换为Arrhenius公式所需的参数形式，以方便程序的编写以及保持计算形式的统一。

对于三温度模型，振动能量方程生成源项Θv和电子能量方程生成源项Θe的计算形式分别为

(3)

两温度模型方程通过振动能量方程与电子能量方程迭加得到，因此两温度模型的能量方程生成源项由Θv和Θe相加得到。多振动温度模型源项及其相关参数模型方法参考文献[12]。

2.3 数值方法

当前，NNW-HyFLOW默认条件下均采用LUSGS隐式时间推进方法进行流场数值迭代求解。LUSGS方法的统一形式[23]为

(4)

式中：Q为守恒量向量；R为右端项；D、L、U分别为通量项分裂得到的对角阵、下三角阵和上三角阵。热化学非平衡计算条件下矩阵D包含源项S的雅可比矩阵Z。为此，NNW-HyFLOW软件针对源项处理提供了全隐式和点隐式两类计算方法。由于点隐式方法编程简单，计算效率高，且稳定性和收敛性也较好[5]，故默认条件主要采用点隐式方法计算非平衡源项。

另外，底层数值算法库可提供Roe、Steger和AUSM类等多种通量数值格式，以及Vanleer、Minmod、Vanalbada等常用限制器用于计算对流通量，黏性通量主要采用中心格式计算。

3 算例验证

完全气体条件的算例测试与验证可参考文献[17]。针对热化学非平衡效应，以下从激波位置、电子数密度、非平衡气动力和气动热环境等非平衡流场特性的计算精准度以及对复杂工程外形的适用性等方面对NNW-HyFLOW软件已完成的求解器功能进行校验和确认。

3.1 HEG风洞试验

针对HEG激波风洞圆柱模型试验[24-25]开展数值模拟测试。计算采用单温度、5组分Dunn-Kang化学模型，壁温取Tw=300 K。计算网格如图10所示，其对称面网格量为129×91，第1层高度取Δh=1.0×10-6m。

图10 HEG圆柱半模型及计算网格Fig.10 Half model of HEG cylinder and its computational grids

图11给出了横截面(z=0.01 m)的压强等值线分布，同时图12和图13分别给出完全催化壁(FCW)和非催化壁(NCW)条件下计算得到的壁面压力分布和热流分布。可以看到，本文计算得到的脱体激波位置与试验和文献结果非常吻合，壁面压力值分布与试验结果一致，同时两种壁面条件计算得到的热流分布趋势符合实际物理机制和理论分析，证明了NNW-HyFLOW软件具有可靠的非平衡流场模拟和激波捕获能力，在高温气体非平衡气动力和气动热预测方面具有较高的计算精度。

图11 激波位置对比Fig.11 Comparison of shock wave positions

图12 壁面压力分布Fig.12 Pressure distribution on wall surface

图13 壁面热流分布Fig.13 Heat flux distribution on wall surface

3.2 RAM-C飞行试验

这里以RAM-C钝锥模型[26-27]为模拟对象，考核NNW-HyFLOW软件常用化学反应模型的正确性，以及对等离子体分布预测的计算精度。测试飞行高度H=71 km，飞行速度U∞=7 650 m/s(Ma∞≈25.914 4)。计算采用单温度模型以及7组分的Dunn-Kang、Gupta和Park等3种化学模型，来流氧气(O2)和氮气(N2)的质量分数分别为0.233和0.767，壁温取Tw=1 500 K。计算网格如图14所示，其对称面网格量为161×64，第1层高度取Δh=1.0×10-5m。

图14 RAM-C钝锥模型及计算网格Fig.14 RAM-C blunt cone and its computational grids

图15给出了轴向x/R=8.1处沿法向(y/R)的电子数密度分布。可以看到，3种化学模型分别计算得到的电子数密度曲线与AEROPH_Flow程序[12]的计算值非常吻合，与飞行试验测量值差异小，处于相同数量级。测试结果证明了NNW-HyFLOW软件所含化学反应模型的正确性和有效性，同时显示了软件在电子数密度分布预测方面具有较高的计算精度。

图15 x/R=8.1处法向电子数密度分布Fig.15 Electron density distribution along normal direction at position of x/R=8.1

3.3 Electre钝锥模型试验

这里选取Electre钝锥模型[12,28]为研究对象，参考飞行试验和Muylaert等的数值分析结果[28]开展模拟测试，考核NNW-HyFLOW软件对热化学非平衡壁面热流的模拟能力。计算采用两温度、7组分Gupta化学模型，来流O2和N2的质量分数分别为0.233和0.767，壁温取Tw=343 K。计算网格采用类似RAM-C外形的网格划分，半模型网格量约72万，第1层网格高度取Δh=1.0×10-6m。

图16和图17分别给出完全催化壁(FCW)和非催化壁(NCW)条件下计算得到的驻点线组分质量分数分布和温度分布，图18给出了采用两种壁面边界条件计算得到的轴向壁面热流分布。图中：T为平动温度；Tv为振动温度。可以看到，本文计算的组分质量分数分布曲线与Muylaert数值结果吻合，温度分布曲线与Muylaert数值结果一致。另外，两种壁面条件计算的热流值曲线与测量值相比分布趋势合理，与AEROPH_Flow的计算值几乎重合，与Muylaert的计算结果符合良好。测试结果证明了NNW-HyFLOW两温度模型求解器及化学反应模型的正确性和可靠性，同时显示了软件在热化学非平衡壁面热流预测方面具有较高的计算精度。

图17 驻点线温度分布Fig.17 Temperature distribution on stagnation line

图18 轴向壁面热流分布Fig.18 Heat flux distribution on wall along axial direction

3.4 航天飞机STS飞行试验

选取“哥伦比亚”号(OV-102)航天飞机飞行试验条件[12,29]进一步开展数值模拟测试。测试飞行高度H=57.2 km，飞行速度U∞=4 560 m/s(Ma∞≈14.25)，攻角α=40.56°。计算采用单温度、5组分Park化学模型，来流O2和N2的质量分数分别为0.233和0.767，壁温Tw=1 500 K。计算网格如图19所示，半模型网格量约652万，第1层网格高度取Δh=1.0×10-6m。

图19 航天飞机OV-102模型及计算网格Fig.19 Space shuttle OV-102 model and its computational grids

图20和图21给出对称面O2、NO的质量分数等值线分布，以及壁面的压强、热流分布云图，同时图22给出了采用非催化壁(NCW)条件计算得到的机身迎风面中心线的热流分布。可以看到，本文计算得到的壁面热流值与STS-2、STS-5等飞行试验测量值[29]符合较好，与AEROPH_Flow计算值基本吻合，显示NNW-HyFLOW软件在工程型号复杂外形化学非平衡流动及其影响的气动热特性模拟方面具有较高的计算可信度，能够适应高超声速复杂飞行器高温非平衡流动数值模拟的需求。

图20 对称面主要组分分布云图Fig.20 Distribution contour of main species on symmetry plane

图21 壁面主要气动参数分布云图Fig.21 Distribution contour of main aerodynamic parameters on surface

图22 机身迎风面中心线热流分布Fig.22 Heat flux distribution on centerline of windward side of fuselage

4 结束语

基于风雷开源软件的基础框架和面向对象的编程思想，开发的NNW-HyFLOW高超声速流动模拟软件具有底层代码复用、功能兼容性好、扩展能力强以及接口灵活等优点。当前测试版本初步实现完全气体、化学非平衡气体和热化学非平衡气体数据结构的一体化设计，集成了常用空气化学反应5组分、7组分和11组分Park、Gupta和Dunn-Kang模型，主要构建了热力学一温度、两温度和三温度能量输运体系，初步具备了多种气体模态热化学非平衡效应数值模拟能力。结合HEG风洞圆柱模型试验、RAM-C飞行试验、Electre钝锥模型试验和航天飞机STS飞行试验条件，开展了高温气体非平衡流数值模拟校验。软件测试表明：该软件测试版本在热化学非平衡效应及其影响的气动力特性、气动热环境和等离子体分布特性预测与评估方面，具有较高的数值计算精准度，初步满足了高超声速复杂飞行器高温非平衡流动数值模拟的需求。

为了进一步提升NNW-HyFLOW软件的综合模拟能力，后续将在以下方面开展优化：

1) 在网格架构多样性方面，基于现有数据结构和功能模块，实现结构/非结构混合和重叠网格扩展，更好地满足复杂外形飞行器局部细微结构精细模拟能力。

2) 在功能的完整性方面，继续集成国内各类物理化学模型，包括但不限于新的地球大气热化学反应体系、火星大气反应体系、多种燃气反应体系以及烧蚀反应体系；进一步完善多物理效应耦合模拟模块，包括但不限于湍流效应、稀薄滑移、气体辐射、表面有限催化及电磁效应的耦合策略和接口设计。

3) 在代码的可扩展性和模块的可移植性方面，继续加强软件各功能模块的封装，加强各级数据之间的分层管理，进一步提升软件移植性。

4) 在程序执行的效率和稳定性方面，持续优化并行底层封装和并行原理构建，结合多种隐式耦合算法，提升软件的计算效率和数值格式的鲁棒性。

致 谢

本研究得到了中国空气动力研究与发展中心计算空气动力研究所风雷课题组和气动物理课题组的大力支持，在此表示感谢。