一种致盲弹辐射能量数值仿真系统及方法

一种致盲弹辐射能量数值仿真系统及方法
【专利摘要】本发明公开了一种致盲弹辐射能量数值仿真系统及方法。所述仿真系统由化学反应动力学计算模块、传热计算模块、各层颗粒吸收和衰减系数计算模块、各层光程计算模块四部分构成。本发明计算了药剂燃烧过程中的7步子反应的能量吸收与释放，计算了燃烧“火球”各层之间能量传递过程。计算获得的燃烧“火球”空间尺度、光谱辐射能量均与实验数据匹配良好，并且得到了连续易于分析的光谱辐射能量数据分布。与传统的应用实验测得“火球”表面温度再用黑体辐射定律计算的方法相比，省去了实验步骤，提高了计算效率、计算精度和可靠性。该方法为强光致盲弹燃烧过程仿真提供了一种新思路，可作为工程计算的一种有效模型和方法。
【专利说明】一种致盲弹辐射能量数值仿真系统及方法

【技术领域】
[0001]本发明属于含能材料燃烧光谱辐射特性计算仿真领域，涉及不同配方的警用致盲弹燃烧过程中的光谱辐射能量数值仿真系统及方法。

【背景技术】
[0002]闪光弹(Flashbang)，又称致盲弹、炫目弹或炫晕弹等，是一种以强光阻碍目标视力功能的一种轻型非致命武器，属于手榴弹的一种，为战术性的辅助工具之一。闪光弹被投掷后，会于数秒内发射刺眼强光，可以致使被攻击目标于短时间内有短暂性失明，使得目标顿时丧失反抗能力。除了以人为目标外，闪光弹亦有被用以投掷至坦克上光学器材的膜层，致使探测器失去了探测能力；另外，闪光弹亦有干扰敌人的战术用途。由于其携带方便，作用广泛，强光弹经常被用于警用、军用等许多方面。
[0003]强光致盲弹药本身蕴藏有很高的化学能，通过燃烧反应部分转变成光能，并以强辐射光的形式释放，达到毁伤敌人目的。这种弹药具有致人眩晕失明失去抵抗能力和使光电探测器材“致盲”不能正常工作的特点。强光致盲剂热辐射诱导下冲击激励发光强度参数的计算，对于指导强光致盲剂配方设计及强光致盲弹弹体设计、装药设计以及点火方式等都有十分重要的意义。热辐射诱导下的发光或各向均匀强光是它与普通爆炸闪光、爆炸激励定向强激光和红外光的发出是有区别的。因此，对强光致盲弹的爆炸过程进行仿真分析是很有必要的。
[0004]传统的强光致盲弹模型的建立均是以实验为基础，将致盲弹燃烧过程近似为黑体辐射，进而通过实验测得致盲弹燃烧产生“火球”表面温度来计算致盲弹的辐射能量，该方法是不准确的而且是非时变的，在某些精密条件下是不适用的。


【发明内容】

[0005]本发明目的是提供一种致盲弹辐射能量数值仿真系统及方法，基于含能材料燃烧爆炸“火球”形成机理研究，根据化学反应动力学理论、燃烧理论、传热理论、黑体辐射理论，考虑强光致盲弹化学药剂燃烧过程微观化学反应建立的燃烧、传热模型，能够精确地计算燃烧过程中化学反应各个步骤释放能量，从而计算得出致盲弹燃烧过程中某一时刻、各个波段的辐射能量及“火球”(辐射源)的空间尺度。
[0006]本发明的目的是通过以下技术方案实现的:
[0007]—种致盲弹辐射能量数值仿真系统，包括化学反应动力学计算模块、传热计算模块、各层颗粒吸收和衰减系数计算模块和各层光程计算模块，其中:
[0008]化学反应动力学计算模块，用于通过化学反应动力学对药剂的燃烧过程进行计算，得出药剂化学反应过程中各个时刻反应物、生成物、中间产物的质量与体积；
[0009]传热计算模块，用于利用化学反应动力学计算模块得到的各个时刻反应物、生成物、中间产物的质量与体积通过对燃烧“火球”各层之间热传递、气体膨胀做功计算得到“火球”各层温度，从而可以计算出“火球”内各层微粒光谱辐射能量；
[0010]各层颗粒吸收和衰减系数计算模块，用于利用化学反应动力学计算模块得到的各个时刻反应物、生成物、中间产物的质量与体积计算燃烧“火球”内部各层的颗粒吸收系数和裳减系数；
[0011]各层光程计算模块，用于计算燃烧“火球”内部某层的某一颗粒向四周发射的各条光线经过“火球”各层的光程，并与各层颗粒吸收和衰减系数计算模块计算得到的燃烧“火球”内部各层的颗粒吸收系数和衰减系数相结合得到各层光谱辐射透过率。最后通过对之前模块计算得到的数据结合计算出致盲弹燃烧过程中某一时刻、各个波段的辐射能量。
[0012]一种利用上述仿真系统进行致盲弹辐射能量数值仿真的方法，包括如下步骤:
[0013]步骤一、用户输入药剂质量，以化学反应药剂微粒为研究对象，将研究对象分为内层、中间层和外层，通过化学反应动力学计算模块对药剂的燃烧过程进行计算，得出药剂化学反应过程中各个时刻反应物、生成物、中间产物的质量与体积；
[0014]步骤二、以燃烧“火球”为研究对象，将“火球”分为Ns层，将步骤一得到的药剂化学反应过程中各个时刻反应物、生成物、中间产物的质量与体积代入传热计算模块，通过下式计算得到“火球”各层温度Tk, i ；
/ — /V — V/ — V
_ _fc—l.j—I k,i—l k,i—l k+?,?—Ι \ j / / γ 7 τ \ tt/ orkk j_i
[0015]AM-------)_ A- (K]_ KJk-ι)_ / —妒=mkck—τ:—
tsxtaAt ,
[0016]其中，Ak为“火球”内部第k层表面积，Tk^1为第k-ι层在1-ι时刻的温度，Tkji为第k层在i时刻的温度，h^, H和Jlri, H分别为第第k-Ι层在1-Ι时刻的扩散热焓和扩散通量，hk, i和Jk, i分别为第k层在i时刻的扩散热焓和扩散通量，Wrad和W分别为辐射和对外做功损失的能量，mk、ck分别为第k层气体质量和比热容；
[0017]由得到的某时刻“火球”内各层温度Tu可以计算出该时刻“火球”内各层微粒光谱辐射能量QkiiU):
[0018]
η I，、 2/?c.2 I
?,? (2) -he5
[0019]其中，h为普朗克常数，c为光速，λ为波长，ks为波尔兹曼常数；
[0020]步骤三、将步骤一得到的药剂化学反应过程中各个时刻反应物、生成物、中间产物的质量与体积代入各层颗粒吸收和衰减系数计算模块，通过以下公式计算出弥散体系中各层颗粒吸收系数ε和衰减系数α，并与各层光程计算模块计算得到的燃烧“火球”内部某层的某一颗粒向四周发射的各条光线经过“火球”各层的光程相结合得到各层光谱辐射透过率TkU):
[0021]ε k(A) = Qabs(A)；
[0022]ak(X) = 2Jir2NQextU)
[0023]
2U1)=Σ 哪[_巧(取]；



k=l
[0024]其中=Qabs(X)为各波段吸收效率因子，Qext(A)为各波段衰减效率因子，Cik(A)
为对应的消光系数，Ik为微粒在第k层经过的光程，N为本征载流子数密度，r为“火球”半径;
[0025]步骤四、将步骤二得到的各层微粒光谱辐射能量Qk, i ( λ )、步骤三得到的各层光谱辐射透过率TkU)结合以下公式得到燃烧“火球”对外辐射能量数据:
[0026]


I Λ*
α外取⑷;



iVft J.=l
[0027]
Q外(争$成外W.-




?=1
[0028]其中:?为粒子所发出的光线数，QuU)为该层粒子光谱辐射能量，Qk# (λ)为该层粒子发出经“火球”吸收后辐射出的能量，Pk为各层粒子浓度，Q#.U)为“火球”对外光谱辐射能量。
[0029]本发明计算了药剂燃烧过程中的7步子反应的能量吸收与释放，计算了燃烧“火球”各层之间能量传递过程。计算获得的燃烧“火球”空间尺度、光谱辐射能量均与实验数据匹配良好，并且得到了连续易于分析的光谱辐射能量数据分布。与传统的应用实验测得“火球”表面温度再用黑体辐射定律计算的方法相比，省去了实验步骤，提高了计算效率、计算精度和可靠性。总之，该方法为强光致盲弹燃烧过程仿真提供了一种新思路，可作为工程计算的一种有效模型和方法。
[0030]本发明具有如下优点:
[0031]1、能够快速仿真获得燃烧过程中“火球”(辐射源)的空间尺度、光谱辐射能量数据。
[0032]2、能够获得连续易于观察的光谱辐射能量数据分布，比起以往计算获得的离散数据，本发明得到更加直观、精确、可靠地燃烧光谱辐射能量数据。
[0033]3、无需实验测量中间数据，采用纯仿真方法获得燃烧光谱辐射能量数据。
[0034]4、本发明可以推广到任意弥散体系燃烧光谱辐射量的计算。

【专利附图】

【附图说明】
[0035]图1为致盲弹燃烧辐射能量计算总体方案；
[0036]图2为化学反应研究对象分层示意图；
[0037]图3为分层示意图；
[0038]图4为坐标系示意图；
[0039]图5为光线传输方向不意图1 ;
[0040]图6为光线传输方向不意图2 ；
[0041]图7为A1+KC104光谱辐射曲线；
[0042]图8为A1+KC104辐射能量时变曲线。

【具体实施方式】
[0043]下面结合附图对本发明的技术方案作进一步的说明，但并不局限如此，凡是对本发明技术方案进行修改或者等同替换，而不脱离本发明技术方案的精神和范围，均应涵盖在本发明的保护范围中。
[0044]1、化学反应动力学计算模块
[0045]强光致盲弹燃烧过程，化学反应过程中还原剂Mg和Al是主要研究对象。下面以A1+KC104配方为研究对象介绍一下该模型的物理机理。
[0046]反应动力学仍以单个颗粒及周围气体为研究对象，并认为“火球”中的每一层温度分布均匀，因此同一层中各个颗粒的化学状态函数是相同的。
[0047]如图2所示，将研究对象分为内层、中间层和外层。内层为液态Al，主要发生蒸发反应，并与中间层发生对流换热；中间层发生各种气相化学反应，生成的最终氧化产物Al2O3凝结于内层斜下方，并释放大量热；外层不发生化学反应，仅与中间层发生传导换热。中间层和外层的主要物质成分均是氧气。
[0048]金属铝颗粒燃烧是介于液滴燃烧和碳颗粒燃烧之间的一种类型，过程十分复杂，生成大量不稳定的中间态产物，再经过多级分解或凝聚反应，最终变成终态的三氧化二铝(Al2O3)。为简化计算，根据铝本身的物性及可观察到的化学反应，燃烧应经历如下过程:
[0049]氧化剂生成反应:
[0050]KClO4 — KCl+02 (Rl)
[0051]药剂中的KClO4分解生成O2作为铝燃烧的氧化剂。由于KClO4分解较容易进行，认为该反应在起始时刻之前已经完成。并假设生成的O2均匀分布于反应核空间。
[0052]表面反应:
[0053]Α1ω—Al(g) (R2)
[0054]液态的铝颗粒表面会蒸发出气态铝分子，并吸收热量。铝蒸发潜热hA1,vap为
11.835K.J/go
[0055]招颗粒表面的蒸发速率(g.τα2.s-1)可由Hertz-Langmuir方程表示为:
M I
[0056](I)
ImRT
[0057]其中，Mai为铝分子量(kg) ；R为标准气体常数(J -moF1 -Γ1) ；Pd为颗粒外铝蒸汽压。
[0058]Pd可通过Kelvin方程来描述:
[0059]

4σ^\\
Wp(命)(2)
[0060]其中，P。为温度T(K)时，平整表面的铝蒸汽压；σ @为熔融状态下铝液表面张力；V1(His)为铝分子体积；d(m)为铝粒子直径；k为波尔兹曼常数。
[0061]蒸汽压P。和表面张力σ @可由如下式子计算得到:
[0062]P0 = exp( 13.07- ~ψ~) (atm)(3)
[0063]σ 培=948-0.202Τ (mN/m) (4)
[0064]气相反应:
[0065]Al,,,.,+^-O2~^AlO(R3)
[0066]A10+^0:^->y\10,(R4)
[0067]其中，k3、k4分别为反应(R3)和反应(R4)的速率系数，是温度的函数，表示为:

SO
[0068]1、=9.76x101' exp(—■—) cm 'mol *.v !( 5 )

100S
[0069]/ = 4.63x 10l4exp(--^-) cm.'morV-1(6)
[0070]反应(R3)中，生成Imol的AlO吸收的热量H3 = 91kJ/mol ;反应(R4)中，生成Imol的AlO2吸收的热量H4 = 38kJ/mol。
[0071]生成物AlO和AlO2是化学性质极不稳定、处于激发态的中间产物，会辐射出一定的能量。
[0072]解离反应:
[0073]AI2O,A1f^-O2(R5)
[0074]当温度较高时，反应(R5)反应速率加快，并吸收大量热量，使体系温度保持在氧化铝沸点(4000K)以下。体系很难达到使金属氧化物蒸发的温度，因为其分解所需热量大于使其温度上升至沸点以上所需的热量。反应吸收热量H5设为变量。速率系数匕写为:
[0075]I5 =9.72xlOi 'exp(- ) cm'mol 1V"1(7)
[0076]凝结反应:
[0077]2 AlO+O2—^Al2O,,,.(R6)
[0078]AlO2+AlO2 A1,0.(；) + ^-O2(R7)
[0079]反应释放出大量热,反应(R6)中，生成Imol的Al2O3释放的热量H6 = 1831.849kJ/mo I ;反应(R7)中，生成 Imol 的 Al2O3 吸收的热量 H7 = 1725.8kJ/mol。
[0080]凝结反应过程假设分为两步进行:
[0081 ] aCm + bCn r-~5? Al2O3^j — ~> Al2O3^j(R8)
[0082]第一步主要是由AlO和AlO2反应生成Al2O3 ;第二步发生气态Al2O3迅速凝固为液态聚团(小液滴)。气态Al2O3的存在是有问题的，所以反应速率主要依赖于气态变成液态氧化物的凝结过程。
[0083]根据质量作用定律，第一步、第二步反应速率可写为:
[0084]~=kr(l('(8)
[0085]o82 = r_Cc (9)
[0086]其中r_是凝结成核速率系数。根据经典成核，rcon可写为:
[0087]




(?ο)

kmT πηι ρ
[0088]Ni*是单位体积内半径等于临界半径的胚团数(当胚团半径大于临界半径时，胚团才能平稳生长，形成稳定的结晶核)。
[0089]
N* = N1 exp[--了平 ^~-](11)
'1 3^3Γ3(1η5)2V J
[0090]其中，m是分子质量；P是液体密度；σ〒是平坦液面表面张力，可由(4)式表示；a j是冷凝系数；v是分子体积；S是过饱和度；队单位体积内是未凝结成胚团的单原子数。过饱和度S可由下式计算:
Σα
[0091]5 = 1 + ^(12)

Pa'
[0092]其中，Pi是种类i的偏压，i = A10, AlO2。
[0093]根据理想气体状态方程PV = η物RT，η物为物质的量，可得:
[P:[ Il
[0094]^ = 1 + ^- = 1 + ^(13)
Pm?Al
[0095]N1可通过如下方法进行估算:
[0096]
/=_2arnTmK(14)

Δ"Δ7.
[0097]其中，σ @为熔融状态下铝液表面张力，用(4)式表达；Tm为Al2O3沸点；VS为分子体积；ΛΗ为分子蒸发热，S卩ΛΗ = mH(H是蒸发一克Al2O3所吸收的热量)；AT是过冷度。
[0098]现在通过两个简单的假设来获得均匀形核的计算公式。假设一，于形核开始后临界胚团数Ni*仍保持平衡时的量值；假设二，一个原子进入临界胚团，即可使超过临界尺寸，并迅速生长。于是，单位体积胚团生长速度I可写为:
[0099]I=NlWf(卜)
[0100]
A =^exp(-(16)


波
[0101]ω*是围绕一个临界胚团的原子数目汀是能越过界面的原子频率数。对于球形胚团，ω*和f可由公式求出:
[0102]ω =~—(17)

cr
[0103]/ = j^(18)

(T
[0104]a是原子间距，D是气态扩散系数。
[0105]Turnbull D和Fisher J C在早期研究蒸汽中液滴形核的工作中，导出了 I的计算表达式:

[0106]
J-Bd^ cxpr l^wXv' 1π9)
Dlm pI 3AH2AT2^T1()
[0107]Dlffl是平衡熔点Tm时液体扩散系数。对于液态金属1，而B可取133的数值。



Ulm
[0108]因此，N1可由下式计算得到:
[0109]N1G^v=133( 20)
[0110]取分子半径l(T9m，Tm = 3800K, H = 1860KJ/mol, ΔΤ = 460K,分子间距 a =2Xl(T8m，扩散系数D = 2Xl(T5m2/s，由上式可计算得到N1 = 1.72 X 121。
[0111]反应(R8)的总速率ω8可写为:

Ir L.r Cl' Cb
m=~-—= f r c°" n(21)
[0112]^ 丄 + 丄 KC:C^ CcTcon
% W
[0113]根据文献，kr与r_成正比，则《8可简化为:
[0114]^=KrcmClCi(22)
[0115]如图3所示，为描述第一阶段单个金属粒子的化学反应过程，将反应区域分为内层、中间层和外层。
[0116]内层:该层是液化的铝颗粒(液滴)，密度和温度分布均匀，蒸发出大量铝蒸汽，并吸收热量。
[0117]中间层:化学反应R3?R8均在此进行，并释放出大量热，作为体系温度上升的能量来源。反应生成的最终氧化产物在中间层与内层的边界处凝固。在反应进行的过程中，中间层的温度始终高于内层和外层。一般认为铝蒸汽的氧化反应进行得较为迅速，中间层主要成分是氧气。
[0118]外层:不发生任何化学反应，与外界和中间层进行能量交换。主要成分是氧气。
[0119]定义变量如下，分别对内层、中间层和外层反应动力学方程。
[0120]pA1-招液滴密度
[0121]Reg——雷诺数
[0122]Eai—招液滴表面每秒蒸发量(g/s)
[0123]pg-氧气密度(g/m3)
[0124]d-液滴直径(m)
[0125]Ug-液滴表面氧气运动速度(m/s)
[0126]A-液滴表面积(m2)
[0127]Ud-招液滴运动速度(m/s)
[0128]pA1-铝液滴密度(g/m3)
[0129]kv—氧气膨胀系数(1/K)
[0130]T1——招液滴温度(K)
[0131]T3-夕卜层氧气温度(K)
[0132]Iii1-招液滴质量(g)
[0133]Cd——阻力系数
[0134]Cai——招比热容(J/ (K.g))
[0135]heva-招蒸发热(J/g)
[0136]μ g——氧气粘度(Pa.s)
[0137]Sa1203 (t)——中间层反应剩余的Al2O3质量(g)
[0138]Cg-氧气比热容(J/ (K.g))
[0139]Nud—努赛尔数
[0140]ε A1——招发射率
[0141]ο-Stefan-Boltzmann 常数
[0142]Pr——普朗特数
[0143]T2——中间层温度(K)
[0144]Sai (t)-中间层反应剩余的Al质量(g)
[0145]1、内层
[0146]内层主要发生液滴表面铝的蒸发反应。Eai表示为:
[0147]Eai = ωΑ1Α (23)
[0148]其中，ωΑ1是蒸发速率，可由式I?4得到；A是液滴表面积。A是时间的函数，选取液滴初始半径为15 μ m，初始质量为3.342X 10_14g，t时刻A和d通过下式计算得到:
"、^ 10.026x10 14
[0149]Α(0 = 4π{-—)3(24)
? ^/10.026xl0~!4-3/:\,/' ；
[0150]d = 2(-生”(25)

4^Pm
[0151]根据文献,招液滴密度Pai可表示为:
[0152]P A1 = 2364000 [1-0.0002268 (Tr933) ] (g.m3) (26)
[0153]内层的热量交换主要来自于与中间层氧气的对流换热，自身热辐射散热，以及蒸发吸收的热量。氧气通常认为不辐射也不吸收热辐射，因此，内层不接收中间层和外层的辐射热。所以，可得到能量传递方程如下:
wIcAl ~77 = Qconv ~ Qrad ~ Qhhmoll
at,一、
[0154](27)
=^^(T2-T})md -εΜσΤ^-EJwa

P、
[0155]其中Hi1 = 3.342Χ 10—8g_tEA1。
[0156]努赛尔数可通过雷诺数和普朗特数表示[2]:

I I
[0157]Y" =2[1 + Re,] Pr^/ 3](28)
[0158]将铝液滴直径作为特征长度，结合液滴与氧气的相对速度，可求得雷诺数:
[0159]Red = d P g I Ug-Ud | / μ g (29)
[0160]氧气的运动速度Ug主要来自于反应核气体的整体热膨胀，铝液滴速度Ud主要来自于气体膨胀过程中对其的推力(以地面为参考系是推力，以液滴为参考系是阻力)。
[0161]反应核发生化学反应时，产生热量，使得反应核中氧气膨胀。由于反应核处于高空自由空间，认为反应核周围压强不变，因此，根据气体热膨胀定律，反应核中气体当前温度为T，初始温度为Tc1、初始体积为Vci时，气体体积(即反应核体积)为:
[0162]V = V0 [ I+kv (T-T0)] (30)
[0163]由体积可得到半径，取半径的时间导数，即可得到以地面为参考系时，氧气的运动速度Ug (取初始体积0.0335m3，初始温度930K，kv = 0.00367):
[0164]u =0,3224^1 + 0-00367(Γ-^1(31)
g di ’
[0165]铝液滴主要受到气体膨胀过程中对其的推力，根据牛顿第二定律:
[0166]Tnl = pgd2CD(u - u )(32)
Qt 8
[0167]阻力系数Cd可通过下式求得:
[0168]Re^1000(33)
i 0.424 Re" >1000
[0169]2、中间层
[0170]内层表面的铝蒸汽蒸发至中间层后，与氧气发生R3反应。根据质量作用定律，Al的消耗速率和AlO的生成速率(g/cm3)为:
[01711 气和=2ikd 气細=39k3CAICa}(34)
[0172]其中CA1、CA1。为 Al 和 AlO 的摩尔浓度(mol/cm3)。
[0173]设t时刻,反应后，中间层Al剩余量为Sai (t),某一时刻Al剩余量的该变量等于该时刻内层蒸发的Al减去R3反应消耗的Al，所以SA1(t)服从方程:
必Al —fT _L.V (￡1-Y
[0174]dt A1 Al(32}(35)
、 ^ (0) = 0
[0175]某一时刻反应R3消耗的O2的质量1? (g)及产生的热量Qk3 (J)为:
[0176]/M02 =1.391(3.342x10 14(36)
[0177]Qe3 = -3956.5 (3.342 X I(T14-Hi1-Sai) (37)
[0178]反应R4、R5的产物消耗速率和产物生成速率仍可通过质量作用定律得到。同理，t时刻，反应R4、R5和R6进行后，AlO剩余量为Sa1⑴服从方程:
O IΠ I
0AlQ _ I V (L±\2 C (1-1Λ2
[0179]dt ~ 3 Al(32} 4 a1o(32}(38)
、5'Alo(0) = 0
[0180]反应生成的AlO2的质量mA1Q2(g)及产生的热量Qk4 (J)为:
[0181]/Wa1O2 =1.41(/wa10-5a10)(39)
[0182]Qe4 = -974.36 (mA10-SA10) (40)
[0183]中间层与外层的热交换形式为传导换热。通过联立以上方程组，即可解出所需要的物理量。
[0184]2、传热计算模块
[0185]从起始时刻起，最外层会与外界发生热交换，从而导致各层出现温度梯度，发生热传导现象；各层温度不同，会导致化学反应进行的程度不同，对反应气体的消耗程度不同，因而各层气体浓度发生变化，出现浓度梯度，发生扩散热交换；同时，“火球”在膨胀，对外做功，并且对外辐射能量，导致各层温度降低。各层传热现象可用以下方程描述:
[0186]- ~Wl) - 4 H - Vh ) - , = ^1 —二
AxAxAt
[0187](41)
[0188]其中，A为该层(第k层)表面积。Th, 为第k_l层在i_l时刻的温度，Tk, i为第k-Ι层在1-Ι时刻的温度。h和J分别为扩散热焓和扩散通量。Wrad和W分别为辐射和对外做功损失的能量。mk、Ck分别为该层气体质量和比热容。
[0189]通过公式(41)计算得到某时刻“火球”内各层温度Tu，由此可以计算出该时刻“火球”内各层微粒光谱辐射能量Qk, i (入)。
[0190]
QkM=(42)


e^TkJ _ I
[0191]其中h为普朗克常数，c为光速，λ为波长，ks为波尔兹曼常数。
[0192]3、各层颗粒吸收和衰减系数计算模块
[0193]燃烧产物在高温下会电离出自由电荷，同时产生从Al或O原子中电离出的空穴。产物连续发射谱归因于本征载流子。因此，可以通过Drude-Lorentz模型来求解燃烧产物的复折射率虚部，进而可以得到发射率。
[0194]类比Lorentz色散模型，自由电荷运动方程可写为:
[0195]m + ηιΓ^^- = -eK{!)(43)
drdiv 7
[0196]方程解为:
r ?'t \ ^? Λ Λ \
[0197]r(!) = ——~——(44)

m (?τ + ?Γω)
[0198]每个电子提供的原子偶极距为P’ = _eXr，则介质宏观极化强度为:
[0199]P = NXP，= ε 0xE (45)
[0200]其中，N为载流子数密度。由此可得电极化率X:
[0201 ] χ( ω) = 一~~(46)

ur + /'Γ ω
[0202]其中，ωρ2 = Ne2/ ε 0m, ε0为真空介电常数，m为载流子质量。进一步可得:
[0203]χ\ω) = -ω' 、χ'\ω) = αξ，ω.(47)

or + 厂-F ar -?- Γ*.
[0204]由 1+ X = ε 得:
[0205]?::{ω) = \-ωΙ , 1 ” ε.”(ω) = ω?「出,(48)
ρ ω +TF or +Τ~
[0206]其中，Γ = τ-1，τ为产物粒子弛豫时间。
[0207]对于粒子复折射率，n2_k2 = ε / r,2nk = ε" y ￡ = (<2+f2)1"。粒子复折射率虚部可写为:
[0208]k虚=[(ε-ε' r)/2]1/2 (49)
[0209]本征载流子数密度N:
[0210]通过固体物理方法的推导，可得本征载流子数密度为:
[0211]
N = 2{^-f \memhr 哪(S)(50)
2πη2k波 I
[0212]其中，ks为玻尔兹曼常数，me、mh分别为电子和空穴有效质量，T为温度，h为约化普朗克常数，将式(50)代入(46)，再代入(48)，即可求解出发射率。带隙能Eg写为:
[0213]Eg = 9.1eV- (0.0011eV/K) X T-Epol (51)
[0214]其中，Epq1 = 2/3eV(eV-电子伏特，K-开尔文)。
[0215]不同温度，不同波长对应不同的Iije值。
[0216]Al2O3复折射率实部k实的计算:
[0217]值可通过以下经验式进行计算，
[0218]
兮_ ? 二 Σ為义2 / (儿2 _几)
A2 = 0.00377588, A =1.023798,門、
H1(52)
[0219]Λ: = 0.0122544, A2 = 1.058264
為2 =321.3616, A3 =5.280792
[0220]通过上式即可得到不同波长下的折射率实部值。
[0221]吸收和衰减效率因子的Mie理论计算:
[0222]根据理论Mie，当光强为Itl，在颗粒周围介质中波长为λ的自然光平行入射到一半径为Y的各向同性球形颗粒上时，在散射角为Θ，距离散射体处的散射光强为:
[0223]




(53)
[0224]其中
[0225]I1 = ,.V, (m，θ, a) X..V,' (m, θ, a)
[0226]h - S1 (m..θ'a)XS;(m、Θ.a)(54)
[0227]S1 (Θ) = ^ + | (α?π? (cos Θ) + hntn (cos Θ))

n=i n(n + l)

ο? 2/1? I
[0228]‘S:狀)=Σ—~-(^r,(cos6>) + /7,^(cos^))
L J?=i η(η + ?)(55)

Τζ / )
[0229](X = lJ-
Λ(56)
?K (cos6?) (J!X
[cosθ) ^~——^—~~—
sin θ a [cosθ)
[0230]! / 、^
洲J(57)
[0231 ] Pm1为m阶勒让德多项式
α _ μην],, {mfDia jAa^ - μ J11 {(Χ)[ιηα]η {ma)^
η Min1J, (ma)[ahip (a)]' - μ Jii^ (a)[ma jr (Jna)X
b.....βχ?η {ma)[xjx (α)? - β-- {a)[m(xj? Qop]'
[0233]? //,./? (Ina)IxiC {a)~1- phlp (a ^rnajil (Iiia)J(58)
[0234]其中，α是尺寸参数，D为粒子球直径，m = k实-1kji为粒子球的复折射率，jn(x)为球贝塞尔函数，hna)(x)为第一类汉克尔函数。Pn(X)为勒让德多项式。
[0235]由计算得到的复折射率，和化学反应模块输出的粒子半径，通过Mie理论即可计算得到各波段吸收效率因子Qabs和衰减效率因子Qrartt5
[0236]Q- = Λ Σ (2〃 + 1Xlflfw Γ+ΚΓ)(59)
w ^=I
[0237]4 (2"十丨)Re(% + ^)(60)

^ --=1
[0238]从而，颗粒吸收系数ε和衰减系数α可分别写为:
[0239]ε (λ) = Qabs ( λ )
[0240](61)
[0241 ] α ( λ ) = 2 π T2NQext ( λ )
[0242]其中:N为本征载流子数密度，r为“火球”半径。
[0243]4、各层光程计算模块
[0244]该模块用于计算“火球”模型内任意一微粒所发出光线在“火球”内各层的光程。由于“火球”的旋转对称性，可以将“火球”模型抽象成二维模型，用一系列同心圆表示，以圆心为坐标原点可建立平面直角坐标系。示意图如图4:
[0245]图中最内层表示爆炸所产生的真空层，该层折射率为1，故不在计算范围内。
[0246]发光微粒可抽象为同心圆内任意位置的一点，由于圆环的对称性可将该点设定在坐标系中的X轴上，设该点为Xtl点，该点所在层数为第k层。该点发出的光线可抽象为以该点为端点的射线。设定一个角度步长，可以对0° -180°内任意角度进行光程计算。由圆环的轴对称性，只计算X轴上方两象限的光线即可，X轴下方的光线可由对称关系得到。
[0247]由于追迹的光线相对于X轴偏角不同，则可以将光程计算分为两种情况:光线与目标点以内的圆环相交的情况和光线与目标点以内的圆环相离的情况，用几何关系很容易表不。两种情况的不意图分别如图5和6:
[0248]当所追迹的光线与目标点以内的圆环相离时，则有光线在目标点以内的光程为0，只需计算光线在目标点以外的光程即可。如上图所示，已知ΟΧρΟΑ、Θ角，则可由余弦定理求出XtlA, X0A即为该光线在第i层内的光程。同理可以求出XtlB，二者之差即为该光线在第k+Ι层内的光程。以此类推可以求出该光线在以外各层内的光程。
[0249]当光线与目标点以内的圆环相交时，除真空层内光程可以忽略不计外，其余层内光程均需要计算。将直线方程与不同半径的圆方程联立即可求得二者交点，由交点坐标即可求出两交点间距离。如图6所示，BC即为真空层内光程，AD-BC为第一层内光程，以此类推，直到第k层。第k层以后的光程可由上一种情况的计算方法计算。
[0250]该模块计算得出的光程Ik与之前计算得到的消光系数经可以计算出“火球”内各层光谱透过率TkU )。
[0251]

【权利要求】
1.一种致盲弹辐射能量数值仿真系统，其特征在于所述仿真系统由化学反应动力学计算模块、传热计算模块、各层颗粒吸收和衰减系数计算模块、各层光程计算模块四部分构成，其中: 化学反应动力学计算模块，用于通过化学反应动力学对药剂的燃烧过程进行计算，得出药剂化学反应过程中各个时刻反应物、生成物、中间产物的质量与体积； 传热计算模块，用于利用化学反应动力学计算模块得到的各个时刻反应物、生成物、中间产物的质量与体积通过对燃烧“火球”各层之间热传递、气体膨胀做功计算得到“火球”各层温度，从而可以计算出“火球”内各层微粒光谱辐射能量； 各层颗粒吸收和衰减系数计算模块，用于利用化学反应动力学计算模块得到的各个时刻反应物、生成物、中间产物的质量与体积计算燃烧“火球”内部各层的颗粒吸收系数和衰减系数； 各层光程计算模块，用于计算燃烧“火球”内部某层的某一颗粒向四周发射的各条光线经过“火球”各层的光程，并与各层颗粒吸收和衰减系数计算模块计算得到的燃烧“火球”内部各层的颗粒吸收系数和衰减系数相结合得到各层光谱辐射透过率，最后通过对之前模块计算得到的数据结合计算出致盲弹燃烧过程中某一时刻、各个波段的辐射能量。
2.—种致盲弹辐射能量数值仿真方法，其特征在于所述仿真方法步骤如下: 步骤一、用户输入药剂质量，以化学反应药剂微粒为研究对象，将研究对象分为内层、中间层和外层，通过化学反应动力学计算模块对药剂的燃烧过程进行计算，得出药剂化学反应过程中各个时刻反应物、生成物、中间产物的质量与体积； 步骤二、以燃烧“火球”为研究对象，将“火球”分为Ns层，将步骤一得到的药剂化学反应过程中各个时刻反应物、生成物、中间产物的质量与体积代入传热计算模块，通过下式计算得到“火球”各层温度Tk, i; A^{k-u\TkJ-1 -- Ak— Vm ) - Wrad-W = mkck TkJ:卜-1 hxhocAt , 其中，Ak为“火球”内部第k层表面积，Tm1为第k-Ι层在1-Ι时刻的温度，Tka为第k层在i时刻的温度，和Jk^1分别为第k-Ι层在1-Ι时刻的扩散热焓和扩散通量，hu和Ju分别为第k层在i时刻的扩散热焓和扩散通量，Wrad和W分别为辐射和对外做功损失的能量，mk、ck分别为第k层气体质量和比热容； 由得到的某时刻“火球”内各层温度Tu可以计算出该时刻“火球”内各层微粒光谱辐射能量QuU):
,Ihc2 IQkj ) - Phe， 其中，h为普朗克常数，c为光速，λ为波长，波尔兹曼常数； 步骤三、将步骤一得到的药剂化学反应过程中各个时刻反应物、生成物、中间产物的质量与体积代入各层颗粒吸收和衰减系数计算模块，通过以下公式计算出弥散体系中各层颗粒吸收系数ε和衰减系数α，并与各层光程计算模块计算得到的燃烧“火球”内部某层的某一颗粒向四周发射的各条光线经过“火球”各层的光程相结合得到各层光谱辐射透过率Tk“):
ε k (入)=Qabs (入)
α k ( λ ) = 2 JI r2NQext (λ)；



k=l 其中:Qabs(x)为各波段吸收效率因子，QraitU)为各波段衰减效率因子，α,(λ)为对应的消光系数，Ik为微粒在第k层经过的光程，N为本征载流子数密度，r为“火球”半径；步骤四、将步骤二得到的各层微粒光谱辐射能量Qk, i ( λ )、步骤三得到的各层光谱辐射透过率TkU)结合以下公式得到燃烧“火球”对外辐射能量数据:QmU) ;
光片Q外;



?=1 其中:Ν|为粒子所发出的光线数，Qu(X)为该层粒子光谱辐射能量，Qk# (λ)为该层粒子发出经“火球”吸收后辐射出的能量，Pk为各层粒子浓度，Q#.(λ)为“火球”对外光谱辐射能量。
3.根据权利要求2所述的致盲弹辐射能量数值仿真方法，其特征在于所述本征载流子数密度N为:N =exp(-^-)，
2πη2k波 I 其中，ks为玻尔兹曼常数，me、mh分别为电子和空穴有效质量，Eg为带隙能，T为温度，A是约化普朗克常数。
4.根据权利要求2所述的致盲弹辐射能量数值仿真方法，其特征在于所述带隙能Eg写为:
Eg = 9.1eV- (0.0OlleV/K) XT-Epol, 其中，Eptjl = 2/3eV, eV-电子伏特，K-开尔文。
5.根据权利要求2所述的致盲弹辐射能量数值仿真方法，其特征在于所述各波段吸收效率因子Qabs写为: 0^=4?^?+1'χΚΓ+Ι6πΓ)，
^ n=l 其中，α为尺寸参数，an、bnSMie散射计算中间数据。
6.根据权利要求2所述的致盲弹辐射能量数值仿真方法，其特征在于所述衰减效率因子Qeixt写为:
Qext = ~Σ(2? +1)Re(-- +b?)。
^ n=\
【文档编号】G06F19/00GK104134009SQ201410392093
【公开日】2014年11月5日 申请日期:2014年8月11日 优先权日:2014年8月11日
【发明者】万子南, 李韬, 陈守谦, 王丽, 党凡阳, 王玉雷, 秦兰琦, 汪东生, 刘振奇, 范志刚 申请人:哈尔滨工业大学