基于动网格技术的固体燃料冲压发动机燃面瞬态退移速率研究①

魏 韬，武晓松

(南京理工大学机械工程学院航空宇航系，南京 210094)

基于动网格技术的固体燃料冲压发动机燃面瞬态退移速率研究①

魏 韬，武晓松

(南京理工大学机械工程学院航空宇航系，南京 210094)

为了研究固体燃料冲压发动机(SFRJ)燃面退移速率在工作过程中的变化特性，基于发动机工作特点及动网格技术，考虑到燃烧流动及燃料表面的对流、辐射换热与燃料热解退移等过程耦合的影响，建立了SFRJ燃面瞬态退移速率预示方法，并对某带补燃室、以聚乙烯(PE)为燃料的试验发动机的燃烧室-喷管统一内流场进行数值计算，得到在移动边界条件下的瞬态流场分布，并分析了内弹道参数云图及其随时间的变化规律。结果表明，燃烧主要发生在当量比函数φ在-2～2之间的区域;随着发动机工作，燃速逐渐降低，且再附点向下游移动，燃料通道出口处流速和温度有降低趋势;此外，在小型发动机工作初期，燃料通道尾部出现类似固体火箭发动机的侵蚀燃烧现象。研究表明，该方法能成功求解发动机复杂的非定常工作过程，较好揭示燃面退移过程。所得结论对发动机设计和试验具有一定指导意义。

固体燃料冲压发动机;瞬态退移速率;动网格;数值仿真

0 引言

在固体燃料冲压发动机工作过程中，燃面瞬态退移速率随时间和轴向位置变化，形成的燃气内流场几何边界十分复杂，流动区域不断变化，且受结构和工作条件影响，难以用统一公式描述。用试验方法获得燃面瞬态退移速率所需系统复杂，成本昂贵，而传统数值计算方法采用准定常假设，显然不能准确描述发动机复杂非定常流场的细节。目前，针对SFRJ燃面瞬态退移速率应用动网格技术和瞬态方法，较详细地反映流场变化规律的研究，国内外很少有报道。

SFRJ燃面退移速率的传统分析方法［1-4］是基于边界层假设的传热理论，来确定向燃面的换热量及燃速。研究者通常采用准定常假设［3-6］，忽略热辐射［1，3-5］，甚 至 假 设 燃 面 温 度 和 有 效 汽 化 热 为 常数［2-4］，这样的简化不能考虑到各种复杂物理化学过程的相互作用。因此，针对SFRJ在运动边界条件下的非定常流场，建立更为精确的燃面瞬态退移速率预示方法具有重要意义和必要性。

针对以上问题，基于动网格技术，将燃烧流动及固体燃料表面的对流、辐射换热与燃料热解退移等过程耦合计算。采用低Re数k-ε湍流模型，并考虑加质对对流换热的影响和温度对燃料有效汽化热的影响，建立SFRJ燃面瞬态退移速率的预示方法，并对代夫特科技大学试验发动机［3］的燃烧室、喷管流场进行一体化计算。研究移动边界条件下的瞬态流场，并分析内弹道参数云图、燃面瞬态退移速率、再附点、燃料通道出口参数及其随时间的变化规律。研究结果可为SFRJ的工程设计和试验提供参考。

1 物理数学模型

1.1 物理模型

图1为试验发动机结构简图［3］。燃料药柱通道直径为45 mm，长度为300 mm，入口直径为15 mm，突扩台阶高度为15 mm，补燃室长度为100 mm，喉部直径为20 mm。固体燃料为聚乙烯(PE)。

图1 试验固体燃料冲压发动机简图Fig.1 Configuration of experimental SFRJ

1.2 基本假设

固体燃料冲压发动机工作过程是一个极其复杂的能量转化过程。其中，燃料热解退移、燃烧、湍流流动和传热等均是复杂的物理化学过程，且相互之间高度耦合。目前，对其进行严格的数学描述尚存在诸多困难。因此，本文假设:(1)发动机为二维轴对称模型;(2)聚乙烯热解产物只含C2H4，燃气为纯气相;(3)化学反应为C2H4-空气两步反应模型;(4)不考虑壁面烧蚀与传热。

1.3 控制方程

1.3.1 流动控制方程

对于具有运动边界的有限体积控制体，采用任意拉格朗日-欧拉(Arbitrary Lagrange-Euler-ALE)有限体积法描述的积分形式守恒型控制方程为

式中ρ为流体密度;φ为通用变量为燃气速度为燃面网格运动速度;Γφ为广义扩散系数;Sφ为广义源项。对流项和扩散项采用二阶迎风格式离散。

1.3.2 气-固相交界面控制方程

固体燃料和燃气的耦合，通过交界面质量和能量守恒方程实现，即

式中为燃料热解产物的质量流率;ρs为固体燃料密度，取926 kg/m3;˙rb为当地燃面退移速率;和分别为对流和辐射换热项=(Tw-T0)为燃料内部温升吸热率，其中cs为固体燃料比热容，取2 142 J/(kg·K)，Tw为当地燃面温度，T0为燃料初温，取300 K;=为燃料热解、气化潜热，其中hv为燃料有效汽化热。

1.4 数学子模型

除以上控制方程外，还须附加数学子模型来描述化学反应、湍流、燃料热解、对流和辐射换热以及燃料有效汽化热。

1.4.1 气相化学反应模型

固体燃料冲压发动机再附点下游的燃料表面上形成湍流边界层，在边界层内形成湍流扩散火焰，发动机中的燃烧主要为扩散燃烧，反应速率受气体扩散过程控制。因此，燃烧模拟采用涡耗散模型(EDM)。用简化的两步总包反应模拟气相化学反应，即

1.4.2 湍流模型

由于燃面存在不均匀的质量和能量输运，在靠近燃面的区域，湍流脉动动能强烈衰减，而耗散率达到最大值［7］。为使计算能从高Re数区域一直进行到燃面(该处Ret=0)，本文使用考虑了近壁面效应的低Re数k-ε湍流模型(low-Re k-εmodel)［8-9］封闭方程组。湍流输运引入的源项为

其中，ηt=cμ∣fμ∣，f1=1.0，f2=1-0.22 ×exp［-(Ret/6)2］，fμ=1-exp(-0.011 5y+)，模型常数cμ=0.09，c1=1.35，c2=1.8，σk=1.0，σε=1.3。

1.4.3 固体燃料热解模型

根据Hadar Ian［10］等人的研究，碳氢型固体燃料的退移速率与燃面温度的关系符合Arrhenius公式，即

式中A、Ea和R分别为指前因子、活化能和通用气体常数。

模型常数通过热分解实验［11］获得，对于PE，A=8.25 ×105mm/s，Ea=133 539.35 J/mol。将式(6)代入式(3)，采用牛顿迭代法解得Tw，进而得到当地燃速。

1.4.4 对流换热模型

燃气向固体燃料对流换热热流密度为

式中T∞为边界层内火焰面温度;h为对流换热系数。

对圆管内壁有质量(即燃料热解气体)加入的湍流流动，h的计算式为［12］

式中cp、ρ和Pr分别为燃面附近燃气定压比热容、密度和普朗特数;uin为入口流速;Red为按uin计算的通道雷诺数为平均燃速。

1.4.5 辐射换热模型

燃气向燃料辐射换热净热流密度为［13］

式中σ为Stefan-Boltzmann常数;ε＇w=(εw+1)/2为燃料的有效发射率，将燃料壁面考虑为漫射的灰体表面［14］，并考虑到炭黑附着，发射率取εw=0.8;εg和αg分别为燃气的发射率和吸收率;Tc为燃气静温。

热力计算表明，火焰区主要气体为CO2、H2O、N2和O2，其中N2和O2为非极性对称型双原子气体，对热射线的发射和吸收能力微弱，可认为是透明体，故不考虑其对辐射换热的贡献，而CO2和H2O等多原子气体是主要辐射源。因此，燃气发射率εg和吸收率αg分别为

查文献［13］中图 7-1到 7-3并计算得，εg=0.024 85，αg=0.175 55。

1.4.6 燃料有效汽化热

燃料的有效汽化热与温度有关，并取决于传热机理，准确的hv对于精确预测燃速非常关键［3］。根据文献［15］，hv可表达为

式中cp为热解产物的定压比热容;hmel为解链反应热;hvap为热解产物的汽化热;hpyr为产生热解产物所需的反应热。

对于 PE，hmel=225 kJ/kg，hvap=485 kJ/kg，hpyr=3 335 kJ/kg。

Rihani-Doraiswamy(RD)基团贡献法在求解SFRJ燃烧室环境中化合物的定压比热容时有更高的精度［6］，它将cp表示为

式中ni为第i种基团的数量;ai、bi、ci和di分别表示不同基团对化合物定压比热容的贡献。

应用RD法所得C2H4定压比热容的基团贡献值见表1。C2H4的定压比热容为

表1 Rihani-Doraiswamy基团贡献法C2H4定压比热容的基团贡献值Table 1 Heat capacity of C2H4organic compounds from group contribution by RD method

1.5 初始计算网格

利用Gambit生成计算网格，局部网格见图2。由于燃面节点退移速率不同，变形后燃面形状复杂。为了动网格的实现，在燃面附近采用加密的非结构网格，以适应不规则的气流通道。为提高计算精度和速度，对计算域进行分区，在远离燃面的区域，采用稍稀疏的四边形结构网格，并对喉部和近壁区进行加密。初始网格数为48 746。

图2 局部计算网格(燃烧室头部)Fig.2 Partial computation mesh

1.6 边界及初始条件

初始条件采用初场赋值方法:首先在定网格条件下，利用DEFINE_PROFILE宏定义燃速，计算定常流场。待求解收敛后，采用动网格技术模拟燃面退移的非定常流场。

边界条件:(1)燃面边界。包括加质和移动边界条件。燃面节点退移速率由式(6)给出。(2)固体壁面边界。绝热、无滑移条件。(3)轴对称边界。沿轴线法向速度为零，各物理量在轴线上的法向梯度为零。(4)空气入口边界。采用质量流率入口边界条件=0.15 kg/s，总温T0=600 K。(5)喷管出口边界。喷管出口为超音速流，流动参数由上游参数二阶外推插值得到。

1.7 动网格技术

为了动态模拟SFRJ工作过程中燃面瞬态退移过程，需使用动网格技术，根据计算得到的当地燃面退移速率实时更新网格。

1.7.1 动网格更新算法

动态网格更新算法通常有动态层铺法、弹性光顺法及局部网格重构法。针对燃面变形复杂和运动幅度较小的特点，本文联合使用弹性光顺法和局部网格重构法。

(1)弹性光顺法。该网格更新算法假定相邻两节点之间的边为理想弹簧。变形前的初始位置为平衡状态，当节点运动时，与其相连的各边将产生正比于边长变化的力。将作用于节点的力写为胡克定律的形式，即

平衡状态下，作用于节点的合力为零，可得迭代方程为

当得到节点位移后，通过对内部节点的Jacobi迭代来求解方程。当迭代收敛后，则按式(17)更新节点位置

式中n+1和n分别表示下一时刻和当前时刻的节点位置。

(2)局部网格重构法。当边界位移较小时，仅靠网格变形就能适应边界运动，则采用弹性光顺法;当边界位移相对当地网格尺寸较大时，网格变形会形成严重扭曲的单元，或发生网格退化，下一时刻的求解可能会产生收敛性问题。局部网格重构法将不符合畸变率或尺寸标准的网格聚集起来，通过插值运算，在局部重新生成网格，以保证计算精度。若新网格满足畸变率和尺寸标准，就对网格进行局部更新，否则放弃新网格。

1.7.2 计算流程

图3为流场计算程序流程图。主程序主要有更新流场参数、流场计算和动网格技术更新网格等模块。其中，更新流场参数模块在每个时间步开始时，利用上一时刻流场参数的计算结果，对当前流场进行更新;流场计算模块为Fluent提供的流场计算程序;动网格模块的作用是在当前时间步下流场计算收敛后，计算并更新燃面在下一时刻的位置，同时添加源项。其中，燃面的运动利用DEFINE_GRID_MOTION宏定义。

图3 流场计算程序流程图Fig.3 Flow chart of CFD program

1.8 模型校验

为验证所建立计算方法的准确性，对文献［3］的试验发动机进行计算，燃面退移速率计算值与试验结果的对比如图4所示。可见，两者基本吻合，表明所建立的预示方法适用于SFRJ燃面退移速率的计算。

图4 燃面退移速率计算值与试验数据对比Fig.4 Comparison between experimental data and computational results

2 计算结果与分析

选用Fluent提供的隐式、非定常耦合求解器，采用UDF编程方式处理燃面退移。选取时间步长为0.05 s，迭代步数为600，内迭代最大步数为100。由于发动机长径比较大，为有效描述内流场特性，将发动机的X轴(轴向)与Y轴(径向)坐标之比设为1/3。

图5为14.0 s时刻发动机中速度及流线分布图。可看出，突扩台阶后及节流板后形成两个回流区，它们将在一定程度上强化空气和燃料的掺混，并延长气体在燃烧室的停留时间，有利于扩散燃烧的充分进行;此外，突扩台阶后的回流区还会加强燃气对燃料的对流换热，这有利于提高局部燃面退移速率。发动机轴线附近的流速在突扩台阶和节流板后，受回流涡旋的影响有一个减小过程;流速沿轴向有较大变化，这是由于燃料热解加质、燃烧及通道结构引起的;同时，燃面附近流速很小，且存在回流区，这正是SFRJ能维持稳定燃烧的原因。

图5 发动机中速度和流线图Fig.5 Velocity and path-line in SFRJ

图6为发动机工作过程中速度沿轴线的变化曲线(喷管加速过程未完全画出)。可见，轴线上的速度变化趋势相同;燃气进入补燃室之前，流速有先减小后增大的过程，减速是由于节流板的阻滞，使流动以较大的径向速度向轴线汇集，加速是由于节流板通道直径较小，流动连续性要求使然;此外，在突扩台阶和节流板后流速有不同程度的降低，这也印证了对图5的速度分析。应当指出，随着发动机工作，燃料通道截面逐渐增大，通道出口处(x=0.29 m)速度逐渐降低，在6 s时刻，速度为130 m/s，到26 s时刻，速度仅为102 m/s;燃烧室入口处的速度基本维持在180 m/s，这是由于空气质量流率保持不变。

为直观地表述燃烧室中燃料和氧化剂的分布状况，引入组分当量比函数φ=lg(Yoxid/Yfuel/φ0)。其中，Yoxid和Yfuel分别为当地氧化剂和燃料的摩尔分数;φ0为氧化剂和燃料的摩尔恰当比。在φ=0的区域，氧化剂和燃料完全反应;在φ＜0的区域，燃料剩余;在φ＞0的区域，氧化剂剩余。

图6 不同时刻沿发动机轴线的速度曲线Fig.6 Curves of velocity along axis at different time

图7为14.0 s时刻发动机中组分当量比函数φ和温度等值线图。可见，燃烧主要集中在燃料表面和轴线之间的一定区域内(-2＜φ＜2)，在此区域外，氧燃比偏离反应恰当比φ0较大，几乎无燃烧发生;轴线前半部及燃面附近氧气与燃料混合较差，故温度较低;火焰面最高温度达2 400 K，沿轴线向后，由于氧气的消耗，燃烧区逐渐扩展至轴线上;此外，补燃室温度达2 200 K以上，说明补燃室可提高燃烧效率。

图7 发动机中当量比函数和温度等值线图Fig.7 Contours of species equivalent ratio function and temperature in SFRJ

图8为不同时刻发动机轴线上温度变化曲线。可看出，温度变化规律相同;结合图7可知，在燃烧室前部，由于氧气充足，燃料未扩散至轴线处即被消耗，因此温度基本等于来流温度;在x=0.2 m处，温度开始迅速上升，这是由于沿轴线向后，随着氧气消耗，燃料与氧气在越来越靠近轴线的区域混合充分，并发生扩散燃烧;此外，补燃室轴线温度相对稳定，气体在喷管中膨胀加速，因而温度迅速下降。应当指出，随着发动机工作，燃料通道出口处温度有降低趋势，在6 s时刻温度为2 435 K，到30 s时刻温度为2 285 K。这是因为燃料通道截面逐渐增大，空气相对集中在轴线附近，燃料和氧气混合更困难，使火焰面远离轴线。

图8 不同时刻沿发动机轴线的温度曲线Fig.8 Curves of temperature along axis at different time

图9为发动机工作过程中燃面瞬态退移速率变化曲线。可发现，燃面退移速率沿轴向有较大变化，且不同时刻燃面退移速率的变化趋势基本相同，都是先增大后减小，这是由于燃面退移速率与热解表面温度之间符合Arrhenius热解速率公式。这造成燃料通道呈现出两头细中间粗的型面，如图5和图7所示;紧靠突扩台阶后面存在一个流动滞止点，燃料难以与氧气充分反应，使局部燃料绝对富余(当量比函数φ为-7左右)，因此燃速很低，约为0.05 mm/s;在此滞止点后，回流区的产生强化了燃料和氧气的掺混，以及燃气向燃料表面的对流换热，使得燃速迅速上升，并在再附点处达到峰值，在6 s时刻再附点(x=104 mm)燃速为0.399 mm/s;峰值后边界层充分发展，随着氧气的消耗，火焰面向轴线靠近，对固体燃料热解表面的传热减弱，因而燃面退移速率逐渐降低。

随着发动机工作，燃面退移速率逐渐降低，这是由于燃料通道截面逐渐增大，使空气相对集中在轴线附近，燃料和氧气在远离燃面的区域发生扩散燃烧，从而传向燃面的热量减少。还可发现，随着发动机工作，再附点逐渐向下游移动。这是由于燃料通道直径增大，使得台阶高度增加引起的。此外，在发动机工作初期(t＜10 s)，燃料通道出现类似于固体火箭发动机中的侵蚀燃烧现象，即燃料通道尾部的局部退移速率呈缓慢增大趋势。这是因为燃烧初期小型发动机的燃气通道较小，平行于燃料表面的流速大于出现侵蚀效应的临界速度。

图9 不同时刻燃面瞬态退移速率曲线Fig.9 Curves of instantaneous regression rate of the solid fuel grain at different times

3 结论

(1)通过对比计算结果与试验数据，验证了所建立的固体燃料冲压发动机燃面瞬态退移速率预示方法的准确性。该方法成功求解了发动机复杂的非定常工作过程，较好地揭示了燃面不规则退移过程。

(2)固体燃料冲压发动机的流场具有明显不均匀性，燃烧主要发生在燃料表面和轴线之间的一定区域，该区域的当量比函数φ在-2～2之间。

(3)在燃烧室入口处，燃速很低，之后形成的回流区强化了燃料和氧气的掺混及对流换热，使局部燃面退移速率迅速升高，并在再附点处达到峰值，峰值后平缓降低。

(4)随着发动机工作，燃面退移速率逐渐降低，且再附点逐渐向下游移动;燃料通道出口处的流速和温度有降低趋势。

(5)在小型发动机工作初期，燃料通道尾部出现类似固体火箭发动机的侵蚀燃烧现象。

［1］Krishnan S，Philmon G.Solid fuel ramjet combustor design［J］.Progress in Aerospace Sciences，1998，34:219-256.

［2］Korting P A O G，Van der Geld C W M，Wijchers T，et al.Combustion of polymethylmethacrylate in a solid fuel ramjet［J］.Journal of Propulsion and Power，1990，6(3):263-270.

［3］Elands P J M，Korting P A O G，Wijchers T，et al.Comparison of combustion experiments and theory in polyethylene solid fuel ramjets［J］.Journal of Propulsion and Power，1990，6(6):732-739.

［4］向敏.固体燃料冲压增程炮弹工作过程仿真及性能分析研究［D］.长沙:国防科技大学，2006.

［5］刘巍.固体燃料冲压发动机燃烧组织技术研究［D］.长沙:国防科技大学，2010.

［6］夏强.固体燃料冲压发动机工作过程研究［D］.南京:南京理工大学，2011.

［7］So R M C，Zhang H S，Speziale C G.Near-wall modeling of the dissipation rate equation［J］.AIAA Journal，1991，29(12):2069-2076.

［8］陶文铨.数值传热学(第2版)［M］.西安:西安交通大学出版社，2001.

［9］Chien K.Predictions of channel and boundary-layer flows with a low-Reynolds-number turbulence model［J］.AIAA Journal，1982，20(1):33-38.

［10］Hadar I，Gany A.Fuel regression mechanism in a solid fuel ramjet［J］.Propellants，Explosives，Pyrotechnics，1992，17(2):70-76.

［11］De Wilde J P.Fuel pyrolysis effects on hybrid rocket and solid fuel ramjet combustor performance［D］.Delft:Delft University of Technology，1991.

［12］王守范.固体火箭发动机燃烧与流动［M］.北京:北京工业学院出版社，1987.

［13］郑亚，陈军，鞠玉涛，等.固体火箭发动机传热学［M］.北京:北京航空航天大学出版社，2006.

［14］杨世铭，陶文铨.传热学(第四版)［M］.北京:高等教育出版社，2006.

［15］De Wilde J P.The heat of gasification of polyethylene and polymethylmethacrylate［M］.Delft:Delft University of Technology，1988.

Study of instantaneous regression rate in solid fuel ramjet based on dynamic mesh

WEI Tao，WU Xiao-song
(Department of Aeronautics and Astronautics，School of Mechanical Engineering，Nanjing University of Science and Technology，Nanjing 210094，China)

In order to study the changes in regression rate of solid grain in solid fuel ramjet(SFRJ)，a method to indicate the instantaneous regression rate of solid grain in SFRJ was developed based on the operational feature，dynamic mesh，as well as coupled simulation of gas-phase combustion，heat transfer and regression of solid grain.Then whole inner flow of an experimental motor with additional chamber and polyethylene(PE)fuel was numerically simulated.The time-dependent flow was obtained with dynamic boundaries，and the internal ballistic parameters were analyzed.The results show that combustion mainly occurs in the area where φ is between-2 and 2;as motor works，the regression rate decreases，the reattachment point moves downstream，and the velocity and temperature reduce at the outlet of solid grain;in addition，there is erosive burning in small SFRJ during initial operation stage.The study shows that this method can simulate the unsteady working process and fuel regression.The conclusions offer instruction for designing and experiment of solid fuel ramjet.

solid fuel ramjet;instantaneous regression rate;dynamic mesh;numerical simulation

V435

A

1006-2793(2012)04-0450-07

2011-06-10;

2012-02-17。

国防预研项目。

魏韬(1986—)，男，硕士，主要从事固体燃料冲压发动机工作过程研究。E-mail:nustweitao@gmail.com

(编辑:崔贤彬)