本帖最后由 六月痕 于 2013-12-7 20:58 编辑
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（论坛可找到)
该贴已经同步到 六月痕的微博
如果在fortran中没有定义组分的顺序，那么组分的顺序就会和aspen中component中的顺序是一样的，这个时候写rates（1），rates（2）......的时候就要特别注意，记得和组分相对应，而不是只看那个反应式，否则会出现质量不守恒的情况。
好像一些简单的编程可以直接在CALCULATOR里面进行fortran语言编程，但是偶不会，连在calculator里设置变量都出问题，现在又急需用到这块，太需要学习学习了
guanxian1990 于 2014-05-10 20:49:20 补充以下内容:
看着都觉得头大啊
楼主，你好，我想问一下，通过fortran子程序可以设置同一反应器内的多个化学反应的动力学参数吗？看了一些例子，都是针对某一个反应的，不知道有不有涉及多个反应的例子？
楼主，我想问下，这个反应动力学代码适不适用于反应精馏塔呢？我要是写反应精馏塔的代码需要用到密度跟相对分子质量，还有体积之类的参数，能不能借用您这里面的相对分子质量，和密度的代码呢？还有计算密度的公式时：STWORK_VL代表什么意思呢？
SPECIFIC VOLUME MODEL VL2RKT WAS CALLED WITH INVALID INPUT：ALL MOLE FRACTIONS ARE ZERO. PROPERTIES WILL NOT BE CALCULATED. 请问这个错误应该从那儿入手改正呢？
USER KINETIC.pptx (274.7 KB) 
reaction.rar (312.38 KB) 
AspenPlusUserModelsV7_3-Ref.pdf (1.51 MB) 
