Aspen Dynamic 自带例题反应精馏MTBE Fortran代码请教
本帖最后由 楚天湘水 于 2013-1-9 12:15 编辑大家好：最近两天看了一下反应精馏ＭＴＢＥ的合成，反应式和动力学方程如下
楚天湘水 发表于 2013-1-9 12:19
C thermodynamic rate constant DKAC =============================== 9010FORMAT(1X,3(G13.6,1X))...
楼主，请教两个问题：
1、你所贴的程序目的是计算反应速率，而计算反应速率需要活度（程序中的ACTIVE变量）。但是为什么在计算活度前先调用了PPMON_FUGLY？这是用来计算逸度系数的。而且程序中计算活度时也没有用到逸度系数（程序中的PHI变量）。
2、计算逸度系数时，程序是这样调用的CALL PPMON_FUGLY(T,P,X(1,1)
+ ,Y,NCOMP,IDX,NBOPST,KDIAG,KPHI,PHI,DPHI,KER)。
其中的参数X（1，1）表示液相中甲醇的摩尔分数。这里为什么用甲醇的，而不是用其他组分的，如IC4、NC4、MTBE，对应的应为X(2,1)、X(3,1)、X(4,1)？ 本帖最后由 楚天湘水 于 2013-1-9 14:34 编辑
C 个人感觉这个值就是每块板所装填的催化剂吧，单位可能是g，所以乘以-3次方C 但代码提供的注释又说是use mass holdup askgcata，不知道二者有什么关系 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_MAX3C好像是释放一个宏的意思，具体不太懂 ，求解释。 END 以上代码红色字体是我加上的注释，同时还加了一张图片。其他部分是aspen提供的例子的源码。现在我再说说我的问题，大家一起讨论一下。1、我感觉这个程序好像是每计算一块存在反应的塔板，都得调用一次这个程序，这个程序每执行一次，都会有一些中间的变量写入文件，但这些变量写入了哪个文件，我不清楚，也查不到这些中间变量的值。2、程序中并没有输入塔板数，温度等变量，是不是程序可以自动调用aspen中的这些值？具体怎么调用的呢？3、IDBG是个什么变量呢，因为我改变他的命名为IDBH，可以通过aspcomp的编译，但aspen运行这个程序不错误，是不是像KER一样，是某个函数的一个输出值呢？查了好几版本的user module ，都没查到这个值或者函数？他是怎样实现又0到1的变化的呢？求解释。4、DMS_GTPRPS是自定义的函数呢？还是aspen user module自身就有的函数呢？5、#undef P_MAX3的意思？网上查了一下，好像是释放宏的意思，是不是程序运行结束，要释放调用过的aspen user module的公用函数呢？没有学过Fortran，求解释。 谢谢大家了，关于Fortran动力学的资料太少了，欢迎大家一起讨论。同时上传这个例子的ＢＫＰ文件和相关的．ｆ文件（因为０６和０７的ｆ文件有一点小小的改变，故两个都上传）这几个ｃａｓｅ是可以在ａｓｐｅｎ安装的目录下找到的： \AspenTech\Aspen Plus Dynamics V7.1\Examples 文件名为ＤｙＭＴＢＥ．ＢＫＰ 和ＲａＭＴＢＥ．ｆ 另外美国里海大学威廉鲁宾教授和台湾大学余政靖教授合著的书Reactive Distillation Design and Control第九章也提到了这个例子，但没具体解释代码。我上传第九章节，有兴趣的可以看一看。
附件已经在一楼了，谢谢大家的关注。 楚天湘水 发表于 2013-1-9 12:20
C 个人感觉这个值就是每块板所装填的催化剂吧，单位可能是g，所以乘以-3次方C 但代码提供的注释又说 ...
回應一下：
1.2.ASPEN計算會將中間的計算值全部存在PLEX裡面，那是一個很大的舉陣。所有的資料以及參數都會儲存在裡面。
我們只要能夠從PLEX抓出我們需要參數即可在FORTRAN中計算。
ASPEN--->FORTRAN CALCULATE--->ASPEN （REPEAT*）
要抓那個參數要去查USER MODEL對每個不同參數的代號的命名，可從手冊上得知。
3. IDBG，KER那個是執行CODE，例如KER＝1代表要求程式計算出你所想知道的值
KER＝0代表不要程式計算任何值
（KER有0-3的變化，各有不同意思，在手冊上有，可直接搜尋KER即可找出）
4. DMS_*一般是指ASPEN裡面的運作，你的例子要求ACTIVITY，在新版中有不同的呼叫方式。
請在手冊上搜尋，DMS_ALIPOFF3(24)，新版2006.5後因應儲存PLEX便利性
要抓ACTIVITY已經改變成ALIPOFF而非書上的那種寫法。
5. 那句可刪掉也可以正常運行，詳細功能不清楚。
建議可以從ASPEN的FORTRAN範例（USER跟KEN）下手研究
6. 關於寫出FORTRAN抓的數值，程式會在相圖目錄下產生一個新的文字檔*.dat。可用WORDPAD打開
可參考OPEN OPEN(10,FILE='fortran.DAT')
WRITE(10,想看的變數）
祝順利
PSESKY 发表于 2013-6-17 17:06
回應一下：
1.2.ASPEN計算會將中間的計算值全部存在PLEX裡面，那是一個很大的舉陣。所有的資料以及參數 ...
不好意思，今天才看到你的回复。
首先谢谢你的热心回复，1,2的回答让我明白了PLEX是什么东西，以后还要多关注这个。
文件的读写过段时间在研究一下。
台湾大学在这方面做的比较好，一段时间我也曾多次登录台湾大学PSE的网站下载资料学习，希望你能在论坛多多发言，共同进步。 IDBG是个什么变量呢，你只赋给了它初值，它中间怎么变的，为什么会大于1
x(1,1)中的x是什么样的一个矩阵，也不是很清楚，求指导 如果IDBG值就是0的话，那接下来的几个输出和计算就不运行了 楼主求指教啊？最近在做这个但是一直找不到动力学数据，楼主qq多少？我的505606742 我也想问楼主个问题，偏移量 offset 是什么玩意。 {:1106_362:} xue xi le ,tai ganxie louzhu 学习了，感谢楼主 DECLARE VARIABLES USED IN DIMENSIONING
变量是如何选取的？请大神们指教 楚天湘水 发表于 2013-7-8 10:14
不好意思，今天才看到你的回复。
首先谢谢你的热心回复，1,2的回答让我明白了PLEX是什么东西，以后还要 ...
请问编程中变量是如何选取的？？ 请问有没有这方面可参考的书籍资料 XUEIXIYIXIA 请问楼主，图中的real值与integer是如何设定的 感觉楼主分享，最近正在学习动力学的应用，感觉难度比较大 代码C***********************************************************************C LICENSED MATERIAL.PROPERTY OFASPEN TECHNOLOGY, INC.TO BE C TREATED AS ASPEN TECH PROPRIETARYINFORMATION UNDER THE TERMS COFTHE ASPEN PLUSSUBSCRIPTION AGREEMENT. *C***********************************************************************C-----------------------------------------------------------------------C COPYRIGHT (C) 1994C ASPEN TECHNOLOGY, INC.C CAMBRIDGE, MAC-----------------------------------------------------------------------CC calculation of MTBE synthesis reaction rateC User Kinetics Subroutine for RADFRAC, BATCHFRAC, RATEFRACC 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)CC VARIABLES IN ARGUMENT LISTCC VARIABLEI/OTYPE DIMENSION DESCRIPTION AND RANGEC NSTAGE I I - STAGE NUMBERC NCOMP I I - NUMBER OF COMPONENTSC NR I I - TOTAL NUMBER OF KINETICC REACTIONSC NRL I I - NUMBER OF LIQUID PHASEC KINETIC REACTIONSC NRV I I - NUMBER OF VAPOR PHASEC KINETIC REACTIONSC 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 BYRATEFRAC **C P I R - STAGE PRESSURE (N/SQ.M)C VF I R - VAPOR FRACTIONC F I R - TOTAL FLOW ON STAGEC (VAPOR+LIQUID)(KMOL/SEC)C X I R NCOMP,3 LIQUID MOLE FRACTIONC Y I R NCOMP VAPOR MOLE FRACTIONC IDX I I NCOMP COMPONENT INDEX VECTORC NBOPST I I 6 OPTION SET BEAD POINTERC KDIAG I I - LOCAL DIAGNOSTIC LEVEL
C STOIC I R NCOMP,NR REACTION STOICHIOMETRYC IHLBAS I I - BASIS FOR LIQUIDC HOLDUP SPECIFICATIONC 1:VOLUME,2:MASS,3:MOLEC HLDLIQ I R - LIQUID HOLDUPC IHLBAS UNITSC 1 CU.M.C 2 KGC 3 KMOLC TIMLIQ I R - LIQUID RESIDENCE TIMEC (SEC)C IHVBAS I I - BASIS FOR VAPORC HOLDUP SPECIFICATIONC 1:VOLUME,2:MASS,3:MOLEC HLDVAP I R - VAPOR HOLDUPC IHVBAS UNITSC 1 CU.M.C 2 KGC 3 KMOLC TIMVAP I R - VAPOR RESIDENCE TIME(SEC)C NINT I I - LENGTH OF INTEGER VECTORC INT I/O I NINT INTEGER VECTORC NREAL I I - LENGTH OF REAL VECTORC REAL I/O R NREAL REAL VECTORC RATES O R NCOMP COMPONENT REACTIONRATESC (KMOL/SEC)C RATEL O R NRL INDIVIDUAL REACTIONRATESC INTHE LIQUID PHASEC (KMOL/SEC)C *USED ONLY BY RATEFRAC **C RATEV O R NRV INDIVIDUAL REACTIONRATESC INTHE VAPOR PHASEC (KMOL/SEC)C * USED ONLY BY RATEFRAC **C NINTB I I - LENGTH OF INTEGER VECTORC (FROM UOS BLOCK)C INTB I/O I NINTB INTEGER VECTORC (FROM UOS BLOCK)C NREALB I I - LENGTH OF REAL VECTORC (FROM UOS BLOCK)C REALB I/O R NREALB REAL VECTORC (FROM UOS BLOCK)C NIWORK I I - LENGTH OF INTEGER WORK 本帖最后由 楚天湘水 于 2013-1-9 14:32 编辑
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
CC DECLARE SYSTEM FUNCTIONSC REAL*8 DLOGCC DECLARE LOCAL VARIABLESC 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)CC#include "dms_maxwrt.cmn"C user model 的解释是write to the Terminal File #include "dms_ipoff3.cmn"C user model 的解释是directly access some calculated or intermediate properties from theplex REAL*8 DMS_GTPRPSC：这一句没看懂，查不到DMS_GTPRPS的意思，CIn version 2006, a projectto reduce memory usage by Aspen Plus has Cresulted in a major changeto the way these calculated properties are stored. CAs a result, labeled commonIPOFF3_IPOFF3 can no longer be used in this #include "dms_plex.cmn" EQUIVALENCE(B(1),IB(1)) CC DATA STATEMENTSC DATA IDBG/0/C IDBG参数设置为0 ，赋值语句，非执行语句C IDBG参数由0到1怎么变化的，不清楚？C 我尝试改变IDBG的命名，比如改成ＩＤＢＨ，可以通过ａｓｐｃｏｍｐ编译，但运C 行ａｓｐｅｎ的时候会出现错误，在user model里面查找IDBG，没有找到这个参数。 本帖最后由 楚天湘水 于 2013-1-9 14:29 编辑
C thermodynamic rate constant DKAC =============================== 9010FORMAT(1X,3(G13.6,1X)) 9000FORMAT(' fugly failed at T=',G12.5,' P=',G12.5,' ker=',I4) 9020FORMAT(' compo ',I3,' mole-frac=',G12.5,' activity=',G12.5) 9030FORMAT(' stage=',I4,' spec-rate=',G12.5,' net-rate=',G12.5)CC BEGIN EXECUTABLE CODEC 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 constantC ====================== 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) ENDIFC这些写入的值最后可以在哪查看呢？没有找到。C calculation of components activitiesC ====================================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) ENDIFC KER 的user module 解释是：Error return code: C ≠0 if an error or warning conditionoccurred =0 otherwiseC set offset to get activity coefficientsC (see vol5, p 11-11 and asp$sor search for 'GAMMAL')C L_GAMMA=IPOFF3_IPOFF3(24) C calculate activities for plex data DO J=1,NCOMPC ACTIV(J)=DEXP(B(L_GAMMA+J))*X(J,1)C aspen2006提供的例子计算活度是用的上面的公式，V7.1的例子改用了下面的公式，DMS_GTPRPS没看明白，求解释，24是活度系数的代码 ACTIV(J)=DEXP(DMS_GTPRPS(24,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 rateC =============C RATE : mol/s/kgcataC Q : equiv/kgcataC DKR: mol/s/equiv Q=4.9D+00 RATE=Q*DKR* * (ACTIV(K_IC4)/ACTIV(K_MEOH) * -ACTIV(K_MTBE)/DKA/ACTIV(K_MEOH)**2)C use mass holdup as kgcataC SCRATCH THAT...THATS WITH VERSION 9, IN 10 USE REALB(1) AS TOTAL MASS. RATNET=RATE*REALB(1)*1.D-03CREALB user module的解释是：Vector of realparameters (from unit operation block UserSubroutine form)。下图中1000就应该是REALB（1）。C
页:
[1]