李映坤，韩珺礼，2，陈 雄，鞠玉涛

(1.南京理工大学 机械工程学院，南京 210094;2.北京机电研究所，北京 100012)

0 引言

固体火箭发动机具有结构简单、使用方便、可靠性好等优点，在国防领域受到广泛重视及应用。与传统的单室单推、单室双推发动机相比，双脉冲固体火箭发动机具有改善固体火箭发动机能量可控性的突出优势，可多次点火，并提供不连续推力，有利于提高导弹武器的作战能力［1］。

级间隔离装置是双脉冲固体火箭发动机的关键部件之一。根据隔离装置承力情况，隔离装置可分为隔板式(硬隔离)和隔层式(软隔离)。其中，隔板式隔离装置又包括了隔塞式隔板、陶瓷舱盖隔板、金属膜片式隔板等，它们都具有结构简单、易于加工、承压以及密闭可靠性好等优点，国内外对类似结构的双脉冲发动机已进行了试验研究［2－3］，但试验测得的推力-时间曲线与计算结果有一定的偏差，难以重现内弹道的性能。这是因为在二级脉冲燃烧室工作期间，一级脉冲燃烧室中存在多个隔塞(或者陶瓷碎片)，隔塞从喷管排出的过程中对喷管的内流场产生了影响，造成发动机瞬时推力的波动，而国内外未见对隔塞在喷管中运动过程仿真的报道。

本文基于AUSM-PW矢通量分裂格式，求解圆柱坐标系下的二维轴对称非定常可压缩雷诺平均Navier-Stokes方程，采用基于格心的有限体积法及双时间步LU-SGS方法，结合运动嵌套网格技术，研究了隔塞沿着喷管轴线运动过程中喷管内流场的演变过程，并分析了隔塞在喷管内运动对发动机推力的影响，为相关设计提供了参考。

1 控制方程及数值求解方法

1.1 流动控制方程

采用ALE(Arbitrary Lagrangian Eulerian)有限体积方法描述的可压缩非定常Navier-Stokes方程为

式中 U为守恒变量的一般形式;Fc为无粘通量;Fv为粘性通量;H为轴对称几何源项;Fc=Fi+Gj;Fv=Fvi+Gvj;(i，j)表示直角坐标系2个坐标方向(x，y)的单位矢量。

式中 ρ为气体密度;u、v为气体运动速度矢量的2个分量;xt、yt分别为网格运动速度矢量的2个分量;T为气体的温度;E为单位体积气体的总能量;τ表示应力张量，其具体形式参考文献［4］。

对于湍流问题，本文采用Menter F R提出的k-ε剪切应力输运(shear-stress-transport)模式［5］，该模型通过混合函数F1将k-ε模型和k-ε模型结合起来，这样充分发挥了k-ε模型对自由流和k-ω模型对壁面受限流动的处理优势。具体描述如下:

式中 k为湍动能;ω为比耗散率;μt为湍流粘性系数，其他它参数的具体形式见参考文献［5］。

1.2 数值求解方法

本文采用基于格心的有限体积法，对方程进行离散。其中，无粘项采用AUSM-PW失通量分裂格式。以i方向上的通量F为例，i+1/2边界上的通量可写为［6］

式中 c为单元界面声速;Φ为守恒通量;p为压力项。

粘性项采用Jameson中心差分法离散，为了提高非定常流动的时间计算精度，同时又要具有较高的计算效率，本文采用Jameson提出的一种双时间步［7］的计算方法，在冻结的真实时刻点上，引入类似牛顿迭代的虚拟时间迭代过程，通过这种内迭代过程来提高LU-SGS隐式算法处理过程中损失的时间精度。

2 嵌套网格技术

隔塞在喷管内运动的过程为典型的移动边界问题。目前，处理移动边界问题较成熟的方法包括非结构动网格技术及动态嵌套网格技术等，但非结构网格在处理具有较大相对位移问题时，网格变形与重构的表现能力较差。因此，本文选择动态嵌套网格方法，其优点是对于复杂几何体，不要求流场各个计算域共享边界，减轻了网格生成难度，适用于运动幅度较大的相对运动。

2.1 洞边界的确定

在嵌套网格系统中，建立人工插值内边界的过程称之为挖洞过程。为了提高确定洞边界的效率，本文采用 Chiu I T和 Meakin R L提出的“Hole-Map”方法［8］，其核心思想是对于给定的嵌套网格体系，若已知其拓扑结构，就能用均匀的笛卡尔网格单元去近似每个曲面边界，从而得到该曲面边界的笛卡尔网格近似。在此基础上，网格点P与曲面边界的位置关系可近似转换为网格点P与笛卡尔网格的关系，而这种关系十分容易确定。

2.2 贡献单元的搜寻

在嵌套网格系统中，通过贡献单元将流场解的信息插值到插值边界面网格上。因此，寻找贡献单元技术的优劣，对整个嵌套网格方法有很大影响，有时甚至是嵌套网格方法成败的关键。本文结合Wang Z J的矢量判别法［9］，提出了一种简单易行的寻找插值点贡献单元的方法，具体步骤如下:

(1)对于任意给定网格单元，寻找距离其最近的网格单元;

(2)采用Wang Z J的矢量判别法，寻找包围该点的贡献单元。对于任意网格单元，如果其所有边界均满足rfc·n≥0，则点C位于网格单元内，从而该网格单元就为点C的贡献单元，其中f为网格单元边界的中心，n为边界外法线矢量;

(3)使用双线性插值方法进行插值，完成该点物理信息的传递。

3 算例验证

3.1 静态NACA0012翼型绕流

为了验证本文静态嵌套网格计算程序，根据文献［10］，计算了攻角 α0=2.05°，来流马赫数 Ma∞=0.755，雷诺数为 9.9×106工况下 NACA0012 翼型绕流。为了便于对比分析，本文分别采用了单块网格和嵌套网格。图1所示的是采用单块网格和嵌套网格计算的翼型表面压力系数分布，并与实验结果进行了对比。从图1中可看出，本文的计算与实验值吻合较好，而且嵌套网格计算的结果与单块网格计算的结果相差不大。因此，证明本文所编制的计算程序可信、可靠。

图1 翼型表面压力系数对比Fig.1 Comparison of steady pressuredistributions for the airfoil

3.2 俯仰振动的NACA0012翼型绕流

以NACA0012翼型绕1/4弦点作简谐俯仰振动为算例，验证本文动态嵌套网格计算程序的准确性，来流马赫数为0.755，翼型攻角的变化规律是

式中 初始攻角(平均攻角)α=0.016°;振幅 αm=2.51°;无量纲角频率 k=wc/u∞;c表示弦长;ω 是角频率。

图2给出了翼型的升力系数随攻角的变化规律，并与Batina J T［11］的结果以及实验结果进行了比较。从图2中可看出，本文的计算结果与文献［11］的计算结果吻合很好，但都与实验值稍有差异。Batina J T认为导致这一差异的原因是实验数据可能在稍大的平均攻角下获得的。

图2 升力系数随攻角的变化Fig.2 Lift coefficient vs angle of attack

图3 计算区域网格(上)和边界条件(下)Fig.3 Computational domain(upper)and prescribed boundary conditions(lower)

4 隔塞运动对喷管内流场的影响及结果分析

4.1 物理模型和边界条件

不考虑发动机工作过程中喷管的形变，计算区域取喷管和一段圆柱形外流场区域。喷管入口直径为220 mm，喷管喉部直径为72 mm，扩张比为2.5，收敛半角和扩张半角分别为30°和15°，一般隔塞式脉冲发动机多使用多个小隔塞，为了重点研究小隔塞在喷管内运动对发动机性能影响的规律，并简化仿真模型，本文研究了一个小隔塞。在喷管内运动的过程，隔塞的直径为30 mm，长度为20 mm。外流场区域径向取5倍的喷管出口直径，轴向取10倍的喷管出口直径。喷管入口的质量流率为18 kg/s，总温为3 300 K，总压为10 MPa。由于喷管入口为亚音速，因此在计算中采用局部准一维特征分析确定入口边界条件，壁面采用无滑移边界条件，绝热壁假设，计算域出口边界条件根据马赫数判定。当出口为超声速时，此时所有物理量外推;当出口为亚声速时，给定环境反压，其它参数由内向外插值。

首先，计算定常状态下喷管的流场;然后，让隔塞从初始位置(此时记时间t=0)沿着轴线运动。根据文献［12］的实验结果，隔塞从隔板飞出后，在一级燃烧室中的运动方向与燃烧室的轴线基本保持平行，隔塞到达喷管后，靠近中心线的一小部分隔塞直接排出喷管［12］，且实验得到隔塞到达喷管之前的飞行速度约为30 m/s。本文重点研究隔塞在喷管内运动过程中喷管内流场结构的演化，暂不考虑隔塞在喷管中的运动规律对流场结构的影响。因此，根据文献［12］的实验，假设隔塞以30m/s的速度在喷管中运动。计算区域网格如图3所示，包围隔塞的贴体网格整体在喷管的网格上运动，总网格数目为49 654。

图4 不同时刻的马赫数云图Fig.4 Mach number contour at different times

图 5 马赫数云图(局部)t=3.1 ms(上)t=6.0 ms(下)Fig.5 Mach number contour(local)at 3.1 ms(upper)and 6.0 ms(lower)

图6 t=3.1 ms时刻轴线上有无隔塞时马赫数分布Fig.6 Axial mach number distributions with plug and without plug at 3.1 ms

4.2 喷管内流场结构演化分析

以马赫数云图来揭示隔塞运动过程中喷管内流场的演变过程，并与不含隔塞的云图进行了对比。图4所示的是隔塞在喷管内运动的过程中，6个时刻的马赫数云图。图4中，上半部分为喷管中含隔塞时不同时刻的马赫数云图;下半部分为喷管中不含隔塞的云图。由图4可见，图4(a)、图4(b)中，隔塞位于喷管的收敛段内，距离喷管的喉部还有一段距离，其流场结构与不含隔塞时的结构相差不大。

在t=3.1 ms时刻(图4(c)所示)，隔塞正好位于喷管的喉部。此时，隔塞前端的燃气流为亚音速，后端的燃气流为超音速，附着在隔塞上的燃气流离开物面，流动出现分离，并在分离区上游引起分离激波;同时，燃气流在分离点后形成涡流，由于此时隔塞附近的流场为轴对称的，因此隔塞后端分离形成的涡流交汇形成尾迹，很大程度上影响了燃气流在喷管扩张段的流场结构，如图5上半部分所示。轴线上有无隔塞时，燃气流马赫数对比如图6所示。由图6可见，由于隔塞的影响，喷管轴线上的马赫数均比无隔塞的时候要小，尤其是在喷管的扩张段，平均要低0.5;同时，隔塞后端的马赫数先增大后减小，有一段很小的波动，这是因为燃气在隔塞的后端形成的涡流所引起的。

图4(d)、图4(e)、图4(f)中，隔塞已经完全通过了喷管喉部，位于喷管的扩张段中，此时喷管扩张段里的燃气流为超音速。因此，在隔塞的前端形成了弓形脱体激波，如图5下半部分所示。随着隔塞的后移，隔塞后端形成的尾迹与喷管出口的马赫锥逐渐相互作用，形成了另一道激波，如图4(e)、4(f)所示。同时，当隔塞位于喷管扩张段上游时，形成的弓形脱体激波还作用在喷管壁面上，造成此处压力的变化，如图7所示。由图7可见，与初始时刻相比，隔塞在喷管内运动的过程中，作用在喷管扩张段壁面上的压力剧烈波动。

图7 不同时刻喷管壁面上的压力Fig.7 Pressure distributions in nozzle wall at different time

4.3 隔塞运动对发动机推力的影响

火箭发动机的推力定义为静止条件下发动机所受内、外壁面压强的合力，发动机内壁面压强的合力由两部分组成，即

式中 FM为作用在燃烧室壳体上的力，FN为作用在喷管上的力。

每种力都可由作用在壁面上的压强和粘性力积分得到，但粘性力仅占到总推力的1%以下。因此，计算中忽略粘性力，同时还假设作用在燃烧室壳体上的力FM和发动机外壁面受到的力保持不变。所以，在隔塞运动的过程中，只有FM变化。

图8中所示的为发动机无量纲推力随时间的变化曲线(取火箭运动的方向为正方向)。从图8中可看出，在隔塞沿着喷管轴线运动的过程中，发动机的瞬时推力先逐渐减小，在t=3.1 ms时刻降到最小，并在接下来的一段时间内，瞬时推力的大小变化不显著，t=3.5 ms之后，瞬时推力开始急剧升高，并在t=4.5 ms时刻推力达到最大值，随后缓慢下降，最终趋于初始时刻的值。隔塞的运动是造成发动机瞬时推力变化的主要原因，隔塞通过喷管喉部后，由于喷管扩张段中燃气流为超音速，因此在隔塞前端形成了弓形脱体激波，并与喷管壁面相互作用，造成喷管壁面上压力的变化，如图7所示。t=3.1 ms时刻，作用在喷管扩张段内的压力急剧减小，而此时收敛段内的压力变化不大，积分后的合力减小。所以，此时瞬时推力降至最小。t=4.5 ms时刻，作用在喷管扩张段内的压强显著增大，积分后的瞬时推力也达到最大。

图8 无量纲推力随时间的变化曲线Fig.8 Normalized thrust vs time

同时，图8中还显示隔塞在喷管收敛段运动时，作用在喷管壁面上的力变化不大，发动机的瞬时推力几乎不变;隔塞位于喷管喉部时，发动机瞬时推力降至最小，相比于初始时刻减小了27.7%;隔塞完全通过喷管喉部后，并位于喷管扩张段上游时，发动机瞬时推力达到最大值，比初始时刻增大了8.9%。一般隔塞式双脉冲发动机二脉冲工作时，一级燃烧室中的隔塞不止一个。因此，当这些隔塞依次排出喷管时，会使发动机瞬时推力发生巨大改变，造成推力-时间曲线的波动，使发动机的内弹道性能难以重现，严重影响发动机的正常工作。因此，有必要针对隔塞的形状及运动规律进行大量数值计算，对其特性进行详细研究，以总结出规律。

5 结论

(1)本文所采用的动态结构嵌套网格方法能处理运动幅度较大的相对运动问题，具有较高的工程应用价值，为数值模拟隔塞在一级燃烧室中的运动轨迹奠定基础。

(2)隔塞通过喷管喉部后，燃气流在隔塞前端形成弓形脱体激波，在隔塞后端形成尾迹，作用在喷管收敛段上的压强变化显著。

(3)第二级燃烧室中的多个隔塞依次从喷管排出的过程中，发动机的瞬时推力逐渐增大，到达最大值后保持一段很小的时间，然后开始下降，最终缓慢增加至初始时刻的值。

(4)发动机正常工作时，当直径为30 mm的隔塞运动到直径为72 mm喷管喉部时，发动机瞬时推力降至最小，比初始时刻减小了27.7%;隔塞完全通过喷管喉部后，并位于喷管扩张段上游时，发动机瞬时推力达到最大值，比初始时刻增大了8.9%。

［1］Naumann K W，Stadler L.Double-pulse solid rocket motor technology-applications and technical solutions［R］.AIAA 2010-6754.

［2］Nishii S，Fukuda K，Kubota N.Combustion tests of two-stage pulse rocket［R］.AIAA 89-2426.

［3］Wang Chang-hui，Liu Yu ，Liu Ya-bing.Design and experimental studies on ceramic port cover for dual pulse motor［J］.Acta.Astronautica，2011，68(11):1881-1890.

［4］张涵信，陈坚强，高树椿.H2/O2燃烧的超声速非平衡流动的数值模拟［J］.宇航学报，1994，15(2):14-23.

［5］Menter F R.Two-equation eddy-viscosity turbulence models for engineering applications［J］.AIAA Journal，1994，32(08):1598-1605.

［6］Liou Meng-sing.Methods for the accurate computations of hypersonic flows［J］.Journal of Computational Physics，2001，174(1):38-80.

［7］Jameson Antony.Time dependent calculations using multigrid，with applications to unsteady flows past airfoils and wings［R］.AIAA 91-1596.

［8］Chiu Ing-Tsau ，Robert L Meakin.On automating domain connectivity for overset grids［R］.AIAA 95-0854.

［9］Wang Z J，Parthasarathy V.A fully automated Chimera methodology for multiple moving body problems［J］.International Journal for Numerical Methods in Fluids，2000，33(7):919-938.

［10］John B McDevitt，Arthur F Okuno.Static and dynamic pressure measurements on a NACA0012 airfoil in the ames high Reynolds number facility［R］.NACA-TP-2485.

［11］John T Batina.Unsteady euler airfoil solutions using unstructured dynamic meshes［J］.AIAA Journal ，1990，28(9):1381-1388.

［12］李江，肖育民，何国强，等.双脉冲固体火箭发动机二次点火内视研究［J］.推进技术，1998，19(3):61-64.