標(biāo)準(zhǔn)試驗(yàn)發(fā)動(dòng)機(jī)燃燒時(shí)間的自動(dòng)計(jì)算方法*

張 榮，師 堯，王 剛，龐昀昭

(中國(guó)航天科技集團(tuán)公司四院41所,西安 710025)

0 引言

固體火箭標(biāo)準(zhǔn)Ф127 mm發(fā)動(dòng)機(jī)試驗(yàn)的主要目的是獲取藥柱燃速，進(jìn)而預(yù)示和評(píng)估全尺寸發(fā)動(dòng)機(jī)的燃燒性能。而燃速由固定參數(shù)肉厚與燃燒時(shí)間確定，因此燃燒時(shí)間求解的準(zhǔn)確性直接關(guān)系到燃速這一重要性能指標(biāo)的正確性，是評(píng)定發(fā)動(dòng)機(jī)內(nèi)彈道性能及其預(yù)示精度的關(guān)鍵。

固體火箭發(fā)動(dòng)機(jī)燃燒時(shí)間tb是藥柱燃燒表面開始點(diǎn)燃(通常規(guī)定燃燒室壓強(qiáng)上升到300 kPa的時(shí)刻)到肉厚燃燒結(jié)束的時(shí)間間隔。目前，雙切線法由于其使用的便捷性是固體發(fā)動(dòng)機(jī)試驗(yàn)中使用最為廣泛的一種燃燒時(shí)間確定方法，具體為：分別在壓強(qiáng)-時(shí)間曲線的平穩(wěn)段與下降段做兩條切線,畫出兩條切線夾角的平分線,其與壓強(qiáng)-時(shí)間曲線交點(diǎn)在時(shí)間軸上的投影定為燃燒終止時(shí)刻，如圖1所示。

雙切線法的優(yōu)點(diǎn)是直觀、簡(jiǎn)便、易操作，但是必須依托人工進(jìn)行操作，可重復(fù)性差，存在較大的人為誤差。尤其對(duì)于圖2中這種沒有明確的平穩(wěn)段和下降段的曲線，難以確定切線繪制方式。如圖3所示，3種切線繪制方式得到了不同的結(jié)果，導(dǎo)致數(shù)據(jù)處理的準(zhǔn)確性降低。

由于標(biāo)準(zhǔn)Ф127 mm發(fā)動(dòng)機(jī)試驗(yàn)量與日俱增，雙切線法必須通過人工取點(diǎn)畫切線的方式，無法實(shí)現(xiàn)試驗(yàn)數(shù)據(jù)處理的自動(dòng)化，影響試驗(yàn)效率的提升。因此文中從數(shù)據(jù)處理角度進(jìn)行探討，提出一種不受人為因素影響的自動(dòng)判定燃燒終止時(shí)間算法。

圖2 3種拐點(diǎn)處特殊的曲線形狀

1 燃燒時(shí)間的自動(dòng)計(jì)算方法

固體火箭發(fā)動(dòng)機(jī)燃燒時(shí)間tb主要通過對(duì)發(fā)動(dòng)機(jī)燃燒室壓強(qiáng)信號(hào)處理得到。自動(dòng)求解燃燒時(shí)間的方法有很多種，例如：曲率最大法[1-3]、燃速公式法、計(jì)算機(jī)掃描法[4]、沖量系數(shù)法[5]、壓強(qiáng)系數(shù)法[6]、線性加權(quán)極值法等。

燃速公式法是根據(jù)公式r=aPn計(jì)算燃速，其中發(fā)動(dòng)機(jī)燃速系數(shù)a和發(fā)動(dòng)機(jī)壓強(qiáng)指數(shù)n需要大量的試驗(yàn)進(jìn)行確定，不僅操作麻煩，而且當(dāng)裝藥不完全燃燒時(shí)，此方法不可行。計(jì)算機(jī)掃描法將地面試驗(yàn)中計(jì)算機(jī)采集的壓強(qiáng)-時(shí)間數(shù)據(jù)進(jìn)行模擬顯示，判斷顯示的壓強(qiáng)-時(shí)間曲線是否有明顯的拐點(diǎn)，該方法對(duì)于拐點(diǎn)不明顯的曲線不可行。沖量系數(shù)法與壓強(qiáng)系數(shù)法需要通過相同型號(hào)以往的大量數(shù)據(jù)獲得相應(yīng)的計(jì)算系數(shù)，在實(shí)際使用中比較麻煩。

圖3 尾部三段型的3種切線繪制方式

曲率最大法是一種較新的燃燒時(shí)間自動(dòng)求取方法，通過對(duì)壓強(qiáng)曲線求取一階導(dǎo)數(shù)、二階導(dǎo)數(shù)以及曲率值，最終曲率最大值所對(duì)應(yīng)的時(shí)間即為燃燒時(shí)間終點(diǎn)。這種方法可以對(duì)絕大部分發(fā)動(dòng)機(jī)壓強(qiáng)信號(hào)進(jìn)行求解，也不需要以往的數(shù)據(jù)作為基礎(chǔ)，方便快捷，適用性強(qiáng)，唯一的缺點(diǎn)是由于真實(shí)的發(fā)動(dòng)機(jī)燃燒室壓強(qiáng)-時(shí)間曲線是由離散數(shù)據(jù)組成，數(shù)據(jù)不連續(xù)也不光滑，不存在一個(gè)連續(xù)的一階導(dǎo)數(shù)。目前有利用三次樣條曲線擬合實(shí)際壓強(qiáng)-時(shí)間曲線的方法，但是由于擬合精度不夠，導(dǎo)致求得的燃燒時(shí)間與目前普遍使用的雙切線法存在一定差距。

文中在曲率最大法的基礎(chǔ)上提出了一種依托實(shí)際壓強(qiáng)-時(shí)間曲線來確定標(biāo)準(zhǔn)Ф127 mm發(fā)動(dòng)機(jī)燃燒時(shí)間的自動(dòng)算法，不需要擬合樣條曲線，直接在實(shí)際壓強(qiáng)-時(shí)間曲線上對(duì)燃燒時(shí)間進(jìn)行求解，并且為了保證以往燃速數(shù)據(jù)結(jié)果的延續(xù)性，對(duì)算法進(jìn)行了一系列改進(jìn)以使其最大程度地接近雙切線法處理結(jié)果。

2 標(biāo)準(zhǔn)Ф127 mm發(fā)動(dòng)機(jī)燃燒時(shí)間的曲率自動(dòng)算法

2.1 算法理論依據(jù)

固體發(fā)動(dòng)機(jī)內(nèi)彈道曲線的形狀取決于燃?xì)馍陕屎团懦雎实年P(guān)系。當(dāng)生成率大于排出率時(shí)，壓強(qiáng)上升；當(dāng)二者平衡時(shí)，壓強(qiáng)穩(wěn)定；當(dāng)燃?xì)馍陕市∮谌細(xì)馀懦雎蕰r(shí)，壓強(qiáng)下降。裝藥燃燒終止瞬間，燃面突然消失或突降，燃?xì)馍陕仕查g變?yōu)榱慊蛲唤担瑝簭?qiáng)曲線必然發(fā)生急劇轉(zhuǎn)折。內(nèi)彈道壓強(qiáng)-時(shí)間曲線信息中包含著終燃點(diǎn)信息，曲線變化最快的地方，即曲率最大值處可認(rèn)定為藥柱燃燒終止時(shí)刻。

2.2 算法思路

針對(duì)標(biāo)準(zhǔn)Ф127 mm固體火箭發(fā)動(dòng)機(jī)燃燒室實(shí)際采集得到的壓強(qiáng)-時(shí)間曲線進(jìn)行求解，首先確定燃燒終止時(shí)間所在的區(qū)間；然后對(duì)該區(qū)間內(nèi)的壓強(qiáng)-時(shí)間曲線進(jìn)行歸一化處理；接著給定一個(gè)合適的時(shí)間步長(zhǎng)，從所確定區(qū)間的起始位置開始逐點(diǎn)向后每隔該時(shí)間步長(zhǎng)進(jìn)行微分求導(dǎo)，如圖4所示，為從點(diǎn)1開始按步長(zhǎng)為4個(gè)采樣間隔進(jìn)行的逐點(diǎn)求導(dǎo)，并按式(1)計(jì)算該處的曲率，以曲率最大位置為燃燒終止時(shí)間點(diǎn)。

(1)

其中：K為某點(diǎn)曲率；P′(t)為某點(diǎn)壓強(qiáng)的一階導(dǎo)數(shù)；P″(t)為某點(diǎn)壓強(qiáng)的二階導(dǎo)數(shù)。

圖4 按時(shí)間步長(zhǎng)為4個(gè)采樣間隔進(jìn)行求導(dǎo)

該算法中的關(guān)鍵在于如何在給定區(qū)間內(nèi)完成壓強(qiáng)-時(shí)間曲線的“歸一化”以及確定一個(gè)合適的“求導(dǎo)時(shí)間步長(zhǎng)”。

2.2.1 壓強(qiáng)-時(shí)間曲線的“歸一化”

壓強(qiáng)-時(shí)間曲線的縱軸代表壓強(qiáng)值，如圖5(a)所示，若按照實(shí)際物理量進(jìn)行求導(dǎo)計(jì)算曲率，計(jì)算結(jié)果都處于較大的量值，不便于進(jìn)行比較，因此，對(duì)給定區(qū)間內(nèi)的縱軸按照橫軸進(jìn)行等比例變換。在歸一化之前需先選定歸一化區(qū)間，該區(qū)間必須包含燃燒終止時(shí)刻，可以肯定燃燒終止時(shí)刻一定處于最大壓強(qiáng)對(duì)應(yīng)時(shí)刻與工作結(jié)束時(shí)刻之間，因此，選定歸一化區(qū)間為[TPmax，Ta]，歸一化處理的具體算法為：

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

(2)

其中：PB為變化后的壓強(qiáng)點(diǎn)值；P為實(shí)際壓強(qiáng)點(diǎn)值；Ta為發(fā)動(dòng)機(jī)工作結(jié)束時(shí)刻點(diǎn)；TPmax為平穩(wěn)段壓強(qiáng)最大值對(duì)應(yīng)的時(shí)刻點(diǎn)；Pmax為平穩(wěn)段壓強(qiáng)最大值。

歸一化之后，縱軸的最大值和最小值與橫軸時(shí)間值一致，范圍完全相同，縱軸的區(qū)間在數(shù)值上也是[TPmax，Ta]，如圖5(b)所示。

圖5 壓強(qiáng)-時(shí)間曲線的歸一化處理

2.2.2 求導(dǎo)時(shí)間步長(zhǎng)

理想情況下，若壓強(qiáng)-時(shí)間曲線是一條光滑連續(xù)、可用具體函數(shù)進(jìn)行表征的曲線，那么根據(jù)2.1節(jié)中的理論依據(jù)，可直接通過計(jì)算曲率最大值的位置來確定燃燒終止時(shí)刻，但實(shí)際壓強(qiáng)-時(shí)間曲線是由離散點(diǎn)組成的，且采樣點(diǎn)由于電信號(hào)的干擾，在小范圍內(nèi)會(huì)有一定的波動(dòng)，既不連續(xù)也不光滑，無法直接對(duì)曲線進(jìn)行求導(dǎo)，如圖6所示。算法引入時(shí)間步長(zhǎng)求導(dǎo)的方法，利用插值法對(duì)曲線上每隔一個(gè)固定的求導(dǎo)時(shí)間步長(zhǎng)的兩個(gè)時(shí)間點(diǎn)根據(jù)式(3)進(jìn)行求導(dǎo)計(jì)算，計(jì)算結(jié)果P′(t)作為[t1，t2]中點(diǎn)時(shí)刻(t1+t2)/2的一階導(dǎo)數(shù)。使用相同的方法根據(jù)式(4)對(duì)一階導(dǎo)數(shù)P′(t)逐點(diǎn)向后進(jìn)行求導(dǎo)計(jì)算，得到二階導(dǎo)數(shù)P″(t)，之后利用式(1)計(jì)算各點(diǎn)曲率及最大值位置。

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

(3)

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

(4)

圖6 離散狀的壓強(qiáng)-時(shí)間曲線

對(duì)于一條光滑連續(xù)的曲線，根據(jù)式(5)一階導(dǎo)數(shù)的定義，求導(dǎo)時(shí)間步長(zhǎng)越短，得到的結(jié)果越準(zhǔn)確。但對(duì)于實(shí)際離散的壓強(qiáng)-時(shí)間曲線，求導(dǎo)的步長(zhǎng)太短或太長(zhǎng)獲得的結(jié)果都與雙切線法的結(jié)果相差較大，須確定一個(gè)合適的時(shí)間步長(zhǎng)來保證處理結(jié)果的準(zhǔn)確性。

(5)

為確定合適的時(shí)間步長(zhǎng)，對(duì)涵蓋高、中、低燃速近20種型號(hào)500發(fā)標(biāo)準(zhǔn)Ф127 mm發(fā)動(dòng)機(jī)試驗(yàn)數(shù)據(jù)進(jìn)行了數(shù)據(jù)分析驗(yàn)證，將其結(jié)果與雙切線法處理結(jié)果進(jìn)行對(duì)比。驗(yàn)證分為兩個(gè)階段，第一階段將時(shí)間步長(zhǎng)分別選取為0.005 s，0.01 s，0.02 s，…，0.1 s，記錄所得燃燒時(shí)間與雙切線法最接近的時(shí)間步長(zhǎng)及其所對(duì)應(yīng)的發(fā)動(dòng)機(jī)個(gè)數(shù)，并做統(tǒng)計(jì)結(jié)果如圖7所示。從圖中并未發(fā)現(xiàn)合適的時(shí)間步長(zhǎng)選取規(guī)律，但是通過進(jìn)一步分析，發(fā)現(xiàn)步長(zhǎng)與工作時(shí)間ta有一定的關(guān)系，時(shí)間步長(zhǎng)在0.02～0.04倍的工作時(shí)間范圍內(nèi)，所得結(jié)果與雙切線法結(jié)果更接近。第二階段將時(shí)間步長(zhǎng)選取為0.005ta，0.01ta，0.015ta，0.02ta，…，0.06ta，記錄所得燃燒時(shí)間與雙切線法最接近的時(shí)間步長(zhǎng)及其所對(duì)應(yīng)的發(fā)動(dòng)機(jī)個(gè)數(shù)如圖8所示。結(jié)果證明選取0.025ta作為求導(dǎo)步長(zhǎng)效果最好。因此初步選取0.025ta作為微分求導(dǎo)步長(zhǎng)進(jìn)行一階導(dǎo)數(shù)與二階導(dǎo)數(shù)求解。

圖7 最佳時(shí)間步長(zhǎng)分布規(guī)律

圖8 最佳時(shí)間步長(zhǎng)占工作時(shí)間比例分布規(guī)律

2.3 算法實(shí)例

對(duì)某標(biāo)準(zhǔn)Ф127 mm發(fā)動(dòng)機(jī)應(yīng)用此算法進(jìn)行燃速計(jì)算，使用PXI測(cè)試系統(tǒng)對(duì)該發(fā)動(dòng)機(jī)燃燒室壓強(qiáng)信號(hào)以2 000點(diǎn)/s的采樣速率進(jìn)行掃描采集，發(fā)動(dòng)機(jī)燃燒室壓強(qiáng)信號(hào)P，點(diǎn)火信號(hào)時(shí)刻點(diǎn)為t0=9.157 0 s，如圖9所示。

圖9 某標(biāo)準(zhǔn)發(fā)動(dòng)機(jī)壓強(qiáng)-時(shí)間曲線

算法流程如下：1)計(jì)算得到壓強(qiáng)最大值Pmax=13 031.17 kPa，對(duì)應(yīng)時(shí)刻點(diǎn)TPmax=9.580 5 s，工作結(jié)束時(shí)刻點(diǎn)為Ta=12.260 5 s，確定計(jì)算區(qū)間為[9.580 5 s，12.260 5 s]；

2)將[9.580 5 s，12.260 5 s]區(qū)間內(nèi)的壓強(qiáng)數(shù)值P進(jìn)行歸一化處理后為PB；

3)選取微分求導(dǎo)的時(shí)間步長(zhǎng)Δt=(12.260 5-9.580 5)×0.025=0.067 s；

4)從9.580 5 s開始，以0.067 s為步長(zhǎng)，對(duì)變換后的PB曲線逐點(diǎn)向后進(jìn)行求導(dǎo)計(jì)算，得到一階導(dǎo)數(shù)P′，對(duì)應(yīng)區(qū)間為[9.614 s，12.227 s]；

5)對(duì)一階導(dǎo)數(shù)P′在區(qū)間[10.691 0 s, 13.030 0 s]內(nèi)再次逐點(diǎn)向后進(jìn)行求導(dǎo)計(jì)算，得到二階導(dǎo)數(shù)P″，對(duì)應(yīng)區(qū)間為[9.647 5 s，12.193 5 s]；

6)利用曲率計(jì)算公式計(jì)算曲率K，得到最大曲率值為51.455 3，該點(diǎn)對(duì)應(yīng)的采集時(shí)刻值為12.046 5 s。

7)得到燃燒時(shí)間tb為12.046 5-9.157 0=2.889 5 s。

3 曲率自動(dòng)算法改進(jìn)

將自動(dòng)算法得到的燃燒時(shí)間tb與雙切線法處理結(jié)果對(duì)比，結(jié)果顯示即使選取效果最好的0.025ta作為求導(dǎo)步長(zhǎng)，兩者的平均時(shí)間差仍有0.013 9 s，結(jié)果相差較大，相較于雙切線法，曲率自動(dòng)算法計(jì)算得到的燃燒時(shí)間偏小。

為了提高曲率自動(dòng)算法處理結(jié)果與以往試驗(yàn)處理結(jié)果的延續(xù)性與一致性，需要對(duì)算法進(jìn)行改進(jìn)，使兩者的處理結(jié)果盡可能接近。通過大量的試驗(yàn)數(shù)據(jù)統(tǒng)計(jì)分析，在自動(dòng)算法中，當(dāng)采用同一時(shí)間步長(zhǎng)，若歸一化區(qū)間設(shè)定不同，得到的最大曲率位置也會(huì)有所差異，且區(qū)間越小，得到的燃燒時(shí)間越大。

目前在應(yīng)用雙切線法進(jìn)行處理時(shí)，通常只選取下降段拐點(diǎn)附近的一小段數(shù)據(jù)，對(duì)其進(jìn)行歸一化處理后繪制切線，該區(qū)間包含燃燒終止時(shí)刻，但其具體范圍是人為隨機(jī)選取的，導(dǎo)致燃燒時(shí)間求取不穩(wěn)定。曲率自動(dòng)算法的歸一化區(qū)間為[TPmax，Ta]，相較于雙切線法較大，所以計(jì)算得到的燃燒時(shí)間偏小。為了減小兩者的時(shí)間差須將曲率自動(dòng)算法的歸一化區(qū)間盡可能接近雙切線法。由于試驗(yàn)數(shù)據(jù)在沒有經(jīng)過分析之前，無法自動(dòng)選取一段包含tb且與雙切線法大致一致的分析區(qū)間，但可以肯定tb一定在最大壓強(qiáng)時(shí)刻與工作結(jié)束時(shí)刻之間，因此，首先確定數(shù)據(jù)的歸一化區(qū)間為[TPmax，Ta]，得到曲率最大位置tb1，然后將歸一化區(qū)間再次設(shè)置為[tb1-3Δt,ta]，其中Δt為求導(dǎo)時(shí)間步長(zhǎng)，得到曲率最大位置tb2，對(duì)計(jì)算結(jié)果進(jìn)行多次迭代，直到計(jì)算結(jié)果不再變化，確定其為tb。根據(jù)經(jīng)驗(yàn)，一般最多迭代3次即可。這種區(qū)間選取方法只與Δt有關(guān)，相較于雙切線法，區(qū)間選取更為確定，處理結(jié)果更加穩(wěn)定，不受人為因素干擾。利用新的區(qū)間求取方法對(duì)原曲率自動(dòng)算法進(jìn)行改進(jìn)，并對(duì)2 600余發(fā)標(biāo)準(zhǔn)Ф127 mm發(fā)動(dòng)機(jī)試驗(yàn)數(shù)據(jù)重新進(jìn)行數(shù)據(jù)分析驗(yàn)證，結(jié)果顯示改進(jìn)后的曲率自動(dòng)算法選取時(shí)間步長(zhǎng)為0.02ta效果最好，與雙切線法結(jié)果的平均時(shí)間差為0.003 6 s，標(biāo)準(zhǔn)差為0.007 5 s，結(jié)果統(tǒng)計(jì)如表1所示。

表1 自動(dòng)算法與雙切線法的燃燒時(shí)間差

4 結(jié)論

提出了一種依托實(shí)際壓強(qiáng)-時(shí)間曲線的曲率最大位置自動(dòng)算法，以此自動(dòng)確定標(biāo)準(zhǔn)Ф127 mm發(fā)動(dòng)機(jī)燃燒時(shí)間。該算法引入迭代算法，通過大量數(shù)據(jù)統(tǒng)計(jì)，當(dāng)選取發(fā)動(dòng)機(jī)工作時(shí)間的2%作為求導(dǎo)時(shí)間步長(zhǎng)，其處理結(jié)果最接近雙切線法，得到的燃燒時(shí)間差在5 ms之內(nèi)，滿足了試驗(yàn)數(shù)據(jù)結(jié)果的延續(xù)性。該算法相比于雙切線法具有可實(shí)現(xiàn)自動(dòng)化、在參數(shù)設(shè)置不變的情況下結(jié)果數(shù)據(jù)唯一的特點(diǎn)，提高了固體標(biāo)準(zhǔn)發(fā)動(dòng)機(jī)試驗(yàn)數(shù)據(jù)處理的可靠性與效率。