ASPEN Fortran 代码求解
本帖最后由 fzm1108 于 2015-9-17 17:27 编辑在Aspen的实例中MTBE反应精馏中有如下代码：C***********************************************************************CLICENSED MATERIAL.PROPERTY OF ASPEN TECHNOLOGY, INC.TO BE *
CTREATED AS ASPEN TECH PROPRIETARY INFORMATION UNDER THE TERMS *
COF THE ASPEN PLUS SUBSCRIPTION AGREEMENT. *
C***********************************************************************
C-----------------------------------------------------------------------
C COPYRIGHT (C) 1994
C ASPEN TECHNOLOGY, INC.
C CAMBRIDGE, MA
C-----------------------------------------------------------------------
C
C calculation of MTBE synthesis reaction rate
C User Kinetics Subroutine for RADFRAC, BATCHFRAC, RATEFRAC
C
SUBROUTINE RAMTBE (NSTAGE, NCOMP, NR, NRL, NRV,
2 T, TLIQ, TVAP, P, VF,
3 F, X, Y, IDX, NBOPST,
4 KDIAG,STOIC, IHLBAS, HLDLIQ,TIMLIQ,
5 IHVBAS, HLDVAP,TIMVAP, NINT, INT,
6 NREAL,REAL, RATES,RATEL, RATEV,
7 NINTB,INTB, NREALB, REALB, NIWORK,
8 IWORK,NWORK, WORK)
C
C VARIABLES IN ARGUMENT LIST
C
C VARIABLEI/OTYPE DIMENSION DESCRIPTION AND RANGE
C NSTAGE I I - STAGE NUMBER
C NCOMP I I - NUMBER OF COMPONENTS
C NR I I - TOTAL NUMBER OF KINETIC
C REACTIONS
C NRL I I - NUMBER OF LIQUID PHASE
C KINETIC REACTIONS
C NRV I I - NUMBER OF VAPOR PHASE
C KINETIC REACTIONS
C T I R - STAGE TEMPERATURE (K)
C TLIQ I R - LIQUID TEMPERATURE (K)
C * USED ONLY BY RATEFRAC **
C TVAP I R - VAPOR TEMPERATURE (K)
C * USED ONLY BY RATEFRAC **
C P I R - STAGE PRESSURE (N/SQ.M)
C VF I R - VAPOR FRACTION
C F I R - TOTAL FLOW ON STAGE
C (VAPOR+LIQUID) (KMOL/SEC)
C X I R NCOMP,3 LIQUID MOLE FRACTION
C Y I R NCOMP VAPOR MOLE FRACTION
C IDX I I NCOMP COMPONENT INDEX VECTOR
C NBOPST I I 6 OPTION SET BEAD POINTER
C KDIAG I I - LOCAL DIAGNOSTIC LEVEL
C STOIC I R NCOMP,NR REACTION STOICHIOMETRY
C IHLBAS I I - BASIS FOR LIQUID
C HOLDUP SPECIFICATION
C 1:VOLUME,2:MASS,3:MOLE
C HLDLIQ I R - LIQUID HOLDUP
C IHLBAS UNITS
C 1 CU.M.
C 2 KG
C 3 KMOL
C TIMLIQ I R - LIQUID RESIDENCE TIME
C (SEC)
C IHVBAS I I - BASIS FOR VAPOR
C HOLDUP SPECIFICATION
C 1:VOLUME,2:MASS,3:MOLE
C HLDVAP I R - VAPOR HOLDUP
C IHVBAS UNITS
C 1 CU.M.
C 2 KG
C 3 KMOL
C TIMVAP I R - VAPOR RESIDENCE TIME (SEC)
C NINT I I - LENGTH OF INTEGER VECTOR
C INT I/O I NINT INTEGER VECTOR
C NREAL I I - LENGTH OF REAL VECTOR
C REAL I/O R NREAL REAL VECTOR
C RATES O R NCOMP COMPONENT REACTION RATES
C (KMOL/SEC)
C RATEL O R NRL INDIVIDUAL REACTION RATES
C IN THE LIQUID PHASE
C (KMOL/SEC)
C * USED ONLY BY RATEFRAC **
C RATEV O R NRV INDIVIDUAL REACTION RATES
C IN THE VAPOR PHASE
C (KMOL/SEC)
C * USED ONLY BY RATEFRAC **
C NINTB I I - LENGTH OF INTEGER VECTOR
C (FROM UOS BLOCK)
C INTB I/O I NINTB INTEGER VECTOR
C (FROM UOS BLOCK)
C NREALB I I - LENGTH OF REAL VECTOR
C (FROM UOS BLOCK)
C REALB I/O R NREALB REAL VECTOR
C (FROM UOS BLOCK)
C NIWORK I I - LENGTH OF INTEGER WORK
C VECTOR
C IWORK I/O I NIWORK INTEGER WORK VECTOR
C NWORK I I - LENGTH OF REAL WORK VECTOR
C WORK I/O R NWORK REAL WORK VECTOR
C
C***********************************************************************
C
IMPLICIT NONE
C
C DECLARE VARIABLES USED IN DIMENSIONING
C
INTEGER NCOMP, NR, NRL, NRV, NINT,
+ NINTB, NREALB,NIWORK,NWORK, N_COMP
C
C DECLARE PARAMETERS & VARIABLES USED IN PARAMETERS
C
INTEGER K_MEOH,K_IC4, K_NC4, K_MTBE
PARAMETER(K_MEOH=1)
PARAMETER(K_IC4=2)
PARAMETER(K_NC4=3)
PARAMETER(K_MTBE=4)
PARAMETER(N_COMP=4)
C component order
C ===============
C this routine assumes that the components are in this order :
C
C DECLARE ARGUMENTS
C
INTEGER IDX(NCOMP), NBOPST(6), INT(NINT),
+ INTB(NINTB),IWORK(NIWORK),NSTAGE,
+ KDIAG, IHLBAS,IHVBAS,NREAL, KPHI,
+ KER, L_GAMMA, J
REAL*8 X(NCOMP,3), Y(NCOMP),
+ STOIC(NCOMP,NR), RATES(NCOMP),
+ RATEL(NRL), RATEV(NRV),
+ REALB(NREALB),WORK(NWORK),B(1),T,
+ TLIQ,TVAP,P, VF, F
REAL*8 HLDLIQ,TIMLIQ,HLDVAP,TIMVAP,TZERO,
+ FT
C
C
C DECLARE SYSTEM FUNCTIONS
C
REAL*8 DLOG
C
C DECLARE LOCAL VARIABLES
C
INTEGER IMISS, IDBG
REAL*8 REAL(NREAL),RMISS, C1, C2, C3,
+ C4, C5, C6, DKA, DKR,
+ Q, RATE,RATNET
REAL*8 PHI(N_COMP)
REAL*8 DPHI(N_COMP)
REAL*8 ACTIV(N_COMP)
C
#include "ppexec_user.cmn"
EQUIVALENCE (RMISS, USER_RUMISS)
EQUIVALENCE (IMISS, USER_IUMISS)
C
C
#include "dms_maxwrt.cmn"
#include "dms_ipoff3.cmn"
#include "dms_lclist.cmn"
INTEGER FN
#include "dms_plex.cmn"
EQUIVALENCE(B(1),IB(1))
FN(J)=J+LCLIST_LBLCLIST
C
C DATA STATEMENTS
C
DATA IDBG/0/
C thermodynamic rate constant DKA
C ===============================
9010 FORMAT(1X,3(G13.6,1X))
9000 FORMAT(' fugly failed at T=',G12.5,' P=',G12.5,' ker=',I4)
9020 FORMAT(' compo ',I3,' mole-frac=',G12.5,' activity=',G12.5)
9030 FORMAT(' stage=',I4,' spec-rate=',G12.5,' net-rate=',G12.5)
C
C BEGIN EXECUTABLE CODE
C
TZERO=298.15D+00
C1=-1493.D+00
C2=-77.4D+00
C3=0.508D+00
C4=-0.913D-03
C5=1.11D-06
C6=-0.628D-09
FT=C1*(1.D+00/T-1.D+00/TZERO)
*+C2*DLOG(T/TZERO)
*+C3*(T-TZERO)
*+C4*(T*T-TZERO*TZERO)
*+C5*(T**3-TZERO**3)
*+C6*(T**4-TZERO**4)
DKA=284D+00*DEXP(FT)
C reaction rate constant
C ======================
DKR=3.67D+12*DEXP(-11110.D+00/T)
IF(IDBG.GE.1)THEN
WRITE(MAXWRT_MAXBUF(1),9010) FT,DKA,DKR
CALL DMS_WRTTRM(1)
ENDIF
C calculation of components activities
C ====================================
C calculate only fugacity coefficient
KPHI=1
C fugacity coefficient of components in the mixture
CALL PPMON_FUGLY(T,P,X(1,1)
+ ,Y,NCOMP,IDX,NBOPST,KDIAG,KPHI,PHI,DPHI,KER)
IF(KER.NE.0)THEN
WRITE(MAXWRT_MAXBUF(1),9000) T,P,KER
CALL DMS_WRTTRM(1)
ENDIF
C set offset to get activity coefficients
C (see vol5, p 11-11 and asp$sor search for 'GAMMAL')
L_GAMMA=IPOFF3_IPOFF3(24)
C calculate activities for plex data
DO J=1,NCOMP
ACTIV(J)=DEXP(B(FN(L_GAMMA)+J))*X(J,1)
END DO
IF(IDBG.GE.1)THEN
DO J=1,NCOMP
WRITE(MAXWRT_MAXBUF(1),9020) J,X(J,1),ACTIV(J)
CALL DMS_WRTTRM(1)
END DO
ENDIF
C reaction rate
C =============
C RATE : mol/s/kgcata
C Q : equiv/kgcata
C DKR: mol/s/equiv
Q=4.9D+00
IF (ACTIV(K_IC4) .GT. 0.D0 .AND. ACTIV(K_MEOH) .GT. 0.D0) THEN
RATE=Q*DKR*
* (ACTIV(K_IC4)/ACTIV(K_MEOH)
* -ACTIV(K_MTBE)/DKA/ACTIV(K_MEOH)**2)
ELSE
RATE = 0.D0
END IF
C use mass holdup as kgcata
CSCRATCH THAT...THATS WITH VERSION 9, IN 10 USE REALB(1) AS TOTAL MASS.
RATNET=RATE*REALB(1)*1.D-03
RATES(K_IC4)=-RATNET
RATES(K_MEOH)=-RATNET
RATES(K_MTBE)=RATNET
RATES(K_NC4)=0.D+00
IF(IDBG.GE.1)THEN
WRITE(MAXWRT_MAXBUF(1),9030) NSTAGE,RATE,RATNET
CALL DMS_WRTTRM(1)
ENDIF
RETURN
#undef P_MAX3
END
在上述的代码中求得反应速率之后用到 REALB(1)，按上述的解释似乎为实型指针，请问此处的意义是什么，（1）的意义又是什么呢，如何选择holdup中的不同条件？
{:1106_362:}{:1106_362:}{:1106_362:} 此外关于Fortran中怎么引用液相摩尔分率呢？
页:
[1]