北京化工大学2003-2004学年第一学期 《计算化学》期末试卷答案 9 1.计算原理（化学原理和计算方法）(20分) 反应 A+B C+D B+C->D+E 其动力学方程式为一个常微分方程组： dIA]=-[A][B]+(CID] At d[B]-[A][B]+k.(CID]-&,[BJC] dIC]=IA]IB]-&(CID]-&,[B]C] d d[D]-k[A]B]-&[CID]+[BC] dt d回=kBC] 利用Run e-Kutta法解此方程组，可得不同时刻[A],B,[C],D,E值 2.程序框图 30分) 开始 输入：[A,B,[C,D,E初值 时间间隔H,计算精度EPS 计算：积分步长H0 调用变步长Runge--Kuta积分法子程序解动力学方程组计算时刻t 的[A,B,[C],D,E值 (其中调用计算动力学方程组的子程序计算DY值) 列表输出：0-200min内每隔10min所对应的[A,B,[C],D],E]值
北京化工大学 2003-2004 学年第一学期 《计算化学》期末试卷答案 A 1．计算原理（化学原理和计算方法）（20 分） 反应 k1 k2 A+B C+D B C D+E 3k + ⎯⎯→ 其动力学方程式为一个常微分方程组： 1 2 1 23 1 23 1 23 3 d[A] [A][B] [C][D] d d[B] [A][B] [C][D] [B][C] d d[C] [A][B] [C][D] [B][C] d d[D] [A][B] [C][D] [B][C] d d[E] [B][C] d k k t kkk t kkk t kkk t k t =− + =− + − = −− = −+ = 利用 Runge-Kutta 法解此方程组，可得不同时刻[A]，[B]，[C]，[D]，[E]值。 2．程序框图（30 分） 开始 输入：[A]，[B]，[C]，[D]，[E]初值 时间间隔 H， 计算精度 EPS 计算：积分步长 H0 调用变步长 Runge-Kutta 积分法子程序解动力学方程组计算时刻 t 的[A]，[B]，[C]，[D] ，[E]值 （其中调用计算动力学方程组的子程序计算 DY 值） 列表输出：0-200min 内每隔 10min 所对应的[A]，[B]，[C]，[D] ，[E]值
结柬 3.源程序(30分) PROGRAM MAIN IMPLICIT REAL*8(A-H.O-Z) DIMENSION YO(10).Y1(10),Y2(10) OPEN(6 FILE-TESTA TXT STATUS-UNKNOWNY DATA Y0/1.0.1.0.0.0.0.0.0.0,5*0.0/ T0=0.0 T1=100.0 H0=10 EPS-1.0E-6 H-HO WRITE(6.11)N.TO.TI.HO.EPS WRITE(6,22)T0,H,(Y00,I-1,N) M=IDINT((T1-TOVHO+10) D08010=1M CALLRK4(T0.Y0.Y1.H0.N.K) D030KK=1,100 K-K+K H=HO/K CALLRK4(T0.Y0.Y2.H.N.K) Es-00 D010J=1,N 10 ES=ES+ABS((Y2(-Y10VY20) IF(ES.LT.EPS)GOTO 50 D020J=1N YI()-Y2() CONTINUE WRITE(6,33) STOP 50 T0=T0+H0 D060J=1,N 60 Y0J)=Y2) WRITE(6.22)TO.H.(Y2(D)I=1.N) 80 CONTINUE WRITE(6.44) FORMATUIX N=122X TO-F6 12X TI=F8 L2X H0=G8.3,1X'EPS=,G10.320XY,e12,3,4N& 1X96(1H-)/5X,t/min,3X,'H/min,6X,'[A/molL-1',& 6X.[Bl/molL-I'.6X.'[C]/molL-I'6X.'[DVmolL-I'.6X.'[E]/molL-1'& 1X,961H-》
3．源程序（30 分） PROGRAM MAIN IMPLICIT REAL *8(A-H,O-Z) DIMENSION Y0(10),Y1(10),Y2(10) OPEN(6,FILE='TESTA.TXT',STATUS='UNKNOWN') DATA Y0/1.0,1.0,0.0,0.0,0.0,5*0.0/ N=5 T0=0.0 T1=100.0 H0=10 EPS=1.0E-6 H=H0 WRITE(6,11)N,T0,T1,H0,EPS WRITE(6,22)T0,H,(Y0(I),I=1,N) M=IDINT((T1-T0)/H0+10) DO 80 I0=1,M K=1 CALL RK4(T0,Y0,Y1,H0,N,K) DO 30 KK=1,100 K=K+K H=H0/K CALL RK4(T0,Y0,Y2,H,N,K) ES=0.0 DO 10 J=1,N 10 ES=ES+ABS((Y2(J)-Y1(J))/Y2(J)) IF(ES.LT.EPS)GOTO 50 DO 20 J=1,N 20 Y1(J)=Y2(J) 30 CONTINUE WRITE(6,33) STOP 50 T0=T0+H0 DO 60 J=1,N 60 Y0(J)=Y2(J) WRITE(6,22)T0,H,(Y2(I),I=1,N) 80 CONTINUE WRITE(6,44) 11 FORMAT(/1X,'N=',I2,2X,'T0=',F6.1,2X,'T1=',F8.1,2X, & 'H0=',G8.3,1X,'EPS=',G10.3//20X,'Y(I),I=1,2,3,4,...,N'//& 1X,96(1H-)/5X,'t/min',3X,'H/min',6X,'[A]/molL-1', & 6X,'[B]/molL-1',6X,'[C]/molL-1',6X,'[D]/molL-1',6X,'[E]/molL-1'& /1X,96(1H-)) 结束
2 F0RMAT1X2F8.2,2X5G15.6) 33 FORMAT(/IX.'FAILED TO FIND STABLE SOLUTION IN MAX IT/) ) SUBROUTINE F(Y,DY) IMPLICIT REAL8(A-H.O-ZY REAL*8K1K2K3 DIMENSION Y(10).DY(10) DATA K1/1.25E-3,K22.6E-4/,K3/8.4E-3 DY(1=-K1*Y(1)*Y(2)+K2*Y(3)Y(4) DY(2=K1*Y(1)*Y(2+K2*Y(3)*Y(4)-K3*Y(2)*Y(3) DY(3)=K1*Y(1)*Y(2K2*Y(3)*Y(4)K3*Y(2)*Y(3) DY)-K1Y(Y2)-K2-Y(3)Y(4)+K3Y()Y() DY(5)=K3*Y(2)*Y(3) RETURN END DIMENSION YO(N).Y(N) DIMENSION Y2(20).DY(20),RK(20) X=XO D05=1,N 5 Y()-YOND) D050L=1,M CALL F(YDY) DO I0ILN RK(D=DY(D 10 Y2=Y+0.5*H*DY CALLF(Y2.DY DO 20 I=1.N RKT=RKI①+2◆DYD 20 Y20Y(D+0.5*HDY(I) CALLF(Y2,DY) DO30I=1.N RK=RKI+2*DY④ 30 Y2I=Y(D+H◆DYID CALL F(Y2.DY) D0401=1N RK(I)-RK(I)+DY() 40 Y(I)-Y(I)+RKHW 50 CONTINUE RETURN
22 FORMAT(1X,2F8.2,2X,5G15.6) 33 FORMAT(/1X,'FAILED TO FIND STABLE SOLUTION IN MAX IT'/) 44 FORMAT(1X,96(1H-)) END SUBROUTINE F(Y,DY) IMPLICIT REAL *8(A-H,O-Z) REAL *8 K1,K2,K3 DIMENSION Y(10),DY(10) DATA K1/1.25E-3/,K2/2.6E-4/,K3/8.4E-3/ DY(1)=-K1*Y(1)*Y(2)+K2*Y(3)*Y(4) DY(2)=-K1*Y(1)*Y(2)+K2*Y(3)*Y(4)-K3*Y(2)*Y(3) DY(3)=K1*Y(1)*Y(2)-K2*Y(3)*Y(4)-K3*Y(2)*Y(3) DY(4)=K1*Y(1)*Y(2) -K2*Y(3)*Y(4)+K3*Y(2)*Y(3) DY(5)=K3*Y(2)*Y(3) RETURN END SUBROUTINE RK4(X0,Y0,Y,H,N,M) IMPLICIT REAL *8(A-H,O-Z) DIMENSION Y0(N),Y(N) DIMENSION Y2(20),DY(20),RK(20) X=X0 DO 5 I=1,N 5 Y(I)=Y0(I) DO 50 L=1,M CALL F(Y,DY) DO 10 I=1,N RK(I)=DY(I) 10 Y2(I)=Y(I)+0.5*H*DY(I) CALL F(Y2,DY) DO 20 I=1,N RK(I)=RK(I)+2*DY(I) 20 Y2(I)=Y(I)+0.5*H*DY(I) CALL F(Y2,DY) DO 30 I=1,N RK(I)=RK(I)+2*DY(I) 30 Y2(I)=Y(I)+H*DY(I) CALL F(Y2,DY) DO 40 I=1,N RK(I)=RK(I)+DY(I) 40 Y(I)=Y(I)+RK(I)*H/6 50 CONTINUE RETURN
END 0T1= 100.0H0=10.0 EPS=.100E-05 Y0.=1,2,3,4N Vmin H/min ICymolL- [DVmolL-I 00 10.00 100000 100000 000000 000000 .000000 1000 987657 987154 118414E-01 12456E01 502113E03 975627 2324520E-0 262947E.01 192133E0 250 3195 022 41366 500 95250 94546 4045660 545352F-0 7039300 50.0 50 94140 .93087 .480618E-01 691260E-01 105321E-01 500 93061 91608 548602E-01 839168E-0 145283E-01 5.00 92011 13821 87118 7 288113E 10000 5.0 890355 85622 755091E0 143780 341355E-0 10.0 5.00 09 84133 793524E-01 15866 96583E-0 5.0 g634 110 电 858139E.0 57034IE4 831 2168 2990E-0 16000 500 8379s7 6900 930259E-01 231000 689868E-0 170.00 5.00 830101 5509 945967E.01 244901 750024E.0 5.00 52436 9%5449E-0 25858 810196E-D 1900 72796 79927E0 370231E.0 20.0 922E 28520 1.计算原理（化学原理和计算方法）(20分) 合成甲醇的反应是可逆 放热和物质的量减少的反应。由于反应的可逆性 通常产品产率较低。设反应过平衡时生成甲醇的物质的量为xmol,有 C0+2H2 ±CH3OH 起始物质的量mol 0
END 4．运行结果。（20 分） N= 5 T0= .0 T1= 100.0 H0=10.0 EPS= .100E-05 Y(I),I=1,2,3,4,...,N --------------------------------------------------------------------------------------------------------------------------------- t/min H/min [A]/molL-1 [B]/molL-1 [C]/molL-1 [D]/molL-1 [E]/molL-1 --------------------------------------------------------------------------------------------------------------------------------- .00 10.00 1.00000 1.00000 .000000 .000000 .000000 10.00 2.50 .987657 .987154 .118414E-01 .128456E-01 .502113E-03 20.00 2.50 .975627 .973705 .224520E-01 .262947E-01 .192133E-02 30.00 2.50 .963910 .959773 .319535E-01 .402268E-01 .413663E-02 40.00 5.00 .952504 .945465 .404566E-01 .545352E-01 .703930E-02 50.00 5.00 .941406 .930874 .480618E-01 .691260E-01 .105321E-01 60.00 5.00 .930611 .916083 .548602E-01 .839168E-01 .145283E-01 70.00 5.00 .920115 .901164 .609337E-01 .988358E-01 .189511E-01 80.00 5.00 .909911 .886179 .663565E-01 .113821 .237321E-01 90.00 5.00 .899994 .871182 .711950E-01 .128818 .288113E-01 100.00 5.00 .890355 .856220 .755091E-01 .143780 .341355E-01 110.00 5.00 .880989 .841331 .793524E-01 .158669 .396583E-01 120.00 5.00 .871888 .826549 .827730E-01 .173451 .453388E-01 130.00 5.00 .863045 .811904 .858139E-01 .188096 .511411E-01 140.00 5.00 .854452 .797418 .885139E-01 .202582 .570341E-01 150.00 5.00 .846102 .783111 .909076E-01 .216889 .629905E-01 160.00 5.00 .837987 .769000 .930259E-01 .231000 .689868E-01 170.00 5.00 .830101 .755099 .948967E-01 .244901 .750024E-01 180.00 5.00 .822436 .741416 .965449E-01 .258584 .810196E-01 190.00 5.00 .814984 .727961 .979927E-01 .272039 .870231E-01 200.00 5.00 .807740 .714740 .992602E-01 .285260 .929998E-01 --------------------------------------------------------------------------------------------------------------------------------- B 1.计算原理（化学原理和计算方法）（20 分） 合成甲醇的反应是可逆、放热和物质的量减少的反应。由于反应的可逆性， 通常产品产率较低。设反应过平衡时生成甲醇的物质的量为 x mol，有 CO+2H2 CH3OH 起始物质的量/mol CO,0 n H2 ,0 n 0
平衡物质的量/mol nco0x”,0-2x 令起始反应物的总物质的量为1mol,即：no0+m,0=1 反应达平衡时系统中总物质的量为：no0-x+nH,0-2x+x=1-2x 故平衡时系统中各物质的分压为： Pco=pYeo =p-("co 1-2x P4=m4=p-(0-3 1-2x Pcuon PYanon p-(1-2x 式中：y为物质的量分数，p为系统总压。 代入平衡关系式得： K= -22p (1) (c 1-2x 1-2x (1)式经整理后得 1+K pimo(incog K p'nco(r0 41+Kp2) 41+Kp） (2) 令S-1+K,Pa4ne+n） Q=KPrcoa( 41+Kp) 41+K,p) (2)式可写为x2-x2+5x-Q=0 用D表示原料气体中C0与H2的体积比，有 nco0+n40=1 (3) 解(3)得：aoo=D1+D) n4.0=1/1+D) 令：fx)=x3-x2+5Sx-Q 调用Newton--Raphson法求解方程f(x)=0 取初值x。=0.00
平衡物质的量/mol CO,0 n -x H2 ,0 n -2x x 令起始反应物的总物质的量为 1mol，即： 1 nCO,0 + nH2 ,0 = 反应达平衡时系统中总物质的量为： n x n 2x x 1 2x CO,0 − + H2 ,0 − + = − 故平衡时系统中各物质的分压为： ) 1 2 ( ) 1 2 ( ) 1 2 ( CH OH CH OH H ,0 H H CO,0 CO CO 3 3 2 2 2 x x p py p x n x p py p x n x p py p − = = ⋅ − − = = ⋅ − − = = ⋅ 式中：y 为物质的量分数，p 为系统总压。 代入平衡关系式得： CO,0 H ,0 2 2 2 ( ) 1 2 ( )( ) 12 12 p x p x K n x n x p p x x − = − − − − （1） （1） 式经整理后得 22 2 2 22 3 2 H ,0 CO,0 H ,0 CO,0 H ,0 2 2 1 (4 ) ( ) 0 4(1 ) 4(1 ) p p p p K pn n n K pn n xx x Kp Kp ⎡ ⎤ + + −+ − = ⎢ ⎥ + + ⎣ ⎦ （2） 令 22 2 2 2 H ,0 CO,0 H ,0 CO,0 H ,0 2 2 1 (4 ) ( ) 4(1 ) 4(1 ) p p p p K p n n n K pn n S Q Kp Kp + + = = + + （2） 式可写为 0 3 2 x − x + Sx − Q = 用 D 表示原料气体中 CO 与 H2 的体积比，有 n n D n n = + = CO,0 H ,0 CO,0 H ,0 2 2 / 1 （3） 解（3）得： 1/(1 ) /(1 ) H ,0 CO,0 2 n D n D D = + = + 令： f x = x − x + Sx − Q 3 2 ( ) 调用 Newton-Raphson 法求解方程 f (x) = 0 取初值 x0 = 0.00