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

摘要

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

说明书

技术领域

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

背景技术

闪光弹(Flashbang)，又称致盲弹、炫目弹或炫晕弹等，是一种 以强光阻碍目标视力功能的一种轻型非致命武器，属于手榴弹的一 种，为战术性的辅助工具之一。闪光弹被投掷后，会于数秒内发射刺 眼强光，可以致使被攻击目标于短时间内有短暂性失明，使得目标顿 时丧失反抗能力。除了以人为目标外，闪光弹亦有被用以投掷至坦克 上光学器材的膜层，致使探测器失去了探测能力；另外，闪光弹亦有 干扰敌人的战术用途。由于其携带方便，作用广泛，强光弹经常被用 于警用、军用等许多方面。

强光致盲弹药本身蕴藏有很高的化学能，通过燃烧反应部分转变 成光能，并以强辐射光的形式释放，达到毁伤敌人目的。这种弹药具 有致人眩晕失明失去抵抗能力和使光电探测器材“致盲”不能正常工 作的特点。强光致盲剂热辐射诱导下冲击激励发光强度参数的计算， 对于指导强光致盲剂配方设计及强光致盲弹弹体设计、装药设计以及 点火方式等都有十分重要的意义。热辐射诱导下的发光或各向均匀强 光是它与普通爆炸闪光、爆炸激励定向强激光和红外光的发出是有区 别的。因此，对强光致盲弹的爆炸过程进行仿真分析是很有必要的。

传统的强光致盲弹模型的建立均是以实验为基础，将致盲弹燃烧 过程近似为黑体辐射，进而通过实验测得致盲弹燃烧产生“火球”表 面温度来计算致盲弹的辐射能量，该方法是不准确的而且是非时变 的，在某些精密条件下是不适用的。

发明内容

本发明目的是提供一种致盲弹辐射能量数值仿真系统及方法，基 于含能材料燃烧爆炸“火球”形成机理研究，根据化学反应动力学理 论、燃烧理论、传热理论、黑体辐射理论，考虑强光致盲弹化学药剂 燃烧过程微观化学反应建立的燃烧、传热模型，能够精确地计算燃烧 过程中化学反应各个步骤释放能量，从而计算得出致盲弹燃烧过程中 某一时刻、各个波段的辐射能量及“火球”(辐射源)的空间尺度。

本发明的目的是通过以下技术方案实现的：

一种致盲弹辐射能量数值仿真系统，包括化学反应动力学计算模 块、传热计算模块、各层颗粒吸收和衰减系数计算模块和各层光程计 算模块，其中：

化学反应动力学计算模块，用于通过化学反应动力学对药剂的燃 烧过程进行计算，得出药剂化学反应过程中各个时刻反应物、生成物、 中间产物的质量与体积；

传热计算模块，用于利用化学反应动力学计算模块得到的各个时 刻反应物、生成物、中间产物的质量与体积通过对燃烧“火球”各层 之间热传递、气体膨胀做功计算得到“火球”各层温度，从而可以计 算出“火球”内各层微粒光谱辐射能量；

各层颗粒吸收和衰减系数计算模块，用于利用化学反应动力学计 算模块得到的各个时刻反应物、生成物、中间产物的质量与体积计算 燃烧“火球”内部各层的颗粒吸收系数和衰减系数；

各层光程计算模块，用于计算燃烧“火球”内部某层的某一颗粒 向四周发射的各条光线经过“火球”各层的光程，并与各层颗粒吸收 和衰减系数计算模块计算得到的燃烧“火球”内部各层的颗粒吸收系 数和衰减系数相结合得到各层光谱辐射透过率。最后通过对之前模块 计算得到的数据结合计算出致盲弹燃烧过程中某一时刻、各个波段的 辐射能量。

一种利用上述仿真系统进行致盲弹辐射能量数值仿真的方法，包 括如下步骤：

步骤一、用户输入药剂质量，以化学反应药剂微粒为研究对象， 将研究对象分为内层、中间层和外层，通过化学反应动力学计算模块 对药剂的燃烧过程进行计算，得出药剂化学反应过程中各个时刻反应 物、生成物、中间产物的质量与体积；

步骤二、以燃烧“火球”为研究对象，将“火球”分为N_层层， 将步骤一得到的药剂化学反应过程中各个时刻反应物、生成物、中间 产物的质量与体积代入传热计算模块，通过下式计算得到“火球”各 层温度T_k，i；

$A_k\lambda (\frac{T_{k-1,i-1}-T_{k,i-1}}{\mathrm{\Delta x}}-\frac{T_{k,i-1}-T_{k+1,i-1}}{\mathrm{\Delta x}})-A_k(h_{k-1}{J}_{k-1}-h_k{J}_{k-1})W_{\mathrm{rad}}-W=m_k{c}_k\frac{T_{k,i}-T_{k,i-1}}{\mathrm{\Delta t}},$

其中，A_k为“火球”内部第k层表面积，T_k-1，i-1为第k-1层在i-1 时刻的温度，T_k，i为第k层在i时刻的温度，h_k-1，i-1和J_k-1，i-1分别为第第 k-1层在i-1时刻的扩散热焓和扩散通量，h_k，i和J_k，i分别为第k层在i 时刻的扩散热焓和扩散通量，W_rad和W分别为辐射和对外做功损失 的能量，m_k、c_k分别为第k层气体质量和比热容；

由得到的某时刻“火球”内各层温度T_k，i可以计算出该时刻“火 球”内各层微粒光谱辐射能量Q_k，i(λ)：

其中，h为普朗克常数，c为光速，λ为波长，k_波为波尔兹曼常 数；

步骤三、将步骤一得到的药剂化学反应过程中各个时刻反应物、 生成物、中间产物的质量与体积代入各层颗粒吸收和衰减系数计算模 块，通过以下公式计算出弥散体系中各层颗粒吸收系数ε和衰减系数 α，并与各层光程计算模块计算得到的燃烧“火球”内部某层的某一 颗粒向四周发射的各条光线经过“火球”各层的光程相结合得到各层 光谱辐射透过率T_k(λ)：

ε_k(λ)＝Q_abs(λ)；

α_k(λ)＝2πr²NQ_ext(λ)

其中：Q_abs(λ)为各波段吸收效率因子，Q_ext(λ)为各波段衰减效率 因子，α_k(λ)为对应的消光系数，l_k为微粒在第k层经过的光程，N为 本征载流子数密度，r为“火球”半径；

步骤四、将步骤二得到的各层微粒光谱辐射能量Q_k，i(λ)、步骤三 得到的各层光谱辐射透过率T_k(λ)结合以下公式得到燃烧“火球”对 外辐射能量数据：

其中：N_光为粒子所发出的光线数，Q_k，i(λ)为该层粒子光谱辐射 能量，Q_k外(λ)为该层粒子发出经“火球”吸收后辐射出的能量，ρ_k为 各层粒子浓度，Q_外·(λ)为“火球”对外光谱辐射能量。

本发明计算了药剂燃烧过程中的7步子反应的能量吸收与释放， 计算了燃烧“火球”各层之间能量传递过程。计算获得的燃烧“火球”空 间尺度、光谱辐射能量均与实验数据匹配良好，并且得到了连续易于 分析的光谱辐射能量数据分布。与传统的应用实验测得“火球”表面温 度再用黑体辐射定律计算的方法相比，省去了实验步骤，提高了计算 效率、计算精度和可靠性。总之，该方法为强光致盲弹燃烧过程仿真 提供了一种新思路，可作为工程计算的一种有效模型和方法。

本发明具有如下优点：

1、能够快速仿真获得燃烧过程中“火球”(辐射源)的空间尺度、 光谱辐射能量数据。

2、能够获得连续易于观察的光谱辐射能量数据分布，比起以往 计算获得的离散数据，本发明得到更加直观、精确、可靠地燃烧光谱 辐射能量数据。

3、无需实验测量中间数据，采用纯仿真方法获得燃烧光谱辐射 能量数据。

4、本发明可以推广到任意弥散体系燃烧光谱辐射量的计算。

附图说明

图1为致盲弹燃烧辐射能量计算总体方案；

图2为化学反应研究对象分层示意图；

图3为分层示意图；

图4为坐标系示意图；

图5为光线传输方向示意图1；

图6为光线传输方向示意图2；

图7为Al+KClO₄光谱辐射曲线；

图8为Al+KClO₄辐射能量时变曲线。

具体实施方式

下面结合附图对本发明的技术方案作进一步的说明，但并不局限 如此，凡是对本发明技术方案进行修改或者等同替换，而不脱离本发 明技术方案的精神和范围，均应涵盖在本发明的保护范围中。

1、化学反应动力学计算模块

强光致盲弹燃烧过程，化学反应过程中还原剂Mg和Al是主要 研究对象。下面以Al+KClO₄配方为研究对象介绍一下该模型的物理 机理。

反应动力学仍以单个颗粒及周围气体为研究对象，并认为“火球” 中的每一层温度分布均匀，因此同一层中各个颗粒的化学状态函数是 相同的。

如图2所示，将研究对象分为内层、中间层和外层。内层为液态 Al，主要发生蒸发反应，并与中间层发生对流换热；中间层发生各种 气相化学反应，生成的最终氧化产物Al₂O₃凝结于内层斜下方，并释 放大量热；外层不发生化学反应，仅与中间层发生传导换热。中间层 和外层的主要物质成分均是氧气。

金属铝颗粒燃烧是介于液滴燃烧和碳颗粒燃烧之间的一种类型， 过程十分复杂，生成大量不稳定的中间态产物，再经过多级分解或凝 聚反应，最终变成终态的三氧化二铝(Al₂O₃)。为简化计算，根据铝 本身的物性及可观察到的化学反应，燃烧应经历如下过程：

氧化剂生成反应：

KClO₄→KCl+O₂    (R1)

药剂中的KClO₄分解生成O₂作为铝燃烧的氧化剂。由于KClO₄分解较容易进行，认为该反应在起始时刻之前已经完成。并假设生成 的O₂均匀分布于反应核空间。

表面反应：

Al_(l)→Al_(g)    (R2)

液态的铝颗粒表面会蒸发出气态铝分子，并吸收热量。铝蒸发潜 热h_Al，vap为11.835K·J/g。

铝颗粒表面的蒸发速率(g·m^-2·s^-1)可由Hertz-Langmuir方程表 示为：

${\omega }_{\mathrm{Al}}={\left(\frac{M_{\mathrm{Al}}}{2\mathrm{\pi RT}}\right)}^{\frac{1}{2}}·P_d---\left(1\right)$

其中，M_Al为铝分子量(kg)；R为标准气体常数(J·mol^-1·K^-1)； P_d为颗粒外铝蒸汽压。

P_d可通过Kelvin方程来描述：

其中，P₀为温度T(K)时，平整表面的铝蒸汽压；σ_熔为熔融状 态下铝液表面张力；v₁(m³)为铝分子体积；d(m)为铝粒子直径；k为 波尔兹曼常数。

蒸汽压P₀和表面张力σ_熔可由如下式子计算得到：

$P_0=\exp (13.07-\frac{36373}{T})\left(\mathrm{atm}\right)---\left(3\right)$

σ_熔＝948-0.202T  (mN/m)    (4)

气相反应：

${\mathrm{Al}}_{\left(g\right)}+\frac{1}{2}{O}_2\stackrel{k_3}{\to }\mathrm{AlO}---\left(R3\right)$

$\mathrm{AlO}+\frac{1}{2}{O}_2\stackrel{k_4}{\to }\mathrm{Al}{O}_2---\left(R4\right)$

其中，k₃、k₄分别为反应(R3)和反应(R4)的速率系数，是温度的 函数，表示为：

$\left(\begin{array}{cc}{k}_3=9.76\times {10}^{13}\exp (-\frac{80}{T})& {\mathrm{cm}}^3{\mathrm{mol}}^{-1}{s}^{-1}---\left(5\right)\end{array}\right)$

$\left(\begin{array}{cc}{k}_4=4.63\times {10}^{14}\exp (-\frac{10008}{T})& {\mathrm{cm}}^3{\mathrm{mol}}^{-1}{s}^{-1}---\left(6\right)\end{array}\right)$

反应(R3)中，生成1mol的AlO吸收的热量H₃＝91kJ/mol；反应(R4) 中，生成1mol的AlO₂吸收的热量H₄＝38kJ/mol。

生成物AlO和AlO₂是化学性质极不稳定、处于激发态的中间产 物，会辐射出一定的能量。

解离反应：

${\mathrm{Al}}_2{O}_3\stackrel{k_5}{\to }\mathrm{AlO}+\frac{1}{2}{O}_2---\left(R5\right)$

当温度较高时，反应(R5)反应速率加快，并吸收大量热量，使体 系温度保持在氧化铝沸点(4000K)以下。体系很难达到使金属氧化 物蒸发的温度，因为其分解所需热量大于使其温度上升至沸点以上所 需的热量。反应吸收热量H₅设为变量。速率系数k₅写为：

$\left(\begin{array}{cc}{k}_5=9.72\times {10}^{15}\exp (-\frac{25000}{T})& {\mathrm{cm}}^3{\mathrm{mol}}^{-1}{s}^{-1}---\left(7\right)\end{array}\right)$

凝结反应：

$2\mathrm{AlO}+\frac{1}{2}{O}_2\stackrel{k_6}{\to }{\mathrm{Al}}_2{O}_{3\left(l\right)}---\left(R6\right)$

$\mathrm{Al}{O}_2+\mathrm{Al}{O}_2\stackrel{k_7}{\to }{\mathrm{Al}}_2{O}_{3\left(l\right)}+\frac{1}{2}{O}_2---\left(R7\right)$

反应释放出大量热，反应(R6)中，生成1mol的Al₂O₃释放的热 量H₆＝1831.849kJ/mol；反应(R7)中，生成1mol的Al₂O₃吸收的热量 H₇＝1725.8kJ/mol。

凝结反应过程假设分为两步进行：

${\mathrm{aC}}_m+{\mathrm{bC}}_m\stackrel{k_r}{\to }{\mathrm{Al}}_2{O}_{3\left(g\right)}\stackrel{r_{\mathrm{con}}}{\to }{\mathrm{Al}}_2{O}_{3\left(l\right)}---\left(R8\right)$

第一步主要是由AlO和AlO₂反应生成Al₂O₃；第二步发生气态 Al₂O₃迅速凝固为液态聚团(小液滴)。气态Al₂O₃的存在是有问题的， 所以反应速率主要依赖于气态变成液态氧化物的凝结过程。

根据质量作用定律，第一步、第二步反应速率可写为：

${\omega }_{\mathrm{8,1}}=k_r{C}_m^a{C}_n^b---\left(8\right)$

ω_8，2＝r_conC_c    (9)

其中r_con是凝结成核速率系数。根据经典成核，r_con可写为：

N_i^*是单位体积内半径等于临界半径的胚团数(当胚团半径大于 临界半径时，胚团才能平稳生长，形成稳定的结晶核)。

其中，m是分子质量；ρ是液体密度；σ_平是平坦液面表面张力， 可由(4)式表示；α_i*是冷凝系数；v是分子体积；S是过饱和度；N₁单位体积内是未凝结成胚团的单原子数。过饱和度S可由下式计算：

$S=1+\frac{\underset{i}{\Sigma }{p}_i}{p_{\mathrm{Al}}}---\left(12\right)$

其中，p_i是种类i的偏压，i＝AlO，AlO₂。

根据理想气体状态方程PV＝n_物RT，n_物为物质的量，可得：

$S=1+\frac{\underset{i}{\Sigma }{p}_i}{p_{\mathrm{Al}}}=1+\frac{\underset{i}{\Sigma }{n}_i}{n_{\mathrm{Al}}}---\left(13\right)$

N₁可通过如下方法进行估算：

其中，σ_熔为熔融状态下铝液表面张力，用(4)式表达；T_m为Al₂O₃沸点；V_s为分子体积；ΔH为分子蒸发热，即ΔH＝mH(H是蒸发一 克Al₂O₃所吸收的热量)；ΔT是过冷度。

现在通过两个简单的假设来获得均匀形核的计算公式。假设一， 于形核开始后临界胚团数N_i^*仍保持平衡时的量值；假设二，一个原 子进入临界胚团，即可使超过临界尺寸，并迅速生长。于是，单位体 积胚团生长速度I可写为：

$I=N_i^{*}{\omega }^{*}f---\left(15\right)$

ω^*是围绕一个临界胚团的原子数目；f是能越过界面的原子频率 数。对于球形胚团，ω^*和f可由公式求出：

${\omega }^{*}=\frac{4\pi r^{*2}}{a^2}---\left(17\right)$

$f=\frac{D}{a^2}---\left(18\right)$

a是原子间距，D是气态扩散系数。

Turnbull D和Fisher J C在早期研究蒸汽中液滴形核的工作中， 导出了I的计算表达式：

D_lm是平衡熔点T_m时液体扩散系数。对于液态金属而B 可取10³³的数值。

因此，N₁可由下式计算得到：

N₁ω^*v＝10³³                          (20)

取分子半径10^-9m，T_m＝3800K，H＝1860KJ/mol，ΔT＝460K，分 子间距a＝2×10^-8m，扩散系数D＝2×10^-5m²/s，由上式可计算得到 N₁＝1.72×10²¹。

反应(R8)的总速率ω₈可写为：

${\omega }_8=\frac{1}{\frac{1}{{\omega }_{\mathrm{8,1}}}+\frac{1}{{\omega }_{\mathrm{8,2}}}}=\frac{C_c{k}_r{r}_{\mathrm{con}}{C}_m^a{C}_n^b}{k_r{C}_m^a{C}_n^b+C_c{r}_{\mathrm{con}}}---\left(21\right)$

根据文献，k_r与r_con成正比，则ω₈可简化为：

${\omega }_8={\mathrm{Kr}}_{\mathrm{con}}{C}_m^a{C}_n^b---\left(22\right)$

如图3所示，为描述第一阶段单个金属粒子的化学反应过程，将 反应区域分为内层、中间层和外层。

内层：该层是液化的铝颗粒(液滴)，密度和温度分布均匀，蒸 发出大量铝蒸汽，并吸收热量。

中间层：化学反应R3～R8均在此进行，并释放出大量热，作为 体系温度上升的能量来源。反应生成的最终氧化产物在中间层与内层 的边界处凝固。在反应进行的过程中，中间层的温度始终高于内层和 外层。一般认为铝蒸汽的氧化反应进行得较为迅速，中间层主要成分 是氧气。

外层：不发生任何化学反应，与外界和中间层进行能量交换。主 要成分是氧气。

定义变量如下，分别对内层、中间层和外层反应动力学方程。

ρ_Al——铝液滴密度

Re_g——雷诺数

E_Al——铝液滴表面每秒蒸发量(g/s)

ρ_g——氧气密度(g/m³)

d——液滴直径(m)

u_g——液滴表面氧气运动速度(m/s)

A——液滴表面积(m²)

u_d——铝液滴运动速度(m/s)

ρ_Al——铝液滴密度(g/m³)

k_V——氧气膨胀系数(1/K)

T₁——铝液滴温度(K)

T₃——外层氧气温度(K)

m₁——铝液滴质量(g)

C_D——阻力系数

c_Al——铝比热容(J/(K·g))

h_eva——铝蒸发热(J/g)

μ_g——氧气粘度(Pa·s)

S_Al2O3(t)——中间层反应剩余的Al₂O₃质量(g)

c_g——氧气比热容(J/(K·g))

Nu_d——努赛尔数

ε_Al——铝发射率

σ——Stefan-Boltzmann常数

P_r——普朗特数

T₂——中间层温度(K)

S_Al(t)——中间层反应剩余的Al质量(g)

1、内层

内层主要发生液滴表面铝的蒸发反应。E_Al表示为：

E_Al＝ω_AlA     (23)

其中，ω_Al是蒸发速率，可由式1～4得到；A是液滴表面积。A 是时间的函数，选取液滴初始半径为15μm，初始质量为3.342×10^-14g， t时刻A和d通过下式计算得到：

$A\left(t\right)=4\pi {\left(\frac{10.026\times {10}^{-14}-3E_{\mathrm{Al}}t}{4{\mathrm{\pi \rho }}_{\mathrm{Al}}}\right)}^{\frac{2}{3}}---\left(24\right)$

$d=2{\left(\frac{10.026\times {10}^{-14}-3E_{\mathrm{Al}}t}{4{\mathrm{\pi \rho }}_{\mathrm{Al}}}\right)}^{\frac{1}{3}}---\left(25\right)$

根据文献，铝液滴密度ρ_Al可表示为：

ρ_Al＝2364000[1-0.0002268(T₁-933)](g·m³)     (26)

内层的热量交换主要来自于与中间层氧气的对流换热，自身热辐 射散热，以及蒸发吸收的热量。氧气通常认为不辐射也不吸收热辐射， 因此，内层不接收中间层和外层的辐射热。所以，可得到能量传递方 程如下：

$\left(\begin{array}{c}{m}_1{c}_{\mathrm{Al}}\frac{{\mathrm{dT}}_1}{\mathrm{dt}}=Q_{\mathrm{conv}}-Q_{\mathrm{rad}}-Q_{\mathrm{Al},\mathrm{melt}}\\ =\mathrm{\pi d}\frac{{\mu }_g{c}_g}{{\mathrm{Pr}}_g}(T_2-T_1){\mathrm{Nu}}_d-{ϵ}_{\mathrm{Al}}\sigma T_1^4-E_{\mathrm{Al}}{h}_{\mathrm{eva}}\end{array}\right)---\left(27\right)$

其中m₁＝3.342×10^-8g-tE_Al。

努赛尔数可通过雷诺数和普朗特数表示^[2]：

${\mathrm{Nu}}_d=2[1+{{\mathrm{Re}}_d}^{\frac{1}{2}}{{\mathrm{Pr}}_g}^{\frac{1}{3}}/3]---\left(28\right)$

将铝液滴直径作为特征长度，结合液滴与氧气的相对速度，可求 得雷诺数：

Re_d＝dρ_g|u_g-u_d|/μ_g     (29)

氧气的运动速度u_g主要来自于反应核气体的整体热膨胀，铝液 滴速度u_d主要来自于气体膨胀过程中对其的推力(以地面为参考系 是推力，以液滴为参考系是阻力)。

反应核发生化学反应时，产生热量，使得反应核中氧气膨胀。由 于反应核处于高空自由空间，认为反应核周围压强不变，因此，根据 气体热膨胀定律，反应核中气体当前温度为T，初始温度为T₀、初始 体积为V₀时，气体体积(即反应核体积)为：

V＝V₀[1+k_v(T-T₀)]     (30)

由体积可得到半径，取半径的时间导数，即可得到以地面为参考 系时，氧气的运动速度u_g(取初始体积0.0335m³，初始温度930K， k_V＝0.00367)：

$u_g=0.3224\frac{d\left[\sqrt[3]{1+0.00367(T_3-930)}\right]}{\mathrm{dt}}---\left(31\right)$

铝液滴主要受到气体膨胀过程中对其的推力，根据牛顿第二定 律：

$m_1\frac{{\mathrm{du}}_d}{\mathrm{dt}}=\frac{\pi }{8}{\rho }_g{d}^2{C}_D{(u_g-u_d)}^2---\left(32\right)$

阻力系数C_D可通过下式求得：

$C_D=\left(\begin{array}{cc}24(1+{{\mathrm{Re}}_d}^{\frac{1}{2}}/6)/{\mathrm{Re}}_d& {\mathrm{Re}}_d\le 1000\\ 0.424& {\mathrm{Re}}_d>1000\end{array}\right)---\left(33\right)$

2、中间层

内层表面的铝蒸汽蒸发至中间层后，与氧气发生R3反应。根据 质量作用定律，Al的消耗速率和AlO的生成速率(g/cm³)为：

${\omega }_{3,\mathrm{Al}}=23k_3{C}_{\mathrm{Al}}{C_{O_2}}^{\frac{1}{2}},{\omega }_{3,\mathrm{AlO}}=39k_3{C}_{\mathrm{Al}}{C_{O_2}}^{\frac{1}{2}}---\left(34\right)$

其中C_Al、C_AlO为Al和AlO的摩尔浓度(mol/cm³)。

设t时刻，反应后，中间层Al剩余量为S_Al(t)，某一时刻Al剩余 量的该变量等于该时刻内层蒸发的Al减去R3反应消耗的Al，所以 S_Al(t)服从方程：

$\left(\begin{array}{c}\frac{{\mathrm{dS}}_{\mathrm{Al}}}{\mathrm{dt}}{E}_{\mathrm{Al}}-k_3{S}_{\mathrm{Al}}{\left(\frac{{\rho }_g}{32}\right)}^{\frac{1}{2}}\\ S_{\mathrm{Al}}\left(0\right)=0\end{array}\right)---\left(35\right)$

某一时刻反应R3消耗的O₂的质量m_O2(g)及产生的热量Q_R3(J) 为：

$m_{O_2}=1.391(3.342\times {10}^{-14}-m_1-S_{\mathrm{Al}})---\left(36\right)$

Q_R3＝-3956.5(3.342×10^-14-m₁-S_Al)    (37)

反应R4、R5的产物消耗速率和产物生成速率仍可通过质量作用 定律得到。同理，t时刻，反应R4、R5和R6进行后，AlO剩余量为 S_AlO(t)服从方程：

$\left(\begin{array}{c}\frac{{\mathrm{dS}}_{\mathrm{AlO}}}{\mathrm{dt}}=k_3{S}_{\mathrm{Al}}{\left(\frac{{\rho }_g}{32}\right)}^{\frac{1}{2}}-k_4{S}_{\mathrm{AlO}}{\left(\frac{{\rho }_g}{32}\right)}^{\frac{1}{2}}\\ S_{\mathrm{AlO}}\left(0\right)=0\end{array}\right)---\left(38\right)$

反应生成的AlO₂的质量m_AlO2(g)及产生的热量Q_R4(J)为：

$m_{{\mathrm{AlO}}_2}=1.41(m_{\mathrm{AlO}}-S_{\mathrm{AlO}})---\left(39\right)$

Q_R4＝-974.36(m_AlO-S_AlO)    (40)

中间层与外层的热交换形式为传导换热。通过联立以上方程组， 即可解出所需要的物理量。

2、传热计算模块

从起始时刻起，最外层会与外界发生热交换，从而导致各层出现 温度梯度，发生热传导现象；各层温度不同，会导致化学反应进行的 程度不同，对反应气体的消耗程度不同，因而各层气体浓度发生变化， 出现浓度梯度，发生扩散热交换；同时，“火球”在膨胀，对外做功， 并且对外辐射能量，导致各层温度降低。各层传热现象可用以下方程 描述：

$A_k\lambda (\frac{T_{k-1,i-1}-T_{k,i-1}}{\mathrm{\Delta x}}-\frac{T_{k,i-1}-T_{k+1,i-1}}{\mathrm{\Delta x}})-A_k(h_{k-1}{J}_{k-1}-h_k{J}_{k-1})-W_{\mathrm{rad}}-W=m_k{c}_k\frac{T_{k,i}-T_{k,i-1}}{\mathrm{\Delta t}}$

$\left(41\right)$

其中，A为该层(第k层)表面积。T_k-1，i-1为第k-1层在i-1时刻 的温度，T_k，i为第k-1层在i-1时刻的温度。h和J分别为扩散热焓和 扩散通量。W_rad和W分别为辐射和对外做功损失的能量。m_k、c_k分别 为该层气体质量和比热容。

通过公式(41)计算得到某时刻“火球”内各层温度T_k，i，由此 可以计算出该时刻“火球”内各层微粒光谱辐射能量Q_k，i(λ)。

其中h为普朗克常数，c为光速，λ为波长，k_波为波尔兹曼常数。

3、各层颗粒吸收和衰减系数计算模块

燃烧产物在高温下会电离出自由电荷，同时产生从Al或O原子 中电离出的空穴。产物连续发射谱归因于本征载流子。因此，可以通 过Drude-Lorentz模型来求解燃烧产物的复折射率虚部，进而可以得 到发射率。

类比Lorentz色散模型，自由电荷运动方程可写为：

$m\frac{{\partial }^2\overrightarrow{r}\left(t\right)}{\partial t^2}+\mathrm{m\Gamma }\frac{\partial \overrightarrow{r}\left(t\right)}{\partial t}=-e\overrightarrow{E}\left(t\right)---\left(43\right)$

方程解为：

$\overrightarrow{r}\left(t\right)=\frac{e}{m}\frac{\overrightarrow{E}\left(\omega \right)}{({\omega }^2+\mathrm{i\Gamma \omega })}---\left(44\right)$

每个电子提供的原子偶极距为P’＝-e×r，则介质宏观极化强度为：

P＝N×P’＝ε₀χE    (45)

其中，N为载流子数密度。由此可得电极化率χ：

$\chi \left(\omega \right)=-\frac{{\omega }_p^2}{{\omega }^2+\mathrm{i\Gamma \omega }}---\left(46\right)$

其中，ω_p²＝Ne²/ε₀m，ε₀为真空介电常数，m为载流子质量。进 一步可得：

${\chi }^{\prime }\left(\omega \right)=-{\omega }_p^2\frac{1}{{\omega }^2+{\Gamma }^2},{\chi }^{\prime \prime }\left(\omega \right)={\omega }_p^2\frac{\Gamma /\omega }{{\omega }^2+{\Gamma }^2}---\left(47\right)$

由1+χ＝ε得：

${{ϵ}_r}^{\prime }\left(\omega \right)=1-{\omega }_P^2\frac{1}{{\omega }^2+{\Gamma }^2},{{ϵ}_r}^{\prime \prime }\left(\omega \right)={\omega }_P^2\frac{\Gamma /\omega }{{\omega }^2+{\Gamma }^2}---\left(48\right)$

其中，Γ＝τ^-1，τ为产物粒子弛豫时间。

对于粒子复折射率，n²-k²＝ε′_r，2nk＝ε″_r，粒子复折 射率虚部可写为：

k_虚＝[(ε-ε′_r)/2]^1/2    (49)

本征载流子数密度N：

通过固体物理方法的推导，可得本征载流子数密度为：

其中，k_波为玻尔兹曼常数，m_e、m_h分别为电子和空穴有效质量， T为温度，h为约化普朗克常数，将式(50)代入(46)，再代入(48)，即 可求解出发射率。带隙能E_g写为：

E_g＝9.1eV-(0.0011eV/K)×T-E_pol    (51)

其中，E_pol＝2/3eV(eV-电子伏特，K-开尔文)。

不同温度，不同波长对应不同的k_虚值。

Al₂O₃复折射率实部k_实的计算：

k_实值可通过以下经验式进行计算，

$\left(\begin{array}{c}{\lambda }_1^2=0.00377588,A_1=1.023798\\ {\lambda }_2^2=0.0122544,A_2=1.058264\\ {\lambda }_3^2=321.3616,A_3=5.280792\end{array}\right)---\left(52\right)$

通过上式即可得到不同波长下的折射率实部k_实值。

吸收和衰减效率因子的Mie理论计算：

根据理论Mie，当光强为I₀，在颗粒周围介质中波长为λ的自然 光平行入射到一半径为γ的各向同性球形颗粒上时，在散射角为θ， 距离散射体r_散处的散射光强为：

其中

$i_1=S_1(m,\theta ,\alpha )\times S_1^{*}(m,\theta ,\alpha )$

$i_2=S_2(m,\theta ,\alpha )\times S_2^{*}(m,\theta ,\alpha )---\left(54\right)$

$S_1\left(\theta \right)=\underset{n=1}{\overset{\infty }{\Sigma }}\frac{2n+1}{n(n+1)}(a_n{\pi }_n\left(\cos \theta \right)+b_n{\tau }_n\left(\cos \theta \right))$

$S_2\left(\theta \right)=\underset{n=1}{\overset{\infty }{\Sigma }}\frac{2n+1}{n(n+1)}(a_n{\tau }_n\left(\cos \theta \right)+b_n{\pi }_n\left(\cos \theta \right))---\left(55\right)$

$\alpha =\frac{\mathrm{\pi D}}{\lambda }---\left(56\right)$

$\left(\begin{array}{c}{\pi }_m\left(\cos \theta \right)=\frac{P_m^1\left(\cos \theta \right)}{\sin \theta }=\frac{d{P}_m\left(\cos \theta \right)}{d\left(\cos \theta \right)}\\ {\tau }_m\left(\cos \theta \right)=\frac{d{P}_m^1\left(\cos \theta \right)}{\mathrm{d\theta }}\end{array}\right)---\left(57\right)$

P_m¹为m阶勒让德多项式

$a_n=\frac{\mu m^2{j}_n\left(\mathrm{m\alpha }\right){\left[\alpha j_x\left(\alpha \right)\right]}^{\prime }-{\mu }_1{j}_n\left(\alpha \right){\left[\mathrm{m\alpha }{j}_n\left(\mathrm{m\alpha }\right)\right]}^{\prime }}{\mu m^2{j}_n\left(\mathrm{m\alpha }\right){\left[\alpha h_n^{\left(1\right)}\left(\alpha \right)\right]}^{\prime }-{\mu }_1{h}_n^{\left(1\right)}\left(\alpha \right){\left[\mathrm{m\alpha }{j}_n\left(\mathrm{m\alpha }\right)\right]}^{\prime }}$

$b_n=\frac{{\mu }_1{j}_n\left(\mathrm{m\alpha }\right){\left[x{j}_x\left(\alpha \right)\right]}^{\prime }-\mu j_n\left(\alpha \right){\left[\mathrm{m\alpha }{j}_n\left(\mathrm{m\alpha }\right)\right]}^{\prime }}{{\mu }_1{j}_n\left(\mathrm{m\alpha }\right){\left[x{h}_n^{\left(1\right)}\left(\alpha \right)\right]}^{\prime }-\mu h_n^{\left(1\right)}\left(\alpha \right){\left[\mathrm{m\alpha }{j}_n\left(\mathrm{m\alpha }\right)\right]}^{\prime }}---\left(58\right)$

其中，α是尺寸参数，D为粒子球直径，m＝k_实-ik_虚为粒子球的 复折射率，j_n(x)为球贝塞尔函数，h_n⁽¹⁾(x)为第一类汉克尔函数。P_n(x) 为勒让德多项式。

由计算得到的复折射率，和化学反应模块输出的粒子半径，通过 Mie理论即可计算得到各波段吸收效率因子Q_abs和衰减效率因子Q_ext。

$Q_{\mathrm{abs}}=\frac{2}{{\alpha }^2}\underset{n=1}{\overset{\infty }{\Sigma }}(2n+1)({\left|a_n\right|}^2+{\left|b_n\right|}^2)---\left(59\right)$

$Q_{\mathrm{ext}}=\frac{2}{{\alpha }^2}\underset{n=1}{\overset{\infty }{\Sigma }}(2n+1)\mathrm{Re}(a_n+b_n)---\left(60\right)$

从而，颗粒吸收系数ε和衰减系数α可分别写为：

ε(λ)＝Q_abs(λ)

                          (61)

α(λ)＝2πr²NQ_ext(λ)

其中：N为本征载流子数密度，r为“火球”半径。

4、各层光程计算模块

该模块用于计算“火球”模型内任意一微粒所发出光线在“火球” 内各层的光程。由于“火球”的旋转对称性，可以将“火球”模型抽 象成二维模型，用一系列同心圆表示，以圆心为坐标原点可建立平面 直角坐标系。示意图如图4：

图中最内层表示爆炸所产生的真空层，该层折射率为1，故不在 计算范围内。

发光微粒可抽象为同心圆内任意位置的一点，由于圆环的对称性 可将该点设定在坐标系中的x轴上，设该点为X₀点，该点所在层数 为第k层。该点发出的光线可抽象为以该点为端点的射线。设定一个 角度步长，可以对0°-180°内任意角度进行光程计算。由圆环的轴对 称性，只计算x轴上方两象限的光线即可，x轴下方的光线可由对称 关系得到。

由于追迹的光线相对于x轴偏角不同，则可以将光程计算分为两 种情况：光线与目标点以内的圆环相交的情况和光线与目标点以内的 圆环相离的情况，用几何关系很容易表示。两种情况的示意图分别如 图5和6：

当所追迹的光线与目标点以内的圆环相离时，则有光线在目标点 以内的光程为0，只需计算光线在目标点以外的光程即可。如上图所 示，已知OX₀、OA、θ角，则可由余弦定理求出X₀A，X₀A即为该 光线在第i层内的光程。同理可以求出X₀B，二者之差即为该光线在 第k+1层内的光程。以此类推可以求出该光线在以外各层内的光程。

当光线与目标点以内的圆环相交时，除真空层内光程可以忽略不 计外，其余层内光程均需要计算。将直线方程与不同半径的圆方程联 立即可求得二者交点，由交点坐标即可求出两交点间距离。如图6所 示，BC即为真空层内光程，AD-BC为第一层内光程，以此类推，直 到第k层。第k层以后的光程可由上一种情况的计算方法计算。

该模块计算得出的光程l_k与之前计算得到的消光系数经可以计 算出“火球”内各层光谱透过率T_k(λ)。

其中：Q_abs(λ)为各波段吸收效率因子，Q_ext(λ)为各波段衰减效率 因子，α_k(λ)为对应的消光系数，l_k为微粒在第k层经过的光程。

通过T_k(λ)和之前计算得到的Q_k，j(λ)可以计算“火球”内部任意小 颗粒辐射到“火球”尺度外的光谱辐射能量，只要光线取的足够密集， 可以计算得出较为准确的光谱辐射能量值，最后通过(64)对“火球” 尺度内的颗粒进行积分即可得出致盲弹燃烧过程中某一时刻、各个波 段的辐射能量。

其中：T_k(λ)为该粒子所发出的光线的光谱透过率，N_层为所划分 的层数，N_光为粒子所发出的光线数，Q_k，i（λ）为该层粒子光谱辐射能量， Q_k外(λ)为该层粒子发出经“火球”吸收后辐射出的能量，ρ_k为各层粒 子浓度，Q_外.(λ)为“火球”对外光谱辐射能量。

通过以上过程可以计算得出致盲弹燃烧过程中某一时刻、各个波 段的辐射能量及“火球”(辐射源)的空间尺度。光谱辐射能量计算结 果如图7-8所示，仿真结果与实验结果比较，误差在和接受范围内， 仿真结果可靠，能反应强光致盲弹燃烧过程的辐射特性。

另外，该仿真方法具有普遍性，可以应用到任意高温稀疏体系辐 射能量仿真中。