C********************************************************************** C* PRMIX.FOR * C* PRMIX.FOR CALCULATES PROPERTIES OF MIXTURE USING PENG-ROBINSON * C* EQUATION OF STATE. * C* 1994.3 R.AKASAKA * C* * C* Modified for ver.11.1 ( Feb, 1999 ) * C* - Add SUBMXT and SUBMXH * C* - Change return value of IPHASE * C* single phase = 1 * C* two phase = 2 * C* - Add O2-CO2 mixture ( KOMBI = 17 ) * C* - Subroutine FMF001 is replaced to original code. * C* * C* Modified for ver.12.1 ( May 7, 2001 ) * C* - return value of IDENTM: 11.1 -> 12.1 * C********************************************************************** C----------------------------------- C FMF001 SOLVES CUBIC EQUATION. C----------------------------------- SUBROUTINE FMF001(A, X, IND) IMPLICIT DOUBLE PRECISION (A-H, O-Z) DIMENSION A(3), X(3), AA(2), XX(2) DATA PI /3.1415926535897932385/ DATA DELTA /5D-2/ C A1 = A(1) A2 = A(2) A3 = A(3) P = A2/3D0 - A1**2/9D0 Q = A3/2D0 - A1*A2/6D0 + A1**3/2.7D1 DECIDE = P**3 + Q**2 C IF (DECIDE.GE.0D0) THEN IND = 1 EX = 1D0/3D0 PQ1 = -Q + SQRT(DECIDE) IF (PQ1.GE.0D0) THEN U = PQ1**EX ELSE U = -((-PQ1)**EX) ENDIF PQ2 = -Q - SQRT(DECIDE) IF (PQ2.GE.0D0) THEN V = PQ2**EX ELSE V = -((-PQ2)**EX) ENDIF Y1 = U + V X(1) = Y1 - A1/3D0 ELSE IND = 3 PSQRT = DSQRT(DABS(P)) PI2 = 2D0/3D0*PI A13 = A1/3D0 PQ3 = Q/P/PSQRT FAI = DACOS(PQ3)/3D0 Y1 = 2D0*PSQRT*DCOS(FAI) Y2 = 2D0*PSQRT*DCOS(FAI + PI2) Y3 = 2D0*PSQRT*DCOS(FAI + 2D0*PI2) X(1) = Y1 - A13 X(2) = Y2 - A13 X(3) = Y3 - A13 ENDIF C IF (IND.EQ.1) RETURN C IF (DABS(X(2)).GT.DABS(X(1))) THEN DAMMY = X(1) X(1) = X(2) X(2) = DAMMY ENDIF IF (DABS(X(3)).GT.DABS(X(1))) THEN DAMMY = X(1) X(1) = X(3) X(3) = DAMMY ENDIF IF (DABS(X(3)).GT.DABS(X(2))) THEN DAMMY = X(2) X(2) = X(3) X(3) = DAMMY ENDIF EPS = DABS(A1)*DELTA IF (DABS(X(3)).GT.EPS) THEN RETURN ELSEIF (DABS(X(2)).GT.EPS.AND.DABS(X(3)).LT.EPS) THEN X(3) = -A3/X(1)/X(2) ELSEIF (DABS(X(2)).LT.EPS.AND.DABS(X(3)).LT.EPS) THEN AAA = A1 + A2/X(1) + A3/X(1)**2 AA(1) = (A2 + A3/X(1))/AAA AA(2) = A3/AAA CALL FMFN01(AA, XX, INDD) X(2) = XX(1) X(3) = XX(2) ENDIF RETURN END C C*************************************************************************** C SUBROUTINE FMFN01(A, X, IND) IMPLICIT DOUBLE PRECISION (A-H, O-Z) DIMENSION A(2), X(2) A1 = A(1) A2 = A(2) DECIDE = A1**2 - 4D0*A2 IF (DECIDE.LT.0D0) THEN IND = 0 RETURN ENDIF IND = 2 IF (A1.GE.0D0) THEN X(1) = (-A1 - DSQRT(DECIDE))/2D0 ELSE X(1) = (-A1 + DSQRT(DECIDE))/2D0 ENDIF X(2) = A2/X(1) RETURN END C C*********************************************************** C C************************************************************* C INTEGER FUNCTION FMF002(T,P,Z) DOUBLE PRECISION T,P,Z,PS,TC,PC,VC,LIM,TC1,TC2,PC1,PC2 INTEGER JC,FMF003,COMBI DOUBLE PRECISION PR(1:11),CP(1:20) COMMON/FMFC/PR,CP,COMBI C LIM=1.0D-10 PC1=PR(3) PC2=PR(8) TC1=PR(2) TC2=PR(7) C C---- PURE SUBSTANCE -------------------------------------- IF(Z.LT.LIM.OR.Z.GT.1.0-LIM) THEN IF(Z.LT.LIM) THEN IF (P.GT.PC2) THEN FMF002=4 RETURN ELSEIF (T.GT.TC2) THEN FMF002=3 RETURN ENDIF CALL FMF059(2,J,T,PS) ELSE IF (T.GT.PC1) THEN FMF002=4 RETURN ELSEIF (T.GT.TC1) THEN FMF002=3 RETURN ENDIF CALL FMF059(1,J,T,PS) ENDIF IF(J.EQ.0) THEN IF(P.GT.PS) FMF002=1 IF(P.LT.PS) FMF002=3 RETURN ELSE FMF002=-1 RETURN ENDIF ENDIF C C---- MIXTURE --------------------------------------------- C C modified by Akasaka, June 1, 1998 C C CALL FMF074(1,JC,Z,TC,PC,VC) C IF(JC.NE.0) THEN C FMF002=-1 C RETURN C ENDIF C IF(P.GT.PC) THEN C FMF002=4 C RETURN C ENDIF C IF (T.GT.TC) THEN C FMF002=3 C RETURN C ENDIF FMF002=FMF003(T,P,Z) RETURN END C C******************************************************** C INTEGER FUNCTION FMF003(T,P,Z) DOUBLE PRECISION T,P,Z,FMF026,FMF031,PB,PD c PB=FMF026(T,Z) PD=FMF031(T,Z) IF(PB.LT.0.0D0.OR.PD.LT.0.0D0) THEN FMF003=-1 RETURN ELSE IF(P.GT.PB) FMF003=1 IF(P.LE.PB.AND.P.GE.PD) FMF003=2 IF(P.LT.PD) FMF003=3 RETURN ENDIF END C C************************************************************ C SUBROUTINE FMF004(IND,A,XR1,XR2,XI) DOUBLE PRECISION A(3),XR1,XR2,XI,X1,X2,DISC INTEGER IND c X1=-A(2)/(2.0*A(1)) DISC=A(2)*A(2)-4.0*A(1)*A(3) IF(DISC.LT.0.0) THEN IND=0 X2=SQRT(-DISC) XR1=X1 XR2=X1 XI=X2/(2.0*A(1)) ELSEIF(DISC.EQ.0.0) THEN IND=1 XR1=X1 XR2=X1 XI=0.0 ELSEIF(DISC.GT.0.0) THEN IND=2 X2=SQRT(DISC)/(2.0*A(1)) XR1=X1+X2 XR2=X1-X2 XI=0.0 ENDIF RETURN END C C********************************************************** C SUBROUTINE FMF005(C4D,XR,XI) DOUBLE PRECISION C4D(5),CB(3),CBX(3),CQ(3),A0,A1 $ ,A2,A3,A,B,C,D,Z,XR(4),XI(4) INTEGER I1,I2 C A3=C4D(2)/C4D(1) A2=C4D(3)/C4D(1) A1=C4D(4)/C4D(1) A0=C4D(5)/C4D(1) C CB(1)=-A2 CB(2)=A1*A3-4.0*A0 CB(3)=A0*(4.0*A2-A3**2)-A1**2 CALL FMF001(CB,CBX,IND) C IF(IND.EQ.1) THEN Z=CBX(1) ELSEIF(IND.EQ.3) THEN Z=MAX(CBX(1),CBX(2),CBX(3)) ENDIF C B=Z/2.0 A=A3/2.0 DISC=B*B-A0 IF(DISC.GT.0.0) THEN D=SQRT(DISC) C=(A*B-A1/2.0)/D ELSE D=0 C=SQRT(A*A-A2+Z) ENDIF C CQ(1)=1.0 CQ(2)=A-C CQ(3)=B-D CALL FMF004(I1,CQ,XR(1),XR(2),XI(1)) CQ(1)=1.0 CQ(2)=A+C CQ(3)=B+D CALL FMF004(I2,CQ,XR(3),XR(4),XI(3)) XI(2)=-XI(1) XI(4)=-XI(3) RETURN END C C************************************************************ C C------------------------------------- C (TC,PC,OMEGA,T,P) --> (A,A',B) C------------------------------------- SUBROUTINE FMF006(TC,PC,OMEGA,T,P,A,AD,B) DOUBLE PRECISION TC,PC,OMEGA,T,P,R,SA,SAD,SB,A,AD,B c R=8.31451 CALL FMF007(TC,PC,OMEGA,T,SA,SAD,SB) A=SA*P/(R*T)**2 AD=SAD*P/(R*T)**2 B=SB*P/(R*T) RETURN END C C*************************************************************** C C------------------------------------- C (TC,PC,OMEGA,T) --> (a,a',b) C------------------------------------- SUBROUTINE FMF007(TC,PC,OMEGA,T,SA,SAD,SB) DOUBLE PRECISION TC,PC,OMEGA,T,R,K,ALFA,SA,SAD,SB c R=8.31451 K=0.37464+1.54226*OMEGA-0.26992*OMEGA**2 ALFA=(1+K*(1-DSQRT(T/TC)))**2 SA=0.45723552892138*(R**2)*(TC**2)/PC*ALFA SAD=-0.45723552892138*(R**2)*(TC**2)/PC*K*DSQRT(ALFA/(T*TC)) SB=0.077796073903888*R*TC/PC RETURN END C C************************************************************ C C----------------------------------- C (A,A',B,Y) --> (Am,Am',Bm) C----------------------------------- SUBROUTINE FMF008(A,AD,B,Y,AM,ADM,BM) DOUBLE PRECISION A(2),AD(2),B(2),Y,AM,ADM,BM,IAP DOUBLE PRECISION PR(1:11),CP(1:20) INTEGER COMBI COMMON/FMFC/PR,CP,COMBI c IAP=PR(11) AM=Y**2*A(1)+2.0*Y*(1.0-Y)*(1.0-IAP)*DSQRT(A(1)*A(2)) & +(1.0-Y)**2*A(2) ADM=Y**2*AD(1)+2.0*Y*(1.0-Y)*(1.0-IAP)*(AD(1)*A(2)+A(1)*AD(2)) & /(2.0*DSQRT(A(1)*A(2)))+(1.0-Y)**2*AD(2) BM=Y*B(1)+(1.0-Y)*B(2) RETURN END C C************************************************************** C C-------------------------------------- C (A,A',B,Y1,Y2) --> (Am,Am',Bm) C-------------------------------------- SUBROUTINE FMF009(A,AD,B,Y1,Y2,AM,ADM,BM) DOUBLE PRECISION A(2),AD(2),B(2),Y1,Y2,AM,ADM,BM,IAP DOUBLE PRECISION PR(1:11),CP(1:20) INTEGER COMBI COMMON/FMFC/PR,CP,COMBI c IAP=PR(11) AM=Y1**2*A(1)+2.0*Y1*Y2*(1.0-IAP)*DSQRT(A(1)*A(2))+Y2**2*A(2) ADM=Y1**2*AD(1)+2.0*Y1*Y2*(1.0-IAP)*(AD(1)*A(2)+A(1)*AD(2)) & /(2.0*DSQRT(A(1)*A(2)))+Y2**2*AD(2) BM=Y1*B(1)+Y2*B(2) RETURN END C C************************************************************* C C------------------------ C (T,Y) --> (am,bm) C------------------------ SUBROUTINE FMF010(T,Y,SAM,SBM) DOUBLE PRECISION T,Y,SAM,SBM,TC1,TC2,PC1,PC2 $ ,OMEGA1,OMEGA2,KIJ,SA(2),SAD(2),SB(2) DOUBLE PRECISION PR(1:11),CP(1:20) INTEGER COMBI COMMON/FMFC/PR,CP,COMBI C TC1=PR(2) PC1=PR(3) OMEGA1=PR(5) TC2=PR(7) PC2=PR(8) OMEGA2=PR(10) KIJ=PR(11) C CALL FMF007(TC1,PC1,OMEGA1,T,SA(1),SAD(1),SB(1)) CALL FMF007(TC2,PC2,OMEGA2,T,SA(2),SAD(2),SB(2)) SAM=Y**2*SA(1)+2.0*Y*(1.0-Y)*(1.0-KIJ)*DSQRT(SA(1)*SA(2)) & +(1.0-Y)**2*SA(2) SBM=Y*SB(1)+(1.0-Y)*SB(2) RETURN END C C************************************************************ C C-------------------------------- C (T,P,Y) --> (Am,Am',Bm) C-------------------------------- SUBROUTINE FMF011(T,P,Y,AM,ADM,BM) DOUBLE PRECISION T,P,Y,R,TC1,PC1,OMEGA1 $ ,TC2,PC2,OMEGA2,A(2),AD(2),B(2),AM,ADM,BM DOUBLE PRECISION PR(1:11),CP(1:20) INTEGER COMBI COMMON/FMFC/PR,CP,COMBI C R=8.31451 TC1=PR(2) PC1=PR(3) OMEGA1=PR(5) TC2=PR(7) PC2=PR(8) OMEGA2=PR(10) C CALL FMF006(TC1,PC1,OMEGA1,T,P,A(1),AD(1),B(1)) CALL FMF006(TC2,PC2,OMEGA2,T,P,A(2),AD(2),B(2)) CALL FMF008(A,AD,B,Y,AM,ADM,BM) RETURN END C C**************************************************************** C SUBROUTINE FMF012(I,J,Z,TC,PC,VC) DOUBLE PRECISION X(3),X2(2),DOM(2) DOUBLE PRECISION TI,VI,TC,VC,VCC,Z,PC,EP,ER,AM,BM,FMF018 INTEGER I,J,KPA,MESS,KSTAN,KAS,COMBI DOUBLE PRECISION PR(1:11),CP(1:20) COMMON/UNIT/KPA,MESS,KSTAN,KAS COMMON/FMFC/PR,CP,COMBI C IF(I.EQ.1) THEN CALL FMF015(Z,TI,VI,PR) X(1)=TI X(2)=VI X(3)=Z DOM(1)=1.0D-2 DOM(2)=1.0D-4 EP=1.0D-1 CALL FMF016(J,I,X,X2,DOM,EP,PR) IF(J.NE.0) THEN J=-1 ER=-1.0D10 TC=ER PC=ER VC=ER RETURN ELSE J=0 TC=X2(1) VC=X2(2) CALL FMF010(TC,Z,AM,BM) VCC=VC*1.0D-3 PC=FMF018(TC,VCC,AM,BM) RETURN ENDIF C ELSEIF(I.EQ.2) THEN Z=0.5D0 CALL FMF013(J,TC,X,PR) IF (J.NE.0) THEN J=-1 ER=-1.0D10 Z=ER PC=ER VC=ER RETURN ELSE J=0 Z=X(1) VC=X(2) CALL FMF010(TC,Z,AM,BM) VCC=VC*1.0D-3 PC=FMF018(TC,VCC,AM,BM) RETURN ENDIF ENDIF END C C ------------------------------------------------------------ SUBROUTINE FMF013(J,T,X,PR) INTEGER J,TIME DOUBLE PRECISION PR(1:11),X(1:2),Y(1:2) + ,EP,LAM1,LAM2,T1,T2,DT,T C EP=1.0D-7 LAM1=1.0D0 LAM2=0.0D0 T1=PR(2) T2=PR(7) DO 10 TIME=1,100 LAM2=(LAM1-LAM2)*(T-T2)/(T1-T2)+LAM2 CALL FMF014(LAM2,Y,PR) T2=Y(1) DT=((T2-T)/T)**2 IF (DT.LT.EP) THEN X(1)=LAM2 X(2)=Y(2) RETURN ENDIF 10 CONTINUE J=-1 RETURN END C ----------------------------------------------------------- SUBROUTINE FMF014(L,Y,PR) INTEGER I,J DOUBLE PRECISION PR(1:11),X(1:3),DOM(1:2),X2(1:2),Y(1:2) + ,L,TI,VI,EP C CALL FMF015(L,TI,VI,PR) C X(1)=TI X(2)=VI X(3)=L DOM(1)=1.0D-2 DOM(2)=1.0D-4 EP=1.0D-1 J=0 I=1 CALL FMF016(J,I,X,X2,DOM,EP,PR) IF (J.NE.0) THEN RETURN ELSE Y(1)=X2(1) Y(2)=X2(2) RETURN ENDIF END C----------------------------------------------------------- SUBROUTINE FMF015(Z,TI,VI,PR) INTEGER MESP,MESL,MESPM,MESLM,I,M,N DOUBLE PRECISION Z,TI,VI,TA,VAL,VALM,FMF017 + ,PR(1:11),U(1:3),TM(1:100),VM(1:100),GIB(1:100,1:100) C U(3)=Z TA=Z*PR(2)+(1.0D0-Z)*PR(7) IF (PR(2).LT.PR(7)) THEN IF (TA.LT.(PR(2)+5.0D0)) THEN DO 10 I=1,50 TM(I)=PR(2)+DBLE(I)*1.0D-1 TM(I+50)=PR(2)+5.0D0+DBLE(I)*5.0D-1 10 CONTINUE ELSEIF (TA.GT.(PR(7)-5.0D0)) THEN DO 20 I=1,50 TM(I)=PR(7)-DBLE(I)*1.0D-1 TM(I+50)=PR(7)-5.0D0-DBLE(I)*5.0D-1 20 CONTINUE ELSEIF (TA.LT.(PR(2)+1.5D1)) THEN DO 11 I=1,50 TM(I)=PR(2)+DBLE(I)*2.0D-1 TM(I+50)=PR(2)+1.0D1+DBLE(I)*6.0D-1 11 CONTINUE ELSEIF (TA.GT.(PR(7)-1.5D1)) THEN DO 21 I=1,50 TM(I)=PR(7)-DBLE(I)*2.0D-1 TM(I+50)=PR(7)-1.0D1-DBLE(I)*6.0D-1 21 CONTINUE ELSEIF (TA.LT.(PR(2)+3.0D1)) THEN DO 12 I=1,50 TM(I)=PR(2)+DBLE(I)*1.0D0 TM(I+50)=PR(2)+3.0D1+DBLE(I)*1.0D0 12 CONTINUE ELSEIF (TA.GT.(PR(7)-3.0D1)) THEN DO 22 I=1,50 TM(I)=PR(7)-DBLE(I)*1.0D0 TM(I+50)=PR(7)-3.0D1-DBLE(I)*1.0D0 22 CONTINUE ELSE DO 30 I=1,100 TM(I)=TA-3.0D1+DBLE(I)*1.0D0 30 CONTINUE ENDIF ELSE IF (TA.LT.(PR(7)+5.0D0)) THEN DO 40 I=1,50 TM(I)=PR(7)+DBLE(I)*1.0D-1 TM(I+50)=PR(7)+5.0D0+DBLE(I)*5.0D-1 40 CONTINUE ELSEIF (TA.GT.(PR(2)-5.0D0)) THEN DO 50 I=1,50 TM(I)=PR(2)-DBLE(I)*1.0D-1 TM(I+100)=PR(2)-1.0D0-DBLE(I)*5.0D-1 50 CONTINUE ELSEIF (TA.LT.(PR(7)+1.5D1)) THEN DO 41 I=1,50 TM(I)=PR(7)+DBLE(I)*2.0D-1 TM(I+50)=PR(7)+1.0D1+DBLE(I)*6.0D-1 41 CONTINUE ELSEIF (TA.GT.(PR(2)-1.5D1)) THEN DO 51 I=1,50 TM(I)=PR(2)-DBLE(I)*2.0D-1 TM(I+50)=PR(2)-1.0D1-DBLE(I)*6.0D-1 51 CONTINUE ELSEIF (TA.LT.(PR(7)+3.0D1)) THEN DO 42 I=1,50 TM(I)=PR(7)+DBLE(I)*6.0D-1 TM(I+50)=PR(7)+3.0D1+DBLE(I)*1.0D0 42 CONTINUE ELSEIF (TA.GT.(PR(2)-3.0D1)) THEN DO 52 I=1,50 TM(I)=PR(2)-DBLE(I)*6.0D-1 TM(I+50)=PR(2)-3.0D1-DBLE(I)*1.0D0 52 CONTINUE ELSE DO 60 I=1,100 TM(I)=TA-3.0D1+DBLE(I)*1.0D0 60 CONTINUE ENDIF ENDIF C DO 70 MESP=1,99 DO 80 MESL=1,100 VM(MESL)=(1.0D-1)+DBLE(MESL)*6.00D-3 U(1)=TM(MESP) U(2)=VM(MESL) GIB(MESP,MESL)=FMF017(1,U,PR) 80 CONTINUE 70 CONTINUE VALM=GIB(1,1) MESPM=0 MESLM=0 DO 90 M=1,99 DO 100 N=1,100 VAL=GIB(M,N) IF (VAL.LT.VALM) THEN VALM=VAL MESPM=M MESLM=N ENDIF 100 CONTINUE 90 CONTINUE TI=TM(MESPM) VI=VM(MESLM) RETURN END C----------------------------------------------------------- C Complex法により2元連立方程式 C d2G/dx2=0,d3G/dx3=0 C の解を求めます。 C N:number of parameters C P:initial data C P2:result C DOM:domain of section C EP:condition of convergence C JJ:error code C----------------------------------------------------------- SUBROUTINE FMF016(JJ,N,P,P2,DOM,EP,PR) INTEGER I,J,JH,JL,JJ,LIM,NN,TIME,M,LAMBDA,MU,IR,N PARAMETER (NN=2) PARAMETER (M=1664501,LAMBDA=1229,MU=351750) DOUBLE PRECISION DOM(1:NN),P(1:NN+1),X(1:2*NN,1:NN) + ,R(1:NN*(NN*2-1)),P1(1:NN+1),G(1:NN+1) + ,P2(1:NN),XN(1:NN+1),Y(1:2*NN),PR(1:11) DOUBLE PRECISION ALFA,ALFA2,BETA,EP,YH,YL,YG,YN,FMF017,INVM C ALFA=1.3D0 BETA=0.5D0 IR=0 INVM=1.0D0/DBLE(M) LIM=100000 TIME=0 G(3)=P(3) P1(3)=P(3) XN(3)=P(3) C Y(1)=FMF017(N,P,PR) X(1,1)=P(1) X(1,2)=P(2) C DO 170 I=1,6 IR=MOD(LAMBDA*IR+MU,M) R(I)=DBLE(IR)*INVM 170 CONTINUE DO 20 J=2,4 X(J,1)=P(1)+DOM(1)*(R(2*J-3)-0.5D0) X(J,2)=P(2)+DOM(2)*(R(2*J-2)-0.5D0) 20 CONTINUE C DO 40 J=2,4 P1(1)=X(J,1) P1(2)=X(J,2) Y(J)=FMF017(N,P1,PR) C WRITE(*,*) Y(J) 40 CONTINUE C 1000 YH=Y(1) JH=1 YL=Y(1) JL=1 DO 60 J=2,4 IF (Y(J).GT.YH) THEN YH=Y(J) JH=J ELSEIF (Y(J).LT.YL) THEN YL=Y(J) JL=J ENDIF 60 CONTINUE C TIME=TIME+1 C WRITE(*,*) ' TIME=',TIME IF (TIME.GT.LIM) THEN JJ=-1 RETURN ENDIF IF (YL.LT.EP) THEN P2(1)=X(JL,1) P2(2)=X(JL,2) JJ=0 RETURN ENDIF C G(1)=0.0D0 G(2)=0.0D0 DO 80 I=1,2 DO 90 J=1,4 G(I)=G(I)+X(J,I) 90 CONTINUE 80 CONTINUE G(1)=(G(1)-X(JH,1))/3.0D0 G(2)=(G(2)-X(JH,2))/3.0D0 YG=FMF017(N,G,PR) C WRITE(*,*) YG C ALFA2=ALFA IF (YG.LT.YH) THEN C 3000 XN(1)=G(1)+ALFA2*(G(1)-X(JH,1)) XN(2)=G(2)+ALFA2*(G(2)-X(JH,2)) C WRITE(*,*) ' 3' YN=FMF017(N,XN,PR) C WRITE(*,*) YN IF (YN.LT.YH) THEN X(JH,1)=XN(1) X(JH,2)=XN(2) Y(JH)=YN GOTO 1000 ELSE ALFA2=ALFA2*BETA GOTO 3000 ENDIF ELSE C 4000 XN(1)=X(JL,1)+ALFA2*(X(JL,1)-X(JH,1)) XN(2)=X(JL,2)+ALFA2*(X(JL,2)-X(JH,2)) C WRITE(*,*) ' 4' YN=FMF017(N,XN,PR) C WRITE(*,*) YN IF (YN.LT.YH) THEN X(JH,1)=XN(1) X(JH,2)=XN(2) Y(JH)=YN GOTO 1000 ELSE ALFA2=ALFA2*BETA GOTO 4000 ENDIF ENDIF C END C---------------------------------------------------------- C 混合媒体の臨界点はギブスの自由エネルギーを C 混合比で2回微分したものと3回微分したものが C 同時に0になるという条件を満たします。 C C この関数はx1=x、x2=1-xとして、xで微分を行っ C たものです。 C---------------------------------------------------------- C DOUBLE PRECISION FUNCTION FMF017(I,Z,PROP) INTEGER I DOUBLE PRECISION Z(1:3),PROP(1:11) DOUBLE PRECISION T,V,X,R,TC1,TC2,PC1,PC2,OMEGA1,OMEGA2,K12 + ,ALPHA,BETA,KAPPA1,KAPPA2,TR1,TR2,RT,A1,A2 + ,B1,B2,A,B,DAX,D2AX,DBX,INVX,D2H1X,D3H1X,INV1X + ,D2H2X,D3H2X,INVVB,D2H3X,D3H3X,INVB,Y0,DY0X + ,D2Y0X,D3Y0X,C1,C2,INVC1,INVC2,Y1,DY1X,D2Y1X + ,D3Y1X,Y2,DY2X,D2Y2X,D3Y2X,INV2,D2H4X,D3H4X + ,D2H5X,D3H5X,D2HX,D3HX,P2,INVP2,DP1X,D2P1X + ,DP1V,D2P1V,D2P1XV,DP2X,D2P2X,DP2V,D2P2V,D2P2XV + ,DP3X,D2P3X,DP3V,D2P3V,D2P3XV,DPX,D2PX,DPV + ,D2PV,D2PXV,D2GX,D3GX C IF (I.EQ.1) THEN T=Z(1) V=Z(2) X=Z(3) ELSEIF (I.EQ.2) THEN X=Z(1) V=Z(2) T=Z(3) ENDIF C R=8.31451D3 TC1=PROP(2) TC2=PROP(7) PC1=PROP(3) PC2=PROP(8) OMEGA1=PROP(5) OMEGA2=PROP(10) K12=PROP(11) ALPHA=4.5723552892138D-1 BETA=7.7796073903888D-2 KAPPA1=(3.7464D-1)+(1.54226D0)*OMEGA1-(2.6992D-1)*OMEGA1**2 KAPPA2=(3.7464D-1)+(1.54226D0)*OMEGA2-(2.6992D-1)*OMEGA2**2 TR1=T/TC1 TR2=T/TC2 RT=R*T C A1=ALPHA*R**2*TC1**2*(1.0D0+KAPPA1*(1.0D0-DSQRT(TR1)))**2/PC1 A2=ALPHA*R**2*TC2**2*(1.0D0+KAPPA2*(1.0D0-DSQRT(TR2)))**2/PC2 B1=BETA*R*TC1/PC1 B2=BETA*R*TC2/PC2 A=A1*X**2+2.0D0*X*(1.0D0-X)*(1.0D0-K12)*DSQRT(A1*A2)+ + A2*(1.0D0-X)**2 B=X*B1+(1.0D0-X)*B2 DAX=2.0D0*A1*X+2.0D0*(1.0D0-2.0D0*X)*(1.0D0-K12)*DSQRT(A1*A2)+ + 2.0D0*A2*(X-1.0D0) D2AX=2.0D0*A1-4.0D0*(1.0D0-K12)*DSQRT(A1*A2)+2.0D0*A2 DBX=B1-B2 INVX=1.0D0/X D2H1X=RT*INVX D3H1X=-RT*INVX**2 INV1X=1.0D0/(1.0D0-X) D2H2X=RT*INV1X D3H2X=RT*INV1X**2 INVVB=1.0D0/(V-B) D2H3X=RT*DBX**2*INVVB**2 D3H3X=2.0D0*RT*DBX**3*INVVB**3 INVB=1.0D0/B Y0=A*INVB DY0X=(DAX*B-DBX*A)*INVB**2 D2Y0X=D2AX*INVB-2.0D0*DAX*DBX*INVB**2+2.0D0*A*DBX**2*INVB**3 D3Y0X=-3.0D0*D2AX*DBX*INVB**2+6.0D0*DAX*DBX**2*INVB**3- + 6.0D0*A*DBX**3*INVB**4 C1=1.0D0+DSQRT(2.0D0) C2=1.0D0-DSQRT(2.0D0) INVC1=1.0D0/(V+C1*B) INVC2=1.0D0/(V+C2*B) Y1=DLOG(V+C1*B) DY1X=C1*DBX*INVC1 D2Y1X=-(C1**2)*DBX**2*INVC1**2 D3Y1X=2.0D0*C1**3*DBX**3*INVC1**3 Y2=DLOG(V+C2*B) DY2X=C2*DBX*INVC2 D2Y2X=-(C2**2)*DBX**2*INVC2**2 D3Y2X=2.0D0*C2**3*DBX**3*INVC2**3 INV2=1.0D0/(2.0D0*DSQRT(2.0D0)) D2H4X=-INV2*(D2Y0X*Y1+2.0D0*DY0X*DY1X+Y0*D2Y1X) D3H4X=-INV2*(D3Y0X*Y1+3.0D0*D2Y0X*DY1X+3.0D0*DY0X*D2Y1X+ + Y0*D3Y1X) D2H5X=INV2*(D2Y0X*Y2+2.0D0*DY0X*DY2X+Y0*D2Y2X) D3H5X=INV2*(D3Y0X*Y2+3.0D0*D2Y0X*DY2X+3.0D0*DY0X*D2Y2X+ + Y0*D3Y2X) D2HX=D2H1X+D2H2X+D2H3X+D2H4X+D2H5X D3HX=D3H1X+D3H2X+D3H3X+D3H4X+D3H5X C P2=V**2+2.0D0*V*B-B**2 INVP2=1.0D0/P2 DP1X=DBX*RT*INVVB**2 D2P1X=2.0D0*DBX**2*RT*INVVB**3 DP1V=-RT*INVVB**2 D2P1V=2.0D0*RT*INVVB**3 D2P1XV=-2.0D0*DBX*RT*INVVB**3 DP2X=2.0D0*V*DBX-2.0D0*B*DBX D2P2X=-2.0D0*DBX**2 DP2V=2.0D0*(V+B) D2P2V=2.0D0 D2P2XV=2.0D0*DBX DP3X=DAX*INVP2-A*DP2X*INVP2**2 D2P3X=D2AX*INVP2-(2.0D0*DAX*DP2X+A*D2P2X)*INVP2**2+ + 2.0D0*A*DP2X**2*INVP2**3 DP3V=-A*DP2V*INVP2**2 D2P3V=-A*D2P2V*INVP2**2+2.0D0*A*DP2V**2*INVP2**3 D2P3XV=-(DAX*DP2V+A*D2P2XV)*INVP2**2+2.0D0*A*DP2X*DP2V*INVP2**3 DPX=DP1X-DP3X D2PX=D2P1X-D2P3X DPV=DP1V-DP3V D2PV=D2P1V-D2P3V D2PXV=D2P1XV-D2P3XV D2GX=D2HX+DPX**2/DPV D3GX=D3HX+3.0D0*DPX*D2PX/DPV-3.0D0*DPX**2*D2PXV/DPV**2+ + DPX**3*D2PV/DPV**3 FMF017=D2GX**2+D3GX**2 C WRITE(*,*) FMF017 RETURN END C C*********************************************************** C DOUBLE PRECISION FUNCTION FMF018(T,V,A,B) DOUBLE PRECISION T,V,A,B,R c R=8.31451 FMF018=R*T/(V-B)-A/(V*(V+B)+B*(V-B)) RETURN END C C************************************************************** C SUBROUTINE FMF019(FUN) CHARACTER FUN*6 INTEGER KPA,MESS,KSTAN,KAS COMMON/UNIT/KPA,MESS,KSTAN,KAS c IF(MESS.NE.0) THEN WRITE(*,100)'**** NO CONVERGENCE AT ',FUN,' ****' ENDIF 100 FORMAT(3A) RETURN END C C*************************************************************** C SUBROUTINE FMF020(FUN) CHARACTER FUN*6 INTEGER KPA,MESS,KSTAN,KAS COMMON/UNIT/KPA,MESS,KSTAN,KAS c IF(MESS.NE.0) THEN WRITE(*,100)'**** OUT OF RANGE AT ',FUN,' ****' ENDIF 100 FORMAT(3A) RETURN END C C**************************************************************** C SUBROUTINE FMF021(FUN) CHARACTER FUN*6 INTEGER KPA,MESS,KSTAN,KAS COMMON/UNIT/KPA,MESS,KSTAN,KAS c IF(MESS.NE.0) THEN WRITE(*,100)'**** NO COMBINATION AT ',FUN,' ****' ENDIF 100 FORMAT(3A) RETURN END C C******************************************************** C DOUBLE PRECISION FUNCTION FMF022(I,XIN) REAL XIN DOUBLE PRECISION FMF025 INTEGER KPA,MESS,KSTAN,KAS,I COMMON/UNIT/KPA,MESS,KSTAN,KAS C IF(I.EQ.1) THEN IF(KPA.EQ.1.OR.KPA.EQ.3) THEN FMF022=DBLE(XIN+273.15) ELSE FMF022=DBLE(XIN) ENDIF RETURN ENDIF C ------------------------------------------------- IF(I.EQ.2) THEN IF(KPA.EQ.1.OR.KPA.EQ.2) THEN FMF022=DBLE(XIN*1.0E5) ELSE FMF022=DBLE(XIN) ENDIF RETURN ENDIF C ------------------------------------------------- IF(I.EQ.3) THEN IF(KAS.EQ.1) THEN FMF022=FMF025(DBLE(XIN)) ELSE FMF022=DBLE(XIN) ENDIF RETURN ENDIF END C C********************************************************** C DOUBLE PRECISION FUNCTION FMF023(I,XIN,Z) DOUBLE PRECISION XKMOL,MW,HSTAN,SSTAN,FMF036 REAL XIN,Z,AKMOL,MOL INTEGER KPA,MESS,KSTAN,KAS,I,COMBI DOUBLE PRECISION PR(1:11),CP(1:20) COMMON/UNIT/KPA,MESS,KSTAN,KAS COMMON/FMFC/PR,CP,COMBI C IF(KAS.EQ.1) THEN MOL=AKMOL(Z) MW=FMF036(AKMOL(Z)) XKMOL=DBLE(XIN)*MW ELSE MOL=Z XKMOL=DBLE(XIN) ENDIF C --------------------------------------------------- IF(I.EQ.1) THEN FMF023=XKMOL RETURN ENDIF C --------------------------------------------------- IF(I.EQ.2) THEN IF(KSTAN.EQ.1) THEN HSTAN= & (DBLE(MOL)*CP(9)*PR(1) & +DBLE(1.0-MOL)*CP(19)*PR(6))*1.0D3 FMF023=XKMOL-HSTAN ELSE FMF023=XKMOL ENDIF RETURN ENDIF C --------------------------------------------------- IF(I.EQ.3) THEN IF(KSTAN.EQ.1) THEN SSTAN= & (DBLE(MOL)*CP(10)*PR(1) & +DBLE(1.0-MOL)*CP(20)*PR(6))*1.0D3 FMF023=XKMOL-SSTAN ELSE FMF023=XKMOL ENDIF RETURN ENDIF END C C*********************************************************** C DOUBLE PRECISION FUNCTION FMF024(X) DOUBLE PRECISION X,MW1,MW2 DOUBLE PRECISION PR(1:11),CP(1:20) INTEGER COMBI COMMON/FMFC/PR,CP,COMBI c MW1=PR(1) MW2=PR(6) FMF024=X*MW1/(X*MW1+(1.0-X)*MW2) RETURN END C C************************************************************ C DOUBLE PRECISION FUNCTION FMF025(X) DOUBLE PRECISION X,MW1,MW2 DOUBLE PRECISION PR(1:11),CP(1:20) INTEGER COMBI COMMON/FMFC/PR,CP,COMBI c MW1=PR(1) MW2=PR(6) FMF025=X/MW1/(X/MW1+(1.0-X)/MW2) RETURN END C C********************************************************** C C-------------------------------------------------- C FMF026-02 CALCULATE PRESSURE AT BUBLE POINT. C-------------------------------------------------- DOUBLE PRECISION FUNCTION FMF026(T,Z) DOUBLE PRECISION T,Z,PB,FMF027,FMF028,LIM c LIM=1.0D-7 IF(Z.LT.LIM) THEN CALL FMF059(2,J,T,PB) ELSEIF(Z.GT.1.0-LIM) THEN CALL FMF059(1,J,T,PB) ELSE PB=FMF027(T,Z) IF(PB.LE.0.0) PB=FMF028(T,Z) ENDIF FMF026=PB RETURN END C C************************************************************ C DOUBLE PRECISION FUNCTION FMF027(T,X) DOUBLE PRECISION T,X,TC1,PC1,OMEGA1,TC2,PC2 $ ,OMEGA2,P,A(2), & AD(2),B(2),P0,Y1,Y2,X1,X2,AMY,ADMY,BMY,ZL,ZV,ZLM,ZVM, & SKX1,SKX2,AMX,ADMX,BMX,K1,K2,FMF035,FL1,FL2,FV1,FV2, & FMF030 INTEGER J,N,NN,COMBI DOUBLE PRECISION PR(1:11),CP(1:20) COMMON/FMFC/PR,CP,COMBI C TC1=PR(2) PC1=PR(3) OMEGA1=PR(5) TC2=PR(7) PC2=PR(8) OMEGA2=PR(10) P0=FMF030(T,X) X1=X X2=1.0-X C P=P0 DO 2000 N=1,2000 CALL FMF006(TC1,PC1,OMEGA1,T,P,A(1),AD(1),B(1)) CALL FMF006(TC2,PC2,OMEGA2,T,P,A(2),AD(2),B(2)) CALL FMF008(A,AD,B,X,AMX,ADMX,BMX) CALL FMF073(AMX,BMX,ZV,ZLM) IF(ZLM.LT.0.0) THEN FMF027=100.0 RETURN ENDIF FL1=FMF035(1,0,P,A,B,AMX,BMX,X1,ZLM) FL2=FMF035(2,0,P,A,B,AMX,BMX,X2,ZLM) C Y1=1.0D-10 Y2=1.0-X1 SKX1=0.0 C DO 1000 NN=1,200 CALL FMF009(A,AD,B,Y1,Y2,AMY,ADMY,BMY) CALL FMF073(AMY,BMY,ZVM,ZL) FV1=FMF035(1,0,P,A,B,AMY,BMY,Y1,ZVM) FV2=FMF035(2,0,P,A,B,AMY,BMY,Y2,ZVM) C K1=FL1/FV1 K2=FL2/FV2 C SKX2=X1*K1+X2*K2 IF(DABS(SKX1/SKX2-1.0).LT.1.0D-7) GOTO 1500 Y1=(X1*K1)/SKX2 Y2=(X2*K2)/SKX2 SKX1=SKX2 1000 CONTINUE C 1500 IF(DABS(SKX2-1.0).LT.1.0D-7) THEN J=0 GOTO 2500 ENDIF P=P*SKX2 2000 CONTINUE J=-1 C 2500 IF(J.EQ.0) THEN FMF027=P ELSE FMF027=-1.0D10 ENDIF RETURN END C C*********************************************************** C DOUBLE PRECISION FUNCTION FMF028(T,Z) DOUBLE PRECISION T,Z,X,Y,TC1,P1,P2,TC,VC,P DOUBLE PRECISION PR(1:11),CP(1:20) INTEGER COMBI COMMON/FMFC/PR,CP,COMBI C TC1=PR(2) IF(T.LE.TC1) THEN CALL FMF059(1,J1,T,P1) ELSE CALL FMF074(1,J1,Z,TC,P1,VC) ENDIF CALL FMF059(2,J2,T,P2) C DO 1000 N=1,200 P=DEXP(DLOG(P1*P2)*0.5) CALL FMF044(JJ,T,P,X,Y) IF(DABS(X/Z-1.0).LT.1.0E-7) THEN J=0 GOTO 2000 ENDIF IF(X.GT.Z) THEN P1=P ELSE P2=P ENDIF 1000 CONTINUE J=-1 2000 IF(J.EQ.0) THEN FMF028=P ELSE FMF028=-1.0D10 ENDIF RETURN END C C************************************************************** C C---------------------------------------------------- C FMF029 CALCULATES TEMPERATURE AT BUBLE POINT. C---------------------------------------------------- DOUBLE PRECISION FUNCTION FMF029(P,Z) DOUBLE PRECISION P,Z,TC1,PC1,TC2,PC2,T1,T2,P1,P2 $ ,T,FMF026,PB,LIM,TC,PC,VC INTEGER J1,J2,J,COMBI DOUBLE PRECISION PR(1:11),CP(1:20) COMMON/FMFC/PR,CP,COMBI C LIM=1.0D-7 IF(Z.LT.LIM) THEN CALL FMF061(2,J,T,P) FMF029=T RETURN ELSEIF(Z.GT.1.0-LIM) THEN CALL FMF061(1,J,T,P) FMF029=T RETURN ENDIF C TC1=PR(2) PC1=PR(3) TC2=PR(7) PC2=PR(8) IF(P.LE.PC1) THEN CALL FMF061(1,J1,T1,P) ELSE T1=TC1 ENDIF CALL FMF074(1,J2,Z,TC,PC,VC) T2=TC C DO 1000 N=1,5 T=0.5*T1+0.5*T2 PB=FMF026(T,Z) IF(PB.GT.P) THEN T2=T ELSE T1=T ENDIF 1000 CONTINUE C P1=FMF026(T1,Z) P2=FMF026(T2,Z) DO 1500 N=1,50 T=T1+(P-P1)*(T2-T1)/(P2-P1) PB=FMF026(T,Z) IF(DABS(PB/P-1.0).LT.1.0D-7) THEN FMF029=T RETURN ENDIF IF(PB.GT.P) THEN T2=T P2=PB ELSE T1=T P1=PB ENDIF 1500 CONTINUE FMF029=-1.0D10 RETURN END C C*********************************************************** C--------------------------------------------------------------------- C FMF030 CALCULATES INITIAL VALUE OF PRESSURE AT DEW OR BUBLE POINT. C--------------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF030(T,Y) DOUBLE PRECISION T,Y,TC1,PC1,OMEGA1,TC2,PC2,OMEGA2, & SA(2),SAD(2),SB(2),SAM,SADM,SBM,V1,V2,P1,P2 DOUBLE PRECISION PR(1:11),CP(1:20) INTEGER COMBI COMMON/FMFC/PR,CP,COMBI c TC1=PR(2) PC1=PR(3) OMEGA1=PR(5) TC2=PR(7) PC2=PR(8) OMEGA2=PR(10) CALL FMF007(TC1,PC1,OMEGA1,T,SA(1),SAD(1),SB(1)) CALL FMF007(TC2,PC2,OMEGA2,T,SA(2),SAD(2),SB(2)) CALL FMF008(SA,SAD,SB,Y,SAM,SADM,SBM) CALL FMF058(T,SAM,SBM,V1,V2,P1,P2) IF(P1.LE.0.0) P1=1.0 FMF030=P1*0.5+P2*0.5 RETURN END C C********************************************************** C C--------------------------------------------- C FMF031-02 CALCULATES PRESSURE AT DEW POINT. C--------------------------------------------- DOUBLE PRECISION FUNCTION FMF031(T,Z) DOUBLE PRECISION T,Z,PD,FMF033,X,LIM c LIM=1.0D-7 IF(Z.LT.LIM) THEN CALL FMF059(2,J,T,PD) ELSEIF(Z.GT.1.0-LIM) THEN CALL FMF059(1,J,T,PD) ELSE CALL FMF032(J,T,PD,X,Z) IF(J.NE.0.0) PD=FMF033(T,Z) ENDIF FMF031=PD RETURN END C C********************************************************* C SUBROUTINE FMF032(J,T,PD,X,Y) DOUBLE PRECISION T,PD,X,Y,TC1,PC1,OMEGA1,TC2,PC2 $ ,OMEGA2,P,A(2), & AD(2),B(2),P0,Y1,Y2,X1,X2,AMY,ADMY,BMY,ZL,ZV,ZLM,ZVM, & SKY1,SKY2,AMX,ADMX,BMX,K1,K2,FMF035,FL1,FL2,FV1,FV2, & FMF030 INTEGER J,N,NN,COMBI DOUBLE PRECISION PR(1:11),CP(1:20) COMMON/FMFC/PR,CP,COMBI C TC1=PR(2) PC1=PR(3) OMEGA1=PR(5) TC2=PR(7) PC2=PR(8) OMEGA2=PR(10) P0=FMF030(T,Y) Y1=Y Y2=1.0-Y C P=P0 DO 2000 N=1,200 CALL FMF006(TC1,PC1,OMEGA1,T,P,A(1),AD(1),B(1)) CALL FMF006(TC2,PC2,OMEGA2,T,P,A(2),AD(2),B(2)) CALL FMF008(A,AD,B,Y,AMY,ADMY,BMY) CALL FMF073(AMY,BMY,ZVM,ZL) FV1=FMF035(1,0,P,A,B,AMY,BMY,Y1,ZVM) FV2=FMF035(2,0,P,A,B,AMY,BMY,Y2,ZVM) C X1=1.0D-7 X2=1.0-X1 SKY1=0.0 C DO 1000 NN=1,200 CALL FMF009(A,AD,B,X1,X2,AMX,ADMX,BMX) CALL FMF073(AMX,BMX,ZV,ZLM) IF(ZLM.LT.0.0) THEN J=0 PD=100.0 X=0.0 RETURN ENDIF FL1=FMF035(1,0,P,A,B,AMX,BMX,X1,ZLM) FL2=FMF035(2,0,P,A,B,AMX,BMX,X2,ZLM) C K1=FL1/FV1 K2=FL2/FV2 C SKY2=Y1/K1+Y2/K2 IF(DABS(SKY1/SKY2-1.0).LT.1.0D-7) GOTO 1500 X1=(Y1/K1)/SKY2 X2=(Y2/K2)/SKY2 SKY1=SKY2 1000 CONTINUE C 1500 IF(DABS(SKY2-1.0).LT.1.0D-7) THEN J=0 GOTO 2500 ENDIF P=P/SKY2 2000 CONTINUE J=-1 C 2500 IF(J.EQ.0) THEN PD=P X=X1 ELSE PD=-1.0D10 X=-1.0D10 ENDIF RETURN END C C C********************************************************** C DOUBLE PRECISION FUNCTION FMF033(T,Z) DOUBLE PRECISION T,Z,X,Y,TC1,P1,P2,TC,VC,P DOUBLE PRECISION PR(1:11),CP(1:20) INTEGER COMBI COMMON/FMFC/PR,CP,COMBI C TC1=PR(2) IF(T.LE.TC1) THEN CALL FMF059(1,J1,T,P1) ELSE CALL FMF074(1,J1,Z,TC,P1,VC) ENDIF CALL FMF059(2,J2,T,P2) C DO 1000 N=1,200 P=DEXP(DLOG(P1*P2)*0.5) CALL FMF044(JJ,T,P,X,Y) IF(DABS(Y/Z-1.0).LT.1.0E-7) THEN J=0 GOTO 2000 ENDIF IF(Y.GT.Z) THEN P1=P ELSE P2=P ENDIF 1000 CONTINUE J=-1 2000 IF(J.EQ.0) THEN FMF033=P ELSE FMF033=-1.0D10 ENDIF RETURN END C C*************************************************************** C C---------------------------------------------------- C FMF034 CALCULATES TEMPERATURE AT DEW POINT. C---------------------------------------------------- DOUBLE PRECISION FUNCTION FMF034(P,Z) DOUBLE PRECISION P,Z,TC1,PC1,TC2,PC2,T1,T2,P1,P2,T, & FMF031,PD,LIM,TC,PC,VC INTEGER J1,J2,J,COMBI DOUBLE PRECISION PR(1:11),CP(1:20) COMMON/FMFC/PR,CP,COMBI C LIM=1.0D-7 IF(Z.LT.LIM) THEN CALL FMF061(2,J,T,P) FMF034=T RETURN ELSEIF(Z.GT.1.0-LIM) THEN CALL FMF061(1,J,T,P) FMF034=T RETURN ENDIF C TC1=PR(2) PC1=PR(3) TC2=PR(7) PC2=PR(8) IF(P.LE.PC1) THEN CALL FMF061(1,J1,T1,P) ELSE T1=TC1 ENDIF CALL FMF074(1,J2,Z,TC,PC,VC) T2=TC C DO 1000 N=1,5 T=0.5*T1+0.5*T2 PD=FMF031(T,Z) IF(PD.GT.P) THEN T2=T ELSE T1=T ENDIF 1000 CONTINUE C P1=FMF031(T1,Z) P2=FMF031(T2,Z) DO 1500 N=1,50 T=T1+(P-P1)*(T2-T1)/(P2-P1) PD=FMF031(T,Z) IF(DABS(PD/P-1.0).LT.1.0D-7) THEN FMF034=T RETURN ENDIF IF(PD.GT.P) THEN T2=T P2=PD ELSE T1=T P1=PD ENDIF 1500 CONTINUE FMF034=-1.0D10 RETURN END C C*********************************************************************** C C------------------------------------------------------------------- C FMF035 CALCULATES FUGACITY OR FUGACITY COEFFICIENT OF MIXTURE. C------------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF035(I,L,P,A,B,AM,BM,Y,Z) DOUBLE PRECISION P,A(2),B(2),AM,BM,Y,Z $ ,IAP,A12,SIGMA,FAI,FUG INTEGER I,L,COMBI DOUBLE PRECISION PR(1:11),CP(1:20) COMMON/FMFC/PR,CP,COMBI C IAP=PR(11) A12=(1.0-IAP)*DSQRT(A(1)*A(2)) SIGMA=2.0*(Y*A(I)+(1.0-Y)*A12) FAI=B(I)/BM*(Z-1.0)-DLOG(Z-BM) FAI=FAI-(AM/(2.0*DSQRT(2.0D0)*BM)) & *(SIGMA/AM-B(I)/BM)*DLOG((Z+(1.0+DSQRT(2.0D0))*BM)/(Z+ & (1.0-DSQRT(2.0D0))*BM)) FAI=DEXP(FAI) FUG=Y*P*FAI IF(L.EQ.0) THEN FMF035=FAI ELSE FMF035=FUG ENDIF RETURN END C C******************************************************************* C DOUBLE PRECISION FUNCTION FMF036(Z0) DOUBLE PRECISION Z,MW1,MW2 REAL Z0 DOUBLE PRECISION PR(1:11),CP(1:20) INTEGER COMBI COMMON/FMFC/PR,CP,COMBI c MW1=PR(1) MW2=PR(6) Z=DBLE(Z0) FMF036=Z*MW1+(1.0-Z)*MW2 RETURN END C C************************************************************** C C---------------------------------------- C FMF037 CALCULATES PARTIAL ENTHALPY. C---------------------------------------- SUBROUTINE FMF037(IP,J,T,P,X,Y,HL,HV) DOUBLE PRECISION T,P,X,Y,HL(2),HV(2) $ ,R,TC1,TC2,PC1,PC2, & OMEGA1,OMEGA2,KIJ,A(2),AD(2),B(2), & EQ(2),C(2,5),HIG(2),T0, & HR(2),HRL(2),HRV(2), & AM,ADM,BM,Z,ZV,ZL, & AMX,ADMX,BMX, & AMY,ADMY,BMY INTEGER J,K,IP,FMF002,COMBI DOUBLE PRECISION PR(1:11),CP(1:20) COMMON/FMFC/PR,CP,COMBI C R=8.31451 TC1=PR(2) PC1=PR(3) OMEGA1=PR(5) TC2=PR(7) PC2=PR(8) OMEGA2=PR(10) KIJ=PR(11) CALL FMF006(TC1,PC1,OMEGA1,T,P,A(1),AD(1),B(1)) CALL FMF006(TC2,PC2,OMEGA2,T,P,A(2),AD(2),B(2)) C C ----- IDEAL GAS ENTHALPY ------------------------------------ EQ(1)=CP(1) C(1,1)=CP(2) C(1,2)=CP(3) C(1,3)=CP(4) C(1,4)=CP(5) C(1,5)=CP(6) EQ(2)=CP(11) C(2,1)=CP(12) C(2,2)=CP(13) C(2,3)=CP(14) C(2,4)=CP(15) C(2,5)=CP(16) T0=298.15 DO 1000 K=1,2 IF(EQ(K).LT.1.1) THEN HIG(K)=C(K,1)*T+C(K,2)*C(K,3)/DTANH(C(K,3)/T) & -C(K,4)*C(K,5)*DTANH(C(K,5)/T) & -C(K,1)*T0-C(K,2)*C(K,3)/DTANH(C(K,3)/T0) & +C(K,4)*C(K,5)*DTANH(C(K,5)/T0) ELSE HIG(K)=C(K,1)*(T-T0)+C(K,2)*(T**2-T0**2)/2.0+C(K,3) & *(T**3-T0**3)/3.0+C(K,4)*(T**4-T0**4)/4.0 ENDIF 1000 CONTINUE C-------------------------------------------------------------------- IP=FMF002(T,P,Y) C C ---- SINGLE PHASE REGION ------------------------------- IF(IP.NE.2) THEN CALL FMF008(A,AD,B,Y,AM,ADM,BM) CALL FMF073(AM,BM,ZV,ZL) IF(IP.EQ.1) THEN Z=ZL ELSE Z=ZV ENDIF C CALL FMF038(T,A,AD,B,Y,Z,HR) C IF(IP.EQ.1) THEN HL(1)=HR(1)+HIG(1) HL(2)=HR(2)+HIG(2) HV(1)=-1.0D20 HV(2)=-1.0D20 ELSE HL(1)=-1.0D20 HL(2)=-1.0D20 HV(1)=HR(1)+HIG(1) HV(2)=HR(2)+HIG(2) ENDIF X=-1.0D20 J=0 RETURN ENDIF C C C ---- TWO-PHASE REGION ---------------------------------- IF(IP.EQ.2) THEN CALL FMF043(J,T,P,X,Y) IF(J.EQ.-1) THEN ER=-1.0D10 HL(1)=ER HL(2)=ER HV(1)=ER HV(2)=ER RETURN ENDIF CALL FMF008(A,AD,B,X,AMX,ADMX,BMX) CALL FMF073(AMX,BMX,Z,ZL) CALL FMF008(A,AD,B,Y,AMY,ADMY,BMY) CALL FMF073(AMY,BMY,ZV,Z) C CALL FMF038(T,A,AD,B,X,ZL,HRL) CALL FMF038(T,A,AD,B,Y,ZV,HRV) C HL(1)=HRL(1)+HIG(1) HL(2)=HRL(2)+HIG(2) HV(1)=HRV(1)+HIG(1) HV(2)=HRV(2)+HIG(2) J=0 RETURN ENDIF END C C*************************************************************** C SUBROUTINE FMF038(T,A,AD,B,Y,Z,H) DOUBLE PRECISION T,A(2),AD(2),B(2),Y,Z,H(2),R,KIJ, & AM,ADM,BM,AY,BY,ZY,ADY,HR,HRY DOUBLE PRECISION PR(1:11),CP(1:20) INTEGER COMBI COMMON/FMFC/PR,CP,COMBI C R=8.31451 KIJ=PR(11) C CALL FMF008(A,AD,B,Y,AM,ADM,BM) AY=2.0*Y*A(1) & +2.0*(1.0-2.0*Y)*(1.0-KIJ)*DSQRT(A(1)*A(2)) & -2.0*(1.0-Y)*A(2) BY=B(1)-B(2) ZY=(Z*(-BY*Z-AY+6.0*BM*BY+2.0*BY)+BM*AY+AM*BY & -2.0*BM*BY-3.0*BM*BM*BY) & /(3.0*Z*Z+2.0*(-1.0+BM)*Z+(AM-3.0*BM*BM-2.0*BM)) C ADY=2.0*Y*AD(1) & +(1.0-2.0*Y)*(1.0-KIJ)*(A(1)*A(2))**(-0.5) & *(AD(1)*A(2)+A(1)*AD(2)) & -2.0*(1.0-Y)*AD(2) C HR=R*T*(Z-1.0) & +R*T*(T*ADM-AM)/(2.0*DSQRT(2.0D0)*BM) & *DLOG((Z+(1.0+DSQRT(2.0D0))*BM)/(Z+(1.0-DSQRT(2.0D0))*BM)) C HRY=R*T*ZY & +(R*T/(2.0*DSQRT(2.0D0)))*((T*ADY-AY)*BM-(T*ADM-AM)*BY) & /(BM**2)*DLOG((Z+(1.0+DSQRT(2.0D0))*BM) & /(Z+(1.0-DSQRT(2.0D0))*BM)) & +R*T*(T*ADM-AM)/(2.0*DSQRT(2.0D0)*BM) & *((ZY+(1.0+DSQRT(2.0D0))*BY)/(Z+(1.0+DSQRT(2.0D0))*BM) & -(ZY+(1.0-DSQRT(2.0D0))*BY)/(Z+(1.0-DSQRT(2.0D0))*BM)) C H(1)=(HR+(1.0-Y)*HRY)*1.0D3 H(2)=(HR-Y*HRY)*1.0D3 C RETURN END C C******************************************************** C C-------------------------------------- C FMF039 CALCULATES PARTIAL ENTROPY. C-------------------------------------- SUBROUTINE FMF039(IP,J,T,P,X,Y,SL,SV) DOUBLE PRECISION T,P,X,Y,SL(2),SV(2) $ ,R,TC1,TC2,PC1,PC2, & OMEGA1,OMEGA2,KIJ,A(2),AD(2),B(2), & EQ(2),C(2,5),SIG(2),T0,P0, & SR(2),SRL(2),SRV(2), & AM,ADM,BM,Z,ZV,ZL, & AMX,ADMX,BMX, & AMY,ADMY,BMY INTEGER J,K,IP,FMF002,COMBI DOUBLE PRECISION PR(1:11),CP(1:20) COMMON/FMFC/PR,CP,COMBI C TC1=PR(2) PC1=PR(3) OMEGA1=PR(5) TC2=PR(7) PC2=PR(8) OMEGA2=PR(10) KIJ=PR(11) CALL FMF006(TC1,PC1,OMEGA1,T,P,A(1),AD(1),B(1)) CALL FMF006(TC2,PC2,OMEGA2,T,P,A(2),AD(2),B(2)) C C ----- IDEAL GAS ENTROPY ------------------------------------ EQ(1)=CP(1) C(1,1)=CP(2) C(1,2)=CP(3) C(1,3)=CP(4) C(1,4)=CP(5) C(1,5)=CP(6) EQ(2)=CP(11) C(2,1)=CP(12) C(2,2)=CP(13) C(2,3)=CP(14) C(2,4)=CP(15) C(2,5)=CP(16) R=8.31451 RR=8314.51 T0=298.15 P0=1.0E5 C DO 1000 K=1,2 IF(EQ(K).LT.1.1) THEN SIG(K)=C(K,1)*DLOG(T) & +C(K,2)*(C(K,3)/T/DTANH(C(K,3)/T) & -DLOG(DABS(DSINH(C(K,3)/T)))) & -C(K,4)*(C(K,5)/T*DTANH(C(K,5)/T) & -DLOG(DCOSH(C(K,5)/T))) & -C(K,1)*DLOG(T0) & -C(K,2)*(C(K,3)/T0/DTANH(C(K,3)/T0) & -DLOG(DABS(DSINH(C(K,3)/T0)))) & +C(K,4)*(C(K,5)/T0*DTANH(C(K,5)/T0) & -DLOG(DCOSH(C(K,5)/T0))) & -RR*DLOG(P/P0) ELSE SIG(K)=C(K,1)*(DLOG(T/T0))+C(K,2)*(T-T0) & +C(K,3)*(T**2-T0**2)/2.0+C(K,4) & *(T**3-T0**3)/3.0-R*DLOG(P/P0) ENDIF 1000 CONTINUE C-------------------------------------------------------------------- IP=FMF002(T,P,Y) C C ---- SINGLE PHASE REGION ------------------------------- IF(IP.NE.2) THEN CALL FMF008(A,AD,B,Y,AM,ADM,BM) CALL FMF073(AM,BM,ZV,ZL) IF(IP.EQ.1) THEN Z=ZL ELSE Z=ZV ENDIF C CALL FMF040(T,A,AD,B,Y,Z,SR) C IF(IP.EQ.1) THEN SL(1)=SR(1)+SIG(1)-RR*DLOG(Y) SL(2)=SR(2)+SIG(2)-RR*DLOG(1.0-Y) SV(1)=-1.0D20 SV(2)=-1.0D20 ELSE SL(1)=-1.0D20 SL(2)=-1.0D20 SV(1)=SR(1)+SIG(1)-RR*DLOG(Y) SV(2)=SR(2)+SIG(2)-RR*DLOG(1.0-Y) ENDIF X=-1.0D20 J=0 RETURN ENDIF C C C ---- TWO-PHASE REGION ---------------------------------- IF(IP.EQ.2) THEN CALL FMF043(J,T,P,X,Y) IF(J.EQ.-1) THEN ER=-1.0D10 SL(1)=ER SL(2)=ER SV(1)=ER SV(2)=ER RETURN ENDIF CALL FMF008(A,AD,B,X,AMX,ADMX,BMX) CALL FMF073(AMX,BMX,Z,ZL) CALL FMF008(A,AD,B,Y,AMY,ADMY,BMY) CALL FMF073(AMY,BMY,ZV,Z) C CALL FMF038(T,A,AD,B,X,ZL,SRL) CALL FMF038(T,A,AD,B,Y,ZV,SRV) C SL(1)=SRL(1)+SIG(1)-RR*DLOG(X) SL(2)=SRL(2)+SIG(2)-RR*DLOG(1.0-X) SV(1)=SRV(1)+SIG(1)-RR*DLOG(Y) SV(2)=SRV(2)+SIG(2)-RR*DLOG(1.0-Y) J=0 RETURN ENDIF END C C********************************************************** C SUBROUTINE FMF040(T,A,AD,B,Y,Z,S) DOUBLE PRECISION T,A(2),AD(2),B(2),Y,Z,S(2),R,KIJ, & AM,ADM,BM,AY,BY,ZY,ADY,SR,SRY DOUBLE PRECISION PR(1:11),CP(1:20) INTEGER COMBI COMMON/FMFC/PR,CP,COMBI C R=8.31451 KIJ=PR(11) C CALL FMF008(A,AD,B,Y,AM,ADM,BM) AY=2.0*Y*A(1) & +2.0*(1.0-2.0*Y)*(1.0-KIJ)*DSQRT(A(1)*A(2)) & -2.0*(1.0-Y)*A(2) BY=B(1)-B(2) ZY=(Z*(-BY*Z-AY+6.0*BM*BY+2.0*BY)+BM*AY+AM*BY & -2.0*BM*BY-3.0*BM*BM*BY) & /(3.0*Z*Z+2.0*(-1.0+BM)*Z+(AM-3.0*BM*BM-2.0*BM)) C ADY=2.0*Y*AD(1) & +(1.0-2.0*Y)*(1.0-KIJ)*(A(1)*A(2))**(-0.5) & *(AD(1)*A(2)+A(1)*AD(2)) & -2.0*(1.0-Y)*AD(2) C SR=R*LOG(Z-BM) & +R*T*ADM/(2.0*DSQRT(2.0D0)*BM) & *DLOG((Z+(1.0+DSQRT(2.0D0))*BM)/(Z+(1.0-DSQRT(2.0D0))*BM)) C SRY=R*(ZY-BY)/(Z-BM)+R/(2.0*DSQRT(2.0D0)) & *(T*ADY*BM-T*ADM*BY) & /(BM**2)*DLOG((Z+(1.0+DSQRT(2.0D0))*BM) & /(Z+(1.0-DSQRT(2.0D0))*BM)) & +(R*T*ADM)/(2.0*DSQRT(2.0D0)*BM) $ *((ZY+(1.0+DSQRT(2.0D0))*BY) & /(Z+(1.0+DSQRT(2.0D0))*BM)-(ZY+(1.0-DSQRT(2.0D0))*BY) & /(Z+(1.0-DSQRT(2.0D0))*BM)) C S(1)=(SR+(1.0-Y)*SRY)*1.0D3 S(2)=(SR-Y*SRY)*1.0D3 C RETURN END C C********************************************************** C C-------------------------------------- C FMF041 CALCULATES PARTIAL VOLUME. C-------------------------------------- SUBROUTINE FMF041(IP,J,T,P,X,Y,VL,VV) DOUBLE PRECISION T,P,X,Y,VL(2),VV(2),R,TC1,TC2,PC1,PC2, & OMEGA1,OMEGA2,KIJ,A(2),AD(2),B(2),AM,ADM,BM, & Z,ZV,ZL,AY,BY,ZY1,ZY2,ZY,AX,BX,AMX,BMX,ZX1, & ZX2,ZX,AMY,BMY,ER INTEGER J,FMF002,COMBI DOUBLE PRECISION PR(1:11),CP(1:20) COMMON/FMFC/PR,CP,COMBI C R=8.31451 TC1=PR(2) PC1=PR(3) OMEGA1=PR(5) TC2=PR(7) PC2=PR(8) OMEGA2=PR(10) KIJ=PR(11) CALL FMF006(TC1,PC1,OMEGA1,T,P,A(1),AD(1),B(1)) CALL FMF006(TC2,PC2,OMEGA2,T,P,A(2),AD(2),B(2)) IP=FMF002(T,P,Y) C C ---- SINGLE PHASE REGION ------------------------------- IF(IP.NE.2) THEN CALL FMF008(A,AD,B,Y,AM,ADM,BM) CALL FMF073(AM,BM,ZV,ZL) IF(IP.EQ.1) THEN Z=ZL ELSE Z=ZV ENDIF AY=2.0*Y*A(1) & +2.0*(1.0-2.0*Y)*(1.0-KIJ)*DSQRT(A(1)*A(2)) & -2.0*(1.0-Y)*A(2) BY=B(1)-B(2) ZY1=Z*(-BY*Z-AY+6.0*BM*BY+2.0*BY) & +BM*AY+AM*BY & -2.0*BM*BY-3.0*BM*BM*BY ZY2=3.0*Z*Z+2.0*(-1.0+BM)*Z & +(AM-3.0*BM*BM-2.0*BM) ZY=ZY1/ZY2 IF(IP.EQ.1) THEN VL(1)=(R*T/P)*(Z+(1.0-Y)*ZY)*1.0D3 VL(2)=(R*T/P)*(Z-Y*ZY)*1.0D3 X=Y Y=-1.0D20 ELSE VV(1)=(R*T/P)*(Z+(1.0-Y)*ZY)*1.0D3 VV(2)=(R*T/P)*(Z-Y*ZY)*1.0D3 X=-1.0D20 ENDIF J=0 RETURN ENDIF C C ---- TWO-PHASE REGION ---------------------------------- IF(IP.EQ.2) THEN CALL FMF043(J,T,P,X,Y) IF(J.EQ.-1) THEN ER=-1.0D10 VL(1)=ER VL(2)=ER VV(1)=ER VV(2)=ER RETURN ENDIF CALL FMF008(A,AD,B,X,AMX,ADM,BMX) CALL FMF073(AMX,BMX,Z,ZL) CALL FMF008(A,AD,B,Y,AMY,ADM,BMY) CALL FMF073(AMY,BMY,ZV,Z) AX=2.0*X*A(1) & +2.0*(1.0-2.0*X)*(1.0-KIJ)*DSQRT(A(1)*A(2)) & -2.0*(1.0-X)*A(2) BX=B(1)-B(2) ZX1=ZL*(-BX*ZL-AX+6.0*BMX*BX+2.0*BX) & +BMX*AX+AMX*BX & -2.0*BMX*BX-3.0*BMX*BMX*BX ZX2=3.0*ZL*ZL+2.0*(-1.0+BMX)*ZL & +(AMX-3.0*BMX*BMX-2.0*BMX) ZX=ZX1/ZX2 AY=2.0*Y*A(1) & +2.0*(1.0-2.0*Y)*(1.0-KIJ)*DSQRT(A(1)*A(2)) & -2.0*(1.0-Y)*A(2) BY=B(1)-B(2) ZY1=ZV*(-BY*ZV-AY+6.0*BMY*BY+2.0*BY) & +BMY*AY+AMY*BY & -2.0*BMY*BY-3.0*BMY*BMY*BY ZY2=3.0*ZV*ZV+2.0*(-1.0+BMY)*ZV & +(AMY-3.0*BMY*BMY-2.0*BMY) ZY=ZY1/ZY2 VL(1)=(R*T/P)*(ZL+(1.0-X)*ZX)*1.0D3 VL(2)=(R*T/P)*(ZL-X*ZX)*1.0D3 VV(1)=(R*T/P)*(ZV+(1.0-Y)*ZY)*1.0D3 VV(2)=(R*T/P)*(ZV-Y*ZY)*1.0D3 J=0 RETURN ENDIF END C C****************************************************** C C------------------------------------------ C FMF042 CALCULATES Q FROM (T,P,Z). C------------------------------------------ SUBROUTINE FMF042(J,T,P,Z,Q) DOUBLE PRECISION T,P,Z,Q,X,Y,FMF024 INTEGER KPA,MESS,KSTAN,KAS,J COMMON/UNIT/KPA,MESS,KSTAN,KAS C CALL FMF043(J,T,P,X,Y) IF(J.EQ.0) THEN IF(KAS.EQ.1) THEN X=FMF024(X) Y=FMF024(Y) Z=FMF024(Z) ENDIF Q=(Z-X)/(Y-X) ELSE Q=-1.0D10 ENDIF RETURN END C C****************************************************** C C----------------------------------------------------------------- C FMF043-2 CALCULATES (X,Y) AT VAPOR-LIQUID EQUILIBLIUM STATE. C----------------------------------------------------------------- SUBROUTINE FMF043(J,T,P,X,Y) DOUBLE PRECISION T,P,X,Y c CALL FMF044(J,T,P,X,Y) IF(J.NE.0) THEN CALL FMF045(J,T,P,X,Y) ENDIF RETURN END C C******************************************************** C SUBROUTINE FMF044(J,T,P,X,Y) DOUBLE PRECISION T,P,X,Y,TC1,PC1,OMEGA1,TC2,PC2,OMEGA2, & A(2),B(2),AD,X1,X2,Y1,Y2,K1(2),K2(2),AMX,AMY,ADMX,ADMY, & BMX,BMY,ZL,ZV,ZLM,ZVM,FMF035,FL1,FL2,FV1,FV2,EPS INTEGER J,COMBI DOUBLE PRECISION PR(1:11),CP(1:20) COMMON/FMFC/PR,CP,COMBI C TC1=PR(2) PC1=PR(3) OMEGA1=PR(5) TC2=PR(7) PC2=PR(8) OMEGA2=PR(10) CALL FMF006(TC1,PC1,OMEGA1,T,P,A(1),AD,B(1)) CALL FMF006(TC2,PC2,OMEGA2,T,P,A(2),AD,B(2)) C X1=1.0D-7 Y1=1.0-X1 X2=1.0-X1 Y2=1.0-Y1 K1(1)=Y1/X1 K2(1)=Y2/X2 C ---------------------------------------------------------- DO 1000 N=1,200 CALL FMF011(T,P,X1,AMX,ADMX,BMX) CALL FMF011(T,P,Y1,AMY,ADMY,BMY) C CALL FMF073(AMX,BMX,ZV,ZLM) CALL FMF073(AMY,BMY,ZVM,ZL) C FL1=FMF035(1,1,P,A,B,AMX,BMX,X1,ZLM) FL2=FMF035(2,1,P,A,B,AMX,BMX,X2,ZLM) FV1=FMF035(1,1,P,A,B,AMY,BMY,Y1,ZVM) FV2=FMF035(2,1,P,A,B,AMY,BMY,Y2,ZVM) C ----------------------------------------------- EPS=DABS(FL1/FV1-1.0)+DABS(FL2/FV2-1.0) IF(EPS.LE.1.0D-7) THEN J=0 X=X1 Y=Y1 RETURN ENDIF C ----------------------------------------------- K1(2)=(FL1/FV1)*K1(1) K2(2)=(FL2/FV2)*K2(1) X1=(K2(2)-1.0)/(K2(2)-K1(2)) Y1=K1(2)*X1 X2=1.0-X1 Y2=1.0-Y1 K1(1)=K1(2) K2(1)=K2(2) 1000 CONTINUE C ---------------------------------------------------------- IF(EPS.LE.1.0D-5) THEN J=10 X=X1 Y=Y1 ELSEIF(EPS.LE.1.0D-3) THEN J=100 X=X1 Y=Y1 ELSE J=-1 X=1.0D-10 Y=1.0D-10 ENDIF RETURN END C C************************************************************ C SUBROUTINE FMF045(J,T,P,X,Y) DOUBLE PRECISION T,P,X,Y,Y1,Y2,PD,PC,VC DOUBLE PRECISION PR(1:11),CP(1:20) INTEGER COMBI COMMON/FMFC/PR,CP,COMBI c TC1=PR(2) Y1=0.0 IF(T.LT.TC1) THEN Y2=1.0 ELSE CALL FMF074(2,JC,Y2,T,PC,VC) ENDIF C DO 1000 N=1,200 Y=0.5*(Y1+Y2) CALL FMF032(J,T,PD,X,Y) IF(J.NE.0) THEN J=-1 X=-1.0D10 Y=-1.0D10 RETURN ENDIF C ------------------------------------------- IF(DABS(PD/P-1.0).LT.1.0D-5) THEN J=0 RETURN ENDIF C ------------------------------------------- IF(PD.LT.P) Y1=Y IF(PD.GT.P) Y2=Y 1000 CONTINUE J=-1 X=-1.0D10 Y=-1.0D10 RETURN END C C**************************************************************** C C-------------------------------------------------------------------- C FMF046-4 CALCULATES PROPERTIES AT VAPOR-LIQUID EQUILIBLIUM STATE. C-------------------------------------------------------------------- SUBROUTINE FMF046(J,T,P,X,Y,VL,VV,HL,HV,SL,SV) DOUBLE PRECISION T,P,X,Y,FMF072,VL,VV,HL,HV,SL,SV,ER C CALL FMF043(J,T,P,X,Y) IF(J.EQ.0) THEN VL=FMF072('V',1,T,P,X) VV=FMF072('V',3,T,P,Y) HL=FMF072('H',1,T,P,X) HV=FMF072('H',3,T,P,Y) SL=FMF072('S',1,T,P,X) SV=FMF072('S',3,T,P,Y) ELSE ER=-1.0D10 VL=ER VV=ER HL=ER HV=ER SL=ER SV=ER ENDIF RETURN END C C********************************************************* C SUBROUTINE FMF047(J,T,P,Z,V,H,S) DOUBLE PRECISION T,P,Z,V,H,S,X,Y,VL,VV,HL,HV,SL,SV,Q INTEGER J c CALL FMF046(J,T,P,X,Y,VL,VV,HL,HV,SL,SV) IF(J.EQ.0) THEN Q=(Z-X)/(Y-X) V=(1.0-Q)*VL+Q*VV H=(1.0-Q)*HL+Q*HV S=(1.0-Q)*SL+Q*SV ELSE ER=-1.0D10 V=ER H=ER S=ER ENDIF RETURN END C C************************************************************** C C----------------------------------------------- C FMF048 CALCULATES (V,H,S) FROM (T,P,Z). C----------------------------------------------- SUBROUTINE FMF048(J,T,P,Z,V,H,S) DOUBLE PRECISION T,P,Z,V,H,S,LIM INTEGER J,IP,FMF002 C LIM=1.0D-7 IF(Z.LT.LIM) THEN CALL FMF063(2,J,T,P,V,H,S) ELSEIF(Z.GT.1.0-LIM) THEN CALL FMF063(1,J,T,P,V,H,S) ELSE IP=FMF002(T,P,Z) IF(IP.EQ.2) THEN CALL FMF047(J,T,P,Z,V,H,S) ELSE CALL FMF068(IP,T,P,Z,V,H,S) J=0 ENDIF ENDIF RETURN END C C****************************************************** C C--------------------------------------------- C FMF049 CALCULATES (T,V,S) FROM (P,Z,H). C--------------------------------------------- SUBROUTINE FMF049(J,T,P0,Z0,V,H0,S) DOUBLE PRECISION T,P0,Z0,V,H,H0,S,TC,PC,VC, & T1,T2,H1,H2,DH,DH1,DH2,ER,HMAX,HMIN,FMF072 INTEGER IT,IT1,IT2,J,FMF003,FMF002,COMBI DOUBLE PRECISION PR(1:11),CP(1:20) COMMON/FMFC/PR,CP,COMBI C J=-1 T1=DMAX1(CP(7),CP(17)) T2=DMIN1(CP(8),CP(18)) IT1=FMF002(T1,P0,Z0) IT2=FMF002(T2,P0,Z0) HMIN=FMF072('H',IT1,T1,P0,Z0) HMAX=FMF072('H',IT2,T2,P0,Z0) IF(H0.GT.HMAX.OR.H0.LT.HMIN) THEN J=-2 ER=-1.0D20 T=ER V=ER S=ER RETURN ENDIF C CALL FMF074(1,J,Z0,TC,PC,VC) IF(P0.GT.PC) THEN IT=3 ELSE DO 1000 N=1,10 T=0.5*(T1+T2) CALL FMF048(J,T,P0,Z0,V,H,S) IF(H.GT.H0) THEN T2=T ELSE T1=T ENDIF IF(T1.GT.TC) THEN IT=3 ELSE IT1=FMF003(T1,P0,Z0) IT2=FMF003(T2,P0,Z0) IF(IT1.EQ.3) THEN IT=3 ELSEIF(IT2.EQ.1) THEN IT=1 ELSEIF(IT1.EQ.2.AND.IT2.EQ.2) THEN IT=2 ELSE IT=0 ENDIF ENDIF IF(IT.NE.0) GOTO 1200 1000 CONTINUE ENDIF C---------------------------------------------------------- 1200 IF(IT.EQ.0) THEN C*******METHOD OF BISECTION USING FMF048. ******** C DO 2000 N=1,100 T=0.5*(T1+T2) CALL FMF048(J,T,P0,Z0,V,H,S) CALL FMF075(J,H,H0,1.0D-6) IF(J.EQ.0) RETURN IF(H.GT.H0) THEN T2=T ELSE T1=T ENDIF 2000 CONTINUE C---------------------------------------------------------- ELSEIF(IT.EQ.2) THEN C*******METHOD OF REGULA-FALSI USING FMF047. ******** C CALL FMF047(J,T1,P0,Z0,V,H1,S) CALL FMF047(J,T2,P0,Z0,V,H2,S) DH1=H1-H0 DH2=H2-H0 DO 2100 N=1,50 T=(T1*DH2-T2*DH1)/(DH2-DH1) CALL FMF047(J,T,P0,Z0,V,H,S) CALL FMF050(J,H,H0) IF(J.EQ.0) RETURN DH=H-H0 IF(DH.GT.0.0) THEN T2=T DH2=DH ELSE T1=T DH1=DH ENDIF 2100 CONTINUE C---------------------------------------------------------- ELSE C*******METHOD OF REGULA-FALSI USING FMF068. ******** C CALL FMF068(IT,T1,P0,Z0,V,H1,S) CALL FMF068(IT,T2,P0,Z0,V,H2,S) DH1=H1-H0 DH2=H2-H0 DO 2200 N=1,50 T=(T1*DH2-T2*DH1)/(DH2-DH1) CALL FMF068(IT,T,P0,Z0,V,H,S) CALL FMF050(J,H,H0) IF(J.EQ.0) RETURN DH=H-H0 IF(DH.GT.0.0) THEN T1=T DH1=DH ELSE T2=T DH2=DH ENDIF 2200 CONTINUE ENDIF C---------------------------------------------------------- ER=-1.0D10 T=ER V=ER S=ER J=-1 RETURN END C C****************************************************** C SUBROUTINE FMF050(J,H,H0) DOUBLE PRECISION H,H0 INTEGER J c J=-1 IF(DABS(H0).LT.1.0) THEN IF(DABS(H-H0).LT.1.0) J=0 ELSE IF(DABS(H/H0-1.0).LT.1.0D-7) J=0 ENDIF RETURN END C C***************************************************** C C--------------------------------------------- C FMF051 CALCULATES (T,V,H) FROM (P,Z,S). C--------------------------------------------- SUBROUTINE FMF051(J,T,P0,Z0,V,H,S0) DOUBLE PRECISION T,P0,Z0,V,H,S,S0,TC,PC,VC, & T1,T2,S1,S2,DS,DS1,DS2,ER,SMAX,SMIN,FMF072 INTEGER IT,IT1,IT2,J,FMF003,FMF002,COMBI DOUBLE PRECISION PR(1:11),CP(1:20) COMMON/FMFC/PR,CP,COMBI C J=-1 T1=DMAX1(CP(7),CP(17)) T2=DMIN1(CP(8),CP(18)) IT1=FMF002(T1,P0,Z0) IT2=FMF002(T2,P0,Z0) SMIN=FMF072('S',IT1,T1,P0,Z0) SMAX=FMF072('S',IT2,T2,P0,Z0) IF(S0.GT.SMAX.OR.S0.LT.SMIN) THEN J=-2 ER=-1.0D20 T=ER V=ER H=ER RETURN ENDIF C CALL FMF074(1,J,Z0,TC,PC,VC) IF(P0.GT.PC) THEN IT=3 GOTO 1200 ENDIF C DO 1000 N=1,10 T=0.5*(T1+T2) CALL FMF048(J,T,P0,Z0,V,H,S) IF(S.GT.S0) THEN T2=T ELSE T1=T ENDIF IF(T1.GT.TC) THEN IT=3 ELSE IT1=FMF003(T1,P0,Z0) IT2=FMF003(T2,P0,Z0) IF(IT1.EQ.3) THEN IT=3 ELSEIF(IT2.EQ.1) THEN IT=1 ELSEIF(IT1.EQ.2.AND.IT2.EQ.2) THEN IT=2 ELSE IT=0 ENDIF ENDIF IF(IT.NE.0) GOTO 1200 1000 CONTINUE C---------------------------------------------------------- 1200 IF(IT.EQ.0) THEN C*******METHOD OF BISECTION USING FMF048. ******** C DO 2000 N=1,100 T=0.5*(T1+T2) CALL FMF048(J,T,P0,Z0,V,H,S) CALL FMF052(J,S,S0) IF(J.EQ.0) RETURN IF(S.GT.S0) THEN T2=T ELSE T1=T ENDIF 2000 CONTINUE C---------------------------------------------------------- ELSEIF(IT.EQ.2) THEN C*******METHOD OF REGULA-FALSI USING FMF047. ******** C CALL FMF047(J,T1,P0,Z0,V,H,S1) CALL FMF047(J,T2,P0,Z0,V,H,S2) DS1=S1-S0 DS2=S2-S0 DO 2100 N=1,50 T=(T1*DS2-T2*DS1)/(DS2-DS1) CALL FMF047(J,T,P0,Z0,V,H,S) CALL FMF052(J,S,S0) IF(J.EQ.0) RETURN DS=S-S0 IF(DS.GT.0.0) THEN T2=T DS2=DS ELSE T1=T DS1=DS ENDIF 2100 CONTINUE C---------------------------------------------------------- ELSE C*******METHOD OF REGULA-FALSI USING FMF068. ******** C CALL FMF068(IT,T1,P0,Z0,V,H,S1) CALL FMF068(IT,T2,P0,Z0,V,H,S2) DS1=S1-S0 DS2=S2-S0 DO 2200 N=1,50 T=(T1*DS2-T2*DS1)/(DS2-DS1) CALL FMF068(IT,T,P0,Z0,V,H,S) CALL FMF052(J,S,S0) IF(J.EQ.0) RETURN DS=S-S0 IF(DS.GT.0.0) THEN T1=T DS1=DS ELSE T2=T DS2=DS ENDIF 2200 CONTINUE ENDIF C---------------------------------------------------------- ER=-1.0D10 T=ER V=ER H=ER RETURN END C C********************************************************** C SUBROUTINE FMF052(J,S,S0) DOUBLE PRECISION S,S0 INTEGER J c J=-1 IF(DABS(S0).LT.1.0) THEN IF(DABS(S-S0).LT.1.0) J=0 ELSE IF(DABS(S/S0-1.0).LT.1.0D-7) J=0 ENDIF RETURN END C C*********************************************************** C C--------------------------------------------- C FMF053 CALCULATES (T,H,S) FROM (P,Z,V). C--------------------------------------------- SUBROUTINE FMF053(J,T,P0,Z0,V0,H,S) DOUBLE PRECISION T,P0,Z0,V,V0,H,S,TC,PC,VC, & T1,T2,V1,V2,DV,DV1,DV2,ER,TC1,TC2,PC1,PC2,R,B1,B2,VMIN INTEGER IT,IT1,IT2,J,FMF003,COMBI DOUBLE PRECISION PR(1:11),CP(1:20) COMMON/FMFC/PR,CP,COMBI C TC1=PR(2) PC1=PR(3) TC2=PR(7) PC2=PR(8) R=8314.5 B1=0.077796073903888*R*TC1/PC1 B2=0.077796073903888*R*TC2/PC2 VMIN=B1*Z0+B2*(1.0-Z0) IF(V0.LE.VMIN) THEN J=-2 ER=-1.0D20 T=ER H=ER S=ER RETURN ENDIF C T1=DMAX1(CP(7),CP(17)) T2=DMIN1(CP(8),CP(18)) CALL FMF074(1,J,Z0,TC,PC,VC) IF(P0.GT.PC) THEN IT=3 GOTO 1200 ENDIF C DO 1000 N=1,10 T=0.5*(T1+T2) CALL FMF048(J,T,P0,Z0,V,H,S) IF(V.GT.V0) THEN T2=T ELSE T1=T ENDIF IF(T1.GT.TC) THEN IT=3 ELSE IT1=FMF003(T1,P0,Z0) IT2=FMF003(T2,P0,Z0) IF(IT1.EQ.3) THEN IT=3 ELSEIF(IT2.EQ.1) THEN IT=1 ELSEIF(IT1.EQ.2.AND.IT2.EQ.2) THEN IT=2 ELSE IT=0 ENDIF ENDIF IF(IT.NE.0) GOTO 1200 1000 CONTINUE C---------------------------------------------------------- 1200 IF(IT.EQ.0) THEN C*******METHOD OF BISECTION USING FMF048. ******** C DO 2000 N=1,100 T=0.5*(T1+T2) CALL FMF048(J,T,P0,Z0,V,H,S) CALL FMF054(J,V,V0) IF(J.EQ.0) RETURN IF(V.GT.V0) THEN T2=T ELSE T1=T ENDIF 2000 CONTINUE C---------------------------------------------------------- ELSEIF(IT.EQ.2) THEN C*******METHOD OF REGULA-FALSI USING FMF047. ******** C CALL FMF047(J,T1,P0,Z0,V1,H,S) CALL FMF047(J,T2,P0,Z0,V2,H,S) DV1=V1-V0 DV2=V2-V0 DO 2100 N=1,50 T=(T1*DV2-T2*DV1)/(DV2-DV1) CALL FMF047(J,T,P0,Z0,V,H,S) CALL FMF054(J,V,V0) IF(J.EQ.0) RETURN DV=V-V0 IF(DV.GT.0.0) THEN T2=T DV2=DV ELSE T1=T DV1=DV ENDIF 2100 CONTINUE C---------------------------------------------------------- ELSE C*******METHOD OF REGULA-FALSI USING FMF068. ******** C CALL FMF068(IT,T1,P0,Z0,V1,H,S) CALL FMF068(IT,T2,P0,Z0,V2,H,S) DV1=V1-V0 DV2=V2-V0 DO 2200 N=1,50 T=(T1*DV2-T2*DV1)/(DV2-DV1) CALL FMF068(IT,T,P0,Z0,V,H,S) CALL FMF054(J,V,V0) IF(J.EQ.0) RETURN DV=V-V0 T1=T2 T2=T DV1=DV2 DV2=DV 2200 CONTINUE ENDIF C---------------------------------------------------------- J=-1 ER=-1.0D10 T=ER H=ER S=ER RETURN END C C********************************************************** C SUBROUTINE FMF054(J,V,V0) DOUBLE PRECISION V,V0,LV,LV0 INTEGER J c J=-1 LV=DLOG(V) LV0=DLOG(V0) IF(DABS(LV0).LT.1.0D-5) THEN IF(DABS(LV).LT.1.0D-7) J=0 ELSE IF(DABS(LV/LV0-1.0).LT.1.0D-7) J=0 ENDIF RETURN END C C********************************************************* C FUNCTION FMF055(I,XOUT) DOUBLE PRECISION XOUT,FMF024 INTEGER KPA,MESS,KSTAN,KAS,I COMMON/UNIT/KPA,MESS,KSTAN,KAS C IF(I.EQ.1) THEN IF(KPA.EQ.1.OR.KPA.EQ.3) THEN FMF055=REAL(XOUT)-273.15 ELSE FMF055=REAL(XOUT) ENDIF RETURN ENDIF C ------------------------------------------------ IF(I.EQ.2) THEN IF(KPA.EQ.1.OR.KPA.EQ.2) THEN FMF055=REAL(XOUT)*1.0E-5 ELSE FMF055=REAL(XOUT) ENDIF RETURN ENDIF C ------------------------------------------------ IF(I.EQ.3) THEN IF(KAS.EQ.1) THEN FMF055=REAL(FMF024(XOUT)) ELSE FMF055=REAL(XOUT) ENDIF RETURN ENDIF END C C************************************************************ C FUNCTION FMF056(I,XOUT,Z) DOUBLE PRECISION MW,HSTAN,SSTAN,XOUT,Z,FMF036 $ ,XOUTT INTEGER KPA,MESS,KSTAN,KAS,I,COMBI DOUBLE PRECISION PR(1:11),CP(1:20) COMMON/UNIT/KPA,MESS,KSTAN,KAS COMMON/FMFC/PR,CP,COMBI C IF(I.EQ.1) THEN IF(KAS.EQ.1) THEN MW=FMF036(REAL(Z)) XOUTT=XOUT/MW ELSE XOUTT=XOUT ENDIF FMF056=REAL(XOUTT) RETURN ENDIF C --------------------------------------------------- IF(I.EQ.2) THEN IF(KSTAN.EQ.1) THEN HSTAN=(Z*CP(9)*PR(1)+(1.0-Z)*CP(19)*PR(6))*1.0D3 XOUT=XOUT+HSTAN ENDIF IF(KAS.EQ.1) THEN MW=FMF036(REAL(Z)) XOUT=XOUT/MW ENDIF FMF056=REAL(XOUT) RETURN ENDIF C --------------------------------------------------- IF(I.EQ.3) THEN IF(KSTAN.EQ.1) THEN SSTAN=(Z*CP(10)*PR(1)+(1.0-Z)*CP(20)*PR(6))*1.0D3 XOUT=XOUT+SSTAN ENDIF IF(KAS.EQ.1) THEN MW=FMF036(REAL(Z)) XOUT=XOUT/MW ENDIF FMF056=REAL(XOUT) RETURN ENDIF END C C********************************************************** C C------------------------------------------------------ C FMF057 CALCULATES FUGACITY OF PURE SUBSTANCE. C------------------------------------------------------ DOUBLE PRECISION FUNCTION FMF057(P,Z,A,B) DOUBLE PRECISION P,Z,A,B,LNF,F c LNF=(Z-1.0)-LOG(Z-B)-(A/(2.0*SQRT(2.0)*B))*LOG((Z+(1.0 & +SQRT(2.0))*B)/(Z+(1.0-SQRT(2.0))*B)) F=P*EXP(LNF) FMF057=F RETURN END C C******************************************************** C SUBROUTINE FMF058(T,A,B,VMIN,VMAX,PMIN,PMAX) DOUBLE PRECISION T,A,B,VMIN,VMAX,PMIN,PMAX $ ,C(5),XR(4),XI(4),FMF018 c R=8.31451 C(1)=-R*T C(2)=2.0*A-4*B*R*T C(3)=-2.0*A*B-2.0*B**2*R*T C(4)=4.0*B**3*R*T-2.0*A*B**2 C(5)=2.0*A*B**3-B**4*R*T CALL FMF005(C,XR,XI) VMIN=DMIN1(XR(3),XR(4)) VMAX=DMAX1(XR(3),XR(4)) PMIN=FMF018(T,VMIN,A,B) PMAX=FMF018(T,VMAX,A,B) RETURN END C C************************************************************** C C---------------------------------------------------------------- C FMF059 CALCULATES SATURATED VAPOR PRESSURE OF PURE SUBSTANCE. C---------------------------------------------------------------- SUBROUTINE FMF059(I,J,T,P) DOUBLE PRECISION TC,PC,OMEGA,T,P,P1,P2,A,AD,B,ZV,ZL, & FV,FL,DF1,DF2,FMF057 INTEGER J,COMBI DOUBLE PRECISION PR(1:11),CP(1:20) COMMON/FMFC/PR,CP,COMBI C IF(I.EQ.1) THEN TC=PR(2) PC=PR(3) OMEGA=PR(5) ELSEIF(I.EQ.2) THEN TC=PR(7) PC=PR(8) OMEGA=PR(10) ENDIF IF(T.GT.TC) THEN J=-2 P=-1.0D20 RETURN ENDIF C CALL FMF060(TC,PC,OMEGA,T,P1) DO 1000 N=1,100 CALL FMF006(TC,PC,OMEGA,T,P1,A,AD,B) CALL FMF073(A,B,ZV,ZL) FV=FMF057(P1,ZV,A,B) FL=FMF057(P1,ZL,A,B) DF1=FL-FV IF(N.EQ.1) P2=P1*FL/FV CALL FMF006(TC,PC,OMEGA,T,P2,A,AD,B) CALL FMF073(A,B,ZV,ZL) FV=FMF057(P2,ZV,A,B) FL=FMF057(P2,ZL,A,B) DF2=FL-FV C -------------------------------------------- IF(ABS(FL/FV-1.0).LT.1.0D-7) THEN J=0 P=P2 RETURN ENDIF C -------------------------------------------- P=P1-DF1/(DF2-DF1)*(P2-P1) P1=P2 P2=P 1000 CONTINUE J=-1 P=-1.0D10 RETURN END C C********************************************************* C C------------------------------------------------------------------ C FMF060 CALCULATES INITIAL VALUE OF SATURATED VAPOR PRESSURE. C------------------------------------------------------------------ SUBROUTINE FMF060(TC,PC,OMEGA,T,P) DOUBLE PRECISION TC,PC,OMEGA,P,T,A,AD,B,C(5),XR(4) $ ,XI(4),V1,V2,V,FMF018 c CALL FMF007(TC,PC,OMEGA,T,A,AD,B) R=8.31451 C(1)=-R*T C(2)=2.0*A-4*B*R*T C(3)=-2.0*A*B-2.0*B**2*R*T C(4)=4.0*B**3*R*T-2.0*A*B**2 C(5)=2.0*A*B**3-B**4*R*T CALL FMF005(C,XR,XI) V1=XR(3) V2=XR(4) V=0.5*V1+0.5*V2 P=FMF018(T,V,A,B) RETURN END C C************************************************************** C C------------------------------------------------------------------- C FMF061 CALCULATES SATURATED VAPOR TEMPERATURE OF PURE SUBSTANCE. C------------------------------------------------------------------- SUBROUTINE FMF061(I,J,T,P) DOUBLE PRECISION P,T,TC,PC,T0,P0,T1,P1,PS,OMEGA,TB,PB INTEGER I,J,JJ,N,COMBI DOUBLE PRECISION PR(1:11),CP(1:20) COMMON/FMFC/PR,CP,COMBI C IF(I.EQ.1) THEN TC=PR(2) PC=PR(3) OMEGA=PR(5) ELSEIF(I.EQ.2) THEN TC=PR(7) PC=PR(8) OMEGA=PR(10) ENDIF C TB=0.7*TC PB=PC*10.0**(-OMEGA-1.0) C ------------------------------------ IF(P.GT.PC) THEN J=-2 T=-1.0D20 RETURN ENDIF C ------------------------------------ T0=TB P0=PB T1=TC P1=PC DO 1000 N=1,100 T=T0+(T1-T0)/(P1-P0)*(P-P0) CALL FMF059(I,JJ,T,PS) C --------------------------------------- IF(DABS(PS/P-1.0).LT.1.0D-5) THEN J=0 RETURN ENDIF C --------------------------------------- T0=T1 P0=P1 T1=T P1=PS 1000 CONTINUE J=-1 T=-1.0D10 RETURN END C C************************************************************ C C------------------------------------------------------ C FMF062 CALCULATES PROPERTIES AT SATURATED STATE. C------------------------------------------------------ SUBROUTINE FMF062(I,J,T,P,VL,VV,HL,HV,SL,SV) DOUBLE PRECISION TC,Y,T,P,VL,VV,HL,HV,SL,SV,FMF072,ER INTEGER I,J,COMBI DOUBLE PRECISION PR(1:11),CP(1:20) COMMON/FMFC/PR,CP,COMBI C IF(I.EQ.1) THEN TC=PR(2) Y=1.0 ELSEIF(I.EQ.2) THEN TC=PR(7) Y=0.0 ENDIF C IF(T.GT.TC) THEN J=-2 ELSE CALL FMF059(I,J,T,P) ENDIF C IF(J.EQ.0) THEN VL=FMF072('V',1,T,P,Y) VV=FMF072('V',3,T,P,Y) HL=FMF072('H',1,T,P,Y) HV=FMF072('H',3,T,P,Y) SL=FMF072('S',1,T,P,Y) SV=FMF072('S',3,T,P,Y) ELSEIF(J.EQ.-1) THEN ER=-1.0D10 P=ER VL=ER VV=ER HL=ER HV=ER SL=ER SV=ER ELSEIF(J.EQ.-2) THEN ER=-1.0D20 P=ER VL=ER VV=ER HL=ER HV=ER SL=ER SV=ER ENDIF RETURN END C C********************************************************** C C----------------------------------------------------- C FMF063 CALCULATES PROPERTIES OF PURE SUBSTANCES. C----------------------------------------------------- SUBROUTINE FMF063(I,J,T,P,V,H,S) DOUBLE PRECISION TC,Y,T,P,PS,V,H,S,FMF072,ER INTEGER I,IP,COMBI DOUBLE PRECISION PR(1:11),CP(1:20) COMMON/FMFC/PR,CP,COMBI C IF(I.EQ.1) THEN TC=PR(2) Y=1.0 ELSEIF(I.EQ.2) THEN TC=PR(7) Y=0.0 ENDIF C IF(T.LT.TC) THEN CALL FMF059(I,J,T,PS) IF(P.GE.PS) THEN IP=1 ELSE IP=3 ENDIF ELSE IP=3 ENDIF C IF(J.EQ.0) THEN V=FMF072('V',IP,T,P,Y) H=FMF072('H',IP,T,P,Y) S=FMF072('S',IP,T,P,Y) ELSE ER=-1.0D10 V=ER H=ER S=ER ENDIF RETURN END C C******************************************************** C SUBROUTINE FMF064(I,II,J,T,P,V,H,S) DOUBLE PRECISION T,P,V,H,S INTEGER I,II c IF(I.EQ.1) THEN CALL FMF063(II,J,T,P,V,H,S) ELSEIF(I.EQ.2) THEN CALL FMF065(II,J,T,P,V,H,S) ELSEIF(I.EQ.3) THEN CALL FMF066(II,J,T,P,V,H,S) ELSEIF(I.EQ.4) THEN CALL FMF067(II,J,T,P,V,H,S) ENDIF RETURN END C C******************************************************* C SUBROUTINE FMF065(I,J,T,P0,V,H0,S) DOUBLE PRECISION T,P0,V,H0,S,PC,TMAX $ ,TMIN,ER,TS,PS, & VL,VV,HL,HV,SL,SV,T1,T2,H,H1,H2,DH,DH1,DH2,Q INTEGER I,J,COMBI DOUBLE PRECISION PR(1:11),CP(1:20) COMMON/FMFC/PR,CP,COMBI C IF(I.EQ.1) THEN PC=PR(3) TMIN=CP(7) TMAX=CP(8) ELSE PC=PR(8) TMIN=CP(17) TMAX=CP(18) ENDIF C ---------------------------------------------------- IF(P0.LT.PC) THEN CALL FMF061(I,J,TS,P0) CALL FMF062(I,J,TS,PS,VL,VV,HL,HV,SL,SV) IF(H0.GE.HL.AND.H0.LE.HV) THEN J=0 T=TS Q=(H0-HL)/(HV-HL) V=(1.0-Q)*VL+Q*VV S=(1.0-Q)*SL+Q*SV RETURN ELSEIF(H0.LT.HL) THEN T1=TMIN T2=TS CALL FMF063(I,J,T1,P0,V,H1,S) DH1=H1-H0 DH2=HL-H0 ELSEIF(H0.GT.HV) THEN T1=TS T2=TMAX CALL FMF063(I,J,T2,P0,V,H2,S) DH1=HV-H0 DH2=H2-H0 ENDIF ENDIF C ---------------------------------------------------- IF(P0.GT.PC) THEN T1=TMIN T2=TMAX CALL FMF063(I,J,T1,P0,V,H1,S) CALL FMF063(I,J,T2,P0,V,H2,S) DH1=H1-H0 DH2=H2-H0 ENDIF C ---------------------------------------------------- IF(DH2.LT.0.0) THEN J=-2 ER=-1.0D20 T=ER V=ER S=ER RETURN ENDIF C----------------------------------------------------------- C**** METHOD OF REGULA-FALSI AND BISECTION USING FMF063. **** DO 1000 N=1,50 IF(N.LE.3) T=0.5*(T1+T2) IF(N.GE.4) T=(T1*DH2-T2*DH1)/(DH2-DH1) CALL FMF063(I,J,T,P0,V,H,S) C --------------------------------------- IF(DABS(H0).LT.1.0D-7) THEN IF(DABS(H).LT.1.0D-7) RETURN ELSE IF(DABS(H/H0-1.0).LT.1.0D-7) RETURN ENDIF C --------------------------------------- DH=H-H0 IF(DH.GT.0.0) THEN T2=T DH2=DH ELSE T1=T DH1=DH ENDIF 1000 CONTINUE J=-1 ER=-1.0D10 T=ER V=ER S=ER RETURN END C C************************************************************* C SUBROUTINE FMF066(I,J,T,P0,V,H,S0) DOUBLE PRECISION T,P0,V,H,S,S0,PC $ ,TMAX,TMIN,ER,TS,PS, & VL,VV,HL,HV,SL,SV,T1,T2,S1,S2,DS,DS1,DS2,Q INTEGER I,J,COMBI DOUBLE PRECISION PR(1:11),CP(1:20) COMMON/FMFC/PR,CP,COMBI C IF(I.EQ.1) THEN PC=PR(3) TMIN=CP(7) TMAX=CP(8) ELSE PC=PR(8) TMIN=CP(17) TMAX=CP(18) ENDIF C ---------------------------------------------------- IF(P0.LT.PC) THEN CALL FMF061(I,J,TS,P0) CALL FMF062(I,J,TS,PS,VL,VV,HL,HV,SL,SV) IF(S0.GE.SL.AND.S0.LE.SV) THEN J=0 T=TS Q=(S0-SL)/(SV-SL) V=(1.0-Q)*VL+Q*VV H=(1.0-Q)*HL+Q*HV RETURN ELSEIF(S0.LT.SL) THEN T1=TMIN T2=TS CALL FMF063(I,J,T1,P0,V,H,S1) DS1=S1-S0 DS2=SL-S0 ELSEIF(S0.GT.SV) THEN T1=TS T2=TMAX CALL FMF063(I,J,T2,P0,V,H,S2) DS1=SV-S0 DS2=S2-S0 ENDIF ENDIF C ---------------------------------------------------- IF(P0.GT.PC) THEN T1=TMIN T2=TMAX CALL FMF063(I,J,T1,P0,V,H,S1) CALL FMF063(I,J,T2,P0,V,H,S2) DS1=S1-S0 DS2=S2-S0 ENDIF C ---------------------------------------------------- IF(DS2.LT.0.0) THEN J=-2 ER=-1.0D20 T=ER V=ER H=ER RETURN ENDIF C----------------------------------------------------------- C**** METHOD OF REGULA-FALSI AND BISECTION USING FMF063. **** DO 1000 N=1,50 IF(N.LE.3) T=0.5*(T1+T2) IF(N.GE.4) T=(T1*DS2-T2*DS1)/(DS2-DS1) CALL FMF063(I,J,T,P0,V,H,S) C --------------------------------------- IF(DABS(S0).LT.1.0D-7) THEN IF(DABS(S).LT.1.0D-7) RETURN ELSE IF(DABS(S/S0-1.0).LT.1.0D-7) RETURN ENDIF C --------------------------------------- DS=S-S0 IF(DS.GT.0.0) THEN T2=T DS2=DS ELSE T1=T DS1=DS ENDIF 1000 CONTINUE J=-1 ER=-1.0D10 T=ER V=ER H=ER RETURN END C C************************************************************* C SUBROUTINE FMF067(I,J,T,P0,V0,H,S) DOUBLE PRECISION T,P0,V,V0,H,S,TC,PC,VC,ER,TS,PS,VMIN, & VL,VV,HL,HV,SL,SV,T1,T2,V2,DV,DV1,DV2,Q,DT INTEGER I,J,COMBI DOUBLE PRECISION PR(1:11),CP(1:20) COMMON/FMFC/PR,CP,COMBI C IF(I.EQ.1) THEN TC=PR(2) PC=PR(3) VC=PR(4) ELSE TC=PR(7) PC=PR(8) VC=PR(9) ENDIF C R=8.31451 VMIN=0.077796073903888*R*TC/PC*1.0D3 IF(V0.LE.VMIN) THEN J=-2 ER=-1.0D20 T=ER H=ER S=ER RETURN ENDIF C ---------------------------------------------------- IF(P0.LT.PC) THEN CALL FMF061(I,J,TS,P0) CALL FMF062(I,J,TS,PS,VL,VV,HL,HV,SL,SV) C ---------------------------------------------------- IF(V0.GE.VL.AND.V0.LE.VV) THEN J=0 T=TS Q=(V0-VL)/(VV-VL) WRITE(*,*) 'Q=',Q H=(1.0-Q)*HL+Q*HV S=(1.0-Q)*SL+Q*SV RETURN C ---------------------------------------------------- ELSEIF(V0.GT.VV) THEN T1=TS T2=TS*1.01 CALL FMF063(I,J,T2,P0,V2,H,S) DV1=VV-V0 DV2=V2-V0 C ******* METHOD OF SECANT USING FMF063. ******* DO 1000 N=1,50 T=(T1*DV2-T2*DV1)/(DV2-DV1) CALL FMF063(I,J,T,P0,V,H,S) IF(DABS(V/V0-1.0).LT.1.0D-7) RETURN DV=V-V0 T1=T2 T2=T DV1=DV2 DV2=DV 1000 CONTINUE C ---------------------------------------------------- ELSEIF(V0.LT.VL) THEN T1=TS DV1=VL-V0 ENDIF ENDIF C ---------------------------------------------------- IF(P0.GT.PC) THEN T1=TC DV1=VC-V0 ENDIF C ---------------------------------------------------- DT=5.0 IF(V0.LT.VC) THEN DO 2000 N=1,100 T=T1-DT CALL FMF063(I,J,T,P0,V,H,S) IF(DABS(V/V0-1.0).LT.1.0D-7) RETURN IF(V.LT.V0) THEN DT=DT*0.5 ELSE T1=T ENDIF 2000 CONTINUE ENDIF C ---------------------------------------------------- IF(V0.GT.VC) THEN DO 3000 N=1,100 T=T1+DT CALL FMF063(I,J,T,P0,V,H,S) IF(DABS(V/V0-1.0).LT.1.0D-7) RETURN IF(V.GT.V0) THEN DT=DT*0.5 ELSE T1=T ENDIF 3000 CONTINUE ENDIF C ---------------------------------------------------- J=-1 ER=-1.0D10 T=ER H=ER S=ER RETURN END C C************************************************************* C C------------------------------------------------------- C FMF068 CALCULATES PROPERTIES AT SINGLE PHASE REGION C------------------------------------------------------- SUBROUTINE FMF068(I,T,P,Z,V,H,S) DOUBLE PRECISION T,P,Z,V,H,S,FMF072 INTEGER I c V=FMF072('V',I,T,P,Z) H=FMF072('H',I,T,P,Z) S=FMF072('S',I,T,P,Z) RETURN END C C***************************************************** C DOUBLE PRECISION FUNCTION FMF069(I,T,P,Y) DOUBLE PRECISION T,P,Y,R,ZV,ZL,Z,AM,BM,ADM INTEGER I C CALL FMF011(T,P,Y,AM,ADM,BM) CALL FMF073(AM,BM,ZV,ZL) IF(I.EQ.1) THEN Z=ZL ELSE Z=ZV ENDIF C R=8.31451 FMF069=(Z*R*T/P)*1.0D3 RETURN END C C---------------------------------------------------------------------- C DOUBLE PRECISION FUNCTION FMF070(I,T,P,Y) DOUBLE PRECISION T,P,Y,R,Z,ZV,ZL,C(2,5) $ ,A,B,AD,HIG(2),HIGM, & HR,T0,LIM,TMAX,TMIN,EQ(2) INTEGER I,K,COMBI DOUBLE PRECISION PR(1:11),CP(1:20) COMMON/FMFC/PR,CP,COMBI C EQ(1)=CP(1) C(1,1)=CP(2) C(1,2)=CP(3) C(1,3)=CP(4) C(1,4)=CP(5) C(1,5)=CP(6) EQ(2)=CP(11) C(2,1)=CP(12) C(2,2)=CP(13) C(2,3)=CP(14) C(2,4)=CP(15) C(2,5)=CP(16) C ------------------------------------------ LIM=1.0D-7 IF(Y.LT.LIM) THEN TMIN=CP(17) TMAX=CP(18) ELSEIF(Y.GT.1.0-LIM) THEN TMIN=CP(7) TMAX=CP(8) ELSE TMIN=DMAX1(CP(7),CP(17)) TMAX=DMIN1(CP(8),CP(18)) ENDIF IF(T.GT.TMAX.OR.T.LT.TMIN) THEN FMF070=-1.0D20 RETURN ENDIF C ------------------------------------------ CALL FMF011(T,P,Y,A,AD,B) CALL FMF073(A,B,ZV,ZL) IF(I.EQ.1) THEN Z=ZL ELSE Z=ZV ENDIF C R=8.31451 T0=298.15 C C-------ENTHALPY AT IDIAL GAS STATE ------------------------------- DO 1000 K=1,2 IF(EQ(K).LT.1.1) THEN HIG(K)=C(K,1)*T+C(K,2)*C(K,3)/DTANH(C(K,3)/T) & -C(K,4)*C(K,5)*DTANH(C(K,5)/T) & -C(K,1)*T0-C(K,2)*C(K,3)/DTANH(C(K,3)/T0) & +C(K,4)*C(K,5)*DTANH(C(K,5)/T0) ELSE HIG(K)=C(K,1)*(T-T0)+C(K,2)*(T**2-T0**2)/2.0+C(K,3) & *(T**3-T0**3)/3.0+C(K,4)*(T**4-T0**4)/4.0 ENDIF 1000 CONTINUE HIGM=Y*HIG(1)+(1.0-Y)*HIG(2) C C------RESIDUAL ENTHALPY ------------------------------------------ HR=R*T*(Z-1.0)+R*T*(T*AD-A)/(2.0*DSQRT(2.0D0)*B) & *DLOG((Z+(1.0+DSQRT(2.0D0))*B)/(Z+(1.0-DSQRT(2.0D0))*B)) C FMF070=HIGM+HR*1.0D3 C RETURN END C C******************************************************** C DOUBLE PRECISION FUNCTION FMF071(I,T,P,Y) DOUBLE PRECISION T,P,Y,R,RR,Z,ZV,ZL,C(2,5) $ ,SIG(2),SIGM,SR,SMIX & ,T0,P0,A,AD,B,LIM,TMAX,TMIN,EQ(2) INTEGER I,K,COMBI DOUBLE PRECISION PR(1:11),CP(1:20) COMMON/FMFC/PR,CP,COMBI C EQ(1)=CP(1) C(1,1)=CP(2) C(1,2)=CP(3) C(1,3)=CP(4) C(1,4)=CP(5) C(1,5)=CP(6) EQ(2)=CP(11) C(2,1)=CP(12) C(2,2)=CP(13) C(2,3)=CP(14) C(2,4)=CP(15) C(2,5)=CP(16) C ------------------------------------------ LIM=1.0D-7 IF(Y.LT.LIM) THEN TMIN=CP(17) TMAX=CP(18) ELSEIF(Y.GT.1.0-LIM) THEN TMIN=CP(7) TMAX=CP(8) ELSE TMIN=DMAX1(CP(7),CP(17)) TMAX=DMIN1(CP(8),CP(18)) ENDIF IF(T.GT.TMAX.OR.T.LT.TMIN) THEN FMF071=-1.0D20 RETURN ENDIF C ------------------------------------------ CALL FMF011(T,P,Y,A,AD,B) CALL FMF073(A,B,ZV,ZL) IF(I.EQ.1) THEN Z=ZL ELSE Z=ZV ENDIF C R=8.31451 RR=8314.51 T0=298.15 P0=1.0E5 C C-------ENTROPY AT IDIAL GAS STATE ------------------------------- DO 1000 K=1,2 IF(EQ(K).LT.1.1) THEN SIG(K)=C(K,1)*DLOG(T) & +C(K,2)*(C(K,3)/T/DTANH(C(K,3)/T) & -DLOG(DABS(DSINH(C(K,3)/T)))) & -C(K,4)*(C(K,5)/T*DTANH(C(K,5)/T) & -DLOG(DCOSH(C(K,5)/T))) & -C(K,1)*DLOG(T0) & -C(K,2)*(C(K,3)/T0/DTANH(C(K,3)/T0) & -DLOG(DABS(DSINH(C(K,3)/T0)))) & +C(K,4)*(C(K,5)/T0*DTANH(C(K,5)/T0) & -DLOG(DCOSH(C(K,5)/T0))) & -RR*DLOG(P/P0) ELSE SIG(K)=C(K,1)*(DLOG(T/T0))+C(K,2)*(T-T0) & +C(K,3)*(T**2-T0**2)/2.0+C(K,4) & *(T**3-T0**3)/3.0-RR*DLOG(P/P0) ENDIF 1000 CONTINUE SIGM=Y*SIG(1)+(1.0-Y)*SIG(2) C C------RESIDUAL ENTROPY ------------------------------------------ SR=R*DLOG(Z-B)+R*T*AD/(2.0*DSQRT(2.0D0)*B) & *DLOG((Z+(1.0+DSQRT(2.0D0))*B)/(Z+(1.0-DSQRT(2.0D0))*B)) C C------MIXING ENTROPY -------------------------------------------- IF(Y.LT.1.0D-10.OR.Y.GT.1.0-1.0D-10) THEN SMIX=0.0 ELSE SMIX=Y*(-R*DLOG(Y))+(1.0-Y)*(-R*DLOG(1.0-Y)) ENDIF C FMF071=SIGM+(SR+SMIX)*1.0D3 RETURN END C C******************************************************** C C------------------------------------------------- C FMF072 CALCULATES V OR H OR S FROM (T,P,Y). C------------------------------------------------- DOUBLE PRECISION FUNCTION FMF072(A,I,T,P,Y) DOUBLE PRECISION T,P,Y,FMF069,FMF070,FMF071 CHARACTER A c IF(A.EQ.'V') THEN FMF072=FMF069(I,T,P,Y) ELSEIF(A.EQ.'H') THEN FMF072=FMF070(I,T,P,Y) ELSEIF(A.EQ.'S') THEN FMF072=FMF071(I,T,P,Y) ENDIF RETURN END C C************************************************************ C C-------------------------------------------------------------- C FMF073 CALCULATES REAL ROOTS OF COMPRESSIBILITY FACTOR,Z. C-------------------------------------------------------------- SUBROUTINE FMF073(A,B,ZV,ZL) DOUBLE PRECISION A,B,C(3),X(3),ZV,ZL INTEGER IND c C(1)=B-1.0 C(2)=A-3.0*B**2-2.0*B C(3)=-A*B+B**2+B**3 CALL FMF001(C,X,IND) IF(IND.EQ.1) THEN ZV=X(1) ZL=X(1) ELSEIF(IND.EQ.3) THEN ZV=MAX(X(1),X(2),X(3)) ZL=MIN(X(1),X(2),X(3)) ENDIF RETURN END C C****************************************************** C SUBROUTINE FMF075(J,H,H0,EP) DOUBLE PRECISION H,H0,EP INTEGER J c J=-1 IF(DABS(H0).LT.1.0) THEN IF(DABS(H-H0).LT.EP) J=0 ELSE IF(DABS(H/H0-1.0).LT.EP) J=0 ENDIF RETURN END C C***************************************************** C c ---------------------------------------------------- c This program calculates the critical point T,P,V c from Z. c ---------------------------------------------------- SUBROUTINE FMF074(I,J,ZC,TC,PC,VC) DOUBLE PRECISION ZC,TC,PC,VC,TA(1:5),PA(1:5),VA(1:5),ZA(1:5) $ ,T01(1:5),P01(1:5),V01(1:5),Z01(1:5) $ ,PR(1:11),CP(1:20) INTEGER I,J,COMBI,N COMMON/FMFC/ PR,CP,COMBI c DATA T01/4.568311D2,-5.115951D1,-2.746161D1,1.653764D0 $ ,-1.049713D1/ DATA P01/3.663638D6,2.541498D6,-4.256569D5,9.467489D5 $ ,-1.751854D6/ DATA V01/3.186666D-1,-1.570070D-1,3.188625D-2,-4.257654D-2 $ ,3.859562D-2/ DATA Z01/-7.113347D1,7.465860D-1,-2.879637D-3,4.950808D-6 $ ,-3.236612D-9/ IF (COMBI.EQ.1) THEN DO 10 N=1,5 TA(N)=T01(N) PA(N)=P01(N) VA(N)=V01(N) ZA(N)=Z01(N) 10 CONTINUE ELSE CALL FMF012(I,J,ZC,TC,PC,VC) RETURN ENDIF IF (I.EQ.1) THEN TC=TA(1)+TA(2)*ZC+TA(3)*ZC**2+TA(4)*ZC**3+TA(5)*ZC**4 PC=PA(1)+PA(2)*ZC+PA(3)*ZC**2+PA(4)*ZC**3+PA(5)*ZC**4 VC=VA(1)+VA(2)*ZC+VA(3)*ZC**2+VA(4)*ZC**3+VA(5)*ZC**4 J=0 RETURN ELSEIF (I.EQ.2) THEN ZC=ZA(1)+ZA(2)*TC+ZA(3)*TC**2+ZA(4)*TC**3+ZA(5)*TC**4 PC=PA(1)+PA(2)*ZC+PA(3)*ZC**2+PA(4)*ZC**3+PA(5)*ZC**4 VC=VA(1)+VA(2)*ZC+VA(3)*ZC**2+VA(4)*ZC**3+VA(5)*ZC**4 J=0 RETURN ELSE J=-2 RETURN ENDIF END C---------------------------------------- C AKG CONVERSES [KMOL/KMOL] TO [KG/KG] C---------------------------------------- REAL FUNCTION AKG(X) REAL X DOUBLE PRECISION XX,FMF024,KG IF(X.LT.0.0.OR.X.GT.1.0) THEN CALL FMF020('AKG ') AKG=-1.0E20 RETURN ENDIF XX=DBLE(X) KG=FMF024(XX) AKG=REAL(KG) RETURN END C C*************************************************************** C C------------------------------------------ C AKMOL CONVERSES [KG/KG] TO [KMOL/KMOL] C------------------------------------------ REAL FUNCTION AKMOL(X) REAL X DOUBLE PRECISION XX,FMF025,KMOL IF(X.LT.0.0.OR.X.GT.1.0) THEN CALL FMF020('AKMOL ') AKMOL=-1.0E20 RETURN ENDIF XX=DBLE(X) KMOL=FMF025(XX) AKMOL=REAL(KMOL) RETURN END C C************************************************************** C C C------------------------------------------------------ C CRPM OUTPUTS CRITICAL PROPERTY OF PURE SUBSTANCE. C------------------------------------------------------ REAL FUNCTION CRPM(I,A) DOUBLE PRECISION TC,PC,VC,HC,SC,Z,FMF072 CHARACTER A DOUBLE PRECISION PR(1:11),CP(1:20) INTEGER COMBI COMMON/FMFC/PR,CP,COMBI C IF(I.EQ.1) THEN Z=1.0 TC=PR(2) PC=PR(3) VC=PR(4) ELSEIF(I.EQ.2) THEN Z=0.0 TC=PR(7) PC=PR(8) VC=PR(9) ELSE CALL FMF020('CRPM ') CRPM=-1.0E20 RETURN ENDIF C IF(A.EQ.'T') THEN CRPM=FMF055(1,TC) ELSEIF(A.EQ.'P') THEN CRPM=FMF055(2,PC) ELSEIF(A.EQ.'V') THEN CRPM=FMF056(1,VC,Z) ELSEIF(A.EQ.'H') THEN HC=FMF072(A,4,TC,PC,Z) CRPM=FMF056(2,HC,Z) ELSEIF(A.EQ.'S') THEN SC=FMF072(A,4,TC,PC,Z) CRPM=FMF056(3,SC,Z) ELSE CALL FMF020('CRPM ') CRPM=-1.0E20 ENDIF RETURN END C C-------------------------------------- C FCM OUTPUTS FUNDAMENTAL CONSTANTS. C-------------------------------------- REAL FUNCTION FCM(I,A) INTEGER I,COMBI CHARACTER A DOUBLE PRECISION PR(1:11),CP(1:20) COMMON/FMFC/PR,CP,COMBI C IF(A.EQ.'M') THEN IF(I.EQ.1) THEN FCM=REAL(PR(1)) ELSEIF(I.EQ.2) THEN FCM=REAL(PR(6)) ELSE CALL FMF020('FCM ') FCM=-1.0E20 ENDIF RETURN C ELSEIF(A.EQ.'R') THEN IF(I.EQ.1) THEN FCM=8314.51/REAL(PR(1)) ELSEIF(I.EQ.2) THEN FCM=8314.51/REAL(PR(6)) ELSE CALL FMF020('FCM ') FCM=-1.0E20 ENDIF RETURN C ELSEIF(A.EQ.'T') THEN IF(I.EQ.1) THEN FCM=FMF055(1,PR(2)) ELSEIF(I.EQ.2) THEN FCM=FMF055(1,PR(7)) ELSE CALL FMF020('FCM ') FCM=-1.0E20 ENDIF RETURN C ELSEIF(A.EQ.'P') THEN IF(I.EQ.1) THEN FCM=REAL(FMF055(2,PR(3))) ELSEIF(I.EQ.2) THEN FCM=REAL(FMF055(2,PR(8))) ELSE CALL FMF020('FCM ') FCM=-1.0E20 ENDIF RETURN C ELSEIF(A.EQ.'V') THEN IF(I.EQ.1) THEN FCM=FMF056(1,PR(4),1.0D0) ELSEIF(I.EQ.2) THEN FCM=FMF056(1,PR(9),0.0D0) ELSE CALL FMF020('FCM ') FCM=-1.0E20 ENDIF RETURN C ELSEIF(A.EQ.'E') THEN IF(I.EQ.1) THEN FCM=REAL(PR(5)) ELSEIF(I.EQ.2) THEN FCM=REAL(PR(10)) ELSE CALL FMF020('FCM ') FCM=-1.0E20 ENDIF RETURN C ELSE CALL FMF020('FCM ') FCM=-1.0E20 RETURN ENDIF END C C************************************************************ C Added by Akasaka, June 2, 1998 C-------------------------------------------------------- C IDENTM RETUENS NAME OF PURE COMPONENT C-------------------------------------------------------- CHARACTER*40 FUNCTION IDENTM( I, A ) CHARACTER PRNAME(2)*16, PRCHEM(2)*16 CHARACTER*1 A COMMON/NAME/PRNAME, PRCHEM IF ( I.EQ.1.OR.I.EQ.2 ) THEN IF ( A.EQ.'C' ) THEN IDENTM = PRCHEM( I ) ELSEIF ( A.EQ.'S' ) THEN IDENTM = PRNAME( I ) ELSE IDENTM = '12.1' ENDIF ELSE IDENTM = 'Out of range' ENDIF RETURN END C C************************************************************ C C-------------------------------------------------------- C IPHASE DECIDES STATE C COMPRESSED WATER , WET VAPOR , SUPERHEATED VAPOR C-------------------------------------------------------- INTEGER FUNCTION IPHASE(T,P,Z) REAL T,P,Z DOUBLE PRECISION TT,PP,ZZ,FMF022 C INTEGER FMF002 C modified Akasaka June 3, 1998 INTEGER FMF002, II c TT=FMF022(1,T) PP=FMF022(2,P) ZZ=FMF022(3,Z) C IPHASE=FMF002(TT,PP,ZZ) C modified Akasaka June 3, 1998 II = FMF002( TT, PP, ZZ ) IF ( II.EQ.1.OR.II.EQ.3 ) THEN IPHASE = 1 ELSE IF ( II.EQ.2 ) THEN IPHASE = 2 ELSE IPHASE = II ENDIF RETURN END C C********************************************************* C SUBROUTINE KPAMES(KPAA,MESSS) INTEGER KPAA,KPA,MESSS,MESS,KSTAN,KAS COMMON/UNIT/KPA,MESS,KSTAN,KAS c KPA=KPAA MESS=MESSS RETURN END C-------------------------------------------------- C MKTABL MAKES TABLE OF FUNDAMETAL CONSTANTS C AND IDEAL GAS HEAT CAPACITY. C-------------------------------------------------- SUBROUTINE MKTABL(J,KOMBI) REAL MW1,MW2,KIJ INTEGER J,KOMBI,COMBI CHARACTER*16 PRNAME(2),PRCHEM(2),NAME1,NAME2,CHEM1,CHEM2 DOUBLE PRECISION PR(1:11),CP(1:20) COMMON/FMFC/PR,CP,COMBI COMMON/NAME/PRNAME,PRCHEM C CALL START1(J,KOMBI) C IF(J.EQ.0) THEN C --- FUNDAMENTAL CONSTANTS ------------------------------------ NAME1=PRNAME(1) NAME2=PRNAME(2) CHEM1=PRCHEM(1) CHEM2=PRCHEM(2) MW1=REAL(PR(1)) MW2=REAL(PR(6)) R1=8314.51/MW1 R2=8314.51/MW2 TC1=REAL(PR(2)) TC2=REAL(PR(7)) PC1=REAL(PR(3))*1.0E-6 PC2=REAL(PR(8))*1.0E-6 VC1=REAL(PR(4)) VC2=REAL(PR(9)) OM1=REAL(PR(5)) OM2=REAL(PR(10)) KIJ=PR(11) WRITE(*,*) 'FUNDAMENTAL CONSTANTS' WRITE(*,100) NAME1,NAME2 WRITE(*,101) & ' CHEMICAL FORMULA ',CHEM1,CHEM2 WRITE(*,110)' MOLECULER WEIGHT [KG/KMOL] ',MW1,MW2 WRITE(*,110)' GAS CONSTANT [J/(KG･K)]',R1,R2 WRITE(*,120)' CRITICAL TEMPERATURE [K] ',TC1,TC2 WRITE(*,130)' CRITICAL PRESSURE [MPA] ',PC1,PC2 WRITE(*,140)' CRITICAL VOLUME [M^3/KMOL]',VC1,VC2 WRITE(*,150)' ACENTRIC FACTOR [-] ',OM1,OM2 WRITE(*,160)' INTERACTION PATAMETER[-] ',KIJ 100 FORMAT(32X,A,3X,A) 101 FORMAT(2A,3X,A) 110 FORMAT(A,F16.3,3X,F16.3) 120 FORMAT(A,F16.2,3X,F16.2) 130 FORMAT(A,F16.4,3X,F16.4) 140 FORMAT(A,F16.6,3X,F16.6) 150 FORMAT(A,F16.4,3X,F16.4) 160 FORMAT(A,17X,F7.4) C NEQ1=INT(CP(1)) A1=REAL(CP(2)) B1=REAL(CP(3)) C1=REAL(CP(4)) D1=REAL(CP(5)) E1=REAL(CP(6)) TMIN1=REAL(CP(7)) TMAX1=REAL(CP(8)) NEQ2=INT(CP(11)) A2=REAL(CP(12)) B2=REAL(CP(13)) C2=REAL(CP(14)) D2=REAL(CP(15)) E2=REAL(CP(16)) TMIN2=REAL(CP(17)) TMAX2=REAL(CP(18)) WRITE(*,*) ' ' WRITE(*,*) 'IDEAL GAS HEAT CAPACITY' WRITE(*,100) NAME1,NAME2 WRITE(*,200)' EQUATION NUMBER ',NEQ1,NEQ2 WRITE(*,210)' A ',A1,A2 WRITE(*,210)' B ',B1,B2 WRITE(*,210)' C ',C1,C2 WRITE(*,210)' D ',D1,D2 WRITE(*,210)' E ',E1,E2 WRITE(*,220)' LOW TEMPERATURE LIMIT [K] ',TMIN1,TMIN2 WRITE(*,220)' HIGH TEMPERATURE LIMIT [K] ',TMAX1,TMAX2 200 FORMAT(A,I16,3X,I16) 210 FORMAT(A,E18.6,1X,E18.6) 220 FORMAT(A,F16.1,3X,F16.1) J=0 RETURN C ELSE J=3 RETURN ENDIF END C C******************************************************* C C------------------------------------------------- C PBT CALCULATES PRESSURE AT BUBLE POINT. C------------------------------------------------- REAL FUNCTION PBT(T,Z) REAL T,Z,FMF055 DOUBLE PRECISION FMF022,TT,ZZ,PP,LIM,TC,PC,VC,FMF026,TC1 INTEGER J,JJ,COMBI DOUBLE PRECISION PR(1:11),CP(1:20) COMMON/FMFC/PR,CP,COMBI C IF(Z.LT.0.0.OR.Z.GT.1.0) THEN CALL FMF020('PBT ') PBT=-1.0E20 RETURN ENDIF C TT=FMF022(1,T) ZZ=FMF022(3,Z) C LIM=1.0D-7 J=0 IF(ZZ.LT.LIM) THEN ZZ=0.0 CALL FMF059(2,J,TT,PP) ELSEIF(ZZ.GT.1.0-LIM) THEN ZZ=1.0 CALL FMF059(1,J,TT,PP) ELSE TC1=PR(2) IF(TT.LE.TC1) THEN PP=FMF026(TT,ZZ) IF(PP.LT.0.0) J=-1 ELSE CALL FMF074(1,JJ,ZZ,TC,PC,VC) IF(TT.LT.TC) THEN PP=FMF026(TT,ZZ) IF(PP.LT.0.0) J=-1 ELSE J=-2 ENDIF ENDIF ENDIF C IF(J.EQ.0) THEN PBT=FMF055(2,PP) ELSEIF(J.EQ.-1) THEN CALL FMF019('PBT ') PBT=-1.0E10 ELSE CALL FMF020('PBT ') PBT=-1.0E20 ENDIF RETURN END C C************************************************************** C C-------------------------------------------- C PDT CALCULATES PRESSURE AT DEW POINT. C-------------------------------------------- REAL FUNCTION PDT(T,Z) REAL T,Z,FMF055 DOUBLE PRECISION FMF022,TT,ZZ,PP,LIM,TC,PC,VC,FMF031,TC1 INTEGER J,JJ,COMBI DOUBLE PRECISION PR(1:11),CP(1:20) COMMON/FMFC/PR,CP,COMBI C IF(Z.LT.0.0.OR.Z.GT.1.0) THEN CALL FMF020('PDT ') PDT=-1.0E20 RETURN ENDIF C TT=FMF022(1,T) ZZ=FMF022(3,Z) C LIM=1.0D-7 J=0 IF(ZZ.LT.LIM) THEN ZZ=0.0 CALL FMF059(2,J,TT,PP) ELSEIF(ZZ.GT.1.0-LIM) THEN ZZ=1.0 CALL FMF059(1,J,TT,PP) ELSE TC1=PR(2) IF(TT.LE.TC1) THEN PP=FMF031(TT,ZZ) IF(PP.LT.0.0) J=-1 ELSE CALL FMF074(1,JJ,ZZ,TC,PC,VC) IF(TT.LT.TC) THEN PP=FMF031(TT,ZZ) IF(PP.LT.0.0) J=-1 ELSE J=-2 ENDIF ENDIF ENDIF C IF(J.EQ.0) THEN PDT=FMF055(2,PP) ELSEIF(J.EQ.-1) THEN CALL FMF019('PDT ') PDT=-1.0E10 ELSE CALL FMF020('PDT ') PDT=-1.0E20 ENDIF RETURN END C---------------------------------------------------------------- C PSTM CALCULATES SATURATED VAPOR PRESSURE OF PURE SUBSTANCE C---------------------------------------------------------------- REAL FUNCTION PSTM(I,T) DOUBLE PRECISION TT,PP,FMF022 REAL T,FMF055 INTEGER J C IF(I.EQ.1.OR.I.EQ.2) THEN TT=FMF022(1,T) CALL FMF059(I,J,TT,PP) IF(J.EQ.0) THEN PSTM=FMF055(2,PP) ELSEIF(J.EQ.-1) THEN CALL FMF019('PSTM ') PSTM=-1.0E10 ELSEIF(J.EQ.-2) THEN CALL FMF020('PSTM ') PSTM=-1.0E20 ENDIF RETURN ELSE CALL FMF020('PSTM ') PSTM=-1.0E20 RETURN ENDIF END C C************************************************************** C C------------------------------------------------------------------- C QMIX CALCULATES QUALITY AT VAPOR-LIQUID EQUILIBRIUM STATE. C------------------------------------------------------------------- REAL FUNCTION QMIX(T,P,Z) REAL T,P,Z DOUBLE PRECISION TT,PP,ZZ,QQ,FMF022,LIM INTEGER J,IP,FMF002 c LIM=1.0D-7 C TT=FMF022(1,T) PP=FMF022(2,P) ZZ=FMF022(3,Z) IF(ZZ.LT.LIM.OR.ZZ.GT.1.0-LIM) THEN J=-2 CALL FMF020('QMIX ') QMIX=-1.0E20 RETURN ENDIF C IP=FMF002(TT,PP,ZZ) C IF(IP.NE.2) THEN J=-2 CALL FMF020('QMIX ') QMIX=-1.0E20 RETURN ENDIF C CALL FMF042(J,TT,PP,ZZ,QQ) IF(J.EQ.-1) THEN CALL FMF019('QMIX ') QMIX=-1.0E10 ELSE QMIX=REAL(QQ) ENDIF RETURN END C C************************************************************ C C---------------------------------------------------------------------- C START1 SET INTO PR AND CP. WHERE PR AND CP ARE C DOUBLE PRECISION ARRAY. C C PR C 1 : MW1 2 : TC1 3 : PC1 4 : VC1 5 : OMEGA1 C 6 : MW2 7 : TC2 8 : PC2 9 : VC2 10 : OMEGA2 C 11 : KIJ C C CP C CP = A + B[(C/T)SINH(C/T)]^2 + D[(E/T)COSH(E/T)]^2 (1) C CP = A + B*T + C*T^2 + D*T^3 + E*T^4 (2) C 1 : EQUATION FORM FOR 1 (1 OR 2) 2 : A FOR 1 3 : B FOR 1 C 4 : C FOR 1 5 : D FOR 1 6 : E FOR 1 7 : TMIN1 8 : TMAX1 C 9 : H01 10 : S01 C 11 : EQUATION FORM FOR2 (1 OR 2) 12 : A FOR 2 13 : B FOR 2 C 14 : C FOR 2 15 : D FOR 2 16 : E FOR 2 17 : TMIN2 18 : TMAX2 C 19 : H02 20 : S02 C C---------------------------------------------------------------------- SUBROUTINE START1(J,KOMBI) DOUBLE PRECISION R22(5), R123(5), R32(5), R125(5), & R134A(5), CH4(5), C2H4(5), C2H6(5), & C3H6(5), C3H8(5), IC4H10(5), NC4H10(5), & IC5H12(5),NC5H12(5), NC6H14(5), C6H6(5), & CC6H12(5),NC7H16(5), NC8H18(5), N2(5), & CO(5), CO2(5), H2S(5), O2(5), & CPR22(10), CPR123(10), CPR32(10), CPR125(10), & C134A(10), CPCH4(10), CPC2H4(10), CPC2H6(10), & CPC3H6(10), CPC3H8(10), CIC4HT(10),CNC4HT(10), & CIC5HW(10),CNC5HW(10),CNC6HF(10),CPC6H6(10), & CCC6HW(10),CNC7HS(10),CNC8HE(10),CPN2(10), & CPCO(10), CPCO2(10), CPH2S(10), CPO2(10) DOUBLE PRECISION KIJ CHARACTER PRNAME(2)*16,PRCHEM(2)*16 INTEGER KOMBI,I,J,COMBI DOUBLE PRECISION PR(1:11),CP(1:20) COMMON/FMFC/PR,CP,COMBI COMMON/NAME/PRNAME,PRCHEM C C C--- THERMODYNAMICS PROPERTIES OF PURE SUBSTANCES (MW,TC,PC,VC,OMEGA) --- C DATA(R22(I),I=1,5)/86.468,369.30,4.971D6,0.168615,0.2192/ C* DATA(R123(I),I=1,5)/152.931,456.86,3.666D6,0.275276,0.2816/ C* DATA(R32(I),I=1,5)/52.024,351.60,5.8302D6,0.12100,0.2763/ C* DATA(R125(I),I=1,5)/120.020,339.4,3.631D6,0.20983,0.306 / C* DATA(R134A(I),I=1,5)/102.030,374.21,4.056D6,0.19812,0.326 / C* DATA(CH4(I),I=1,5)/16.043,190.58,4.6043D6,0.09925,0.0108/ C* DATA(C2H4(I),I=1,5)/28.054,282.36,5.0318D6,0.12907,0.0852/ C* DATA(C2H6(I),I=1,5)/30.070,305.42,4.8801D6,0.14792,0.0990/ C* DATA(C3H6(I),I=1,5)/42.081,364.76,4.6126D6,0.18100,0.1424/ C* DATA(C3H8(I),I=1,5)/44.096,369.82,4.2492D6,0.20288,0.1518/ C* DATA(IC4H10(I),I=1,5)/58.123,408.14,3.6480D6,0.26270,0.1770/ C* DATA(NC4H10(I),I=1,5)/58.123,425.18,3.7969D6,0.25490,0.1993/ C* DATA(IC5H12(I),I=1,5)/72.150,460.43,3.3812D6,0.30583,0.2275/ C* DATA(NC5H12(I),I=1,5)/72.150,433.78,3.1992D6,0.30358,0.1964/ C* DATA(NC6H14(I),I=1,5)/86.177,507.43,3.0123D6,0.36990,0.3046/ C* DATA(C6H6(I),I=1,5)/78.114,562.16,4.8980D6,0.25894,0.2108/ C* DATA(CC6H12(I),I=1,5)/84.161,553.54,4.0748D6,0.30788,0.2118/ C* DATA(NC7H16(I),I=1,5)/100.204,540.26,2.7358D6,0.43190,0.3511/ C* DATA(NC8H18(I),I=1,5)/114.231,568.83,2.4863D6,0.49205,0.3962/ C* DATA(N2(I),I=1,5)/28.014,126.10,3.3944D6,0.09010,0.0403/ C* DATA(CO(I),I=1,5)/28.010,132.92,3.4988D6,0.09310,0.0663/ C* C DATA(CO2(I),I=1,5)/44.010,304.19,7.3815D6,0.09400,0.2276/ DATA(CO2(I),I=1,5)/44.010,304.21,7.3830D6,0.09400,0.2276/ C* DATA(H2S(I),I=1,5)/34.082,373.53,8.9629D6,0.09849,0.0827/ C* DATA(O2(I),I=1,5)/31.999,154.58,5.0430D6,0.07340,0.0218/ C C*--- IDEAL GAS HEAT CAPACITY ----------------------- C DATA(CPR22(I),I=1,10) & /1.0D0,3.4600D4,7.0840D4,9.8750D2,4.3500D4,4.3500D2, & 100.0,1500,430.06567,1.9847835/ c & 100.0,1500,430.06567,1.98457026/ C* DATA(CPR123(I),I=1,10) & /2.0D0,6.02126D3,0.50517D3,-7.6316D-1,5.1834D-4,0.0D0, c $ 150.0,500.0,395.86546,1.7135793/ & 150.0,500.0,395.86546,1.65377557/ C* DATA(CPR32(I),I=1,10) & /1.0D0,3.4060D4,7.1440D4,1.4220D3,3.9400D4,6.7700D2, & 100.0,1500.0,565.2668,2.650756/ C* DATA(CPR125(I),I=1,10) & /2.0D0,2.36022D4,2.83723D2,-1.23028D-1,-5.67252D-5,0.0D0, & 150.0,500.0,363.0407,1.584097/ C* DATA(C134A(I),I=1,10) & /2.0D0,1.67813D4,2.86357D2,-2.27337D-1,1.133121D-4,0.0D0, & 150.0,500.0,426.3636,1.819680/ C* DATA(CPCH4(I),I=1,10) & /1.0D0,3.3295D4,8.0295D4,2.1018D3,4.2130D4,9.9510D2, & 50.0,1500.0,0.0D0,0.0D0/ C* DATA(CPC2H4(I),I=1,10) & /1.0D0,3.3380D4,9.4790D4,1.5960D3,5.5100D4,7.4080D2, & 60.0,1500.0,0.0D0,0.0D0/ C* DATA(CPC2H6(I),I=1,10) & /1.0D0,3.5650D4,1.3520D5,1.4300D3,6.1800D4,6.1200D2, & 100.0,1500.0,613.6744,3.308224/ C DATA(CPC3H6(I),I=1,10) & /1.0D0,4.1300D4,1.5250D5,1.3520D3,7.4400D4,5.7800D2, & 90.0,1500.0,0.0D0,0.0D0/ C DATA(CPC3H8(I),I=1,10) & /1.0D0,4.4000D4,1.93800D5,1.3690D3,9.8000D4,5.8300D2, & 100.0,1500.0,0.0D0,0.0D0/ C DATA(CIC4HT(I),I=1,10) & /1.0D0,6.5490D4,2.4776D5,1.5870D3,1.5750D5,-7.0699D2, & 200.0,1500.0,0.0D0,0.0D0/ C DATA(CNC4HT(I),I=1,10) & /1.0D0,7.1340D4,2.4300D5,1.6300D3,1.5033D5,-7.3042D2, & 200.0,1500.0,0.0D0,0.0D0/ C DATA(CIC5HW(I),I=1,10) & /1.0D0,7.6400D4,3.3030D5,1.6190D3,2.0100D5,-6.9040D2, & 200.0,1500.0,0.0D0,0.0D0/ C DATA(CNC5HW(I),I=1,10) & /1.0D0,7.6200D4,3.6870D5,1.5550D3,2.1200D5,-6.3290D2, & 200.0,1500.0,0.0D0,0.0D0/ C DATA(CNC6HF(I),I=1,10) & /1.0D0,1.0440D5,3.5230D5,1.6946D3,2.3690D5,-7.6160D2, & 200.0,1500.0,0.0D0,0.0D0/ C DATA(CPC6H6(I),I=1,10) & /1.0D0,4.4420D4,2.3205D5,1.4946D3,1.7213D5,-6.7815D2, & 200.0,1500.0,0.0D0,0.0D0/ C DATA(CCC6HW(I),I=1,10) & /1.0D0,4.3200D4,3.7350D5,1.1920D3,1.6350D5,-5.3010D2, & 100.0,1500.0,0.0D0,0.0D0/ C* DATA(CNC7HS(I),I=1,10) & /1.0D0,1.2015D5,4.0010D5,1.6766D3,2.7400D5,-7.5640D2, & 200.0,1500.0,0.0D0,0.0D0/ C* DATA(CNC8HE(I),I=1,10) & /1.0D0,1.3554D5,4.4310D5,1.63560D3,3.0540D5,7.4640D2, & 200.0,1500.0,0.0D0,0.0D0/ C* DATA(CPN2(I),I=1,10) & /1.0D0,2.9105D4,8.6149D3,1.7016D3,1.0347D2,9.0979D2, & 50.0,1500.0,0.0D0,0.0D0/ C* DATA(CPCO(I),I=1,10) & /1.0D0,2.9108D4,8.7730D3,3.0851D3,8.4553D3,1.5382D3, & 60.0,1500.0,0.0D0,0.0D0/ C* DATA(CPCO2(I),I=1,10) & /1.0D0,2.9370D4,3.45400D4,-1.4280D3,2.6400D4,5.8800D2, & 50.0,5000.0,0.0D0,0.0D0/ C* DATA(CPH2S(I),I=1,10) & /1.0D0,3.3288D4,2.6086D4,9.1340D2,-1.7979D4,9.4940D2, & 100.0,1500.0,0.0D0,0.0D0/ C* DATA(CPO2(I),I=1,10) & /1.0D0,2.9103D4,1.0040D4,2.5265D3,9.3560D3,1.1538D3, & 50.0,1500.0,0.0D0,0.0D0/ C C--- SET CONSTANTS INTO PR AND CP ---------- C COMBI=KOMBI C --- R22-R123 ----------------------------- IF(KOMBI.EQ.1) THEN J=0 KIJ=0.0D0 DO 10 I=1,5 PR(I)=R22(I) 10 CONTINUE DO 11 I=1,5 PR(I+5)=R123(I) 11 CONTINUE PR(11)=KIJ DO 12 I=1,10 CP(I)=CPR22(I) 12 CONTINUE DO 13 I=1,10 CP(I+10)=CPR123(I) 13 CONTINUE PRNAME(1)=' R22' PRNAME(2)=' R123' PRCHEM(1)=' CHClF2' PRCHEM(2)=' CHCl2CF3' RETURN C C --- R32-R125 ----------------------------- ELSEIF(KOMBI.EQ.2) THEN J=0 KIJ=0.013 DO 20 I=1,5 PR(I)=R32(I) 20 CONTINUE DO 21 I=1,5 PR(I+5)=R125(I) 21 CONTINUE PR(11)=KIJ DO 22 I=1,10 CP(I)=CPR32(I) 22 CONTINUE DO 23 I=1,10 CP(I+10)=CPR125(I) 23 CONTINUE PRNAME(1)=' R32' PRNAME(2)=' R125' PRCHEM(1)=' CH2F2' PRCHEM(2)=' CF3CHF2' RETURN C C --- R32-R134A ----------------------------- ELSEIF(KOMBI.EQ.3) THEN J=0 KIJ=-0.01816 DO 30 I=1,5 PR(I)=R32(I) 30 CONTINUE DO 31 I=1,5 PR(I+5)=R134A(I) 31 CONTINUE PR(11)=KIJ DO 32 I=1,10 CP(I)=CPR32(I) 32 CONTINUE DO 33 I=1,10 CP(I+10)=C134A(I) 33 CONTINUE PRNAME(1)=' R32' PRNAME(2)=' R134a' PRCHEM(1)=' CH2F2' PRCHEM(2)=' CF3CH2F' RETURN C C--- CH4-C2H4 ------------------------------- ELSEIF(KOMBI.EQ.4) THEN J=0 KIJ=0.022 DO 40 I=1,5 PR(I)=CH4(I) 40 CONTINUE DO 41 I=1,5 PR(I+5)=C2H4(I) 41 CONTINUE PR(11)=KIJ DO 42 I=1,10 CP(I)=CPCH4(I) 42 CONTINUE DO 43 I=1,10 CP(I+10)=CPC2H4(I) 43 CONTINUE PRNAME(1)=' METHANE' PRNAME(2)=' ETHYLENE' PRCHEM(1)=' CH4' PRCHEM(2)=' C2H4' RETURN C C--- CH4-C2H6 ------------------------------- ELSEIF(KOMBI.EQ.5) THEN J=0 KIJ=-0.003 DO 50 I=1,5 PR(I)=CH4(I) 50 CONTINUE DO 51 I=1,5 PR(I+5)=C2H6(I) 51 CONTINUE PR(11)=KIJ DO 52 I=1,10 CP(I)=CPCH4(I) 52 CONTINUE DO 53 I=1,10 CP(I+10)=CPC2H6(I) 53 CONTINUE PRNAME(1)=' METHANE' PRNAME(2)=' ETHANE' PRCHEM(1)=' CH4' PRCHEM(2)=' C2H6' RETURN C C--- CH4-C3H6 ------------------------------- ELSEIF(KOMBI.EQ.6) THEN J=0 KIJ=0.033 DO 60 I=1,5 PR(I)=CH4(I) 60 CONTINUE DO 61 I=1,5 PR(I+5)=C3H6(I) 61 CONTINUE PR(11)=KIJ DO 62 I=1,10 CP(I)=CPCH4(I) 62 CONTINUE DO 63 I=1,10 CP(I+10)=CPC3H6(I) 63 CONTINUE PRNAME(1)=' METHANE' PRNAME(2)=' PROPYLENE' PRCHEM(1)=' CH4' PRCHEM(2)=' C3H6' RETURN C C--- CH4-C3H8 ------------------------------- ELSEIF(KOMBI.EQ.7) THEN J=0 KIJ=0.016 DO 70 I=1,5 PR(I)=CH4(I) 70 CONTINUE DO 71 I=1,5 PR(I+5)=C3H8(I) 71 CONTINUE PR(11)=KIJ DO 72 I=1,10 CP(I)=CPCH4(I) 72 CONTINUE DO 73 I=1,10 CP(I+10)=CPC3H8(I) 73 CONTINUE PRNAME(1)=' METHANE' PRNAME(2)=' PROPANE' PRCHEM(1)=' CH4' PRCHEM(2)=' C3H8' RETURN C C--- CH4-i-C4H10 ------------------------------- ELSEIF(KOMBI.EQ.8) THEN J=0 KIJ=0.026 DO 80 I=1,5 PR(I)=CH4(I) 80 CONTINUE DO 81 I=1,5 PR(I+5)=IC4H10(I) 81 CONTINUE PR(11)=KIJ DO 82 I=1,10 CP(I)=CPCH4(I) 82 CONTINUE DO 83 I=1,10 CP(I+10)=CIC4HT(I) 83 CONTINUE PRNAME(1)=' METHANE' PRNAME(2)=' ISOBUTANE' PRCHEM(1)=' CH4' PRCHEM(2)=' i-C4H10' RETURN C C--- CH4-n-C4H10 ------------------------------- ELSEIF(KOMBI.EQ.9) THEN J=0 KIJ=0.019 DO 90 I=1,5 PR(I)=CH4(I) 90 CONTINUE DO 91 I=1,5 PR(I+5)=NC4H10(I) 91 CONTINUE PR(11)=KIJ DO 92 I=1,10 CP(I)=CPCH4(I) 92 CONTINUE DO 93 I=1,10 CP(I+10)=CNC4HT(I) 93 CONTINUE PRNAME(1)=' METHANE' PRNAME(2)=' n-BUTANE' PRCHEM(1)=' CH4' PRCHEM(2)=' n-C4H10' RETURN C C--- CH4-i-C5H12 ------------------------------- ELSEIF(KOMBI.EQ.10) THEN J=0 KIJ=-0.006 DO 100 I=1,5 PR(I)=CH4(I) 100 CONTINUE DO 101 I=1,5 PR(I+5)=IC5H12(I) 101 CONTINUE PR(11)=KIJ DO 102 I=1,10 CP(I)=CPCH4(I) 102 CONTINUE DO 103 I=1,10 CP(I+10)=CIC5HW(I) 103 CONTINUE PRNAME(1)=' METHANE' PRNAME(2)=' ISOPENTAN' PRCHEM(1)=' CH4' PRCHEM(2)=' i-C5H12' RETURN ENDIF C C--- CH4-n-C5H12 ------------------------------- IF(KOMBI.EQ.11) THEN J=0 KIJ=0.026 DO 110 I=1,5 PR(I)=CH4(I) 110 CONTINUE DO 111 I=1,5 PR(I+5)=NC5H12(I) 111 CONTINUE PR(11)=KIJ DO 112 I=1,10 CP(I)=CPCH4(I) 112 CONTINUE DO 113 I=1,10 CP(I+10)=CNC5HW(I) 113 CONTINUE PRNAME(1)=' METHANE' PRNAME(2)=' NEOPENTAN' PRCHEM(1)=' CH4' PRCHEM(2)=' n-C5H12' RETURN C C--- CH4-n-C6H14 ------------------------------- ELSEIF(KOMBI.EQ.12) THEN J=0 KIJ=0.040 DO 120 I=1,5 PR(I)=CH4(I) 120 CONTINUE DO 121 I=1,5 PR(I+5)=NC6H14(I) 121 CONTINUE PR(11)=KIJ DO 122 I=1,10 CP(I)=CPCH4(I) 122 CONTINUE DO 123 I=1,10 CP(I+10)=CNC6HF(I) 123 CONTINUE PRNAME(1)=' METHANE' PRNAME(2)=' n-HEXANE' PRCHEM(1)=' CH4' PRCHEM(2)=' n-C6H14' RETURN C C--- CH4-C6H6 ------------------------------- ELSEIF(KOMBI.EQ.13) THEN J=0 KIJ=0.055 DO 130 I=1,5 PR(I)=CH4(I) 130 CONTINUE DO 131 I=1,5 PR(I+5)=C6H6(I) 131 CONTINUE PR(11)=KIJ DO 132 I=1,10 CP(I)=CPCH4(I) 132 CONTINUE DO 133 I=1,10 CP(I+10)=CPC6H6(I) 133 CONTINUE PRNAME(1)=' METHANE' PRNAME(2)=' BENZENE' PRCHEM(1)=' CH4' PRCHEM(2)=' C6H6' RETURN C C--- CH4-c-C6H12 ------------------------------- ELSEIF(KOMBI.EQ.14) THEN J=0 KIJ=0.039 DO 140 I=1,5 PR(I)=CH4(I) 140 CONTINUE DO 141 I=1,5 PR(I+5)=CC6H12(I) 141 CONTINUE PR(11)=KIJ DO 142 I=1,10 CP(I)=CPCH4(I) 142 CONTINUE DO 143 I=1,10 CP(I+10)=CCC6HW(I) 143 CONTINUE PRNAME(1)=' METHANE' PRNAME(2)=' CYCLOHEXANE' PRCHEM(1)=' CH4' PRCHEM(2)=' c-C6H12' RETURN C C--- CH4-n-C7H16 ------------------------------- ELSEIF(KOMBI.EQ.15) THEN J=0 KIJ=0.035 DO 150 I=1,5 PR(I)=CH4(I) 150 CONTINUE DO 151 I=1,5 PR(I+5)=NC7H16(I) 151 CONTINUE PR(11)=KIJ DO 152 I=1,10 CP(I)=CPCH4(I) 152 CONTINUE DO 153 I=1,10 CP(I+10)=CNC7HS(I) 153 CONTINUE PRNAME(1)=' METHANE' PRNAME(2)=' n-HEPTANE' PRCHEM(1)=' CH4' PRCHEM(2)=' n-C7H16' RETURN C C--- CH4-n-C8H18 ------------------------------- ELSEIF(KOMBI.EQ.16) THEN J=0 KIJ=0.050 DO 160 I=1,5 PR(I)=CH4(I) 160 CONTINUE DO 161 I=1,5 PR(I+5)=NC8H18(I) 161 CONTINUE PR(11)=KIJ DO 162 I=1,10 CP(I)=CPCH4(I) 162 CONTINUE DO 163 I=1,10 CP(I+10)=CNC8HE(I) 163 CONTINUE PRNAME(1)=' METHANE' PRNAME(2)=' n-OCTANE' PRCHEM(1)=' CH4' PRCHEM(2)=' n-C8H18' RETURN C C--- O2-CO2 ------------------------------- ELSEIF(KOMBI.EQ.17) THEN J=0 KIJ=0.09D0 DO 170 I=1,5 PR(I)=O2(I) 170 CONTINUE DO 171 I=1,5 PR(I+5)=CO2(I) 171 CONTINUE PR(11)=KIJ DO 172 I=1,10 CP(I)=CPO2(I) 172 CONTINUE DO 173 I=1,10 CP(I+10)=CPCO2(I) 173 CONTINUE PRNAME(1)=' OXYGEN' PRNAME(2)=' CARBON DIOXIDE' PRCHEM(1)=' O2' PRCHEM(2)=' CO2' RETURN C ELSE J=3 CALL FMF021('START1') RETURN ENDIF END C C******************************************************** C SUBROUTINE START2(J,PR1,PR2,KKIJ,CP1,CP2) REAL PR1(5),PR2(5),CP1(8),CP2(8),A(5),B(8) REAL KKIJ INTEGER I,COMBI DOUBLE PRECISION PR(1:11),CP(1:20) COMMON/FMFC/PR,CP,COMBI C IF(PR1(2).GT.PR2(2)) THEN J=-1 DO 1000 I=1,5 A(I)=PR2(I) 1000 CONTINUE DO 1100 I=1,5 PR2(I)=PR1(I) 1100 CONTINUE DO 1200 I=1,5 PR1(I)=A(I) 1200 CONTINUE DO 1300 I=1,8 B(I)=CP2(I) 1300 CONTINUE DO 1400 I=1,8 CP2(I)=CP1(I) 1400 CONTINUE DO 1500 I=1,8 CP1(I)=B(I) 1500 CONTINUE ELSE J=0 ENDIF C C--- SET CONSTANTS INTO PR -------------- DO 2000 I=1,5 PR(I)=DBLE(PR1(I)) 2000 CONTINUE DO 2100 I=1,5 PR(I+5)=DBLE(PR2(I)) 2100 CONTINUE PR(11)=DBLE(KKIJ) C C--- SET CONSTANTS INTO CPR ------------- CP(1)=DBLE(CP1(1)) CP(2)=DBLE(CP1(2)) CP(3)=DBLE(CP1(3)) CP(4)=DBLE(CP1(4)) CP(5)=DBLE(CP1(5)) CP(6)=DBLE(CP1(6)) CP(7)=DBLE(CP1(7)) CP(8)=DBLE(CP1(8)) CP(9)=0.0D0 CP(10)=0.0D0 CP(11)=DBLE(CP2(1)) CP(12)=DBLE(CP2(2)) CP(13)=DBLE(CP2(3)) CP(14)=DBLE(CP2(4)) CP(15)=DBLE(CP2(5)) CP(16)=DBLE(CP2(6)) CP(17)=DBLE(CP2(7)) CP(18)=DBLE(CP2(8)) CP(19)=0.0D0 CP(20)=0.0D0 C RETURN END C C******************************************************* C SUBROUTINE STNKAS(KSTANN,KASS) INTEGER KASS,KAS,KSTANN,KSTAN,KPA,MESS COMMON/UNIT/KPA,MESS,KSTAN,KAS c KSTAN=KSTANN KAS=KASS RETURN END C C********************************************************* C C----------------------------------------------------- C SUBCRT CALCULATES PROPERTIES AT CRITICAL STATE. C----------------------------------------------------- SUBROUTINE SUBCRT(I,J,TC,PC,Z,VC,HC,SC) REAL TC,PC,Z,VC,HC,SC DOUBLE PRECISION LIM,ZZ,FMF022,TTC,PPC,VVC,HHC $ ,SSC,FMF072,TC1,TC2 DOUBLE PRECISION PR(1:11),CP(1:20) INTEGER COMBI COMMON/FMFC/PR,CP,COMBI C IF(I.EQ.1) THEN LIM=1.0D-5 IF(Z.LT.0.0.OR.Z.GT.1.0) THEN J=-2 CALL FMF020('SUBCRT') ER=-1.0E20 TC=ER PC=ER VC=ER HC=ER SC=ER RETURN ENDIF C ZZ=FMF022(3,Z) IF(ZZ.LT.LIM) THEN J=0 ZZ=0.0 TTC=PR(7) PPC=PR(8) VVC=PR(9) ELSEIF(ZZ.GT.1.0-LIM) THEN J=0 ZZ=1.0 TTC=PR(2) PPC=PR(3) VVC=PR(4) ELSE CALL FMF074(1,J,ZZ,TTC,PPC,VVC) ENDIF C IF(J.EQ.0) THEN HHC=FMF072('H',4,TTC,PPC,ZZ) SSC=FMF072('S',4,TTC,PPC,ZZ) TC=FMF055(1,TTC) PC=FMF055(2,PPC) VC=FMF056(1,VVC,ZZ) HC=FMF056(2,HHC,ZZ) SC=FMF056(3,SSC,ZZ) RETURN ELSE J=-1 CALL FMF019('SUBCRT') ER=-1.0E10 TC=ER PC=ER VC=ER HC=ER SC=ER RETURN ENDIF ENDIF C-------------------------------------------------------------- C IF(I.EQ.2) THEN TC1=PR(2) TC2=PR(7) TTC=FMF022(1,TC) IF(TTC.LT.TC1.OR.TTC.GT.TC2) THEN J=-2 CALL FMF020('SUBCRT') ER=-1.0E20 Z=ER PC=ER VC=ER HC=ER SC=ER RETURN ENDIF CALL FMF074(2,J,ZZ,TTC,PPC,VVC) IF(J.EQ.0) THEN HHC=FMF072('H',4,TTC,PPC,ZZ) SSC=FMF072('S',4,TTC,PPC,ZZ) Z=FMF055(3,ZZ) PC=FMF055(2,PPC) VC=FMF056(1,VVC,ZZ) HC=FMF056(2,HHC,ZZ) SC=FMF056(3,SSC,ZZ) RETURN ELSE J=-1 CALL FMF019('SUBCRT') ER=-1.0E10 Z=ER PC=ER VC=ER HC=ER SC=ER RETURN ENDIF ENDIF C-------------------------------------------------------------- C IF(I.EQ.3) THEN J=-2 CALL FMF020('SUBCRT') ER=-1.0E20 Z=ER TC=ER PC=ER VC=ER HC=ER SC=ER RETURN ENDIF END C C********************************************************** C C---------------------------------------------- C SUBMIX CALCULATES PROPERTIES OF MIXTURE. C---------------------------------------------- SUBROUTINE SUBMIX(I,J,T,P,Z,V,H,S) REAL T,P,Z,V,H,S,ER,FMF055,FMF056 DOUBLE PRECISION TT,PP,ZZ,VV,HH,SS,FMF022,FMF023,LIM C IF(Z.LT.0.0.OR.Z.GT.1.0) THEN CALL FMF020('SUBMIX') T=-1.0E20 P=-1.0E20 V=-1.0E20 H=-1.0E20 S=-1.0E20 RETURN ENDIF C LIM=1.0D-15 C IF(I.EQ.1) THEN TT=FMF022(1,T) PP=FMF022(2,P) ZZ=FMF022(3,Z) IF(ZZ.LT.LIM) THEN CALL FMF064(1,2,J,TT,PP,VV,HH,SS) ELSEIF(ZZ.GT.1.0-LIM) THEN CALL FMF064(1,1,J,TT,PP,VV,HH,SS) ELSE CALL FMF048(J,TT,PP,ZZ,VV,HH,SS) ENDIF IF(J.EQ.0) THEN V=FMF056(1,VV,ZZ) H=FMF056(2,HH,ZZ) S=FMF056(3,SS,ZZ) ELSE CALL FMF019('SUBMIX') ER=-1.0E10 V=ER H=ER S=ER ENDIF RETURN ENDIF C---------------------------------------------------------- IF(I.EQ.2) THEN PP=FMF022(2,P) ZZ=FMF022(3,Z) HH=FMF023(2,H,Z) IF(ZZ.LT.LIM) THEN CALL FMF064(2,2,J,TT,PP,VV,HH,SS) ELSEIF(ZZ.GT.1.0-LIM) THEN CALL FMF064(2,1,J,TT,PP,VV,HH,SS) ELSE CALL FMF049(J,TT,PP,ZZ,VV,HH,SS) ENDIF IF(J.EQ.0) THEN T=FMF055(1,TT) V=FMF056(1,VV,ZZ) S=FMF056(3,SS,ZZ) ELSEIF(J.EQ.-1) THEN CALL FMF019('SUBMIX') ER=-1.0E10 T=ER V=ER S=ER ELSE CALL FMF020('SUBMIX') ER=-1.0E20 T=ER V=ER S=ER ENDIF RETURN ENDIF C---------------------------------------------------------- IF(I.EQ.3) THEN PP=FMF022(2,P) ZZ=FMF022(3,Z) SS=FMF023(3,S,Z) IF(ZZ.LT.LIM) THEN CALL FMF064(3,2,J,TT,PP,VV,HH,SS) ELSEIF(ZZ.GT.1.0-LIM) THEN CALL FMF064(3,1,J,TT,PP,VV,HH,SS) ELSE CALL FMF051(J,TT,PP,ZZ,VV,HH,SS) ENDIF IF(J.EQ.0) THEN T=FMF055(1,TT) V=FMF056(1,VV,ZZ) H=FMF056(2,HH,ZZ) ELSEIF(J.EQ.-1) THEN CALL FMF019('SUBMIX') ER=-1.0E10 T=ER V=ER H=ER ELSE CALL FMF020('SUBMIX') ER=-1.0E20 T=ER V=ER H=ER ENDIF RETURN ENDIF C---------------------------------------------------------- IF(I.EQ.4) THEN PP=FMF022(2,P) ZZ=FMF022(3,Z) VV=FMF023(1,V,Z) IF(ZZ.LT.LIM) THEN CALL FMF064(4,2,J,TT,PP,VV,HH,SS) ELSEIF(ZZ.GT.1.0-LIM) THEN CALL FMF064(4,1,J,TT,PP,VV,HH,SS) ELSE CALL FMF053(J,TT,PP,ZZ,VV,HH,SS) ENDIF IF(J.EQ.0) THEN T=FMF055(1,TT) H=FMF056(2,HH,ZZ) S=FMF056(3,SS,ZZ) ELSEIF(J.EQ.-1) THEN CALL FMF019('SUBMIX') ER=-1.0E10 T=ER H=ER S=ER ELSE CALL FMF020('SUBMIX') ER=-1.0E20 T=ER H=ER S=ER ENDIF RETURN C---------------------------------------------------------- ELSE J=-2 CALL FMF020('SUBMIX') ER=-1.0E20 T=ER P=ER Z=ER V=ER H=ER S=ER RETURN ENDIF END C C***************************************************************** C SUBROUTINE SUBMXH and SUBMXT were written by Chen. C ( added by Akasaka June 3, 1998 ) C C***************************************************************** C ---------------------------------------------------------- C ** SUBROUTINE SUBMXH ** C ** Thermal Properties of Mixtures of Mixed by ** C ** Two Different Compositions ** C ** [Input] ** C ** PP : Pressure [Pa],[bar] ** C ** ZA : Total Composition of A ** C ** ZB : Total Composition of B ** C ** HA : Enthalpy of A [J/kmol],[J/kg] ** C ** HB : Enthalpy of B [J/kmol],[J/kg] ** C ** W : Fraction of A in Mixing ** C ** [Output] ** C ** J : Error Detection Code ** C ** TT : Temperature after Mixing ** C ** Z : Total Composition after Mixing ** C ** V : Volume after Mixing ** C ** H : Enthalpy after Mixing ** C ** S : Entropy after Mixing ** C ---------------------------------------------------------- SUBROUTINE SUBMXH(J,T,P,ZA,ZB,Z,V,HA,HB,H,S,W) INTEGER J,KPA,MESS,KSTAN,KAS REAL T,P,ZA,ZB,Z,V,HA,HB,H,S,W CHARACTER*6 PRNAME COMMON /UNIT/KPA,MESS,KSTAN,KAS PRNAME='SUBMXH' J=0 Z=W*ZA+(1-W)*ZB H=W*HA+(1-W)*HB CALL SUBMIX(2,J,T,P,Z,V,HH,S) RETURN END C ---------------------------------------------------------- C ** SUBROUTINE SUBMXT ** C ** Thermal Properties of Mixtures of Mixed by ** C ** Two Different Compositions ** C ** [Input] ** C ** PP : Pressure [Pa],[bar] ** C ** ZA : Total Composition of A ** C ** ZB : Total Composition of B ** C ** TA : Temperature of A [C],[K] ** C ** TB : Temperature of B [C],[K] ** C ** W : Fraction of A in Mixing ** C ** [Output] + ** C ** J : Error Detection Code ** C ** TT : Temperature after Mixing ** C ** Z : Total Composition after Mixing ** C ** V : Volume after Mixing ** C ** H : Enthalpy after Mixing ** C ** S : Entropy after Mixing ** C ---------------------------------------------------------- SUBROUTINE SUBMXT(J,T,P,ZA,ZB,Z,V,TA,TB,H,S,W) INTEGER J,KPA,MESS,KSTAN,KAS REAL T,P,ZA,ZB,Z,V,TA,TB,H,S,W,HA,HB,SA,SB CHARACTER*6 PRNAME COMMON /UNIT/KPA,MESS,KSTAN,KAS PRNAME='SUBMXT' J=0 Z=W*ZA+(1-W)*ZB CALL SUBMIX(1,J,TA,P,ZA,VA,HA,SA) CALL SUBMIX(1,J,TB,P,ZB,VB,HB,SB) H=W*HA+(1-W)*HB CALL SUBMIX(2,J,T,P,Z,V,H,S) RETURN END C********************************************************** C C---------------------------------------------- C SUBPAR CALCULATES PARTIAL PROPERTIES. C---------------------------------------------- SUBROUTINE SUBPAR(A,J,T,P,Z,XL,XV) REAL T,P,Z,XL(2),XV(2) DOUBLE PRECISION TT,PP,ZZ,XX,YY,VL(2),VV(2),HL(2),HV(2) $ ,SL(2),SV(2), & FMF022,LIM INTEGER KPA,MESS,KSTAN,KAS,COMBI CHARACTER A DOUBLE PRECISION PR(1:11),CP(1:20) COMMON/UNIT/KPA,MESS,KSTAN,KAS COMMON/FMFC/PR,CP,COMBI C TT=FMF022(1,T) PP=FMF022(2,P) ZZ=FMF022(3,Z) C LIM=1.0D-10 IF(ZZ.LT.LIM) ZZ=LIM IF(ZZ.GT.1.0-LIM) ZZ=1.0-LIM C ---------------------------------------------------- IF(A.EQ.'V') THEN YY=ZZ CALL FMF041(IP,J,TT,PP,XX,YY,VL,VV) IF(J.EQ.-1) THEN CALL FMF019('SUBPAR') XL(1)=-1.0E10 XL(2)=-1.0E10 XV(1)=-1.0E10 XV(2)=-1.0E10 RETURN ENDIF IF(IP.EQ.1) THEN IF(KAS.EQ.1) THEN XL(1)=REAL(VL(1)/PR(1)) XL(2)=REAL(VL(2)/PR(6)) ELSE XL(1)=REAL(VL(1)) XL(2)=REAL(VL(2)) ENDIF XV(1)=-1.0E20 XV(2)=-1.0E20 ELSEIF(IP.EQ.2) THEN IF(KAS.EQ.1) THEN XL(1)=REAL(VL(1)/PR(1)) XL(2)=REAL(VL(2)/PR(6)) XV(1)=REAL(VV(1)/PR(1)) XV(2)=REAL(VV(2)/PR(6)) ELSE XL(1)=REAL(VL(1)) XL(2)=REAL(VL(2)) XV(1)=REAL(VV(1)) XV(2)=REAL(VV(2)) ENDIF ELSE XL(1)=-1.0E20 XL(2)=-1.0E20 IF(KAS.EQ.1) THEN XV(1)=REAL(VV(1)/PR(1)) XV(2)=REAL(VV(2)/PR(6)) ELSE XV(1)=REAL(VV(1)) XV(2)=REAL(VV(2)) ENDIF ENDIF RETURN ENDIF C IF(A.EQ.'H') THEN YY=ZZ CALL FMF037(IP,J,TT,PP,XX,YY,HL,HV) IF(J.EQ.-1) THEN CALL FMF019('SUBPAR') XL(1)=-1.0E10 XL(2)=-1.0E10 XV(1)=-1.0E10 XV(2)=-1.0E10 RETURN ENDIF IF(IP.EQ.1) THEN IF(KAS.EQ.1) THEN XL(1)=REAL(HL(1)/PR(1)) XL(2)=REAL(HL(2)/PR(6)) IF(KSTAN.EQ.1) THEN XL(1)=XL(1)+REAL(CP(9)*1.0D3) XL(2)=XL(2)+REAL(CP(19)*1.0D3) ENDIF ELSE XL(1)=REAL(HL(1)) XL(2)=REAL(HL(2)) IF(KSTAN.EQ.1) THEN XL(1)=XL(1)+REAL(PR(1)*CP(9)*1.0D3) XL(2)=XL(2)+REAL(PR(6)*CP(19)*1.0D3) ENDIF ENDIF XV(1)=-1.0E20 XV(2)=-1.0E20 C ELSEIF(IP.EQ.2) THEN IF(KAS.EQ.1) THEN XL(1)=REAL(HL(1)/PR(1)) XL(2)=REAL(HL(2)/PR(6)) XV(1)=REAL(HV(1)/PR(1)) XV(2)=REAL(HV(2)/PR(6)) IF(KSTAN.EQ.1) THEN XL(1)=XL(1)+REAL(CP(9)*1.0D3) XL(2)=XL(2)+REAL(CP(19)*1.0D3) XV(1)=XV(1)+REAL(CP(9)*1.0D3) XV(2)=XV(2)+REAL(CP(19)*1.0D3) ENDIF ELSE XL(1)=REAL(HL(1)) XL(2)=REAL(HL(2)) XV(1)=REAL(HV(1)) XV(2)=REAL(HV(2)) IF(KSTAN.EQ.1) THEN XL(1)=XL(1)+REAL(PR(1)*CP(9)*1.0D3) XL(2)=XL(2)+REAL(PR(6)*CP(19)*1.0D3) XV(1)=XV(1)+REAL(PR(1)*CP(9)*1.0D3) XV(2)=XV(2)+REAL(PR(6)*CP(19)*1.0D3) ENDIF ENDIF C ELSE IF(KAS.EQ.1) THEN XV(1)=REAL(HV(1)/PR(1)) XV(2)=REAL(HV(2)/PR(6)) IF(KSTAN.EQ.1) THEN XV(1)=XV(1)+REAL(CP(9)*1.0D3) XV(2)=XV(2)+REAL(CP(19)*1.0D3) ENDIF ELSE XV(1)=REAL(HV(1)) XV(2)=REAL(HV(2)) IF(KSTAN.EQ.1) THEN XV(1)=XV(1)+REAL(PR(1)*CP(9)*1.0D3) XV(2)=XV(2)+REAL(PR(6)*CP(19)*1.0D3) ENDIF ENDIF XL(1)=-1.0E20 XL(2)=-1.0E20 ENDIF RETURN ENDIF C IF(A.EQ.'S') THEN YY=ZZ CALL FMF039(IP,J,TT,PP,XX,YY,SL,SV) IF(J.EQ.-1) THEN CALL FMF019('SUBPAR') XL(1)=-1.0E10 XL(2)=-1.0E10 XV(1)=-1.0E10 XV(2)=-1.0E10 RETURN ENDIF IF(IP.EQ.1) THEN IF(KAS.EQ.1) THEN XL(1)=REAL(SL(1)/PR(1)) XL(2)=REAL(SL(2)/PR(6)) IF(KSTAN.EQ.1) THEN XL(1)=XL(1)+REAL(CP(10)*1.0D3) XL(2)=XL(2)+REAL(CP(20)*1.0D3) ENDIF ELSE XL(1)=REAL(SL(1)) XL(2)=REAL(SL(2)) IF(KSTAN.EQ.1) THEN XL(1)=XL(1)+REAL(PR(1)*CP(10)*1.0D3) XL(2)=XL(2)+REAL(PR(6)*CP(20)*1.0D3) ENDIF ENDIF XV(1)=-1.0E20 XV(2)=-1.0E20 C ELSEIF(IP.EQ.2) THEN IF(KAS.EQ.1) THEN XL(1)=REAL(SL(1)/PR(1)) XL(2)=REAL(SL(2)/PR(6)) XV(1)=REAL(SV(1)/PR(1)) XV(2)=REAL(SV(2)/PR(6)) IF(KSTAN.EQ.1) THEN XL(1)=XL(1)+REAL(CP(10)*1.0D3) XL(2)=XL(2)+REAL(CP(20)*1.0D3) XV(1)=XV(1)+REAL(CP(10)*1.0D3) XV(2)=XV(2)+REAL(CP(20)*1.0D3) ENDIF ELSE XL(1)=REAL(SL(1)) XL(2)=REAL(SL(2)) XV(1)=REAL(SV(1)) XV(2)=REAL(SV(2)) IF(KSTAN.EQ.1) THEN XL(1)=XL(1)+REAL(PR(1)*CP(10)*1.0D3) XL(2)=XL(2)+REAL(PR(6)*CP(20)*1.0D3) XV(1)=XV(1)+REAL(PR(1)*CP(10)*1.0D3) XV(2)=XV(2)+REAL(PR(6)*CP(20)*1.0D3) ENDIF ENDIF C ELSE IF(KAS.EQ.1) THEN XV(1)=REAL(SV(1)/PR(1)) XV(2)=REAL(SV(2)/PR(6)) IF(KSTAN.EQ.1) THEN XV(1)=XV(1)+REAL(CP(10)*1.0D3) XV(2)=XV(2)+REAL(CP(20)*1.0D3) ENDIF ELSE XV(1)=REAL(SV(1)) XV(2)=REAL(SV(2)) IF(KSTAN.EQ.1) THEN XV(1)=XV(1)+REAL(PR(1)*CP(10)*1.0D3) XV(2)=XV(2)+REAL(PR(6)*CP(20)*1.0D3) ENDIF ENDIF XL(1)=-1.0E20 XL(2)=-1.0E20 ENDIF RETURN ENDIF END C C********************************************************* C C------------------------------------------------- C SUBPB CALCULATES PROPERTIES AT BUBLE POINT. C------------------------------------------------- SUBROUTINE SUBPB(J,T,P,Z,V,H,S) REAL T,P,Z,V,H,S,FMF055,FMF056 DOUBLE PRECISION FMF022,TT,PP,ZZ,VV,HH,SS,LIM,TC,PC,VC $ ,FMF026,FMF072 INTEGER J,JJ C IF(Z.LT.0.0.OR.Z.GT.1.0) THEN CALL FMF020('SUBPB ') P=-1.0E20 V=-1.0E20 H=-1.0E20 S=-1.0E20 RETURN ENDIF C TT=FMF022(1,T) ZZ=FMF022(3,Z) C LIM=1.0D-7 IF(ZZ.LT.LIM) THEN ZZ=0.0 CALL FMF059(2,J,TT,PP) ELSEIF(ZZ.GT.1.0-LIM) THEN ZZ=1.0 CALL FMF059(1,J,TT,PP) ELSE C--------------------------------------------------------- C modified by Akasaka at June 3, 1998. C CALL FMF074(1,JJ,ZZ,TC,PC,VC) C IF(JJ.EQ.0) THEN C IF(TT.LT.TC) THEN PP=FMF026(TT,ZZ) IF(PP.GT.0.0) THEN J=0 ELSE J=-1 ENDIF C ELSE C J=-2 C ENDIF C ELSE C J=-1 C ENDIF ENDIF C IF(J.EQ.0) THEN VV=FMF072('V',1,TT,PP,ZZ) HH=FMF072('H',1,TT,PP,ZZ) SS=FMF072('S',1,TT,PP,ZZ) P=FMF055(2,PP) V=FMF056(1,VV,ZZ) H=FMF056(2,HH,ZZ) S=FMF056(3,SS,ZZ) ELSEIF(J.EQ.-1) THEN CALL FMF019('SUBPB ') P=-1.0E10 V=-1.0E10 H=-1.0E10 S=-1.0E10 ELSE CALL FMF020('SUBPB ') P=-1.0E20 V=-1.0E20 H=-1.0E20 S=-1.0E20 ENDIF RETURN END C C********************************************************* C C---------------------------------------------- C SUBPD CALCULATES PROPERTIES AT DEW POINT. C---------------------------------------------- SUBROUTINE SUBPD(J,T,P,Z,V,H,S) REAL T,P,Z,V,H,S,FMF055,FMF056 DOUBLE PRECISION FMF022,TT,PP,ZZ,VV,HH,SS,LIM,TC,PC,VC $ ,FMF031,FMF072 INTEGER J,JJ C IF(Z.LT.0.0.OR.Z.GT.1.0) THEN CALL FMF020('SUBPD ') P=-1.0E20 V=-1.0E20 H=-1.0E20 S=-1.0E20 RETURN ENDIF C TT=FMF022(1,T) ZZ=FMF022(3,Z) C LIM=1.0D-7 IF(ZZ.LT.LIM) THEN ZZ=0.0 CALL FMF059(2,J,TT,PP) ELSEIF(ZZ.GT.1.0-LIM) THEN ZZ=1.0 CALL FMF059(1,J,TT,PP) ELSE C--------------------------------------------------------- C modified by Akasaka at June 3, 1998. C CALL FMF074(1,JJ,ZZ,TC,PC,VC) C IF(JJ.EQ.0) THEN C IF(TT.LT.TC) THEN PP=FMF031(TT,ZZ) IF(PP.GT.0.0) THEN J=0 ELSE J=-1 ENDIF C ELSE C J=-2 C ENDIF C ELSE C J=-1 C ENDIF ENDIF C IF(J.EQ.0) THEN VV=FMF072('V',3,TT,PP,ZZ) HH=FMF072('H',3,TT,PP,ZZ) SS=FMF072('S',3,TT,PP,ZZ) P=FMF055(2,PP) V=FMF056(1,VV,ZZ) H=FMF056(2,HH,ZZ) S=FMF056(3,SS,ZZ) ELSEIF(J.EQ.-1) THEN CALL FMF019('SUBPD ') P=-1.0E10 V=-1.0E10 H=-1.0E10 S=-1.0E10 ELSE CALL FMF020('SUBPD ') P=-1.0E20 V=-1.0E20 H=-1.0E20 S=-1.0E20 ENDIF RETURN END C C********************************************************* C C--------------------------------------------------------------------- C SUBPST CALCULATES PROPERTIES AT VAPOR-LIQUID EQUILIBRIUM STATE C OF PURE SUBSTACES. C--------------------------------------------------------------------- SUBROUTINE SUBPST(I,J,T,P,VL,VV,HL,HV,SL,SV) DOUBLE PRECISION Z,TT,FMF022,PP,VVL,VVV,HHL,HHV,SSL,SSV REAL T,P,VL,VV,HL,HV,SL,SV,FMF055,FMF056,ER INTEGER I,J C IF(I.EQ.1) THEN Z=1.0 ELSEIF(I.EQ.2) THEN Z=0.0 ELSE CALL FMF020('SUBPST') ER=-1.0E20 P=ER VL=ER VV=ER HL=ER HV=ER SL=ER SV=ER RETURN ENDIF C TT=FMF022(1,T) CALL FMF062(I,J,TT,PP,VVL,VVV,HHL,HHV,SSL,SSV) C IF(J.EQ.0) THEN P=FMF055(2,PP) VL=FMF056(1,VVL,Z) VV=FMF056(1,VVV,Z) HL=FMF056(2,HHL,Z) HV=FMF056(2,HHV,Z) SL=FMF056(3,SSL,Z) SV=FMF056(3,SSV,Z) ELSEIF(J.EQ.-1) THEN CALL FMF019('SUBPST') ER=-1.0E10 P=ER VL=ER VV=ER HL=ER HV=ER SL=ER SV=ER ELSEIF(J.EQ.-2) THEN CALL FMF020('SUBPST') ER=-1.0E20 P=ER VL=ER VV=ER HL=ER HV=ER SL=ER SV=ER ENDIF RETURN END C C********************************************************* C C------------------------------------------------------ C SUBPUR CALCULATES PROPERTIES OF PURE SUBSTACES. C------------------------------------------------------ SUBROUTINE SUBPUR(I,J,T,P,V,H,S) DOUBLE PRECISION Z,TT,PP,VV,HH,SS,FMF022 REAL T,P,V,H,S,FMF056,ER INTEGER I,J,COMBI DOUBLE PRECISION PR(1:11),CP(1:20) COMMON/FMFC/PR,CP,COMBI C IF(I.EQ.1) THEN Z=1.0 ELSEIF(I.EQ.2) THEN Z=0.0 ELSE CALL FMF020('SUBPUR') ER=-1.0E20 V=ER H=ER S=ER RETURN ENDIF C TT=FMF022(1,T) PP=FMF022(2,P) CALL FMF063(I,J,TT,PP,VV,HH,SS) C IF(J.EQ.0) THEN V=FMF056(1,VV,Z) H=FMF056(2,HH,Z) S=FMF056(3,SS,Z) ELSE CALL FMF019('SUBPUR') ER=-1.0E10 V=ER H=ER S=ER ENDIF RETURN END C C******************************************************* C C-------------------------------------------------- C SUBTB CALCULATES PROPERTIES AT BUBLE POINT. C-------------------------------------------------- SUBROUTINE SUBTB(J,T,P,Z,V,H,S) REAL T,P,Z,V,H,S,FMF055,FMF056 DOUBLE PRECISION TT,PP,ZZ,VV,HH,SS,LIM,FMF022 $ ,PC1,FMF029,FMF072 DOUBLE PRECISION PR(1:11),CP(1:20) INTEGER COMBI COMMON/FMFC/PR,CP,COMBI C IF(Z.LT.0.0.OR.Z.GT.1.0) THEN CALL FMF020('SUBTB ') T=-1.0E20 V=-1.0E20 H=-1.0E20 S=-1.0E20 RETURN ENDIF C PP=FMF022(2,P) ZZ=FMF022(3,Z) C LIM=1.0D-7 IF(ZZ.LT.LIM) THEN ZZ=0.0 CALL FMF061(2,J,TT,PP) ELSEIF(ZZ.GT.1.0-LIM) THEN ZZ=1.0 CALL FMF061(1,J,TT,PP) ELSE PC1=PR(3) IF(PP.LT.PC1) THEN TT=FMF029(PP,ZZ) IF(TT.GT.0.0) THEN J=0 ELSE J=-1 ENDIF ELSE J=-2 ENDIF ENDIF C IF(J.EQ.0) THEN VV=FMF072('V',1,TT,PP,ZZ) HH=FMF072('H',1,TT,PP,ZZ) SS=FMF072('S',1,TT,PP,ZZ) T=FMF055(1,TT) V=FMF056(1,VV,ZZ) H=FMF056(2,HH,ZZ) S=FMF056(3,SS,ZZ) ELSEIF(J.EQ.-1) THEN CALL FMF019('SUBTB ') T=-1.0E10 V=-1.0E10 H=-1.0E10 S=-1.0E10 ELSE CALL FMF020('SUBTB ') T=-1.0E20 V=-1.0E20 H=-1.0E20 S=-1.0E20 ENDIF RETURN END C C******************************************************** C C--------------------------------------------------- C SUBTD CALCULATES PROPERTIES AT BUBLE POINT. C--------------------------------------------------- SUBROUTINE SUBTD(J,T,P,Z,V,H,S) REAL T,P,Z,V,H,S,FMF055,FMF056 DOUBLE PRECISION TT,PP,ZZ,VV,HH,SS,LIM,FMF022 $ ,PC1,FMF034,FMF072 DOUBLE PRECISION PR(1:11),CP(1:20) INTEGER COMBI COMMON/FMFC/PR,CP,COMBI C IF(Z.LT.0.0.OR.Z.GT.1.0) THEN CALL FMF020('SUBTD ') T=-1.0E20 V=-1.0E20 H=-1.0E20 S=-1.0E20 RETURN ENDIF C PP=FMF022(2,P) ZZ=FMF022(3,Z) C LIM=1.0D-7 IF(ZZ.LT.LIM) THEN ZZ=0.0 CALL FMF061(2,J,TT,PP) ELSEIF(ZZ.GT.1.0-LIM) THEN ZZ=1.0 CALL FMF061(1,J,TT,PP) ELSE PC1=PR(3) IF(PP.LT.PC1) THEN TT=FMF034(PP,ZZ) IF(TT.GT.0.0) THEN J=0 ELSE J=-1 ENDIF ELSE J=-2 ENDIF ENDIF C IF(J.EQ.0) THEN VV=FMF072('V',3,TT,PP,ZZ) HH=FMF072('H',3,TT,PP,ZZ) SS=FMF072('S',3,TT,PP,ZZ) T=FMF055(1,TT) V=FMF056(1,VV,ZZ) H=FMF056(2,HH,ZZ) S=FMF056(3,SS,ZZ) ELSEIF(J.EQ.-1) THEN CALL FMF019('SUBTD ') T=-1.0E10 V=-1.0E10 H=-1.0E10 S=-1.0E10 ELSE CALL FMF020('SUBTD ') T=-1.0E20 V=-1.0E20 H=-1.0E20 S=-1.0E20 ENDIF RETURN END C C************************************************************* C C------------------------------------------------------------------- C SUBXY CALCULATES PROPERTIES AT VAPOR-LIQUID EQUILIBRIUM STATE C OF MIXTURE. C------------------------------------------------------------------- SUBROUTINE SUBXY(J,T,P,X,Y,VL,VV,HL,HV,SL,SV) REAL T,P,X,Y,VL,VV,HL,HV,SL,SV,FMF055,FMF056,ER DOUBLE PRECISION TT,PP,FMF022,P1,P2,ZC,PC,VC $ ,XX,YY,VVL,VVV, & HHL,HHV,SSL,SSV INTEGER J1,J2,J,COMBI DOUBLE PRECISION PR(1:11),CP(1:20) COMMON/FMFC/PR,CP,COMBI C TT=FMF022(1,T) PP=FMF022(2,P) C TC1=PR(2) TC2=PR(7) IF(TT.GT.TC2) J=-2 IF(TT.LE.TC2) THEN IF(TT.GE.TC1) THEN C CALL FMF074(2,J1,ZC,TT,PC,VC) C P1=PC P1=PR(3) CALL FMF059(2,J2,TT,P2) ELSEIF(TT.LT.TC1) THEN CALL FMF059(1,J1,TT,P1) CALL FMF059(2,J2,TT,P2) ENDIF IF(PP.GT.P1.OR.PP.LT.P2) J=-2 ENDIF C IF(J.NE.2) THEN CALL FMF046(J,TT,PP,XX,YY,VVL,VVV,HHL,HHV,SSL,SSV) ENDIF C IF(J.EQ.0) THEN X=FMF055(3,XX) Y=FMF055(3,YY) VL=FMF056(1,VVL,XX) VV=FMF056(1,VVV,YY) HL=FMF056(2,HHL,XX) HV=FMF056(2,HHV,YY) SL=FMF056(3,SSL,XX) SV=FMF056(3,SSV,YY) ELSEIF(J.EQ.-1) THEN CALL FMF019('SUBXY ') ER=-1.0E10 X=ER Y=ER VL=ER VV=ER HL=ER HV=ER SL=ER SV=ER ELSE CALL FMF020('SUBXY ') ER=-1.0E20 X=ER Y=ER VL=ER VV=ER HL=ER HV=ER SL=ER SV=ER ENDIF RETURN END C C********************************************************* C C------------------------------------------------ C TBP CALCULATES TEMPERATURE AT BUBLE POINT. C------------------------------------------------ REAL FUNCTION TBP(P,Z) REAL P,Z,FMF055 DOUBLE PRECISION TT,PP,ZZ,LIM,FMF022,PCMAX,FMF029 DOUBLE PRECISION PR(1:11),CP(1:20) INTEGER COMBI COMMON/FMFC/PR,CP,COMBI C IF(Z.LT.0.0.OR.Z.GT.1.0) THEN CALL FMF020('TBP ') TBP=-1.0E20 RETURN ENDIF C PP=FMF022(2,P) ZZ=FMF022(3,Z) C LIM=1.0D-7 IF(ZZ.LT.LIM) THEN ZZ=0.0 CALL FMF061(2,J,TT,PP) ELSEIF(ZZ.GT.1.0-LIM) THEN ZZ=1.0 CALL FMF061(1,J,TT,PP) ELSE IF(PR(3).GT.PR(8)) PCMAX=PR(3) IF(PR(3).LT.PR(8)) PCMAX=PR(8) IF(PP.LT.PCMAX) THEN TT=FMF029(PP,ZZ) IF(TT.GT.0.0) THEN J=0 ELSE J=-1 ENDIF ELSE J=-2 ENDIF ENDIF C IF(J.EQ.0) THEN TBP=FMF055(1,TT) ELSEIF(J.EQ.-1) THEN CALL FMF019('TBP ') TBP=-1.0E10 ELSE CALL FMF020('TBP ') TBP=-1.0E20 ENDIF RETURN END C C********************************************************* C C---------------------------------------------- C TDP CALCULATES TEMPERATURE AT DEW POINT. C---------------------------------------------- REAL FUNCTION TDP(P,Z) REAL P,Z,FMF055 DOUBLE PRECISION TT,PP,ZZ,LIM,FMF022,PCMAX,FMF034 DOUBLE PRECISION PR(1:11),CP(1:20) INTEGER COMBI COMMON/FMFC/PR,CP,COMBI C IF(Z.LT.0.0.OR.Z.GT.1.0) THEN CALL FMF020('TDP ') TDP=-1.0E20 RETURN ENDIF C PP=FMF022(2,P) ZZ=FMF022(3,Z) C LIM=1.0D-7 IF(ZZ.LT.LIM) THEN ZZ=0.0 CALL FMF061(2,J,TT,PP) ELSEIF(ZZ.GT.1.0-LIM) THEN ZZ=1.0 CALL FMF061(1,J,TT,PP) ELSE IF(PR(3).GT.PR(8)) PCMAX=PR(3) IF(PR(3).LT.PR(8)) PCMAX=PR(8) IF(PP.LT.PCMAX) THEN TT=FMF034(PP,ZZ) IF(TT.GT.0.0) THEN J=0 ELSE J=-1 ENDIF ELSE J=-2 ENDIF ENDIF C IF(J.EQ.0) THEN TDP=FMF055(1,TT) ELSEIF(J.EQ.-1) THEN CALL FMF019('TDP ') TDP=-1.0E10 ELSE CALL FMF020('TDP ') TDP=-1.0E20 ENDIF RETURN END C C********************************************************** C C------------------------------------------------------------------ C TSPM CALCULATES SATURATED VAPOR TEMPERATURE OF PURE SUBSTANCE. C------------------------------------------------------------------ REAL FUNCTION TSPM(I,P) DOUBLE PRECISION TT,PP,FMF022 REAL P,FMF055 INTEGER J C IF(I.EQ.1.OR.I.EQ.2) THEN PP=FMF022(2,P) CALL FMF061(I,J,TT,PP) IF(J.EQ.0) THEN TSPM=FMF055(1,TT) ELSEIF(J.EQ.-1) THEN CALL FMF019('TSPM ') TSPM=-1.0E10 ELSEIF(J.EQ.-2) THEN CALL FMF020('TSPM ') TSPM=-1.0E20 ENDIF RETURN ELSE CALL FMF020('TSPM ') TSPM=-1.0E20 RETURN ENDIF END C C********************************************************** C