正癸烷/空气旋转爆轰波形成与传播数值模拟

邵晓峰, 赵宁波, 刘世铮, 孟庆洋, 李雅军, 郑洪涛

(1.哈尔滨工程大学 动力与能源工程学院,黑龙江 哈尔滨 150001; 2.中国船舶集团有限公司第七〇三研究所,黑龙江 哈尔滨 150078)

由于等压燃烧过程熵增较大,现代布雷顿循环热力发动机的循环效率严重,受限于新材料技术的发展。旋转爆轰燃烧作为一种近等容燃烧技术,具有熵增小、释热快、自增压等优势,在航空、航天以及舰船动力推进领域均有广泛的应用前景[1-2]。

燃料属性是影响旋转爆轰波形成传播与燃烧室性能的关键因素之一。近年来,国内外有关该方向的研究主要集中于氢气[3]、丙烷[4]、乙烯[5]等高活性气态燃料方面,而对工程实际应用中的液态燃料研究相对较少。Bykovskii等[6-8]对两相旋转爆轰波的形成与稳定传播进行了实验探索,通过在氧化剂中添加氢气助燃的方式成功得到了航空煤油/空气旋转爆轰波,并得出在没有助燃剂加入的情况下,航空煤油和空气在常温下无法起爆。Kindracki[9]对液态煤油喷入燃烧室的冷态过程进行了实验研究,其研究结果表明,煤油液滴的索特尔平均直径约33～38 μm,能够在高温环境下迅速蒸发。在冷态实验基础上,Kindracki[10]通过添加少量氢气助燃的方法实现了煤油/空气旋转爆轰,并研究了硝酸异丙酯对提高煤油/空气两相爆轰灵敏度的影响。郑权等[11-13]以汽油为燃料,富氧空气为氧化剂,在环形阵列式旋转爆轰燃烧室内形成了旋转爆轰波,详细分析了燃烧室长度[11]、燃料喷注压力[12]和当量比[13]对爆轰波传播模态、传播速度、推力和比冲的影响。

为进一步了解气液两相旋转爆轰波的形成与传播机理,Meng等[14]对正庚烷/空气两相旋转爆轰进行了二维数值模拟,采用在正庚烷液滴中添加预蒸发正庚烷蒸气的方式来提高燃料活性,发现液滴粒径对爆轰波传播速度的影响较大,并且随着液滴粒径的增大,爆轰波速度会衰减约5%～30%。此外,Meng等[15]还研究了进气温度对正庚烷/空气两相旋转爆轰波传播特性的影响,指出随着进气温度的升高,液滴的蒸发率逐渐增加并在填充区与缓燃层之间形成高浓度的正庚烷蒸气层。Hayashi等[16]以JP-10为燃料,空气为氧化剂,数值模拟研究了不同液滴粒径(1～10 μm)对爆轰波传播演变规律的影响,结果表明,当液滴粒径过大时,爆轰波后会存在未燃尽的液滴而导致爆轰波熄灭。Sun等[17]研究了进口总温和燃料进口间距对辛烷/空气两相旋转爆轰波形成与传播过程的影响,指出爆轰波能够稳定传播的临界温度和燃料进口间距分别为600 K和10 mm。徐高等[18]在汽油/富氧空气两相旋转爆轰数值模拟中考虑了液滴的雾化破碎过程,并研究了不同燃料喷注压力和氧化剂填充比对爆轰波传播模态的影响,发现不同工况下爆轰波会呈现稳定单波、稳定双波、不稳定双波和不稳定单波4种传播模态。

针对上述研究现状和近年来液态燃料旋转爆轰发动机的快速发展需求,本文以液态正癸烷-空气为研究对象,采用三维数值模拟方法研究非预混旋转爆轰波的形成与传播过程,重点分析两相旋转爆轰波的形成演变特点和自持传播机理。

1 物理模型及数值计算方法

1.1 计算模型

旋转爆轰燃烧室(rotating detonation chamber, RDC)采用环缝-喷孔形式,其中空气由渐缩-渐扩型环缝进入RDC,燃料由均匀分布在环缝扩张段的90个燃料孔进行喷注,如图1(a)所示。此外,基于正癸烷/空气的胞格尺寸,RDC内半径设为Rin=68 mm,外半径设为Rout=78 mm,轴向长度设为L=85 mm,如图1(b)所示。

图1 非预混旋转爆轰燃烧室结构示意

1.2 计算方法

采用商业软件ANSYS Fluent进行数值模拟,气相部分的控制方程为:

(1)

(2)

(3)

(4)

p=ρRT

(5)

对于液相部分,假设液滴在破碎和蒸发过程均保持球形且液滴的内部温度保持一致。由于液滴的体积分数小于1‰,故忽略液滴之间的相互作用。液相部分的控制方程为:

(6)

dud/dt=(u-ud)/τr

(7)

(8)

液滴的破碎过程采用KH-RT破碎模型[19]进行模拟,该模型考虑了液滴表面的瑞利-泰勒(Rayleigh-Taylor, RT)不稳定性对液滴破碎过程的影响。液滴表面的最快扰动频率为:

(9)

式中:gt为液滴传播方向的加速度;ρg和ρl分别为气相和液相的密度;σ为液滴的表面张力系数。液滴表面相应的波数为:

(10)

液滴的破碎时间为:

tRT=Cτ/ΩRT

(11)

式中Cτ为RT破碎时间常数,值为0.5。破碎后的子液滴半径为:

rc=πCRT/KRT

(12)

式中CRT为破碎半径常数,值为0.1。

在液相的计算中,需要与气相进行耦合求解,SY为液相组分时SY=Sm,其他组分时为0。耦合求解的源项Sm、Su、Se在每个网格中分别表示为:

(13)

(14)

(15)

式中:Vc为每个网格的体积;Nd为网格数。

在求解过程中,对于气相部分,对流项采用三阶单调守恒迎风格式(monotonic upwind scheme for conservation laws,MUSCL)格式进行离散,物理通量采用对激波捕捉效果较好的迎风型矢通量分裂法(advection upstream splitting method,AUSM)进行分解,时间推进方式采用四阶龙格-库塔法。湍流模型采用标准k-ε模型[20],燃烧模型采用层流有限速率模型,化学反应机理采用Jia等[21]验证的总包反应,化学反应速率为反应速率常数和反应物浓度的乘积,其中反应速率常数由阿伦尼乌斯公式给出:

(16)

式中:A为指前因子,值为2.587×109s-1;b为温度指数,值为0;Ea为活化能,值为1.256×108J/kmol;F和O分别为燃料和氧化剂的浓度;α和β为对应的化学反应级数,值分别为0.25和1.5。

对于液相部分,式(6)～(8)通过一阶隐式欧拉法进行求解,式(13)～(15)与式(1)～(4)中的源项进行耦合求解。

1.3 边界条件

空气采用质量流量进口,液滴喷注过程利用两相DPM方法进行模拟,具体的边界参数如表1所示。燃烧室出口为压力出口,内外壁面均为绝热无滑移壁面,且假设液滴与壁面只发生反射而不被壁面捕捉。

表1 边界参数

1.4 点火方式

采用如图2(a)所示的燃料初始分布方式进行点火起爆,雷知迪等[22]验证了该方法的可行性并指出该方法能够实现单次点火建立稳定自持的旋转爆轰波。RDC被分割为2个区域(即RDC1和RDC2)但两者相互连通,高能点火区域(半径为5 mm、压力为2 MPa、温度为3 000 K的球形)位于RDC2一侧。点火过程时序如图2(b)所示,首先对RDC1和RDC2进行空气喷注,然后将燃料从RDC2一侧的燃料孔喷入燃烧室,并且Δt1=300 μs时RDC2点火,Δt2=60 μs时RDC1开始喷入燃料。

图2 点火方式示意和时序

1.5 网格无关性检验及模型验证

为确保数值结果的有效性,首先对3种不同网格尺度(0.2、0.4和0.6 mm)下的网格进行无关性检验。图3给出了监测点处(x=20 mm,y=0 mm,z=73 mm)压力随时间的变化曲线,从图中可以看出在网格尺度为0.4 mm时已经达到无关性要求。此外,在本文的条件下,正癸烷/空气的胞格尺寸为14.8 mm。网格尺寸为0.4 mm时,在一个胞格内存在超过1 000个网格,网格分辨率足够捕捉流场信息。因此本文选择网格尺度0.4 mm用于后续相关模拟分析。

图3 网格无关性检验

为验证所选数值方法的准确性,表2所示为数值模拟结果与试验值和理论值vCJ为1 830 m/s的对比。其中理论值由SD Toolbox[23]计算得到。从表2中可以看出,由于液滴蒸发破碎,燃料掺混不均匀等因素的影响,实验中两相爆轰波速度衰减达到30%以上。而与文献[24-25]数值模拟结果相对比,本文所得到的爆轰波速度衰减仅为10.7%,处于合理范围内,可以说明所采用的数值方法是可行的。

表2 数值模拟计算值与实验值对比结果

2 爆轰结果与讨论

本文基于前期对旋转爆轰的研究基础,拟从起爆、对撞、稳定传播3个阶段详细剖析正癸烷/空气非预混旋转爆轰波的形成与传播特性。

2.1 起爆阶段

图4给出了点火后不同时刻下,燃烧室x=20 mm截面处的压力和温度分布云图。由于本文所得到的爆轰波平均高度为37.5 mm,因此本文选择靠近燃烧室进口的x=20 mm截面即可以保证爆轰波形成过程的一致性。从图中可以看出,t=8 μs时,燃烧室内形成了一道逆时针传播的爆轰波(压力和温度峰值可达2 MPa、3 000 K)、一道顺时针传播的弱压力波;t=16 μs时,压力波传播至RDC1,由于缺乏燃料供应,其压力波阵面与火焰面发生明显解耦;t=36 μs时,顺时针传播的压力波穿过RDC1区域,并在RDC2区域新鲜混合物的支持下强度有所提高。同时RDC1区域内还观察到另一道逆时针传播的压力波,其对流场的扰动使得爆轰产物的压力和温度进一步提高。

图4 起爆后压力与温度分布云图

图5为t=36 μs时,燃烧室在x=20 mm截面处的压力与温度沿中环线(内外半径间的圆环线,θ)的变化。其中0°所在位置为RDC1和RDC2区域的右侧交界线,角度旋转方向为逆时针。从图中可以明显看出,燃烧室内只在270°附近存在压力峰值(约1.3 MPa)并与温度峰值(约2 800 K)强耦合现象,虽然在90°位置处压力和温度曲线同样存在间断,但压力峰值和温度峰值并未耦合,判断该处未形成爆轰波。此外,在2道波面之间(0°～90°和270°～360°),压力和温度均存在复杂波动,这主要是回传压力波对爆轰产物的扰动所导致。

图5 t=36 μs时刻沿中环线的压力和温度分布

2.2 对撞阶段

经2.1节起爆过程之后,爆轰波与压力波继续沿各自的传播方向运动,并在t=76 μs发生首次对撞。图6所示为爆轰波对撞阶段x=20 mm平面上的压力和温度分布。从图中可以看出,此时RDC中存在3道波面(其中1表示逆时针传播的爆轰波;2表示压力波;3表示回传压力波),并且在爆轰波1和压力波2的对撞点附近产生局部高压高温区,压力峰值达到2.8 MPa,温度峰值超过3 600 K。

图6 t=76 μs时刻的压力和温度分布云图

图7为t=120 μs时刻下x=20 mm平面上的压力、温度、气态正癸烷质量分数和反应速率云图。从图7(a)中可以看出,爆轰波1和压力波2在经过t=76 μs时刻的对撞后仍沿原方向继续传播,Bluemner等[29-30]分别在实验和数值模拟中发现过类似的现象。从图7中还发现,压力波2与回传压力波3发生对撞时,对撞点处压力低于1 MPa,温度低于2 000 K,并且化学反应速率低,因此判断此处未形成爆轰波。对于爆轰波1,从气态正癸烷质量分数云图中可以看出,爆轰波前存在大量的可燃气体且在爆轰波面处化学反应速率较高,因此认为可燃气体燃烧释放能量是维持爆轰波传播的主要原因。

图7 t=120 μs时刻的不同参数分布云图

图8为t=120 μs时刻x=20 mm平面上沿中环线的压力与温度分布。此时燃烧室中仅存在一个压力峰值(约为2 MPa)且与温度峰值耦合。压力波2与回传压力波3对撞点处的压力变化不大且未出现压力峰值与温度峰值耦合,即没有出现火焰面跟随压力波传播的现象。由此可以推断,此时燃烧室中仅存在一道爆轰波。

图8 t=120 μs时刻沿中环线的压力和温度分布

图9为不同时刻下爆轰波1与压力波2的对撞情况。发现回传压力波3在t=120 μs对撞结束后逐渐衰减并消失,燃烧室中只存在爆轰波1和压力波2。此外,爆轰波1与压力波2共经历2次对撞,分别发生在t=176 μs和t=260 μs;压力波2在对撞后逐渐衰减,并在t=364 μs时完全消失,此时爆轰波1以单波模态传播。

图9 不同时刻的压力与温度云图

2.3 稳定传播阶段

图10所示为监测点处(x=15 mm,y=0 mm,z=69 mm)压力随时间(600 μs

图10 监测点处的压力变化

图11对比了爆轰波稳定传播时燃烧室内壁面的压力、温度和气态正癸烷质量分数。值得指出的是,两相旋转爆轰波与气相旋转爆轰波结构的最大区别在于爆轰三叉点处存在未反应气流区[31],未反应气流区中有大量未燃尽的气态正癸烷,这会降低爆轰波面和爆轰波后部分区域温度。

图11 不同的流场参数云图

为了深入揭示未反应气流区的形成原因,图12进一步给出了爆轰波后的压力和马赫数云图并叠加液滴直径D分布。从图中可以看出,在环缝扩张段出现了压力波回传现象,并且在压力波的影响下,环缝扩张段出口背压提高,使得局部马赫数降低。假设液滴射流角度定义为θ,从图12可以看出,液滴在来流高速气流和回传压力波的作用下迅速破碎,大部分液滴直径小于5 μm可以快速参与燃烧,但液滴射流角度在回传压力波的作用下可能会发生改变,进而使得燃料与氧化剂的掺混过程更加复杂[32]。

图12 压力、马赫数以及液滴直径分布云图

图13进一步分析了环缝扩张段周向截面上液滴射流角度与马赫数之间的关系,其中爆轰波面处于阶段1和阶段3的交界面,特征为来流马赫数和射流角度均有明显的突变。从图13中可以看出,在爆轰波前后,马赫数与射流角度存在动态调节机制。阶段1为压力波回传造成空气进气阻塞,扩张段马赫数降低,为亚音速状态,射流角度较大,爆轰波后为逆向喷注状态(θ>90°);阶段2,随着爆轰波后的等熵膨胀过程,压力迅速下降,空气进气逐渐恢复,表现为扩张段马赫数逐渐上升,射流角度逐渐下降,恢复为顺喷状态(θ<90°);阶段3,空气进气完全恢复,表现为马赫数和射流角度基本稳定。

图13 沿周向的射流角度和马赫数变化

从以上分析可以看出,压力波回传会改变来流马赫数,从而影响液滴射流角度,使得燃料掺混不均匀,燃料无法被完全消耗;同时,液态燃料在燃烧之前需要经过破碎蒸发,化学反应区较宽,这也会导致部分燃料难以在较小的空间内充分反应。由此可以推断,燃料不均匀掺混及其不完全燃烧共同促进了未反应气流区的形成。

图14给出了684 μs≤t≤712 μs不同时刻下RDC内壁面爆轰波的压力云图。从图14中可以看出,在t=684 μs时,未反应气流区已经基本形成;由于未反应气流区存在大量未燃反应物,其可在t=688 μs时被点燃形成局部爆炸点,并且随后爆炸点沿轴向发展成为轴向爆轰波2和3;然后轴向爆轰波2和3分别与爆轰波1和斜激波对撞,爆轰波3传播的方向由于缺乏可燃气体支撑,在对撞后很快熄灭,燃烧室中只存在爆轰波1和轴向爆轰波2;随着轴向爆轰波2不断向上游传播,爆轰波1不断被压缩,所占比重逐渐降低,最终在燃烧室中只存在轴向爆轰波2。此外,在轴向爆轰波2发展的同时,原未反应气流区中的燃料被消耗,并逐渐消失。燃料填充区内的新鲜燃料继续向三叉点附近聚集,当燃料聚集到一定程度时,流场局部温度和压力出现明显下降,形成新的未反应气流区,如图14(h)所示,未反应气流区形成时间的平均周期为28 μs。由此可以看出,未反应气流区的周期性起爆是两相爆轰波能够自持传播的重要原因。

图14 不同时刻的压力云图

3 结论

1) 在爆轰波的形成过程中,存在多次爆轰波与压力波之间的对撞,在该过程中压力波逐渐消失,爆轰波以单波模态传播。

2) 在爆轰波自持传播过程中,燃料的不均匀掺混以及液态燃料的不完全燃烧促进了未反应气流区的形成,未反应气流区的周期性起爆是两相旋转爆轰波能够自持传播的重要原因。

3) 在旋转爆轰燃烧室中获得了自持传播的两相旋转爆轰波,模拟结果为两相旋转爆轰发动机的发展提供了理论支撑。