三角形网格CE／SE方法在带锥形尾喷管两相脉冲爆轰发动机流场计算中的应用

王研艳，翁春生

（南京理工大学 瞬态物理国家重点实验室，江苏 南京 210094）

0 引 言

脉冲爆轰发动机（Pulse Detonation Engine，简称PDE）是一种利用脉冲爆轰波来产生周期性冲量的新概念非定常推进系统，具有热效率高、结构简单、单位燃料消耗率低、工作范围广等优点。但如何充分利用爆轰产物的高温高压能量来提高脉冲爆轰发动机推力仍为得到解决，方法之一就是在发动机尾部安装尾喷管。常规发动机尾喷管增加推进效能的原理是设计一个最佳喷管构型，将燃烧室内排出的高温高压燃气的热能尽可能多地转化为动能，以增加发动机的推力。而脉冲爆轰发动机尾喷管中的流动是周期性的、非定常的，存在复杂激波和膨胀波，甚至存在爆轰波，在同一周期的不同时刻喷管内流场不尽相同，因此很难找到一个能适合所有时刻的喷管构型，这就使得建立在传统定常流动基础上的喷管设计理论和概念在这里不再适用［1］。Nobuyuki T等人［2］利用 Strang算子分裂法和迎风TVD格式数值模拟了多循环下带收敛扩张喷管的PDE的内外流场，研究了以H2为燃料时喷管进出口面积比和燃料填充率对性能的影响，得出在不同飞行马赫数条件下存在不同的最佳喷管面积比的结论。秦亚欣等人［3］对带不同结构喷管的PDE进行了二维轴对称单循环数值研究，指出常规喷管中扩张喷管可以产生最大的瞬时推力峰值和比冲。王志武等人［4］试验研究了锥形喷管的进出口面积比和长径比对吸气式脉冲爆轰发动机推力的影响，指出不同频率下各类型喷管对发动机的平均推力影响不同，但总体来讲收敛型尾喷管的增推效果最好。李建中等人［5］采用基元反应模型和显式迎风TVD差分格式模拟了带共用尾喷管多管脉冲爆轰发动机各爆轰管内以及共用尾喷管内的流场。Caldwell等人［6］实验研究了不同填充率下引射器扩张角度对脉冲爆轰发动机流场和推进性能的影响。总体而言脉冲爆轰发动机尾喷管的最佳设计与PDE的频率、飞行马赫数、喷管类型、高度、燃料填充率等工作参数紧密相关［7－8］，所以至今还未得出非定常脉冲爆轰发动机尾喷管的设计准则。

CE／SE方法最先由Chang提出，它将时间与空间统一处理，从守恒律积分方程出发，设立守恒元和（求解元，使求解格式局部和整体都保证物理意义上的守恒。由于CE／SE方法无需黎曼分解，具有计算格式简单、精度高、捕捉间断能力强等优点，所以近年来CE／SE方法成为国内外学者研究的焦点之一。Chang等人［9］在最初提出的CE／SE方法的基础上进行了改善，提高了它的计算精度和稳定性。Gary C Cheng等人［10］在Chang的基础上将CE／SE方法进一步发展，针对粘性流采用长扁型网格发展了CE／SE的几种计算格式，并进行了算例验证和效率分析。Wang等人［11］运用矩形网格CE／SE方法对两相脉冲爆轰发动机内流场进行了数值研究。李昕等人［12］运用矩形网格CE／SE方法求解流体动力学N－S方程与麦克斯韦方程，实现了对电磁轨道炮中等离子体电枢的数值研究。刘建文等人［13］考虑到流场的不规则性运用非结构网格CE／SE方法对单相多管脉冲爆轰发动机模型进行了数值研究。

本文在Chang的基础上推导了带源项的结构化三角形网格CE／SE方法计算格式，这方面的内容在国内还未见报道。而此结构化三角形网格CE／SE方法与矩形网格CE／SE方法相比无需坐标变换，网格分布均匀；与非结构网格CE／SE方法相比推导过程简化，计算格式简单。本文应用此CE／SE方法数值研究带锥形扩张尾喷管的两相脉冲爆轰发动机的内外流场。分析尾喷管内不同填充状态时PDE的内外流场及推进性能。

1 理论模型

计算采用轴对称两相燃烧转爆轰控制方程［14－15］：

式中，α＝0为二维流动；α＝1为轴对称流动。其中下标g、l分别表示气相和液相；φg、φl为体积分数比，满足φg＋φl＝1；ρ、u、v、p分别是密度、轴向速度、径向速度、压力。E是单位总能，E＝e＋（u2＋v2）。Id、Fd、Qd、Qc的含义和表达式见参考文献［15］。

2 三角形网格CE／SE方法

2.1 三角形CE／SE方法中守恒元与求解元的确定

CE／SE方法［9－10，16］将整个空间－时间计算区域划分为若干个求解元。在每个求解元内，假设流场的变量是连续的，并可以用Taylor级数展开，穿过相邻求解元的边界，流场的变量可以不连续。而在每个网格点对应的守恒元上，空间－时间的积分通量是守恒的。考虑如图1（a）所表示的交替网格，每个三角形的中心点用实心圆和空心圆分别表示其交替网格点。令空间Ω表示E3上的所有网格点M（r，s，n）的集合。令空心圆点表示tn时刻，实心点表示tn＋1／2时刻。如图1（b）所示，GMBFDHC表示tn时刻，G′M′B′F′D′H′C′表示tn＋1／2时刻，G″M″B″F″D″H″C″表示tn＋1时刻。那么点G′对应的求解元S（G′）为GBB′G′，GDD′G′、GCC′G′和M′B′F′D′H′C′；对应的守恒元为四棱柱GBFDG′B′F′D′，GDHCG′D′H′C′和GCMBG′C′M′B′。

2.2 计算格式

在方程（1）中，根据Chang对时间和空间统一处理 的思想［9－10，14，16］，令x1＝x，x2＝y，x3＝t作为三维欧氏空间的3个坐标，且令S＝R－；h＝（F，G，U），利用高斯定理对其积分有：

图1 二维三角形网格CE／SE方法中时空几何图形Fig.1 Spatial scheme of two－dimensional CE／SE method using congruent triangular meshes

式中，S（V）是V在E3上任意时空区域的边界，h＝（F，G，U）是时间－空间通量流密度矢量。

引入雅克比系数矩阵A、B：A＝因为F和G均是U的齐次函数，所以F＝AU，G＝BU。这样h进一步表示为：h＝（AU，BU，U）。

根据求解元和守恒元的定义，对于任意的（x，y，t）∈SE（r，s，n），U（x，y，t）、F（x，y，t）、G（x，y，t）可用其相应的离散量U＊（x，y，t；r，s，n）、F＊（x，y，t；r，s，n）、G＊（x，y，t；r，s，n）近似代替。利用泰勒级数展开有：

所以方程（2）用离散通量表示如下：

为了方便计算引入一个新的坐标系，如图2。新坐标系下取单位长度：Δη＝MH＝和Δξ，其中w＝

图2 （x，y）和（ξ，η）坐标系Fig.2 Cartesian coordinates and new coordinates

在各CE上对式（3）进行求解，得到三角形网格CE／SE方法求解格式如下：

当r＝4j＋1，s＝2k＋1，或当r＝4j＋3，s＝2k（s＝1，2，3…；k＝1，2，3…）时：

当r＝4j，s＝2k＋1，或当r＝4j＋2，s＝2k（s＝1，2，3…；k＝1，2，3…）时：

上式中，I是单位矩阵。

和Uη＋求解如下［14］：

式中：α是一个可调的参数，通常取1或2。

式（6）中：

2.3 源项的处理

方程（1）中等号右侧出现了化学反应源项，由于化学反应的特征时间远小于对流的特征时间，其源项是刚性的，因此需要采用龙格－库塔法对方程（1）中的化学反应源项进行求解［14］。具体方法是（1）先考虑除化学反应之外的其他源项，在CE／SE格式下求出（2）将作为初值，用龙格－库塔法来求解方程＝R。龙格－库塔法的时间步长为：Δt＝rk，N通常取5～20。

3 计算结果及分析

采用三角形网格CE／SE方法对带扩张尾喷管PDE的内外流场进行耦合计算，其轴对称计算区域如图3，AEFGHA为PDE爆轰管和尾喷管平面计算区域，尾喷管扩张角度为30°。EDCBHGFE为PDE外流场计算区域。

初值条件：模拟的初始时刻PDE有两种填充方案：方案1：爆轰管和锥形管作为整体，均填充化学当量比的汽油空气混合物；方案2：将锥形管仅作喷管使用，即爆轰管内填充化学当量比的汽油空气混合物，但锥形管内只填充空气。同时外流场区域则布满均匀的空气。管内外气体速度均为0m／s，温度为298K，压力为0.1MPa，管内液相汽油液滴半径为50μm。

为了验证运用三角形网格CE／SE方法模拟PDE流场的准确性，先对不带尾喷管的PDE内外流场进行数值计算，并与实验结果进行对比验证。比较发现模拟得到的图线变化趋势和实验结果吻合得很好（图4）。数值计算的爆轰波波阵面很陡，只有2～4个网格，说明此三角形网格CE／SE方法能够精确捕获爆轰波等强间断。

图3 PDE模型计算区域示意图Fig.3 Scheme of the PDE computational domain

图4 爆轰管尾部压力随时间变化曲线Fig.4 Pressure contribution vs.time at the end of PDE

图5是采用填充方案1时PDE的内外流场压力云图，图6为对应的PDE中心轴线上压力随x变化曲线。联系图6分析图5，由图5（a）知PDE内已形成稳定的爆轰波，爆轰波后是从封闭端发出的Taylor膨胀波扇，以满足封闭端速度为零的条件。在封闭端和Taylor膨胀波波尾之间是均匀的平台压力区。从图6中曲线5、6可以看出形成稳定爆轰时爆轰压力为1.97MPa，均匀区平台压力为0.48MPa。图5（b）中爆轰波传播到锥形管内，燃料燃烧所释放的热量维持爆轰波在管内继续传播。锥形管的扩张作用使得爆轰压力略有下降，爆轰波波面出现向右弯曲现象。由图6曲线7知，至t＝0.58ms时，中心轴线上爆轰压力值降为1.85MPa。图5（c）中爆轰波传出锥形管，失去能量支持后退化为球形激波在外流场扩散。锥形管出口处产生一系列膨胀波向上游传播，管内受到外流场的干扰压力降低。图5（d），外流场球形波向左传播至0.61m，向右传播到1.44m，锥形管出口外侧出现明显的涡环区域，锥形管外侧背压发生变化，发动机推进性能受到影响。发动机管内膨胀波传播至0.36m，管内压力降低，平台压力区缩短，由图6曲线16可知随着排气过程的进行，t＝2.91ms时发动机头部压力降低到0.14MPa。

图5 采用填充方案1时PDE内外流场压力分布云图Fig.5 Pressure distribution in and out of PDE with fuel filled in the conical tube

图6 采用填充方案1时PDE内外流场轴线上压力分布曲线Fig.6 Pressure distribution on the symmetry axis of PDE with fuel filled in conical tube

为了进一步分析锥形管内外的流动情况，图7给出了采用填充方案1时PDE内外流场的速度矢量图。图7（a）中速度沿锥形管壁面分布解释了扩张管内爆轰压力降低的原因。图7（b）气流流出扩张管后在出口外侧出现速度涡环，这是因为当曲面爆轰波上下边界传出扩张管后，管内欠膨胀的高温高压气体尾随激波进入大气环境中，其斜压作用使得管口外侧出现涡环。

图7 采用填充方案1时PDE内外流场不同时刻速度矢量图Fig.7 Velocity vector distribution in and out of PDE with fuel filled in conical tube

图8是采用填充方案2时PDE内外流场的压力分布云图，图9为对应的PDE中心轴线上的压力随x变化曲线。联系两者分析：图8（a）中爆轰波传播到喷管后退化为无化学反应的主激波，受喷管扩展壁面产生的膨胀波影响，在靠近壁面处发生弯曲，形成曲面激波。此时锥形管内流场比方案一中更加复杂。从图9的曲线7、8、9可以看出，随着主激波在扩张喷管中的传播，激波越来越弱，波后气体的压力越来越小。至t＝0.59ms时，波后压力值为1.58MPa。图8（b）中主激波传出喷管，以球形激波向外流场扩散，喷管内已存在的膨胀波系和外流场新产生的膨胀波共同作用，降低管内压力。对比图6和图9发现，方案二的排气过程比方案一时快。具体体现在，t＝1.27ms时，方案一的膨胀波降低PDE内部平台压力至x＝0.34m处；而方案二中PDE内部x＝0.30m处的平台压力已经开始下降。

图8 采用填充方案2时PDE内外流场不同时刻压力分布云图Fig.8 Pressure distribution in and out of PDE at different times，without fuel in conical tube

图9 采用填充方案2时PDE内外流场轴线上压力分布曲线Fig.9 Pressure distribution on the symmetry axis of PDE at different times，without fuel in conical tube

图10是爆轰波传入锥形管后中心轴线上的马赫数分布曲线，比较（a、b）图发现当锥形管内不填充燃料时，爆轰波传入锥形管内没有化学反应的支持退化为激波，但是波峰位置处的马赫数反而变大，分析其原因：在锥形管内不填充燃料时波峰位置处压力变小，速度值也变小，但前者的变化幅度远大于后者，导致当地声速变小，且其减小量大于速度减小量。

图10 爆轰波传入锥形管后中心轴线上马赫数分布曲线Fig.10 Mach number distribution on the symmetry axis of PDE at different moments

考虑到脉冲爆轰发动机非定常特性，论文中采用表面力积分来分析带锥形尾喷管PDE的推进性能［17］。图11为不同填充方案时PDE瞬时推力F、冲量I、燃料比冲Ispf随时间变化曲线。图中曲线1对应填充方案1；曲线2对应填充方案2。分析三幅图发现，当锥形扩张管内不填充燃料时，它起到喷管的作用，使爆轰管内排出的高温高压燃气在锥形管内得到膨胀，转化为动能，从而增加了发动机的推力，比较图11（c）图中两条曲线发现，至t＝4.8ms时此方案发动机的燃料比冲比另一方案的增加了15.5%，说明这种方案提高了燃料的利用率；但当锥形扩张管内填充燃料时，它作为发动机爆轰的一部分，充分利用了脉冲爆轰发动机的空间，增加了发动机的瞬时推力和冲量，比较图11（b）图中两条曲线发现，至t＝4.8ms时此方案发动机的冲量比另一方案增加了11.1%。显然，当我们考虑在每个周期中增加PDE的推进性能时，需要利用所有可以利用的空间，此时可以考虑将其锥形尾喷管作为发动机的一部分，在其中填充或者部分填充燃料，从而增加发动机的总体推力。

图11 锥形管不同方案时PDE的瞬时推力、冲量和燃料比冲比较Fig.11 The comparison of thrust，impulse and fuel－specific impulse of PDEs with different filling conditions in conical tube

4 结 论

本文推导了带源项三角形网格CE／SE方法，数值求解带锥形尾喷管脉冲爆轰发动机的内外流场。计算结果表明：

（1）该三角形网格CE／SE方法可以有效捕捉到直管和变截面管中的爆轰波、激波等强间断，很好地刻画了PDE内外流场的细节。

（2）将锥形扩张管作为脉冲爆轰发动机的喷管而不填充燃料时，扩张管中压力峰值大幅下降，但马赫数峰值增加，发动机的排气速度增加，此时PDE可以得到更大的燃料比冲，发动机对燃料的利用率提高。

（3）而将锥形扩张管与爆轰管当成整体全部填充燃料时，锥形管中流场更简单，压力峰值比爆轰管中略有下降，此时发动机充分利用所有空间，瞬时推力和冲量增大，总体推进性能提高。

［1］BROPHY C M，DAUSEN D F，SMITH L R，et al.Fluidic nozzles for pulse detonation combustors［R］.AIAA Paper，2012－1035.

［2］NOBUYUKI T，YUICHIRO K，HAYASHI A K，et al.Numerical study and performance evaluation for pulse detonation engine with exhaust nozzle［R］.AIAA Paper，2009－5315.

［3］秦亚欣，于军力，高歌.脉冲爆震发动机喷管性能数值分析［J］.航空动力学报，2010，25（2）：366－372.

［4］WANG Z W，YAN C J.Experimental investigation of nozzle effects on a two－phase valveless air－breathing pulse detonation engine［R］.AIAA Paper，2008－991.

［5］李建中，王家骅，王春，等.共用尾喷管多管脉冲爆震发动机数值模拟研究［J］.空气动力学学报，2008，26（1）：96－100.

［6］CALDWELL N，GUTMARK E，HOKE J，et al.Investigation of fundamental processes leading to pulse detonation engine／ejector thrust augmentation ［R］.AIAA Paper，2008－116.

［7］曾昊，何立明，章雄伟，等.喷管收敛－扩张角对爆震发动机性能影响分析［J］.推进技术，2011，32（1）：97－102.

［8］KAILASANATH K.A review of research on pulse detonation engine nozzles［R］.AIAA Paper，2001－3932.

［9］CHANG S C.A new approach for constructing highly stable high order CESE schemes［R］.AIAA Paper，2010－543.

［10］CHENG G C，VENKATACHARI B S，CHANG C L，et al.Comparative study of different numerical approaches in space－time CE／SE framework for high－fidelity flow simulations［J］.ComputersandFluids，2011，45：47－54.

［11］WANG G，ZHANG D L，LIU K X，et al.An improved CE／SE scheme for numerical simulation of gaseous and two－phase detonations［J］.ComputersandFluids，2010，39：168－177.

［12］LI X，WENG C S.2－D viscous MHD simulation of plasma armatures by the CE／SE method［J］.Chinese ScienceBulletin，2009，54（10）：1641－1647.

［13］刘建文，赵书苗，钟诚文，等.CE／SE方法在多管爆轰流场并行计算中的应用［J］.爆炸与冲击，2008，28（3）：229－235。

［14］翁春生，王浩.计算内弹道学［M］.北京：国防工业出版社，2006：317－326.

［15］马丹花，翁春生.爆震管内扰流片对爆震波影响的数值分析［J］.推进技术，2011，32（3）：425－430.

［16］CHANG S C，WANG X Y，CHOW C Y.The spacetime conservation element and solution element method：a new high－resolution and genuinely multidimensional paradigm for solving conservation laws［J］.JournalofComputationalPhysics，1999，156：89－136.

［17］VENKAT E T，ANTHONY J D，NOBUYUKI T，et al.Performance of a pulse detonation engine under subsonic and supersonic flight conditions［R］.AIAA Paper，2007－1245.