第三章 化学反应动力学的计算

第三章化学反应动力学的计算化学反应的速度各不相同，有的反应速度极快，只要几个毫微秒就达到平衡(接近扩散速度，如无机酸碱中和)，有的反应速度极慢，几乎看不到变化(如自然界的某些变化)。大部分有机化学反应可用常规方法测量，对某些快速反应则可用停留法、驰豫法等测量。不论反应速度的快慢，动力学方程都是类似的。一、化学反应动力学方程反应物浓度随时间的变化绝大部分不是线性关系，而是一条曲线，见图3-1。反应速度公式可用微分方程来表示。具有简单级数的化学反应的反应速度公式可用积分式表示：一级如：0AA1AdcAC=a,-=kcdt生成物：，㏑CA=㏑a–K1t二级A+A→产物CA0=a2A2A2Adc11-kC,=+ktdtca对于反应1-1kkAB这一可逆反应生成物反应物初始条件t=0a0时间t时t=ta-xx达到平衡时，B的浓度为Xe，则可逆反应的速度积分式为：级数：1-11-10kAAee1A-1Bk0e0C=adcxxAB=-kC+kC:=ktdtax-xC=0ln1-21-10Ak0AeeeB1A-1BCkee0CC=adcxax+x(a-x)AB+CC=0=-kC+kCC:=ktdt2a-xa(x-x)C=0ln二、常微分方程的解化学反应动力学方程是用微分方程表示的，对于简单的反应，可直接求得微分方程的解。微分方程：()(1)(,,,......)......(1)nnyfxyyy在区间axb的解，是指()yx,这样一个函数，在所述区间内存在导数()(),(),......()nxxx。且对于区间axb内的每一个x，等式(1)都成立。()yx就叫微分方程的解。（1）对于一级反应，kA2t2A1t1ABd[A]k[A]dtd[A]kdt[A]d[A]kdt[A]1212[A]lnk(tt)[A]t1为0，初始浓度为A0，则上式为：-kt00[A]lnktA=Ae[A]　　或　　(2)对于一级可逆反应:12kk12ABd[A]=-k[A]+k[B]dt开始时只有[A0]。反应时若A的浓度为[A]，则B的浓度为[A0]-[A]，那么：1202012212012d[A]=-k[A]+k([A]-[A])dt=k[A]-(k+k)[A]k=-(k+k)([A]-[A])k+k平衡时：k1[Ae]=k2[Be]即e0e1ee2[B][A]-[A]k==[A][A]k则2e012k[A]=[A]k+k12ed[A]=-(k+k)([A]-[A])dt0At12A0ed[A]-=(k+k)dt[A]-[A]0e12e[A]-[A]ln=(k+k)t[A]-[A]对于一阶线性微分方程dy+P(x)y=Q(x)dx，其通解可表示成：-PdxPdxy=e(eQdx)+c......(2)所以，对于一阶线性微分方程可用（2）式求通解。但求一阶非线性微分方程的解析解非常困难甚至不可能。例如22dy=x+ydx，就无法求其解析解。在化学反应中，某些稍稍复杂一点的反应，如甲基磺酸盐胆甾醇酯在甲醇和三氯甲烷溶液中与硫脲反应，其解析解的求解非常困难。而在原油裂解、催化成环等反应中可以考虑的动力学方程可达几十个，因此求方程的解析解更是无法实现，这就不得不借助于电子计算机来求微分方程的数值解。二、常微分方程的数值解（Rung-Kutta法）1．泰勒公式（泰勒Taylor中值定理）在计算机的近似方法处理中经常遇到要使用泰勒公式，在此讨论一下该公式。任意一个函数f(x)，如果在某区间内，一直到(n+1)阶导数都存在，而且连续。设点xa是这区间的一个内点，则当x取这区间的任何值时，()fx可以依()xa的方幂展开：()2()()()()()()()......()()1!2!!nnfafafafxfaxaxaxaRnxn其中()Rnx是余项。令(),xhxxa代入上式，则()2()()()()......()1!2!!nnfxfxffxhfxhhhRnxhn二元函数的泰勒展开：一个二元函数(,)fxy，把(,)fxy在00(,)xy点附近展开成泰勒级数，只取二阶导数项，忽略高次项：''00000000(,)(,)(,)()(,)()xyfxyfxyfxyxxfxyyy22100000002[(,)()(,)()()xyxfxyxxfxyxxyy220000000(,)()()(,)()]yxyfxyxxyyfxyyy令00,xxxyyy则上式改写成：''00000000(,)(,)(,)(,)xyfxxyyfxyfxyxfxyy222210000002[(,)2(,)(,)]xyxyfxyxfxyxyfxyy2．Eular法方向场方程222xyc(1)是一个圆的方程，圆心在原点，半径为c。(1)式的微分方程为：dyxdxy(2)(1)与(2)式之间的关系为：在圆周上的点作一短小切线，这样的切线一定符合方程(2)，而这些短小切线所构成的图形符合方程(1)，(1)和(2)式在几何上是有明确的意义的。对一个微分方程(,)dyfxydx，设函数(,)fxy在xy平面的定义域是G，在区域G内每一点(,)xy画出斜率等于(,)fxy的小线段，因而在G内每一点确定一个方向。我们说，这个微分方程在区域G内定义一个方向场，经过G内若干个点，画出具有已斜率的小线段，即可得方向场的略图(图3-2)从方向场各点的略图可以推出微分方程的原函数图形。例如画出微分方程dyydxx的方向场略图，其解为：在xy平面上（除原点外）的若干个点，画出斜率等于yx的小线段，于是得到dyydxx的方向场略图图3－2(图3－3)。由图可见，dyydxx的解是直线族ycx（c是任意常数，这是微分方程的通解）。如果给定的初始条件(指定其一点)就是特解，比如给定初始条件003,3xy，则其解必定为,1yxc。图3－3画出微分方程dyxdxy的方向场略图。在xy平面内若干点(,)xy画出斜率等于xy的小线段，它垂直于该点至原点的连线，所得方向场如图3－4，其解是以原点为中心的圆族：222xyc(通解)从几何观点来看，所谓微分方程的解，就是这样的曲线，即在曲线上每一点的切线方向都是一致的，我们称这样的曲线为微分方程的积分曲线。实际问题中遇到的微分方程(,)dyfxydx，其解析解往往不能求得。但我们总可以画出它的方向场，从而图3－4对解的几何性态可以有所了解。微分方程的数值解是根据方向场概念而建立的。xy对于化学反应动力学的微分方程来说，只要体系是真实的，表达的方程是准确的，那就一定存在数值解。微分方程(,)dyfxydx的右端函数(,)fxy，是在xy平面上的带形区域R0(,)axby内确定了一个方向场求解初值的问题。(,)dyfxydx(0)y其几何意义相当于确定一条通过初始点00(,)xy的曲线()yyx，使曲线上每一点的切线方向与已给定的方向一致。常微分方程数值解使用的是离散方法，即找出一个有效的数值方法，计算自变量x的离散点012,,,...,nxxxx和对应()yx的近似值012,,,...,nyyyy。微分方程(,)dyfxydx从起点000(,)pxy开始(初始条件)，沿着该点的方向线段移动一段距离10()xx到新点111(,)pxy，然后改变方向，从点111(,)pxy的方向移动距离21()xx到第二个新点222(,)pxy，等等，如图3－5所示，就这样构成了一个“折线”函数。通常离散点x取等距离，即令10211...nnxxxxxxh图3－5Eular法图P0(x0,y0)P1(x1,y1)P2(x2,y2)P3P4yxx0x1x2x3x4常数h称为步长。微分方程(,)dyfxydx给定初值00(,)xy，代入到微分方程(,)dyfxydx，就可求得该点00(,)xy的斜率。移动h，在1x处求得1y，以点111(,)pxy再代入微分方程，求得点222(,)pxy。如此重复直到要求的点为止。这种方法称为Eular方法或折线法。它简单易懂，有明确的几何解释，但这一方法误差太大，故实用价值受到限制。Eular法是用递推公式1(,)nnnnyyhfxy，从00(,)xy点出发求(,)nnxy的过程。以dyxydx为例进行说明。初始条件000,0xy，方程的准确解为1xyex。用Eular方法和准确解分别计算当0.1,0.2,0.3,...,1.0x时对应的()yx。选择步长h=0.1,用Eular公式1(,)nnnnyyhfxy计算。1000(0.1)(,)00.1(00)0yyhfxy2111(0.2)(,)00.1(0.10)0.01yyhfxy3222(0.3)(,)0.010.1(0.20.01)0.031yyhfxy.............................................表3－4方程dyxydx的Eular法和准确解的比较xy(Eular)y(准确解)0000.100.00517090.20.010.021400.30.03100.04990.40.06410.09130.50.11050.14370.60.17160.22210.70.24370.31380.80.34360.42550.90.45790.55961.00.59370.7183从表3－4的计算结果来看，Eular法与准确解还有不小的误差，并且随着计算次数的增加，积累误差增大。如果缩小步长h，取切线很小的范围，就能逼近函数的真实解。但随着步长的缩小，计算次数将增加，又会引起舍入误差。仍以dyxydx00(0,0)xy为例计算(0.1)y的值，观察步长对结果的影响。如果步长h=0.1(0.1)0y如果步长h=0.05000(0.05)(,)0yyhfxy(0.1)0.0025y如果步长h=0.025(0.025)0y(0.05)0.000625y(0.075)0.00189y(0.1)0.00381y表3－5列出不同步长h时所计算(0.1)y的结果。表3－5不同步长h计算(0.1)y值h(0.1)y0.10.00.050.00250.0250.003810.01250.004490.006250.004830.0031250.005000.00156250.005080.000781250.005130.0003906250.005150.00019531250.00516从表3－4查得(0.1)y的准确解值为0.005171。通常选择步长的方法是先用h作步长，算出0xx处的近似值()0hy，再用2h作步长计算0xx处的近似值2()0hy。如果两个0y值之间的绝对值不超过容许的误差范围：2()()0hhayy则认为选择步长h对()hy是满足了所要求的精度内的解。Eular方法产生误差的原因是由于函数()fx为曲线。用泰勒公式在nx处展开式子：2()()'()()...2nnnnhyxhyxhyxyxEular法只取了前二项而忽略了高次项，所以产生了误差。3．Runge-Kutta方法3.1常微分方程的Runge-Kutta方法Runge-Kutta方法是建立在泰勒公式基础上的一种方法。通常采用的是四阶R－K公式，即考虑了泰勒公式中四次项，而Eular公式只取了一次项。故R-K公式比Eular公式有了很大的改进。R-K方法在求解范围大、精度要求主的情况下是一种比较好的方法，并且计算工作量不算太大，所以在化学化工中应用颇多。四阶的R-K公式为：1123412132431(22)6(,)11(,)2211(,)22(,)nnnnnnnnnnyyKKKKKhfxyKhfxhyKKhfxhyKKhfx