化学动力学是物理化学^[1]、化学动力学^[2]、反应工程、无机化学^[3]等课程的主要内容之一。直观动态地显示各类反应系统反应组分的浓度随时间的关系及温度的影响，有助于了解和掌握相关反应的动力学规律，这通常需要求解反应动力学微分方程组。然而，这类方程组一般无函数形式的通解，因此需借助数值积分方法得到数值解。反应动力学微分方程形式多样，编写能应用于各种情形的通用软件亦比较困难。到目前为止，国内已研制了多个适用于特定情形的软件^[4–8]，国外亦研制了多个软件^{[9, 10]}，但这些软件很少涉及温度的影响。我们经过多年在物理化学和化学动力学课程方面的教学研究基础上^{[11, 12]}，编写出适用于各类反应系统的动力学模拟计算通用软件，可以研究不同温度下各组分浓度随时间的关系，有助于学生在教学过程中了解时间和温度对反应的影响规律。

设反应系统有m个组分，分别记为A_i，i = 1, 2, …, m，存在如下n个反应：

(1) $\sum\limits_{i}{{{(-{{v}_{i}}{{{\rm{A}}}_{i}})}_{s}}}\underset{{{k}_{-, s}}}{\overset{{{k}_{+, s}}}{\longleftrightarrow}}\sum\limits_{j}{{{({{v}_{j}}{{{\rm{A}}}_{j}})}_{s}}}, \;\;s=1, 2, \cdots , n$

式中(v_i)_s为第s个反应中组分A_i的化学计量数，k_{+, s}和k_{–, s}分别为第s个反应正向和逆向的速率常数。速率常数k与温度T的关系常用推广的阿伦尼乌斯方程描述：

(2) $k=A{{T}^{a}}\exp \left( -E{\rm{a}}/RT \right)$

其中A、a、E_a、R分别为指前因子、温度系数、活化能和摩尔气体常数。相应地，恒温时反应系统的动力学微分方程可写为：

(3) $\text{d}{{c}_{{{\text{A}}_{j}}}}/\text{d}t=\sum\limits_{s\text{=1}}^{n}{{{(v_{j}^{+})}_{s}}}{{k}_{+, s}}\prod\limits_{i}{{{(c_{{{\text{A}}_{i}}}^{n_{i}^{+}})}_{s}}}\text{+}\sum\limits_{s\text{=1}}^{n}{{{(v_{j}^{-})}_{s}}}{{k}_{-, s}}\prod\limits_{i}{{{(c_{{{\text{A}}_{i}}}^{n_{i}^{-}})}_{s}}}=\\\ {{f}_{j}}(t\text{, }\ {{c}_{{{\text{A}}_{1}}}}\text{, }\ {{c}_{{{\text{A}}_{2}}}}\text{, }\ \cdots \text{ , }\ {{c}_{{{\text{A}}_{m}}}}), \ j=1, 2, \cdots , m$

对于第s个反应，式中$ {(v_j^ + )_s} $和$ {(v_j^ - )_s} $分别为正向、逆向中组分A_j的化学计量数；$ n_i^ + $和$ n_i^ - $对应为正向、逆向中组分A_i的分级数，为0时相当于该组分不参与反应。

式(3)移项积分可得：

(4) ${c}_{{{\rm{A}}}_{j}}(t)={c}_{{{\rm{A}}}_{j}}(t=0)+{\displaystyle {\int }_{t{\rm{=}}0}^{t}{\rm{ }}{f}_{j}{\rm{(}}t，{c}_{{{\rm{A}}}_{1}}，{c}_{{{\rm{A}}}_{2}}，\cdots ，{c}_{{{\rm{A}}}_{m}}{\rm{)}}}{\rm{ d}}t, \;j=1, 2, \cdots , m$

若知道t = 0时各组分的初浓度，即：

(5) ${c_{{{{\rm{A}}}_j}}}{{\rm{(}}}t = {{\rm{0)}}} = {c_{{{{\rm{A}}}_j}}}({{\rm{0}}}) , \;j=1, 2, \cdots , m$

便可知道以后任意时刻(t)任意组分的浓度。但式(4)难以通过直接积分的方法得到各组分浓度与时间的函数关系，常用4阶Runge-Kunta数值积分法进行求解^{[11, 12]}。设计算的最长反应时间为t_max，将反应时间从0–t_max分成p个等分，则每个等分(即积分步长) h为：

(6) $h={{t}_{\max }}/p$

于是t_k+1 = (k + 1)h时各组分的浓度为：

(7) ${c_{{{{\rm{A}}}_j}}}{{\rm{(}}}{t_{k + 1}}) = {c_{{{{\rm{A}}}_j}}}{{\rm{(}}}{t_k}) + ({f_{j1}} + 2{f_{j2}} + 2{f_{j3}} + {f_{j4}})/6 , \;j=1, 2, \cdots , m$

$ {f}_{j{\rm{1}}}=h{f}_{j}({t}_{k}，{c}_{{{\rm{A}}}_{1}}({t}_{k})，{c}_{{{\rm{A}}}_{2}}({t}_{k})，\cdots ，{c}_{{{\rm{A}}}_{n}}({t}_{k})) $

$ {f}_{j{\rm{2}}}=h{f}_{j}({t}_{k}+h/2，{c}_{{{\rm{A}}}_{1}}({t}_{k})+{f}_{j1}/2，{c}_{{{\rm{A}}}_{2}}({t}_{k})+{f}_{j1}/2，\cdots ，{c}_{{{\rm{A}}}_{n}}({t}_{k})+{f}_{j1}/2) $

$ {f}_{j{\rm{3}}}=h{f}_{j}({t}_{k}+h/2，{c}_{{{\rm{A}}}_{1}}({t}_{k})+{f}_{j2}/2，{c}_{{{\rm{A}}}_{2}}({t}_{k})+{f}_{j2}/2，\cdots ，{c}_{{{\rm{A}}}_{n}}({t}_{k})+{f}_{j2}/2) $

$ {f}_{j{\rm{4}}}=h{f}_{j}({t}_{k}+h，{c}_{{{\rm{A}}}_{1}}({t}_{k})+{f}_{j3}，{c}_{{{\rm{A}}}_{2}}({t}_{k})+{f}_{j3}，\cdots ，{c}_{{{\rm{A}}}_{n}}({t}_{k})+{f}_{j3})$

利用k = 0即初时刻t₀ = 0对应的初浓度(5)，使用式(7)可计算以后任意时刻t_k (k = 1, 2, …, p)的各组分浓度。

微分方程组的数值求解是利用不同时刻的浓度值近似描述一系列连续的浓度–时间函数。显然，时间等分数越多，求解越准确，但计算量越大，因此可通过逐渐增大时间等分数即p值的方式找到合理值。考虑到四阶Runge-Kutta法具有较高的计算精度，p值只要选取10000以上使计算所得的浓度合理时即可达到精度要求。

本软件使用Microsoft Visual Basic 6.0语言在Windows系统下研制，已转成可执行程序。使用前需安装和注册VB插件msflxgrd.ocx。运行软件后根据提示单击“任意动力学方程”按钮(其中“一级反应网络”功能可看参考^[12])，即可出现如图 1所示的界面，分为功能按钮、计算参数设置、反应参数输入和绘图四大区域。

实际使用时，在反应参数输入区逐行输入实际的反应式和各组分的分级数，正、逆反应的速率常数或指前因子、活化能及温度系数，一行一个反应。反应组分可用实际的分子式或用非数值开始的代号，若是数值开始的分子式如4-甲基苯酚，必须加括号即(4-甲基苯酚)或用A1等符号表示组分，否则数值4会被认为是化学计量数。若组分的分级数不是1时，则在其组分符后面加箭头号“^”和级数值，如图 1中2A^2表示组分A的分级数为2，化学计量数为−2。反应符指定为英文格式的等号“=”。反应式与速率常数的各数值之间使用英文格式的分号“; ”区分。所有反应均认为正向和逆向可同时进行，因此处理实际机理时反应式可适当减少。按正向和逆向指前因子或速率常数、正向和逆向活化能、正向与逆向温度系数的次序输入速率常数与温度的关系。若活化能和温度系数没有输入或为0时，相当于只输入速率常数，此时不考虑温度的影响。速率常数为0时表示不存在该方向的反应。反应组分的符号和数量根据输入的反应式自动确定。

反应式输入完后再输入各组分的初浓度，采用组分符号与初浓度数值加英文格式的冒号“: ”分隔，可一行一个组分或一行多个组分的方式输入，一行多个初浓度值时用英文格式的分号“; ”分隔，初浓度为0的值可不用输入。亦可输入不同时刻多种组分实际浓度的实验或准确值，此时先输入符号t加冒号及不同时间值，各值间用分号分隔，整体作为一行；接着一行是组分符加冒号及不同时间对应的浓度值，浓度值之间亦使用分号。这时可直接显示出不同时刻的实际浓度，以便与计算结果进行比较，图 1即属此情形。

输入完反应式和初浓度后，在计算参数区修改反应温度、浓度单位和绘图浓度范围、反应时间和时间单位，然后单击“绘制浓度时间曲线”按钮，即可得到各组分的浓度与时间的关系曲线。计算得到的数据可以以纯文本格式的“Txt”文件保存，以便使用其他软件如Origin或Excel等进一步处理，亦可读取纯文本格式的数据进行重新绘图。

需要注意的是，输入的数值可以是负数或带e格式(即10的幂次)的数，但不能是计算式或带有括号。为便于阅读理解，反应参数输入区可输入空格和说明行，说明行是不带有英文格式的等号和冒号的内容。输入区可进行常规的编辑操作，亦可用复制或粘贴命令。计算时，反应时间取值较大时会使微分不能用数值差分近似而异常退出，建议以最小的反应速率常数值的倒数作为反应时间选取的参考，从较小值开始逐渐增大计算。当然，亦可增大时间步长值，如从一万到几万甚至几十万，但这样会显著增加计算时间。通过与实验值比较或计算的浓度是否符合质量守恒定律等判断计算结果是否合理。

程序设定单个方向反应的最大组分数为10个，最多可研究同时存在500个反应的情形。

在介绍各种级数的反应动力学特征时，不同级数对反应完全时间的影响规律是一个重要知识点。可以考虑同时存在单组分级数n为0、1、1.5、2、3级和二组分2级的如下反应：

A0$ \xrightarrow{{{k_0}}} $B0，0级；A1$ \xrightarrow{{{k_1}}} $B1，1级；A15$ \xrightarrow{{{k_{1.5}}}} $B15，1.5级；A2$ \xrightarrow{{{k_2}}} $B2，2级；

A3$ \xrightarrow{{{k_3}}} $B3，3级；A21 + A22$ \xrightarrow{{{k_{11}}}} $B21，反应组分均为1级

设反应正向的速率常数均为1 (mol·L⁻¹)¹⁻ⁿ·s⁻¹，逆向速率常数均为0，即不存在逆反应。若各反应的反应物的初浓度均为2 mol·L⁻¹，具体的反应参数和计算参数选取输入及计算结果如图 2所示。从计算得到的曲线可知，产物B0、B1、B15、B2 (= B21)和B3的浓度曲线对应图中从右上方到中间依序的5条曲线，而反应物A0、A1、A15、A2 (= A21 = A22)和A3的浓度曲线则对应图中从右下方到中间依序的5条曲线。此时2个二级反应的3个反应物及2个产物B2和B21的浓度保持相等。但与t < 1.3 s之前的曲线次序刚好相反，原因是此时的反应物浓度数值较高且主要大于1，级数越大反应速率越大，所以级数大的产物浓度增加快而处于上方；然而，当浓度数值小于1时，级数越大反应速率越小，级数大的产物浓度增长慢反而处于下方。从图中还可知道，当反应级数大于或等于1时，反应物的浓度趋于0对应的时间∞即为反应完全时间。对级数为0，反应完全时间2 s。对于小于1的非整数或负级数反应，单独绘制所得的曲线趋势与0级类似，一定时间(即完全反应时间)后反应物浓度由正变为零，然后变为负。这时计算机自动取零处理，与实际情况对应。总的结论：其他条件相同时，级数越小，反应完全时间越小，级数等于或大于1的反应完全时间趋于无穷大。这个结论亦可从另一角度得到证明。对单一反应物A的n级单向反应：A$ \xrightarrow{{{k_n}}} $B，动力学方程为：dc_A/dt = −k_nc_Aⁿ，移项积分结果为：c_A¹⁻ⁿ = c_A0¹⁻ⁿ + (n − 1)k_nt。c_A = 0时t = (c_A0¹⁻ⁿ − 0¹⁻ⁿ)/[(1 − n)k_n]，即为反应完全时间。显然，n ≥ 1时，t = ∞；n < 1时t = c_A0¹⁻ⁿ/[(1 − n)k_n]，c_A0和k_n固定时，t随n减小而变小。另外，此时计算的浓度与准确的浓度值最大误差为±5 × 10⁻⁷ mol·L⁻¹，完全达到精度要求。

放热的对行反应有最佳温度，即转化率一定时反应速率达到最大对应的温度，在物理化学教材^[13]中常以反应2SO₂(g) + O₂(g) = 2SO₃(g)作为实例，但无实质的相关数据说明，学生难以理解。该反应正、逆方向的指前因子及活化能^[14]分别为1.63 × 10⁵ L²·mol⁻²·s⁻¹、7.20 × 10⁹ L·mol⁻¹·s⁻¹、75.32 kJ·mol⁻¹和169.0 kJ·mol⁻¹，设SO₂(g)和O₂(g)的初浓度均为0.1 mol·L⁻¹，反应参数区输入如下数据：

2SO2^2 + O2 = 2SO3^2；1.63e5；7.20e9；75320；169000

SO2: 0.1；O2: 0.1

温度分别取500、600、700和800 K时计算得到三个组分的浓度与时间关系如图 3所示。温度在500–700 K之间时转化率均可达到75%，但需要的时间分别为95000，4700和750 s，而800 K时的平衡转化率仅为59%。显然，转化率为75%时的较佳温度为700 K。这里的计算虽然难以得到准确的最佳温度，但通过分析亦可得到大致的最佳温度，使学生对产生反应最佳温度的原因得到较好的了解。

H₂和Cl₂的热反应是介绍链反应动力学常用的例子，为何自由原子H和Cl的浓度可用稳态法处理，是讨论的重点。此问题可通过具体的计算得以说明。其机理相关的反应式和相应正、逆方向的指前因子A₊、A₋及活化能E_{a, +}、E_{a, −}依次取值如下^[15]：

Cl₂ = 2Cl；8.51 × 10¹⁵ L·mol⁻¹·s⁻¹，1.25 × 10¹⁵ L·mol⁻¹·s⁻¹；234 kJ·mol⁻¹，6.820 kJ·mol⁻¹

Cl + H₂ = HCl + H；2.65 × 10¹³ L·mol⁻¹·s⁻¹，1.69 × 10¹³ L·mol⁻¹·s⁻¹；21.37 kJ·mol⁻¹，17.29 kJ·mol⁻¹

H + Cl₂ = HCl + Cl；4.822 × 10¹³ L·mol⁻¹·s⁻¹，2.84 × 10¹⁴ L·mol⁻¹·s⁻¹；3.46 kJ·mol⁻¹，198 kJ·mol⁻¹

假设H₂和Cl₂的初浓度均为10⁻⁵ mol·L⁻¹，温度为1500 K，所得的计算结果如图 4所示。从图 4中可知，反应时间为10⁻⁶ s时，H₂的转化率已达到80%以上。反应过程中H和Cl浓度最高约为0.009 × 10⁻⁵和0.068 × 10⁻⁵ mol·L⁻¹，是H₂初浓度的6.8%以下，即维持很小值。当温度降低到1000 K时反应(积分步长需选取200000，图略)，反应时间为10⁻³ s时H₂的转化率为94.7%，H和Cl的浓度最高分别为H₂初浓度的0.004%和0.064%；但当反应温度升至2000 K时(积分步长选取20000即可，图略)，反应时间为10⁻⁶ s时H₂的转化率达92.6%，H和Cl的浓度最高分别达到H₂初浓度的60%和2%。显然，温度低于1500 K时反应过程H和Cl浓度能维持较低值，所以可用稳态法处理；但温度接近2000 K时这两种中间物的浓度较高，稳态法不再能够使用。

软件还提供了复合二级反应、复合一级反应、指前因子及活化能改变对反应时间的影响，稳态近似法与平衡态近似法比较、振荡反应，酶催化和其他催化反应等许多实例，教学时可根据实际要求随时调用。亦可自行输入一些实例，丰富教学内容。

本反应动力学模拟软件具有输入方便、通用性强、界面友好、操作方便、计算快速(一般不超过10 s)等特点，同时提供了许多用于物理化学和化学动力学等相关课程的教学实例。教学实践表明，本软件提供的反应动力学实例直观明了，有助于学生掌握课程的知识点和了解动力学的原理和规律。另外，本软件亦可作为科学研究参考之用。