aspen外挂fortran动力学方程探讨

六月痕

C    User Kinetics Subroutine for RCSTR, RPLUG, RBATCH, PRES-RELIEF

C

      SUBROUTINEUSRKIN (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 以上定义子程序变量

      IMPLICITNONE   !fortran 语句 清除隐含变量

C

C    DECLARE VARIABLES USED IN DIMENSIONING

C

      INTEGERNSUBS, NINT,  NPO,   NIWORK,NWORK,

     +        NC,    NR,   NTCAT, NTSSAT,NCOMP,

     +        NRALL, NUSERV,NINTR,NREALR,NIWR,

     +        NWR

C

#include "ppexec_user.cmn"                ！公用区，存储运行时间控制标志

      EQUIVALENCE(RMISS, USER_RUMISS)    ！EQUIVALENCE 是fortran中的等价语句

      EQUIVALENCE(IMISS, USER_IUMISS     ！USER_RUMISS 实型补遗代码  

                                            USER_IUMISS 实型补遗代码

C

C

C

C

C.....RCSTR...

#include "rcst_rcstri.cmn"        ！公共模块包含 RCSTR 模块的整型结构参数

#include "rxn_rcstrr.cmn"         ！公共模块包含 RCSTR 模块的实型结构参数

C

C.....RPLUG...

#include "rplg_rplugi.cmn"        !公共模块包括RPLUG模块的整型结构参数

#include "rplg_rplugr.cmn"        !公共模块包括RPLUG模块的实型结构参数

      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"     ！公共区 实型性质（例如温度。压强等）P71

      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"                      ！SHS_STWORK 是物流闪蒸工作偏移量公用区 P134

      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)

#include "pputl_ppglob.cmn"    ！由所有模拟程序性质程序使用的物性全局公用区

#include "dms_ncomp.cmn"       ！包括各种类别组分的数目和与组分相关的流段的长度

#include "dms_plex.cmn"        ！DMS_PLEX含有Plex的长的内核存储区域，存储fortran调用所需的物性参数

      EQUIVALENCE(IB(1), B(1))  ！B 实型的Plex区域  IB 整型的Plex区域

C

C    DECLARE ARGUMENTS

C

      INTEGERIDXSUB(NSUBS),ITYPE(NSUBS), INT(NINT),

     +       IDS(2),NBOPST(6,NPO),IWORK(NIWORK),

     +        IDX(NCOMP),   INTR(NINTR), IWR(NIWR),

     +        NREAL, KCALL, KFAIL,KFLASH,I

      REAL*8SOUT(1),      WORK(NWORK),

     +       STOIC(NC,NSUBS,NR),  RATES(1),

     +       FLUXM(1),     FLUXS(1),     RATCAT(NTCAT),

     +       RATSSA(NTSSAT),      Y(NCOMP),

     +       X(NCOMP),     X1(NCOMP),    X2(NCOMP)

      REAL*8RATALL(NRALL),USERV(NUSERV),

     +       REALR(NREALR),WR(NWR),      XCURR

C

C

C    DECLARE LOCAL VARIABLES

C

      INTEGERIMISS, MKBAS, MKPHAS,MTAPP, MKBASS,

     +        MTAPPS,LMW

      REAL*8REAL(NREAL), XL(3), B(1),  RMISS, XLEN,

     +       DIAM,  TEMP, PRES,  VFRAC, BETA,

     +       VVAP,  VLIQ, VLIQS, SSALT, VSALT,

     +       FSALT, RTEMP, DENLIQ,CONAA,CONACT,

     +       FTERM, STERM, RATE,  SUMR

      INTEGERXMW

C

C    INITIALIZE RATES

C

C

C STATEMENT FUNCTIONS FOLLOW

C

C

C    BEGIN EXECUTABLE CODE

C

      XMW(I) =LMW + I         ！组分分子量，The MW values start from B(LMW + 1)  P42

C

C SET PLEX OFFSETS

C

      LMW = IPOFF1_IPOFF1(306)   ！??? 没找到依据，猜测确定组分数据区域的DMS_PLEX偏移量

      DO100 I = 1, NC

        RATES(I) = 0D0

        XL(I)    = 0D0

100 CONTINUE

C                               ！以上排数偏移量的影响

      DO101 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)         : Totalflowrates(kg-moles/sec)

C    SOUT(NCC+2)         :Temperature(K)

C    SOUT(NCC+3)         :Pressure(N/SQM)

C    SOUT(NCC+4)         : Massenthalpy(J/KG)

C    SOUT(NCC+5)         : Molar vaporfraction

C    SOUT(NCC+6)         : Molar liquidfraction

C    SOUT(NCC+7)         : Massentropy(J/KG-K)

C    SOUT(NCC+8)         : Massdensity(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

      DENLIQ = 1.0 / STWORK_VL         ！求密度

C

C    Compute Concentration for reactants:

C

C    CONC(i) = X(i) * RHO

C

      CONAA = XL(1) * DENLIQ          ！计算Allyl-Alcohol浓度

      CONACT = XL(2) * DENLIQ          ！计算Acetone浓度

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                                                           ！动力学方程主代码

      FTERM = REALR(1) * EXP( -REALR(2) / PPGLOB_RGAS / RTEMP )   ！计算k值，其中real(1), real(2)均为subroutine中的可变整型量

      STERM = conaa * conact**0.5

      RATE  = FTERM * STERM * RCSTRR_VOLRC       ！ ？？不知道为什么要乘以反应器体积，难道求速率方程式时是以单位体积为基准？

      RATES(1) = -rate                           ！每个组分必须均指定速率（非常重要否则计算出错）

      RATES(2) = -rate

      RATES(3) = rate

C

C    Checking Mass Balance for RATES Vector       ！核对质量是否守恒

C

      SUMR = 1.0

      DO102 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,*)

     +          'CalculatedRates Not in Mass Balance!'

         WRITE(USER_NHSTRY,*)

     +          'Errorsin User Kinetics Subroutines!'

      ENDIF

       write(USER_NHSTRY,*)'volrc = ', RCSTRR_VOLRC

      write(USER_NHSTRY,*)'vfrrc = ', RCSTRR_VFRRC

       RETURN

#undef P_NPOFF1

      END

问题：

1  LMW = IPOFF1_IPOFF1(306)   ！??? 没找到依据，猜测确定组分数据区域的DMS_PLEX偏移量

2  #include "dms_ipoff1.cmn"        ！？？ aspen内部之间相互调用的存储空间

3 RATE  = FTERM *STERM * RCSTRR_VOLRC       ！ ？？不知道为什么要乘以反应器体积，难道求速率方程式时是以单位体积为基准？

体会：

  首先感谢论坛提供的资料，以上实例来自论坛，大家可以搜索下载。上面注释中的页码对应ASPEN_PLUS_10.0_用户模型.pdf。对于aspen调用fortran动力学方

程，正处于学习阶段。由于论坛此方面的资料较少，特开此贴希望大家积极参与讨论，如果上述有不妥之处，也希望能者不吝指教。

  对于程序的书写，谈一下个人感受。Aspen自带有模板D:（安装盘）\Program Files (x86)\AspenTech\Aspen Plus V7.3\Engine\User下的usrkin.f。可以在此基础上添加主代码即可。

Top: 1 fortran 程序名字不要超过6字符，程序内均采用双精度型变量

     2 合理调用aspen内部物性参数（难点），注意偏移量的影响

     3 对方程式每个组分指定反应速率（曾经漏掉过，一直报错）

     4 最好核对质量守恒，曾经见过一个程序没有此段代码。否则容易出错（曾经写过一程序，一直提示质量不守恒，不知是不是缺少此段代码的过）

以下隐藏为bkp文件及相应ppt（论坛可找到)

                               
登录/注册后可看大图

该贴已经同步到 六月痕的微博

下载地址

• 

  USER KINETIC.pptx

  2013-12-7 20:54 大小：274.7 KB 下载次数：1112
• 

  reaction.rar

  2013-12-7 20:54 大小：312.38 KB 下载次数：895

风云天下

不会很伤啊

resgroupinc

请问，能外挂c++的库吗

zjwmcl

都有什么啊？

六月痕

zjwmcl 发表于 2013-12-7 21:06
都有什么啊？

附件是讲解用户动力学程序的ppt，压缩包是bkp文件和fortran文件

华歆

这确实是很需要探讨的问题！尤其对反应器来说很有必要！

Wordsworth

看不懂呀

微博评论

六月痕,亲，可不可以关注下我哦[抱抱]\n\n

                               
登录/注册后可看大图

来自 Fennie梓珊 的新浪微博

高天云

好资料，不学不能说会ASPEN

yeyu

看不懂，代码控

liuxin

限售条件股份

supercaosong

下下来学习学习

lingerking

[udu_catsoul=5]给力！[/udu_catsoul]高手啊

SuperDJ

没整过！围观！

fg2009

好资料，正好需要哦

liuyanfei

是关于Fortran的ppt吗，先看看

bkqcycyqm

一提到编程，头就大了。

dragon747005

正需要fortran资料，学习学习，楼主辛苦了。

tianshanxue

{:1106_362:}
楼主厉害

840882292

正在做动力学的问题，希望有所帮助

qjhnfy

正需要fortran资料，学习学习，楼主辛苦了。

roger309-l0

謝謝樓主分享