1 引言

双模态超燃冲压发动机燃烧室内部流场结构复杂，存在激波串与附面层相互作用、流动分离再附等复杂的物理现象，为发动机的性能分析带来挑战。三维数值仿真虽然能够对这些复杂的物理现象提供详细的流场信息，但计算过程耗时耗力；而一维分析方法计算量小，可以通过合适的计算模型获得较高精度的计算结果^{[1, 2]}，因而被广泛应用于燃烧室性能评估^[3~6]。

目前，超燃燃烧室一维分析方法主要分为两类：(1)基于掺混和化学反应经验公式，获得燃烧室沿程释热分布，从而来计算其他一维流动参数。该方法需要采用合适的激波串模型，进而计算获得较为准确的压力、马赫数的分布^[7]。(2)基于实验获得的壁面压力分布和芯流面积拟合来计算其他一维流动参数。本文选用的是第二种方法，根据核心流再附的经验性假设，将超燃燃烧室内的流动分为芯流和分离区两部分，其中芯流反映了燃烧室内核心区流体的流动情况，可用于燃烧室的性能分析，其流动参数主要采用一维欧拉方法求解^[8]。从William Heiser等^[9]的一维分析结果可以看出，在燃烧室入口气体流速为Ma=2的工况下，回流区截面积与燃烧室几何截面积之比可高达25%。计算过程中如果直接将燃烧室截面积带入一维模型，会使流动分离再附过程的计算结果出现较大偏差^[10]，因此需要对芯流面积进行有效的评估。William Heiser等^[9]给出了芯流面积的经验性公式，但公式中的经验常数变化范围较大，难以根据实际情况给出合理数值；文献[9, 11, 12]根据等压燃烧假设，将燃烧室分为三段：绝热压缩段由气体动力学公式计算芯流面积，等压燃烧段通过差值或拟合的方式获得芯流面积，膨胀释热段由多变过程关系式“PAⁿ=常数”计算芯流面积，但发动机实际工作过程中，燃料的释热过程难以维持等压过程^[13]，多变关系式也难以满足任意工况的壁面压力数据^[9]；张鹏^[14]、程柳维^[15]等放弃等压燃烧假设，通过确定释热起始点和芯流再附点位置，通过拟合的方式获得芯流面积，但释热点与再附点位置是根据实验经验确定的，在一定程度上存在人的主观因素。

研究表明，为获得流动分离再附过程中的芯流面积，给出某个变量的变化规律是一个比较好的方法^[14]。因此，本文基于实验获得的壁面压力数据，引入燃烧室释热分布模型，对芯流面积进行求解。求解过程中，放弃了选取再附点的方法，认为燃料释热过程满足瑞利概率密度分布，释热模型中的经验常数由燃烧室构型参数和壁面压力反算获得，在一定程度上降低了人为干扰因素，提高了该模型的通用性。

2 物理模型和计算方法

在超燃冲压发动机燃烧室工作过程中，由于燃料喷注、释热等因素，导致燃烧室内的流动过程存在逆压力梯度，出现回流区，使流动分为芯流和分离区两部分，如图 1所示。在$ {{x}_{\text{sep}}}<x<{{x}_{\text{re}}}$段，芯流依次发生分离、加质、释热、再附过程。

2.1 控制方程

本文的研究基于定常流动假设，考虑芯流面积变化、摩擦、释热、加质的影响，构建了一维流动分析模型，控制方程组如下

状态方程

$ \frac{{{\rm{d}}\rho }}{\rho } + \frac{{{\rm{d}}T}}{T} = \frac{{{\rm{d}}p}}{p} $

(1)

连续性方程

$ \frac{{{\rm{d}}\rho }}{\rho } + \frac{{{\rm{d}}V}}{V} + \frac{{{\rm{d}}{A_{\rm{c}}}}}{{{A_{\rm{c}}}}} = \frac{{{\rm{d}}\dot m}}{{\dot m}} $

(2)

动量方程方程

$ kM{a^2}\frac{{{\rm{d}}V}}{V} = - E - F - \frac{{{\rm{d}}p}}{p} $

(3)

速度/马赫数方程

$ \frac{{2{\rm{d}}Ma}}{{Ma}} + \frac{{{\rm{d}}T}}{T} = \frac{{2{\rm{d}}V}}{V} $

(4)

总温方程

$ \frac{{{\rm{d}}{T_{\rm{t}}}}}{{{T_{\rm{t}}}}} = \frac{{{\rm{d}}T}}{T} + \frac{{\left( {k - 1} \right)M{a^2}}}{{1 + \frac{{\left( {k - 1} \right)}}{2}M{a^2}}}\frac{{{\rm{d}}p}}{p} $

(5)

式中p，$ \mathit{\dot m}$为壁面静压、质量流量，由实验测得；ρ，V，Ma，A_c，T，T_t，k分别为密度、速度、马赫数、芯流面积、静温、总温、比热比；$ E = \frac{{{C_{\rm{f}}}kM{a^2}}}{2}\frac{{{\rm{d}}{A_{\rm{w}}}}}{{{A_0}}} $，为摩擦项，C_f为壁面摩擦系数^[16]，A_w为湿表面积，A₀为燃烧室截面积；$ F = kM{a^2}\left( {1 - \bar y} \right)\frac{{{\rm{d}}\dot m}}{{\dot m}} $，为添质项，y为燃料在主流方向分速度与燃料喷注速度之比。

由公式(1)~(5)可得到Ma和芯流面积A_c的一维微分方程

$ \begin{array}{*{20}{l}} {\frac{{{\rm{d}}Ma}}{{Ma}} = - \left( {1 + \frac{{k - 1}}{2}M{a^2}} \right) \times }\\ {\left( {\frac{1}{{kM{a^2}}}\frac{{{\rm{d}}p}}{p} + \frac{1}{2}\frac{{{\rm{d}}{T_{\rm{t}}}}}{{{T_{\rm{t}}}}}} \right) - E' - F''} \end{array} $

(6)

$ \begin{array}{*{20}{l}} {\frac{{{\rm{d}}{A_{\rm{c}}}}}{{{A_{\rm{c}}}}} = \left( {\frac{1}{{kM{a^2}}} - \frac{1}{k}} \right)\frac{{{\rm{d}}p}}{p} + }\\ {\left( {1 + \frac{{k - 1}}{2}M{a^2}} \right)\frac{{{\rm{d}}{T_{\rm{t}}}}}{{{T_{\rm{t}}}}} - E'' - F''} \end{array} $

(7)

式中E'，E"为摩擦项，F'，F"为添质项。可以看出，马赫数是总温和壁面压力的函数，而芯流面积是马赫数、总温、壁面压力的函数。可以看出，在其他参数不变的情况下，释热越快，总温抬升越快，对应芯流面积增加得也就越快，因此选取合适的释热分布模型，即可通过四阶Runge-Kutta法计算得到燃烧室芯流面积变化情况。

2.2 释热分布模型

William Heiser等^[9]曾采用总温分布函数来描述燃烧室内的释热过程，但该模型的释热分布曲线为单调递减函数，释热最快的位置发生在释热起始点，这一结论与实际情况存在较大差异；近年来，本课题组^[17]通过瑞利概率密度分布来描述燃烧室内的释热分布情况，获得了较好的计算结果。图 2给出了瑞利概率密度分布曲线和Micka等^[18]通过实验获得的OH基团的分布。可以看出，二者的曲线具有较好的相似性，而OH基团的分布则反映了燃烧室内释热分布情况，故采该模型具有一定的可行性。

采用瑞利概率密度描述燃烧室释热分布的具体关系式为

$ q\left( x \right) = \frac{{\left( x \right)}}{{{\mu ^2}}} \times {{\rm{e}}^{ - \left( {\frac{{{x^2}}}{{2{\mu ^2}}}} \right)}} $

(8)

$ \frac{{{\rm{d}}Q\left( x \right)}}{{{\rm{d}}x}} = {\dot m_{{\rm{fuel}}}}q\left( x \right){E_{\rm{u}}}{\eta _{\rm{c}}} $

(9)

式(8)为瑞利概率密度函数，μ和η_c为经验常数，μ为反应释热分布形态的唯一参数，η_c为燃烧效率；式(9)为释热分布函数，Q(x)为燃料释热量，$ {\dot m_{{\rm{fuel}}}} $为燃料质量流率，E_u为燃料热值。燃料释热前，总温保持不变，始终等于燃烧室入口总温；燃料释热过程中满足式(9)所示关系。因此，通过积分，可得到超燃燃烧室总温分布

$ {T_{\rm{t}}}\left( x \right) = \left\{ \begin{array}{l} {T_{t, {\rm{in}}}} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left( {x \le {x_0}} \right)\\ {T_{{\rm{t}}, {\rm{in}}}} + \frac{{Q\left( {x - {x_0}} \right)}}{{{c_p}}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;(x > {x_0}) \end{array} \right. $

(10)

$ Q\left( x \right) = {m_{{\rm{fuel}}}}{E_{\rm{u}}}{\eta _{\rm{c}}}\left( {1 - {{\rm{e}}^{ - \left( {\frac{{{x^2}}}{{2{\mu ^2}}}} \right)}}} \right) $

(11)

式中x₀为释热起始位置；T_{t, in}为燃烧室入口总温，c_p为定压比热容。

2.3 计算方法

将式(10)代入式(6)、(7)中，即可对燃烧室芯流马赫数、芯流面积进行求解，经验常数μ，η_c取不同数值时，可得到不同的芯流面积分布。因此，选取适合的经验常数μ，η_c，即可得到适应实际工况的释热分布，从而求解芯流面积。

在超燃燃烧室芯流完成再附过程后，忽略附面层等因素的影响，可以认为芯流面积等于燃烧室的截面积，如图 1所示。因此，在再附点x_re后选取x₁，x₂两个截面作为参考截面，令x₁，x₂处计算得到的芯流面积A_c与参考截面的燃烧室截面积A₀相等，即可通过迭代计算获得适用当前工况的常数μ，η_c值；同时，流动再附后分离区消失，芯流面积等于燃烧室截面积，因此在x＞x₂段令A_c=A₀，芯流马赫数Ma为

$ \begin{array}{*{20}{l}} {\frac{{{\rm{d}}Ma}}{{Ma}} = \frac{1}{2}[\left( { - \frac{E}{{kM{a^2}}}} \right) + \left( {\frac{{{\rm{d}}\dot m}}{{\dot m}} - \frac{F}{{kM{a^2}}}} \right) - }\\ {\;\;\;\;\;\;\;\;\;\;\;\;\left( {1 + \frac{1}{{kM{a^2}}}} \right)\frac{{{\rm{d}}p}}{p} - \frac{{{\rm{d}}{A_0}}}{{{A_0}}}]} \end{array} $

(12)

综上所述，超燃燃烧室芯流计算流程如图 3所示。

3 算例验证及结果分析

本次研究分别对气动斜坡^[19]和凹腔^[20]两种超燃燃烧室构型进行了验证，实验结果来自于北京航空航天大学直连式超燃冲压发动机试验台，仿真结果由CFD++进行计算，湍流模型为k-Ω (SST)模型，燃烧模型采用有限速率模型。通过对比发现，一维计算与仿真结果具有较好的一致性。

3.1 验证算例一：基于气动斜坡的超燃燃烧室

基于气动斜坡的超燃燃烧室采用矩形截面设计，通过气动斜坡-燃气发生器结构促进火焰稳定，燃烧室入口马赫数为Ma=2，流量约为0.96kg/s，总压0.85MPa，总温1200K，燃料为乙烯，燃烧室结构如图 4所示，仿真选用1/2模型，并采用十组分十步化学反应机理进行计算。

图 5给出了当量比为0.3，0.37和0.5三种工况下实验和仿真的壁面静压分布，二者体现出较好的一致性，表明仿真结果能够较好地反映超燃燃烧室内部真实的流场结构，可以作为一维计算结果的验证算例。

图 6给出了三种工况下，x=0~0.6m区间内燃烧室芯流面积的计算结果，选取的参考截面为x₁=0.55m，x₂=0.6m。可以看出，流动过程的分离过程主要发生在x=0.2~0.45m区间内，且当量比越高，流动发生分离的位置越靠近上游，分离区内的芯流面积越小；在x≤0.2m段，芯流面积的大小存在波动，这是因为流动过程存在膨胀波，使壁面压力出现波动，进而影响到了芯流面积的计算结果。

图 7(a)，(b)，(c)给出了分离区内芯流面积和释热分布的一维计算结果和三维仿真结果，其中三维仿真的芯流面积是由燃烧室截面积减去回流区截面积获得的。整体而言，一维计算和三维仿真得到的释热分布具有较好的一致性，一维计算得到的芯流面积也能够较为准确地反应流动的分离再附过程；在数值上，一维计算得到的芯流面积略小于三维仿真结果，这是因为一维模型对燃烧室内的流动过程进行了简化，不能对附面层、流动不均匀等复杂的流动过程进行捕捉，而是将其影响反映到了芯流面积的计算结果中。

图 8给出了x方向速度u=0m/s等值面及x=0.3m，0.35m，0.4m，0.45m截面的密流(即密度与流向速度之积)分布云图。可以看出，回流区主要分布在x=0.3~0.4m段，在回流区下游，流密分布具有明显的不均匀性。

图 9(a)，(b)，(c)分别给出了Φ=0.3，0.37和0.5三种工况下，超燃燃烧室内马赫数分布的一维结果和三维结果。可以看出，两种工况下的一维计算结果均与三维计算结果具有较好的一致性。

为进一步分析参考截面x₁，x₂的位置对释热分布的影响，研究中选取的多组参考截面进行了计算，结果如表 1所示。考虑x=0.65m处存在压缩波，而压缩波、膨胀波会对芯流面积的计算结果产生一定影响，因此参考截面均选取在x=0.65m上游位置。从表中可以看出，参考截面的位置对经验常数的计算结果影响较小，表征燃烧效率的经验常数η_c均在0.9附近。在三维仿真中，由于燃气发生器对乙烯燃烧具有较强的促进作用，因此获得了较高的燃烧效率——95.7%。通过对比可以发现，一维计算结果与三维仿真结果的误差在7%之内。因此可以认为在避免膨胀波、压缩波的前提下，参考截面位置对一维计算结果的影响相对不大。

3.2 验证算例二：基于凹腔的超燃燃烧室

基于凹腔的超燃燃烧室构型采用矩形截面设计，燃烧室入口马赫数为Ma=2，流量约为1.5kg/s，总压2.0MPa，总温1900K，燃料为煤油，燃烧室结构如图 10所示，仿真选用1/2模型，并采用16组分23步化学反应机理进行计算。

图 11给出了实验和仿真的壁面静压分布，二者体现出较好的一致性，表明仿真结果可以用于一维计算的验证算例。

图 12和图 13给出了燃烧室内芯流面积和马赫数分布的一维计算与三维仿真的结果，选取x₁=0.95m，x₂=1.0m两截面作为参考截面。可以看出，虽然芯流面积的一维计算结果在x=0~0.3m区间内仍然受膨胀波、附面层的影响而出现波动，但整体而言，一维计算结果较为准确地捕捉到了燃烧室在x=0.3m~0.75m区间内流动分离再附的过程。与此同时，一维计算得到的马赫数分布也与仿真结果体现出较好的一致性，表明该一维方法可以用于以煤油为燃料、以凹腔为火焰稳定装置的超燃燃烧室芯流的计算。

4 结论

本文对双模态超燃燃烧室一维计算方法进行了研究，通过与两种燃烧室构型仿真结果的对比，可以得出以下结论：

(1) 通过对两个验证算例的分析，可以认为本文提出的一维分析模型能用于超燃燃烧室芯流面积的计算，对燃烧室芯流的分离再附过程和马赫数分布情况起到较好的预测作用。

(2) 通过选取参考截面获得超燃燃烧室释热分布情况的方法是可行的，并且通过该方法，可以使本文提出的一维分析方法适用于不同燃料、不同构型的超燃燃烧室芯流分析。

(3) 超燃燃烧室内的膨胀波、压缩波会对芯流面积的计算结果产生一定影响，因此参考截面应当在再附点后、膨胀波或压缩波出现前的位置选取。由验证算例一的结果表明，在避免膨胀波和压缩波的前提下，参考截面的位置对释热分布计算结果的影响不大，燃烧效率的波动在4%之内，与仿真结果的误差在7%之内。

后续的研究将针对如何降低膨胀波、压缩波对芯流面积计算结果的影响展开，从而进一步降低本文所述一维分析方法对参考截面选取位置的约束，使该一维分析方法具有更好的适用性。