正交配置求解问题.docx
《正交配置求解问题.docx》由会员分享，可在线阅读，更多相关《正交配置求解问题.docx（34页珍藏版）》请在冰豆网上搜索。
正交配置求解问题
正交配置求解问题：
运用正交配置法求解有轴向扩散的固定床反应器中催化反应的温度和浓度分布。
柱形固体床
反应器中催化反应的温度和浓度方程为：
T_
Z
1T
r——
rrr
+'R（c,T）
c
1c
r+
R（c,T）
Z
rrr
T,
c.
Ir=1一——|r=0=0
r
r
T
c,
-|r=
：
1一Biw[T（1,z）-T
w（z）],-1r=1=0
r
r
T（r,0）=T
o,c（r,0）=c
0
T
T
1
T
'R（c,T）
-r
+
Z
Z
rr
r
c
1
c
Z
r+
R（c,T）
r
rr
T
c.
1r=1一
1r=0=0
r
r
T
c,
-
|r=1=Biw[T（1,z）-T
w（Z）],-
|r=1=0
r
r
T（r,0）=To
c（r,0）=c
0
其中R（c,T）为催化反应的速率方程，其形式为R（c,T）=1ce
1T
'-r+'R（c,T）
rrr
c
1
c
—r+R（c,1）
Z
r
rr
T,
c
1
r=1=-
|r=0=0
r
r
T
|r=1=Biw[T（1,z）-Tw（z）],-—|r=1=0
r
r
T（r,O）=To,c（r,O）=c0
其中R（c,T）为催化反应的速率方程，其形式为R（c,T）=1ce
解题思路：
应用对称的正交配置法，有下面的方程和初始条件：
dT.N1
-='BT+'（1-c）eT
dZiijiij
dcjn1
j=B
dZi1
jici+'（1-cj）eT
jiij
Tj（0）=T0,c
j（0）=c0
N1
N1
边界条件为：
-
AN+1,iTi=Biw（TN+1-Tw）,
An+1,ici=0
i1i1
将温度和浓度的边界条件代入微分方程，消去边界值，可得2N个常微分方程，而将两边界
条件的代数方程同2N个常微分方程组联合，就组成2N+2个微分代数方程组。
结合正交配置系数的计算程序与常微分方程组或微分方程组求解程序，可得到反应器中的温度和浓度分布。
具体做法如下：
一、利用对称的正交配置格式：
1、对称常微分方程程序：
（COLLAB.FORQLSODE.FOR）
主程序：
IMPLICITREAL*8（A-H,O-Z）
EXTERNALFEX,JEX
DIMENSIONAS（19,19）,BS（19,19）,Q（19,19）,XS（19）,WS（19）
DIMENSIONDIF1（19）,DIF2（19）,DIF3（19）,ROOT（19）,V1（19）,V2（19）
DIMENSIONY（99）,ATOL（99）,RW0RK（10920）,IWORK（120）
DOUBLEPRECISIONYN1,YN2
COMMON/AB/N,AS,BS
COMMON/BC/YN1,YN2
CN---FORSYMMETRICCOLLOCATIONUSEDFORPARTICLEAND
CM---FORASYMMETRICCOLLOCATIONUSEDFORCOLUMNN=7
IW=1
IS=2
CALLCOLL（AS,BS,Q,XS,WS,19,N,IW,IS）
NS=N+1
WRITE（*,*）'*SYMMETRICSITUATION:
*'
WRITE（*,*）'*POLYNOMIALROOTS*'
WRITE（*,*）（XS（I）,I=1,NS）
WRITE（*,*）WRITE（*,*）'*A-MATRIX*'
DO20I=1,NS
20WRITE（*,*）（AS（I,J）,J=1,NS）WRITE（*,*）WRITE（*,*）'*B-MATRIX*'DO30I=1,NS
30WRITE（*,*）（BS（I,J）,J=1,NS）WRITE（*,*）
WRITE（*,*）'*W-MATRIX*'
WRITE（*,*）（WS（J）,J=1,NS）
BE
CCALCULATINGTHEPARAMETERSOFTHEPROBLEM,WHICHWILLUSED
CFORTHEDIMENSIONLESSFORMOFANDDEFININGOFTHEPROBLEM.
NEQ=2*N
LRW=22+9*NEQ+NEQ**2
LIW=20+NEQ
CINITIALCONDITIONSDO201I=1,NY（I）=1.D0
Y（N+I）=0.D0
201CONTINUEYN1=1.0D0YN2=0.D0
T=0.D0
DT=5.D-2
ITOL=2
RTOL=1.D-6
DO203I=1,NEQATOL（I）=1.D-6
203CONTINUEITASK=1ISTATE=1IOPT=0
MF=22DO240IOUT=1,20TOUT=DT*DFLOAT（IOUT）
CALLLSODE（FEX,NEQ,Y,T,TOUT,ITOL,RTOL,ATOL,ITASK,ISTATE,
1IOPT,RWORK,LRW,IWORK,LIW,JEX,MF）OPEN（2,FILE='LW_S_ODE.OUT'）
WRITE（2,'（''Z:
'',F8.4）'）TWRITE（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:
'CDO205I=1,N
WRITE（2,'（10（4X,D11.5））'）（Y（N+I）,I=1,N）,YN2
CWRITE（2,*）
C205CONTINUE
C220FORMAT（7HATT=,D12.4,6HY=,3D15.7）
IF（ISTATE.LT.0）GOTO280
240CONTINUEWRITE（3,260）IWORK（11）,IWORK（12）,IWORK（13）
260FORMAT（/12HNO.STEPS=,I4,11HNO.F-S=,I4,11HNO.J-S=,I4）
STOP
280WRITE（3,290）ISTATE
290FORMAT（///22HERRORHALT..ISTATE=,I3）STOPEND子程序：
SUBROUTINEFEX（NEQ,T,Y,YP）
CVARIABLESINTHEORDEROFY
（1）TOY（M）FORTHECOLUMNCOLLOCATION
CPOINTSOF2TOM+1,Y（M+1）TOY（M+N）FORPARTICLECOLLOCATIONPOINTS
CFROM1TONINCOLUMNCOLLOCATIONPOINT1,ANDY（M+N*（M+1）+1）TO
CY（M+N*（M+2））FORPARTICLECOLLOCATIONPOINTSFROM1TONINCOLUMN
CCOLLOCATIONPOINTM+2.
IMPLICITREAL*8（A-H,O-Z）
DIMENSIONY（19*19）,YP（19*19）,AS（19,19）,BS（19,19）
DOUBLEPRECISIONYN1,YN2COMMON/AB/N,AS,BSCOMMON/BC/YN1,YN2
T0=0.D0
TT0=0.D0
DO990I=1,N
T0=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）
DO99J=1,N
T1=BS（J,N+1）*YN1
T2=BS（J,N+1）*YN2
DO991I=1,N
T1=T1+BS（J,I）*Y（I）
T2=T2+BS（J,I）*Y（I+N）
991CONTINUE
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））
99CONTINUE
RETURN
END
SUBROUTINEJEX（NEQ,T,Y,ML,MU,PD,NRPD）DOUBLEPRECISIONPD,T,Y
DIMENSIONY（NEQ）,PD（NRPD,NEQ）
RETURN
END
输出结果：
Z:
.0500
R
.14993D+00
.33864D+00
.51555D+00
.67294D+00
.80460D+00
.90541D+00
T-
.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-01
Z:
.1000
D
R
.14993D+00
.33864D+00
.51555D+00
.67294D+00
.80460D+00
.90541D+00
T-
.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
.67294D+00
.80460D+00
.90541D+00
T・
.97146D+00
.10000D+01
T:
.13623D+01
.13464D+01
.13210D+01
.12905D+01
.12599D+01
.12338D+01
.12155D+01
.12046D+01
C:
.10000D+01
.10000D+01
.10000D+01
.10000D+01
.10000D+01
.10000D+01
.10000D+01
.10000D+01
Z:
1.0000
D
R
.14993D+00
.33864D+00
.51555D+00
.67294D+00
.80460D+00
.90541D+00
T-
.97146D+00
.10000D+01
T:
.13290D+01
.13142D+01
.12906D+01
.12624D+01
.12341D+01
.12099D+01
.11930D+01
.11825D+01
C:
.10000D+01
.10000D+01
.10000D+01
.10000D+01
.10000D+01
.10000D+01
.10000D+01
.10000D+01
NO.J-S=29
NO.STEPS=210NO.F-S=677
2、对称代数微分方程程序：
主程序：
PROGRAMDEMO
C
-WECOUNTFUNCTIONANDJACOBIANEVALUATIONS-
INTEGER
NF,NJ
COMMON/SCORE/NF,NJ
C
-FORPROBLEMS-
INTEGER
N,NB
PARAMETER
（NB=100）
DOUBLEPRECISION
T,H,TOUT,TREP,Y（NB）,B（NB,NB）,ATOLER,RTOLER
C
-WORKSPACEDECLARATION-
INTEGER
NWORK
PARAMETER
（NWORK=10000）
INTEGER
KOUNTI,KOUNTD,IWORK（NWORK）
DOUBLEPRECISION
DWORK（NWORK）
C
-FORBESIRK-
INTEGER
INFO,RESULT,MAXIT,ITER,I,STANDU
DOUBLEPRECISION
TREPRT,HNEXT,HMIN,HMAX,NSTEP,
+
ABSTOL（NB）,RELTOL（NB）
LOGICAL
NUMJAC,ALGEQN,STANDT
EXTERNAL
STANDU,STANDT
C
-PROBLEMROUTINE-
INTEGERREPT
EXTERNALFUNC,REPT
INTEGERM,NS
DOUBLEPRECISIONAS（19,19）,BS（19,19）,XS（19）
COMMON/AB/M,NS,AS,BS,XS
C-STEPSIZEISZERO->BESIRKFINDSIT-
H=0.0D0
C-INITIALIZEPROBLEM:
-
CALLLW（）
WRITE（*,*）'MAIN'
WRITE（*,*）（AS（1,J）,J=1,NS）,M
CALLINIT（N,T,H,TOUT,TREP,Y,B,NB,NUMJAC,ATOLER,RTOLER）C-INITIALIZECOUNTERS-
NF=0
NJ=0
C-SETINFOLEVEL-
INFO=0
C-SETNEXTREPORTTIME-
TREPRT=TREP
C-SETMINIMUMANDMAXIMUMSTEPSIZE-
HMIN=1.0D-10
HMAX=TOUT
C-SETNUMERICALJACOBIANSTEP-
NSTEP=1.0D-3
C-MAXIMUMNUMBEROFITERATIONS-
MAXIT=500
C-SETTOLERANCES（FOREACHVARIABLE）-
DOI=1,NABSTOL（I）=ATOLERRELTOL（I）=RTOLER
ENDDO
ALGEQN=.FALSE.
C-SETWORKSPACECOUNTERS-
KOUNTI=NWORKKOUNTD=NWORK
C-SOLVE-
CALLINIT3
CALL
BESIRK（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,KOUNTD,KOUNTI）
C-REPORT-
WRITE（*,*）
WRITE（*,*）'BESIRKRETURNS=',RESULT
WRITE（*,*）'FUNCTIONEVALUATIONS=',NF
WRITE（*,*）'JACOBIANEVALUATIONS=',NJ
WRITE（*,*）
C-TERMINATE-
STOP
END
子程序：
SUBROUTINELW（）
IMPLICITREAL*8（A-H,O-Z）
CEXTERNALFEX,JEX
DIMENSIONAS（19,19）,BS（19,19）,Q（19,19）,XS（19）,WS（19）
CDIMENSIONAU（19,19）,BU（19,19）,XA（19）,WA（19）
DIMENSIONDIF1（19）,DIF2（19）,DIF3（19）,ROOT（19）,V1（19）,V2（19）
DIMENSIONY（99）,ATOL（99）,RWORK（10920）,IWORK（120）
COMMON/AB/N,NS,AS,BS,XS
CCOMMON/BC/Y1,YM2,YN1
CCOMMON/PARA/PSAI,SITA,PE,SAI,PAK
CN---FORSYMMETRICCOLLOCATIONUSEDFORPARTICLEAND
CM---FORASYMMETRICCOLLOCATIONUSEDFORCOLUMN
N=7
IW=1
IS=2
CALLCOLL（AS,BS,Q,XS,WS,19,N,IW,IS）
NS=N+1
WRITE（*,*）'*SYMMETRICSITUATION:
*'
WRITE（*,*）'*POLYNOMIALROOTS*'
WRITE（*,*）（XS（I）,I=1,NS）
WRITE（*,*）
WRITE（*,*）'*A-MATRIX*'
DO20I=1,NS
20WRITE（*,*）（AS（I,J）,J=1,NS）
WRITE（*,*）
WRITE（*,*）'*B-MATRIX*'
DO30I=1,NS
30WRITE（*,*）（BS（I,J）,J=1,NS）
WRITE（*,*）
WRITE（*,*）'*W-MATRIX*'
WRITE（*,*）（WS（J）,J=1,NS）
WRITE（*,*）'ORXO'
WRITE（*,*）（AS（1,I）,I=1,NS）,N
RETURN
END
SUBROUTINEFUNC（N,INFO,F,Y,DFDY,T,CALCJ,ERROR,+DW,IW,KOUNTD,KOUNTI）
INTEGERN,INFO,ERROR,KOUNTD,KOUNTI,IW（*）
DOUBLEPRECISIONF（N）,Y（N）,DFDY（N,N）,T,DW（*）LOGICALCALCJ
INTEGERNF,NJ
COMMON/SCORE/NF,NJ
INTEGERM,NS
DOUBLEPRECISIONAS（19,19）,BS（19,19）,XS（19）
COMMON/AB/M,NS,AS,BS,XS
NF=NF+1
CWRITE（*,*）'FUNC'
CWRITE（*,*）（AS（1,J）,J=1,NS）,MM
CPAUSE
DO99J=1,M
T1=0.0D0
DO991I=1,NS
T1=T1+BS（J,I）*Y（I）
991CONTINUE
T2=0.D0
DO992L=1,NS
T2=T2+BS（J,L）*Y（NS+L）
992CONTINUE
F（J）=T1+0.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））
99CONTINUE
T0=0.D0
TT0=0.D0
DO990I=1,NS
T0=T0+AS（NS,I）*Y（I）
TT0=TT0+AS（NS,I）*Y（NS+I）
990CONTINUE
F（NS）=（Y（NS）-0.92）+T0
F（2*NS）=TT0
RETURN
END
SUBROUTINEINIT（N,T,H,TOUT,TREP,Y,B,ND,NUMJAC,ATOLER,RTOLER）
INTEGER
DOUBLEPRECISION
N,ND
T,H,TOUT,TREP,Y（ND）,B（ND,ND）,ATOLER,RTOLER
LOGICAL
INTEGER
DOUBLEPRECISION
NUMJAC
I
ZERO,ONE
PARAMETER
（ZERO=0.0D0,ONE=1.0D0）
INTEGERM,NS
DOUBLEPRECISIONAS（19,19）,BS（19,19）,XS（19）
COMMON/AB/M,NS,AS,BS,XS
N=2*NS
T=ZERO
TOUT=1.0D0
TREP=0.1D0
DOI=1,NSY（I）=1.0D0Y（NS+I）=0.0D0ENDDOCALLMAT0（B,ND,N）DOI=1,MB（I,I）=ONEB（NS+I,NS+I）=ONE
ENDDO
NUMJAC=.TRUE.ATOLER=1.0D-6RTOLER=0.0D0
RETURN
END
INTEGERFUNCTIONREPT（N,ITER,INFO,T,TREPRT,Y,E）
INTEGERN,ITER,INFO
DOUBLEPRECISIONT,Y（N）,E（N）,TREPRT
INTEGERM,NS
DOUBLEPRECISIONAS（19,19）,BS（19,19）,XS（19）
COMMON/AB/M,NS,AS,BS,XS
OPEN（2,FILE='JJM_S_DAE.OUT'）
WRITE（3,*）'Z=',T
WRITE（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,*）'T'
WRITE（3,'（9（4X,D11.5））'）
WRITE（3,*）'C'
WRITE（3,'（9（4X,D11.5））'）
REPT=0
RETURN
END
C
升级会员 