1 引言

大多数液体火箭发动机在研制过程中, 都会碰到燃烧不稳定性的问题^{[1, 2]}。燃烧不稳定性按频率可以分为低频、中频和高频燃烧不稳定性^[1]。高频燃烧不稳定性是由燃烧放热和具有声学振型特征的压力振荡耦合造成的, 通常具有声学特征, 而又被称之为高频声学燃烧不稳定性, 因其会破坏壁面换热条件而导致发动机推力室瞬时烧毁或爆炸而具有极大的危害性。因此这种燃烧不稳定性得到了人们广泛和持续的关注和研究。对于一个确定构型的推力室而言, 人们和工程最为关心的是该构型最容易激发的高频声学振型及其构型参数的影响规律, 其对燃烧不稳定性的控制具有重要意义。

燃烧室的声学特性主要采用理论、实验和数值模拟等方法进行研究。在理论研究中, 求解基于小扰动假设建立的声学方程, 可获得燃烧室的声学振型频率^[2], 但由于燃烧室结构特殊性和内部物理过程的复杂性, 理论推导中存在很多与实际相差很大的假设。在冷态声学实验中, 通常在模型燃烧室壁面安装扬声器或二次喷管来激发燃烧室的压力振荡, 并对其进行FFT分析, 获得燃烧室声学振型频率及其衰减率等声学特性^[3~7], 但扬声器和二次喷管往往产生的是特定频率的扰动, 激发单一振型的压力振荡, 并且姿轨控发动机燃烧室的尺寸在厘米量级, 增大了传感器的安装和测量的难度, 实验中得到的数据有限。在数值模拟中, 可以获得丰富的数据, 可显示压力波的时空演化形成声学振型的过程。因此, 数值模拟成为燃烧室声学特性研究的首选。

在数值模拟中, 有基于线性和非线性声学方程的有限元方法（FEM）^[7~9]和基于N-S方程的CFD方法^{[10, 11]}。通常在燃烧室中引入质量^[12]、能量^[13]、压力^[14~18]或速度扰动^{[19, 20]}等外部激励来激发燃烧室的压力振荡, 进一步通过频谱分析方法获得燃烧室声学特性。这些研究的主要目的是对于给定声学振型, 确定其频率及其衰减特性；或者通过具有特定频率的小扰动来激发燃烧室自激高频声学燃烧不稳定。在这些施加激励的方式中, 质量扰动通常是通过燃烧室入口流量的振荡实现的, 而其他几种扰动方式则是通过燃烧室局部相关量的人为振荡实现的。在后几种扰动方式中, 人为引入不同形式的能量, 但却没有相应的质量引入来承载, 从而使其激励幅值相对较小, 否则会导致低密度或高温度等非物理异常, 甚至使计算崩溃。因此, 这些施加激励的方式难以用来确定给定构型燃烧室的最容易激发的声学振型。

要确定给定构型推力室最容易激发的声学振型, 则必须在推力室中施加非特定频率大幅值激励, 如一些实验研究中采用的定容弹激励^[1]。另外在液体火箭发动机中, 燃料和氧化剂均是通过头部面板上的喷嘴引入, 从而在头部会形成一些剧烈燃烧的空间离散点源, 其可能形成燃烧室的激励源而诱发声学振型燃烧不稳定性发生^{[2, 21, 22]}。为此, 本文提出了施加非特定频率强激励的“数值定容弹”模型, 并在给定构型燃烧室头部易形成剧烈燃烧的局部施加该数值定容弹激励, 获得了具有多种声学振型特征的压力振荡。基于该压力振荡确定出了最容易激发的声学振型, 并进一步研究了推力室收缩比对最容易激发声学振型的影响。为工程中抑制燃烧不稳定性发生提供依据和参考。

2 数学物理模型 2.1 基本方程

以自燃推进剂姿轨控发动机推力室为对象, 基于数值定容弹模型在冷态条件下开展了其声学特性研究。对真实工作条件下的推力室, 在喉部和扩张段分别达到声速和超声速状态, 使得喉部下游区域的扰动无法传到上游, 因此喷管扩张段及喉部以下区域的流动对燃烧室的声学特性没有影响, 喉部可看为声学闭合边界条件。鉴于此, 在本文推力室声学特性的数值模拟研究中, 只考虑头部面板、圆柱段、收敛段和喉部截面围成的区域, 如图 1所示。圆柱段长度和直径分别为7.1mm和32mm, 收敛段长度和收缩比分别为23.3mm和14。

基于k-ε双方程湍流模型, 描述燃烧室中流动的URANS控制方程如下^[21]

$ \frac{{\partial \rho }}{{\partial t}} + \frac{{\partial \left( {\rho {u_j}} \right)}}{{\partial {x_j}}} = 0 $

(1)

$ \frac{{\partial \left( {\rho {u_i}} \right)}}{{\partial t}} + \frac{{\partial \left( {\rho {u_i}{u_j}} \right)}}{{\partial {x_j}}} =-\frac{{\partial \rho }}{{\partial {x_i}}}-\frac{\partial }{{\partial {x_i}}}\left( {\frac{2}{3}\rho k} \right) + \frac{{\partial {\tau _{ij}}}}{{\partial {x_j}}} $

(2)

$ \frac{{\partial \left( {\rho e} \right)}}{{\partial t}} + \frac{{\partial \left( {\rho e{u_j}} \right)}}{{\partial {x_j}}} =-p\frac{{\partial {u_j}}}{{\partial {x_j}}}-\frac{\partial }{{\partial {x_j}}}\left( {-K\frac{{\partial T}}{{\partial {x_j}}} - \rho Dh\left( T \right)} \right) + \rho \varepsilon $

(3)

$ \frac{{\partial \left( {\rho k} \right)}}{{\partial t}} + \frac{{\partial \left( {\rho {u_j}k} \right)}}{{\partial {x_j}}} =-\frac{2}{3}\rho k\frac{{\partial {u_j}}}{{\partial {x_j}}} + {\tau _{ij}}\frac{{\partial {u_i}}}{{\partial {x_j}}} + \frac{\partial }{{\partial {x_j}}}\left( {\frac{{{\mu _{{\rm{eff}}}}}}{{{{Pr }_k}}}\frac{{\partial k}}{{\partial {x_j}}}} \right)-\rho \varepsilon $

(4)

$ \begin{array}{l} \frac{{\partial \left( {\rho \varepsilon } \right)}}{{\partial t}} + \frac{{\partial \left( {\rho {u_j}\varepsilon } \right)}}{{\partial {x_j}}} =- \left( {\frac{2}{3}{C_{\varepsilon 1}}- {C_{\varepsilon 3}}} \right)\rho \varepsilon \frac{{\partial {u_j}}}{{\partial {x_j}}} + \\ \frac{\partial }{{\partial {x_j}}}\left( {\frac{{{\mu _{{\rm{eff}}}}}}{{{{Pr }_\varepsilon }}}\frac{{\partial \varepsilon }}{{\partial {x_j}}}} \right) + \frac{\varepsilon }{k}\left[{{C_{\varepsilon 1}}{\tau _{ij}}\frac{{\partial {u_i}}}{{\partial {x_j}}}-{C_{\varepsilon 2}}\rho \varepsilon } \right] \end{array} $

(5)

式中ρ, e, p, T, k, ε分别代表燃烧室密度、内能、压力、温度、湍动能及其耗散率。u_j和x_j分别代表j方向的速度和位置坐标。μ_eff和τ_ij分别为有效粘性系数和雷诺应力。

该推力室放置在无限大静止环境之中, 其不工作。在预测推力室声学特性时, 通常不考虑喷注面板的局部细节并将之看成为声学封闭边界, 因此在本文推力室冷态声学特性研究中, 喷注面板处理为声学闭合边界。采用自主开发的三维非定常湍流两相流动与燃烧数值模拟程序, 对图 1所示构型内由数值定容弹所激发的流动与压力波传播过程进行了三维非定常模拟。网格尺寸为1mm。燃烧室喷注面板和侧壁面都采用无滑移和绝热壁面边界条件, 燃烧室出口采用连续出流边界条件, 给定出口压力为0.1MPa。对流项采用二阶迎风格式, 扩散项采用中心差分格式, 时间项采用一阶向前欧拉格式, 并采用SIMPLE算法^[23]来求解方程。

2.2 数值定容弹模型

在燃烧室声学特性实验研究中, 通过采用爆炸弹、脉冲枪、二次喷管等^[1]人为激励来激发压力振荡。在数值仿真中, 在燃烧室中引入压力^[14~18]、质量^[12]或能量扰动^[13]等方法来激发压力振荡。在这些扰动模型中, 通常需要预先设定激励频率或脉冲型扰动幅值不能太大^[14]。对一个压力脉冲扰动, 迅速增加了加载区的压力, 但由于加载时间很短来不及膨胀, 从而必然导致加载区温度同倍数增加。对液体火箭发动机推力室, 其温度通常高达3000K以上, 因此如在温度为3400K的区域中加载幅值为1.5倍室压的压力脉冲, 则加载区域的温度会超过5000K, 从而导致计算发散。这就是目前压力脉冲幅值较小或采用特定频率的压力扰动作为激励源的原因。由于衰减的原因, 较低压力脉冲幅值很难激发多模态声学振型, 从而难以确定给定推力室构型的最容易激发声学振型, 在热态燃烧情况下尤为困难。类比于实验中定容弹, 提出了数值定容弹模型, 即通过高压工质极短时间导入加载区域来实现压力脉冲。假定燃烧室任意大小的子区域Ω为压力脉冲加载区域, 加载前其压力、温度和密度分别为p₀, T₀和ρ₀, 某一瞬时开始在该区域导入压力、温度和密度分别为γp₀, αT₀和βρ₀的高压气体, 二者具有相同的组分且瞬时混合, 并以混合后的状态维持一段时间, 从而形成压力脉冲。基于加载前后加载区气体均满足状态方程以及加载是一个定容过程, 可以获得基于该数值定容弹模型所实现的压力脉冲的状态参数如下

$ \begin{array}{l} p\left( {x, y, z, t} \right) = {p_0}\left( {x, y, z} \right) + \gamma {p_0}\left( {x, y, z} \right)\\ \rho \left( {x, y, z, t} \right) = {\rho _0}\left( {x, y, z} \right) + \beta {\rho _0}\left( {x, y, z} \right)\\ T\left( {x, y, z, t} \right) = \frac{{1 + \alpha \beta }}{{1 + \beta }}{T_0}\left( {x, y, z} \right)\\ \gamma = \alpha \beta \\ x, y, z \in \mathit{\Omega} ;{t_0} \le t \le {t_0} + \Delta t \end{array} $

(6)

式中t₀为数值定容弹加载的起始时间, Δt为其实现的压力脉冲的脉宽。可见, γ为本模型所实现的压力脉冲的过压系数。和已有的通过仅仅改变加载区压力来实现脉冲压力的模型不同, 本模型中的压力脉冲是通过一定质量（βρ₀）的介质来承载的, 从而可以实现更大的脉冲压力。已有模型的过压系数大于1.5时, 就会出现过高温度而使计算发散, 本文提出的数值定容弹模型可以获得不小于19的过压系数, 为激发多模态的压力振荡奠定了基础。从而能够观察给定构型推力室中最容易激发的声学振型, 并确定各种影响因素的作用规律。

3 计算结果与讨论 3.1 计算结果验证

首先对基于本文提出的数值定容弹模型研究推力室声学特性的方法进行了验证。为此选取直径R和长度L分别为32mm和40mm, 且两端封闭的圆柱段作为算例进行了模拟, 其中充满压力和温度分别为0.1MPa和298K的静止空气。基于方程（6）所建立的数值定容弹模型, 在距离推力室头部面板3mm和圆柱1/2半径处的半径为2mm、长度为3mm的圆柱形区域施加γ和β均为19的数值定容弹, 持续0.1μs。在施加定容弹的径向反方向一侧近壁面位置设置观测点, 其压力随时间的振荡过程如图 2所示, 相应的FFT分析结果如图 3所示。可见, 在计算模型中激发出了显著的压力振荡, 且具有多种特征频率, 四个最强突频分别为4.6kHz, 6.7kHz, 8.1kHz和11.199kHz。

燃烧室横截面六个瞬时时刻的压力分布如图 4所示。图 4（a）中, 压力呈左右两侧高低压力分布, 压力波沿着圆周壁面向左侧的低压区传播, 形成两个近壁端高压区, 右侧高压区压力降低, 如图 4（b）所示。压力波继续向左传播到达左侧, 左侧变为高压区, 两个近壁端高压区变为低压区, 如图 4（c）所示。压力波从左右侧向两近壁端低压区传播, 形成左侧高压右侧低压的压力分布, 如图 4（d）所示。右传的压力波继续沿壁面向右侧低压区传播, 形成了两个近壁端的高压区, 左侧因压力波传出变成低压区, 如图 4（e）所示, 最后到达右侧, 右侧变成高压区, 形成左右两侧高低压力分布, 如图 4（f）所示。图 4（a）, 4（d）和4（f）是典型的一阶切向振型的压力分布, 其时间间隔约为1/2的一阶切向振型的周期。图 4（b）, 4（c）和4（e）是典型的二阶切向振型的压力分布, 其时间间隔约为1/2的二阶切向振型的周期。可见燃烧室中出现了一阶切向振型和二阶切向振型。

为确认图 3所示的压力振荡形态是否具有声学特征, 给出了模拟构型的理论声学特征频率^[2], 如下式所示

$ {f_{lmn}} = \frac{c}{{2\pi }}\sqrt {\frac{{\lambda _{mn}^2}}{{{R^2}}} + \frac{{{l^2}{\pi ^2}}}{{{L^2}}}}, l, m, n = 0, 1, 2, \cdots $

(7)

式中λ_mn为横向本征值^[2], l, m, n分别代表纵向振型阶数、切向振型阶数和径向振型阶数, c为声速, 取值为347m/s。该模拟构型的一阶切向（1T）、二阶切向（2T）、一阶纵向（1L）和一切-一纵（1T1L）声学振型对应的理论特征频率如表 1所示, 并和基于数值定容弹预测的结果进行了对比。可见, 预测的突频能够被该构型的声学振荡模态所辨识。图 3中四个峰值最大的突频依次可以辨识为一阶切向、二阶切向、一阶纵向和一切-一纵振型, 且预测的这些振型的特征频率与相应的理论声学振型频率相符。因此在模拟算例中激发的压力振荡是一个多模态的声学振荡, 且预测的声学振型频率与理论预测结果相符, 偏差在7%以内。表明本文建立的数值方法能够合理准确地模拟给定构型的声学特性。

3.2 推力室声学振荡特性分析

在图 1所示的推力室中, 充满压力和温度分别为0.1MPa和298K静止空气。实际推力室中, 在喷注面板附近的撞击对区域因为剧烈燃烧而容易出现压力扰动^{[2, 21, 22]}, 因此在位于燃烧室1/2半径并距离燃烧室头部3mm处的长度和半径分别为3mm和2mm的圆柱区域施加γ和β均为19的数值定容弹, 持续0.1μs。在施加定容弹的径向反方向一侧近壁面位置设置观测点, 其压力随时间的振荡过程如图 5所示, 相应的FFT分析结果如图 6所示。可见施加的数值定容弹激发出了显著的压力振荡, 而且具有多种特征频率, 分别为能量较大的7.3kHz和11.699kHz, 与能量较小的15.899kHz。

图 7给出了在局部施加定容弹激励之后, 依次八个瞬时推力室的压力分布, 能够展示出其压力振荡的时空演化过程, 红色代表高压, 蓝色代表低压。圆柱状定容弹施加在燃烧室1/2半径处, 并在此形成了高压区, 压力波以圆弧形式向外传播, 如图 7（a）所示。圆形压力波左半部分继续以圆弧型向左侧方向传播, 此时最左侧波面已接近推力室中心轴线区域, 而右半部分压力波在向右侧近端壁面传播, 其最右侧已经与壁面碰撞, 并沿近壁区域周向双向传播, 如图 7（b）所示。此时压力波到达右壁面和燃烧室轴心, 这些区域成为高压区, 而初始脉冲压力形成区由于压力波的传出而变成相对低压区。沿近壁区域周向传播的压力波和新近到达壁面的右传的初始压力波在近壁区汇聚, 并形成了如图 7（c）所示的近壁高压区, 而左传的初始压力波继续以圆弧状向左传播。随后汇聚在近壁区的高压区由于壁面反射向中心轴线区域传播, 而左传的初始压力波继续左传, 并整体同时到达远端壁面, 如图 7（d）所示。而后整体同时到达远端壁面的压力波聚焦到远端壁面与轴线之间的中间区域并与近端壁面反射过来的压力波交汇于此, 形成了一半区域高压, 一半区域低压的压力分布形态, 如图 7（e）所示, 具备了一阶切向振型的雏形。进一步演化出现了如图 7（f）~7（h）所示的典型的一阶切向和二阶切向的压力分布形态。

图 8为推力室压力分布的瞬态纵剖面图, 与图 7一致呈现出典型的切向振型的压力分布形态, 且在推力室整个纵向都呈现出了强弱程度不同的切向振型压力分布形态。

为了将图 6所显示的这些突频与图 7所展示的压力分布形态进行关联和一致性确认, 进一步基于推力室构型理论声学振型频率对这些突频进行了辨认。对图 1所示的推力室构型, 通常将其等效成两端封闭的圆柱来估算其理论声学频率。等效长度为燃烧室圆柱段长度加上二分之一的收敛段长度^[2]。采用公式（7）进行计算, 声速取为347m/s, 该构型的理论声学特征频率如表 2所示。将图 6的突频与理论声学振型频率相对照, 7.3kHz, 11.699kHz和15.899kHz分别可辨识为一阶切向、二阶切向和三阶切向, 可见与图 7和图 8中观察的形态是一致的。因此对如图 1给定构型的推力室, 基于本文提出的数值定容弹模型, 能够激发出多模态的声学振型。一阶切向、二阶切向和三阶切向的振幅分别为3.788kPa, 1.5871kPa和0.5659kPa。所以, 对如图 1给定的构型, 在头部偏心激励作用下, 最容易激发的声学振型是一阶切向振型。这是由于该激励在以圆弧型波面传播过程中受圆柱形壁面的碰撞、反射和聚焦, 从而容易形成周向传播的压力波或推力室周向高低压分区的分布形态, 从而易于激发出切向振型。切向振型是对发动机可靠稳定工作危害最大的一种燃烧不稳定形式, 也是工程最为关心的一种振型。

3.3 收缩比对燃烧室声学特性的影响

工程中, 在燃烧室直径和特征长度不变的情况下, 改变收敛段收缩比, 形成表 3所示的一组推力室构型, 其中, ζ, D, L_ch, L_cv和L^*分别为收缩比、圆柱段直径、圆柱段长度、收敛段长度和特征长度。收缩比为14的构型和上节的构型一样, 其反映了在特征长度不变的情况下推力室构型由细长到短粗的变化。为此, 针对这一组构型, 采用数值定容弹模型从最容易激发的声学振型角度研究了构型变化的影响。在如上节相同位置施加相同的激励, 并在相同的观测点记录了相应的压力振荡。

对四个工况均观察到了类似图 5的显著的压力振荡, 其频谱分析结果如图 9所示。基于这些构型的理论声学振型频率对图 9中的突频进行了辨识, 相应振型如图 9所标注。可见, 对表 3所示的一组推力室构型, 均激发出了多模态的声学振型的压力振动, 但所激发的声学振型及其振荡幅值随收缩比有明显的变化。随着收缩比的增加, 一阶纵向和一纵-一切振型削弱, 并在收缩比大于11之后没有出现；而一阶切向和二阶切向振型在各构型中均会出现, 并随收缩比的增加而不断增强。在圆柱段直径和特征长度不变的情况下, 较大的收缩比使推力室具有短粗的形状特征（如收缩比为14时, 圆柱段只有7.1mm）, 压力的纵向分布更加均匀和一致, 压力波主要通过周向传播来耗散, 所以形成了很强的切向声学振型, 而没有出现纵向声学振型；而较小的收缩比使推力室具有细长的构型特征, 不存在纵向相距很近的两个壁面使其间的压力分布更均匀的效应, 从而压力波在周向和纵向都可以充分的发展, 因此形成了切向振型、纵向振型和二者的复合振型。可见, 相对短粗的推力室构型, 细长的推力室构型更有利于抑制切向声学振型。

图 10和图 11给出了所激发的声学振型频率和幅值随收缩比的变化。由于推力室圆柱段直径没有改变, 所以收缩比对切向振型频率的影响很小, 基本保持不变。但由于圆柱段长度随收缩比有很大变化, 因此对纵向振型的频率影响显著。切向振型幅值随着收缩比增大而增大, 并且一阶切向和二阶切向振幅差也越来越大。而一阶纵向和一切-一纵振幅随收缩比增加而减小, 收缩比大于11时, 一阶纵向和一切-一纵振型消失。

4 结论

通过本文研究, 得到以下结论:

（1）在圆柱体中进行激励, 计算所得声学振型频率值与理论声学特征频率值相差在7%以内, 表明本文建立的数值方法是可靠的。

（2）在典型的姿轨控发动机燃烧室1/2半径处施加定容弹时, 出现了一阶切向、二阶切向和三阶切向振型的压力振荡, 其中一阶切向振型的压力振荡最剧烈, 振幅最大, 为最容易激发振型。

（3）在推力室特征长度和圆柱段直径不变的情况下, 随收缩比的增加, 纵向声学振型得到削弱, 甚至不再出现, 而切向声学振型得到增强。收缩比对切向声学振型的频率影响较小, 而对纵向声学振型的频率影响较大。和短粗推力室构型相比, 细长构型有利于抑制对发动机危害很大的切向声学振型。