固体火箭发动机水下超音速射流数值研究

王利利， 刘影， 李达钦， 吴钦， 王国玉

(北京理工大学 机械与车辆学院, 北京 100081)

0 引言

潜射导弹具有隐蔽性好、系统简便灵活、体积小、质量轻等优点，作为可靠的二次核打击力量，越来越受到各国的重视[1]。固体火箭发动机由于结构简单、可靠性高、快速反应能力强等优点，逐渐成为潜射导弹的首选动力源[2]。然而，由于发动机水下工作环境介质密度远大于空气，且处于高温、高压、超音速工作状态，其射流过程十分复杂[3]：一方面，大量的高温燃气喷入水中，会在导弹尾部形成瞬态燃气泡，高温燃气泡与水之间存在气- 水掺混、换热和相变等复杂物理现象并伴有激烈的流体动力干扰[4]；另一方面，燃气泡中存在的复杂激波结构[5-6]以及激波结构的运动[7]均会使流场压力剧烈变化，从而导致推力性能难以预测。近年来，为了提高发动机动力系统的推力性能，其采用的喷管扩张比越来越大，导致其稳定性大大降低。因此，深入研究大扩张比水下发动机的工作特性及其相关两相流动机理，对水下发射技术发展具有重要意义。

国内外学者对水下超音速气体射流过程开展了大量研究。王柏懿等[8]、Shi等[9]、施红辉等[10-11]通过试验对扩张比为1.2、马赫数为1～2的水下气体射流过程进行了观察研究，并利用探针排测量了尾流场及喷管壁面上的压力脉动情况，试验结果表明，气体中出现的大幅值压力脉冲与在喷管出口两倍直径处的射流形貌存在关联性，喷口端面处的回击现象与流场气相区中的压力脉动相关，射流的胀鼓与回击均可引起压力脉动,回击频率为胀鼓频率的2～3倍。Fronzeo等[12]针对扩张比为1.7的喷管在常温下不同密度液相环境下的射流过程进行了数值研究，研究发现液体密度是制约射流贯穿距的主要因素，气- 液界面的湍流水平随着环境液体密度的增大而急剧增加，导致混合层变化剧烈、激波胞室结构较早破裂。在此基础上， Li等[13]结合试验与数值计算方法，对含有蒸汽的气体射流提出了一种预测射流贯穿距动量平衡模型。李婷婷等[14]用高速摄像机观测了不同流量下环形喷管的射流夹断过程，结果表明，颈缩导致气体内部压力突增，从而引发射流夹断现象，夹断的频率随着气体流量的增加而降低。Wilson等[15]对1 800 K气体、马赫数分别为1和2的水下超音速射流进行了试验研究，分别考察了圆形、矩形、椭圆形和三角形喷嘴的射流情况，研究表明在使用椭圆形及矩形喷管时，气体射流的传播速度及气- 水之间的混合要比使用圆形喷管时快，而与周围低温水环境之间强烈的混合作用会使射流总温降得较低。由于试验中所用的射流入口压强较低(111～237 Pa)，且没有对速度、压力等参数进行测量，无法充分表征发动机水下工作特性。

水下超音速气体射流过程的推力特性及影响因素一直是研究的重点与难点。王宝寿等[16]通过试验测量了喷管扩张比为3、滞止压强为9 MPa、气体初始温度为2 700 K、水深为10～40 m条件下固体火箭发动机水下工作时的推力和侧向力,获得了推力、侧向力、推力矢量角随发动机喷口堵塞面积比、水深的变化规律。张有为[3]利用球形气泡模型对发动机水下推力特征的影响因素进行了研究，结果表明，推力峰值随着燃气泡初始半径的增大而迅速降低，随着燃面面积的减小而降低。

已有研究表明，水下超音速射流往往发生周期性的变化，表现为形态上的周期性变化，压力场、速度场等出现循环振荡，该现象对推力将产生较大影响。然而现有研究大多关注流场整体特征的变化，且对大扩张比喷管在超高温、超高压工作条件下的超音速水下射流研究较少。本文采用计算流体力学(CFD)方法，针对高温高压下两种扩张比喷管的射流开展模拟研究，揭示射流流场变化与推力变化机理,以获得扩张比对流场结构与推力的影响规律,所得研究成果对研究水下气体射流与固体火箭发动机水下推进具有重要意义。

1 数值计算方法

1.1 流体体积方法及其控制方程

将气体假设为理想气体，水为不可压缩相，采用流体体积方法对气- 液运动界面进行追踪。连续方程、动量方程、能量方程、体积输运方程和状态方程分别如下：

(1)

(2)

式中：ρm为汽- 液混合密度；μm为汽- 液混合黏性系数；ρm、μm可以统一用物质属性φm来表示，φm=φgαg+φwαw，φg、φw分别表示气相密度或黏性系数、液相密度或黏性系数，αg、αw分别表示气相体积分数、液相体积分数，ui、uj为速度分量；p为压强；Fi为体积力，忽略重浮力作用下Fi=0.

(3)

式中：气- 液两相体积分数关系为αg+αw=1；能量En和温度Tt为各项的质量平均变量；γeff为有效热传导率；Sh为源项。

(4)

(5)

式中：R为理想气体常数。

1.2 基于滤波器的修正湍流模型

本文采用基于滤波器的修正湍流模型[17]，在标准k-ε湍流模型中加入滤波函数，对尺度小于滤波器尺寸的湍流采用标准k-ε模型(见(6)式)；对尺度大于滤波器尺寸的湍流结构则采用直接计算方法求解(见(7)式)。

(6)

(7)

式中：ρ为流体密度；k为湍动能；ε为动能耗散率；Pk为湍动能生成项；模型常数σk=1.0，Cε1=1.44，Cε2=1.92，σε=1.30；μt为湍流黏性系数，

(8)

模型常数Cμ=0.09，F为滤波函数，由滤波器尺寸λ和湍流长度比值的大小决定，

(9)

式中：模型常数C3=1.0.

由(9)式可知，当湍流尺度较大时，湍流黏性系数表达为

(10)

1.3 计算模型及边界条件

图1 边界条件及几何示意Fig.1 Boundary conditions and geometry

图2 入口无量纲压力与入口温度函数曲线Fig.2 Curves of inlet pressure and temperature

根据流动物理特征分析，在壁面、喷管喉部及射流核心区等物理量变化梯度较大的区域进行网格加密，其局部放大图如图3所示，计算时间步长Δt=10-6s.

图3 网格划分结果及其放大图Fig.3 Mesh generation and its enlarged view

扩张比η为出口截面面积与喉部截面面积之比，本文以扩张比η=14与η=3.4的喷管为例分析流场和推力变化规律，二者受缩段、喉部段尺寸一致，控制出口直径de以获得不同的扩张比。工作参数如表1所示，其中Ma为设计工况出口马赫数。

表1 喷管工作参数

1.4 数值结果验证

本文根据文献[9]试验结果进行数值方法验证，喷管结构尺寸与工作条件均与文献[9]一致。

表2为射流发展过程形貌对比。从表2中可以看出，数值模拟结果与试验结果呈现出相同的周期性变化特征，数值计算结果对颈缩位置、胀鼓泡和回击流的捕捉均与文献[9]较为一致。

表2 射流发展形貌的数值模拟结果与试验结果[9]对比

Tab.2 Comparison of numerically simulated and experimental results of jet pattern

图4 尾流场静压均方根数值模拟结果与试验结果[9] 对比Fig.4 Comparison of numerically simulated and experimental results[9] about flow field static pressure RMS

1.5 推力计算方法

固体火箭发动机水下所受推力示意图如图5所示，其中As为发动机外壳横截面积，pb为发动机尾部背压，ve为出口排气速度，me为出口排气质量。由图5可见，固体火箭发动机在水下工作时，由于整个工作过程为非稳态湍流过程，发动机外壁所受到的外部介质静压作用非均匀分布，且发动机尾部所受到的背压不等于环境压强，即pb≠pa. 因此，固体水下火箭发动机的推力不可按空气中的简化推力公式[18]来计算。

图5 推力示意图Fig.5 Schematic diagram of thrust

考虑背压pb后的发动机推力计算公式为

(11)

由(9)式可知，推力由3部分组成：

2) 压差推力Ae(pe-pa)，由喷管排气面上的压强pe与外界压强pa之间的不平衡产生，取决于排气截面的燃气参数pe、Ae和发动机所处的环境压强pa.

3) 压差推力(As-Ae)(pb-pa)，由作用于发动机前后端面上的环境压强pa与背压pb之间的压力差产生。

由于不同发动机的外壁尺寸As差异较大，As取不同值时对推力公式中的最后一项影响较大，故省略(11)式最后一项，化为积分形式：

(12)

式中：re为喷管出口半径；ρe为出口截面上流体的密度；rs为出口截面半径。

2 结果与讨论

2.1 不同扩张比喷管射流瞬态演化过程

表3给出了不同扩张比喷管在初始膨胀阶段和充分发展阶段的典型水气云图与压力云图，p*=(p-pa)/0.5Aeρrefu2，以入口中心处为坐标零点，横坐标x*=x/de, 纵坐标y*=y/de.

表3 射流发展过程云图

Tab.3 Water-vapor contour of jet evolution process

由表3可见：

初始膨胀时期，在气囊推开液相作用下，气- 液界面压力略有升高(如虚线框所示)；由于气囊处于膨胀阶段，一部分流线辐射向外，另一部分流线在靠近管口处形成回流。

在充分发展时期，压力波动集中在喷管出口外(如虚线框所示)，随着射流向下游发展，压力逐渐与环境一致；射流下游形成了大尺寸旋涡，这是因为膨胀气体随主射流通道向下游喷射，造成水与气体裹挟与掺混。对比两种喷管射流形态可知，η=14.0的喷管射流轴向贯穿距远大于3.4的喷管，这是因为扩张段的大扩张比增加了气体膨胀惯性。

为了进一步研究射流发展的瞬态演变规律，根据充分发展时期射流特征周期性变化特点，将其分为颈缩、胀鼓和回击3个阶段，表4、表5给出了两种喷管在各阶段的压力云图与水气轮廓线、流线图与速度云图、涡流结构判别标准Q云图，Q以其最大值进行无量纲化处理。Q定义为Q=1/2(|Ω|2-|S|2)，其中Ω为涡度张量，S为应变张量变化率。在颈缩阶段，气流在垂直射流的方向上被压缩，射流通道变窄，由于主射流通道阻塞，喷管口形成了较大范围的高压区，导致喷管壁压力突增。喷管口流线呈辐射状，未出现气- 液倒灌现象，表明高压高温条件有效遏制了流动分离现象。在胀鼓阶段，气体阻塞在颈缩点上游，发生膨胀，且由于过膨胀气体推开水的能力较弱，胀鼓位置紧贴喷管口外壁；随着胀鼓泡内压力能的释放，胀鼓泡膨胀，出口速度增加，喷管内部高压区消失；在胀鼓泡内部，围绕主射流通道形成较高速度的回流，出现了Q的最高值，如表4、表5中胀鼓阶段Q云图所示。在回击阶段，胀鼓泡内累积能量与液相能量平衡最终冲破颈缩点，射流通道打开，形成了稳定的高速射流，此时出口速度达到设计马赫数。由表4中回击阶段Q云图可以看到，形成了大范围的回流，回流气体来自胀鼓泡内被主射流带动形成负速度的气体。同时，喷管口附近形成低压区，如表4、表5中压力云图所示。由于回流对主射流通道有剪切作用、高压水环境对射流通道的挤压作用，射流通道在喷管口附近再次发生颈缩，上述过程将周期性地重复。

表4η=14.0时射流发展过程云图各阶段水气分布图、压力云图、Q云图(局部放大)

Tab.4 Partially enlarged contours of water vapor distribution, pressure andQforη=14.0

表5 η=3.4时射流发展过程云图各阶段水气分布图、压力云图、Q云图(局部放大)

对比扩张比η=14.0与η=3.4的喷管射流形貌可知：η=3.4的喷管在水气轮廓图上未发生较大变化，颈缩阶段没有出现明显的主射流通道变窄现象，这是因为扩张比小，气体离开喷管后延径向扩张的惯性较小，故仍在轴向保持较大的动能与压力能；η=3.4的喷管出口维持稳定的超音速射流，除了在颈缩点出现小范围高压区，出口一直维持高速低压状态，而η=14.0喷管在颈缩阶段出口速度为亚音速；η=3.4喷管出口处没有观察到紊乱的流线，且该处气体径向膨胀的趋势小，故流线未呈辐射状。胀鼓阶段在水气轮廓图上并未出现明显胀鼓泡，出口速度增加，由高压转为低压，回流产生在液相区。在回击阶段出口速度增加，出现大范围低压区；从流线图上可以看到，回流并未直接回击到管壁，表明该喷管回击带来的压力振荡与腐蚀较小。高速气体射流气- 液边界的相对运动较稳定，气- 液掺混较为缓和。

2.2 不同扩张比喷管的推力特性

为了进一步研究不同扩张比喷管推力特性，给出了喷管出口中心处监测点A(见图6)的瞬态压力变化，如图7所示，其中横坐标为时间无量纲化结果t*=t/T. 由图7可知：η=14.0时，在初始膨胀时期呈现高频随机性震荡，在充分发展时期，由于激波推入管内，发生明显的脉冲振荡，平均周期为0.010 12 s；对比发现，η=3.4的压力曲线基本呈现一条直线，未出现脉冲，这是因为扩张比较小，出口压力较大，形成了稳定射流，激波一直未进入喷管内。

图6 监测点A位置Fig.6 Location of Point A

图7 监测点A压力变化曲线Fig.7 Pressure at Point A

图8 单个周期内动量推力与压差推力曲线Fig.8 Momentum thrust and differential pressure thrust curves in a single cycle

图9 喷管口附近压力云图Fig.9 Pressure contours near nozzle

2.3 不同扩张比喷管背压振荡特性

作为推力的影响参数，背压振荡与推力振荡有相同的变化规律。为了研究背压振荡规律，在出口侧壁上布置4个监测点B1、B2、B3、B4如图10所示，每个监测点相距10 mm. 图11和图12分别给出了B1点压力振荡曲线及相应快速傅里叶变换频谱图；表6给出了各点的主频。从表6中可知，各点主频相近，η=14.0的喷管均在100 Hz左右，η=3.4的喷管均在200 Hz左右。由于高频噪声在水下衰减更快，η=3.4的喷管在水下推进时隐蔽性更好。

图10 监测点位置及放大图Fig.10 Locations of monitoring points

图11 监测点B1压力曲线Fig.11 Pressure curve of B1

图12 B1点频谱图Fig.12 Spectrogram of B1

表7为两种扩张比喷管的推力主要参数。从表7中可以看到：两个扩张比喷管均呈现周期性特征；η=14.0的喷管振荡频率低于η=3.4的喷管，后者周期短，压力振荡频率高；η=14.0的喷管呈现明显的胀鼓现象，出现大范围的高压区；η=3.4的喷管周期性特点较弱，胀鼓程度小，因此脉动幅值小，射流稳定；η=3.4的喷管负推力小，平均推力大。

表6 B1、B2、B3、B4点主频

表7 η=3.4与η=14.0推力主要参数

3 结论

本文对比了两种扩张比喷管的流场结构与推力变化，并分析了二者之间的关系。得到结论如下：

1) 固体火箭发动机水下射流分为初始膨胀时期和周期性发展时期，其中周期性发展时期包括颈缩、胀鼓、回击3个阶段：颈缩阶段，激波由喷管外向管内移动，喷管口处压力在颈缩末期达到峰值，此时动量推力最小；胀鼓阶段，激波向外移动，喷口速度增加，动量推力上升，压差推力下降； 出口处达到稳定的高速射流；回击阶段，出口压差推力不变，动量推力小幅度震荡。

2) 不同扩张比的喷管射流形态不同。η=3.4的喷管各阶段射流气体形貌变化较弱，η=14.0的喷管各阶段形貌特性明显；η=3.4的喷管出口维持着稳定的超音速射流，而η=14.0的喷管在颈缩阶段出口速度为亚音速；η=3.4的喷管高速气体射流气液边界的相对运动较稳定，气- 液掺混较为缓和。

3) 不同扩张比的喷管推力特性不同。η=3.4的喷管胀鼓程度小，因此脉动幅值小，射流稳定；η=3.4的喷管负推力小，平均推力大；η=3.4的喷管荡频率高于η=14.0的喷管，因此工作噪声衰减更快、隐蔽性更好。故η=3.4的喷管比η=14.0的喷管更有利于水下推进。