那旭东，夏智勋，马立坤，颜小婷

(国防科技大学 空天科学学院，长沙 410073)

0 引言

固体火箭发动机是火箭的动力装置，由于其具有结构简单、使用方便、可靠性高、成本低等特点，常作为一级或多级的动力装置被广泛应用于运载火箭、战术武器以及战略武器等方面[1]。固体火箭发动机的设计涉及推进剂的点火与燃烧，推进剂、壳体、绝热层及喷管的结构力学，内部流场及外部羽流的流体动力学，含能材料物理化学特性，组件的老化损坏以及不同潜在失效模型分析等一系列学科[2]。此外，固体火箭发动机工作过程涉及不同空间(μm～m)、时间(10-6～102s)尺度的燃烧及流动[3]，在燃烧室内点火器点火产生出高温高压的气体点燃固体推进剂，推进剂燃烧的气体产物生成压力作用于固体药柱表面，促使药柱变形，另外气体产物向药柱传递热量，促进推进剂的快速燃烧。同时，药柱的变形与快速燃烧又影响着燃烧室内气体流动与压力分布等[4]。因此，固体火箭发动机的工作过程是一个复杂的、强耦合的物理化学过程。

自20世纪70年代以来，采用计算机仿真技术对固体火箭发动机进行辅助设计及工作过程模拟已成为固体火箭发动机研究的一个重要手段。近年来，随着计算机性能的提升以及CFD技术的发展，各国对固体火箭发动机工作过程的数值模拟也变得越来越精细，从零维到多维、从不可压到可压、从定常到非定常、从无粘到有粘、从单相到多相等各个方面不断发展[5]。各国在对固体火箭发动机工作过程进行数值仿真时，除了采用Fluent、Star-CD及ABAQUS等商业软件外，还开发了许多相应的设计仿真程序和软件[6-12]，并取得了一系列可喜的成果，其中以美国先进火箭仿真中心CSAR(Center for Simulation of Advanced Rockets)开发的Rocstar软件包最具有代表性[13]。

Rocstar能够对固体火箭发动机从点火到稳态燃烧工作过程进行三维、全耦合及高精度仿真模拟，求解还包括铝颗粒、烟尘和各种化学成分的流场[14-15]。本文主要就CSAR的由来，Rocstar的发展历程、程序架构组成、求解模块及对应的算法、程序验证和求解流程几方面对Rocstar进行介绍。

1 Rocstar开发背景、机构组成及发展历程

1.1 开发背景和机构组成

20世纪90年代，美国为了保持战略平衡，隶属于能源部的国家核安全管理局发起了先进模拟与计算ASC (Advanced Simulation and Computing Program)项目[16-17]。ASC项目的目标是提升先进的计算建模和仿真能力来实现对复杂、多组件系统的大规模并行仿真模拟，以便在没有地下核试验的情况下，通过数值模拟试验来检测其核武器的安全性与可靠性。

伊利诺斯大学为响应ASC计划，于1997年成立了先进火箭仿真中心CSAR[18]。CSAR希望开发一套固体火箭发动机虚拟样机工具，用于开展固体火箭发动机主要组成部件及部件间动态相互作用的详细建模和仿真，从而对整个固体火箭发动机工作状态进行精确的预测[19-20]。

为完成固体火箭发动机全系统仿真的这一终极目标，CASR将开发人员分成了燃烧和含能材料组、计算科学组、流体力学组、结构和材料组4个大的研究组以及为了完成特定问题而建立的其他小组。从1997年创建CSAR到2006年底，将近8年时间，CSAR完成了立项之时所设定的目标[4]。此后，Illinois Rocstar LLC作为一个独立的商业部门和IR&D公司开始负责Rocstar的开发和技术支持，但仍由美国能源部提供资金支持[13,21]。

1.2 发展历程及可实现的功能

CSAR在开发Rocstar程序过程中采用了逐步深入的开发策略[22]，图1给出了Rocstar各个版本所能处理的几何及物理模型的复杂程度。随着版本的不断提高，其所能处理物理问题也越来越精细，几何模型也越来越复杂。

当版本发展到Rocstar 3时，程序中已包含有许多物理应用求解器，这些求解器通过一个强大的耦合框架整合在一起。所有求解器都可进行高效且大规模的并行求解，从而可以对异常复杂的问题进行精细化的仿真计算[23]。图2为Rocstar 3中一系列可供用户选择的计算模型和算法。

图1 Rocstar各版本功能Fig.1 Functions of Rocstar in each version

图2 Rocstar 3中可用的模型和算法Fig.2 Models and algorithms in Rocstar 3

在Illinois Rocstar LLC接手Rocstar 3后，将其最新一代的多物理场仿真软件包Rocstar 3代码开源[24-25]，并将其命令为Rocstar仿真套件(Rocstar Simulation Suite)。

2 Rocstar架构及模块功能

2.1 Rocstar开发思想及架构组成

固体火箭发动机是一种性能优越的火箭动力装置，其结构简单，其主要由点火器、壳体、推进剂、绝热层和喷管等部件组成，如图3所示。

固体火箭发动机的工作本质是燃烧流动过程进行的两次能量转换，其机理复杂、环节众多：如在燃烧室多相流动过程中就涉及金属分散燃烧、凝相破碎/聚合、熔渣沉积、热结构烧蚀、多相辐射对流热传导等过程，如图3所示。此外，组成发动机各部件间的耦合是强非线性的，工作过程涉及不同空间、时间尺度的燃烧及流动，因此其工作过程是一个复杂的、强耦合的物理化学过程，涉及大量的学科和科学问题。

图3 固体火箭发动机基本组成部件及相关物理化学过程Fig.3 Basic components and physicochemical processes of SRM

从基本原理出发，进行详细、全部件耦合的固体火箭发动机仿真模拟，需要开发代码来实现每个火箭部件的数学模型以及它们之间的相互作用。这不仅需要更大的计算能力，而且要求程序中各求解模块的数值正确，且相互保持一致的物理作用。因此，进行固体火箭发动机全系统的建模和仿真极为复杂。

CSAR在开发Rocstar时，多物理场系统集成的方法是开发一个主程序和一套接口代码，以及一系列可单独运行的求解程序，然后通过接口代码实现主程序对相应求解程序的调用[4]。Rocstar仿真应用程序的体系结构见图4，包括了大量用于对流体-固体-燃烧三者之间相互作用、耦合的组件程序。

图4 Rocstar仿真套件体系结构Fig.4 Rocstar simulation suite architecture

从图4可看出，Rocstar主要分为前处理、物理模块求解器、集成架构以及后处理4个功能区域，其中集成框架是整个Rocstar的核心所在，包含了许多用于支持多物理场耦合仿真的服务及协调模块。

在Rocstar开发过程中，采用了面向对象的程序开发思想，服务及协调模块采用C++语言编写，物理应用求解器采用Fortran 90语言编写并采用MPI实现并行仿真。Rocstar可在集群平台上进行固体火箭发动机一体化、全系统耦合的大规模仿真计算[22]。

2.2 协调及服务模块

Roccom 是Rocstar最底层的框架。Roccom可实现不同程序模块中数据及函数交互、内部程序的输入输出(Rocin,Rocout)、网格之间的通信操作(Rocsurf,Rocmap)等功能，其中程序模块可支持使用不同程序语言，主要为C++和F90。通过调用Roccom，物理模块相互之间可以进行通信，同时也可以访问一系列服务，包括网格修改、数据传递、数据输入输出、性能分析及负载平衡等。

Rocman是Rocstar的主要驱动和协调服务模块，主要执行和管理不同组件程序间的相关作用，控制物理应用的执行过程，包含数值仿真的初始化、耦合时间步、交界面处的跳跃条件以及输出存储和仿真终止等。Rocman最新的版本为Rocman 3，采用C++编写，其实现的功能比以前更加完善和通用。

Rocprop、Rocmop、Rocrem及Rocface用于耦合界面演化与数据传递。

Rocprop是不同域之间接触面演化服务模块，主要为Rocstar提供拉格朗日面追踪及演化服务，用于计算推进剂表面由于燃烧所带来的退移过程。Rocprop既可以用来进行多域耦合仿真也可以单独进行流体仿真或结构力学仿真。Rocprop 实现了两种分界面演化算法：(1)早期采用的粒子标记方法；(2)更先进的面位移算法。面位移算法(FOM)对于处理边角附近的面的运动时很有优势[26]。

对于仿真过程涉及诸如固体药柱燃面退移、药柱裂纹动态发展等几何构型发生改变的情况，为保证网格的质量，在Rocstar中有两种网格修改方案来处理不同等级的网格失效情况。非结构网格的光顺(不改变网格节点数量及节点连接情况)是通过Rocmop调用Mesquite程序包完成的，Mesquite[27]是由美国桑迪亚国家实验室开发的网格算法程序。每一个分区同时调用Mesquite程序，用于实时提供真实和虚拟(在计算域外)的节点。如果光顺不能提高网格的质量，则采用Rocrem模块对网格进行局部的修改或重新进行网格生成，Rocrem通过调用商业软件包Simmetrix来完成这一工作。当网格重新生成后，采用Rocface模块将计算结果映射到新的网格上继续进行计算。Rocface[28]是一个数据映射服务模块，可将独立离散面之间的求解数据进行传送。

2.3 Rocstar前处理

对于几何建模，只要CAD软件可生成被网格划分软件读取的几何模型就可以采用。CSAR开发人员采用Pro/Engineer来生成流体及固体域的CAD模型，并输出为通用几何模型接口IGES格式。对于网格划分，CSAR开发人员指定了Rocstar前处理中的网格及边界条件的格式。目前，对于流体力学程序采用Gridgen软件来生成网格及边界条件，对于结构力学程序采用Partran或Truegrid来生成对应的网格及边界条件。当网格和边界条件以指定格式输出后，在输入数据集准备工具Rocprep的辅助下为每个物理仿真应用创建数据输入集以及为并行运算进行分区划分。

2.4 求解器

2.4.1 流体力学求解器

Rocstar中流体动力学求解器主要有两个，分别为Rocflo和Rocflu，主要用于求解可压缩多相流。两个求解器均采用基于体心的有限体积离散方法求解。

Rocflo采用任意拉格朗日-欧拉ALE (Arbitrary Lagrangian Eulerian)算法在多块结构网格上进行求解[28]，ALE可以考虑由于推进剂变形或燃烧带来的边界运动情况。Rocflo求解器空间离散格式采用中心格式或二阶TVD格式。时间离散，除了显示的Runge-Kutta方法，Rocflo还可以使用一个双时间步算法来处理比CFL条件限制下的时间步长更大的时间步长。Rocflo的并行划分是基于块拓扑结构，每个处理器包含1个或多个网格块。

Rocflu可在非结构四面体网格或混合四面体/六面体/金字塔/棱柱网格上进行计算，该求解器所实现的动网格(ALE算法)技术可用来模拟推进剂燃烧及变形带来的拓扑结构改变的情况。在采用有限体积法进行方程离散时，采用了高阶的类WENO法，即HLLC法，该方法可处理类似点火流动等强瞬态过程。时间积分方法可通过3阶或4阶显式多步Runge-Kutta时间步进算法。

这两个流体动力学求解器都可以选择是否包含有湍流(Rocturb)、拉格朗日粒子追踪(Rocpart)、欧拉烟雾追踪(Rocsmoke)、化学反应(Rocspecies)以及辐射(Rocrad)等模型。

2.4.2 结构求解器

Rocstar主要有两个有限元结构力学求解器，Rocfrac和Rocsolid[13]。这两个求解器中都包含有ALE算法，因此都可以考虑固体推进剂由固相到气相的变化过程。此外，还可以处理大的应变、旋转过程、三维热传导等问题。

Rocfrac是一个计算结构力学和瞬态热传导求解器，该求解器采用非结构网格，可进行显式和隐式求解。Rocfrac支持多种材料力学模型，主要为线弹性模型、Arruda-Boyce模型及Neoh ookian模型。Rocfrac在普通的四面体单元之间还可以包含有结合力/体有限单元，采用这种单元可以对药柱裂纹的发展过程进行仿真计算。

Rocsolid是一个有限应变结构求解器，它采用隐式时间积分方案，并且包含有多重网格和稳定双共轭梯度在内的多个高性能稀疏迭代求解器，可以进行高效并行线性系统的求解。对可压缩或不可压缩材料Rocsolid均包含多种求解单元类型。Rocsolid可以采用广义有限元法(GFEM)来捕捉不连续性，例如运动前缘或裂纹，而不需要繁琐的几何问题计算。除此之外，Rocsolid可利用多个基于微观力学的模型来解释微观结构的演变，包括裂纹形成和生长的影响，损伤(脱湿)以及应变硬化。

2.4.3 燃烧求解器

Rocburn可仿真推进剂加热、点火及燃烧情况下燃面的退移速率或化学反应。Rocburn有3个燃烧求解模型用于计算推进剂燃率，物理模型在形式上是一维(沿燃面法向方向)，但是在推进剂燃面上的每个单元面都独立的施加，使其相当于三维模型。这3个模型中最简单的模型称为RocburnAPN，采用著名的维也里公式r=apn，其燃速大小正比于当地燃气压强的n次方。另外两个燃速模型为动态燃速模型，通过求解一维非定常热传导方程以获得温度场从而捕捉点火瞬间过程。其中一个动态燃速模型称为RocburnZN，基于Zeldovich-Novozhilov方法开发；另一个称为RocburnPY，采用一个简单的分解定律来计算燃速。RocburnPY可以计算点火器产生的高温气体对推进剂的热传导过程以及推进剂的点火过程。在一个全系统仿真里，RocburnPY可以调用由CSAR开发的详细的三维推进剂燃烧仿真代码来计算推进剂的燃速，此时可考虑推进剂配方对燃速的影响。

Rocfire 用于模拟尺度在1～10 μm的异质复合推进剂的燃烧，异质推进剂颗粒填充几何模型采用Rocpack生成。Rocfire考虑固体和气相之间详细的物理耦合，可考虑的现象包括燃面的非均匀退移、固体中的非稳态热传导、金属颗粒的存在、粘合剂/氧化剂的分解与火焰。通过Rocfire的模拟结果可提供气相火焰结构、燃速特性、涡度和湍流的产生的信息。

2.5 后处理

CSAR开发了一款功能强大的科学可视化后处理软件包Rocketeer[29]。Rocketeer基于VTK开发，可在AIX、Linux、Windows及Solaris等系统上运行。Rocketeer有一个非常友好的图形界面，可用其读取HDF文件格式或CGNS文件格式数据。Rocstar还开发了Rocketeer的批处理并行版本Voyager，以及客户端/服务器版本Apollo/Houston。

除使用CSAR开发的专用后处理Rocketeer外，还可以使用商业软件Tecplot和开源软件ParaView对仿真结果进行处理分析。

3 耦合方案

Rocstar的特点之一就是可处理三维情况下具有共同交接面的多域耦合问题。Rocstar要求各域可变形或移动但不能分开，也不能重叠。在交接面处，质量、动量及能量等保持守恒，一些物理量，如密度等被设置为跳跃边界条件。

Rocstar采用分域计算策略，不同计算域单独进行网格划分，对于不同的物理求解域都有对应的求解器，比如流体域的Rocflo和固体域的Rocfrac。每个域的计算过程是独立的，只在每个全局时间步长计算后，将交接面的数据通过Rocface模块传递给相邻域求解器。例如流体域在计算一个全局时间步长后，将压力分布加载到交接面上并传递给固体域，固体域进行一个全局时间步长的计算，得到交接面的变形及运动速度，又将其传递给流体域，整个过程一直重复直到计算结束为止，如图5所示。交接面处可考虑燃烧或化学反应过程，由燃烧或变形造成的计算域几何形状发生极端变化时，需要对网格进行光顺或重构处理。

图5 Rocstar多物理接口耦合Fig.5 Rocstar multiphysics interface coupled

Rocstar 采用“分割”方法来处理多域耦合时间步进[23，30]。在这个方法中，一个全局时间步长Δtn内，流体域的时间步长Δtf和固体域的时间步长Δts是不同的，由CFL条件来控制其大小。当前域物理求解器在达到了下一全局时刻tn+1时，会通过交界面与相邻域之间进行数据传递。只要全局时间不比任何一个物理计算求解器的时间步长大太多，整个系统仍保持紧密的耦合。图6为Rocstar的显式时间步进耦合方案。

图6 Rocstar耦合时间步进方案Fig.6 Rocstar coupled time stepping scheme

如图6所示，在Δtn时刻计算首先从固体域开始，当固体域求解器经若干时间步Δts的计算到达下一全局时刻tn+1，将燃面的位置、速度以及由于燃烧带来的质量流量传递给流体域求解器，并采用Rocprop模块更新燃面位置。流体域求解器经过多个时间步Δtf的计算从tn时刻到达tn+1时刻。此时，流体域会将燃面处的受力加载传给固体域，并将燃面处的压强及温度传递给燃烧求解器，燃烧求解器计算出的燃速值也会传回给固体域求解器，这样就完成了一个全局时间步的计算。

对于结构求解，如果采用隐式固相结构求解器Rocsolid代替Rocfrac，还可以采用基于预测-修正的隐式时间步算法。

4 Rocstar程序的验证

CSAR在开发Rocstar过程中，选择航天飞机可重复使用固体助推器作为其主要仿真对象，NASA向CSAR提供该发动机的各种数据，以供CSAR的研究人员对Rocstar进行各类验证[22]。

在开发Rocstar过程中，为了对各模块及物理求解器进行验证，CSAR开发人员进行了相应的数值实验对质量守恒、网格精度、系统时间步和时间精度，以及ALE算法的可扩展性进行了研究，在开发物理求解器过程中还将计算结果与其他求解器进行对比以保证Rocstar的计算准确性与精度[31-33]。

除了航天飞机固体助推器发动机，CSAR还对多台战略火箭发动机进行了仿真计算，证明了其仿真结果具有很高的可信度。例如，CSAR开发人员使用Rocstar对1991年美国“大力神4”助推固体火箭发动机(Titan IV SRMU)在首发试车时发生爆炸的原因进行了数值模拟[34-35]。通过计算发现，该发动机是由于燃气沿药柱通道流动存在压强分布，导致后段药柱变形，在连接槽处推进剂发生塌陷，从而减少了燃气流动的有效气流面积，加速了压强分布的不均匀性，而这又加速导致药柱的变形，如图7所示，仿真结果与理论推断相吻合[36]。

图7 变形推进剂中空隙的分布Fig.7 Distribution of voids in the deformed propellant

5 Rocstar仿真求解流程

Rocstar的仿真流程与多数CAE仿真软件(如CFX、Fluent和OpenFOAM等)的仿真流程类似，分为前处理、运算求解和后处理三部分。本节以Rocstar自带的姿态控制发动机ACM(Attitude Control Motor)算例为例介绍采用Rocstar进行固体火箭发动机仿真时，其仿真求解流程及相关模型、参数的设置方法。

该姿控发动机的构型如图8所示，其总长度约为55 mm，是一个小型的固体火箭发动机[33]。采用Rocstar对该发动机进行非定常计算，在考虑燃面动态退移情况下，开展发动机内流场燃烧流动耦合仿真计算。将发动机内部流动简化无粘流动，在推进剂燃面采用质量流量入口，推进剂燃速采用燃速指数公式r=apn进行求解。

图8 姿态控制发动机构型Fig.8 Attitude control motor geometry

5.1 前处理

采用几何造型软件Pro/Engineer绘制发动机几何模型并导出通用几何模型接口IGES格式文件。在Gridgen中，采用多块结构六面体网格对发动机进行划分(如图9所示)，并输出对应的网格文件及边界条件文件。

图9 发动机多块结构网格Fig.9 Block-structured hexahedral mesh of motor

图10 本地数据存档目录结构Fig.10 Set up NDA scheme

Rocstar没有图形操作界面，相关设置主要通过修改对应的配置文件进行。建立如图10所示的NDA(Native Data Archive)目录结构[37]，其中每一个矩形表示一个文件夹，箭头指向为当前文件夹的子文件夹。将相应的网格文件(ACM-PLOT3D.grd)、边界条件文件(ACM-PLOT3D.inp /ACM.bc)、边界条件映射文件(ACM-PLOT3D.bcmp)、物理模型和数值参数等信息配置文件(ACM.inp)以及各求解模块的参数文件(*.txt)放入相应的位置。

当上述工作准备完成后采用2.3节中介绍的前处理工具Rocprep，为每个物理仿真应用创建数据输入集和进行并行运算分区划分。最终生成的姿控发动机算例目录结构如图11所示，其中ACM40为姿控发动机算例求解的主目录。在ACM40目录下包含3个子文件夹：Rocflo文件夹，即RocburnAPN文件夹和Rocman文件夹，用于存放对应模块的配置参数文件、相关数据信息文件和计算结果输出文件等。

图11 姿态控制发动机算例目录结构Fig.11 Directory structure of ACM case

5.2 运算求解

前处理设置完成后，在运算求解阶段主要通过对图11中各求解模块的配置文件进行参数设置。例如在RocstarControl.txt文件中，可以对耦合方案、流体求解模块、燃烧求解模块以及时间步长、自动保存步数等进行设置。RocburnAPN文件夹内的RocburnAPNControl.txt文件用于设置推进剂燃速指数公式r=apn中的系数n和压强指数a，以及绝热火焰温度等参数。

Rocstar采用MPI进行大规模并行求解计算，在求解模型及参数设定完成后采用如下命令开始求解计算：mpirun-np核数rocstar。计算结果保存在对应模块文件夹内的Rocout文件夹。

5.3 后处理

计算完成后，采用Rocketeer后处理软件对仿真结果进行处理分析。图12为选取计算结果中4个时刻下发动机几何和网格拓扑变化，可看出随着推进剂的燃烧，燃面发生了退移。图13为z=0平面发动机内流场压强云图，可以看出随着推进剂的燃烧，发动机内流道变大且发动机内部压强逐渐升高。

图12 几何构型及网格拓扑变化Fig.12 Topology change of geometry and mesh

图13 Z=0平面压强变化Fig.13 Pressure change at Z=0 plane

6 结束语

Rocstar多物理场仿真程序提供了先进的固体火箭发动机仿真模拟能力，在固体火箭发动机全三维模型和精细物理现象捕捉、多物理场非线性耦合、复杂几何模型动态变化、空间和时间尺度上的级大差异、材料属性以及物理过程的复杂性以及大规模并行实现等方面都实现了突破。Rocstar能够对固体火箭发动机工作过程进行全三维、多域耦合、高精度且大规模并行的数值模拟，可用于预示固体火箭发动机的性能和工作效果，分析潜在可能导致发动机工作失效的问题，以及为事故原因分析和发动机结构设计的合理改进提供较准确的技术支持。

Rocstar程序目前已经开源，除采用Rocstar对固体火箭发动机开展已有模型算法的仿真计算，还可以对其源程序代码进行二次开发。例如，Rocstar中的两相流中并没有考虑颗粒之间的碰撞过程，如需在仿真过程中添加该过程，则可以对Rocpart拉格朗日颗粒追踪模块进行修改。此外，得益于其灵活、通用的集成框架和众多的求解模块，还可以采用Rocstar对固体火箭冲压发动机或固体火箭超燃冲压发动机的工作性能开展仿真计算。