李 峥，向红军

(北京航空航天大学宇航学院，北京 100191)

0 引言

固体火箭发动机燃烧室和喷管中的流场是不可分割的统一流场。燃烧室中流速较低，压力变化幅度小，可按不可压流处理，燃气在喷管中的加速膨胀过程则为可压缩流动［1］。目前，大多数固体火箭发动机药柱几何形状都呈现明显的三维特征，使燃面边界及相应的内流场亦具有三维性质。与此同时，发动机在工作过程中会出现过载情况。文献［2－3］指出，在高过载情况下，燃烧室内部Al2O3颗粒偏向过载方向，并在药柱表面沉积，受到影响的药柱燃烧部位燃速增加，药柱燃面平行层燃烧规律被打破，药柱不同部位的燃速受颗粒运动的影响差异很大，使得流场、燃面以及燃烧不再保持轴对称性。因此，有必要对固体火箭发动机燃烧室喷管进行全三维数值模拟。

以往对于求解固体火箭发动机内部流动问题，大多采用基于压力的方法，将连续方程改写为压力修正方程，解出压力场与速度场，这就是著名的SIMPLE算法［4］。由于固体火箭发动机内部流动包含不可压和可压缩问题，压力修正方程的迭代矩阵往往呈现病态，这增加了求解难度。早期工程上多采用直接求解法以及Krylov子空间迭代法求解压力修正方程，受限于其计算时间复杂度［5］。此类方法大多应用于二维或网格数较少的伪三维情况，难以求解全三维大规模实际问题。代数多重网格方法(Algebraic Multigrid—AMG)具有最优的时间复杂度，即迭代收敛速度与问题规模无关，计算量线性依赖于未知量个数。目前，该算法已成功用于求解计算流体力学问题［6－7］，但在固体火箭发动机模拟方面的相关报道较少。

本文将AMG方法引入原有基于SIMPLE算法的求解器，用于求解压力修正方程。在介绍控制方程和离散方法的基础上，给出了AMG方法的主要内容;然后，将Gauss消元法、不完全LU分解BiCGStab方法及AMG方法用于双燃速星型药柱固体发动机三维两相流场仿真，对比了3种方法在不同网格规模下的计算效率;最后，给出了AMG方法下发动机内流场计算结果，分析了颗粒相与过载对发动机内流场的影响。

1 数学模型

表1给出了SIMPLE算法求解步骤及单步循环各步所用时间。其中，求解压力修正方程的时间占用了总时间的86%，为SIMPLE算法的核心步骤。

表1 SIMPLE算法单次循环时间分布(AMG，100×40×20)Table 1 Time distribution of a single SIMPLE cycle(AMG，100×40×20)

1.1 气相控制方程及离散

稳态气相通用控制方程组

式中 φ表示气相通用变量，在三维坐标下分别为速度u、v、w及温度T;φ为1时代表连续方程;Γ表示扩散系数;Sφ为源项;U为速度矢量;ρ为密度，由完全气体状态方程 p=ρRT 确定［3］。

采用基于同位网格的有限体积法离散气相通用控制方程(1)，最终得到线性代数方程组:

其中的系数依赖于离散化过程中所使用的具体格式。对于可压缩流动，压力不仅出现在动量方程中，而且通过状态方程直接与密度以及能量方程相关联。密度的变化也会影响压力变化，所以必须建立压力-速度-密度的耦合关系式，即可压缩形式的压力校正方程:

基于SIMPLE算法的可压缩流动求解步骤与不可压缩几乎完全相同。其中:

从中可看出，该方程的系数矩阵不再具有对角占优特性(aP≠n=∑nban)，且为病态方程组，其求解出现

(i)了困难。

1.2 颗粒相求解方法

采用拉格朗日颗粒轨道模型模拟Al2O3颗粒与连续相之间的动量、能量交换，不考虑颗粒的燃烧、蒸发、碰撞、聚合等，只考虑破碎效应。采用Crowe C T［8］所提出的PSIC方法(单元中的颗粒源项法)，进行气粒两相间的耦合计算。

(1)颗粒破碎模型

当颗粒在气相的牵引下加速时，会逐渐变形直至破碎，采用韦伯数(Weber)［9］来控制颗粒破碎是否发生。文中简单地认为，当颗粒保持为液态，且韦伯数超过临界值时，颗粒破碎成直径为原来一半的粒子。

式中 Dk为颗粒直径;ρ为气相密度;为气相速度;为颗粒速度;σ为表面张力，Al2O3颗粒在2 300 K时的表面张力为0.69 N/m。

当韦伯数超过4时，球形液滴开始变形，在We=10～20范围内粒子发生破碎。

(2)颗粒尺寸分布

采用 Braithwaite［10］和 Salita 等［11］提出的颗粒分布函数，描述推进剂中Al2O3颗粒的直径分布，颗粒尺寸分布为对数正态分布。当采用质量平均直径为5 μm、标准差为0.46的8组颗粒群计算时，平均直径及所代表的尺寸范围和质量分数见表2。

表2 颗粒群平均直径Table 2 Average diameter of particle cloud

2 代数多重网格方法

多重网格方法的思想起源于20世纪60年代Brakhage H［12］、Fedorenko R［13］和 Bachvalov N［14］的工作，最早出现的是几何多重网格(Geometric Multigrid-GMG)方法。GMG方法通过在不同疏密网格下求解原方程组，平顺高、低频误差，以达到加速收敛的效果。

AMG方法是在GMG方法的思想上延伸建立起来的，它不需要利用求解问题的物理和几何背景，仅利用方程组稀疏矩阵中的代数信息，构造多重网格组成部件。这一特点使得AMG方法既保持了GMG方法良好的算法可扩展能力，又可作为具有即插即用性质的黑盒子解法器，大大增强了其鲁棒性［15］。

2.1 基本原理

对于式(3)，可写成如下的线性系统［16］:

式中 A 为 n×n矩阵;u，f∈Rn。

插值算子P将粗网格向量空间转换为细网格向量空间Rnc→Rn，为n×nc矩阵;限制算子R将细网格向量空间转换为粗网格向量空间，为插值算子的转置R=PT。

尽管实际的多重网格格式多采用几层网格，为简便起见，这里仅讨论二层网格(即细网格Ωh和粗网格Ω2h)的V循环，并先介绍最简单的粗网格修正格式:

(1)在细网格上求解

迭代μ1步后，计算所得近似值的残余:

(2)在粗网格上求解误差方程:

(3)进行粗网格修正:

(4)以此为初值，在细网格上再迭代μ2步，然后回到(1)，开始下一个V循环，直到达到一定的收敛标准为止。

经典Ac构建按Galerkin方法得到，即Ac=PTAP。该方法为能量范数形式，可更好地减少校正后的误差。

2.2 粗网格的选取与插值算子［17］

在大多数AMG方法中，嵌套的网格层由细网格变量集合与粗网格变量集合组成，从细变量集合中筛选粗变量的过程为粗网格的选取。粗网格的选取对于AMG方法性能有很大影响，其他部件的构造往往直接依赖于网格粗化的结果。本文以C-AMG方法为例，介绍网格粗化的一般方法。

假设aij为系数矩阵A中的元素，若

则称i点与j点是强连接，粗网格点由强连接方向上选出，同时文献［18］还给出了粗化过程应满足2个原则:(1)对于当前网格上不属于粗网格的节点，它的每个强连接点都应在粗网格上，或与粗网格上的至少一点有强连接;(2)粗网格点集应为细网格点集的极大独立集，即细网格点集中的任意2个粗点之间不存在强连接关系，在此前提下使得粗网格点集网格结点尽量多。文献［16］中详细描述了选取粗网格结点的具体过程，本文不再赘述。

到目前为止，算法的核心问题变为定义插值算P:

假设j点为细网格上的节点，而i为粗网格上的结点。若j点与i点重合，则令j点的误差校正量等于粗网格上i点的值，若j点与i点不重合，可由与j点相邻的粗网格点插值得到。

3 计算结果

算例所模拟的发动机由星形药柱、圆柱药柱、过渡段及长尾喷管组成(图1)，以两种燃速推进剂实现单室双推力。其中，后部的第1级为星形药柱;第2级包括前端的圆柱药柱及过渡段药柱。第1级装高燃速推进剂，第2级装低燃速推进剂。为满足导弹总体需求，采用了长尾喷管。图2给出了圆柱药柱和星形药柱横截面的网格划分。

图1 发动机结构Fig.1 Configuration of solid rocket motor

图2 网格划分Fig.2 Computational grid

推进剂由 65%AP、8%HTPB、13%RDX、10%铝粉及4%癸二酸二辛脂组成，燃烧室平衡压强为6.0 MPa，由热力计算得出燃烧室温度为3 302 K，只计算发动机初始时刻的流场情况。为考虑固体火箭发动机三维特性，对其施加一径向载荷30 g。

3.1 代数多重网格方法数值效率对比

首先，对AMG方法数值计算效率进行测试，将上述算例在不同网格密度下展开，分别对比了直接求解法(高斯消元法-GE)、不完全LU分解BiCGStab方法以及AMG方法的计算效率。

图3给出了二维网格下3种方法单次循环的时间消耗，收敛精度为10－8。从中可看出，对于网格数较少的情况，AMG方法求解效率不明显，3种方法耗时相当;随着网格数增加，3种方法的耗时均出现了不同程度增加，AMG方法的时间增长率远小于高斯消元法和BiCGStab方法;在200×80网格下，AMG方法取代直接求解法，成为耗时最短的算法。图4分别给出了BiCGStab方法与AMG方法单次循环迭代收敛次数。可看出，随网格数增加，AMG方法收敛次数基本不变，而BiCGStab方法收敛次数大大增加，增加了该算法计算时间。

表3给出了三维网格下3种方法单次循环的计算效率，收敛精度为10－5。从中可看出，直接求解法并不适用于三维流场计算，而AMG方法不受计算规模的影响。在三维网格下，其收敛效率大大提高，收敛时间仅为BiCGStab方法的1/20，且仍能保持10次左右迭代收敛。

图3 二维网格下单次循环计算时间对比Fig.3 Comparison of single cycle computational cost in 2D

图4 二维网格下单次循环迭代次数对比Fig.4 Comparison of single cycle iteration in 2D

表3 三维情况下单次循环计算效率对比Table 3 Comparison of single cycle computational efficiency in 3D

3.2 双燃速星型药柱固体火箭发动机内部流场分析

该算例的网格划分为200×40×20，采用AMG算法，可大大提高收敛效率，总计算时长不到BiCGStab方法的1/20。

图5给出了发动机内部流场纵截面的温度分布，图6给出了燃烧室后部及长尾喷管内马赫数的分布。固体推进剂药柱燃烧生成的高温高压燃气填充整个燃烧室，随后在长尾喷管作用下，转化为高速燃气喷出，温度降低，马赫数上升。在喷管扩张段靠近壁面区域出现高速低温区，这是颗粒相作用的结果。颗粒的物理特性决定了颗粒不能通过膨胀加速，将热能转换为动能。受惯性的影响，喷管扩张段对颗粒的作用降低，出现无颗粒区。在颗粒为主流的核心区域，颗粒对气体传热，减缓气体温度降低的趋势，并阻碍气体流动，同时将其挤压到无颗粒区;在以燃气为主的无颗粒区，气体温度较低，流速较高。

图5 温度分布Fig.5 Temperature distrbution

图6 燃烧室后部及长尾喷管内马赫数分布Fig.6 Mach contour at the rear parts of chamber and long-tail nozzle

图7给出了4个不同位置上的横截面云图，分别位于二级圆柱药柱(Z=0.9 m)、一级星型药柱(Z=1.6 m)、长尾喷管喉部(Z=1.88 m)及喷管出口位置(Z=2.0 m)。在星形药柱及长尾喷管内，流动均呈现明显的三维特征，颗粒分布保持星角对称性，这是由星边燃面的周向加质造成的，颗粒的分布使得流场内温度的分布亦成星角型;对于过载问题，在圆柱药柱燃烧室内轴向速度与粒子分布明显地偏向于过载作用方向(图7(a))，而喷管内的流场分布受过载影响不大。在燃烧室内部，燃气与颗粒相的流动较慢，颗粒受过载作用较大，轴向速度较小，更容易发生偏转;在喷管内部，燃气流动较快，颗粒主要受气相作用，轴向流动速度较大，受过载的影响较少。

4 结论

(1)高斯消元法和BiCGStab方法在网格数较少的计算时具有简单高效的特性，但在三维计算中耗时较多;AMG方法计算收敛时间随网格数增加线性增加，在数值模拟三维固体火箭发动机内部流动问题时，收敛时间远小于高斯消元法和BiCGStab方法;

(2)对于过载问题，在圆柱药柱燃烧室内轴向速度与粒子分布明显地偏向于过载作用方向，而喷管内的流场分布受过载影响不大;

(3)颗粒的加入使得星形药柱固体火箭发动机喷管内流动呈星角对称，颗粒对周围燃气有挤压加热作用，燃气在该区域温度升高，速度降低;在以燃气为主的喷管壁面附近，燃气温度较低，流速较高。

图7 不同横截面上的颗粒密度、温度及燃气轴向速度分布Fig.7 Particle density distribution，temperature and gas axial velocity at different cross sections

［1］向红军，方国尧，崔济亚.固体火箭发动机燃烧室喷管统一流场计算［J］.固体火箭技术，2000，22(3).

［2］Langhenry M T，Marietta M.Acceleration effects in solid propellant rocket motors［R］.AIAA 86-1577.

［3］郭颜红，梁晓庚，陈斌.双燃速星孔药柱长尾喷管发动机三维两相流场数值模拟［J］.固体火箭技术，2007，30(3).

［4］蔡国彪，王慧玉.固体火箭发动机燃烧室三维流场数值模拟［J］.推进技术，1995(5).

［5］郭飞宵，杨力，刘荣，等.大型稀疏法方程组的代数多重网格解法［J］.测绘科学技术学报，2012，29(1).

［6］MacLachlan S P，Tang J M，Vuik C.Fast and robust solvers for pressure-correction in bubbly flow problems［J］.Journal of Computational Physics，2008.

［7］Darwish M，Sraj I，Moukalled F.A coupled finite volume solver for the solution of incompressible flows on unstructured grids［J］.Journal of Computational Physics，2009.

［8］周力行.湍流两相流动与燃烧的数值模拟［M］.北京:清华大学出版社，1991.

［9］Clift R，Grace J R.Bubbles，drops and particles［M］.New York:Academic Press，1978.

［10］Braithwaite P C.Quench bomb investigation of aluminum oxide formation from solid rocket propellants(part I):experimental methodology［C］//Proceeding of the 25th JANNAF Combustion Meeting，CPIA 498，1988.

［11］Salita.Quench bomb investigation of Al2O3formation from solid rocket propellants(part II):analysis of data［C］//25th JANNAF Combustion Meeting，CPIA 498，1988.

［12］Brakhage H.Uber die numerishe behandlung von integralgleichungen nach der quadraturformelmethode［J］.Numer.Math.，1960.

［13］Fedorenko R P.A relaxation method for solving elliptic difference equations［J］.USSR Comput.Math.Phys.，1961.

［14］Bachvalov N S.On the convergence of a relaxation method with natural constraints on the elliptic operator［J］.USSR Comput.Math.Phys.

［15］徐小文.可扩展并行代数多重网格算法研究［D］.中国工程物理研究院，2007.

［16］Falgout R D.An introduction to algebraic multigrid［R］.UCRL-JRNL-220851.

［17］Ashby S F，Falgout R D.A parallel multigrid preconditioned conjugate gradient algorithm for groundwater flow simulations［R］.UCRL-JC-122359，1995.

［18］Ruge J W，Stuben K.Algebraic multigrid［M］.In Multigrid Methods，Frontiers Appl.Math.3，McCormick S F editor，SIAM，PhiladelPhia，1987:73-130.