北京化工大学2006一一2007学年第一学期 《计算化学》期末考试试卷标准答案 第一部分填空题（共30分，每空一分） 1.数值计算、非数值计算。 2.建立数学模型。 3.二分法：①函数x)在区间[a,b]内连续，②函数在a、b两点之值a)与b) 异号：牛顿迭代法：①∫”x)有明确表示式且容易计算②容易确定根的初值xo: 二分法：压力对化学反应中反应物平衡转化率的影响。 4.线性方程组的系数矩阵为对角占优势。 5.插值法。 6.一阶导数，由c~1数据求化学反应速率。 7.y=+h+e;y=6+hx+b3++bxn+s。最小二乘法。 8.x)函数形式己知，但其积分不能表示成初等函数的闭合形式，x)函数形 式未知，但其离散数据表已给出。(1)用插值程序求任意点的函数值。(2)用 Simpson求积程序计算[a,b]区间中离散点下的面积。 9.Euler法、Runge-.Kutta法。化学动力学。 10.量子化学/结构化学。 11,不用求导，甚至没有目标函数表达式时也可使用。 I2.Monte Carlo法、分子动力学方法。 13.Chemwindow/Chem Office Origin,Gaussian,Materials studio. 14.稳态模拟、动态模拟
1 北京化工大学 2006——2007 学年第一学期 《计算化学》期末考试试卷标准答案 第一部分 填空题（共 30 分，每空一分） 1．数值计算、非数值计算。 2．建立数学模型。 3．二分法：①函数 f(x)在区间[a,b]内连续，②函数在 a、b 两点之值 f(a)与 f(b) 异号；牛顿迭代法：①f ’(x)有明确表示式且容易计算②容易确定根的初值 x0； 二分法：压力对化学反应中反应物平衡转化率的影响。 4．线性方程组的系数矩阵为对角占优势。 5．插值法。 6．一阶导数，由 c~t 数据求化学反应速率。 7．y=ax+b+ε； = + + + ⋅⋅⋅+ + ε m m y b b x b x b x 0 1 1 2 2 。最小二乘法。 8．f(x)函数形式已知，但其积分不能表示成初等函数的闭合形式，f(x)函数形 式未知，但其离散数据表已给出。（1）用插值程序求任意点的函数值。（2）用 Simpson 求积程序计算[a,b]区间中离散点下的面积。 9．Euler 法、Runge-Kutta 法。化学动力学。 10．量子化学／结构化学。 11．不用求导，甚至没有目标函数表达式时也可使用。 12．Monte Carlo 法、分子动力学方法。 13．Chemwindow/Chem Office ，Origin，Gaussian，Materials studio。 14．稳态模拟、动态模拟
第二部分上机题（共70分） A 1.计算原理（化学原理和计算方法）(15分) 带式板应图装型东：业=5式大业 式中，Fo为反应物A的进气流速，”A为以浓度表示的反应速率。 反应管长为 L=4F 1 dx 其中反应物A的进气流速FA0=xAoF%=1/1+0.5)F=2F/3(5分) kAL-A2 2AKA-1 XA=0 no 0 n。/2 XA no(1-xA)noxA/2 no/2 可得P4=20-2p 3-XA pa 故，-=p-P (1) TAe =TAp /RT
2 第二部分 上机题（共 70 分） A 1． 计算原理（化学原理和计算方法）（15 分） 等温管式反应器模型为： A 0 A, R A0 A d 1 V x r F x c ∫ = 式中， FA0 为反应物 A 的进气流速， c rA, 为以浓度表示的反应速率。 反应管长为 A 0 A, 2 A0 A d 4 1 x D r F L x c ∫ = π 其中反应物 A 的进气流速 FA0 = xA0F0 = 1/(1+ 0.5)F0 = 2F0 / 3 （5 分） 2A A2 kA,1 kA,-1 0 xA = n0 0 n0 / 2 Ax (1 ) 0 A n − x n0 xA / 2 n0 / 2 可得 p x x p A A A 3 - 2(1− ) = p x x p A A A 3 - 2 = 故 p x x p k x x k k p k p t p r p A A A,-1 2 2 A 2 A A,1 A,-1 A 2 A,1 A A A, (3 - ) 3 - 4(1 ) d d 2 − − = = − = − （1） rA,c = rA, p / RT I
所以反应管长L πD20r (2) πD2 (3)(5分 由于被积函数 L有确定的形式（式(1）的倒数)，所以可直接利用Si即som求积 分的方法计算可得到(3)式S的值，其中调用自编子程序计算函数， 1的值，然后代入(2) 式，可得不同的转化率时反应管长。因积分下限0，为了防止出错，计算时取0.0001。(5 分) 程序框图(20分) 开始 输入：数据点数F0,EPS,P,T,XA0,XA①，(I=1,N)(5分) A=0.000L,B=XA①，(=1，N) (5分) 调用Simpson法子程序计算积分S (用自编子程序计算函数值)，L=4F,RTS/πD2 (5分) 输出：不同的转化率时反应管长 (5分) 结束 3. 源程序(20分) PROGRAM MAIN (10分) DOUBLE PRECISION EPS.P.XA.A.B.S2.L.D.FO.FAO.XAO DIMENSION XA(5) EXTERNALF OPEN(6.FILE-2007TESTAOUT TXU STATUS-UNKNOWN)
3 所以反应管长 A 0 A, 2 A0 A 0 A, 2 A0 A A d 4 1 d 4 1 x D r F RT x D r F L x p x c ∫ ∫ = = π π （2） 令 A 0 A, A d 1 x r S x p ∫ = （3） （5 分） 由于被积函数 p rA, 1 有确定的形式（式（1）的倒数），所以可直接利用 Simpson 求积 分的方法计算可得到（3）式 S 的值，其中调用自编子程序计算函数 p rA, 1 的值，然后代入（2） 式，可得不同的转化率时反应管长。因积分下限 0，为了防止出错，计算时取 0.0001。（5 分） 2． 程序框图（20 分） 3． 源程序（20 分） PROGRAM MAIN （10 分） DOUBLE PRECISION EPS,P,XA,A,B,S2,L,D,F0,FA0,XA0 DIMENSION XA(5) EXTERNAL F OPEN(6,FILE='2007TESTAOUT.Txt',STATUS='UNKNOWN') 开始 输入：数据点数 F0,EPS，P，T，XA0，XA (I),（I=1，N） （5 分） A=0.0001, B= XA (I),（I=1，N） （5 分） 调用 Simpson 法子程序计算积分 S （用自编子程序计算函数值）， 2 A0 L = 4F RTS /πD （5 分） 输出：不同的转化率时反应管长 L （5 分） 结束
EPS-10E-6 ATA XA0.10.020,0.30,0.40.0.50/ 0) 00 FAO XA0 R=8.314 D010=l.5 A=0.0001 R-YACD A B EPS S2.F) S2FA043.14159265D/DRT WRITE(6 I XAOL 11 FORMAT(X.XA-.F12.4.4X.L-F12.4) 10 CONTINUE END FUNCTION FCXA) (5分) DOUBLE PRECISION XA.K1.K2.P D-1 0E5 DATA K1/8565E-5/K2/6.833/ F=1K1*4*(1-XA)*(1-XA)产p*p-K2*XA*3-XA)P(3-XA)*2) RETURN END SUBROUTINE SIMP(A.B.EPS.S2.F) (5分) DOUBLE PRECISION A.B.H.EPS.S2.A1 H=B-A T1=0.5H*(F(A+F(B)) N=l 91=0 10 8=0 D020K=1,N A1=A+K-0.5*H S=S+E(A1 T2=0.5*(T1+H*S) S2=T2+T2.T1V30 IF(ABS(S2).LE.1.0)D=ABS(S2-S1) If(ABS(S2)GT LONTHEN D=ABS((S2-S1VS2) END IE IF(D.LT.EPS)GOTO50 NEN+N H=0.5H T1=T2 S1=S2 GOTO 10 ￥
4 EPS=1.0E-6 DATA XA/0.10,0.20,0.30,0.40,0.50/ D=0.10 P=1.0E5 F0=10.0 XA0=1/(1+0.5) FA0=F0*XA0 R=8.314 T=911.0 DO 10 I=1,5 A=0.0001 B=XA(I) IF(B.EQ.0)STOP CALL SIMP(A,B,EPS,S2,F) L=S2*FA0*4/3.14159265/D/D*R*T WRITE(6,11) XA(I),L 11 FORMAT(3X,'XA=',F12.4,4X,'L=',F12.4) 10 CONTINUE END FUNCTION F(XA) （5 分） DOUBLE PRECISION XA,K1,K2,P P=1.0E5 DATA K1/8.565E-5/,K2/6.833/ F=1/((K1*4*(1-XA)*(1-XA)*P*P-K2*XA*(3-XA)*P)/((3-XA)**2)) RETURN END SUBROUTINE SIMP(A,B,EPS,S2,F) （5 分） DOUBLE PRECISION A,B,H,EPS,S2,A1 H=B-A T1=0.5*H*(F(A)+F(B)) N=1 S1=0 10 S=0 DO 20 K=1,N A1=A+(K-0.5)*H 20 S=S+F(A1) T2=0.5*(T1+H*S) S2=T2+(T2-T1)/3.0 IF(ABS(S2).LE.1.0)D=ABS(S2-S1) IF(ABS(S2).GT.1.0)THEN D=ABS((S2-S1)/S2) END IF IF(D.LT.EPS)GOTO 50 N=N+N H=0.5*H T1=T2 S1=S2 GOTO 10
50 RETURN END 运行结果(15分) 1.8765(每个3分) m 42 4000 13.5401 5000 58596 x≤0.4 当x>0.4以后，所需反应器管长剧增。因此设计反应器时应选取 ◇ 1.计算原理（化学原理和计算方法）(15分) 设反应达平衡时生成C0的物质的量为xmol,有 CH4+H2O CO+3H2 起始物质的量/mol 12 00 平衡物质的量mol 1-x2-x x 3x 反应达平衡时系统中总物质的量为：1-x+2-x+x+3x=3+2x 故平衡时系统中各物质的组成为： 1-x 2-x yal=3+2x' yH0=3+2x x 3x yc0=3+2x m=3+2x 式中：y为物质的量分数。 (5分) 代入平衡关系式得： K”=KP x(3x)3「 =1.247 (1) (1-x)2-x)(3+2x)p° (1)式经整理后得 x.(3x)3 -1.247=0 (2) 1-x)2-x)(3+2x)p° 令：f(x)= x-(3x)3 p -1.247 (5分) (1-x)2-x)(3+2x)p° 调用二分法求解方程f(x)=0
5 50 RETURN END 4． 运行结果（15 分） xA= .1000 L/m= 1.8765（每个 3 分） xA= .2000 L/m= 4.2871 xA= .3000 L/m= 7.6705 xA= .4000 L/m= 13.5401 xA= .5000 L/m= 58.5966 由计算结果可知：当 0.4 xA > 以后，所需反应器管长剧增。因此设计反应器时应选取 0.4 xA ≤ 。 B 1．计算原理（化学原理和计算方法）（15 分） 设反应达平衡时生成 CO 的物质的量为 x mol，有 CH4+H2O CO+3H2 起始物质的量/mol 1 2 0 0 平衡物质的量/mol 1-x 2-x x 3x 反应达平衡时系统中总物质的量为： 1− x + 2 − x + x + 3x = 3 + 2x 故平衡时系统中各物质的组成为： x x y x x y x x y x x y 3 2 3 , 3 2 3 2 2 , 3 2 1 2 4 2 CO H CH H O + = + = + − = + − = 式中：y 为物质的量分数。 （5 分） 代入平衡关系式得： 1.247 (1 )(2 ) (3 2 ) (3 ) ( ) 2 3 = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − − + ⋅ = = ∑ θ θ θ x p p x x x x n p p K K B n （1） （1） 式经整理后得 1.247 0 (1 )(2 ) (3 2 ) (3 ) 2 3 − = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − − + ⋅ θ x p p x x x x （2） 令： 1.247 (1 )(2 ) (3 2 ) (3 ) ( ) 2 3 − ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − − + ⋅ = θ x p p x x x x f x （5 分） 调用二分法求解方程 f (x) = 0
取值范围：0.0~1.0 (5分) 2.程序框图(20分) 开始 输入：温度T,压力P,标准平衡常数K,EPS N)(4分) 输入二分区间A=0.0,B=1.0(4分) 调用二分法子程序求解方程f(x)=0(4分) (其中调用自编子程序计算f(x)值)得C0的生成物质的量x 计算平衡时各组分的摩尔分数YC1,YH1,YC2,YH2(4分) 输出：T,P,x,YC1,YHl,YC2,YH2值(4分) 结束 3.源程序(20分) PROGRAM MAIN (10分) REALT.P.K.EPS.YC1.YH1,YC2.YH2.X EXTERNALF COMMONP OPEN(6,FILE-2007TESTB.TXT,STATUS-UNKNOWN) DATA T/900/ DATA EPS/IE-6/ WRITE(,*YPlease input PRESSURE/Pa:' READ(,*)P K=1.247 A=0.0 6
6 取值范围 ：0.0～1.0 （5 分） 2.程序框图（20 分） 3.源程序（20 分） PROGRAM MAIN （10 分） REAL T,P,K,EPS,YC1,YH1,YC2,YH2,X EXTERNAL F COMMON P OPEN(6,FILE='2007TESTB.TXT',STATUS='UNKNOWN') DATA T/900/ DATA EPS/1E-6/ WRITE(*,*)'Please input PRESSURE/Pa:' READ(*,*)P K=1.247 A=0.0 开始 输入：温度 T，压力 P，标准平衡常数 K，EPS N）（4 分） 调用二分法子程序求解方程 f (x) = 0（4 分） （其中调用自编子程序计算 f (x) 值）得 CO 的生成物质的量 x 输出：T，P，x，YC1，YH1，YC2，YH2 值（4 分） 结束 输入二分区间 A＝0.0，B＝1.0（4 分） 计算平衡时各组分的摩尔分数 YC1，YH1，YC2，YH2（4 分）
B=1.0 EPS-1.0E-6 CALLHALF(A,B.EPS,F,X) YC1=1-X)/3+2*X) YH1=2-X3+2X YC2=X3+2*X0 YH2=3*X/3+2*X WRITE(6,11)T.P.X,YC1.YHI,YC2.YH2 11 FORMAT(/5X,T=,F10.4.'K'/5X.'P=,G12.4,'Pa'l& /5X,'X=,F6.4//RESULT:'/5X,'y(CH4)=,F6.4& /5XvH20=F6.4/5XvCO=F6.4/5X.vH2=F6.4) END FUNCTION F(X) (5分) REAL X,P COMMON P F=(X(3x)*3)1-X/2-X)*(P3+2*X/1.0E5)*21247 RETURN END SUBROUTINE HALF(A.B.EPS,F.X) (5分) REALA.B.EPS YO=F(A) K 10 X-A1+B1)0.5 Y=F(X) IF(Y+YO.GT.0.0)AI-X IFY◆Y0.LE.0.0)B1= IF(ABS(B1-AIX.GT.EPS)GOTO 10 K=K+1 IF(K.GT.50)THEN WRITE(,YNO RESOLUTION 0I030 END I B1)*0.5 30 RN END 4运行结果(15分) (1)T=900.0000K p= .1000E+06Pa X=.7414 RESULT: y(CH4)=.0577 yH20)=.2808 y(C0)=.1654 yH2)F4962 (5分) T=900.0000K
7 B=1.0 EPS=1.0E-6 CALL HALF(A,B,EPS,F,X) YC1=(1-X)/(3+2*X) YH1=(2-X)/(3+2*X) YC2=X/(3+2*X) YH2=3*X/(3+2*X) WRITE(6,11)T,P,X,YC1,YH1,YC2,YH2 11 FORMAT(/5X,'T=',F10.4,'K'/5X,'P=',G12.4,'Pa'/& /5X,'X=',F6.4//'RESULT:'/5X,'y(CH4)=',F6.4 & /5X,'y(H2O)=',F6.4 /5X,'y(CO)=',F6.4 /5X,'y(H2)=',F6.4) END FUNCTION F(X) （5 分） REAL X,P COMMON P F=(X*((3*x)**3)/(1-X)/(2-X))*((P/(3+2*X)/1.0E5)**2)-1.247 RETURN END SUBROUTINE HALF(A,B,EPS,F,X) （5 分） REAL A,B,EPS Y0=F(A) A1=A B1=B K=0 10 X=(A1+B1)*0.5 Y=F(X) IF(Y*Y0.GT.0.0)A1=X IF(Y*Y0.LE.0.0)B1=X IF(ABS(B1-A1)/X.GT.EPS)GOTO 10 K=K+1 IF(K.GT.50) THEN WRITE(*,*)'NO RESOLUTION' GOTO 30 END IF X=(A1+B1)*0.5 30 RETURN END 4.运行结果（15 分） （1）T= 900.0000K P= .1000E+06Pa X= .7414 RESULT: y(CH4)= .0577 y(H2O)= .2808 y(CO)= .1654 y(H2)= .4962 （5 分） T= 900.0000K
P=.1000E+07P X=,2911 RESULT: C=0813 (H2F2438 (5分) (2)第一种情况更有利于反应产物的生成。(5分) Q 1.计算原理（化学原理和计算方法）(15分) 对于气相反应： CHCHO- →CH4+CO 1=0 PAO 0 0 1=i1 PA PAO-PA PAO-PA 1时刻时系统中总压为：p=PA十2(PA0~PA)=2pA0PA=PA0十△p 故时刻反应物的分压为P4=P40-P (5分) 化学反应速率方程： r=- 业=p” d 式中，r为化学反应速率，p为反应物的分压，1为时间，n为反应级数， k为反应速率常数，dp仙为反应物压力随时间的变化率。 两边取对数：nr=ln-史)=lnk+nlnp d y=Inr,x=In p,a=Ink,b=n,y=a+bx (5分) 将实验所测得不同时间1的反应物浓度c数据通过插值和中心差商法求出反应速率r。 计算lnr和lp后，用线性回归子程序计算反应级数n和反应速率常数k。 (5分) 2.程序框图(20分) 开始 输入：数据点数N, 反应温度TO,PO,步长H 时间t与delta-p的实验数据T①，PI① (1分) 8
8 P= .1000E+07Pa X= .2911 RESULT: y(CH4)= .1979 y(H2O)= .4771 y(CO)= .0813 y(H2)= .2438 （5 分） （2）第一种情况更有利于反应产物的生成。（5 分） C 1．计算原理（化学原理和计算方法）（15 分） 对于气相反应： CH3CHO CH4+CO t＝0 pA0 0 0 t＝t1 pA pA0- pA pA0- pA t1 时刻时系统中总压为： p 总＝pA＋2（pA0- pA）＝2pA0- pA＝pA0＋Δp 故 t1时刻反应物的分压为 p A = p A − Δp 0 （5 分） 化学反应速率方程： n kp t p r = − = d d 上式中， r 为化学反应速率，p 为反应物的分压，t 为时间， n 为反应级数， k 为反应速率常数，dp/dt 为反应物压力随时间的变化率。 上两边取对数： k n p t p r ) ln ln d d ln = ln(− = + 令 y = ln r, x = ln p, a = ln k,b = n,得y = a + bx （5 分） 将实验所测得不同时间 t 的反应物浓度 c 数据通过插值和中心差商法求出反应速率 r 。 计算 ln r 和 lnp 后，用线性回归子程序计算反应级数 n 和反应速率常数 k。 （5 分） 2.程序框图（20 分） 开始 输入：数据点数 N, 反应温度 TO, P0, 步长 H 时间 t 与 delta-p 的实验数据 T(I), P1(I) （I（1 分）
PIFP0-P1①(=1，N)(2分) X(I)=T(D),Y(I)=P(I)(I=I,N)( 调用中心差商子程序计算反应速率R(W(dpd)(4分) X①=HnPD1,YI=Hn-R①]=1，N)(4分) 调用一元线性回归子程序计算A,B (反应级数S=B,速率常数KS=EXP(A) (4分) 输出：A,B,SKS (1分) 结束 3.源程序(20分) PROGRAM MAIN (10分) OPEN(6.FILE-2007TESTC.TXT.STATUS-UNKNOWN) .0E-6 P0=48.4 50 CALLCF(N.H.X.Y.V) S-B SS=EXP(A) 7RD181,5XR=D1&1& END SUBROUTINE CFON.HXYZ) (5分)
9 3.源程序（20 分） PROGRAM MAIN （10 分） DIMENSION X(5),Y(5),V(5) DOUBLE PRECISION X,Y,V,A,B,R,S,SS,P0,H OPEN(6,FILE='2007TESTC.TXT',STATUS='UNKNOWN') N=5 H=1.0E-6 P0=48.4 DATA X/42.0,105.,242.,840.,1440./ DATA Y/4.53,9.86,17.86,32.53,37.86/ DO 50 I=1,N Y(I)=P0-Y(I) 50 CONTINUE CALL CF(N,H,X,Y,V) DO 100 I=1,N X(I)=LOG(Y(I)) Y(I)=LOG(ABS(V(I))) 100 CONTINUE CALL PK(N,X,Y,A,B,R) S=B SS=EXP(A) WRITE(6,11)A,B,R,S,SS 11 FORMAT(/5X,'A=',D18.10//5X,'B=',D18.10, 5X,'R=',D18.10// & 5X,'n=',D18.10,5X,'k=',D18.10/) END SUBROUTINE CF(N,H,X,Y,Z) （5 分） P(I)=P0-P1(I) （I=1,N） （2 分） 调用中心差商子程序计算反应速率 R(W)（dp/dt） （4 分） X(I)=ln[P(I)],Y(I)=ln[-R(I)] (I=1,N) （4 分） 调用一元线性回归子程序计算 A，B （反应级数 S=B，速率常数 KS=EXP(A)） （4 分） 输出： A,B, S,KS （1 分） 结束 X(I)=T(I),Y(I)=P(I) (I=1,N) （4 分）
DOUBLE PRECISIONHT.YA.YB GRG2CX YNTYA) DO 10L-N-I CALLLGRG2(X,Y.N.T.YA) GRG2(X.Y.N.T.YB) (YB-YAY 10 CALLLGRG2(X,Y.N,TYA) (2分) DIMENSION X(N).Y(N) SXX=0.0 KV-00 D010=1,N SX=SX+XI +X1*X YY+Y 10 CONTINUE DXX=SXX-SX*SX/N B-DXY/DXX R=DXY/SORT(DXX*DYY) ENDETURN (3分) DOUBLE PRECISION X.Y.TZ.S 59N N END IF IF(N.EQ RETURN EXG)-YG)(EX() 0
10 DIMENSION X(N),Y(N),Z(N) DOUBLE PRECISION X,Y,Z,H,T,YA,YB T=X(1)+H CALL LGRG2(X,Y,N,T,YA) T=X(1)+2*H CALL LGRG2(X,Y,N,T,YB) Z(1)=(4*YA-3*Y(1)-YB)/(2*H) DO 10 I=2,N-1 T=X(I)-H CALL LGRG2(X,Y,N,T,YA) T=X(I)+H CALL LGRG2(X,Y,N,T,YB) Z(I)=(YB-YA)/(2*H) 10 CONTINUE T=X(N)-H CALL LGRG2(X,Y,N,T,YA) T=X(N)-2*H CALL LGRG2(X,Y,N,T,YB) Z(N)=(3*Y(N)-4*YA+YB)/(2*H) RETURN END SUBROUTINE PK(N,X,Y,A,B,R) （2 分） DOUBLE PRECISION X,Y,A,B,R DIMENSION X(N),Y(N) SX=0.0 SY=0.0 SXX=0.0 SYY=0.0 SXY=0.0 DO 10 I=1,N X1=X(I) Y1=Y(I) SX=SX+X1 SY=SY+Y1 SXX=SXX+X1*X1 SYY=SYY+Y1*Y1 SXY=SXY+X1*Y1 10 CONTINUE DXX=SXX-SX*SX/N DYY=SYY-SY*SY/N DXY=SXY-SX*SY/N B=DXY/DXX A=(SY-B*SX)/N R=DXY/SQRT(DXX*DYY) RETURN END SUBROUTINE LGRG2(X,Y,N,T,Z) （3 分） DIMENSION X(N),Y(N) DOUBLE PRECISION X,Y,T,Z,S Z=0.0 IF(N.LE.0)RETURN IF(N.EQ.1)THEN Z=Y(1) RETURN END IF IF(N.EQ.2)THEN Z=(Y(1)*(T-X(2))-Y(2)*(T-X(1)))/(X(1)-X(2)) RETURN