Aspen-Fortran L-H动力学反应精馏
求帮忙，这是我编写的一份反应精馏的L-H型动力学方程Fortran文件，能够成功导入ASPEN中，不过运行之后不收敛，并且在反应精馏塔的profile中发现reaction中都是0（没有反应），请大家帮忙看看我编写的Fortran有没有什么问题？- 本文出自马后炮化工-让天下没有难学的化工技术，原文地址：https://meng.horse/thread-162263-1-1.html
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 DIB (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_IB,K_DIB, K_TI, K_IP, K_TBA
PARAMETER(K_IB=1)
PARAMETER(K_DIB=2)
PARAMETER(K_TI=3)
PARAMETER(K_IP=4)
PARAMETER(K_TBA=5)
PARAMETER(N_COMP=5)
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
C
C
C DECLARE SYSTEM FUNCTIONS
C
REAL*8 DLOG
C
C DECLARE LOCAL VARIABLES
C
INTEGER IMISS, IDBG
REAL*8 REAL(NREAL),RMISS,DKR(2),
+ RATE(2),RATNET(2),KA,KB
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
DKR(1)=(2.278D+10*DEXP(-30000.D+00/8.314/T))
DKR(2)=(1.81D+9*DEXP(-1800.D+00/8.314/T))
KA=7.D+00
KB=1.2D-01
IF(IDBG.GE.1)THEN
WRITE(MAXWRT_MAXBUF(1),9010) DKR(1),DKR(2),KA,KB
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 use mass holdup as kgcata
CSCRATCH THAT...THATS WITH VERSION 9, IN 10 USE REALB(1) AS TOTAL MASS.
IF (ACTIV(K_IB) .GT. 0.D0) THEN
RATE(1)=(DKR(1)*(ACTIV(K_IB))**2.d0/(ACTIV(K_IB)+KA*ACTIV(K_TBA)
*+KB*ACTIV(K_IP))**2.d0)
RATE(2)=(DKR(2)*ACTIV(K_IB)*ACTIV(K_DIB)/(ACTIV(K_IB)+KA
**ACTIV(K_TBA)+KB*ACTIV(K_IP))**3.d0)
ELSE
RATE(1) = 0.D+00
RATE(2) = 0.D+00
END IF
RATNET(1)=RATE(1)*REALB(1)*1.D-03
RATNET(2)=RATE(2)*REALB(1)*1.D-03
RATES(K_IB)=-RATNET(1)-RATNET(1)-RATNET(2)
RATES(K_DIB)=RATNET(1)-RATNET(2)
RATES(K_TI)=RATNET(2)
RATES(K_TBA)=0.D+00
RATES(K_IP)=0.D+00
IF(IDBG.GE.1)THEN
WRITE(MAXWRT_MAXBUF(1),9030) NSTAGE,RATE(1),RATE(2),RATNET(1),
+ RATNET(2)
CALL DMS_WRTTRM(1)
ENDIF
RETURN
#undef P_MAX3
END
C:\Users\Administrator\Desktop\1561704934(1) 编程我一窍不通，感觉很高深啊！ 學習學習!~ 厉害了，完全看不懂 可以參考MTBE的例題文件 楼主解决了吗 我也是出现了相同的问题
学习学习，谢谢分享 {:1110_549:}{:1110_549:}
页:
[1]