标准试验发动机燃烧时间的自动计算方法*

张 荣，师 尧，王 刚，庞昀昭

(中国航天科技集团公司四院41所,西安 710025)

0 引言

固体火箭标准Ф127 mm发动机试验的主要目的是获取药柱燃速，进而预示和评估全尺寸发动机的燃烧性能。而燃速由固定参数肉厚与燃烧时间确定，因此燃烧时间求解的准确性直接关系到燃速这一重要性能指标的正确性，是评定发动机内弹道性能及其预示精度的关键。

固体火箭发动机燃烧时间tb是药柱燃烧表面开始点燃(通常规定燃烧室压强上升到300 kPa的时刻)到肉厚燃烧结束的时间间隔。目前，双切线法由于其使用的便捷性是固体发动机试验中使用最为广泛的一种燃烧时间确定方法，具体为：分别在压强-时间曲线的平稳段与下降段做两条切线,画出两条切线夹角的平分线,其与压强-时间曲线交点在时间轴上的投影定为燃烧终止时刻，如图1所示。

双切线法的优点是直观、简便、易操作，但是必须依托人工进行操作，可重复性差，存在较大的人为误差。尤其对于图2中这种没有明确的平稳段和下降段的曲线，难以确定切线绘制方式。如图3所示，3种切线绘制方式得到了不同的结果，导致数据处理的准确性降低。

由于标准Ф127 mm发动机试验量与日俱增，双切线法必须通过人工取点画切线的方式，无法实现试验数据处理的自动化，影响试验效率的提升。因此文中从数据处理角度进行探讨，提出一种不受人为因素影响的自动判定燃烧终止时间算法。

图2 3种拐点处特殊的曲线形状

1 燃烧时间的自动计算方法

固体火箭发动机燃烧时间tb主要通过对发动机燃烧室压强信号处理得到。自动求解燃烧时间的方法有很多种，例如：曲率最大法[1-3]、燃速公式法、计算机扫描法[4]、冲量系数法[5]、压强系数法[6]、线性加权极值法等。

燃速公式法是根据公式r=aPn计算燃速，其中发动机燃速系数a和发动机压强指数n需要大量的试验进行确定，不仅操作麻烦，而且当装药不完全燃烧时，此方法不可行。计算机扫描法将地面试验中计算机采集的压强-时间数据进行模拟显示，判断显示的压强-时间曲线是否有明显的拐点，该方法对于拐点不明显的曲线不可行。冲量系数法与压强系数法需要通过相同型号以往的大量数据获得相应的计算系数，在实际使用中比较麻烦。

图3 尾部三段型的3种切线绘制方式

曲率最大法是一种较新的燃烧时间自动求取方法，通过对压强曲线求取一阶导数、二阶导数以及曲率值，最终曲率最大值所对应的时间即为燃烧时间终点。这种方法可以对绝大部分发动机压强信号进行求解，也不需要以往的数据作为基础，方便快捷，适用性强，唯一的缺点是由于真实的发动机燃烧室压强-时间曲线是由离散数据组成，数据不连续也不光滑，不存在一个连续的一阶导数。目前有利用三次样条曲线拟合实际压强-时间曲线的方法，但是由于拟合精度不够，导致求得的燃烧时间与目前普遍使用的双切线法存在一定差距。

文中在曲率最大法的基础上提出了一种依托实际压强-时间曲线来确定标准Ф127 mm发动机燃烧时间的自动算法，不需要拟合样条曲线，直接在实际压强-时间曲线上对燃烧时间进行求解，并且为了保证以往燃速数据结果的延续性，对算法进行了一系列改进以使其最大程度地接近双切线法处理结果。

2 标准Ф127 mm发动机燃烧时间的曲率自动算法

2.1 算法理论依据

固体发动机内弹道曲线的形状取决于燃气生成率和排出率的关系。当生成率大于排出率时，压强上升；当二者平衡时，压强稳定；当燃气生成率小于燃气排出率时，压强下降。装药燃烧终止瞬间，燃面突然消失或突降，燃气生成率瞬间变为零或突降，压强曲线必然发生急剧转折。内弹道压强-时间曲线信息中包含着终燃点信息，曲线变化最快的地方，即曲率最大值处可认定为药柱燃烧终止时刻。

2.2 算法思路

针对标准Ф127 mm固体火箭发动机燃烧室实际采集得到的压强-时间曲线进行求解，首先确定燃烧终止时间所在的区间；然后对该区间内的压强-时间曲线进行归一化处理；接着给定一个合适的时间步长，从所确定区间的起始位置开始逐点向后每隔该时间步长进行微分求导，如图4所示，为从点1开始按步长为4个采样间隔进行的逐点求导，并按式(1)计算该处的曲率，以曲率最大位置为燃烧终止时间点。

(1)

其中：K为某点曲率；P′(t)为某点压强的一阶导数；P″(t)为某点压强的二阶导数。

图4 按时间步长为4个采样间隔进行求导

该算法中的关键在于如何在给定区间内完成压强-时间曲线的“归一化”以及确定一个合适的“求导时间步长”。

2.2.1 压强-时间曲线的“归一化”

压强-时间曲线的纵轴代表压强值，如图5(a)所示，若按照实际物理量进行求导计算曲率，计算结果都处于较大的量值，不便于进行比较，因此，对给定区间内的纵轴按照横轴进行等比例变换。在归一化之前需先选定归一化区间，该区间必须包含燃烧终止时刻，可以肯定燃烧终止时刻一定处于最大压强对应时刻与工作结束时刻之间，因此，选定归一化区间为[TPmax，Ta]，归一化处理的具体算法为：

PB=P×((Ta-TPmax)/Pmax)+TPmax

(2)

其中：PB为变化后的压强点值；P为实际压强点值；Ta为发动机工作结束时刻点；TPmax为平稳段压强最大值对应的时刻点；Pmax为平稳段压强最大值。

归一化之后，纵轴的最大值和最小值与横轴时间值一致，范围完全相同，纵轴的区间在数值上也是[TPmax，Ta]，如图5(b)所示。

图5 压强-时间曲线的归一化处理

2.2.2 求导时间步长

理想情况下，若压强-时间曲线是一条光滑连续、可用具体函数进行表征的曲线，那么根据2.1节中的理论依据，可直接通过计算曲率最大值的位置来确定燃烧终止时刻，但实际压强-时间曲线是由离散点组成的，且采样点由于电信号的干扰，在小范围内会有一定的波动，既不连续也不光滑，无法直接对曲线进行求导，如图6所示。算法引入时间步长求导的方法，利用插值法对曲线上每隔一个固定的求导时间步长的两个时间点根据式(3)进行求导计算，计算结果P′(t)作为[t1，t2]中点时刻(t1+t2)/2的一阶导数。使用相同的方法根据式(4)对一阶导数P′(t)逐点向后进行求导计算，得到二阶导数P″(t)，之后利用式(1)计算各点曲率及最大值位置。

P′(t)=(P(t2)-P(t1))/(t2-t1)

(3)

P″(t)=(P′(t2)-P′(t1))/(t2-t1)

(4)

图6 离散状的压强-时间曲线

对于一条光滑连续的曲线，根据式(5)一阶导数的定义，求导时间步长越短，得到的结果越准确。但对于实际离散的压强-时间曲线，求导的步长太短或太长获得的结果都与双切线法的结果相差较大，须确定一个合适的时间步长来保证处理结果的准确性。

(5)

为确定合适的时间步长，对涵盖高、中、低燃速近20种型号500发标准Ф127 mm发动机试验数据进行了数据分析验证，将其结果与双切线法处理结果进行对比。验证分为两个阶段，第一阶段将时间步长分别选取为0.005 s，0.01 s，0.02 s，…，0.1 s，记录所得燃烧时间与双切线法最接近的时间步长及其所对应的发动机个数，并做统计结果如图7所示。从图中并未发现合适的时间步长选取规律，但是通过进一步分析，发现步长与工作时间ta有一定的关系，时间步长在0.02～0.04倍的工作时间范围内，所得结果与双切线法结果更接近。第二阶段将时间步长选取为0.005ta，0.01ta，0.015ta，0.02ta，…，0.06ta，记录所得燃烧时间与双切线法最接近的时间步长及其所对应的发动机个数如图8所示。结果证明选取0.025ta作为求导步长效果最好。因此初步选取0.025ta作为微分求导步长进行一阶导数与二阶导数求解。

图7 最佳时间步长分布规律

图8 最佳时间步长占工作时间比例分布规律

2.3 算法实例

对某标准Ф127 mm发动机应用此算法进行燃速计算，使用PXI测试系统对该发动机燃烧室压强信号以2 000点/s的采样速率进行扫描采集，发动机燃烧室压强信号P，点火信号时刻点为t0=9.157 0 s，如图9所示。

图9 某标准发动机压强-时间曲线

算法流程如下：1)计算得到压强最大值Pmax=13 031.17 kPa，对应时刻点TPmax=9.580 5 s，工作结束时刻点为Ta=12.260 5 s，确定计算区间为[9.580 5 s，12.260 5 s]；

2)将[9.580 5 s，12.260 5 s]区间内的压强数值P进行归一化处理后为PB；

3)选取微分求导的时间步长Δt=(12.260 5-9.580 5)×0.025=0.067 s；

4)从9.580 5 s开始，以0.067 s为步长，对变换后的PB曲线逐点向后进行求导计算，得到一阶导数P′，对应区间为[9.614 s，12.227 s]；

5)对一阶导数P′在区间[10.691 0 s, 13.030 0 s]内再次逐点向后进行求导计算，得到二阶导数P″，对应区间为[9.647 5 s，12.193 5 s]；

6)利用曲率计算公式计算曲率K，得到最大曲率值为51.455 3，该点对应的采集时刻值为12.046 5 s。

7)得到燃烧时间tb为12.046 5-9.157 0=2.889 5 s。

3 曲率自动算法改进

将自动算法得到的燃烧时间tb与双切线法处理结果对比，结果显示即使选取效果最好的0.025ta作为求导步长，两者的平均时间差仍有0.013 9 s，结果相差较大，相较于双切线法，曲率自动算法计算得到的燃烧时间偏小。

为了提高曲率自动算法处理结果与以往试验处理结果的延续性与一致性，需要对算法进行改进，使两者的处理结果尽可能接近。通过大量的试验数据统计分析，在自动算法中，当采用同一时间步长，若归一化区间设定不同，得到的最大曲率位置也会有所差异，且区间越小，得到的燃烧时间越大。

目前在应用双切线法进行处理时，通常只选取下降段拐点附近的一小段数据，对其进行归一化处理后绘制切线，该区间包含燃烧终止时刻，但其具体范围是人为随机选取的，导致燃烧时间求取不稳定。曲率自动算法的归一化区间为[TPmax，Ta]，相较于双切线法较大，所以计算得到的燃烧时间偏小。为了减小两者的时间差须将曲率自动算法的归一化区间尽可能接近双切线法。由于试验数据在没有经过分析之前，无法自动选取一段包含tb且与双切线法大致一致的分析区间，但可以肯定tb一定在最大压强时刻与工作结束时刻之间，因此，首先确定数据的归一化区间为[TPmax，Ta]，得到曲率最大位置tb1，然后将归一化区间再次设置为[tb1-3Δt,ta]，其中Δt为求导时间步长，得到曲率最大位置tb2，对计算结果进行多次迭代，直到计算结果不再变化，确定其为tb。根据经验，一般最多迭代3次即可。这种区间选取方法只与Δt有关，相较于双切线法，区间选取更为确定，处理结果更加稳定，不受人为因素干扰。利用新的区间求取方法对原曲率自动算法进行改进，并对2 600余发标准Ф127 mm发动机试验数据重新进行数据分析验证，结果显示改进后的曲率自动算法选取时间步长为0.02ta效果最好，与双切线法结果的平均时间差为0.003 6 s，标准差为0.007 5 s，结果统计如表1所示。

表1 自动算法与双切线法的燃烧时间差

4 结论

提出了一种依托实际压强-时间曲线的曲率最大位置自动算法，以此自动确定标准Ф127 mm发动机燃烧时间。该算法引入迭代算法，通过大量数据统计，当选取发动机工作时间的2%作为求导时间步长，其处理结果最接近双切线法，得到的燃烧时间差在5 ms之内，满足了试验数据结果的延续性。该算法相比于双切线法具有可实现自动化、在参数设置不变的情况下结果数据唯一的特点，提高了固体标准发动机试验数据处理的可靠性与效率。