aspen RPLUG 的fortran动力学子程序纠错

kele1350 · 发表于 1970-1-1 08:00:00

这是我仿照论坛上的示例修改的子程序，请高手帮忙看一下，从语法上看看有没什么错误，多谢！编译后无法运行。
示例文件在这里http://meng.horse/forum.php?mod=viewthread&tid=76210

C    User Kinetics Subroutine for Urea synthensis, 2NH3+CO2=CARB; CARB=UREA+H2O
C
   SUBROUTINE USURA (SOUT, NSUBS,  IDXSUB, ITYPE,  NINT,
   2                INT, NREAL,  REAL,    IDS, NPO,
   3                NBOPST, NIWORK, IWORK, NWORK,  WORK,
   4                NC,    NR,    STOIC, RATES,  FLUXM,
   5                FLUXS,  XCURR,  NTCAT, RATCAT, NTSSAT,
   6                RATSSA, KCALL,  KFAIL, KFLASH, NCOMP,
   7                IDX, Y,    X,       X1,    X2,
   8                NRALL,  RATALL, NUSERV, USERV,  NINTR,
   9                INTR, NREALR, REALR, NIWR, IWR,
   1                NWR, WR)
C
   IMPLICIT NONE
C
C    DECLARE VARIABLES USED IN DIMENSIONING
C
   INTEGER NSUBS, NINT,  NPO, NIWORK,NWORK,
   +       NC, NR, NTCAT, NTSSAT,NCOMP,
   +       NRALL, NUSERV,NINTR, NREALR,NIWR,
   +       NWR
C
#include "pputl_ppglob.cmn"
#include "ppexec_user.cmn"
   EQUIVALENCE (RMISS, USER_RUMISS)
   EQUIVALENCE (IMISS, USER_IUMISS)
C
C
C
C
C.....RCSTR...
#include "rcst_rcstri.cmn"
#include "rxn_rcstrr.cmn"
C
C.....RPLUG...
#include "rplg_rplugi.cmn"
#include "rplg_rplugr.cmn"
   EQUIVALENCE (XLEN, RPLUGR_UXLONG)
   EQUIVALENCE (DIAM, RPLUGR_UDIAM)
C
C.....RBATCH...
#include "rbtc_rbati.cmn"
#include "rbtc_rbatr.cmn"
C
C.....PRES-RELIEF...
#include "rbtc_presrr.cmn"
C
C.....REACTOR (OR PRES-RELIEF VESSEL) PROPERTIES...
#include "rxn_rprops.cmn"
   EQUIVALENCE (TEMP, RPROPS_UTEMP)
   EQUIVALENCE (PRES, RPROPS_UPRES)
   EQUIVALENCE (VFRAC, RPROPS_UVFRAC)
   EQUIVALENCE (BETA, RPROPS_UBETA)
   EQUIVALENCE (VVAP, RPROPS_UVVAP)
   EQUIVALENCE (VLIQ, RPROPS_UVLIQ)
   EQUIVALENCE (VLIQS, RPROPS_UVLIQS)
C
#include "shs_stwork.cmn"
   EQUIVALENCE (MKBAS, STWORK_NDUM)
   EQUIVALENCE (MKPHAS, STWORK_NBLM)
   EQUIVALENCE (MTAPP, STWORK_NCOVAR)
   EQUIVALENCE (MKBASS, STWORK_NWR)
   EQUIVALENCE (MTAPPS, STWORK_NIWR)
   EQUIVALENCE (SSALT, STWORK_RDUM1)
   EQUIVALENCE (VSALT, STWORK_RDUM2)
   EQUIVALENCE (FSALT, STWORK_FFSALT)
C
#include "dms_ncomp.cmn"
#include "dms_plex.cmn"
   EQUIVALENCE (IB(1), B(1))
C
C    DECLARE ARGUMENTS
C
   INTEGER IDXSUB(NSUBS),ITYPE(NSUBS), INT(NINT),
   +       IDS(2),NBOPST(6,NPO),IWORK(NIWORK),
   +       IDX(NCOMP), INTR(NINTR),  IWR(NIWR),
   +       NREAL, KCALL, KFAIL, KFLASH,I
   REAL*8 SOUT(1),    WORK(NWORK),
   +    STOIC(NC,NSUBS,NR),  RATES(5),
   +    FLUXM(1),    FLUXS(1),    RATCAT(NTCAT),
   +    RATSSA(NTSSAT),    Y(NCOMP),
   +    X(NCOMP),    X1(NCOMP), X2(NCOMP)
   REAL*8 RATALL(NRALL),USERV(NUSERV),
   +    REALR(NREALR),WR(NWR),    XCURR
C
C
C    DECLARE LOCAL VARIABLES
C
   INTEGER IMISS, MKBAS, MKPHAS,MTAPP, MKBASS,
   +       MTAPPS,LMW
   REAL*8 REAL(NREAL),  XL(7), B(1),  RMISS, XLEN,
   +    DIAM,  TEMP,  PRES,  VFRAC, BETA,
   +    VVAP,  VLIQ,  VLIQS, SSALT, VSALT,
   +    FSALT, RTEMP, DENLIQ,KKA, KKB,
   +    AREA, RTA, RTB, SUMR
   INTEGER XMW
#include "dms_ipoff1.cmn"
C
C    INITIALIZE RATES
C
C
C  STATEMENT FUNCTIONS FOLLOW
C
C
C    BEGIN EXECUTABLE CODE
C
   XMW(I) = LMW + I
C
C  SET PLEX OFFSETS
C
   LMW = IPOFF1_IPOFF1(306)
   DO 100 I = 1, NC
      RATES(I) = 0D0
      XL(I) = 0D0
100  CONTINUE
C
   DO 101 I = 1, NCOMP
      XL(IDX(I)) = X(I)
101  CONTINUE
C
C    The structure in the array SOUT is as follows:
C
C    SOUT(1) - SOUT(NCC) : Component flowrates(kg-moles/sec)
C    SOUT(NCC+1)       : Total flowrates(kg-moles/sec)
C    SOUT(NCC+2)       : Temperature(K)
C    SOUT(NCC+3)       : Pressure(N/SQM)
C    SOUT(NCC+4)       : Mass enthalpy(J/KG)
C    SOUT(NCC+5)       : Molar vapor fraction
C    SOUT(NCC+6)       : Molar liquid fraction
C    SOUT(NCC+7)       : Mass entropy(J/KG-K)
C    SOUT(NCC+8)       : Mass density(KG/CUM)
C    SOUT(NCC+9)       : Molecular Weight
C
C    Details of this infomation can be found in the User Guide
C    on page 17-10 for 8.5-6 or Reference Manual Vol. 6, User Models
C    for Release 9.
C
C    Set Reactor Temperature
C
   RTEMP = SOUT(NCOMP_NCC+2)
C
C    Get liquid mixture density
C
C    DENLIQ = 1.0 / STWORK_VL
C
C    Compute Concentration for reactants:
C
C    CONC(i) = X(i) * RHO
C
C    CONAA = XL(1) * DENLIQ
C    CONACT = XL(2) * DENLIQ
C
C
C
C    CALCULATE RATES OF REACTIONS
C
C    The rates must be set for every single component exists in the
C    system(both reactants and products) according to stoichiometry.
C
C    RATES - The rates for all components
C
C    Acetone + Allyl Alcohol -> N-Propyl Propionate
C    XL1=UREA,XL2=CARB,XL3=NH3,XL4=CO2,XL5=H2O
C    FTERM = REALR(1) * EXP( -REALR(2) / PPGLOB_RGAS / RTEMP )
C    STERM = conaa * conact**0.5
C    RATE  = FTERM * STERM * RCSTRR_VOLRC
C    RATES(1) = -rate
C    RATES(2) = -rate
C    RATES(3) = rate
C    AREA=DFLOAT(RPLUGI_NTUBE)*3.141592654/4.*DIAM**2
   KKA=XL(2)/XL(4)/XL(4)/XL(3)
   KKB=XL(1)*XL(5)/XL(2)
   RTA=1.0e12*EXP(-1E07/PPGLOB_RGAS/RTEMP)*(XL(4)**2*XL(3)-XL(2)/KKA)
   RTB=1.5e09*EXP(-1E08/PPGLOB_RGAS/RTEMP)*(XL(2)-XL(1)*XL(5)/KKB)
   RATES(1)=AREA*RTB
   RATES(2)=AREA*(RTA-RTB)
   RATES(3)=AREA*RTA
   RATES(4)=-2*RTA*AREA
   RATES(5)=RTB*AREA
C    Checking Mass Balance for RATES Vector
C
   SUMR = 1.0
   DO 102 I=1,NCOMP_NCC
      SUMR = SUMR + RATES(I)*B(XMW(I))
102  CONTINUE
   IF (DABS(1.0 - SUMR) .GT. 1E-6) THEN
      WRITE(USER_NHSTRY,*)
   +       'Calculated Rates Not in Mass Balance!'
      WRITE(USER_NHSTRY,*)
   +       'Errors in User Kinetics Subroutines!'
   ENDIF
   write(USER_NHSTRY,*) 'volrc = ', RCSTRR_VOLRC
   write(USER_NHSTRY,*) 'vfrrc = ', RCSTRR_VFRRC
   RETURN
#undef P_NPOFF1
   END

LeeMoon0816 · 发表于 1970-1-1 08:00:00

本帖最后由 LeeMoon0816 于 2020-3-16 11:19 编辑

根据方程式：（A）2NH3+CO2=CARB;（B） CARB=UREA+H2O 和 XL1=UREA,XL2=CARB,XL3=NH3,XL4=CO2,XL5=H2O
来看RATES(3)=AREA*RTA 和RATES(4)=-2*RTA*AREA是错误的，应为RATES(3)=-2*RTA*AREA、RATES(4)=-AREA*RTA。
此外需要提醒的是XL1=UREA,XL2=CARB,XL3=NH3,XL4=CO2,XL5=H2O 还需要与ASPEN中的component ID 对应上，也是就components -specifications| component ID列表中从上到下的顺序应为UREA、CARB、NH3、CO2、H2O，否则物料不守恒就无法运行

好久不见H · 发表于 1970-1-1 08:00:00

LeeMoon0816 发表于 2020-3-16 11:17
根据方程式：（A）2NH3+CO2=CARB;（B） CARB=UREA+H2O 和 XL1=UREA,XL2=CARB,XL3=NH3,XL4=CO2,XL5=H2O
来 ...

您好，请教一下，这个XL具体代表的是反应物的什么呢？

aspen RPLUG 的fortran动力学子程序 纠错

相关帖子

aspen RPLUG 的fortran动力学子程序纠错