Aspen Dynamic 自带例题反应精馏MTBE Fortran代码请教

楚天湘水

大家好：最近两天看了一下反应精馏ＭＴＢＥ的合成，反应式和动力学方程如下

731271733

楚天湘水 发表于 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)？

楚天湘水

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_MAX3

C  好像是释放一个宏的意思，具体不太懂 ，求解释。

     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第九章也提到了这个例子，但没具体解释代码。我上传第九章节，有兴趣的可以看一看。

附件已经在一楼了，谢谢大家的关注。

PSESKY

楚天湘水 发表于 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的网站下载资料学习，希望你能在论坛多多发言，共同进步。

damoguyan

IDBG是个什么变量呢，你只赋给了它初值，它中间怎么变的，为什么会大于1

damoguyan

x(1,1)中的x是什么样的一个矩阵，也不是很清楚，求指导

damoguyan

如果IDBG值就是0的话，那接下来的几个输出和计算就不运行了

xinquanliuyun

楼主求指教啊？最近在做这个但是一直找不到动力学数据，楼主qq多少？我的505606742

hbeden

我也想问楼主个问题，偏移量 offset 是什么玩意。

生命之怒放

{:1106_362:} xue xi le ,tai ganxie louzhu

liuhui8803

学习了，感谢楼主

18562317750

DECLARE VARIABLES USED IN DIMENSIONING
变量是如何选取的？请大神们指教

18562317750

楚天湘水 发表于 2013-7-8 10:14
不好意思，今天才看到你的回复。
首先谢谢你的热心回复，1,2的回答让我明白了PLEX是什么东西，以后还要 ...

请问编程中变量是如何选取的？？

liangyf2010

请问有没有这方面可参考的书籍资料

试剑阁阁主

XUEIXIYIXIA

蔡育群

请问楼主，图中的real值与integer是如何设定的

carbon002

感觉楼主分享，最近正在学习动力学的应用，感觉难度比较大

楚天湘水

代码

C***********************************************************************

C LICENSED MATERIAL.  PROPERTY OFASPEN TECHNOLOGY, INC.  TO BE      

C TREATED AS ASPEN TECH PROPRIETARYINFORMATION UNDER THE TERMS      

C  OFTHE ASPEN PLUSSUBSCRIPTION 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      VARIABLE  I/O  TYPE    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 BYRATEFRAC **

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 REACTIONRATES

C                                             (KMOL/SEC)

C      RATEL      O    R        NRL         INDIVIDUAL REACTIONRATES

C                                             INTHE LIQUID PHASE

C                                            (KMOL/SEC)

C                                             *USED ONLY BY RATEFRAC **

C      RATEV      O    R        NRV         INDIVIDUAL REACTIONRATES

C                                             INTHE 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"

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_GTPRPS

C：这一句没看懂，查不到DMS_GTPRPS的意思，

C  In version 2006, a projectto reduce memory usage by Aspen Plus has

C  resulted in a major changeto the way these calculated properties are stored.

C  As a result, labeled commonIPOFF3_IPOFF3 can no longer be used in this

#include "dms_plex.cmn"

     EQUIVALENCE(B(1),IB(1))

C

C     DATA STATEMENTS

C

     DATA IDBG/0/

C   IDBG参数设置为0 ，赋值语句，非执行语句

C   IDBG参数由0到1怎么变化的，不清楚？

C   我尝试改变IDBG的命名，比如改成ＩＤＢＨ，可以通过ａｓｐｃｏｍｐ编译，但运

C    行ａｓｐｅｎ的时候会出现错误，在user model里面查找IDBG，没有找到这个参数。

楚天湘水

C    thermodynamic rate constant DKA

C    ===============================

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)

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  这些写入的值最后可以在哪查看呢？没有找到。

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    KER 的user module 解释是：Error return code:

C            ≠0 if an error or warning conditionoccurred      =0 otherwise

C    set offset to get activity coefficients

C    (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,NCOMP

C       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 rate

C    =============

C    RATE : mol/s/kgcata

C    Q    : equiv/kgcata

C    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 kgcata

C SCRATCH THAT...THATS WITH VERSION 9, IN 10 USE REALB(1) AS TOTAL MASS.

     RATNET=RATE*REALB(1)*1.D-03

C  REALB user module的解释是：Vector of realparameters (from unit operation block UserSubroutine form)。下图中1000就应该是REALB（1）。

C