化学动力学方程的数值解法

1、增刊化学世界105化学动力学方程的数值解法陆兆仁(华东大学,上海200051)在化学动力学的领域中,动力学方程的求解并非是简单的事。通常提出一个反应机理需要用实验来验证,时间的关系;,至没有解析解,法来作近似处理,使用是有条件的1;,反应机理各种各样,求微分方程组的解往往不仅是很麻烦的事,而且需很多数学上的技巧,以制备胆甾醇异硫脲的甲苯磺酸盐(CholesterylIsothiuroniumTosylate)的反应为例,R.GPearson提出反应机理:A+C,无:龙格2库塔(Runge22卡罗(MonteCarlo)法。龙格-库塔(Runge-Kutta)法3设微分方程为dC dt=f(c,

2、t);初始条件是t=0时,C=C0;把所需计算的时间t分为若干个格点,(最简单可等间距)相邻两点时间间距为t;对应的浓度C由于第一点(t=0;C=C0)已知,以后各点则由前一点的数值逐个根据以下公式推算出来:(第i个格点ti,时为Ci:由前一点ti-1时的Ci-1计算)6;Ci=Ci-1+(K1+2K2+2K3+K4)t式中:K1=f(Ci-1,ti-1);2,ti-1+t 2);K2=f(Ci-1+K1t2,ti-1+t 2);K3=f(Ci-1+K2tK4=f(Ci-1+K3t,ti-1+t);其动力学方程:设t=0时A=A0;B=C=D=0;dA dt=-(k1+k3)A;dB dt=d

3、C dt=k1A-k2BC=2k1A-k2B;2dD dt=k2BC+k3A=k2B+k3A对于微分方程组,在计算时,须把浓度C和函数f都视为向量;同样可逐点计算各组分浓度3。图1为上述制备(CholesterylIsothiuroniumTosylate)反应;-1A0=0.004425mol L;k1=0.0044min;k2=即便是这表面看来简单的反应;虽然可得到解析解2:1 2( B=A0)k1+k32.96min-1mol-1L;k3=0.0021min-1。图中曲J0(2i1(1()11)+iH0(2i)H0+i-(k1+k3)t其中:=e;=iJ1(2i) H1(2i);(k1+

4、k3)2;=k1k2A0J0和J1为虚宗量的贝塞尔(Bessel)函数;H0和H1为虚宗量的汉克(Hankel)函数;i是单位虚数。由于它不能仅用初等函数描述;(其中有含复数的贝塞尔函数和汉克函数);要得到某时刻的浓度,仍需求助计算机。为此,提出了多种数值解法。此地,为开设课程,培养学生结合物修稿日期:2000208220线为解析解,为Runge2Kutta法计算值。由图1可见,符合的很好。用Runge2Kutta法解微分方程运算快,程序相对有也较短,但由于涉及向量,需一些线性代数的知识,同学解微分方程时有点困难,容易出错。2蒙特-卡罗(MonteCarlo)法蒙特2卡罗法是一种抽样统计法,常

5、用于计算机模拟,从1963年L.JChaad最早引入动力学方程求解至今已近半个世纪,其方法随着计算机的发展而不断改进,文献4介绍的模拟方法由于把所取的随机数与速率常数比较而决定是否反应,仅适用于一级反应。本文作了如下的改进:将所取的在0和1之间的随机数与含速率常数和各组分浓度的函数比较106化学世界2000年率;)如其数值大于零,说明它的数量增加,如小于零则说明它减少;因此可按下面步骤模拟反应过程:把一个单位时间分为M(2000)等分,每一等分作为模拟中的一步,(每一步相当1 2000单位时间)设第i步结束时A,B,D的浓度分别为a,b,d,它们分别都用N(10000)个球来代表;每一个球分别

6、对应于a N;b N和d N;(i+1)浓度时:A,之间的随机数,统计1+S;(dt=1 2000)则反应Sa N,经过此步留下未反应的A1-S N)a再计算B,令k=(k1A B-k2B)dt1来决定是否反应,。方法如下:(1)简单n级反应:n-dC dt=kC;t=0时,C=a。根据速率常数中使用的时间单位(如分或秒等)确定若k>0说明经过此步B的浓度增加;取N个0和1之间的随机数,统计其中小于k的个数S;B增加的浓度为Sb N;所以经过此步后B的浓度为(1+S N)b;若k<0;说明经过此步B的浓度减少;令k=-(k1A B-k2B)dt;同样取N个0和1之间的随机数,统计其

7、中小于k的个数S;经过此步后B的浓度等于(1-S N)b;如k=0则经单位时间,如1分(或1秒)定为1个单位时间;将1个单位时间分为M等分(如2000等分)计算每个间隔点的浓度C;如t=i 2000时的浓度为Ci;设t=(i-1) 2000时(前一点)的浓度为Ci-1;将Ci-1用N个完全相同的球(即样本)代表;每个球代表的浓度为Ci-1 N;(在实际计算中取N应足够大;如N=10000);由于每个球完全相同,其反应的概率当然也相同,就等于整体的反应概率。将微分方n程-dC dt=kC两边同除C再乘dt:n-1-dC C=kCdt;-dC为浓度C减少数量;-dC C即减少的百分数,也就是反应体

8、系整体的反应概率;从而也是每个球的反应概率;故可模拟如下:第i步计算时Ci过此步B的浓度不变。仍为b。然后再计算D;因D是最终产物,随时间增加而增加;令k=(k3A D+k2B2 D)dt同样取N个随机数,统计其中小于k的个数S,经过此步后D的浓度为(1+S N)d;如此一步步计算直至计算到所需要计算的时间为止。图1为上述制备CholesterylIsothiuronium取N=10000个在0和1之间的随机数;统计其中小于kCn-1dt的随机数个数S;(式中dt即相邻两点时间的间隔dt=1 2000;kCn-1dt当n=1就等于kdt;此时和文献4一致)留下未反应的球为N-S;Ci=(N-S

9、)Ci-1 N=(1-S N)Ci-1;再同法进行下一步模拟计算。直至最后一点。图2为计算结果,初浓度为1.5mol L。(2)复杂反应以上述制备CholesterylIsothiuroniumTosy2late的反应为例;将动力学方程组改写为:A=-(k1+k3)dtdAB=(k1A B-k2B)dtdBD=(k3A D+k2B2 D)dtdDA,dB B,dD D分别代表A,dAB,D的变化百分数;(绝对值即该物质反应的概图2图中圆点为模拟值,从上而下依次为0,1,2,3级反应。k=0.5(下转第184页)184化学世界ua2000年=O的吸收峰面积;AA1600为脲C=O的吸收峰面积;为

10、苯环C=C的吸收峰面积。计算结果表明随着降解时间的延长,(Xb)CO和Ch都在逐渐减小,也就是说硬段之间的氢键在降解氨酯软段和硬段的相分离。在经过微生物降解之后,随着氢键的减弱,这种软硬段的相分离程度也有所减弱。而这则会使微生物对聚氨酯的作用变得更为容易。(下转第87页)之后减弱了。有关资料显示,在硬段含量达到一定数值后,大部分氢键存在于硬段内部,并因此产生了聚(上接第106页)Tosylate的反应,初始条件t=0时A=0.mol L;B=C=D=0。2(MonteCarlo)前2塔(Runge2Kutta),和解析解(曲线)符合很好图3为连续二级反应:A+Bk1,;运行也更快一些。而用蒙特

11、,(所以称其为模拟法)但运算时间稍长一些。对于化学类学生来讲,用蒙特2卡罗法来解方程组更妥当;至少在结合专业学习程序编制方面而言不失为较好内容。(2)计算的误差:对于龙格2库塔法,格点分割得越细,误差就越小,运行时间也长。对于蒙特2卡罗法代表组分浓度的球N越多,把单位时间分割的越细(M越大),误差也越小,当然运行时间也越长。为了便于与文献值比较,本文N取10000,M取2000,同样在龙格2库塔法中也把一单位时间分为2000个格点,计算结果所有点和解析解比较,相对误差都小于2%;已满足一般动力学需要;如需更精确,可格点分割的更细一些。(3)用蒙特2卡罗法时,如速率常数数值太小;有可能N个随机数全部比反应的概率小(因k很小,dt=1 M也很小),在这种情况下,只要把时间的单位换成大一点就行了(如把分换成小时)否则就有可能运行很长时间,浓度不改变。参考文献:1穆尔JW,等著.孙承锷,等译.化学动力学和历程2匀+Dk2E初始条件:t=0时:A=B=1.5mol L;C=D=E=0;-1-1-1k1=7.5molLmin;k2=5.2molLmin-1;(解析解可根据模拟结果和解析解符合也很好。文献5方法解得)相化学反应研究(第二版).北京:科学出版社,1987.361,369.2ChaadLJ,AmJ.Chem.Soc,1963,(85):3588.图3