恒溫式量熱計測爆熱的穩(wěn)健辨識方法

楊 杰，賀元吉，趙宏偉，陳 華，韓秀鳳，占 君

(中國人民解放軍96901部隊，北京 100094)

引 言

目前，炸藥爆熱的精確測量方法主要有恒溫法與絕熱法[1-8]，兩種方法均需測量內桶水的終點溫度。其中，恒溫法外桶溫度無需跟蹤內桶溫度，但修正溫升需要計算補償；而絕熱法外桶溫度需要跟蹤內桶溫度，但修正溫升無需計算補償。

由于恒溫法避免了絕熱法外桶跟蹤內桶溫度在炸藥點火后短時間精度不足的問題，因此是測量質量較大炸藥爆熱的主要手段，但由于炸藥總放熱量大，經典恒溫法測量時間比較長。GJB772A-97方法701.1“爆熱-恒溫法和絕熱法”中規(guī)定了恒溫法初期與末期時間均為11min，而范桂榮[9]提出將初期與末期時間延長為21min、主期時間為80min時可有效減少偶然誤差。如炸藥起爆后出現系統故障(系統斷電、異常加溫、軟件崩潰等)而導致無法獲取內桶水的準確終點溫度，則該次試驗宣告失敗，將造成炸藥樣品報廢。

為降低系統故障導致測量失敗的風險，本研究提出了基于故障前測量數據利用系統辨識方法計算爆熱的解決方案。根據這一思想，首先基于量熱計傳熱方程推導分析內桶水溫度溫升曲線，然后基于系統辨識的思想給出炸藥爆熱的辨識算法，并分析系統辨識算法的收斂特性與辨識誤差，最后結合西安火炸藥一級計量站的試驗數據驗證本方法的準確性與可行性,以期為恒溫法測爆熱提供有益補充。

1 爆熱測量原理與算法

量熱計結構見圖1，包括內桶系統與外桶系統。內桶系統由爆熱彈、量熱桶組成，主要用于起爆炸藥并容納爆炸產物、吸收爆炸產物與外桶的傳熱量并測量炸藥的爆熱；外桶系統由外桶水、套夾、溫控系統等組成，用于產生一個理想的外部溫度環(huán)境。內外桶通過絕熱良好的支撐連接。

圖1 量熱計結構示意圖Fig.1 Schematic diagram of calorimeter structure

恒溫法測爆熱經典計算公式如下[1,7,9]：

(1)

ΔtJ=TI(τn)-TI(τ0)-Δθ

(2)

(3)

式中：ΔtJ、Δt為經典方法、所有方法的內桶水修正溫升，℃；Δτ為采樣時間間隔， min；TI(τ0)、TI(τn)分別為內桶主期初溫和主期末溫的數值，℃；θ0、θn分別為擬合得到初期與末期的起點溫度，℃；v0、vn分別為擬合得到初期與末期的溫升速率，℃/min；Δθ為溫度修正值，℃；n為內桶的主期溫度讀數個數；Q1為被測炸藥爆熱，J/g；C為量熱計系統熱容，J/℃；M2為傳爆藥質量，g；Q2為傳爆藥爆熱，J/g；q為雷管爆熱，J；M1為被測炸藥質量，g。

1.1 量熱計的傳熱學分析

爆熱彈金屬殼與內桶水容器外殼均為金屬結構，導熱良好，因此，可忽略金屬殼體內外表面溫差，即認為整個內桶系統溫度是一致的，則爆炸產物向內桶系統的傳熱可視為對流傳熱[10]，對應熱流量ΦN1(τ)為：

式中：QN1(τ)為炸藥點火時刻τM至τ爆炸產物向內桶系統的導熱熱量；hN為爆炸產物的對流表面?zhèn)鳠嵯禂?；AN為爆熱彈金屬殼的內表面積，kNZ=hNAN；TE(τ)、TI(τ)分別為爆炸產物、內桶水在時刻τ的測量溫度。

外桶系統向內桶系統的傳熱包括絕熱支撐的熱傳導與熱輻射等兩種形式，在內外桶溫差較小時可認為熱輻射的熱流量與內外桶溫差成正比，因此，外桶系統向內桶系統傳熱的熱流量ΦI2(τ)為：

(5)

式中：QI2(τ)為時刻τC至τ外桶系統向內桶系統的導熱熱量(為方便計算，初期階段傳熱學方程中，τC取初期起始時刻τS；主末期階段傳熱學方程中，τC取炸藥點火時刻τM)；TW為外桶水升溫結束且穩(wěn)定后的溫度；kNW為綜合熱傳導與熱輻射的導熱系數。

由于攪拌內桶水會產熱，因此，調節(jié)內外桶熱平衡時，內桶系統溫度須高于外桶系統溫度(該溫差記為TC)，以便將攪拌熱傳導至外桶，即攪拌熱QC(τ)須滿足：

(6)

初期階段是指炸藥點火前外桶水升溫結束且已穩(wěn)定的時間階段。根據能量守恒定律，內桶系統的溫升由外桶系統向內桶系統傳熱以及攪拌熱的總熱量確定，則有：

QI2(τ)+QC(τ)=cNmN(TI(τ)-Tl(τS))

(7)

式中：cN、mN分別為內桶系統的比熱容與質量。

綜合式(5)、式(6)、式(7)可得初期階段的傳熱學方程如下：

(8)

主期與末期階段分別指炸藥點火至爆炸產物放熱基本結束和之后的時間段。根據能量守恒定律，內桶系統的溫升由爆炸產物與外桶系統向內桶系統傳熱以及攪拌熱的總熱量確定，則有：

QN1(τ)+QI2(τ)+QC(τ)=
cNmN(TI(τ)-TI(τM))

(9)

式中：τM表示主期階段的起始時刻，即炸藥點火時刻，因此，τ≥τM。

這里假設炸藥點火后的產熱是瞬態(tài)的，點火后爆炸產物一直處于放熱過程，則有：

QN1(τ)=cZmZ(TE(τM)-TE(τ))

(10)

式中:cZ、mZ、TE(τM)表示爆炸產物的比熱容、質量與初始溫度，且τ≥τM。

綜合式(4)、(5)、(6)、(9)、(10)，可得主期與末期階段傳熱學方程為：

(11)

1.2 量熱計傳熱學方程解算結果

根據初期階段的傳熱學方程，初期階段內桶水溫度曲線計算結果為：

(12)

可見，初期階段內桶水溫度變化服從指數分布規(guī)律，并且收斂于TW+TC。

根據主期與末期階段的傳熱學方程，主期與末期階段內桶水溫度曲線計算結果為：

(13)

(14)

(15)

(16)

(17)

根據式(13)～式(17)，易證，λ1、λ2均為負數，且λ1>λ2；當測量時間足夠長，內桶水溫度將收斂于TW+TC。

2 爆熱計算的穩(wěn)健辨識方法

炸藥爆熱系統辨識的基本思路是：首先，根據初期階段內桶水溫度曲線辨識出參數kNW/(cNmN)，然后根據主期與末期階段內桶水溫度曲線辨識出參數u1、λ1、λ2。再根據中間參數kNW/(cNmN)、u1、λ1、λ2辨識得到TE(τM)、kNZ/cNmN、kNZ/cZmZ等3個目標參數。最后，根據目標參數辨識得到內桶水的修正溫升，再結合經典爆熱計算公式(3)得到炸藥的爆熱計算結果。

2.1 初期階段的參數辨識

由于內外桶絕熱良好，內桶系統的比熱容與質量的乘積是一個非常大的值，則：

(18)

根據泰勒展開公式，內桶水溫度曲線(見式(12))在整個初期階段具有較好的線性，即：

(19)

初期階段由于炸藥未點火，如出現系統故障，量熱計尚可調整與修復。因此，認為初期階段的數據可全部利用?；诔跗陔A段內桶水溫度采樣數據進行線性擬合得到：

TI(τ)=kS1τ+kS2+ε(τ)

(20)

式中：kS1、kS2分別為擬合得到的一次項與零次項系數；ε(τ)為擬合殘差。

比較式(19)與式(20)可得到中間參數kNW/(cNmN)的辨識結果為：

(21)

2.2 主期與末期階段的參數辨識

主期與末期階段由于炸藥已點火，如出現系統故障，量熱計調整修復難度很大，后續(xù)測量數據面臨不可信甚至丟失的風險。

針對這一情況，采用曲線擬合思想，基于起爆至測量時間τ的內桶水溫度數據辨識參數u1、λ1、λ2的目標函數為：

(22)

根據式(15)可知，kNW/(cNmN)近似為0時，λ1也十分接近0。即λ1反映的是內桶水溫的慢變過程，但慢變趨勢通常不顯著，易被快變過程與隨機擾動干擾，進而造成λ1的辨識值出現較大誤差，使得λ1的辨識值向真實值收斂速度很慢。多次仿真表明，如不對λ1進行約束，當測量時間較短時，λ1的辨識值對溫度長期的慢變過程預測能力很差。因此，為防止λ1陷入局部最優(yōu)解，加快收斂速度，求解目標函數式(22)最小值時，需對λ1進行約束。

綜合式(15)與式(18)分析可得：

(23)

考慮到kNZ?kNW、cNmN?cZmZ，以及kNW/(cNmN)本身的辨識誤差，結合式(23)，本研究給中間參數λ1的約束條件為：

(24)

利用最優(yōu)化方法求解有約束條件式(24)的目標函數f(u1,λ1,λ2)的最小值，即可辨識得到測量時間τ對應的中間參數u1、λ1、λ2。然后，根據式(14)～式(17)可得TE(τM)、kNZ/cNmN、kNZ/cZmZ等3個目標參數在測量時間τ的辨識方法見式(25)：

(25)

式中:(λ1)τ表示參數λ1在時刻τ的系統辨識值;依此類推，其他參數若由小括號及τ的下標構成，則表示該參數在τ時刻的系統辨識值。

2.3 炸藥爆熱的計算

解算主期與末期階段的傳熱學方程(11)得到爆炸產物向內桶系統傳熱量時間函數QN1(τ)，可證明，測量時間無限長時，有：

QN1(+∞)=cZmZ(TE(τM)-TC-TW)

(26)

根據修正溫升的定義可知，內桶水修正溫升在測量時刻τ的辨識方法為：

(27)

因此，代入式(25)目標參數TE(τM)、kNZ/cNmN、kNZ/cZmZ在時刻τ的辨識值，即可得到內桶水在測量時間τ的修正溫升，再根據公式(3)即可計算得到炸藥爆熱。

2.4 辨識誤差分析

2.4.1 收斂精度分析

(28)

比較式(28)與式(27)可知，穩(wěn)健系統辨識方法與經典方法得到的內桶水的修正溫升理論值相同，因此爆熱理論值也相同。

2.4.2 收斂穩(wěn)健性分析

中間參數u1、λ2辨識值較大，收斂比較穩(wěn)定；λ1辨識值近似為0，雖然通過約束式(24)提高了對溫度慢變項的長期預測精度，但溫度隨機擾動與快變項的干擾易引起λ1辨識值振蕩，見圖2。圖中τN表示以炸藥起爆為計時零點的測量時間。因此，本節(jié)主要分析λ1辨識結果對各參數辨識誤差的影響。

圖2 某次試驗u1、λ1、λ2辨識值變化情況Fig.2 Change of identified value of u1, λ1, λ2 in an experiment

將λ1在測量時刻τ的辨識值表示為：

(29)

根據式(24)可知，δτ的值域為[0,0.1]。將式(29)代入式(25)可得：

(30)

模擬結果表明，δτ大部分辨識時間計算結果在10-6至10-3范圍內，但收斂較差，即使爆熱辨識值已收斂至3%誤差范圍內，δτ也常常出現劇烈振蕩，最大值與最小值相差也通?？蛇_數百倍，見圖3。圖4為目標參數與修正溫升變化情況。

圖3 某次試驗δτ的變化情況Fig.3 Change of δτ in an experiment

圖4 某次試驗目標參數與修正溫升變化情況Fig.4 Change of target parameters and modified temperature raise in an experiment

根據這些特點，結合式(23)、式(29)、式(30)，可得到如下結論：

(1)目標參數TE(τM)、kNZ/cNmN對λ1辨識值微小波動異常敏感，見圖4。δτ的劇烈振蕩致其難以收斂，其辨識結果也不可信。

(2)由于δτ很小，內桶水修正溫升ΔtS與kNZ/cZmZ的辨識值具有較快的收斂速度、較高的辨識精度以及較好的可信度，見圖4。再比較式(28)與式(27)可知，爆熱辨識值理論上收斂于經典值，因此，本算法是穩(wěn)健的。

(3)中間參數λ1的辨識值近似收斂于kNW/(cNmN)，因此，增加約束條件(24)雖會出現λ1辨識值小幅振蕩的現象，卻有利于爆熱計算值的快速收斂。

(4)從減少辨識誤差的角度講，增大外桶對內桶傳熱的影響，即增大kNW/(cNmN)有利于提高λ1的辨識精度，并減少溫度擾動的影響，但外桶溫度的波動將會影響到內桶溫度，這顯然是得不償失的。

(5)從故障診斷角度講，如果辨識得到的參數kNW/(cNmN)或λ1較大，則說明外桶至內桶的絕熱性能下降或故障，因此，本方法也具有對量熱計進行故障診斷的能力。

2.5 模型誤差分析

由于非理想炸藥爆轟后存在金屬粉氧化的二次反應，以及在含氧氣環(huán)境下的燃燒反應，而氧化/燃燒時間要遠大于爆轟時間，與本研究模型的瞬態(tài)假設存在一定的偏離。因此需要分析非理想炸藥爆轟后燃燒/氧化產熱對辨識誤差的影響?？紤]燃燒/氧化產熱影響后，式(10)變?yōu)椋?/p>

(31)

如果可燃物的消耗速率是先快后慢，則更有利于模型關于反應過程是瞬態(tài)的假設，為了最大化地分析模型誤差，這里給出不利于瞬態(tài)性的假設，即：可燃物的消耗是勻速的，因此，有：

(32)

為了與瞬態(tài)模型的總產熱量對應，得到總的燃燒/氧化產熱量為：

(33)

由于燃燒/氧化時間較短，而攪拌產熱與外桶傳熱速率較慢，因此，結合式(4)、式(31)、式(32)、式(9)，并忽略攪拌產熱與外桶傳熱的影響，可得燃燒/氧化階段爆炸產物向內桶系統的傳熱方程為：

(34)

結合式(33)，解算式(34)可知，如果

(35)

則燃燒/氧化階段結束時，爆炸產物與內桶水的溫度近似解為：

(36)

根據文獻[11]中TNT炸藥在純氧環(huán)境爆轟試驗溫度曲線可知，炸藥燃燒/氧化時間不超過1/30min，結合前文kNZ/cZmZ、kNZ/cNmN的辨識情況來看，式(35)是成立的。

燃燒/氧化階段結束后的傳熱模型與瞬態(tài)過程的傳熱模型是一致的，只是初始狀態(tài)變?yōu)槭?36)，因此，解算得到燃燒/氧化階段結束后內桶水的溫升曲線為：

(37)

式中：λ1、λ2、Δ的取值與式(15)、式(16)、式(17)相同。

(38)

為了與瞬態(tài)模型中內桶水溫升表達式對應，將式(37)改寫為：

(39)

對比瞬態(tài)模型式(13)與燃燒/氧化模型式(37)，再結合cZmZ?cNmN，kNW?kNZ，可得式(39)中的參數為：

(40)

式中：ΔTA表示試驗預設溫升；ΔTF表示燃燒/氧化引起的修正溫升，且有：

(41)

式(40)表明勻速燃燒/氧化模型相對于瞬態(tài)模型的參數偏移量。為了說明參數偏移量的上限，這里假設炸藥產熱過程全部為勻速燃燒/氧化過程，燃燒/氧化時間τB=1/30min，根據表2中8個樣本的辨識情況計算得到u1的最大偏移量為1.54×10-4℃，TI(τM)的最大偏移量為1.59×10-2℃?？梢?，即使在最差的假設條件下，模型誤差造成的參數偏移量也是可以忽略不計的，因此，本算法是足夠穩(wěn)健的，可以滿足無氧/有氧環(huán)境下的爆熱測量精度要求。

3 試驗仿真分析

對爆熱值為4000～9000J/g的8個典型樣本進行了恒溫法爆熱測量，結果見表2。試驗步驟參照國軍標GJB772A-97中方法701.1爆熱-恒溫法，為減少隨機誤差，初期、主期、末期分別延長為21、100、20min，其他時間主要為溫度調整時間，采樣時間間隔為1min[9]。

定義測量時刻τ爆熱辨識值相對誤差為：

(42)

式中:Q1為目標爆熱值(這里取測量結束后經典法得到的爆熱值)，J/g；(Q1)τ為炸藥爆熱在測量時刻τ的系統辨識值，J/g。

定義相對誤差上限水平(記為UN)為起爆后測量時間τN之后的最大辨識相對誤差；定義極限收斂水平(記為LN)為相對誤差上限水平的最小值。

針對其中某次試驗進行具體分析，內外桶溫度及擬合溫度隨采樣時間變化情況見圖5。由圖5可見，炸藥起爆前的初期階段與爆炸產物放熱結束后的末期階段，內桶水溫度曲線具有極好的線性；系統辨識法得到的整個主末期內桶水溫度曲線具有極好的擬合度。

圖5 內外桶實際溫度及各階段的擬合曲線Fig.5 Actual temperature of inner and outer barrels and fitting curve of each stage

為了檢驗系統辨識方法對內桶水長期溫升的預測能力，取起爆后第30、60、90min時刻辨識得到的溫度曲線參數，分別代入內桶水溫升公式(13)，得到擬合溫升曲線與實際溫升曲線，結果見圖6。

從圖6可知，當測量時間達到起爆后第60min時就可以相當精確地預測內桶水溫的變化趨勢。

圖6 不同測量點辨識得到內桶水的理論與實際溫升曲線的比較Fig.6 Comparison of the theoretical and practical temperature rise curves of the inner barrel water obtained from the identif ication of different measuring points

爆熱系統辨識仿真從起爆后第5min開始，取測量結束后用經典法計算得到的爆熱值4874.52J/g作為目標值。仿真結果表明，爆熱辨識值的相對誤差隨起爆后測量時間的增加而迅速減小，在起爆后第29min、31min分別收斂到了3%、2%的相對誤差水平(對應爆熱辨識值4746J/g、4781J/g)，在第32min達到1.72%的相對誤差水平后，相對誤差在此誤差水平內緩慢變化，并于第119min達到了1.66%的極限收斂水平，分別見表1、圖7和圖8。其中，圖中時間從辨識誤差第一次收斂至3%開始顯示。

圖7 爆熱辨識值的收斂情況Fig.7 Convergence rule to identifying value of explosion heat

圖8 爆熱辨識值相對誤差上限水平的變化情況Fig.8 Change of UN of identified value of explosion heat

表1 典型測量時刻爆熱辨識值及誤差情況Table 1 Identified value and error of explosion heat in typical measuring time

為了驗證本算法的普適性，試驗選用了由HMX、Al粉、黏結劑等組成的E-1、E-2樣本，由CL-20、AP、Al粉、黏結劑等組成的E-3、E-4樣本，由TATB、Al粉、黏結劑等組成的E-5～E-8樣本，針對8個典型樣本利用穩(wěn)健系統辨識方法得到爆熱辨識值的收斂特性見表2。

表2 爆熱辨識值的收斂特性Table 2 Convergence characteristics of identification value of explosion heat

由表2可見，爆熱辨識值收斂至3%的相對誤差水平測量時間不超過49min，收斂至2%的相對誤差水平測量時間也不超過6min(約占主末期的1/2)，其極限收斂水平均未超過2%，且大部分樣本極限收斂水平未超過0.5%，充分說明了本算法的穩(wěn)健性、精確性和快速性。

多次仿真表明，本方法得到的爆熱辨識值相對誤差近似按指數規(guī)律快速收斂至不超過3%的相對誤差水平，然后在該范圍內緩慢變化，收斂時間一般不超過主末期階段的1/2，其辨識結果具有相當高的精度與可信度。而且，本算法辨識得到的中間參數kNW/(cNmN)可作為評判量熱計的絕熱性、一致性與穩(wěn)定性等傳熱性能的重要依據。

為了在實際測量中判斷爆熱辨識值的收斂時刻，本研究基于修正溫升辨識值時間曲線在進入收斂階段后具有波動小、變化慢的特點，在總結表2中8個樣本的修正溫升辨識值時間曲線數據的基礎上，提出可供參考的爆熱辨識值收斂時刻的試驗判據法如下：對修正溫升辨識時間曲線的時刻τJ及其前15min范圍數據進行線性擬合，如果擬合得到的直線斜率絕對值不超過5.5×10-3℃/min，擬合殘差均方根不超過6.0×10-3℃，則時刻τJ對應爆熱辨識值相對誤差上限水平不超過3%。

根據該試驗判據，對表2中8個樣本的修正溫升辨識值曲線進行搜索，得到滿足3%誤差上限水平的第一個測量時刻數據見表3。

表3 3%相對誤差上限水平時間點的搜索結果Table 3 Search result to the time of 3% upper limit level of relative error

從表3可知，本研究提出的試驗判據可有效指導實際測量中爆熱辨識值收斂時間的判斷,但由于該試驗判據是基于小樣本試驗的統計結果，所以還需要在實踐中完善。

4 結 論

(1)穩(wěn)健辨識方法得到的爆熱辨識值能在不超過主末期階段1/2時間內快速、可靠性收斂至3%相對誤差上限水平內，顯著降低了因試驗故障導致測量失敗的風險。

(2)誤差分析表明，爆熱辨識值與經典法計算值的收斂值相等，即系統辨識方法與經典方法理論上具有相同的測量精度；雖然量熱計測爆熱的傳熱學模型基于“炸藥產熱是瞬態(tài)過程”的假設建立，但氧化、燃燒等非瞬態(tài)產熱過程對模型造成的偏移是可以忽略的，因此，本方法是足夠穩(wěn)健的。

(3)本研究提出的收斂時間試驗判據能可靠地判斷爆熱辨識值的收斂時間，但由于該試驗判據是基于小樣本試驗的統計值，所以還需要在實踐中完善。本方法還可應用于爆熱預估、內桶水溫升趨勢預測、故障診斷等領域，是恒溫法爆熱測量經典方法的有益補充。