预混气体爆轰胞格结构定性仿真建模*
2013-04-24曲忠伟颜事龙马宏昊
曲忠伟，颜事龙，马宏昊
(1.安徽理工大学理学院，安徽 淮南 232001；2.中国科学技术大学近代力学系，安徽 合肥 230027)
胞格结构能够有效反映可燃气管道内爆轰、爆燃的过程，烟迹捕捉技术是研究胞格结构的重要手段。利用胞格结构分析手段， Thomas等[1]研究了三波点迹线随楔角角度增大而倾斜的实验现象。Ohyagi等[2]研究了在管道中放置后向梯形块时气体爆轰波的绕射现象。Bartlma等[3]研究了爆轰波演变与入射爆轰波横波间距的关系。Edwards等[4]利用烟迹技术分析了气相爆轰波在90°弯管的管道中的绕射过程。国内张德良[5]实验研究了爆轰波通过管道狭缝时胞格的演变过程。黄中伟等[6]研究了胞格尺寸与超压的关系。夏昌敬、周凯元等[7-9]研究了丙烷/氧/空气混合气非稳定爆轰波传播特性和通过90°弯管时的传播特性，研究了90°矩形弯管中氢/氧/氩混合气爆轰波的胞格结构。
胞格尺寸作为爆轰波的一种基本特征尺寸，与气体的可爆性、最小起爆能量、爆轰传播的临界直径有直接联系[10-11]。Moen等[12]研究了氢-氧-稀释剂混合气体的胞格状况，讨论了初始压力和稀释剂对胞格尺寸的影响，同时给出了预混胞格尺寸的经验公式。Bull等[13]在矩形截面激波管里研究了氢气、乙烯、乙炔、乙烷与空气混合物的胞格结构，结果发现：随着初压和温度的升高，胞格尺寸逐渐减小。
初始条件对可燃气体非稳定爆轰波胞格结构影响很大，初始压力、密度、温度、惰性气体混合等都会影响爆轰波的胞格尺寸，下面主要是通过定量的实验结合定性的分析初始条件对爆轰胞格结构的影响，验证定性仿真在预测爆轰胞格结构试验中的准确性和可行性。
1 实验装置设计
使用高压放电系统在封闭端点火，使用压力传感器分别测量火焰传播速度、冲击波压力。传感器的布置范围应涵盖火焰传播初始阶段、火焰缓慢传播阶段、火焰传播加速阶段、爆燃转爆轰(DDT)阶段。使用高速纹影系统拍摄丙烷/空气爆燃、爆轰火焰的演化过程。利用烟迹技术捕捉爆燃波、爆轰波产生、传播的全过程，分析其胞格结构演化过程，实验设计见图1。
图1 直管道系统实验示意图(以3组传感器示意)A.点火电极;B.真空阀;C.水银压力计;D.真空泵;E.混气罐;F.压力传感器;G.波形存贮器;H.电荷放大器;J.熏烟钢片Fig.1 Straight piping systems experimental diagram(3 sensor schematic)A.ignition electrode;B.vacuum valve;C.mercury manometer;D. vacuum pump; E. mixed gas cylinders; F. pressure sensor; G. waveform memory; H. charge amplifier; J. fume steel
利用本实验室现有的管道设备(50 mm), 设计具有一个、两个弯道的管道系统开展实验，见图2。使用高压放电系统在封闭端点火，使用光电传感器、压力传感器分别测量火焰传播速度、冲击波压力。传感器的布置范围应涵盖火焰传播初始阶段、火焰缓慢传播阶段、火焰传播加速阶段、爆燃转爆轰(DDT)阶段。其中，爆燃火焰的加速过程以及爆轰火焰在通过弯管时的变化规律是研究的重点。使用高速纹影系统(图3)拍摄丙烷/空气爆轰火焰通过弯管时的特征图像。利用烟迹技术捕捉爆燃波加速过程以及爆轰波在弯管处的变化规律，分析其胞格结构演化过程。
图2 弯管传播规律实验示意图A. 点火电极; B.真空阀; C.水银压力计; D.真空泵; E.混气罐; F.压力传感器; H.电荷放大器; G.波形存储器; J.熏烟钢片; N.光电传感器Fig.2 The laws of propagation experiments schematic diagram in bend pipelineA.ignition electrode;B.vacuum valve;C.mercury manometer;D. vacuum pump; E. mixed gas cylinders; F. pressure sensor;H. charge amplifier;G. waveform memory; J. fume steel; N. photoelectric sensor
图3 高速纹影系统示意图A.He-Ne激光器 B.光强调节器 C.光束调节器 D.扩束镜 E.反射镜 F.反射凹面镜 G.观察窗H.观察区域 J.刀口 K.滤光片 L.高速摄像机Fig.3 High-speed schlieren system diagramA. He-Ne laser light B. stressed section C. the beam adjuster D. beam expander mirror E. mirrorF. reflective concave mirror G. observation window H.observed region J.incision K.lighter filter L.speed cameras
实验手段：
1)实验前，制作熏烟钢片，调试设备，对管道抽真空，送入预先混合配置好的丙烷/空气混合气体；
2)开启各种探测记录设备和高速纹影系统，采用高压放电系统对丙烷/空气混合气体点火，产生爆燃、爆轰过程。
2 仿真建模
2.1 QSIM算法
QSIM算法中研究的函数都属于可推理的函数[14-15]。下面介绍可推理函数的定义：
定义1：对[a,b]⊆R*，f：[a,b]→R*是可推理函数的充分必要条件是：
1)f在 [a，b] 区间上连续；
2)f在(a，b)区间上连续可微；
3)f在任意有界区间上，都会有有限的一阶导数等于0的临界点，即f含有有限奇数点；
4)在R*上存在着单边的极限值limt->a+f′(t)和limt->b-f′(t)，然后用这两个极限值去分别定义边界值f′(a)和f′(b)。
定义2：在可推理函数行为、状态上的重要点集合被称为路标值。路标值集合中都包含有[a，b] 区间边界点上的f(a)、f(b)和函数f等于0的点；此外，在定性的推理过程当中，可根据推理过程的需要向路标值集合当中增加新的路标值，路标值所在的时间点也称作可区分时间点。
定义3：路标值的集合，即可推理函数随着系统仿真时间而变化的定性值，是l1
(1)
定义4：一个系统可以由可推理函数F=f1,f2,…,fm组成，它们具有各自的路标值集合、可区分时间点。F的可区分时间点集合，是f的可区分时间点上并集，系统的定性状态是根据单个的函数定性状态而组成的一个m元组：
QS(S,ti)=[QS(f1,ti),…,QS(fm,ti)]
QS(S,ti,ti+1)=[QS(f1,ti,ti+1),…,QS(fm,ti,ti+1)]
(2)
定义5：可推理函数F定性状态序列能够构成一个系统在时间区间[a，b]上的定性行为，QS(F,t0)，QS(F,t0,t1)，QS(F,t1)，QS(F,t1,t2)，…，QS(F,tn)。
2.2 QSIM算法的约束条件(Q2)
在定性仿真过程当中，通常描述一个系统的结构为参数之间的约束集合，并且根据约束来对系统的行为进行推理[14]，从这里不难发现仿真的主要依据是约束。在QSIM算法当中也是通过约束(代数约束和定性约束)来表示系统变量间的关系。
常见的代数约束有加约束、取反约束、微分约束、乘约束。
在函数关系当中，最常见也是最重要的关系是单调的函数关系，下面介绍的M+表示的就是函数之间单调增关系，而M-表示的就是函数之间单调减关系，它们的定义分别如下：
定义6：M+关系为真的充分必要条件是，存在函数H(t)，它的定义域为g，值域是f，有H′(t)>0；即f和g成正比关系并且单调性一致。
定义7：M-关系为真的充分必要条件是，存在函数H(t)，它的定义域为g，值域是g，有H′(t)<0；即f和g成反比关系并且单调性相反。
2.3 基于QSIM的预混气体爆轰波胞格结构建模
考虑到胞格尺寸依赖于ZND诱导区的程度，及λ =AΔI，其中A是常数但是却随着初始条件的变化而变化，ΔI为诱导区长度。白方舟等[15]考虑到A值是依赖于初始条件下的参数，因此得到了特征参数A是关于预混气体的当量比(φ)，氩气稀释浓度(XAr),初始压力(p1/p0)的函数，如下所示：
A=Cφα(1-XAr)β(p1/p0)γ
其中，C为常数，α、β、γ分别为指数。
在初始条件确定的情况下，爆轰胞格与其诱导区长度之间的比例系数A是常数，诱导区长度ΔI可以定量计算的，这样爆轰胞格尺寸就完全由初始状态决定了。在这里定义爆轰胞格特征曲线
λ=f(φ,XAr,p)ΔI
其中φ表示的是预混气体的当量比,XAr是氩气稀释浓度,p表示的是初始压力。本文主要研究预混气体爆轰胞格尺寸(λ)与混合气体当量比(φ)、氩气稀释浓度(XAr)和初始压力(P)等因素之间的关系。
下面采用基于定性约束和量空间的知识对爆轰胞格结构进行定性的描述和建模。
在某一时刻ti爆轰波的定性状态QS(S,ti)可以根据若干个单个的函数定性状态而组成的一个m元组，如下表示：
QS(S,ti)=[QS(f1,ti),…,QS(fm,ti)]
QS(S,ti,ti+1)=[QS(f1,ti,ti+1),…,QS(fm,ti,ti+1)]
此外，可以通过定性建模的方法表示在时刻t，胞格尺寸的定性度量和定性变化方向：
qval(λ,t) =
上式中qmag称为定性度量，而qdir表示的是定性变化方向。
综上所述，我们可以假定在某个时刻爆轰胞格尺寸(λ)与混合气体当量比(φ)、氩气稀释浓度(XAr)和初始压力(P)定性约束关系如下表示：
QS(λ,ti)=[QS(φ,ti),QS(XAr,ti),QS(P,ti)]
显而易见，上式中任何时刻的爆轰胞格结构与此刻的混合气体当量比(φ)、氩气稀释浓度(XAr)和初始压力(P)都有着密切的关系。在以前的文章中都会有学者对某些影响因素开展过研究，也得到了一些结论，本文采取的定性建模仿真方法综合分析初始条件的影响作用，下面进行一系列相关的实验对建模进行分析和验证。
3 爆轰胞格尺寸与初始条件的关系
在经过大量的实验后取得以下的数据，下图中的每一个实验室点都是经过多次实验室数据的算术平均值得到的，实验当日为1个大气压，温度23 ℃，湿度58%。
3.1 爆轰胞格尺寸与混合气体当量比的关系
可燃性气体在不同的当量比时的胞格尺寸与当量形成“U”形的关系。一般情况下，胞格的最小值位于理论配比附近，也即φ= 1，在缺氧一侧(即φ>1)，胞格尺寸的增长幅度小于富氧一侧(即φ<1)。下面的图4显示的是丙烷与空气混合气体在不同当量比的胞格尺寸。
图4 丙烷-空气混合气体在不同当量比例的胞格尺寸(p0 = 10 kPa)Fig.4 Cell size of propane-air mixture in different equivalent ratio (p0 = 10 kPa)
3.2 爆轰胞格尺寸与氩气稀释浓度的关系
惰性气体能明显增加胞格尺寸，随着惰性气体稀释浓度的增加胞格尺寸也相信增大。下图5显示的是丙烷-空气混合气体在加入Ar气体稀释后的胞格尺寸。
图5 丙烷-空气混合气体在不同Ar稀释浓度下的胞格尺寸(p0 = 10 kPa)Fig.5 Cell size of propane-air mixture under different Ar dilution detonation (p0 = 10 kPa)
3.3 爆轰胞格尺寸与初始压力的关系
胞格尺寸随着初始压力的升高而降低，呈指数衰减，并且胞格和初始压力基本上呈现的是线性关系。下面图6显示的是丙烷-空气混合气体在不同初始压力下的胞格尺寸。
图6 丙烷-空气混合气体在不同初始压力下的胞格尺寸Fig.6 Cell size of propane-air mixture under different initial pressure
4 实验结果与分析
通过大量的实验得到上述的实验结果，与其他的实验室(Murray、Knystsutas、Lee、Moen、Bull等)所做的实验结果基本吻合。下面采用定性推理的方式综合分析爆轰胞格尺寸随着各影响因素的变化情况，首先介绍两种定性的状态转换。
在定性模拟中有两种类型的状态转换：一种称为P-转换，即从显著时间点上到显著时间点之间的定性状态转换；另一种是I-转换，即从显著时间点之间到显著时间点上的定性状态转换。两种定性状态的转换可分别描述为：
P-转换：QS(f,ti)→QS(f,ti,ti+1)
I-转换：QS(f,ti,ti+1) →QS(f,ti)
下面我们解释P-与I-转换：
P-转换：首先我们假设f是物理系统的一个参数，并且f在时间点ti的定性状态为
I-转换：f在时间区间(ti,ti+1)的定性状态为
从(lj,lj+1)增加到lj+1，QDIR也从增加改变为保持不变。
定性状态的划分，是对可推理函数连续可微性的充分利用。在系统时间被表示为一个点与区间相隔出现的时间序列前提下，系统各参数的状态可按照P-转换，I-转换，P-转换，……的顺序进行转换，从而实现定性推理。
由上述的两种转换组成的通用函数状态转换如下表1所示，这个表格所列举的连续可微函数是所有可能的状态转换，根据这个表格可以推导出系统中个参数的后记状态，进行可以推测出整个系统的行为目的。
根据上面所介绍的P-和I-转换，再结合流定性相空间的定义，我们把爆轰波胞格尺寸随着气体当量变化进行了定性的处理，对于每个点的变化趋势进行了逐个的分析和推理，得到了爆轰胞格尺寸的流定性相空间，如下图7所示。
图7 爆轰尺寸关于气体当量变化的定性相空间Fig.7 Detonation size of gas equivalent in the qualitative phase space
所谓的流定性相空间，指的是从复杂的动力学系统观察，根据观察所得到的物流状态点，然后构造出系统行为定性的相空间图，系统每个可能的状态都有一相对应的相空间的点。通过P转换和I转换的计算和推理，得到了图7所示的爆轰胞格尺寸关于气体当量的定性相空间。通过P6转换，在初始点胞格尺寸是下降减少(dec)的趋势，在相邻的区间内胞格尺寸的定性状态的dec下降的，再结合I6变换，在初始点和相邻区间爆轰胞格都是减少dec的趋势，结合I6变换得知在相邻区间的边界点也是下降dec的趋势。通过不断的P和I定性转换，再结合相应的定量计算，可以得到如下的定性相空间。
表1 通用函数状态转换表
注：其中I8和I9转换发现了新的路标值l*，lj< l*< lj+1
可以看出随着气体定量的不断增加，爆轰胞格尺寸一开始是逐渐减少的，这个时候丙烷气体含量较少，而氧气含量较高，胞格尺寸与当量是呈现出线性递减的关系。当达到最佳配比的时候，胞格尺寸也到达最小值，这个时候丙烷气体和氧气的比例达到最佳理论值。当丙烷当量超过这个最佳比例继续增加的时候，爆轰胞格尺寸开始呈现出线性递增的趋势。
同样的，我们分析了爆轰胞格尺寸与Ar气体稀释浓度之间的定性关系如下图8所示。
图8 丙烷爆轰胞格尺寸与Ar气体浓度间的定性相空间Fig.8 Qualitative phase space of propane detonation cell size and Ar gas concentrations
从上图我们可以看出，丙烷爆轰胞格尺寸随着Ar气体浓度的不断增加呈现出不断增大的趋势。从初始的Ar浓度开始，胞格尺寸呈现出的是上升递增趋势inc，首先从P4转换，我们可以得到在相邻的区间内丙烷胞格尺寸是随着Ar浓度增加而递增；再结合I4转换和I5转换，可以得到在初始点相邻的连续区间内和区间的边界点胞格尺寸都是随着Ar气体浓度的增加而不断增加的。不断的运用P转化和I转换，就能够得到上图所示的定性相空间图。
上面仅仅是每次都针对某一个因素进行的定性分析，当多个因素共同变化一起对爆轰胞格尺寸产生影响的时候，就需要综合的考虑各个因素产生的叠加效果。下面的定性矩阵的分析能够很方便的处理多因素的共同作用效果。
结合上文的分析，爆轰胞格尺寸随着丙烷气体当量的变化趋势是先减少到最佳比例，定义λ= Fun((0，φ)-)，然后开始随着当量的增加而递增，定义λ= Fun((φ，+∞)+)；爆轰胞格尺寸随着Ar气体浓度的增加而不断增加，定义λ= Fun(Ar+)，；随着初始压力的增加而逐渐减少，定义为λ= Fun(p-)，其中(+)表示λ与影响因素有一个定性正比的关系；(-)表示λ与影响因素有一个定性反比的关系。将这些定性的变化趋势体现到下面的定性矩阵。下面矩阵中的“1”表示的是该初始条件是递增的、“0”表示的是保持不变，“-1”表示的是减少(见表2)。
表2 初始压力p和Ar对爆轰胞格尺寸λ的综合影响Table 2 Initial pressure p and Ar detonation cell size λλ = Fun(p-)⊕Fun(Ar+)
表2中“inc”表示增加，“inc*”表示多个因素的作用都是inc的，产生的效果会叠加增强，“dec”表示减少，“dec*”表示多个因素的作用都是dec的，“steady”表示“保持稳定不变”，“？”表示的是两种相异影响因素共同作用导致影响的结果难以判断，无法直接获得爆轰胞格尺寸是直接增加还是减少。不过为了判断两个相异影响效果的因素共同作用结果，我们可以给每种影响因素添加一个变量——权值strength()，表示的是各个影响因素作用于爆轰波胞格尺寸的影响大小，例如：如果设定strength(Fun(p-))>strength ( Fun(Ar+))，那么上面的表2可以进一步计算成下面的表3所示。
表3 初始压力p和Ar对爆轰胞格尺寸λ的综合影响1)Table 3 Initial pressure p and Ar detonation cell size λλ = Fun(p-)⊕Fun(Ar+)
1)strength(Fun(p-))>strength(Fun(Ar+))
从上表就可以看出初始压力和Ar气体的稀释共同作用于爆轰胞格尺寸的效果，由于strength(Fun(p-))>strength (Fun(Ar+))，在表2当中无法确定的几个元素在表3中都得到了相应的求解值。当初始压力增大时，胞格尺寸呈现出减少的趋势，但是同时Ar稀释浓度增加却对胞格尺寸产生增大的趋势影响，两者结合则由于初始压力p影响strength大而总体呈现出dec的趋势。在表3中当初始压力减少时，爆轰胞格尺寸会逐渐增大，如果此时Ar的稀释浓度增加，同样也会促使爆轰胞格尺寸增大，两个因素共同作用的效果都会促使爆轰胞格尺寸增加，这个时候产生的效果会是叠加增强的“inc*”。
当影响的初始条件比较多的时候，可以采用两两结合的方式，将计算的结果再与剩下的变量计算。当出现某些影响因素之间的strength值无法明显的比较大小时，可以采取定量的实验办法来确定在相同等级的变化幅度内，作用于胞格尺寸的影响效果强弱。
定性矩阵的分析，得到的结果虽然体现出来的都是dec、inc和steady的定性趋势，但是同样可以通过定量的实验室验证，将这些定性的结果转变成定量的实验数据。例如，inc产生的效果是100的数据级，那么inc*可以看成是200的数量级别。当然，具体的定性与定量间转化，需要实验验证。
当研究者希望得到爆轰胞格尺寸在某个范围内变化，或者说希望看到爆轰胞格结构在某个范围内变化的时候，可以通过定性矩阵的分析结果，适当调整影响 爆轰胞格尺寸的变量因素，将胞格尺寸调整到可控制范围内，例如“研究者希望将胞格尺寸从当前的10 mm调整到30 mm，可以采取降低初始压强，或者增加Ar的稀释浓度，也可以调整丙烷的当量，当然研究者可以同时调控多个影响因素，共同作用于爆轰胞格尺寸，使得胞格尺寸从当前的10 mm增加到30 mm。
5 总 结
预混气体爆轰胞格尺寸与多个因素相关，如何清楚的了解爆轰胞格结构的特点，就需要分析这些因素作用在胞格尺寸的不同效果。本文不仅从定量的实验中得到数据结果，同时结合到定性的分析研究出爆轰胞格尺寸与各个影响因素的定性关系，利用定性矩阵分析多个因素的叠加效果，得到的效果可以清楚地判断出胞格尺寸的变化是呈现出增加还是减少的趋势。通过定性与定量的实验结合验证，爆轰胞格结构的建模可以借助定性相空间的构造，推理和验证爆轰胞格尺寸与各影响因素之间的关系，这在理论上是个很大的创新，特别是当实验室环境下很难为实验室构造相应的实验条件时，定性仿真建模的方法很大程度上帮助研究者了解爆轰波的传播特性。
在爆轰的过程中预混气体的密度、压力和温度都是动态变化的，例如反应前期，随着化学反应的进行和能量的不断释放，化学反应区内的温度是不断提高的，但是反应后期随着化学反应速度放缓，放出的热量也在不断的减少，温度也肯定会相应的降低。对于这样的一种时刻动态的变化过程，采取定性的分析方法比定量的计算更加便捷，实验人员能够通过定性的分析实验环境下不同的初始和变化条件，进而解释爆轰波传播规律。
参考文献：
[1] THOMAS G O, WILLIAMS R L I. Detonation interaction with wedges and bends[J]. Shock Waves, 2002, 11: 481-492.
[2] OHYAGI S, OBARA T, NAKATA F, et al. A numerical simulation of reflection process of a detonation wave on a wedge[J]. Shock Waves, 2000, 10: 185-196.
[3] BARTLMA F, SCHRODER K. The diffraction of a plane detonation wave at an convex corner [J]. Combustion and Flame, 1986, 66: 237-248.
[4] EDWARDS D H, THOMAS G O. The diffraction of detonation waves in channels with 90° bends [J]. Combustions, 1983,3(1): 123-126.
[5] 张德良，谢巍，胡湘渝.气相爆轰波马赫反射的实验和数值研究 [C].第十届全国计算流体力学会议论文集, 2000: 381-386.
[6] HUANG Zhongwei, XU Bin. Relation between cell size and detonation parameters in overdriven gaseous detonation[J]. Combustion and Flame, 1987, 67: 95-102.
[7] 夏昌敬, 周凯元, 沈兆武, 等.可燃气体非稳定爆轰波通过90°圆弯管传播特性的实验研究 [J].实验力学, 2002, 17(4): 438-443.
[8] 夏昌敬, 周凯元.气相爆轰波在90°矩形弯管中传播时胞格结构的演化 [J].爆炸与冲击, 2005, 25(2): 151-156.
[9] 夏昌敬, 周凯元, 沈兆武.初始条件影响气体非稳定爆轰波在弯管中传播特性的实验研究[J].中国科学技术大学学报, 2004, 34(1): 92-97.
[10] LEE J H.Annual rev fluid mech[M].Pittsburge:the University of Pittsburge Press,1984:311-336.
[11] LEE J H,KNGSTAUTUS R, CHAN C K. Turbulent flame propagation in obstacle-filled tubes[C]∥20th Symposium (International) on Combustion,Pittsburgh:the Combustion Institute,1984:663-671.
[12] MOEN I O.Detonation in fuel-air explosive clouds[J].Hazard Mater,1993,(33): 159-192.
[13] BULL D C,ElSWORTH J E.Detonation cell structure in fuel-air mixture [J].Combust Flame, 1982, 45: 7-14.
[14] KUIPERS B. Qualitative simulation [J]. Artificial Intelligence, 1986，29(33): 289-338.
[15] 白方周，张雷. 定性仿真导论[M]. 合肥：中国科学技术大学出版社, 1998: 63-66.
