正交配置求解问题.doc

正交配置求解问题：运用正交配置法求解有轴向扩散的固定床反应器中催化反应的温度和浓度分布柱形固体床反应器中催化反应的温度和浓度方程为：＝+ =+︱r=1 =|r=0=0 -|r=1=Biw[T(1,z)-Tw(z)], -|r=1=0T(r,0)=T0 , c (r,0)=c0＝+ =+︱r=1 =|r=0=0 -|r=1=Biw[T(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=Biw[T(1,z)-Tw(z)], -|r=1=0T(r,0)=T0 , c (r,0)=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、对称常微分方程程序：（COLLAB.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--- 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 WRITE(*,*) (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)C CALCULATING THE PARAMETERS OF THE PROBLEM, WHICH WILL BE USEDC FOR THE DIMENSIONLESS FORM OF AND DEFINING OF THE PROBLEM. NEQ=2*N LRW=22+9*NEQ+NEQ**2 LIW=20+NEQC INITIAL CONDITIONS DO 201 I=1,N Y(I) = 1.D0 Y(N+I)=0.D0 201 CONTINUE YN1=1.0D0 YN2=0.D0 T = 0.D0 DT = 5.D-2 ITOL = 2 RTOL = 1.D-6 DO 203 I=1,NEQ ATOL(I) = 1.D-6 203 CONTINUE ITASK = 1 ISTATE = 1 IOPT = 0 MF = 22 DO 240 IOUT = 1,20 TOUT = DT*DFLOAT(IOUT) CALL LSODE(FEX,NEQ,Y,T,TOUT,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)') T WRITE(2,*) 'R' WRITE (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 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 FEX (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 COLLOCATION 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,YN2 T0=0.D0 TT0=0.D0 DO 990 I=1,N T0=T0+AS(N+1,I)*Y(I) TT0=TT0+AS(N+1,I)*Y(N+I)990 CONTINUE YN1=(0.92-T0)/(1+AS(N+1,N+1)) YN2=-TT0/AS(N+1,N+1) DO 99 J=1,N T1=BS(J,N+1)*YN1 T2=BS(J,N+1)*YN2 DO 991 I=1,N T1=T1+BS(J,I)*Y(I) T2=T2+BS(J,I)*Y(I+N)991 CONTINUE YP(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))99 CONTINUE RETURN END SUBROUTINE JEX (NEQ, T, Y, ML, MU, PD, NRPD) DOUBLE PRECISION PD, T, Y DIMENSION Y(NEQ), PD(NRPD,NEQ) RETURN END输出结果：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 。