PROII用户自定义热力学方法求助? 我试图通过用户自定义热力学方法添加FORTRAN子程序,链接到PROII中,但是结果运行不对,麻烦各位前辈帮忙看一下我下面的子程序问题出在了哪儿,不胜感激! 程序介绍:通过实验数据拟合NRTL模型中的分子作用参数 aij和bij ,将拟合出的参数值通过子程序写入到PROII中,其中 α设为0.2, τij=aij+bij/T,此程序为计算物质的K值程序,参数值通 过PROII图形界面中的热力学方法设定时写入,再通过CUDATA模块读入。 PS:我知道可以直接将参数通过PROII图形界面的热力学方法设定时写入作用参数值,不用编程,我只是想学会编程的过程,多谢! SUBROUTINE UKHS1 PARAMETER (MAXCN = 300) INCLUDE 'FACTOR.CMN' INCLUDE 'UTHERX.CMN' INCLUDE 'CUDATA.CMN' DIMENSION AF(MAXCN,MAXCN),TO1(MAXCN,MAXCN),TO2(MAXCN,MAXCN), * G(MAXCN,MAXCN),G1(MAXCN,MAXCN),GX(MAXCN),GXT(MAXCN), * GXA(MAXCN),GXB(MAXCN),A(MAXCN),B(MAXCN),C(MAXCN), * D(MAXCN),VPS(MAXCN),GAMMA(MAXCN),DGAMMA(MAXCN), * STREAM(MAXCN+10) C KDCALL=1 ZXX1=0. ZXX2=0. C DO 15 I=1,MAXCN+10 C STREAM(I)=0.0 C 15 CONTINUE ******************************************************************** * DECIDING THE CALCULATION TYPE * ******************************************************************** JXXFLG=0 IF(IXXFLG .NE. -1) GO TO 10 JXXFLG=3 GO TO 999 10 CONTINUE ************************************************************** * INITIALIZE & INPUT DATA FROM STORAGE BLOCK /CUDATA/ * * WATCH OUT FOR THE COMPONENT ORDER * ************************************************************** TDERIV=0.0 C DO 30 I=1,MAXCN+10 C STREAM(I)=0.0 C 30 CONTINUE C DO 31 I=1,MAXCN C GAMMA(I)=0.0 C DGAMMA(I)=0.0 C GX(I)=0.0 C GXT(I)=0.0 C GXA(I)=0.0 C GXB(I)=0.0 C A(I)=0.0 C B(I)=0.0 C C(I)=0.0 C D(I)=0.0 C VPS(I)=0.0 C DO 32 J=1,MAXCN C AF(I,J)=0.0 C TO1(I,J)=0.0 C TO2(I,J)=0.0 C G(I,J)=0.0 C G1(I,J)=0.0 C 32 CONTINUE C 31 CONTINUE K=0 M=0 DO 11 J=1,NOCZ CALL COINUM(J, IS, JC) DO 12 I=1,NOCZ CALL COINUM(I, IS, IC) K=K+1 AF(IC,JC)=TDATA(K) TO1(IC,JC)=TDATA(10+K) TO2(IC,JC)=TDATA(20+K) C 1001 FORMAT('TO1(IC,JC)= ','IC= ','JC= ','NOCZ= ',5F10.4, C * 5F10.4,5F10.4,5F10.4) C WRITE(IOUT,1001) TO1(IC,JC), IC, JC, NOCZ 12 CONTINUE M=M+10 A(JC)=TDATA(21+M) B(JC)=TDATA(22+M) C(JC)=TDATA(23+M) D(JC)=TDATA(24+M) 11 CONTINUE IF(IXXFLG .EQ. 0) THEN GO TO 999 ELSE IF(IXXFLG .NE. 1) THEN JXXFLG=99 GO TO 999 ENDIF **************************************************************** * CALCULATING VAPOR PRESSURE USING ANTONIE EQUATION * * OR FUNCTION TTPROP * **************************************************************** STREAM(1)=1. STREAM(2)=TEMPX DO 61 I=1,NOCZ STREAM(I+10)=1. 61 CONTINUE STREAM(3)=PRESX DO 62 I=4,10 STREAM(I)=0.0 62 CONTINUE C ITYPE=1 C IFAZE=1 C CALL TTPROP(ITYPE,IFAZE,VPS,STREAM,IERR) TK=(TEMPX+TCONVT)*TFAC/YKTOR PKA=(PRESX+PCONVT)*PFAC/ATM/101.325 DO 60 I=1,NOCZ VPS(I)=EXP(A(I)+B(I)/TK+C(I)*log(TK)+D(I)*TK**2) 60 CONTINUE ***************************************************************************** * CALCULATING GAMMA-VALUE USING NRTL EQUATION * ***************************************************************************** DO 20 I=1,NOCZ DO 20 J=1,NOCZ G(I,J)=EXP(-AF(I,J)*(TO1(I,J)+TO2(I,J)/TK)) G1(I,J)=0.2*(TO1(I,J)+TO2(I,J)/TK) 20 CONTINUE DO 50 I=1,NOCZ GX(I)=0.0 GXT(I)=0.0 GXA(I)=0.0 GXB(I)=0.0 DO 40 M=1,NOCZ GX(I)=GX(I)+XLL(M)*G(M,I) GXT(I)=GXT(I)+XLL(M)*G(M,I)*(TO1(M,I)+TO2(M,I)/TK) GXA(I)=GXA(I)+XLL(M)*G(M,I)*AF(M,I)*TO2(M,I) GXB(I)=GXB(I)+XLL(M)*TO2(M,I)*G(M,I)*(1.0-G1(M,I)) 40 CONTINUE 50 CONTINUE DO 100 I=1,NOCZ GAMMA(I)=GXT(I)/GX(I) DGAMMA(I)=-GXB(I)/GX(I)-GXT(I)*GXA(I)/(GX(I))**2 DO 80 M=1,NOCZ GAMMA(I)=GAMMA(I)+XLL(M)*G(I,M)* * ((TO1(I,M)+TO2(I,M)/TK)-GXT(M)/GX(M))/GX(M) DGAMMA(I)=DGAMMA(I)+XLL(M)*G(I,M)*TO2(I,M)*(G1(I,M)-1.0)/GX(M) * -XLL(M)*G(I,M)*(AF(I,M)*TO2(I,M)*GXT(M)+ * GXA(M)*(TO1(I,M)+TO2(I,M)/TK)-GXB(M))/(GX(M))**2 * +2*XLL(M)*G(I,M)*GXA(M)*GXT(M)/(GX(M))**3 80 CONTINUE GAMMA(I)=EXP(GAMMA(I)) DGAMMA(I)=GAMMA(I)*DGAMMA(I)/TK**2 100 CONTINUE ***************************************************************** * CALCULATING K-VALUE * ***************************************************************** C PTOL=PKA PTOL=PRESX IF (KDECNT .GT. 0) THEN C CALL UH2OP(4,1,TEMPX,PRESX,TDELTX,VWPS,TDERIV) C VPS(KDECNT)=(VWPS+PCONVT)*PFAC/ATM/101.325 C PTOL=PKA*(1.0-XVV(KDECNT)) CALL UH2OP(4,1,TEMPX,PRESX,TDELTX,VPS(KDECNT),TDERIV) PTOL=PRESX*(1.0-XVV(KDECNT)) ENDIF DO 70 I=1,NOCZ XKV(I)= GAMMA(I)*VPS(I)/PTOL 70 CONTINUE IF (KDECNT .GT. 0) THEN YXSAT=VPS(KDECNT)/PTOL !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! TF=(TEMPX+TCONVT)*TFAC-FTOR TF=AMAX1(TF,0.01) XXSOL=10**(-9.0306+3.072*LOG10(TF)) XKV(KDECNT)=YXSAT/XXSOL ENDIF ******************************************************** * CALCULATING dK/dT * ******************************************************** IF (TDELTX .NE.0) THEN DO 90 I=1,NOCZ DRVX(I)=DGAMMA(I)*VPS(I)/PTOL 90 CONTINUE IF (KDECNT .GT. 0) THEN DRVX(KDECNT)=XKV(KDECNT)*TDERIV ENDIF ENDIF 999 CONTINUE RETURN END查看更多1个回答 . 1人已关注