《正交配置求解问题.doc》由会员分享，可在线阅读，更多相关《正交配置求解问题.doc（22页珍藏版）》请在金锄头文库上搜索。
1、正交配置求解问题：运用正交配置法求解有轴向扩散的固定床反应器中催化反应的温度和浓度分布。柱形固体床反应器中催化反应的温度和浓度方程为：+ =+r=1 =|r=0=0 -|r=1=BiwT(1,z)-Tw(z), -|r=1=0T(r,0)=T0 , c (r,0)=c0+ =+r=1 =|r=0=0 -|r=1=BiwT(1,z)-Tw(z), -|r=1=0T(r,0)=T0 , c (r,0)=c0其中R（c，T）为催化反应的速率方程，其形式为 R(c,T)=+ =+r=1 =|r=0=0 -|r=1=BiwT(1,z)-Tw(z), -|r=1=0T(r,0)=T0 , c (r,0)=
2、c0其中R（c，T）为催化反应的速率方程，其形式为 R(c,T)=解题思路：应用对称的正交配置法，有下面的方程和初始条件：=+ (1- ) = + (1-) Tj(0)=T0 , cj(0)=c0边界条件为：-AN+1,iTi=Biw(TN+1-Tw) , AN+1,ici=0将温度和浓度的边界条件代入微分方程，消去边界值，可得2N个常微分方程，而将两边界条件的代数方程同2N个常微分方程组联合，就组成2N+2个微分代数方程组。结合正交配置系数的计算程序与常微分方程组或微分方程组求解程序，可得到反应器中的温度和浓度分布。具体做法如下：一、 利用对称的正交配置格式：1、对称常微分方程程序：（COL
3、LAB.FOR,DLSODE.FOR） 主程序：IMPLICIT REAL*8 (A-H,O-Z) EXTERNAL FEX, JEX DIMENSION AS(19,19),BS(19,19),Q(19,19),XS(19),WS(19) DIMENSION DIF1(19),DIF2(19),DIF3(19),ROOT(19),V1(19),V2(19) DIMENSION Y(99), ATOL(99), RWORK(10920), IWORK(120) DOUBLE PRECISION YN1,YN2 COMMON /AB/ N,AS,BS COMMON /BC/YN1,YN2C N-
4、 FOR SYMMETRIC COLLOCATION USED FOR PARTICLE AND C M- FOR ASYMMETRIC COLLOCATION USED FOR COLUMN N=7 IW=1 IS=2 CALL COLL(AS,BS,Q,XS,WS,19,N,IW,IS) NS=N+1 WRITE(*,*) * SYMMETRIC SITUATION: * WRITE(*,*) * POLYNOMIAL ROOTS * WRITE(*,*) (XS(I),I=1,NS) WRITE(*,*) WRITE(*,*) * A-MATRIX * DO 20 I=1,NS20 WR
5、ITE(*,*) (AS(I,J),J=1,NS) WRITE(*,*) WRITE(*,*) * B-MATRIX * DO 30 I=1,NS30 WRITE(*,*) (BS(I,J),J=1,NS) WRITE(*,*) WRITE(*,*) * W-MATRIX * WRITE(*,*) (WS(J),J=1,NS)CCALCULATING THE PARAMETERS OF THE PROBLEM, WHICH WILL BE USEDC FOR THE DIMENSIONLESS FORM OF AND DEFINING OF THE PROBLEM. NEQ=2*NLRW=22
6、+9*NEQ+NEQ*2LIW=20+NEQC INITIAL CONDITIONS DO 201 I=1,N Y(I) = 1.D0Y(N+I)=0.D0 201 CONTINUE YN1=1.0D0YN2=0.D0 T = 0.D0 DT = 5.D-2ITOL = 2 RTOL = 1.D-6DO 203 I=1,NEQ ATOL(I) = 1.D-6 203 CONTINUE ITASK = 1 ISTATE = 1 IOPT = 0 MF = 22 DO 240 IOUT = 1,20TOUT = DT*DFLOAT(IOUT) CALL LSODE(FEX,NEQ,Y,T,TOUT
7、,ITOL,RTOL,ATOL,ITASK,ISTATE, 1 IOPT,RWORK,LRW,IWORK,LIW,JEX,MF)OPEN(2,FILE=LW_S_ODE.OUT) WRITE (2,(Z:,F8.4) TWRITE(2,*) RWRITE (2,(10(4X,D11.5) (XS(I),I=1,N+1) WRITE (2,*) T: WRITE (2,(10(4X,D11.5) (Y(I),I=1,N),YN1 WRITE (2,*) C:C DO 205 I=1,N WRITE (2,(10(4X,D11.5) (Y(N+I),I=1,N),YN2C WRITE(2,*)C
8、205 CONTINUEC 220 FORMAT(7H AT T =,D12.4,6H Y =,3D15.7) IF (ISTATE .LT. 0) GO TO 280 240 CONTINUE WRITE(3,260)IWORK(11),IWORK(12),IWORK(13) 260 FORMAT(/12H NO. STEPS =,I4,11H NO. F-S =,I4,11H NO. J-S =,I4) STOP 280 WRITE(3,290)ISTATE 290 FORMAT(/22H ERROR HALT. ISTATE =,I3) STOP END子程序： SUBROUTINE F
9、EX (NEQ, T, Y, YP)C VARIABLES IN THE ORDER OF Y(1) TO Y(M) FOR THE COLUMN COLLOCATION C POINTS OF 2 TO M+1, Y(M+1) TO Y(M+N) FOR PARTICLE COLLOCATION POINTSC FROM 1 TO N IN COLUMN COLLOCATION POINT 1, AND Y(M+N*(M+1)+1) TO C Y(M+N*(M+2) FOR PARTICLE COLLOCATION POINTS FROM 1 TO N IN COLUMN C COLLOCA
10、TION POINT M+2. IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(19*19),YP(19*19),AS(19,19),BS(19,19)DOUBLE PRECISION YN1,YN2 COMMON /AB/ N,AS,BS COMMON /BC/ YN1,YN2T0=0.D0TT0=0.D0DO 990 I=1,NT0=T0+AS(N+1,I)*Y(I)TT0=TT0+AS(N+1,I)*Y(N+I)990CONTINUE YN1=(0.92-T0)/(1+AS(N+1,N+1)YN2=-TT0/AS(N+1,N+1) DO 99 J=1,NT1=
11、BS(J,N+1)*YN1T2=BS(J,N+1)*YN2DO 991 I=1,NT1=T1+BS(J,I)*Y(I)T2=T2+BS(J,I)*Y(I+N)991CONTINUEYP(J)=T1+0.2*(1-Y(N+J)*EXP(20-20/Y(J)YP(N+J)=T2+0.3*(1-Y(N+J)*EXP(20-20/Y(J)99CONTINUE RETURN END SUBROUTINE JEX (NEQ, T, Y, ML, MU, PD, NRPD) DOUBLE PRECISION PD, T, Y DIMENSION Y(NEQ), PD(NRPD,NEQ) RETURN END
12、输出结果：Z: .0500 R .14993D+00 .33864D+00 .51555D+00 .67294D+00 .80460D+00 .90541D+00 .97146D+00 .10000D+01 T: .10109D+01 .10104D+01 .10088D+01 .10056D+01 .10008D+01 .99565D+00 .99154D+00 .98955D+00 C: .16524D-01 .16408D-01 .16087D-01 .15522D-01 .14890D-01 .14453D-01 .14290D-01 .14878D-01Z: .1000 R .14993D+00 .33864D+00 .51555D+00 .67294D+00 .80460D+00 .90541D+00 .97146D+00 .10000D+01 T: .10215D+01 .10192D+01 .10150D+01 .10091D+01 .10025D+01 .99642D+00 .99205D+00 .99015D+00 C: .35719D-01 .34822D-01 .33344D-01 .31668D-01 .30296D-01 .29520D-01 .Z: .9500 R .14993D+00 .33864D+00 .51555D+00
