恒温式量热计测爆热的稳健辨识方法

2020-09-01贺元吉赵宏伟韩秀凤

火炸药学报 2020年4期

杨杰，贺元吉，赵宏伟，陈华，韩秀凤，占君

(中国人民解放军96901部队，北京 100094)

引言

目前，炸药爆热的精确测量方法主要有恒温法与绝热法[1-8]，两种方法均需测量内桶水的终点温度。其中，恒温法外桶温度无需跟踪内桶温度，但修正温升需要计算补偿；而绝热法外桶温度需要跟踪内桶温度，但修正温升无需计算补偿。

由于恒温法避免了绝热法外桶跟踪内桶温度在炸药点火后短时间精度不足的问题，因此是测量质量较大炸药爆热的主要手段，但由于炸药总放热量大，经典恒温法测量时间比较长。GJB772A-97方法701.1“爆热-恒温法和绝热法”中规定了恒温法初期与末期时间均为11min，而范桂荣[9]提出将初期与末期时间延长为21min、主期时间为80min时可有效减少偶然误差。如炸药起爆后出现系统故障(系统断电、异常加温、软件崩溃等)而导致无法获取内桶水的准确终点温度，则该次试验宣告失败，将造成炸药样品报废。

为降低系统故障导致测量失败的风险，本研究提出了基于故障前测量数据利用系统辨识方法计算爆热的解决方案。根据这一思想，首先基于量热计传热方程推导分析内桶水温度温升曲线，然后基于系统辨识的思想给出炸药爆热的辨识算法，并分析系统辨识算法的收敛特性与辨识误差，最后结合西安火炸药一级计量站的试验数据验证本方法的准确性与可行性,以期为恒温法测爆热提供有益补充。

1 爆热测量原理与算法

量热计结构见图1，包括内桶系统与外桶系统。内桶系统由爆热弹、量热桶组成，主要用于起爆炸药并容纳爆炸产物、吸收爆炸产物与外桶的传热量并测量炸药的爆热；外桶系统由外桶水、套夹、温控系统等组成，用于产生一个理想的外部温度环境。内外桶通过绝热良好的支撑连接。

图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为爆炸产物的对流表面传热系数；AN为爆热弹金属壳的内表面积，kNZ=hNAN；TE(τ)、TI(τ)分别为爆炸产物、内桶水在时刻τ的测量温度。

外桶系统向内桶系统的传热包括绝热支撑的热传导与热辐射等两种形式，在内外桶温差较小时可认为热辐射的热流量与内外桶温差成正比，因此，外桶系统向内桶系统传热的热流量ΦI2(τ)为：

(5)

式中：QI2(τ)为时刻τC至τ外桶系统向内桶系统的导热热量(为方便计算，初期阶段传热学方程中，τC取初期起始时刻τS；主末期阶段传热学方程中，τC取炸药点火时刻τM)；TW为外桶水升温结束且稳定后的温度；kNW为综合热传导与热辐射的导热系数。

由于搅拌内桶水会产热，因此，调节内外桶热平衡时，内桶系统温度须高于外桶系统温度(该温差记为TC)，以便将搅拌热传导至外桶，即搅拌热QC(τ)须满足：

(6)

初期阶段是指炸药点火前外桶水升温结束且已稳定的时间阶段。根据能量守恒定律，内桶系统的温升由外桶系统向内桶系统传热以及搅拌热的总热量确定，则有：

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

(7)

式中：cN、mN分别为内桶系统的比热容与质量。

综合式(5)、式(6)、式(7)可得初期阶段的传热学方程如下：

(8)

主期与末期阶段分别指炸药点火至爆炸产物放热基本结束和之后的时间段。根据能量守恒定律，内桶系统的温升由爆炸产物与外桶系统向内桶系统传热以及搅拌热的总热量确定，则有：

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

(9)

式中：τM表示主期阶段的起始时刻，即炸药点火时刻，因此，τ≥τM。

这里假设炸药点火后的产热是瞬态的，点火后爆炸产物一直处于放热过程，则有：

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

(10)

式中:cZ、mZ、TE(τM)表示爆炸产物的比热容、质量与初始温度，且τ≥τM。

综合式(4)、(5)、(6)、(9)、(10)，可得主期与末期阶段传热学方程为：

(11)

1.2 量热计传热学方程解算结果

根据初期阶段的传热学方程，初期阶段内桶水温度曲线计算结果为：

(12)

可见，初期阶段内桶水温度变化服从指数分布规律，并且收敛于TW+TC。

根据主期与末期阶段的传热学方程，主期与末期阶段内桶水温度曲线计算结果为：

(13)

(14)

(15)

(16)

(17)

根据式(13)～式(17)，易证，λ1、λ2均为负数，且λ1>λ2；当测量时间足够长，内桶水温度将收敛于TW+TC。

2 爆热计算的稳健辨识方法

炸药爆热系统辨识的基本思路是：首先，根据初期阶段内桶水温度曲线辨识出参数kNW/(cNmN)，然后根据主期与末期阶段内桶水温度曲线辨识出参数u1、λ1、λ2。再根据中间参数kNW/(cNmN)、u1、λ1、λ2辨识得到TE(τM)、kNZ/cNmN、kNZ/cZmZ等3个目标参数。最后，根据目标参数辨识得到内桶水的修正温升，再结合经典爆热计算公式(3)得到炸药的爆热计算结果。

2.1 初期阶段的参数辨识

由于内外桶绝热良好，内桶系统的比热容与质量的乘积是一个非常大的值，则：

(18)

根据泰勒展开公式，内桶水温度曲线(见式(12))在整个初期阶段具有较好的线性，即：

(19)

初期阶段由于炸药未点火，如出现系统故障，量热计尚可调整与修复。因此，认为初期阶段的数据可全部利用。基于初期阶段内桶水温度采样数据进行线性拟合得到：

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

(20)

式中：kS1、kS2分别为拟合得到的一次项与零次项系数；ε(τ)为拟合残差。

比较式(19)与式(20)可得到中间参数kNW/(cNmN)的辨识结果为：

(21)

2.2 主期与末期阶段的参数辨识

主期与末期阶段由于炸药已点火，如出现系统故障，量热计调整修复难度很大，后续测量数据面临不可信甚至丢失的风险。

针对这一情况，采用曲线拟合思想，基于起爆至测量时间τ的内桶水温度数据辨识参数u1、λ1、λ2的目标函数为：

(22)

根据式(15)可知，kNW/(cNmN)近似为0时，λ1也十分接近0。即λ1反映的是内桶水温的慢变过程，但慢变趋势通常不显著，易被快变过程与随机扰动干扰，进而造成λ1的辨识值出现较大误差，使得λ1的辨识值向真实值收敛速度很慢。多次仿真表明，如不对λ1进行约束，当测量时间较短时，λ1的辨识值对温度长期的慢变过程预测能力很差。因此，为防止λ1陷入局部最优解，加快收敛速度，求解目标函数式(22)最小值时，需对λ1进行约束。

综合式(15)与式(18)分析可得：

(23)

考虑到kNZ≫kNW、cNmN≫cZmZ，以及kNW/(cNmN)本身的辨识误差，结合式(23)，本研究给中间参数λ1的约束条件为：

(24)

利用最优化方法求解有约束条件式(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)可知，稳健系统辨识方法与经典方法得到的内桶水的修正温升理论值相同，因此爆热理论值也相同。

2.4.2 收敛稳健性分析

中间参数u1、λ2辨识值较大，收敛比较稳定；λ1辨识值近似为0，虽然通过约束式(24)提高了对温度慢变项的长期预测精度，但温度随机扰动与快变项的干扰易引起λ1辨识值振荡，见图2。图中τN表示以炸药起爆为计时零点的测量时间。因此，本节主要分析λ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)可知，爆热辨识值理论上收敛于经典值，因此，本算法是稳健的。

(3)中间参数λ1的辨识值近似收敛于kNW/(cNmN)，因此，增加约束条件(24)虽会出现λ1辨识值小幅振荡的现象，却有利于爆热计算值的快速收敛。

(4)从减少辨识误差的角度讲，增大外桶对内桶传热的影响，即增大kNW/(cNmN)有利于提高λ1的辨识精度，并减少温度扰动的影响，但外桶温度的波动将会影响到内桶温度，这显然是得不偿失的。

(5)从故障诊断角度讲，如果辨识得到的参数kNW/(cNmN)或λ1较大，则说明外桶至内桶的绝热性能下降或故障，因此，本方法也具有对量热计进行故障诊断的能力。

2.5 模型误差分析

由于非理想炸药爆轰后存在金属粉氧化的二次反应，以及在含氧气环境下的燃烧反应，而氧化/燃烧时间要远大于爆轰时间，与本研究模型的瞬态假设存在一定的偏离。因此需要分析非理想炸药爆轰后燃烧/氧化产热对辨识误差的影响。考虑燃烧/氧化产热影响后，式(10)变为：

(31)

如果可燃物的消耗速率是先快后慢，则更有利于模型关于反应过程是瞬态的假设，为了最大化地分析模型误差，这里给出不利于瞬态性的假设，即：可燃物的消耗是匀速的，因此，有：

(32)

为了与瞬态模型的总产热量对应，得到总的燃烧/氧化产热量为：

(33)

由于燃烧/氧化时间较短，而搅拌产热与外桶传热速率较慢，因此，结合式(4)、式(31)、式(32)、式(9)，并忽略搅拌产热与外桶传热的影响，可得燃烧/氧化阶段爆炸产物向内桶系统的传热方程为：

(34)

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

(35)

则燃烧/氧化阶段结束时，爆炸产物与内桶水的温度近似解为：

(36)

根据文献[11]中TNT炸药在纯氧环境爆轰试验温度曲线可知，炸药燃烧/氧化时间不超过1/30min，结合前文kNZ/cZmZ、kNZ/cNmN的辨识情况来看，式(35)是成立的。

燃烧/氧化阶段结束后的传热模型与瞬态过程的传热模型是一致的，只是初始状态变为式(36)，因此，解算得到燃烧/氧化阶段结束后内桶水的温升曲线为：

(37)

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

(38)

为了与瞬态模型中内桶水温升表达式对应，将式(37)改写为：

(39)

对比瞬态模型式(13)与燃烧/氧化模型式(37)，再结合cZmZ≪cNmN，kNW≪kNZ，可得式(39)中的参数为：

(40)

式中：ΔTA表示试验预设温升；ΔTF表示燃烧/氧化引起的修正温升，且有：

(41)

式(40)表明匀速燃烧/氧化模型相对于瞬态模型的参数偏移量。为了说明参数偏移量的上限，这里假设炸药产热过程全部为匀速燃烧/氧化过程，燃烧/氧化时间τB=1/30min，根据表2中8个样本的辨识情况计算得到u1的最大偏移量为1.54×10-4℃，TI(τM)的最大偏移量为1.59×10-2℃。可见，即使在最差的假设条件下，模型误差造成的参数偏移量也是可以忽略不计的，因此，本算法是足够稳健的，可以满足无氧/有氧环境下的爆热测量精度要求。

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个典型样本利用稳健系统辨识方法得到爆热辨识值的收敛特性见表2。

表2 爆热辨识值的收敛特性Table 2 Convergence characteristics of identification value of explosion heat

由表2可见，爆热辨识值收敛至3%的相对误差水平测量时间不超过49min，收敛至2%的相对误差水平测量时间也不超过6min(约占主末期的1/2)，其极限收敛水平均未超过2%，且大部分样本极限收敛水平未超过0.5%，充分说明了本算法的稳健性、精确性和快速性。

多次仿真表明，本方法得到的爆热辨识值相对误差近似按指数规律快速收敛至不超过3%的相对误差水平，然后在该范围内缓慢变化，收敛时间一般不超过主末期阶段的1/2，其辨识结果具有相当高的精度与可信度。而且，本算法辨识得到的中间参数kNW/(cNmN)可作为评判量热计的绝热性、一致性与稳定性等传热性能的重要依据。

为了在实际测量中判断爆热辨识值的收敛时刻，本研究基于修正温升辨识值时间曲线在进入收敛阶段后具有波动小、变化慢的特点，在总结表2中8个样本的修正温升辨识值时间曲线数据的基础上，提出可供参考的爆热辨识值收敛时刻的试验判据法如下：对修正温升辨识时间曲线的时刻τJ及其前15min范围数据进行线性拟合，如果拟合得到的直线斜率绝对值不超过5.5×10-3℃/min，拟合残差均方根不超过6.0×10-3℃，则时刻τJ对应爆热辨识值相对误差上限水平不超过3%。

根据该试验判据，对表2中8个样本的修正温升辨识值曲线进行搜索，得到满足3%误差上限水平的第一个测量时刻数据见表3。