正交配置求解问题.docx

1、正交配置求解问题正交配置求解问题：运用正交配置法求解有轴向扩散的固定床反应器中催化反应的温度和浓度分布。 柱形固体床反应器中催化反应的温度和浓度方程为：T _Z,1 T rr r r+ R(c,T)c1 cr +R(c,T)Zr r rT ,c.I r=1 一 | r=0=0rrTc ,- | r=：1一BiwT(1,z)-Tw(z), - 1 r=1=0rrT(r,0)=To , c (r,0)=c0TT,1TR(c,T)- r+ZZr rrc1cZr +R(c,T)rr rTc.1 r=1 一 1 r=0=0rrTc ,-| r=1=BiwT(1,z)-Tw(Z),-| r=1=0rrT(

2、r,0)=T o,c (r,0)=c0其中R( c, T)为催化反应的速率方程，其形式为 R(c,T)= 1 c e1 T- r + R(c,T)r r rc1cr + R(c, 1 )Zrr rT ,c1r=1 =-| r=0=0rrT| r=1=BiwT(1,z)-T w(z), - | r=1 =0rrT(r,O)=T o , c (r,O)=c 0其中R( c, T)为催化反应的速率方程，其形式为 R(c,T)= 1 c e解题思路：应用对称的正交配置法，有下面的方程和初始条件：dT . N 1-= B T + (1- c ) e TdZ i i ji i jdc j n 1j= Bd

3、Z i 1jici + (1-cj)e Tji i jT j(0)=T 0 , cj(0)=c 0N 1N 1边界条件为：-AN+1 , iTi=Bi w(T N+1 -T w),An+1 ,ici=0i 1 i 1将温度和浓度的边界条件代入微分方程，消去边界值，可得 2N个常微分方程，而将两边界条件的代数方程同 2N个常微分方程组联合，就组成2N+2个微分代数方程组。结合正交配置 系数的计算程序与常微分方程组或微分方程组求解程序，可得到反应器中的温度和浓度分 布。具体做法如下：一、利用对称的正交配置格式：1、对称常微分方程程序： (COLLAB.FORQLSODE.FOR )主程序：IMPL

4、ICIT REAL*8 (A-H,O-Z)EXTERNAL FEX, JEXDIMENSION 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), RW0RK(10920), IWORK(120)DOUBLE PRECISION YN1,YN2COMMON /AB/ N,AS,BSCOMMON /BC/YN1,YN2C N- FOR SYMMETRIC COLLOCATION USED FOR

5、 PARTICLE ANDC M- FOR ASYMMETRIC COLLOCATION USED FOR COLUMN N=7IW=1IS=2CALL COLL(AS,BS,Q,XS,WS,19,N,IW,IS)NS=N+1WRITE(*,*) * SYMMETRIC SITUATION: *WRITE(*,*) * POL YNOMIAL 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

6、(*,*) * B-MATRIX * DO 30 I=1,NS30 WRITE(*,*) (BS(I,J),J=1,NS) WRITE(*,*)WRITE(*,*) * W-MATRIX * WRITE(*,*) (WS(J),J=1,NS)BEC CALCULATING THE PARAMETERS OF THE PROBLEM, WHICH WILL USEDC FOR THE DIMENSIONLESS FORM OF AND DEFINING OF THE PROBLEM.NEQ=2*NLRW=22+9*NEQ+NEQ*2LIW=20+NEQC INITIAL CONDITIONS D

7、O 201 I=1,N Y(I) = 1.D0Y(N+I)=0.D0201 CONTINUE YN1=1.0D0 YN2=0.D0T = 0.D0DT = 5.D-2ITOL = 2RTOL = 1.D-6DO 203 I=1,NEQ ATOL(I) = 1.D-6203 CONTINUE ITASK = 1 ISTATE = 1 IOPT = 0MF = 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,IWO

8、RK,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),YN1WRITE (2,*) C: C DO 205 I=1,NWRITE (2,(10(4X,D11.5) (Y(N+I),I=1,N),YN2C WRITE(2,*)C 205 CONTINUEC 220 FORMA T(7H AT T =,D12.4,6H Y =,3D

9、15.7)IF (ISTATE .LT. 0) GO TO 280240 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)STOP280 WRITE(3,290)ISTA TE290 FORMAT(/22H ERROR HALT. ISTA TE =,I3) STOP END 子程序：SUBROUTINE FEX (NEQ, T, Y , YP)C VARIABLES IN THE ORDER OF Y(1) TO

10、 Y(M) FOR THE COLUMN COLLOCATIONC 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) TOC Y(M+N*(M+2) FOR PARTICLE COLLOCATION POINTS FROM 1 TO N IN COLUMNC COLLOCATION POINT M+2.IMPLICIT REAL*8 (A-H,O-Z)DIMENSION Y(19*1

11、9),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)990CONTINUEYN1=(0.92-T0)/(1+AS(N+1,N+1)YN2=-TT0/AS(N+1,N+1)DO 99 J=1,NT1=BS(J,N+1)*YN1T2=BS(J,N+1)*YN2DO 991 I=1,NT1=T1+BS(J,I)*Y(I)T

12、2=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)99 CONTINUERETURNENDSUBROUTINE JEX (NEQ, T, Y , ML, MU, PD, NRPD) DOUBLE PRECISION PD, T, YDIMENSION Y(NEQ), PD(NRPD,NEQ)RETURNEND输出结果：Z: .0500R.14993D+00.33864D+00.51555D+00.67294D+00.80460D+00

13、.90541D+00T-.97146D+00.10000D+01T:.10109D+01.10104D+01.10088D+01.10056D+01.10008D+01.99565D+00.99154D+00.98955D+00C:.16524D-01.16408D-01.16087D-01.15522D-01.14890D-01.14453D-01.14290D-01.14878D-01Z: .1000DR.14993D+00.33864D+00.51555D+00.67294D+00.80460D+00.90541D+00T-.97146D+00.10000D+01T:.10215D+01

14、.10192D+01.10150D+01.10091D+01.10025D+01.99642D+00.99205D+00.99015D+00C:.35719D-01.34822D-01.33344D-01.31668D-01.30296D-01.29520D-01Z: .9500R.14993D+00.33864D+00.51555D+00.67294D+00.80460D+00.90541D+00T.97146D+00.10000D+01T:.13623D+01.13464D+01.13210D+01.12905D+01.12599D+01.12338D+01.12155D+01.12046

15、D+01C:.10000D+01.10000D+01.10000D+01.10000D+01.10000D+01.10000D+01.10000D+01.10000D+01Z: 1.0000DR.14993D+00.33864D+00.51555D+00.67294D+00.80460D+00.90541D+00T-.97146D+00.10000D+01T:.13290D+01.13142D+01.12906D+01.12624D+01.12341D+01.12099D+01.11930D+01.11825D+01C:.10000D+01.10000D+01.10000D+01.10000D

16、+01.10000D+01.10000D+01.10000D+01.10000D+01NO. J-S = 29NO. STEPS = 210 NO. F-S = 6772、对称代数微分方程程序： 主程序： PROGRAM DEMOC- WE COUNT FUNCTION AND JACOBIAN EV ALUATIONS -INTEGERNF,NJCOMMON /SCORE/ NF,NJC- FOR PROBLEMS -INTEGERN,NBPARAMETER(NB=100)DOUBLE PRECISIONT,H,TOUT,TREP,Y(NB),B(NB,NB),ATOLER,RTOLERC-

17、 WORKSPACE DECLARA TION -INTEGERNWORKPARAMETER(NWORK=10000)INTEGERKOUNTI,KOUNTD,IWORK(NWORK)DOUBLE PRECISIONDWORK(NWORK)C- FOR BESIRK -INTEGERINFO,RESULT,MAXIT,ITER,I,STANDUDOUBLE PRECISIONTREPRT,HNEXT,HMIN,HMAX,NSTEP,+ABSTOL(NB),RELTOL(NB)LOGICALNUMJAC,ALGEQN,STANDTEXTERNALSTANDU,STANDTC- PROBLEM R

18、OUTINE -INTEGER REPTEXTERNAL FUNC,REPTINTEGER M,NSDOUBLE PRECISION AS(19,19),BS(19,19),XS(19)COMMON /AB/ M,NS,AS,BS,XSC - STEPSIZE IS ZERO - BESIRK FINDS IT -H = 0.0D0C - INITIALIZE PROBLEM: -CALLLW()WRITE(*,*) MAINWRITE(*,*) (AS(1,J),J=1,NS),MCALL INIT(N,T,H,TOUT,TREP,Y ,B,NB,NUMJAC,ATOLER,RTOLER)

19、C - INITIALIZE COUNTERS -NF = 0NJ = 0C - SET INFO LEVEL -INFO = 0C - SET NEXT REPORT TIME -TREPRT = TREPC - SET MINIMUM AND MAXIMUM STEPSIZE -HMIN = 1.0D-10HMAX = TOUTC - SET NUMERICAL JACOBIAN STEP -NSTEP = 1.0D-3C - MAXIMUM NUMBER OF ITERA TIONS -MAXIT = 500C - SET TOLERANCES (FOR EACH V ARIABLE)

20、-DO I = 1, N ABSTOL(I) = ATOLER RELTOL(I) = RTOLERENDDOALGEQN=.FALSE.C - SET WORKSPACE COUNTERS -KOUNTI = NWORK KOUNTD = NWORKC - SOLVE -CALL INIT3CALLBESIRK(N,INFO,Y ,B,NB,T,TOUT,TREPRT,H,HNEXT,HMIN,HMAX,RESULT,+ NSTEP,FUNC,STANDU,STANDT,REPT,NUMJAC,MAXIT,ITER, +ABSTOL,RELTOL,ALGEQN,DWORK,IWORK,KOU

21、NTD,KOUNTI)C - REPORT -WRITE(*,*)WRITE(*,*) BESIRK RETURNS =,RESULTWRITE(*,*) FUNCTION EV ALUATIONS=,NFWRITE(*,*) JACOBIAN EV ALUATIONS=,NJWRITE(*,*)C - TERMINATE -STOPEND子程序：SUBROUTINELW()IMPLICIT REAL*8 (A-H,O-Z)C EXTERNAL FEX, JEXDIMENSION AS(19,19),BS(19,19),Q(19,19),XS(19),WS(19)C DIMENSION AU(

22、19,19),BU(19,19),XA(19),WA(19)DIMENSION DIF1(19),DIF2(19),DIF3(19),ROOT(19),V1(19),V2(19)DIMENSION Y(99), ATOL(99), RWORK(10920), IWORK(120)COMMON /AB/ N,NS,AS,BS,XSC COMMON /BC/ Y1,YM2,YN1C COMMON /PARA/ PSAI,SITA,PE,SAI,PAKC N- FOR SYMMETRIC COLLOCATION USED FOR PARTICLE ANDC M- FOR ASYMMETRIC COL

23、LOCATION USED FOR COLUMNN=7IW=1IS=2CALL COLL(AS,BS,Q,XS,WS,19,N,IW,IS)NS=N+1WRITE(*,*) * SYMMETRIC SITUATION: *WRITE(*,*) * POL YNOMIAL 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(

24、*,*) (BS(I,J),J=1,NS)WRITE(*,*)WRITE(*,*) * W-MATRIX * WRITE(*,*) (WS(J),J=1,NS)WRITE(*,*) ORXOWRITE(*,*)(AS(1,I),I=1,NS),NRETURNENDSUBROUTINE FUNC(N,INFO,F,Y ,DFDY ,T,CALCJ,ERROR, + DW,IW,KOUNTD,KOUNTI)INTEGER N,INFO,ERROR,KOUNTD,KOUNTI,IW(*)DOUBLE PRECISION F(N),Y(N),DFDY(N,N),T,DW(*) LOGICAL CALC

25、JINTEGER NF,NJCOMMON /SCORE/ NF,NJINTEGER M,NSDOUBLE PRECISION AS(19,19),BS(19,19),XS(19)COMMON /AB/ M,NS,AS,BS,XSNF=NF+1C WRITE(*,*) FUNCC WRITE(*,*) (AS(1,J),J=1,NS),MMC PAUSEDO 99 J=1,MT1=0.0D0DO 991 I=1,NST1=T1+BS(J,I)*Y(I)991 CONTINUET2=0.D0DO 992 L=1,NST2=T2+BS(J,L)*Y(NS+L)992CONTINUEF(J)=T1+0

26、.2*(1-Y(NS+J)*EXP(20-20/Y(J)F(NS+J)=T2+0.3*(1-Y(NS+J)*EXP(20-20/Y(J)99 CONTINUET0=0.D0TT0=0.D0DO 990 I=1,NST0=T0+AS(NS,I)*Y(I)TT0=TT0+AS(NS,I)*Y(NS+I)990 CONTINUEF(NS)=(Y(NS)-0.92)+T0F(2*NS)=TT0RETURNENDSUBROUTINE INIT(N,T,H,TOUT,TREP,Y,B,ND,NUMJAC,ATOLER,RTOLER)INTEGERDOUBLE PRECISIONN,NDT,H,TOUT,T

27、REP,Y(ND),B(ND,ND),ATOLER,RTOLERLOGICALINTEGERDOUBLE PRECISIONNUMJACIZERO,ONEPARAMETER(ZERO=0.0D0,ONE=1.0D0)INTEGER M,NSDOUBLE PRECISION AS(19,19),BS(19,19),XS(19)COMMON /AB/ M,NS,AS,BS,XSN = 2*NST = ZEROTOUT = 1.0D0TREP =0.1D0DO I=1,NS Y(I)=1.0D0 Y(NS+I)=0.0D0 ENDDO CALL MAT0(B,ND,N) DO I=1,M B(I,I

28、) = ONE B(NS+I,NS+I)=ONEENDDONUMJAC = .TRUE. ATOLER = 1.0D-6 RTOLER = 0.0D0RETURNENDINTEGER FUNCTION REPT(N,ITER,INFO,T,TREPRT,Y,E)INTEGER N,ITER,INFODOUBLE PRECISION T,Y(N),E(N),TREPRTINTEGER M,NSDOUBLE PRECISION AS(19,19),BS(19,19),XS(19)COMMON /AB/ M,NS,AS,BS,XSOPEN(2,FILE=JJM_S_DAE.OUT)WRITE(3,*)Z=,TWRITE(3,*) R(XS(I),I=1,NS)(Y(I),I=1,NS)(Y(NS+I),I=1,NS)WRITE(3,(9(4X,D11.5)WRITE(3,*) TWRITE(3,(9(4X,D11.5)WRITE(3,*) CWRITE(3,(9(4X,D11.5)REPT=0RETURNENDC