过载条件下固体发动机声能共振规律研究①

2023-04-26宋儒儒张翔宇甘晓松赵天泉

固体火箭技术 2023年1期

宋儒儒,张翔宇,甘晓松,赵天泉

(1.中国航天科技集团有限公司四院四十一所,西安 710025;2.中国航天科技集团有限公司第四研究院,西安 710025)

0 引言

固体火箭发动机用于大型运载火箭助推器和各类导弹武器已超过半个世纪,但不稳定燃烧现象在各种尺度的战术导弹发动机及分段式发动机上时常发生,是发动机研制中不可忽略的问题[1-6]。发动机的不稳定燃烧本质上是其内部各种增益/阻尼因素之间相互作用导致的压强周期性振荡。粒子阻尼作为其中的一项阻尼因素,占总阻尼的1/3以上,远远大于结构阻尼和气相阻尼[1]。

粒子阻尼理论研究中使用最广的是微粒松弛理论,当惰性微粒的动力弛豫时间与声振荡的特征时间相等时,惰性微粒的阻尼效应最大。CAI等[7]基于欧拉-拉格朗日法对固体火箭发动机内两相流与声能的交换机理进行了数值模拟,充分考虑了声场、湍流散布、粒子碰撞与聚合等复杂流场条件下气相和粒子相之间的耦合作用,指出粒子的动力驰豫时间、热驰豫时间和声场的特征时间之间的关系对发动机内流场两相流的耦合起着非常重要的作用。西北工业大学刘洋、何国强等[8]提出了一种可收集燃烧室中凝聚态粒子的方法,设计了粒子收集试验装置,通过改变试验装置的收敛半角和试验状态参数可以模拟真实发动机中的粒子聚集状态。彭小波等[9]改进了燃烧室中凝相燃烧产物的试验收集装置,并开展了不同聚集状态下凝相燃烧产物的燃烧特性分析。金秉宁等[10]针对不同初始铝粉粒度的含铝复合推进剂,对其燃烧产物进行了粒子阻尼特性和铝粒子分布燃烧响应特性的试验研究,得到了关于粒子阻尼的大小及预估、粒子对分布燃烧的影响方面的规律。

随着对粒子阻尼的深入研究,可以从发动机飞行状态变化和燃烧室声腔声压振荡等动态的角度理解粒子阻尼。固体火箭在做机动飞行时,横向过载往往会改变燃烧室内凝相粒子的空间分布[11-12],从而影响粒子的阻尼特性,甚至会导致不稳定燃烧。西北工业大学刘佩进课题组[13]对此展开研究,研究结果表明在大过载条件下,燃烧室内大尺寸凝相颗粒会在极短的时间内在局部形成聚集带和真空带,导致粒子阻尼减小,这对发动机稳定性是十分不利的。同时,在极限环振荡的发展过程中,由于声对粒子的操纵作用,可能导致粒子沿发动机纵向方向重新分布,引起粒子阻尼特性的改变,并参与压强振荡非线性增长过程中的增益与阻尼博弈[14]。北京理工大学王宁飞团队[15-16]使用脉冲衰减法对粒子阻尼进行数值仿真研究,结果表明,由脉冲衰减法计算所得的结果与T型燃烧器的实验结果基本吻合,验证了该仿真方法在计算粒子阻尼时的合理性与有效性。

本文从声学角度出发,用声场介质的声衰减系数代替粒子阻尼,探究固体火箭发动机在不同飞行过载条件下声学特性的变化规律,以声学响应传递函数来评估发动机的声不稳定性。

1 粒子阻尼理论

粒子阻尼作用是由弛豫现象引起的。由于气体的粘性和粒子的不可压缩性,使粒子与气体达到相同的速度、温度需要一定的时间,所以存在动力弛豫过程与热力弛豫过程。

定义动力弛豫时间τv,即粒子速度达到与气体速度平衡所需要的时间。

(1)

式中ms为单个粒子的质量,kg;μ为气体的动力粘度系数,kg/(m·s);σ为粒子半径,μm;ρs为粒子自身密度,kg/m3。

定义热力弛豫时间τt,是粒子温度达到与气体温度平衡所需要的时间。

(2)

式中Pr为普朗特数,Pr=μcp/λg;cs为粒子的比热容,J/(kg·K);cp为气体比定压热容,J/(kg·K)。

用压强振幅的时间衰减率(即衰减系数)αp来表示粒子阻尼的大小,粒子阻尼系数可表示为

(3)

以上各式中动力粘度系数μ由Sutherland定律确定:

式中μ0为参考粘度,1.716×10-5kg/(m·s);T为燃气静温;T0为参考温度,273.11 K;S为Sutherland常数,110.56 K。

2 声能共振数值计算方法

2.1 总体思路

本文通过将粒子阻尼转化为声场介质的声衰减系数,以声学响应传递函数评估发动机声不稳定性,研究过程中,不考虑其他阻尼因素影响,思路如下:

(1)对发动机几何模型划分计算域,并进行声腔模态提取及喷管阻尼计算,以喷管阻尼计算相应粒子阻尼范围,确定凝相粒子粒径与浓度的初始参数;

(2)进行两相流计算,并统计各个计算域的粒子浓度,将各计算域浓度结果带入式(3)计算相应粒子阻尼值;

(3)用声场介质的声衰减系数代替粒子阻尼值,表征粒子阻尼对声腔声学特性的影响;

(4)在各计算域添加相应声衰减系数,计算发动机声腔的声学特性,分析声学响应传递函数的变化规律。

2.2 计算模型

本文选用某发动机工作后期的简化模型作为研究对象。由于喷管喉部下游为超声速流动,声波不再反射,因此在进行仿真分析时选取喷管喉部上游部分结构,几何模型及主要参数如图1所示。为统计过载条件下发动机内凝相粒子浓度分布,对模型进行计算域划分,并参照笛卡尔坐标系方向对各区域编号,如图2所示。为便于表述,将所有计算域划分为4部分:第一部分为坐标系(+Y+Z)区域,包含1～37号计算域,记作A;第二部分为坐标系(+Y-Z)区域,包含2～38号计算域,记作B;第三部分为坐标系(-Y-Z)区域,包含3～39号计算域,记作C;第四部分为坐标系(-Y+Z)区域,包含4～40号计算域,记作D。

图1 某发动机工作后期几何构型

图2 不同计算域编号

提取声腔模态,模态分布如图3所示,模态数值见表1。参考稳态波衰减法计算该模型的喷管阻尼,在发动机头部添加单极子声源,声源频率设置为与声腔固有频率相同,本文取162.094 Hz,幅值为0.1 kg/s,并在声腔内建立起稳定的驻波后关闭声源。声腔内压力振荡幅值以指数形式衰减,将衰减过程的振荡幅值绘制在半对数坐标系中取包络线并进行线性拟合,所得拟合直线斜率即为衰减系数,声腔内声压衰减过程以及拟合直线斜率如图4所示。

(a)Sound pressure attenuation process

表1 声腔各阶模态值

(a)First natural frequency

对声压衰减过程进行直线拟合,斜率k=-25.22,即喷管阻尼αN=25.22 s-1,分析发动机内部各种声能阻尼占比[1],αN∶αp≈3∶2,因此该发动机模型中粒子阻尼在16 s-1左右。经反复校核,在粒径为48 μm,浓度为15%时,通过公式(3)计算所得的粒子阻尼为16.8 s-1,符合两种阻尼在总阻尼中的占比,可进行后续计算。

2.3 计算工况参数

气相控制方程采用N-S方程,颗粒相运动轨迹则由拉格朗日方程确定,获得粒子速度和浓度的空间分布。稳态气相流场计算采用质量流量入口,入口流量为6.072 kg/s,介质温度为3235.88 K,出口为压力出口边界,喉部出口静压为4 035 328 Pa,静温为300 K,壁面选择无滑移边界条件,使用基于压力的求解器,湍流模型选用SSTk-ε模型,对控制方程进行二阶迎风格式离散,计算结果如图5所示。

(a)Gas phase pressure distribution

由固体发动机原理[17]可知,燃烧室末端燃气速度系数一般为0.2～0.5,即燃气速度为200～500 m/s。以长度为3 m的发动机为例,据此估算粒子在燃烧室内最长停留时间为30 ms,因此本文计算时间取30 ms。凝相粒子的初始速度取为气相入口速度的0.4倍[18],壁面条件为弹性反弹,凝相粒子粒径、浓度及过载参数见表2。

表2 凝相粒子初始参数

2.4 粒子阻尼的转换

声场中介质的阻尼效应通常使用声衰减系数表征,即在声速表达式中赋加虚部值,如式(4)所示。

(4)

通过拟合声衰减系数Y与粒子阻尼系数αp的关系,可以将粒子阻尼转化为对应的声衰减系数。对Y分别取10组值,利用稳态波衰减法计算声腔内阻尼值,计算方法同2.2节,计算所得阻尼值减去喷管阻尼值即为粒子阻尼αp,拟合结果如图6所示。

图6 Y与αp 的关系拟合

粒子阻尼αp与声衰减系数Y关系式如式(5):

Y=0.152 6αp+0.253 128

(5)

2.5 声学仿真计算

采用直接频率响应分析方法,在发动机头部添加球状声源,声源幅值为2000 Pa,并在每个计算域中心位置添加一个声压频率监测点,同时定义声学响应传递函数Pf来表征压力振荡程度,见式(6),依据该参数来评估发动机出现声不稳定的趋势。

(6)

式中p为测点声压值,Pa;ps为声源声压值,Pa。

3 结果及分析

3.1 各计算域粒子浓度分布规律

设置轴向过载为15g,横向过载分别为3g、6g、9g、12g、15g,过载施加时间同2.3节仿真计算时间一致,即从0 s开始直至仿真计算结束。探究30 ms时刻下粒子空间分布规律以及各工况下声学响应情况,同时以未加过载即横向过载与轴向过载均为0g的工况作为对照,各工况下粒子浓度的空间分布如图7、图8所示。

(a)Transverse overload=0g (b)Transverse overload=3g

(a)Transverse overload=3g (b)Transverse overload=6g (c)Transverse overload=9g

采用侧向加质,轴向及横向过载方向分别与X、Y轴正方向保持一致。结合图7～图9,在施加轴向过载的工况下发动机内所有粒子的速度方向与X轴相同,不断向喷管喉部运动。发动机头部基本不存在粒子轨迹,形成粒子真空区,喷管喉部的粒子不断聚集,形成粒子聚集区。

在横向过载的影响下,粒子不断沿Y轴矢量方向在发动机近壁面区域聚集,如图8所示。并随着横向过载由3g不断增加到15g的过程中,粒子聚集区域中心沿Y轴的偏移量逐渐增大。当横向过载加至12g和15g时,粒子聚集区域缩小,粒子空间分布更加集中,且聚集区中心沿Y轴偏移较大,同时在Y轴负方向发动机近壁面区域形成较大的粒子真空区,如图8(d)、(e)所示。

在横向、轴向过载的共同作用下,发动机喷管附近A、B区域粒子浓度明显高于C、D区域,如图9所示。

(a)Domain 1～37 (b)Domain 2～38

3.2 网格无关性验证

为检验仿真计算的合理性,对模型绘制800 000、1 000 000、1 200 000数目的网格,进行网格无关性验证,网格模型如图10所示。以表2中case 6为例,分别计算各套网格模型中Domain 1～37计算域的粒子阻尼,结果如图11所示。

(a)80×104

由图11可知,不同网格数目下,粒子阻尼值随计算域监测点的变化趋势基本一致,即阻尼计算结果对网格依赖性极小。在计算误差允许范围内,可证明网格数量对数值计算结果无影响。

图11 不同网格数目下的阻尼值

3.3 粒子浓度分布对声学响应的影响

将3.1节各区域计算所得的粒子浓度带入式(3)计算对应粒子阻尼值,计算所得阻尼值结果带入式(5),可计算相应声衰减系数,将其添加至声场中进行声学响应计算。

对未加过载的工况进行频率响应分析,取37、38、39、40等四个区域中心添加监测点,绘制频率响应曲线如图12所示。可以看出四个监测点的频率响应曲线重合,计算结果一致,且发动机内模态频率分布以一阶频率为主,故后续对不同工况下各监测点声学响应的分析主要在一阶频率下进行。

图12 零过载工况下各测点频率响应曲线

发动机模型具有很好的对称性,故只分析A、D区域测点的声学响应。不同工况下各计算域一阶振荡振幅分布如图13所示,发动机头部和喷管喉部测点的声压振幅最大,发动机中间位置振幅较小,符合一阶轴向振型分布特点,且发动机头部声学响应略高于喷管喉部。

(a)Domain 1～37 (b)Domain 4～40

对发动机头部与喷管喉部区域进行声学响应分析,如图14所示,可见所有监测点的声学响应随横向过载的增加不断增大。图14(a)是Domain 1和Domain 37的监测点声学响应分布情况,Domain 37为粒子聚集区,并随着横向过载的不断增加,该区域的粒子浓度不断增加,但声学响应也呈不断增大的趋势,并未出现由于粒子聚集使该区域的声学响应减小,这是因为横向过载导致粒子的整体空间分布极其不均匀,对声振荡的阻尼作用减小,导致声学响应增大。对比Domain 1和Domain 37的声学响应,横向过载相同时,Domain 1的声学响应高于Domain 37,符合粒子阻尼理论,即粒子浓度增加,粒子阻尼值增加,对声不稳定的抑制作用越强。图14(b)所示为Domain 4和Domain 40的监测点声学响应分布情况,在横向过载由3g不断增加到15g的过程中,声学响应呈现出不断增加的趋势,且相同过载情况下Domain 4的声学响应高于Domain 40。