[摘要]模仿复杂化学反应的数值法，不仅正被人们用来提供一种根据基本物理和化学过程的知识预言化学体系行为的能力，而且被用来深入了解这些体系的反应机制。本文不仅评述了这方面的科技水平，而且预测了将来很有可能的发展。
化学反应不仅是造成地球上生命的决定性因素，而且已经被人们所广泛开发利用，对于我们现今高度发达的文明来讲，它们是必不可少的。许多世纪的实验和新近的理论研究已使我们不仅积累下有关它们的丰富知识，而且使我们能实际应用这些过程。因此，如果我们仍然仅停留在用一般的术语来理解这些复杂体系中的某些体系那就显得有些不相称了。其原因是不难发现的：在诸如燃烧、大气干扰、生命循环等现象中同时存在的大量物理和化学过程，甚至还有化学工业操作都难以定量描述的整个过程，尽管我们可能很详细地了解单个组元的行为。
数学是一种定量语言，任何定量描述都必须用它来写出；可是分析方法解决这样复杂问题的能力已被远远超出。随着计算机容量变大，速度更快，利用它们来解决这些问题的可能性已经增长，而且计算机模拟已进入到物理、生物、社会科学、工程和商业界的许多领域。各项的模型建立和模拟包含着广泛的含义，其范围从测量结果简单地拟合为某一数学函数（但是是任意的）到从因果关系为基础的描述，这种描述还结合着对高度复杂系统多个坦元的基本原理的描述。化学动力学的模拟技术在过去十年里成熟起来了。它们作为一种相关和预测工具的应用已经扩大，而且实现其深入了解反应机制的能力在继续发展之中。
质量作用动力学；均匀系统
化学动力学的正规数学理论是在十九世纪提出的，并且以质量作用定律的阐述达到顶峰。简单地说，质量作用定律给出了基本反应
式中nA是物类A的分子数，方括号表示浓度，速度常数K被认为是该反应的一个基本性质。该反应体系的时间演变可由这类ODE（ODE指常微分方程式，下同）的解给出。这些考虑因素曾经成功地解释了有关反应速率的许多数据，同时也刺激在广泛的化学体系中展开实验。不可避免地发现了一些并不符合这些简单速率定律的情况；这些情况随着反应机制概念的出现逐渐得到阐明，反应机制指出同时或连续地发生着许多反应。这样一种机制的每一组元都被认为是服从方程式（2）给出的它自己的速率定律，而复杂系统的行为是所有过程同时起作用的一种结果。应用质量作用原理给出的联立微分方程组仅能对几种比较简单的情况用解析法求解处理比较复杂情况的一些方法发展起来了，它们是层层采取归并、简化和近似的方法。结果，导出了大量化学机制的文献，虽然这些机制看起来表面上是有道理的，但在细节上，它常常比之真正的化学过程更是随可解数学的选择而定的。
第二次世界大战之后，人们开始可以得到高速计算机，这使许多工作者研究将其用来解质量作用微分方程的可能性，而运用解析法则是不可能做到的。这些努力获得了不同程度的成功，这是因为虽然在那时已经很好地发展了ode的数值解，但是化学速率方程呈现出一些困难，这些方程似乎与处理其它问题得到的经验相矛盾的。自动的ODE解算器因为调整作为变化速率这个增量而逐步及时地推进了上述这问题，并且变动解的曲率，结果保证了一个预定的准确度。计算一个化学反应体系的预期行为是：较短时间阶跃总是用在浓度变化很快的开头段，而随着反应变慢下来并接近平衡时间阶跃增大。但经验证明却正好相反时间阶跃在浓度变化是极小（图1）的稳态条件下去得非常之小。
这种表面上反常行为的起因首先是Curtiss和Hirschfelder解释的，他们说这是由于在由相反反应平衡造成的稳态下的化学系统中固有着很大的负反馈。对这种条件的任何稍微的偏移，随之都将会很快地回复到稳态。为了使解保持稳定，也就是防止误差一步步地传播和增长，相应方程数值积分的时间阶跃必须限定在与这种快速复原力相当的值。与机械伺服机构相类似，人们把显示这种性质的微分方程体系统称为“刚性的”。
求解显示出特征刚性的动力学问题需要大量的时间，它使得在五十年代和六十年代比较慢的机器上处理它们不切实际，因而这个领域曾经一度失去生气。也曾有过偶尔的应用，特别是在化学工程方面，在这个领域，是用把在该组中起决定作用的那些微分方程加以稳态表述的权宜办法来克服这种困难的。如果所得到的答案不是原来问题的正确解，在许多方程中那也是没有什么关系的，因为工程场合的目标通常是找到一条模拟工艺过程数据的途径，以便根据实验测量结果进行插入和外推。只要得出模型参数的实验程序与随后连续模拟中的相同，数学解是错的也无关紧要，对于这个问题来说，模型本身是否精确甚至也是不重要的。但是，当在基本动力学研究中尝试使用这相同的途径时，速度常数分配的矛盾开始显现了，不幸的是，直到可以得到更好的数值解，有时是若干年后，人们才识到许多矛盾是这些误差的结果。
在这期间，随着为了环境和其他目的而模拟大型大气系统的需要的增长，对这些刚性微分方程的数学兴趣得到增强。这些求解方法的稳定特征主要是通过Dahlquist的工作而变得清楚了，而且事情表明，需要有一种隐式的方法，它包含在每一时间阶跃时一组联立代数方程在偏导数（雅各比矩阵）条件下的解。Gear研究了准确性和稳定性之间的权衡，他曾发表过以一种经济的方法进行这种计算的一种实用代码（图1）。自此之后，这种技术作了改进，并想出了其他一些方法。Warner评述了这个领域的进展，而最近Byrne评价了一些可用代码。
建立化学动力学过程模型另有一种可供选用的方法，它考虑了在每一时间间隔里发生不连续反应事件的几率。这种“随机”途径已经表明，它基本上相当于求解连续体系微分方程的“确定”技术。所需的蒙特——卡罗数字计算当然较之ODE解执行起来简单得多，但它收敛较慢，因而对大体系变得昂贵起来。由于这个原因，它们在反应模拟中几乎不大使用。
最早应用模拟来研究化学动力学的工作之一是由Allara和Edelson做出的，他们在复杂机制的研究中引进了一种新的基本原理。他们研究过的工艺过程是恒容恒温热解低链烯烃。这些体系已从实验上研究多年，而且得出大批有关反应速度和产物分布的资料。图2给出的是一个观察其行为的典型例子，其中丙烷的分解速率是作为时间的一个函数示出的，它是根据出现产物之一甲烷来监测的。大多数烷烃的一个共同特性是在反应之初热解速率开始降低。曾经提出过关于这种自抑制效应的许多解释。虽然一些解释包含着具有多达15或20个反应的机制，但它们只是被近似地加以处理，而且开展过有关这种效应的真正起因的争论。似非而可能是的是，已经积累了有关基本自由基反应速率的大量丰富资料，真正的反应机制必须根据这些反应构成。Allara和Edelson提出，要利用所有的一切这种信息，并利用数字方法，以获得大量产生的ODE的解。他们模拟过机制为几百个化学反应的工艺过程。尽管这些反应中有许多的速率常数不过是粗糙的估计，但它们还是得到了与实验值十分相符的结果；图2中约直线是计算结果。产物的分布也表明模仿得不错。
这种通过计算机模拟解决复杂机制动力学的途径，最初人们是以相当大怀疑态度欢迎它的。例如，期待图2中相当简单的曲线来反映包含在数百个速率常数中的信息，人们认为这是愚蠢可笑的。这就误解了这些结果的重要意义，因为推出这些速率常数的独立实验证据的广泛背景也必须加以考虑。适当的结论是，产生这些结果的机制是与这个信息总体一致的。作者们导出过机制行为所处的一个范围，导出的基础是单个反应速率的误差估计。但是，声称提出的机制是极好的是不可能。对基团中间产物计算浓度和单个反应通量的检验意味着，自抑制效应的一个可能起因是在后期反应性较差的不饱和基团的重要性在增长。这个结论虽然是定性的，并且带有点主观，但它是最早的征兆之一，它表明模仿途径可以用来深入了解一个复杂反应机制中的控制途径。
尽管化学界的大部分人一开始不愿意承认模拟是一种有效的解析复杂机制问题的途径，在其后不几年，人们就看到对了解反应体系有所贡献的领域数在逐步增加。1977年3月举行的一次学术讨论会，回顾了范围广泛如电化学，燃烧，大气化学和生物化学等领域互异的成功事例。一次特别有趣的应用是用到振荡反应中，在这种情况下模拟用于阐明复杂的Belousov-Zhabotinsky反应补充了分叉理论在这些振荡中的纯数学应用。它表明，在纯化学体系中振荡是有可能的，此外它还表明，向振荡态的跃迁是质量作用动力学的结果，因而它并不要求在一个多元稳态系统中出现波动。
非均匀系统；偏微分方程
从早期的国防应用还发展出把作为一种预言工具的模拟用于大规模的体系，这些应用与上层大气中核武器效应和导弹返回大气层有关。大气环境引起人们巨大的关切，因为人们关心臭氧层中可能发生的变化以及气候的改变，引起这些变化的因素有超声运输的喷出物或氟化碳气溶胶的弥散。随着有关适宜控制策略的争论盛行，机动车辆的喷出物及其通过产生氧化氮、臭氧、刺激物和烟雾等对空气质量的影响，曾是模型制作者的首要目标。与此有关的是燃烧研究，因为发动机和燃烧炉是主要污染物源。燃料工业经济的快速变化，并且有必要在环境影响最小的前提下最有效的提取能量进一步推动了这方面的研究。在这一切努力中，曾使用模拟来使决定这个问题的许多物理和化学过程的巨大基本信息库归并成一个复杂体系的综合描述。它们的共同目标是预言在一些条件下的体系行为，在这些条件下，直接的实验往往昂贵无比，根本不切实际，或者兼而有之（或者如在核效应的情况中那样是由政治考虑而排除的）。
人们将会立即认识到，这些真实体系（与实验室实验的理想恒定容积充分混合的反应器相反）具有一种不同阶次的数学复杂性，因为它们不仅包含化学物类浓度的时间依赖性，而且有空间依赖性，而（像热传递和质量传递等）物理过程也必须像化学反应一样包括在内。这些过程本身常常是很复杂的，以致模型制作者愿意在全然不考虑化学过程的情况下解决它们。有时可以把化学过程与体系隔开并可单独分开解决，如像在大气问题中那样，在这样的情况下传递是靠外部能量输入（例如阳光）驱动的，化学过程常常不过是一种比较轻微的扰动而已。最复杂的体系也许是燃料，在燃烧的情况下，正是化学过程提供了驱动传递过程的能量，而传递过程又影响着化学反应的速率。为了使这些问题易于处理，化学过程的细节常常被缩简为最小一组反应，缩简的方法是“归并”——也就是把相同类型的反应合并成一个反应——以及规定总的物类浓度有一平均或总的速率常数。在极端情况下，所有化学过程可以以一个方程式表达：
燃料+氧化剂→产物+热 （3）
显然，化学过程的细节失去了。对于燃烧中的许多问题，例如稳态炉锅的模型建立，这种途径也许足以描述该体系的温度和流动特征。其中希望的描述涉及到化学细节的其他问题，例如点火或污染，不可能以这种方式来处理。
在数学上，这些问题包括一类空间和时间的偏微分方程（PDE表示偏微分方程，下同），这与充分搅拌的反应器是大不相同的，充分搅拌的反应器仅是时间的一个常微分方程。解决这些问题的计算技术远不及ODE的技术那么先进。被广泛应用的一大类方法被称之为线段法，进一步细分为有限差和有限元技术。这种空间领域细分为更小区域的格子或“网孔”。在有限差程序中，空间依赖性在各子区间中接近一个常数值，结果以一种分步表示代替了连续函数。把空间区域划分得足够精细是需要的，以便分步表示法把解接近于希望的准确度。写出各区间的一组动力学微分方程，在它们上边增添代表邻近单元间的质量（和能量）流动的一些项。因此，模拟就从一个偏微分方程问题转化成了一个常微分方程问题，这些方程中有着许许多多数目等于网孔单元数的因变量。上述刚性的ODE技术可以用来求解这个问题。但是，既然所需的计算机磁芯存贮库容粗略地以变量数的平方形式出现，而且介于平方和立方之间的这种运算增加，所以完成这些计算的费用增长得很快。
有一种替代程序，即有限元技术，是用一组基础函数，一般是多项仿样函数代表着每个网孔单元的解，这些仿样函数不得不通过网孔边界连续接合起来。应用Rayleigh-Ritz-Galerkin最小化技术来寻找最佳拟合把问题转化成一组多项系数的常微分方程（ODE'S）。所需的网孔单元比相当的有限差方法的要少，但每个变量有更多的系数；对于有轻度空间依赖性的一个性能良好的问题，两种方法计算费用已表明是大体相同的。
当这些空间相关问题变得越来越大时，保持最经济使用计算机设备的问题就变得比较严重起来。如像火焰和震荡这样一些类型问题，情况尤其是这样，在这类问题中，大多数化学过程发生在一个很狭窄的区域里，因此模拟显示出浓度梯度十分陡峭。因此，人们有必要利用一个精细网孔来相当精确地表示其解，但是由于经济的原因，人们不希望在陡峭梯度之外再继续作这种细分（图3）。如果陡峭梯度在运动中，则人们可以选择解决这个问题，要么一个运动参照系中的一个固定格子上解，要么用一个固定坐标系统中的一个运动中的格子解（或两者的某种结合）。对于自动适应可变网孔方法的应用，有限元技术特别适宜。最近发展出的一种方法试图尽力缩小解的多项表示法中的误差，其途径是调节网孔边界的位置（图4）。这些网孔坐标于是在一组辅助微分方程中成为可变的，这组方程与主问题同时求解。另一途径是把真实问题的一个可变网孔变换成一个固定网孔，该问题再在此网孔中求解，同时要这种变换不得不遵循解的预期特征（图5）。这后一技术特别适宜于具有运动的断续函数，例如相边界体系，这些相边界的位置必须遵循而且位于一个网孔边界上。由于这类的几个问题具有重要意义，人们正在积极发展其他的求解途径。
在模拟这些大型问题中的计算机时间，大部分花在解出隐式微分方程求解技术要求的线性代数系统上。在这些系统中，遇到的矩阵常常是稀疏的——也就是它们具有不多的非零项，而且也可能具有一定的结构。例如，一维流动和化学问题的有限差表述具有是块状三对角的雅各比矩阵（图6）。有些专门的程序可用来处理这种“带状”问题，它们使人们没有必要带有大量零元素的计算机计算或存贮大量零元素，因而造成相当大的节约。尽管没有特殊结构，代表化学相互作用的矩阵本身可能是稀疏的（图7）。可以应用重新安排技术来把这些矩阵变换成具有可以用来获益的结构矩阵，在其它情况下，一般稀疏矩阵法也可能适用。这些技术的应用仍然尚处于其初期阶段，但人们可以预期它们是会增长的。二维或三维空间中的问题相应地有着更大的矩阵要处理，因此使计算费用迅速增加。在这些情况下，结构变得特别复杂，且目前根本没有什么专门的方法处理它们。模拟这些较高维数问题的科技水平的提高将与发展求解这些矩阵方程的有效途径密切相联。
灵敏度分析
利用数学来从因果律模型得到物理或化学体系的行为得到这样一种描述，其中的变量通过含有代表模型各个不同组成部分参数的函数彼此发生关系。如果这一描述能够得到分析形式，则行为对这些参数的依赖性一目了然，而且计算其中各个变化的影响或了解模型各个部分的相对意义是一个简单的问题。但是，当必须用计算机获得这些解时，输入和输出之间的这种关系就失去了，而且模型越复杂，导出造成某些模拟特征的部分机制或了解变化将如何影响这种行为就越困难。当然，对参数中的许多变量作重复计算是可能的。对于在r个反应中涉及到的n个物类的纯化学机制，这往往需要对各个速率参数的一个变值依次求解一组n个ODE，也就是（r+1）次。对于大的机制，例如图2的热解体系，当n可能是50左右，而r可能是几百时，这种“强力”式途径可能变得冗长而又很昂贵。比较正规的数学处理这种“灵敏性分析”问题得到一些可能在理智上比较令人满意的结果，但其大多数仍然和要获得的工作一样。但是，有条途径是使用格林（Green's）函数变换，它把问题的主要部分简化为求解n组（而不是r组）n个ODE。此外，求解各组给出一种物类相对于所有速率常数的灵敏性（而不是所有物类相当于一个速率常数，而后一种情况是其他方法的特点）。认识到在一复杂机制中只能观察到极其少的物类，人们显然可见，求解极其少的方程组通常将足以给出一切有用的灵敏性信息。这已导致在实际可用的计算范围内进行广泛的灵敏性分析。
从灵敏性分析可以搜集到比改变参数效果的正式结果多得多的资料。分析前述大批热解机制生动地说明了这一点。对于广泛范围的反应条件，所有反应相对灵敏性的顺序足以用来找到灵敏的组元和应该受到进一步详细研究的速率常数。作为一种必然结果，最不重要的反应可能被除去。尽管必须避免认为这些项是同义的趋向，但在某种意义上讲，灵敏性是反应重要性的一种量度。尽管如此，在某些情况下可以应用这种概念，例如，热解研究，在这种情况下，为了减少重要的某事被省去的可能性，其中包括进了比必需多得多的反应。因而灵敏性信息为选择最小一组反应提供了一个客观的准则，当计算机设备非常需要时，如果化学过程和输送问题相联，这种选择就成了决定性因素。
灵敏性分析的最重要结果可能是增加了对机制本身作用的了解。图8表示的是热解机制中最重要反应的灵敏性对时间的依赖性。某些反应的灵敏性在其他反应灵敏性降低的同时及时地增加，可是这种行为与反应类型之间没有关系。但是，比较深入的摸索表明、及时增加的那些反应灵敏性全都涉及到不饱和基团，而降低的那些则不是这样。这提供了一个客观的证据，证明不饱和基团在热解过程的自抑制中起了作用，这一点人们已经作了假定，但是不可能得到定量的证据。
有关机制的其他重要结论是可以取得的。造成体系各个特征的机制组成部分是直接模拟物类浓度的复杂函数，它们也许是可以识别的。例如，分析浓度 - 时间依赖关系的重复发现了控制振荡化学体系周期的反应。参数 - 参数相互影响也是可以估计的，它们可以回答这样一个问题，即已知参数的误差是如何影响用机制来确定未知参数的能力的。与此类似，浓度测量误差对这些数值确定的影响也是可以估计的。最富探索性的问题是任何重要事情是否已被除去，可以一种受限制的方式加以研究。完整的技术正在快速发展之中，无疑它将会成为反应模拟的一个很重要的方面。
未来的展望
影响所有这些努力的一个共同因素就是可计算性。在早先的讨论中，模拟技术的进展是与发展求解微分方程的适当算法相联的。这种注重软件的情况，肯定一直是这些进展中的主要推动力，但是硬件的发展成果绝对不可忽视。过去的十年，从大容量计算机机组到专用小型和中型计算机，而最近又回复到具有矢量处理能力的大容量超型计算机的来回摆动中已经看到，对于这些大容量计算——边界程序喜欢执行。数字数学已经进展到这样的程度，即在原则上可以令人满意地模仿很复杂的系统；剩下的问题是多少设备和费用可被提交给某一特定问题来用。计算费用的主要部分是用于求解为隐式积分法所要求的线性方程。为了特别好地完成这类计算而设计了一些超型计算机，而且早已了解到大容量计算，例如灵敏性分析，实用性的主要变化。随着这些新计算机结构变得普遍起来，而且通过软件的专门发展和最优化充分利用其能力，今后不多几年有希望看到在模拟应用于反应体系方面将有极大的进展。
[Science 27，1981年11月214卷4524期]