c ---------------------------------------------------------- c subroutine FMF000 c set characteristic properties of 3 components mixture c combi : combination of components c types(n) : the form of equation of isobaric heat c capacity of ideal gas c 1 : 1st component c 2 : 2nd component c 3 : 3rd component c =1 : c_p=a+bT+cT^2+dT^3 c =2 : c_p=a+b{(c/T)sinh(c/T)}^2 c +d{(e/T)cosh(e/T)}^2 c 4 : Number of components c const(n) : fundamental constants c 1 : molecular weight [kg/kmol] c 2 : critical temperature [K] c 3 : critical pressure [Pa] c 4 : critical volume [m^3/kmol] c 5 : acentric factor [-] c 6-9: Coefficients of Saturated curve c P_{sat}=A0+A1*T+A2*T^2+A3*T^3 c 11 : Minimun temperature [K] c 12 : Maximum temperature [K] c 13 : Enthalpy difference between c Standard states [J/kmol] (no need) c 14 : Entropy difference between c Standard states [J/(kmol K)] (no need) c 21- : 2nd component c 41- : 3rd component c 61 : Gas constant R [kJ/(kmol K)] c 62 : Interaction parameter f11 c 63 : Interaction parmaeter f12 c 64 : Interaction parmaeter f13 c 65 : Interaction parmaeter f21 c 66 : Interaction parmaeter f22 c 67 : Interaction parmaeter f23 c 68 : Interaction parmaeter f31 c 69 : Interaction parmaeter f32 c 70 : Interaction parmaeter f33 c coeff(n) : the coefficients of CSD c equation of state c 1-3 : the coefficients of CSD eq. of st. c ai [kJ m^3/kmol^2] c 4-6 : the coefficients of CSD eq. of st. c bi [m^3/kmol] c 7- : the coefficients of isobaric heat c capacity of ideal gas [kJ/(kmol K)] c 21- : the coefficients of 2nd component c 41- : the coefficients of 3rd component c NANES(n) : name of component, etc. c 1 : name of 1st component c 2 : structural formula of 1st component c 3 : chemical abstracts name of 1st component c 4 : name of 2nd component c 5 : structural formula of 2nd component c 6 : chemical abstracts name of 2nd component c 7 : name of 3rd component c 8 : structural formula of 3rd component c 9 : chemical abstracts name of 3rd component c 10 : version number c ---------------------------------------------------------- SUBROUTINE FMF000(J,COMBI,CONST,COEFF,NAMES,TYPES) DOUBLE PRECISION CONST(1:80),COEFF(1:60) $ ,PR32(1:20),PR134A(1:20),PR22(1:20),PR123(1:20) $ ,CR32(1:20),CR134A(1:20),CR22(1:20),CR123(1:20) $ ,PR11(1:20),PR12(1:20),PR13(1:20),PR13B1(1:20) $ ,CR11(1:20),CR12(1:20),CR13(1:20),CR13B1(1:20) $ ,PR23(1:20),PR114(1:20),PR152A(1:20) $ ,CR23(1:20),CR114(1:20),CR152A(1:20) $ ,KIJ1(1:3,1:3),KIJ2(1:3,1:3),KIJ3(1:3,1:3) $ ,KIJ4(1:3,1:3),KIJ5(1:3,1:3),KIJ6(1:3,1:3) $ ,KIJ7(1:3,1:3),KIJ8(1:3,1:3),KIJ9(1:3,1:3) $ ,KIJ10(1:3,1:3) c $ ,PR14(1:20),PR113(1:20),PR142B(1:20) c $ ,CR14(1:20),CR113(1:20),CR142B(1:20) INTEGER I,J,K,COMBI,TYPES(1:10) $ ,TR32,TR134A,TR22,TR123,TR11,TR12,TR13,TR13B1,TR23 $ ,TR114,TR152A c $ ,TR14,TR113,TR142B CHARACTER*40 NAMES(1:10) $ ,NR32(1:3),NR134A(1:3),NR22(1:3),NR123(1:3) $ ,NR11(1:3),NR12(1:3),NR13(1:3),NR13B1(1:3) $ ,NR23(1:3),NR114(1:3),NR152A(1:3) c $ ,NR14(1:3),NR113(1:3),NR142B(1:3) c c ----- R32 ----- c fundamental constants from ??? DATA (PR32(I),I=1,20) /20*0.0D0/ c coefficients of CSD eq. of st. from ??? DATA(CR32(I),I=1,20)/20*0.0D0/ c form of Cp TR32=1 c name of component, etc. DATA(NR32(I),I=1,3)/'R32','NUL' $ ,'NUL'/ c c ----- R134A ----- c fundamental constants from ??? DATA(PR134A(I),I=1,20)/20*0.0D0/ c coefficients of CSD eq. of st. from ??? DATA(CR134A(I),I=1,20)/20*0.0D0/ c form of Cp TR134A=1 c name of component, etc. DATA(NR134A(I),I=1,3)/'R134A','NUL' $ ,'NUL'/ c c ----- R22 ----- c fundamental constants from NIST DATA(PR22(I),I=1,20)/86.468D0,369.30D0,4.9710D6 $ ,0.168615D0,0.2192D0 c $ ,10229.9963858457,-98.09555191286287,0.2281508492856921D0 c $ ,-2.707161D-4,1.186741D-6 $ ,-7.658857D3,1.138805D2,-5.668891D-1,9.487912D-4,0.0D0 $ ,1.65D2,4.0D2 c $ ,0.4044233305507963D5,0.1351617511619961D3 $ ,0.2314873305507963D+05,0.4869375116199605D+02 $ ,6*0.0D0/ c coefficients of CSD eq. of st. from NIST DATA(CR22(I),I=1,20)/2624.6155D0,-2.6730483D-3,-1.3323785D-6 $ ,0.11739546D0,-1.4027169D-4,-0.52160854D-7 $ ,21.98390D0,0.1277439D0,-4.788723D-5,11*0.0D0/ c form of Cp TR22=1 c name of component, etc. DATA(NR22(I),I=1,3)/'R22','CHCLF2' $ ,'CHLORODIFLUOROMETHANE'/ c c ----- R123 ----- c fundamental constants from NIST DATA(PR123(I),I=1,20)/152.931D0,456.86D0,3.666D6 $ ,0.275276D0,0.2816D0 c $ ,8235.918965282606,-62.54124991662765,0.115400051937141D0 c $ ,-6.967085D-5,3.777315D-7 $ ,-6.714577D3,7.870805D1,-3.085707D-1,4.061361D-4,0.0D0 $ ,2.1D2,5.0D2 c $ ,.6478951790857697D5,.2232364862142001D3 $ ,0.3420331790857697D+05,0.7030548621420012D+02 $ ,6*0.0D0/ c coefficients of CSD eq. of st. from NIST DATA(CR123(I),I=1,20)/6033.2949D0,-2.3789106D-3,-0.84728055D-6 $ ,0.19954886D0,-1.8949345D-4,-0.67680325D-7 $ ,17.01154D0,0.404631D0,-4.644803D-4,2.347418D-7 $ ,10*0.0D0/ c form of Cp TR123=1 c name of component, etc. DATA(NR123(I),I=1,3)/'R123','CHCL2-CF3' $ ,'1,1-DICHLORO-2,2,2-TRIFLUOROETHANE'/ c c ----- R11 ----- c fundamental constants from D.D. DATA(PR11(I),I=1,20)/137.368D0,471.20D0,4.4076D6,0.24800D0 $ ,0.279D0 $ ,-5.517673D3,6.641110D1,-2.672594D-1,3.612618D-4,0.0D0 $ ,2.0D2,4.0D2 $ ,0.3172783545539200D+05,0.6210547583629349D+02 $ ,6*0.0D0/ c coefficients of CSD eq. of st. from ASHRAE and JANAF DATA(CR11(I),I=1,20)/4971.54D0,-2.24669D-3,-0.51194D-6 $ ,0.176659D0,-1.74531D-4,-3.49717D-8 $ ,22.0418D0,0.260895D0,-2.45319D-4,11*0.0D0/ c form of Cp TR11=1 c name of component, etc. DATA(NR11(I),I=1,3)/'R11','CCL3F' $ ,'TRICHLOROFLUOROMETHANE'/ c c ----- R12 ----- c fundamental constants from D.D. DATA(PR12(I),I=1,20)/120.913D0,384.95D0,4.1249D6,0.21700 $ ,0.1796 $ ,-6.236706D3,8.687112D1,-4.098606D-1,6.572207D-4,0.0D0 $ ,2.0D2,4.0D2 $ ,0.2397776768074632D+05,0.4986764676386134D+02 $ ,6*0.0D0/ c coefficients of CSD eq. of st. from ??? DATA(CR12(I),I=1,20)/3524.12D0,-2.77230D-3,-0.67318D-6 $ ,0.153755D0,-1.84195D-4,-5.03644D-8 $ ,17.5387D0,0.248546D0,-2.16271D-4,11*0.0D0/ c form of Cp TR12=1 c name of component, etc. DATA(NR12(I),I=1,3)/'R12','CCL2F2' $ ,'DICHLORODIFLUOROMETHANE'/ c c ----- R13 ----- c fundamental constants from D.D. DATA(PR13(I),I=1,20)/104.459D0,301.96D0,3.9460D6,0.18028 $ ,0.1800D0 $ ,-6.796813D3,1.162571D2,-6.793426D-1,1.359273D-3,0.0D0 $ ,2.0D2,4.0D2 $ ,0.1595598499707058D+05,0.3496337457993984D+02 $ ,6*0.0D0/ c coefficients of CSD eq. of st. from ??? DATA(CR13(I),I=1,20)/2298.13D0,-3.41828D-3,-1.52430D-6 $ ,0.128141D0,-1.84474D-4,-10.7951D-8 $ ,13.9300D0,0.232181D0,-1.82929D-4,11*0.0D0/ c form of Cp TR13=1 c name of component, etc. DATA(NR13(I),I=1,3)/'R13','CCLF3' $ ,'CHLOROTRIFLUOROMETHANE'/ c c ----- R13B1 ----- c fundamental constants from D.D. DATA(PR13B1(I),I=1,20)/148.910D0,340.15D0,3.9719D6,0.20000 $ ,0.1727 $ ,-7.187330D3,1.086872D2,-5.605256D-1,9.890945D-4,0.0D0 $ ,2.0D2,4.0D2 $ ,0.1985685773294219D+05,0.4250336857062973D+02 $ ,6*0.0D0/ c coefficients of CSD eq. of st. from ??? DATA(CR13B1(I),I=1,20)/2728.10D0,-2.79791D-3,-1.50848D-6 $ ,0.139949D0,-1.82428D-4,-7.75898D-8 $ ,19.9537D0,0.216394D0,-1.70241D-4,11*0.0D0/ c form of Cp TR13B1=1 c name of component, etc. DATA(NR13B1(I),I=1,3)/'R13B1','CBRF3' $ ,'BROMOTRIFLUOROMETHANE'/ c c ----- R14 ----- c fundamental constants from D.D. c DATA(PR14(I),I=1,20)/88.005D0,227.50D0,3.7389D6,0.14000D0 c $ ,0.1855D0,5*0.0D0,2.0D2,4.0D2,8*0.0D0/ c coefficients of CSD eq. of st. from ??? c DATA(CR14(I),I=1,20)/1393.60D0,-4.81985D-3,-1.89167D-3 c $ ,0.100601D0,-1.94974D-4,-13.5408D-8 c $ ,11.0629D0,0.209740D0,-1.40992D-4,11*0.0D0/ c form of Cp c TR14=1 c name of component, etc. c DATA(NR14(I),I=1,3)/'R14','CF4' c $ ,'CARBON TETRAFLUORIDE'/ c c ----- R23 ----- c fundamental constants from D.D. DATA(PR23(I),I=1,20)/70.014D0,298.89D0,4.8362D6,0.13330D0 $ ,0.2672D0 $ ,-1.544310D4,2.430819D2,-1.307795D0,2.2412726D-3,0.0D0 $ ,2.0D2,4.0D2 $ ,0.1601776059390253D+05,0.3484407190490592D+02 $ ,6*0.0D0/ c coefficients of CSD eq. of st. from ??? DATA(CR23(I),I=1,20)/2025.93D0,-4.68206D-3,0.99552D-6 $ ,0.103137D0,-2.29653D-4,15.5760D-8 $ ,20.4760D0,0.106183D0,-0.12189D-4,11*0.0D0/ c form of Cp TR23=1 c name of component, etc. DATA(NR23(I),I=1,3)/'R23','CHF3' $ ,'TRIFLUOROMETHANE'/ c c ----- R113 ----- c fundamental constants from D.D. c DATA(PR113(I),I=1,20)/187.375D0,487.25D0,3.4147D6,0.32530D0 c $ ,0.2552D0,5*0.0D0,2.4D2,4.2D2,8*0.0D0/ c coefficients of CSD eq. of st. from ??? c DATA(CR113(I),I=1,20)/7489.98D0,-2.34255D-3,-0.50363D-6 c $ ,0.233100D0,-2.03232D-4,-8.2063D-8 c $ ,76.2637D0,0.119641D0,0.71879D-4,11*0.0D0/ c form of Cp c TR113=1 c name of component, etc. c DATA(NR113(I),I=1,3)/'R113','C2CL3F3' c $ ,'1,1,2-TRICHLOROTRIFLUOROETHANE'/ c c ----- R114 ----- c fundamental constants from D.D. DATA(PR114(I),I=1,20)/170.921D0,418.85D0,3.2600D6,0.29400D0 $ ,0.2520 $ ,-8.571122D3,1.005603D2,-4.004197D-1,5.432180D-4,0.0D0 $ ,2.4D2,4.7D2 $ ,0.3144909457642187D+05,0.6922841029102383D+02 $ ,6*0.0D0/ c coefficients of CSD eq. of st. from ??? DATA(CR114(I),I=1,20)/9771.35D0,-5.85557D-3,3.99413D-6 $ ,0.306318D0,-7.96444D-4,78.1059D-8 $ ,20.7005D0,0.464035D0,-4.17589D-4,11*0.0D0/ c form of Cp TR114=1 c name of component, etc. DATA(NR114(I),I=1,3)/'R114','C2CL2F4' $ ,'1,2-DICHLOROTETRAFLUOROETHANE'/ c c ----- R142b ----- c fundamental constants from D.D. c DATA(PR142B(I),I=1,20)/100.495D0,410.20D0,4.1239D6,0.23100 c $ ,0.2368D0,5*0.0D0,2.5D2,5.0D2,8*0.0D0/ c coefficients of CSD eq. of st. from ??? c DATA(CR142B(I),I=1,20)/2990.00D0,-0.54056D-3,-4.12642D-6 c $ ,0.146006D0,-0.89250D-4,-18.0562D-8 c $ ,23.7611D0,0.231706D0,-1.06534D-4,11*0.0D0/ c form of Cp c TR142B=1 c name of component, etc. c DATA(NR142B(I),I=1,3)/'R142B','C2H3CLF2' c $ ,'1-CHLORO-1,1-DIFLUOROETHANE'/ c c ----- R152a ----- c fundamental constants from D.D. DATA(PR152A(I),I=1,20)/66.051D0,386.60D0,4.4988D6,0.18100D0 $ ,0.2629D0 $ ,-1.323144D4,1.701294D2,-7.357102D-1,1.075109D-3,0.0D0 $ ,2.0D2,4.0D2 $ ,0.2541295584467842D+05,0.5221990583774483D+02 $ ,6*0.0D0/ c coefficients of CSD eq. of st. from ??? DATA(CR152A(I),I=1,20)/2254.37D0,-0.58778D-3,-4.37432D-6 $ ,0.116521D0,-0.90488D-4,-11.4563D-8 $ ,22.2804D0,0.154009D0,-0.03067D-4,11*0.0D0/ c form of Cp TR152A=1 c name of component, etc. DATA(NR152A(I),I=1,3)/'R152A','C2H4F2' $ ,'1,1-DIFLUOROETHANE'/ c c ------------------------------------------------ c Interaction parameter k_{ij} c ------------------------------------------------ c ------ R22-R123 ------------ DATA(KIJ1(1,I),I=1,3)/3*0.0D0/ DATA(KIJ1(2,I),I=1,3)/3*0.0D0/ DATA(KIJ1(3,I),I=1,3)/3*0.0D0/ c ------ R13B1-R152A --------- DATA(KIJ2(1,I),I=1,3)/0.0D0,0.089D0,0.0D0/ DATA(KIJ2(2,I),I=1,3)/0.089D0,2*0.0D0/ DATA(KIJ2(3,I),I=1,3)/3*0.0D0/ c ------ R22-R12 ------------- DATA(KIJ3(1,I),I=1,3)/0.0D0,0.041D0,0.0D0/ DATA(KIJ3(2,I),I=1,3)/0.041,2*0.0D0/ DATA(KIJ3(3,I),I=1,3)/3*0.0D0/ c ------ R23-R13 ------------- DATA(KIJ4(1,I),I=1,3)/0.0D0,0.089D0,0.0D0/ DATA(KIJ4(2,I),I=1,3)/0.089D0,2*0.0D0/ DATA(KIJ4(3,I),I=1,3)/3*0.0D0/ c ------ R13-R12 ------------- DATA(KIJ5(1,I),I=1,3)/0.0D0,0.035D0,0.0D0/ DATA(KIJ5(2,I),I=1,3)/0.035D0,2*0.0D0/ DATA(KIJ5(3,I),I=1,3)/3*0.0D0/ c ------ R12-R152a ----------- DATA(KIJ6(1,I),I=1,3)/0.0D0,0.035D0,0.0D0/ DATA(KIJ6(2,I),I=1,3)/0.035D0,2*0.0D0/ DATA(KIJ6(3,I),I=1,3)/3*0.0D0/ c ------ R22-R114 ------------ DATA(KIJ7(1,I),I=1,3)/0.0D0,0.03D0,0.0D0/ DATA(KIJ7(2,I),I=1,3)/0.03D0,2*0.0D0/ DATA(KIJ7(3,I),I=1,3)/3*0.0D0/ c ------ R23-R12 ------------- DATA(KIJ8(1,I),I=1,3)/0.0D0,0.088D0,0.0D0/ DATA(KIJ8(2,I),I=1,3)/0.088D0,2*0.0D0/ DATA(KIJ8(3,I),I=1,3)/3*0.0D0/ c ------ R22-R11 ----------- DATA(KIJ9(1,I),I=1,3)/0.0D0,0.037D0,0.0D0/ DATA(KIJ9(2,I),I=1,3)/0.037D0,2*0.0D0/ DATA(KIJ9(3,I),I=1,3)/3*0.0D0/ c ------ R32-R134a ----------- DATA(KIJ10(1,I),I=1,3)/3*0.0D0/ DATA(KIJ10(2,I),I=1,3)/3*0.0D0/ DATA(KIJ10(3,I),I=1,3)/3*0.0D0/ c c ----------------------------------------------- c version number c ----------------------------------------------- NAMES(10)='12.1' c c ----------------------------------------------- c gas constant [kJ/(kg K)] c ----------------------------------------------- CONST(61)=8.314510D0 c c ----------------------------------------------- c set characteristc constants of componets c ----------------------------------------------- DO 10 I=1,19 CONST(61+I)=0.0D0 10 CONTINUE IF (COMBI.EQ.1) THEN DO 40 I=1,20 CONST(I)=PR22(I) CONST(20+I)=PR123(I) CONST(40+I)=0.0D0 COEFF(I)=CR22(I) COEFF(20+I)=CR123(I) COEFF(40+I)=0.0D0 40 CONTINUE DO 50 I=1,3 NAMES(I)=NR22(I) NAMES(3+I)=NR123(I) NAMES(6+I)='NUL' 50 CONTINUE TYPES(1)=TR22 TYPES(2)=TR123 TYPES(3)=0 TYPES(4)=2 DO 70 I=1,6 TYPES(4+I)=0 70 CONTINUE DO 100 I=1,3 DO 110 K=1,3 CONST(58+3*I+K)=KIJ1(I,K) 110 CONTINUE 100 CONTINUE J=0 ELSEIF (COMBI.EQ.2) THEN DO 210 I=1,20 CONST(I)=PR13B1(I) CONST(20+I)=PR152A(I) CONST(40+I)=0.0D0 COEFF(I)=CR13B1(I) COEFF(20+I)=CR152A(I) COEFF(40+I)=0.0D0 210 CONTINUE DO 220 I=1,3 NAMES(I)=NR13B1(I) NAMES(3+I)=NR152A(I) NAMES(6+I)='NUL' 220 CONTINUE TYPES(1)=TR13B1 TYPES(2)=TR152A TYPES(3)=0 TYPES(4)=2 DO 230 I=1,6 TYPES(4+I)=0 230 CONTINUE DO 240 I=1,3 DO 250 K=1,3 CONST(58+3*I+K)=KIJ2(I,K) 250 CONTINUE 240 CONTINUE J=0 ELSEIF (COMBI.EQ.3) THEN DO 310 I=1,20 CONST(I)=PR22(I) CONST(20+I)=PR12(I) CONST(40+I)=0.0D0 COEFF(I)=CR22(I) COEFF(20+I)=CR12(I) COEFF(40+I)=0.0D0 310 CONTINUE DO 320 I=1,3 NAMES(I)=NR22(I) NAMES(3+I)=NR12(I) NAMES(6+I)='NUL' 320 CONTINUE TYPES(1)=TR22 TYPES(2)=TR12 TYPES(3)=0 TYPES(4)=2 DO 330 I=1,6 TYPES(4+I)=0 330 CONTINUE DO 340 I=1,3 DO 350 K=1,3 CONST(58+3*I+K)=KIJ3(I,K) 350 CONTINUE 340 CONTINUE J=0 ELSEIF (COMBI.EQ.4) THEN DO 410 I=1,20 CONST(I)=PR23(I) CONST(20+I)=PR13(I) CONST(40+I)=0.0D0 COEFF(I)=CR23(I) COEFF(20+I)=CR13(I) COEFF(40+I)=0.0D0 410 CONTINUE DO 420 I=1,3 NAMES(I)=NR23(I) NAMES(3+I)=NR13(I) NAMES(6+I)='NUL' 420 CONTINUE TYPES(1)=TR23 TYPES(2)=TR13 TYPES(3)=0 TYPES(4)=2 DO 430 I=1,6 TYPES(4+I)=0 430 CONTINUE DO 440 I=1,3 DO 450 K=1,3 CONST(58+3*I+K)=KIJ4(I,K) 450 CONTINUE 440 CONTINUE J=0 ELSEIF (COMBI.EQ.5) THEN DO 510 I=1,20 CONST(I)=PR13(I) CONST(20+I)=PR12(I) CONST(40+I)=0.0D0 COEFF(I)=CR13(I) COEFF(20+I)=CR12(I) COEFF(40+I)=0.0D0 510 CONTINUE DO 520 I=1,3 NAMES(I)=NR13(I) NAMES(3+I)=NR12(I) NAMES(6+I)='NUL' 520 CONTINUE TYPES(1)=TR13 TYPES(2)=TR12 TYPES(3)=0 TYPES(4)=2 DO 530 I=1,6 TYPES(4+I)=0 530 CONTINUE DO 540 I=1,3 DO 550 K=1,3 CONST(58+3*I+K)=KIJ5(I,K) 550 CONTINUE 540 CONTINUE J=0 ELSEIF (COMBI.EQ.6) THEN DO 610 I=1,20 CONST(I)=PR12(I) CONST(20+I)=PR152A(I) CONST(40+I)=0.0D0 COEFF(I)=CR12(I) COEFF(20+I)=CR152A(I) COEFF(40+I)=0.0D0 610 CONTINUE DO 620 I=1,3 NAMES(I)=NR12(I) NAMES(3+I)=NR152A(I) NAMES(6+I)='NUL' 620 CONTINUE TYPES(1)=TR12 TYPES(2)=TR152A TYPES(3)=0 TYPES(4)=2 DO 630 I=1,6 TYPES(4+I)=0 630 CONTINUE DO 640 I=1,3 DO 650 K=1,3 CONST(58+3*I+K)=KIJ6(I,K) 650 CONTINUE 640 CONTINUE J=0 ELSEIF (COMBI.EQ.7) THEN DO 710 I=1,20 CONST(I)=PR22(I) CONST(20+I)=PR114(I) CONST(40+I)=0.0D0 COEFF(I)=CR22(I) COEFF(20+I)=CR114(I) COEFF(40+I)=0.0D0 710 CONTINUE DO 720 I=1,3 NAMES(I)=NR22(I) NAMES(3+I)=NR114(I) NAMES(6+I)='NUL' 720 CONTINUE TYPES(1)=TR22 TYPES(2)=TR114 TYPES(3)=0 TYPES(4)=2 DO 730 I=1,6 TYPES(4+I)=0 730 CONTINUE DO 740 I=1,3 DO 750 K=1,3 CONST(58+3*I+K)=KIJ7(I,K) 750 CONTINUE 740 CONTINUE J=0 ELSEIF (COMBI.EQ.8) THEN DO 810 I=1,20 CONST(I)=PR23(I) CONST(20+I)=PR12(I) CONST(40+I)=0.0D0 COEFF(I)=CR23(I) COEFF(20+I)=CR12(I) COEFF(40+I)=0.0D0 810 CONTINUE DO 820 I=1,3 NAMES(I)=NR23(I) NAMES(3+I)=NR12(I) NAMES(6+I)='NUL' 820 CONTINUE TYPES(1)=TR23 TYPES(2)=TR12 TYPES(3)=0 TYPES(4)=2 DO 830 I=1,6 TYPES(4+I)=0 830 CONTINUE DO 840 I=1,3 DO 850 K=1,3 CONST(58+3*I+K)=KIJ8(I,K) 850 CONTINUE 840 CONTINUE J=0 ELSEIF (COMBI.EQ.9) THEN DO 910 I=1,20 CONST(I)=PR22(I) CONST(20+I)=PR11(I) CONST(40+I)=0.0D0 COEFF(I)=CR22(I) COEFF(20+I)=CR11(I) COEFF(40+I)=0.0D0 910 CONTINUE DO 920 I=1,3 NAMES(I)=NR22(I) NAMES(3+I)=NR11(I) NAMES(6+I)='NUL' 920 CONTINUE TYPES(1)=TR22 TYPES(2)=TR11 TYPES(3)=0 TYPES(4)=2 DO 930 I=1,6 TYPES(4+I)=0 930 CONTINUE DO 940 I=1,3 DO 950 K=1,3 CONST(58+3*I+K)=KIJ9(I,K) 950 CONTINUE 940 CONTINUE J=0 ELSEIF (COMBI.EQ.10) THEN DO 20 I=1,20 CONST(I)=PR32(I) CONST(20+I)=PR134A(I) CONST(40+I)=0.0D0 COEFF(I)=CR32(I) COEFF(20+I)=CR134A(I) COEFF(40+I)=0.0D0 20 CONTINUE DO 30 I=1,3 NAMES(I)=NR32(I) NAMES(3+I)=NR134A(I) NAMES(6+I)='NUL' 30 CONTINUE TYPES(1)=TR32 TYPES(2)=TR134A TYPES(3)=0 TYPES(4)=2 DO 60 I=1,6 TYPES(4+I)=0 60 CONTINUE DO 80 I=1,3 DO 90 K=1,3 CONST(58+3*I+K)=KIJ10(I,K) 90 CONTINUE 80 CONTINUE J=0 ELSE J=-2 CALL FMF062(J,'FMF000','THERE IS NO COMBINATION') ENDIF RETURN END c ----------------------------------------------------- c Temp,cmpnum => a_i c ----------------------------------------------------- DOUBLE PRECISION FUNCTION FMF001(CMPNUM,TEMP) DOUBLE PRECISION TEMP INTEGER CMPNUM DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C FMF001=COEFF(20*CMPNUM-19)*DEXP(COEFF(20*CMPNUM-18)*TEMP $ +COEFF(20*CMPNUM-17)*TEMP**2) RETURN END c --------------------------------------------- c Temperature,Componet number => b_i c --------------------------------------------- DOUBLE PRECISION FUNCTION FMF002(CMPNUM,TEMP) DOUBLE PRECISION TEMP INTEGER CMPNUM DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C FMF002=COEFF(20*CMPNUM-16)+COEFF(20*CMPNUM-15)*TEMP $ +COEFF(20*CMPNUM-14)*TEMP**2 RETURN END c ------------------------------------------- c Temperature,Mole fraction => a_{mix} c ------------------------------------------- DOUBLE PRECISION FUNCTION FMF003(MOLFR,TEMP) DOUBLE PRECISION MOLFR(1:3),TEMP,AMIX,FMF001 INTEGER I,J DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C AMIX=0.0D0 DO 10 I=1,3 DO 20 J=1,3 AMIX=AMIX+MOLFR(I)*MOLFR(J)*(1.0D0-CONST(58+3*I+J)) $ *(FMF001(I,TEMP)*FMF001(J,TEMP))**0.5 20 CONTINUE 10 CONTINUE FMF003=AMIX RETURN END c ---------------------------------------------- c Temperature,Mole fraction => b_{mix} c ---------------------------------------------- DOUBLE PRECISION FUNCTION FMF004(MOLFR,TEMP) DOUBLE PRECISION MOLFR(1:3),TEMP,BMIX,FMF002 INTEGER CMPNUM DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C BMIX=0.0D0 DO 10 CMPNUM=1,3 BMIX=BMIX+MOLFR(CMPNUM)*FMF002(CMPNUM,TEMP) 10 CONTINUE FMF004=BMIX RETURN END c ------------------------------------------------------- c Temperature,Componet number => Enthalpy of ideal gas c ------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF005(TEMP,CMPNUM) INTEGER CMPNUM DOUBLE PRECISION TEMP,H0,TEMP0,HIG,FMF060 DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C H0=0.0D0 TEMP0=298.15D0 IF (TYPES(CMPNUM).EQ.1) THEN HIG=COEFF(20*CMPNUM-13)*(TEMP-TEMP0)+0.5D0*COEFF(CMPNUM $ *20-12)*(TEMP**2-TEMP0**2)+COEFF(CMPNUM*20-11) $ *(TEMP**3-TEMP0**3)/3.0D0+0.25D0*COEFF(CMPNUM*20 $ -10)*(TEMP**4-TEMP0**4)+H0 ELSEIF (TYPES(CMPNUM).EQ.2) THEN HIG=FMF060(TEMP,CMPNUM)-FMF060(TEMP0,CMPNUM)+H0 ELSE CALL FMF062(-2,'FMF005','TYPES MISSING') HIG=-1.0D20 ENDIF FMF005=HIG RETURN END c ------------------------------------------------------------- c Temperature,Component number => Entropy of ideal gas (p=p0) c ------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF006(TEMP,CMPNUM) INTEGER CMPNUM DOUBLE PRECISION TEMP,S0,SIG,TEMP0,FMF061 DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C IF (TEMP.LE.0.0D0) THEN CALL FMF062(-2,'FMF006','TEMPERATURE IS LESS TAHN 0') FMF006=-1.0D20 RETURN ENDIF S0=0.0D0 TEMP0=298.15D0 IF (TYPES(CMPNUM).EQ.1) THEN SIG=COEFF(CMPNUM*20-13)*(DLOG(TEMP)-DLOG(TEMP0)) $ +COEFF(CMPNUM*20-12)*(TEMP-TEMP0)+0.5D0*COEFF(CMPNUM $ *20-11)*(TEMP**2-TEMP0**2)+COEFF(CMPNUM*20-10) $ *(TEMP**3-TEMP0**3)/3.0D0+S0 ELSEIF (TYPES(CMPNUM).EQ.2) THEN SIG=TEMP*(FMF061(TEMP,CMPNUM)-FMF061(TEMP0,CMPNUM))+S0 ELSE CALL FMF062(-2,'FMF006','TYPES MISSING') SIG=-1.0D20 ENDIF FMF006=SIG RETURN END c ------------------------------------------------ c Temperature,Volume,Mole fraction => Pressure c ------------------------------------------------ DOUBLE PRECISION FUNCTION FMF007(MOLFR,TEMP,VOL) DOUBLE PRECISION TEMP,VOL,MOLFR(1:3),A,B,Y,Y1,VB $ ,FMF003,FMF004 DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C A=FMF003(MOLFR,TEMP) B=FMF004(MOLFR,TEMP) IF (VOL.LE.1.0D-15) THEN CALL FMF062(-2,'FMF007','DIVIDED BY 0.(VOL)') FMF007=-1.0D20 RETURN ENDIF Y=B*0.25D0/VOL Y1=1.0D0-Y VB=VOL+B IF (Y1.EQ.0.0D0) THEN CALL FMF062(-2,'FMF007','DIVIDED BY 0.(Y1)') FMF007=-1.0D20 RETURN ELSEIF (VB.EQ.0.0D0) THEN CALL FMF062(-2,'FMF007','DIVIDED BY 0.(B+VOL)') FMF007=-1.0D20 RETURN ENDIF FMF007=(1.0D0+Y+Y**2-Y**3)*CONST(61)*TEMP/(VOL*Y1**3) $ -A/(VOL*(VOL+B)) RETURN END c -------------------------------------------------------- c Temperature,Volume,Mole fraction => (dp/dv)_{T,N} c -------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF008(MOLFR,TEMP,VOL) DOUBLE PRECISION TEMP,VOL,MOLFR(1:3),A,B,VB,VB4 $ ,FMF003,FMF004 DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C A=FMF003(MOLFR,TEMP) B=FMF004(MOLFR,TEMP) VB=B+VOL VB4=B-4.0D0*VOL IF (VOL.LE.1.0D-15) THEN CALL FMF062(-2,'FMF008','DVICED BY 0(VOL).') FMF008=-1.0D20 RETURN ELSEIF ((VB.EQ.0.0D0).OR.(VB4.EQ.0.0D0)) THEN CALL FMF062(-2,'FMF008','DVICED BY 0(VB OR VB4).') FMF008=-1.0D20 RETURN ENDIF C FMF008=(A*B**5-1.4D1*A*B**4*VOL+6.4D1*A*B**3*VOL**2-6.4D1 $ *A*B**2*VOL**3-2.56D2*A*B*VOL**4+5.12D2*A*VOL**5 $ -B**6*CONST(61)*TEMP+1.4D1*B**5*CONST(61)*TEMP*VOL $ -3.3D1*B**4*CONST(61)*TEMP*VOL**2 $ -3.68D2*B**3*CONST(61)*TEMP*VOL**3-8.32D2*B**2*CONST(61) $ *TEMP*VOL**4-7.68D2*B*CONST(61)*TEMP*VOL**5-2.56D2 $ *CONST(61)*TEMP*VOL**6)/(VB**2*VB4**4*VOL**2) RETURN END c ----------------------------------------------------------- c Temperature,Volume,Component number c => Gibbs free energy of ideal gas c ----------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF009(TEMP,VOL,CMPNUM) DOUBLE PRECISION TEMP,VOL,A,B,BETA,P0 $ ,FMF001,FMF002,FMF005,FMF006,VB,VBE INTEGER CMPNUM DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C IF (TEMP.LE.1.0D-15) THEN CALL FMF062(-2,'FMF009','TEMPERATURE IS LESS THAN 0') FMF009=-1.0D20 RETURN ELSEIF (VOL.LE.1.0D-15) THEN CALL FMF062(-2,'FMF009','VOLME IS LESS THAN 0') FMF009=-1.0D20 RETURN ENDIF P0=1.0D2 A=FMF001(CMPNUM,TEMP) B=FMF002(CMPNUM,TEMP) BETA=B*0.25D0 VBE=VOL-BETA VB=VOL+B IF (B.EQ.0.0D0) THEN CALL FMF062(-2,'FMF009','DIVIDED BY 0(B).') FMF009=-1.0D20 RETURN ELSEIF (VB.EQ.0.0D0.OR.VBE.EQ.0.0D0) THEN CALL FMF062(-2,'FMF009','DIVIDED BY 0(VB OR VBE).') FMF009=-1.0D20 RETURN ENDIF FMF009=FMF005(TEMP,CMPNUM)-TEMP*FMF006(TEMP $ ,CMPNUM)+CONST(61)*TEMP*(DLOG(CONST(61)*TEMP) $ -DLOG(P0*VOL))-A*(DLOG(VOL+B)-DLOG(VOL))/B+CONST(61)*TEMP $ *BETA*(8.0D0*VOL**2-9.0D0*VOL*BETA+3.0D0*BETA**2) $ /(VBE**3)-A/VB RETURN END c -------------------------------------------------------------- c Component number,Temperature,Volume,Mole fraction c => Chemical potential c -------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF010(MOLFR,TEMP,VOL,CMPNUM) DOUBLE PRECISION MOLFR(1:3),TEMP,VOL $ ,AI,A,AIJ,PRESS0,BI,BETAI,B,BETA,MUIG,VB,VBE $ ,FMF001,FMF002,FMF003,FMF004,FMF005,FMF006 INTEGER CMPNUM,I DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C IF (TEMP.LE.1.0D-15) THEN CALL FMF062(-2,'FMF010','TEMPERATURE IS LESS THAN 0') FMF010=-1.0D20 RETURN ELSEIF (VOL.LE.1.0D-15) THEN CALL FMF062(-2,'FMF010','VOLUME IS LESS THAN 0') FMF010=-1.0D20 RETURN ELSEIF (MOLFR(CMPNUM).LE.1.0D-15) THEN CALL FMF062(-2,'FMF010','MOLFR IS LESS THAN 0') FMF010=-1.0D20 RETURN ENDIF PRESS0=1.0D2 C AI=FMF001(CMPNUM,TEMP) A=FMF003(MOLFR,TEMP) AIJ=0.0D0 DO 10 I=1,3 AIJ=AIJ+MOLFR(I)*(1.0D0-CONST(58+3*CMPNUM+I))*(FMF001( $ CMPNUM,TEMP)*FMF001(I,TEMP))**0.5 10 CONTINUE BI=FMF002(CMPNUM,TEMP) BETAI=BI*0.25D0 B=FMF004(MOLFR,TEMP) BETA=B*0.25D0 VBE=VOL-BETA VB=VOL+B IF (B.EQ.0.0D0.OR.VBE.EQ.0.0D0) THEN CALL FMF062(-2,'FMF010','DIVIDED BY 0(B OR VBE).') FMF010=-1.0D20 RETURN ELSEIF (VB.EQ.0.0D0) THEN CALL FMF062(-2,'FMF010' $ ,'DIVIDED BY 0 OR LOG DOMEIN ERROR(VB).') FMF010=-1.0D20 RETURN ENDIF C MUIG=FMF005(TEMP,CMPNUM)-TEMP*FMF006(TEMP $ ,CMPNUM)+CONST(61)*TEMP*DLOG(MOLFR(CMPNUM)) FMF010=MUIG+CONST(61)*TEMP*(DLOG(CONST(61)*TEMP)-DLOG(PRESS0 $ *VOL))+CONST(61)*TEMP*BETA*(4.0D0*VOL-3.0D0*BETA) $ /(VBE**2)+CONST(61)*TEMP*BETAI*(4.0D0*VOL**2 $ -2.0D0*VOL*BETA)/(VBE**3)+A*BI*(DLOG(VB) $ -DLOG(VOL))/B**2-A*BI/(B*(VOL+B))+2.0D0*AIJ*(DLOG(VOL) $ -DLOG(VOL+B))/B RETURN END c ----------------------------------------------------------- c Temperature,Volume,Mole fraction => (d2p/dv2)_{T,N} c ----------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF011(MOLFR,T,V) DOUBLE PRECISION MOLFR(1:3),T,V,A,B,R,VB,VB4 $ ,FMF003,FMF004 DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C A=FMF003(MOLFR,T) B=FMF004(MOLFR,T) R=CONST(61) VB=V+B VB4=B-4.0D0*V IF (V.LE.1.0D-15) THEN CALL FMF062(-2,'FMF011','DIVIDED BY 0(V).') FMF011=-1.0D20 RETURN ELSEIF (VB.EQ.0.0D0.OR.VB4.EQ.0.0D0) THEN CALL FMF062(-2,'FMF011','DIVIDED BY 0(VB OR VB4).') FMF011=-1.0D20 RETURN ENDIF C FMF011=-1.0D0*(2.0D0*(A*B**7-17.0D0*A*B**6*V+103.0D0*A*B**5*V**2 $ -220.0D0*A*B**4*V**3-160.0D0*A*B**3*V**4+896.0D0*A*B**2*V**5 $ +768.0D0*A*B*V**6-3072.0D0*A*V**7-B**8*R*T+17.0D0*B**7 $ *R*T*V-103.0D0*B**6*R*T*V**2+219.0D0*B**5*R*T*V**3 $ +3252.0D0*B**4*R*T*V**4+8160.0D0*B**3*R*T*V**5+9088.0D0 $ *B**2*R*T*V**6+4864.0D0*B*R*T*V**7+1024.0D0*R*T*V**8)) $ /(VB**3*VB4**5*V**3) RETURN END c --------------------------------------------------------------- c Temperature => (da_i/dT) c --------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF012(CMPNUM,TEMP) DOUBLE PRECISION TEMP INTEGER CMPNUM DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C FMF012=COEFF(20*CMPNUM-19)*(COEFF(20*CMPNUM-18) $ +2.0D0*COEFF(20*CMPNUM-17)*TEMP) $ *DEXP(COEFF(20*CMPNUM-18)*TEMP+COEFF(20*CMPNUM-17)*TEMP**2) RETURN END c --------------------------------------------------------------- c Temperature => (db_i/dT) c --------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF013(CMPNUM,TEMP) DOUBLE PRECISION TEMP INTEGER CMPNUM DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C FMF013=COEFF(20*CMPNUM-15)+2.0D0*COEFF(20*CMPNUM-14)*TEMP RETURN END c --------------------------------------------------------------- c Temperature,Mole fraction => (da/dT) c --------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF014(MOLFR,TEMP) DOUBLE PRECISION MOLFR(1:3),TEMP,AI(1:3),DAI(1:3),DA,AIJ $ ,FMF001,FMF012 INTEGER I,J DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C DO 10 I=1,3 AI(I)=FMF001(I,TEMP) DAI(I)=FMF012(I,TEMP) 10 CONTINUE DA=0.0D0 DO 20 I=1,2 DO 30 J=1,2 AIJ=AI(I)*AI(J) IF (AIJ.EQ.0.0D0) THEN CALL FMF062(-2,'FMF014','DIVIDED BY 0(AIJ).') FMF014=-1.0D20 RETURN ENDIF DA=DA+0.5D0*MOLFR(I)*MOLFR(J)*(1.0D0-CONST(58+3*I+J)) $ *(DAI(I)*AI(J)+AI(I)*DAI(J))/AIJ**0.5D0 30 CONTINUE 20 CONTINUE FMF014=DA RETURN END c --------------------------------------------------------------- c Temperature,Mole fraction => (db/dT) c --------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF015(MOLFR,TEMP) DOUBLE PRECISION MOLFR(1:3),TEMP,DBI(1:3),DB $ ,FMF013 INTEGER I DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C DO 10 I=1,3 DBI(I)=FMF013(I,TEMP) 10 CONTINUE DB=0.0D0 DO 20 I=1,3 DB=DB+MOLFR(I)*DBI(I) 20 CONTINUE FMF015=DB RETURN END c --------------------------------------------------------------- c Temperature => (d2a_i/dT2) c --------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF016(CMPNUM,TEMP) DOUBLE PRECISION TEMP INTEGER CMPNUM DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C FMF016=(2.0D0*COEFF(20*CMPNUM-19)*COEFF(20*CMPNUM-17) $ +COEFF(20*CMPNUM-19)*(COEFF(20*CMPNUM-18) $ +2.0D0*COEFF(20*CMPNUM-17)*TEMP)**2) $ *DEXP(COEFF(20*CMPNUM-18)*TEMP+COEFF(20*CMPNUM-17)*TEMP**2) RETURN END c --------------------------------------------------------------- c (d2b_i/dT2) c --------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF017(CMPNUM) INTEGER CMPNUM DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C FMF017=2.0D0*COEFF(20*CMPNUM-14) RETURN END c --------------------------------------------------------------- c Temperature,Mole fraction => (d2a/dT2) c --------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF018(MOLFR,TEMP) DOUBLE PRECISION MOLFR(1:3),TEMP,AI(1:3),DAI(1:3) $ ,DDAI(1:3),DDA,AIJ $ ,FMF001,FMF012,FMF016 INTEGER I,J DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C DO 10 I=1,3 AI(I)=FMF001(I,TEMP) DAI(I)=FMF012(I,TEMP) DDAI(I)=FMF016(I,TEMP) 10 CONTINUE DDA=0.0D0 DO 20 I=1,3 DO 30 J=1,3 AIJ=AI(I)*AI(J) IF (AIJ.EQ.0.0D0) THEN CALL FMF062(-2,'FMF018','DIVIDED BY 0(AIJ).') FMF018=-1.0D20 RETURN ENDIF DDA=DDA+0.5D0*MOLFR(I)*MOLFR(J)*(1.0D0-CONST(58+3*I+J)) $ *((DDAI(I)*AI(J)+2.0D0*DAI(I)*DAI(J)+AI(I)*DDAI(J)) $ -0.5D0*(DAI(I)*AI(J)+AI(I)*DAI(J))**2/AIJ) $ /AIJ**0.5D0 30 CONTINUE 20 CONTINUE FMF018=DDA RETURN END c --------------------------------------------------------------- c Mole fraction => (d2b/dT2) c --------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF019(MOLFR) DOUBLE PRECISION MOLFR(1:3),DDBI(1:3),DDB $ ,FMF017 INTEGER I DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C DO 10 I=1,3 DDBI(I)=FMF017(I) 10 CONTINUE DDB=0.0D0 DO 20 I=1,3 DDB=DDB+MOLFR(I)*DDBI(I) 20 CONTINUE FMF019=DDB RETURN END c --------------------------------------------------------- c Temperature,Volume,Mole fraction => Enthalpy c --------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF020(MOLFR,TEMP,VOL) DOUBLE PRECISION MOLFR(1:3),TEMP,VOL $ ,A,B,DA,DB,BETA,DBETA,HIG,H1,H2,H3,VB,VBE $ ,FMF003,FMF004,FMF014,FMF015,FMF005 INTEGER I DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C A=FMF003(MOLFR,TEMP) B=FMF004(MOLFR,TEMP) DA=FMF014(MOLFR,TEMP) DB=FMF015(MOLFR,TEMP) BETA=B*0.25D0 DBETA=DB*0.025D0 HIG=0.0D0 DO 10 I=1,TYPES(4) HIG=HIG+MOLFR(I)*FMF005(TEMP,I) 10 CONTINUE VB=VOL+B VBE=VOL-BETA IF (VOL.LE.1.0D-15.OR.VB.EQ.0.0D0) THEN CALL FMF062(-2,'FMF020','LOG DOMAIN ERROR(VOL OR VB).') FMF020=-1.0D20 RETURN ELSEIF (VBE.EQ.0.0D0) THEN CALL FMF062(-2,'FMF020','DIVIDED BY 0(VBE).') FMF020=-1.0D20 RETURN ENDIF H1=(DA*B*TEMP-A*DB*TEMP-A*B)*(DLOG(VB)-DLOG(VOL))/B**2 H2=(A*DB*TEMP-A*B)/(B*VB) H3=CONST(61)*TEMP*(4.0D0*VOL**2-2.0D0*VOL*BETA) $ *(BETA-DBETA*TEMP)/(VBE**3) FMF020=HIG+H1+H2+H3 RETURN END c --------------------------------------------------------- c Temperature,Volume,Mole fraction => Entropy c --------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF021(MOLFR,TEMP,VOL) DOUBLE PRECISION MOLFR(1:3),TEMP,VOL $ ,A,B,DA,DB,BETA,DBETA,SIG,S1,S2,S3,S4,VB,VBE $ ,FMF003,FMF004,FMF014,FMF015,FMF006 INTEGER I DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C A=FMF003(MOLFR,TEMP) B=FMF004(MOLFR,TEMP) DA=FMF014(MOLFR,TEMP) DB=FMF015(MOLFR,TEMP) BETA=B*0.25D0 DBETA=DB*0.025D0 SIG=0.0D0 DO 10 I=1,TYPES(4) IF (MOLFR(I).LE.1.0D-15) THEN CALL FMF062(-2,'FMF021','MOLFR IS LESS THAN 0.') FMF021=-1.0D20 RETURN ENDIF SIG=SIG+MOLFR(I)*FMF006(TEMP,I) $ -CONST(61)*MOLFR(I)*DLOG(MOLFR(I)) 10 CONTINUE VB=VOL+B VBE=VOL-BETA IF (B.EQ.0.0D0.OR.VBE.EQ.0.0D0) THEN CALL FMF062(-2,'FMF021','DIVIDED BY 0(B OR VBE).') FMF021=-1.0D20 RETURN ELSEIF (VB.EQ.0.0D0.OR.VOL.EQ.0.0D0) THEN CALL FMF062(-2,'FMF021','LOG DOMAIN ERROR(VB OR VOL).') FMF021=-1.0D20 RETURN ENDIF S1=(DA*B-A*DB)*(DLOG(VB)-DLOG(VOL))/B**2 S2=(A*DB)/(B*VB) S3=CONST(61)*BETA*(4.0D0*VOL-3.0D0*BETA)/(VBE**2) S4=CONST(61)*TEMP*DBETA*(4.0D0*VOL**2-2.0D0*VOL*BETA) $ /(VBE**3) FMF021=SIG+S1+S2-S3-S4 RETURN END c ---------------------------------------------------------- c (dp/dv)_{T,N}=0 => Volume(1:2) c ---------------------------------------------------------- SUBROUTINE FMF022(J,MOLFR,TEMP,VOL) DOUBLE PRECISION MOLFR(1:3),TEMP $ ,VOL(1:2),EP1,EP2,V0,V1,V2,DP0,DP1,DP2,VFL,VFLG,DPFL $ ,AL,BE,DPFLG $ ,FMF004,FMF008 INTEGER I,J DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C EP1=1.0D-14 EP2=1.0D-2 V0=FMF004(MOLFR,TEMP)*0.25D0+1.0D-7 V1=V0+1.0D-7 DP0=FMF008(MOLFR,TEMP,V0) DP1=FMF008(MOLFR,TEMP,V1) DO 10 I=1,100 V2=V1-DP1*(V0-V1)/(DP0-DP1) DP2=FMF008(MOLFR,TEMP,V2) VFL=V2**2 IF (VFL.GT.1.0D0) THEN VFLG=((V1-V2)/V2)**2 DPFL=DP2**2 IF ((VFLG.LT.EP1).AND.(DPFL.LT.EP2)) THEN J=0 VOL(1)=V2 GOTO 20 ENDIF ELSE VFLG=(V2-V1)**2 DPFL=DP2**2 IF ((VFLG.LT.EP1).AND.(DPFL.LT.EP2)) THEN J=0 VOL(1)=V2 GOTO 20 ENDIF ENDIF V0=V1 V1=V2 DP0=DP1 DP1=DP2 10 CONTINUE J=-1 CALL FMF062(J,'FMF022','VOL(1), NO CONVERGENCE') RETURN 20 AL=1.0D-5 50 V0=VOL(1)+AL DP0=FMF008(MOLFR,TEMP,V0) IF (DP0.LE.0.0D0) THEN AL=AL*0.9D0 GOTO 50 ENDIF BE=1.0D0 30 V1=BE DP1=FMF008(MOLFR,TEMP,V1) IF (DP1.GE.0.0D0) THEN BE=BE*5.0D0 GOTO 30 ENDIF DO 40 I=1,1000 V2=(V0+V1)*0.5D0 DP2=FMF008(MOLFR,TEMP,V2) DPFL=DP2*DP0 IF (DPFL.GT.0.0D0) THEN V0=V2 DP0=DP2 ELSE V1=V2 DP1=DP2 ENDIF V2=V1-DP1*(V0-V1)/(DP0-DP1) DP2=FMF008(MOLFR,TEMP,V2) DPFL=DP2*DP0 IF (DPFL.GT.0.0D0) THEN V0=V2 DP0=DP2 ELSE V1=V2 DP1=DP2 ENDIF VFL=V2**2 IF (VFL.GT.1.0D0) THEN VFLG=((V1-V0)/V1)**2 DPFLG=DP2**2 IF ((VFLG.LT.EP1).AND.(DPFLG.LT.EP2)) THEN J=0 VOL(2)=V2 RETURN ENDIF ELSE VFLG=(V1-V0)**2 DPFLG=DP2**2 IF ((VFLG.LT.EP1).AND.(DPFLG.LT.EP2)) THEN J=0 VOL(2)=V2 RETURN ENDIF ENDIF 40 CONTINUE J=-1 CALL FMF062(J,'FMF022','VOL(2), NO CONVERGENCE') RETURN END c ---------------------------------------------------------------- c (d2p/dv2)_{T,N}=0 ====> Volume c ---------------------------------------------------------------- SUBROUTINE FMF023(J,MOLFR,TEMP,VOL) INTEGER I,J DOUBLE PRECISION MOLFR(1:3),TEMP,VOL $ ,EP1,EP2,VOL0,VOL1,VOL2,VOLEP1,VOLEP2,D2P0,D2P1,D2P2 $ ,D2PEP $ ,FMF004,FMF011 DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C EP1=1.0D-14 EP2=1.0D-2 VOL0=FMF004(MOLFR,TEMP)+1.0D-7 VOL1=VOL0+1.0D-7 D2P0=FMF011(MOLFR,TEMP,VOL0) D2P1=FMF011(MOLFR,TEMP,VOL1) DO 10 I=1,100 VOL2=VOL1-D2P1*(VOL0-VOL1)/(D2P0-D2P1) D2P2=FMF011(MOLFR,TEMP,VOL2) VOLEP1=VOL2**2 IF (VOLEP1.GT.1.0D0) THEN VOLEP2=((VOL2-VOL1)/VOL2)**2 D2PEP=D2P2**2 IF ((VOLEP2.LE.EP1).AND.(D2PEP.LE.EP2)) THEN J=0 VOL=VOL2 RETURN ENDIF ELSE VOLEP2=(VOL2-VOL1)**2 D2PEP=D2P2**2 IF ((VOLEP2.LE.EP1).AND.(D2PEP.LE.EP2)) THEN J=0 VOL=VOL2 RETURN ENDIF ENDIF VOL0=VOL1 VOL1=VOL2 D2P0=D2P1 D2P1=D2P2 10 CONTINUE J=-1 CALL FMF062(J,'FMF023','VOL, NO CONVERGENCE') VOL=FMF004(MOLFR,TEMP)*0.25D0+1.0D-2 RETURN END c --------------------------------------------------------------- c p=p =====> vr(1:2) c type=1 : liquid region c 2 : gas region c --------------------------------------------------------------- SUBROUTINE FMF024(J,MOLFR,TEMP,PRESS,VOLPOL,TYPE,VOL) DOUBLE PRECISION MOLFR(1:3),TEMP,PRESS $ ,VOLPOL,VOL(1:2),EP1,EP2,V0,V1,V2,P0,P1,P2,VFL,VFLG,PFL $ ,PFLG,BE,VE $ ,FMF004,FMF007,FMF008 INTEGER I,J,TYPE DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C EP1=1.0D-14 EP2=1.0D-2 IF (TYPE.EQ.2) THEN VOL(1)=-1.0D40 GOTO 20 ENDIF V1=FMF004(MOLFR,TEMP)*0.25D0+1.0D-6 P1=FMF007(MOLFR,TEMP,V1)-PRESS DO 10 I=1,100 V2=V1-P1/FMF008(MOLFR,TEMP,V1) P2=FMF007(MOLFR,TEMP,V2)-PRESS VE=V2**2 IF (VE.LT.1.0D0) THEN VFL=(V2-V1)**2 ELSE VFL=((V2-V1)/V2)**2 ENDIF IF (VFL.LT.EP1) THEN PFL=P2**2 IF (PFL.LT.EP2) THEN VOL(1)=V2 GOTO 25 ENDIF ENDIF V1=V2 P1=P2 10 CONTINUE J=-1 CALL FMF062(J,'FMF024','VOL(1), NO CONVERGENCE') RETURN 25 IF (TYPE.EQ.1) THEN J=0 VOL(2)=-1.0D40 RETURN ENDIF 20 V0=VOLPOL P0=FMF007(MOLFR,TEMP,V0)-PRESS c IF (P0.LT.0) THEN c CALL FMF062(-2,'FMF024','PRESS(MAX) IS LESS THAN 0.') c VOL(2)=-1.0D40 c RETURN c ENDIF BE=1.2D0 60 V1=CONST(61)*TEMP/PRESS*BE P1=FMF007(MOLFR,TEMP,V1)-PRESS IF (P1.GE.0.0D0) THEN BE=BE*1.5D0 GOTO 60 ENDIF DO 30 I=1,1000 V2=(V0+V1)*0.5D0 P2=FMF007(MOLFR,TEMP,V2)-PRESS PFL=P2*P0 IF (PFL.GT.0.0D0) THEN V0=V2 P0=P2 ELSE V1=V2 P1=P2 ENDIF V2=V1-P1*(V0-V1)/(P0-P1) P2=FMF007(MOLFR,TEMP,V2)-PRESS PFL=P2*P0 IF (PFL.GT.0.0D0) THEN V0=V2 P0=P2 ELSE V1=V2 P1=P2 ENDIF VFL=V1**2 IF (VFL.GT.1.0D0) THEN VFLG=((V1-V0)/V1)**2 PFLG=P2**2 IF ((VFLG.LT.EP1).AND.(PFLG.LT.EP2)) THEN J=0 VOL(2)=V2 RETURN ENDIF ELSE VFLG=(V1-V0)**2 PFLG=P2**2 IF ((VFLG.LT.EP1).AND.(PFLG.LT.EP2)) THEN J=0 VOL(2)=V2 RETURN ENDIF ENDIF 30 CONTINUE J=-1 CALL FMF062(J,'FMF024','VOL(2), NO CONVERGENCE') RETURN END c -------------------------------------------------------------- c Temperature,Volume,Component number => Pressure c -------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF025(CMPNUM,TEMP,VOL) DOUBLE PRECISION TEMP,VOL,A,B,Y $ ,FMF001,FMF002,VB,Y1 INTEGER CMPNUM DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C A=FMF001(CMPNUM,TEMP) B=FMF002(CMPNUM,TEMP) IF (VOL.LE.1.0D-15) THEN CALL FMF062(-2,'FMF025','DEVIDED BY 0(VOL)') FMF025=-1.0D20 RETURN ENDIF Y=B*0.25D0/VOL Y1=1.0D0-Y VB=VOL+B IF (VB.EQ.0.0D0.OR.Y1.EQ.0.0D0) THEN CALL FMF062(-2,'FMF025','DEVIDED BY 0(VB OR Y1).') FMF025=-1.0D20 RETURN ENDIF FMF025=(1.0D0+Y+Y**2-Y**3)*CONST(61)*TEMP/(VOL*Y1**3) $ -A/(VOL*VB) RETURN END c ---------------------------------------------------------------- c Temperature,Volume,Component number => (dv/dp) c ---------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF026(CMPNUM,T,V) DOUBLE PRECISION T,V,R,A,B $ ,FMF001,FMF002,VB,VB4 INTEGER CMPNUM DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C R=CONST(61) A=FMF001(CMPNUM,T) B=FMF002(CMPNUM,T) VB=V+B VB4=B-4.0D0*V IF (V.LE.1.0D-15.OR.VB.EQ.0.0D0) THEN CALL FMF062(-2,'FMF026','DIVIDED BY 0(V OR VB).') FMF026=-1.0D20 RETURN ELSEIF (VB4.EQ.0.0D0) THEN CALL FMF062(-2,'FMF026','DIVIDED BY 0(VB4).') FMF026=-1.0D20 RETURN ENDIF C FMF026=(A*B**5-1.4D1*A*B**4*V+6.4D1*A*B**3*V**2 $ -6.4D1*A*B**2*V**3 $ -2.56D2*A*B*V**4+5.12D2*A*V**5-B**6*R*T+1.4D1*B**5*R*T*V $ -3.3D1*B**4*R*T*V**2-3.68D2*B**3*R*T*V**3 $ -8.32D2*B**2*R*T*V**4-7.68D2*B*R*T*V**5-2.56D2*R*T*V**6) $ /(VB**2*VB4**4*V**2) RETURN END c -------------------------------------------------------------- c (dp/dv)_T=0 =====> Volume(1:2) c -------------------------------------------------------------- SUBROUTINE FMF027(J,CMPNUM,TEMP,VOL) DOUBLE PRECISION TEMP,VOL(1:2) $ ,EP1,EP2,V0,V1,V2,DP0,DP1,DP2,VFL,VFLG,DPFL,DPFLG,AL,BE $ ,FMF002,FMF026,FMF045 INTEGER I,J,CMPNUM DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C EP1=1.0D-14 EP2=1.0D-2 V1=FMF002(CMPNUM,TEMP)+1.0D-7 DP1=FMF026(CMPNUM,TEMP,V1) DO 10 I=1,1000 V2=V1-DP1/FMF045(CMPNUM,TEMP,V1) DP2=FMF026(CMPNUM,TEMP,V2) VFL=V2**2 IF (VFL.GT.1.0D0) THEN VFLG=((V2-V1)/V2)**2 DPFL=DP2**2 IF ((VFLG.LT.EP1).AND.(DPFL.LT.EP2)) THEN J=0 VOL(1)=V2 GOTO 20 ENDIF ELSE VFLG=(V2-V1)**2 DPFL=DP2**2 IF ((VFLG.LT.EP1).AND.(DPFL.LT.EP2)) THEN J=0 VOL(1)=V2 GOTO 20 ENDIF ENDIF V1=V2 DP1=DP2 10 CONTINUE J=-1 CALL FMF062(J,'FMF027','VOL(1), NO CONVERGENCE') RETURN 20 AL=1.0D-4 50 V0=VOL(1)+AL DP0=FMF026(CMPNUM,TEMP,V0) IF (DP0.LE.0.0D0) THEN AL=AL*0.9D0 GOTO 50 ENDIF BE=1.0D0 30 V1=BE DP1=FMF026(CMPNUM,TEMP,V1) IF (DP1.GE.0.0D0) THEN BE=BE*5.0D0 GOTO 30 ENDIF DO 40 I=1,1000 V2=(V0+V1)*0.5D0 DP2=FMF026(CMPNUM,TEMP,V2) DPFL=DP2*DP0 IF (DPFL.GT.0.0D0) THEN V0=V2 DP0=DP2 ELSE V1=V2 DP1=DP2 ENDIF V2=V1-DP1*(V0-V1)/(DP0-DP1) DP2=FMF026(CMPNUM,TEMP,V2) DPFL=DP2*DP0 IF (DPFL.GT.0.0D0) THEN V0=V2 DP0=DP2 ELSE V1=V2 DP1=DP2 ENDIF VFL=V1**2 IF (VFL.GT.1.0D0) THEN VFLG=((V0-V1)/V1)**2 DPFLG=DP2**2 IF ((VFLG.LT.EP1).AND.(DPFL.LT.EP2)) THEN J=0 VOL(2)=V2 RETURN ENDIF ELSE VFLG=(V0-V1)**2 DPFLG=DP2**2 IF ((VFLG.LT.EP1).AND.(DPFL.LT.EP2)) THEN J=0 VOL(2)=V2 RETURN ENDIF ENDIF 40 CONTINUE J=-1 CALL FMF062(J,'FMF027','VOL(2), NO CONVERGENCE') RETURN END c -------------------------------------------------------------- c Pressure,Component number => Volume(1:2) c -------------------------------------------------------------- SUBROUTINE FMF028(J,CMPNUM,TEMP,PRESS,VOLPOL,TYPE,VOL) DOUBLE PRECISION TEMP,PRESS,VOLPOL $ ,VOL(1:2),EP1,EP2,V0,V1,V2,P0,P1,P2,VFL,VFLG,PFL,PFLG $ ,BE $ ,FMF002,FMF025,FMF026 INTEGER I,J,CMPNUM,TYPE DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C EP1=1.0D-14 EP2=1.0D-2 IF (TYPE.EQ.2) THEN VOL(1)=-1.0D40 GOTO 20 ENDIF V1=FMF002(CMPNUM,TEMP)*0.25D0+1.0D-7 P1=FMF025(CMPNUM,TEMP,V1)-PRESS DO 10 I=1,100 V2=V1-P1/FMF026(CMPNUM,TEMP,V1) P2=FMF025(CMPNUM,TEMP,V2)-PRESS VFL=((V2-V1)/V2)**2 IF (VFL.LT.EP1) THEN PFL=P2**2 IF (PFL.LT.EP2) THEN VOL(1)=V2 GOTO 25 ENDIF ENDIF V1=V2 P1=P2 10 CONTINUE J=-1 CALL FMF062(J,'FMF028','VOL(1), NO CONVERGENCE') RETURN 25 IF (TYPE.EQ.1) THEN J=0 VOL(2)=-1.0D40 RETURN ENDIF 20 V0=VOLPOL IF (V0.LT.0.0D0) THEN V0=FMF002(CMPNUM,TEMP)*0.25+1.0D-7 ENDIF P0=FMF025(CMPNUM,TEMP,V0)-PRESS BE=1.2D0 60 V1=CONST(61)*TEMP/PRESS*BE P1=FMF025(CMPNUM,TEMP,V1)-PRESS IF (P1.GE.0.0D0) THEN BE=BE*1.5D0 GOTO 60 ENDIF DO 30 I=1,1000 V2=(V0+V1)*0.5D0 P2=FMF025(CMPNUM,TEMP,V2)-PRESS PFL=P2*P0 IF (PFL.GT.0.0D0) THEN V0=V2 P0=P2 ELSE V1=V2 P1=P2 ENDIF V2=V1-P1*(V0-V1)/(P0-P1) P2=FMF025(CMPNUM,TEMP,V2)-PRESS PFL=P2*P0 IF (PFL.GT.0.0D0) THEN V0=V2 P0=P2 ELSE V1=V2 P1=P2 ENDIF VFL=V1**2 IF (VFL.GT.1.0D0) THEN VFLG=((V1-V0)/V1)**2 PFLG=P2**2 IF ((VFLG.LE.EP1).AND.(PFLG.LE.EP2)) THEN J=0 VOL(2)=V2 RETURN ENDIF ELSE VFLG=(V1-V0)**2 PFLG=P2**2 IF ((VFLG.LE.EP1).AND.(PFLG.LE.EP2)) THEN J=0 VOL(2)=V2 RETURN ENDIF ENDIF 30 CONTINUE J=-1 CALL FMF062(J,'FMF028','VOL(2), NO CONVERGENCE') RETURN END c ---------------------------------------------------------------- c Temperature,Component number c => Saturated pressure,Liquid volume,Gas volume c if (dp/dv) does not converge,seting error code J=-2. c ---------------------------------------------------------------- SUBROUTINE FMF029(J,CMPNUM,TEMP,PRESS,VL,VV) INTEGER I,J,CMPNUM,FLG DOUBLE PRECISION TEMP,PRESS $ ,VL,VV $ ,EP1,EP2,A,B,VOLP(1:2),PMAX,PMIN,VMAX(1:2),VMIN(1:2) $ ,VMIN1(1:2),P0,P1,P2,DV0(1:2),DV1(1:2),DV2(1:2) $ ,G0,G1,G2,PFL,PFLG,GFL,GFLG $ ,FMF001,FMF002,FMF009,FMF025 DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES c c ----- the condition of convergence ----- EP1=1.0D-14 EP2=1.0D-2 c ---------- FLG=0 A=FMF001(CMPNUM,TEMP) B=FMF002(CMPNUM,TEMP) c ----- set maximum and minimum pressure ----- CALL FMF027(J,CMPNUM,TEMP,VOLP) IF (J.EQ.-1) THEN CALL FMF062(J,'FMF029','FMF027, NO CONVERGENCE') RETURN ENDIF PMAX=FMF025(CMPNUM,TEMP,VOLP(2)) PMIN=FMF025(CMPNUM,TEMP,VOLP(1)) IF (PMIN.LE.1.0D-3) THEN PMIN=1.0D-3 FLG=1 ENDIF c ----- vl and vv on maximun pressure ----- CALL FMF028(J,CMPNUM,TEMP,PMAX,VOLP(2),1,VMAX) IF (J.NE.0) THEN CALL FMF062(J,'FMF029','IN FMF028-1') RETURN ENDIF VMAX(2)=VOLP(2) c ----- vl and vv on minimum pressure ----- CALL FMF028(J,CMPNUM,TEMP,PMIN,VOLP(2),2,VMIN) IF (J.NE.0) THEN CALL FMF062(J,'FMF029','IN FMF028-2') RETURN ENDIF IF (FLG.EQ.0) THEN VMIN(1)=VOLP(1) ELSEIF (FLG.EQ.1) THEN CALL FMF028(J,CMPNUM,TEMP,PMIN,VOLP(2),1,VMIN1) IF (J.NE.0) THEN CALL FMF062(J,'FMF029','IN FMF028-3') RETURN ENDIF VMIN(1)=VMIN1(1) ENDIF c ----- main ----- P0=PMAX P1=PMIN DV0(1)=VMAX(1) DV0(2)=VMAX(2) DV1(1)=VMIN(1) DV1(2)=VMIN(2) G0=FMF009(TEMP,DV0(1),CMPNUM) $ -FMF009(TEMP,DV0(2),CMPNUM) G1=FMF009(TEMP,DV1(1),CMPNUM) $ -FMF009(TEMP,DV1(2),CMPNUM) GFL=G0*G1 IF (GFL.GT.0.0D0) THEN J=-2 CALL FMF062(J,'FMF029','GIBBS(LIQUID)*GIBBS(GAS)>0') RETURN ENDIF DO 80 I=1,100 P2=(P0+P1)*0.5D0 CALL FMF028(J,CMPNUM,TEMP,P2,VOLP(2),0,DV2) IF (J.NE.0) THEN CALL FMF062(J,'FMF029','IN FMF028(LOOP 80-1)') RETURN ENDIF G2=FMF009(TEMP,DV2(1),CMPNUM) $ -FMF009(TEMP,DV2(2),CMPNUM) GFL=G0*G2 IF (GFL.GT.0.0D0) THEN P0=P2 G0=G2 DV0(1)=DV2(1) DV0(2)=DV2(2) ELSE P1=P2 G1=G2 DV1(1)=DV2(1) DV1(2)=DV2(2) ENDIF P2=P1-G1*(P0-P1)/(G0-G1) CALL FMF028(J,CMPNUM,TEMP,P2,VOLP(2),0,DV2) IF (J.NE.0) THEN CALL FMF062(J,'FMF029','IN FMF028(LOOP 80-2)') RETURN ENDIF G2=FMF009(TEMP,DV2(1),CMPNUM) $ -FMF009(TEMP,DV2(2),CMPNUM) GFL=G0*G2 IF (GFL.GT.0.0D0) THEN P0=P2 G0=G2 DV0(1)=DV2(1) DV0(2)=DV2(2) ELSE P1=P2 G1=G2 DV1(1)=DV2(1) DV1(2)=DV2(2) ENDIF PFL=P2**2 IF (PFL.GT.1.0D0) THEN PFLG=((P1-P0)/P1)**2 GFLG=G2**2 IF ((PFLG.LE.EP1).AND.(GFLG.LE.EP2)) THEN J=0 PRESS=P2 VL=DV2(1) VV=DV2(2) RETURN ENDIF ELSE PFLG=(P1-P0)**2 GFLG=G2**2 IF ((PFLG.LE.EP1).AND.(GFLG.LE.EP2)) THEN J=0 PRESS=P2 VL=DV2(1) VV=DV2(2) RETURN ENDIF ENDIF 80 CONTINUE J=-1 CALL FMF062(J,'FMF029','NO CONVERGENCE') RETURN END c --------------------------------------------------------- c Temperature,Volume,Component number => Enthalpy c --------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF030(CMPNUM,TEMP,VOL) DOUBLE PRECISION TEMP,VOL $ ,A,B,DA,DB,BETA,DBETA,HIG,H1,H2,H3,VB,VBE $ ,FMF001,FMF002,FMF012,FMF013,FMF005 INTEGER CMPNUM DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C IF (VOL.LE.1.0D-15) THEN CALL FMF062(-2,'FMF030','VOLUME IS LESS THAN 0.') FMF030=-1.0D20 RETURN ENDIF A=FMF001(CMPNUM,TEMP) B=FMF002(CMPNUM,TEMP) DA=FMF012(CMPNUM,TEMP) DB=FMF013(CMPNUM,TEMP) BETA=B*0.25D0 DBETA=DB*0.025D0 VB=VOL+B VBE=VOL-BETA IF (VB.EQ.0.0D0.OR.VBE.EQ.0.0D0) THEN CALL FMF062(-2,'FMF030','DIVIDED BY 0(VB OR VBE).') FMF030=-1.0D20 RETURN ELSEIF (B.EQ.0.0D0) THEN CALL FMF062(-2,'FMF030','DIVIDED BY 0(B).') FMF030=-1.0D20 RETURN ENDIF HIG=FMF005(TEMP,CMPNUM) H1=(DA*B*TEMP-A*DB*TEMP-A*B)*(DLOG(VB)-DLOG(VOL))/B**2 H2=(A*DB*TEMP-A*B)/(B*VB) H3=CONST(61)*TEMP*(4.0D0*VOL**2-2.0D0*VOL*BETA) $ *(BETA-DBETA*TEMP)/(VBE**3) FMF030=HIG+H1+H2+H3 RETURN END c --------------------------------------------------------- c Temperature,Volume,Component number => Entropy c --------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF031(CMPNUM,TEMP,VOL) DOUBLE PRECISION TEMP,VOL $ ,A,B,DA,DB,BETA,DBETA,SIG,S1,S2,S3,S4,VB,VBE $ ,FMF001,FMF002,FMF012,FMF013,FMF006 INTEGER CMPNUM DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C IF (VOL.LE.1.0D-15) THEN CALL FMF062(-2,'FMF031','VOLUME IS LESS THAN 0.') FMF031=-1.0D20 RETURN ENDIF A=FMF001(CMPNUM,TEMP) B=FMF002(CMPNUM,TEMP) DA=FMF012(CMPNUM,TEMP) DB=FMF013(CMPNUM,TEMP) BETA=B*0.25D0 DBETA=DB*0.025D0 VB=VOL+B VBE=VOL-BETA IF (B.EQ.0.0D0.OR.VB.EQ.0.0D0) THEN CALL FMF062(-2,'FMF031','DIVIDED BY 0(B OR VB).') FMF031=-1.0D20 RETURN ELSEIF (VBE.EQ.0.0D0) THEN CALL FMF062(-2,'FMF031','DIVIDED BY 0(VBE).') FMF031=-1.0D20 RETURN ENDIF SIG=FMF006(TEMP,CMPNUM) S1=(DA*B-A*DB)*(DLOG(VB)-DLOG(VOL))/B**2 S2=(A*DB)/(B*VB) S3=CONST(61)*BETA*(4.0D0*VOL-3.0D0*BETA)/(VBE**2) S4=CONST(61)*TEMP*DBETA*(4.0D0*VOL**2-2.0D0*VOL*BETA) $ /(VBE**3) FMF031=SIG+S1+S2-S3-S4 RETURN END c ---------------------------------------------------------------- c Temperature,Pressure,Component number => Volume,Enthalpy,Entropy c (Saturated region => saturated liquid) c ---------------------------------------------------------------- SUBROUTINE FMF032(J,CMPNUM,TEMP,PRESS,VOL,ENTHAL,ENTROP) DOUBLE PRECISION TEMP,PRESS,VOL $ ,ENTHAL,ENTROP,VOLP(1:2),VOLM(1:2),PSAT,VL,VV $ ,FMF030,FMF031 INTEGER J,CMPNUM DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C J=0 CALL FMF029(J,CMPNUM,TEMP,PSAT,VL,VV) IF (J.NE.0) THEN CALL FMF062(J,'FMF032','IN FMF029-1') RETURN ENDIF IF (PRESS.LT.PSAT) THEN CALL FMF027(J,CMPNUM,TEMP,VOLP) IF (J.NE.0) THEN CALL FMF062(J,'FMF032','IN FMF027-1') RETURN ENDIF CALL FMF028(J,CMPNUM,TEMP,PRESS,VOLP(2),2,VOLM) IF (J.NE.0) THEN CALL FMF062(J,'FMF032','IN FMF028-1') RETURN ENDIF VOL=VOLM(2) ENTHAL=FMF030(CMPNUM,TEMP,VOL) ENTROP=FMF031(CMPNUM,TEMP,VOL) ELSEIF (PRESS.EQ.PSAT) THEN VOL=VL ENTHAL=FMF030(CMPNUM,TEMP,VOL) ENTROP=FMF031(CMPNUM,TEMP,VOL) ELSE CALL FMF028(J,CMPNUM,TEMP,PRESS,1.0D0,1,VOLM) IF (J.NE.0) THEN CALL FMF062(J,'FMF032','IN FMF028-2') RETURN ENDIF VOL=VOLM(1) ENTHAL=FMF030(CMPNUM,TEMP,VOL) ENTROP=FMF031(CMPNUM,TEMP,VOL) ENDIF RETURN END c --------------------------------------------------------------- c Pressure,Temperature => x,y c --------------------------------------------------------------- SUBROUTINE FMF033(J,TEMP,PRESS,MOLX,MOLY,VOLL,VOLV) DOUBLE PRECISION TEMP,PRESS,MOLX(1:3) $ ,MOLY(1:3),MOLX0(1:3),MOLX1(1:3),MOLX2(1:3),MOLY0(1:3) $ ,MOLY1(1:3),MOLY2(1:3),MOLFL,DMUFL,DMU0,DMU1,DMU2,EP1,EP2 $ ,VOLL0,VOLL1,VOLL2,VOLV0,VOLV1,VOLV2,VOLL,VOLV INTEGER I,J DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C J=0 EP1=1.0D-14 EP2=1.0D-2 MOLX0(1)=1.0D-7 MOLX0(2)=1.0D0-MOLX0(1) MOLX0(3)=0.0D0 CALL FMF035(J,MOLX0,TEMP,PRESS,DMU0,MOLY0,VOLL0,VOLV0) IF (J.NE.0) THEN WRITE(*,*) ' ERROR IN FMF033(FMF035)-1 J=',J RETURN ENDIF MOLX1(1)=2.0D-7 MOLX1(2)=1.0D0-MOLX1(1) MOLX1(3)=0.0D0 CALL FMF035(J,MOLX1,TEMP,PRESS,DMU1,MOLY1,VOLL1,VOLV1) IF (J.NE.0) THEN WRITE(*,*) ' ERROR IN FMF033(FMF035)-2 J=',J RETURN ENDIF DO 10 I=1,1000 MOLX2(1)=MOLX1(1)-DMU1*(MOLX0(1)-MOLX1(1))/(DMU0-DMU1) MOLX2(2)=1.0D0-MOLX2(1) MOLX2(3)=0.0D0 CALL FMF035(J,MOLX2,TEMP,PRESS,DMU2,MOLY2,VOLL2,VOLV2) IF (J.NE.0) THEN WRITE(*,*) ' ERROR IN FMF033(FMF035)-3 J=',J RETURN ENDIF MOLFL=(MOLX2(1)-MOLX1(1))**2 DMUFL=DMU2**2 IF ((MOLFL.LT.EP1).AND.(DMUFL.LT.EP2)) THEN J=0 VOLL=VOLL2 VOLV=VOLV2 MOLX(1)=MOLX2(1) MOLX(2)=MOLX2(2) MOLX(3)=MOLX2(3) MOLY(1)=MOLY2(1) MOLY(2)=MOLY2(2) MOLY(3)=MOLY2(3) RETURN ENDIF MOLX0(1)=MOLX1(1) MOLX0(2)=MOLX1(2) MOLX0(3)=MOLX1(3) MOLX1(1)=MOLX2(1) MOLX1(2)=MOLX2(2) MOLX1(3)=MOLX2(3) DMU0=DMU1 DMU1=DMU2 10 CONTINUE J=-1 RETURN END c ----------------------------------------------------------------- c Temperature,Pressure,Chemical potential => Mole fraction c ( \mu_2^v => y_{\mu_2^v} ) c ----------------------------------------------------------------- SUBROUTINE FMF034(J,TEMP,PRESS,CHMPOT,MOLFR,VOLV) INTEGER I,J DOUBLE PRECISION TEMP,PRESS,CHMPOT $ ,MOLFR(1:3),EP1,EP2,MOLFR0(1:3),MOLFR1(1:3),MOLFR2(1:3) $ ,VOLP0(1:2),VOLP1(1:2),VOLP2(1:2),VOL0(1:2),VOL1(1:2) $ ,VOL2(1:2),MU2V0,MU2V1,MU2V2,MOLFL,MUFL,VOLV $ ,FMF010 DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C J=0 EP1=1.0D-14 EP2=1.0D-2 MOLFR0(1)=1.0D0-1.0D-7 MOLFR0(2)=1.0D0-MOLFR0(1) MOLFR0(3)=0.0D0 MOLFR1(1)=1.0D0-2.0D-7 MOLFR1(2)=1.0D0-MOLFR1(1) MOLFR1(3)=0.0D0 CALL FMF022(J,MOLFR0,TEMP,VOLP0) IF (J.NE.0) THEN WRITE(*,*) ' ERROR IN FMF034(FMF022)-1 J=',J RETURN ENDIF CALL FMF022(J,MOLFR1,TEMP,VOLP1) IF (J.NE.0) THEN WRITE(*,*) ' ERROR IN FMF034(FMF022)-2 J=',J RETURN ENDIF CALL FMF024(J,MOLFR0,TEMP,PRESS,VOLP0(2),2,VOL0) IF (J.NE.0) THEN WRITE(*,*) ' ERROR IN FMF034(FMF024)-1 J=',J RETURN ENDIF CALL FMF024(J,MOLFR1,TEMP,PRESS,VOLP1(2),2,VOL1) IF (J.NE.0) THEN WRITE(*,*) ' ERROR IN FMF034(FMF024)-2 J=',J RETURN ENDIF MU2V0=FMF010(MOLFR0,TEMP,VOL0(2),2)-CHMPOT MU2V1=FMF010(MOLFR1,TEMP,VOL1(2),2)-CHMPOT DO 10 I=1,1000 MOLFR2(1)=MOLFR1(1)-MU2V1*(MOLFR0(1)-MOLFR1(1))/(MU2V0-MU2V1) MOLFR2(2)=1.0D0-MOLFR2(1) MOLFR2(3)=0.0D0 CALL FMF022(J,MOLFR2,TEMP,VOLP2) IF (J.NE.0) THEN WRITE(*,*) ' ERROR IN FMF034(FMF022)-3 J=',J RETURN ENDIF CALL FMF024(J,MOLFR2,TEMP,PRESS,VOLP2(2),2,VOL2) IF (J.NE.0) THEN WRITE(*,*) ' ERROR IN FMF034(FMF024)-3 J=',J RETURN ENDIF MU2V2=FMF010(MOLFR2,TEMP,VOL2(2),2)-CHMPOT MOLFL=(MOLFR1(1)-MOLFR2(1))**2 MUFL=MU2V2**2 IF ((MOLFL.LT.EP1).AND.(MUFL.LT.EP2)) THEN J=0 VOLV=VOL2(2) MOLFR(1)=MOLFR2(1) MOLFR(2)=MOLFR2(2) MOLFR(3)=MOLFR2(3) RETURN ENDIF MOLFR0(1)=MOLFR1(1) MOLFR0(2)=MOLFR1(2) MOLFR0(3)=MOLFR1(3) MU2V0=MU2V1 MOLFR1(1)=MOLFR2(1) MOLFR1(2)=MOLFR2(2) MOLFR1(3)=MOLFR2(3) MU2V1=MU2V2 10 CONTINUE J=-1 RETURN END c --------------------------------------------------------- c Temperature,Pressure,x => \mu_1^l-\mu_1^v c --------------------------------------------------------- SUBROUTINE FMF035(J,MOLX,TEMP,PRESS,DMU,MOLY,VOLL,VOLV) DOUBLE PRECISION MOLX(1:3),TEMP $ ,PRESS,DMU,MOLY(1:3),VOLP(1:2),VOL(1:2),MU1L,MU2L,MU1V $ ,FMF010,VOLL,VOLV INTEGER J DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C J=0 CALL FMF022(J,MOLX,TEMP,VOLP) IF (J.NE.0) THEN WRITE(*,*) ' ERROR IN FMF035(FMF022)-1 J=',J RETURN ENDIF CALL FMF024(J,MOLX,TEMP,PRESS,VOLP(2),1,VOL) IF (J.NE.0) THEN WRITE(*,*) ' ERROR IN FMF035(FMF024)-1 J=',J RETURN ENDIF MU1L=FMF010(MOLX,TEMP,VOL(1),1) MU2L=FMF010(MOLX,TEMP,VOL(1),2) VOLL=VOL(1) CALL FMF034(J,TEMP,PRESS,MU2L,MOLY,VOLV) IF (J.NE.0) THEN WRITE(*,*) ' ERROR IN FMF035(FMF034)-1 J=',J RETURN ENDIF MU1V=FMF010(MOLY,TEMP,VOLV,1) DMU=MU1L-MU1V RETURN END c --------------------------------------------------------------- c Pressure,Component number => Saturated temperature c --------------------------------------------------------------- SUBROUTINE FMF036(J,CMPNUM,PRESS,TEMP,VL,VV) DOUBLE PRECISION PRESS,TEMP $ ,VL,VV,EP1,EP2,T0,T1,T2,P0,P1,P2,PP0,PP1,PP2 $ ,TFL,PFL,FMF059 INTEGER I,J,CMPNUM DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C J=0 EP1=1.0D-14 EP2=1.0D-14 T0=FMF059(CMPNUM,PRESS) c T0=(-CONST(20*CMPNUM-13)+DSQRT(CONST(20*CMPNUM-13)**2 c $ -4.0D0*CONST(20*CMPNUM-12)*(CONST(20*CMPNUM-14)-PRESS))) c $ *0.5D0/CONST(20*CMPNUM-12) c T0=DEXP(CONST(20*CMPNUM-9)+CONST(20*CMPNUM-8)*DLOG(PRESS)) T1=T0-1.0D-4 CALL FMF029(J,CMPNUM,T0,PP0,VL,VV) IF (J.NE.0) THEN CALL FMF062(J,'FMF036','IN FMF029-1') RETURN ENDIF CALL FMF029(J,CMPNUM,T1,PP1,VL,VV) IF (J.NE.0) THEN CALL FMF062(J,'FMF036','IN FMF029-2') RETURN ENDIF P0=PP0-PRESS P1=PP1-PRESS DO 10 I=1,1000 T2=T1-P1*(T0-T1)/(P0-P1) CALL FMF029(J,CMPNUM,T2,PP2,VL,VV) IF (J.NE.0) THEN CALL FMF062(J,'FMF036','IN FMF029-3') RETURN ENDIF P2=PP2-PRESS TFL=((T1-T2)/T2)**2 IF (TFL.LT.EP1) THEN PFL=P2**2 IF (PFL.LT.EP2) THEN TEMP=T2 RETURN ENDIF ENDIF T0=T1 P0=P1 T1=T2 P1=P2 10 CONTINUE J=-1 RETURN END c --------------------------------------------------------- c Temperature,Mole fraction => Bubble point pressure c --------------------------------------------------------- SUBROUTINE FMF037(J,TEMP,MOLFR,PRESS,VOL) DOUBLE PRECISION TEMP,MOLFR(1:3),PRESS $ ,EP1,EP2,VL,VV,PSAT1,PSAT2,P0,P1,P2 $ ,MOLX0(1:3),MOLX1(1:3),MOLX2(1:3),MX0,MX1,MX2,MFL,MFLG $ ,PFL,MOLY2(1:3),VL2,VV2,VOL INTEGER I,J DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C J=0 EP1=1.0D-14 EP2=1.0D-2 CALL FMF029(J,1,TEMP,PSAT1,VL,VV) IF (J.NE.0) THEN CALL FMF062(J,'FMF037','IN FMF029-1') RETURN ENDIF CALL FMF029(J,2,TEMP,PSAT2,VL,VV) IF (J.NE.0) THEN CALL FMF062(J,'FMF037','IN FMF029-2') RETURN ENDIF P0=PSAT2 P1=PSAT1 MOLX0(1)=0.0D0 MOLX0(2)=1.0D0-MOLX0(1) MOLX0(3)=0.0D0 MOLX1(1)=1.0D0 MOLX1(2)=1.0D0-MOLX1(1) MOLX1(3)=0.0D0 MOLX2(3)=0.0D0 MX0=MOLX0(1)-MOLFR(1) MX1=MOLX1(1)-MOLFR(1) DO 10 I=1,1000 P2=(P0+P1)*0.5D0 CALL FMF046(J,TEMP,P2,MOLX2,MOLY2,VL2,VV2) IF (J.NE.0) THEN CALL FMF062(J,'FMF037','IN FMF046-1') RETURN ENDIF MX2=MOLX2(1)-MOLFR(1) MFL=MX2*MX0 IF (MFL.GT.0.0D0) THEN MX0=MX2 MOLX0(1)=MOLX2(1) MOLX0(2)=MOLX2(2) P0=P2 ELSE MX1=MX2 MOLX1(1)=MOLX2(1) MOLX1(2)=MOLX2(2) P1=P2 ENDIF P2=P1-MX1*(P0-P1)/(MX0-MX1) CALL FMF046(J,TEMP,P2,MOLX2,MOLY2,VL2,VV2) IF (J.NE.0) THEN CALL FMF062(J,'FMF037','IN FMF046-2') RETURN ENDIF MX2=MOLX2(1)-MOLFR(1) PFL=((P1-P0)/P1)**2 IF (PFL.LT.EP2) THEN MFLG=MX2**2 IF (MFLG.LT.EP1) THEN J=0 VOL=VL2 PRESS=P2 RETURN ENDIF ENDIF MFL=MX2*MX0 IF (MFL.GT.0.0D0) THEN MX0=MX2 MOLX0(1)=MOLX2(1) MOLX0(2)=MOLX2(2) P0=P2 ELSE MX1=MX2 MOLX1(1)=MOLX2(1) MOLX1(2)=MOLX2(2) P1=P2 ENDIF 10 CONTINUE J=-1 RETURN END c --------------------------------------------------------- c Temperature,Mole fraction => Dew point pressure c --------------------------------------------------------- SUBROUTINE FMF038(J,TEMP,MOLFR,PRESS,VOL) DOUBLE PRECISION TEMP,MOLFR(1:3),PRESS $ ,EP1,EP2,VL,VV,PSAT1,PSAT2,P0,P1,P2 $ ,MOLY0(1:3),MOLY1(1:3),MOLY2(1:3),MY0,MY1,MY2,MFL,MFLG $ ,PFL,MOLX2(1:3),VL2,VV2,VOL INTEGER I,J DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C J=0 EP1=1.0D-14 EP2=1.0D-2 CALL FMF029(J,1,TEMP,PSAT1,VL,VV) IF (J.NE.0) THEN CALL FMF062(J,'FMF038','IN FMF029-1') RETURN ENDIF CALL FMF029(J,2,TEMP,PSAT2,VL,VV) IF (J.NE.0) THEN CALL FMF062(J,'FMF038','IN FMF029-2') RETURN ENDIF P0=PSAT2 P1=PSAT1 MOLY0(1)=0.0D0 MOLY0(2)=1.0D0-MOLY0(1) MOLY0(3)=0.0D0 MOLY1(1)=1.0D0 MOLY1(2)=1.0D0-MOLY1(1) MOLY1(3)=0.0D0 MOLY2(3)=0.0D0 MY0=MOLY0(1)-MOLFR(1) MY1=MOLY1(1)-MOLFR(1) DO 10 I=1,1000 P2=(P0+P1)*0.5D0 CALL FMF046(J,TEMP,P2,MOLX2,MOLY2,VL2,VV2) IF (J.NE.0) THEN CALL FMF062(J,'FMF038','IN FMF046-1') RETURN ENDIF MY2=MOLY2(1)-MOLFR(1) MFL=MY2*MY0 IF (MFL.GT.0.0D0) THEN MY0=MY2 MOLY0(1)=MOLY2(1) MOLY0(2)=MOLY2(2) P0=P2 ELSE MY1=MY2 MOLY1(1)=MOLY2(1) MOLY1(2)=MOLY2(2) P1=P2 ENDIF P2=P1-MY1*(P0-P1)/(MY0-MY1) CALL FMF046(J,TEMP,P2,MOLX2,MOLY2,VL2,VV2) IF (J.NE.0) THEN CALL FMF062(J,'FMF038','IN FMF046-2') RETURN ENDIF MY2=MOLY2(1)-MOLFR(1) PFL=((P1-P0)/P1)**2 IF (PFL.LT.EP2) THEN MFLG=MY2**2 IF (MFLG.LT.EP1) THEN J=0 VOL=VV2 PRESS=P2 RETURN ENDIF ENDIF MFL=MY2*MY0 IF (MFL.GT.0.0D0) THEN MY0=MY2 MOLY0(1)=MOLY2(1) MOLY0(2)=MOLY2(2) P0=P2 ELSE MY1=MY2 MOLY1(1)=MOLY2(1) MOLY1(2)=MOLY2(2) P1=P2 ENDIF 10 CONTINUE J=-1 RETURN END c ------------------------------------------------------------------ c Temperature,Pressure,Mole fraction => Volume,Enthalpy,Entropy c ------------------------------------------------------------------ SUBROUTINE FMF039(J,MOLFR,TEMP,PRESS,VOLM,ENTHAL,ENTROP $ ,MOLX,MOLY) DOUBLE PRECISION TEMP,PRESS,MOLFR(1:3) $ ,VOLM,ENTHAL,ENTROP,MOLX(1:3),MOLY(1:3),PBUB,PDEW,VOLP(1:3) $ ,VOL(1:2),VOLL,VOLV,HL,HV $ ,SL,SV,QUALTY,FMF020,FMF021,VL,VV INTEGER I,J DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C J=0 CALL FMF037(J,TEMP,MOLFR,PBUB,VL) IF (J.NE.0) THEN CALL FMF062(J,'FMF039','IN FMF037-1') RETURN ENDIF CALL FMF038(J,TEMP,MOLFR,PDEW,VV) IF (J.NE.0) THEN CALL FMF062(J,'FMF039','IN FMF038-1') RETURN ENDIF IF (PRESS.LE.PDEW) THEN DO 10 I=1,3 MOLX(I)=0.0D0 MOLY(I)=0.0D0 10 CONTINUE CALL FMF022(J,MOLFR,TEMP,VOLP) IF (J.NE.0) THEN CALL FMF062(J,'FMF039','IN FMF022-1') RETURN ENDIF CALL FMF024(J,MOLFR,TEMP,PRESS,VOLP(2),2,VOL) IF (J.NE.0) THEN CALL FMF062(J,'FMF039','IN FMF024-1') RETURN ENDIF VOLM=VOL(2) ENTHAL=FMF020(MOLFR,TEMP,VOL(2)) ENTROP=FMF021(MOLFR,TEMP,VOL(2)) RETURN ELSEIF ((PRESS.GT.PDEW).AND.(PRESS.LT.PBUB)) THEN CALL FMF046(J,TEMP,PRESS,MOLX,MOLY,VOLL,VOLV) IF (J.NE.0) THEN CALL FMF062(J,'FMF039','IN FMF046-1') RETURN ENDIF QUALTY=(MOLFR(1)-MOLX(1))/(MOLY(1)-MOLX(1)) VOLM=(1.0D0-QUALTY)*VOLL+QUALTY*VOLV HL=FMF020(MOLX,TEMP,VOLL) HV=FMF020(MOLY,TEMP,VOLV) SL=FMF021(MOLX,TEMP,VOLL) SV=FMF021(MOLY,TEMP,VOLV) ENTHAL=(1.0D0-QUALTY)*HL+QUALTY*HV ENTROP=(1.0D0-QUALTY)*SL+QUALTY*SV RETURN ELSE DO 20 I=1,3 MOLX(I)=0.0D0 MOLY(I)=0.0D0 20 CONTINUE CALL FMF024(J,MOLFR,TEMP,PRESS,1.0D0,1,VOL) IF (J.NE.0) THEN CALL FMF062(J,'FMF039','IN FMF024-2') RETURN ENDIF VOLM=VOL(1) ENTHAL=FMF020(MOLFR,TEMP,VOL(1)) ENTROP=FMF021(MOLFR,TEMP,VOL(1)) RETURN ENDIF END c --------------------------------------------------------------- c Mole fraction,Pressure,Enthalpy => Temperature,Volume,Entropy c --------------------------------------------------------------- SUBROUTINE FMF040(J,MOLFR,TEMP,PRESS,VOLM,ENTHAL,ENTROP $ ,MOLX,MOLY) INTEGER I,J DOUBLE PRECISION MOLFR(1:3),TEMP,PRESS $ ,VOLM,ENTHAL,ENTROP,MOLX(1:3),MOLY(1:3) $ ,TBUB,TDEW,VBUB,VDEW,HBUB,HDEW,T0,T1,T2,V1(1:2),V2(1:2) $ ,H0,H1,H2,VP(1:2),TE,HE,EP1,EP2,HH2,S2,HF,VM2 $ ,FMF020,FMF021 DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C EP1=1.0D-14 EP2=1.0D-14 DO 100 I=1,3 MOLX(I)=0.0D0 MOLY(I)=0.0D0 100 CONTINUE C CALL FMF043(J,PRESS,MOLFR,TBUB,VBUB) IF (J.NE.0) THEN CALL FMF062(J,'FMF040','IN FMF043') RETURN ENDIF CALL FMF044(J,PRESS,MOLFR,TDEW,VDEW) IF (J.NE.0) THEN CALL FMF062(J,'FMF040','IN FMF044') RETURN ENDIF C HBUB=FMF020(MOLFR,TBUB,VBUB) HDEW=FMF020(MOLFR,TDEW,VDEW) IF (ENTHAL.LE.HBUB) THEN H0=HBUB-ENTHAL T0=TBUB T1=TBUB*0.99D0 CALL FMF022(J,MOLFR,T1,VP) IF (J.NE.0) THEN CALL FMF062(J,'FMF040','IN FMF022-1') RETURN ENDIF CALL FMF024(J,MOLFR,T1,PRESS,VP(2),1,V1) IF (J.NE.0) THEN CALL FMF062(J,'FMF040','IN FMF024-1') RETURN ENDIF H1=FMF020(MOLFR,T1,V1(1))-ENTHAL DO 10 I=1,100 T2=T1-H1*(T1-T0)/(H1-H0) CALL FMF022(J,MOLFR,T2,VP) IF (J.NE.0) THEN CALL FMF062(J,'FMF040','IN FMF022-2') RETURN ENDIF CALL FMF024(J,MOLFR,T2,PRESS,VP(2),1,V2) IF (J.NE.0) THEN CALL FMF062(J,'FMF040','IN FMF024-2') RETURN ENDIF H2=FMF020(MOLFR,T2,V2(1))-ENTHAL TE=((T2-T1)/T2)**2 HE=H2**2 IF (TE.LT.EP1.OR.HE.LT.EP2) THEN J=0 TEMP=T2 VOLM=V2(1) ENTROP=FMF021(MOLFR,T2,V2(1)) RETURN ENDIF H0=H1 T0=T1 H1=H2 T1=T2 10 CONTINUE J=-1 RETURN ELSEIF (ENTHAL.GT.HBUB.AND.ENTHAL.LT.HDEW) THEN H0=HBUB-ENTHAL H1=HDEW-ENTHAL T0=TBUB T1=TDEW DO 20 I=1,100 T2=(T1+T2)*0.5D0 CALL FMF039(J,MOLFR,T2,PRESS,VM2,HH2,S2,MOLX,MOLY) IF (J.NE.0) THEN CALL FMF062(J,'FMF040','IN FMF039-1') RETURN ENDIF H2=HH2-ENTHAL HF=H2*H0 IF (HF.GT.0.0D0) THEN H0=H2 T0=T2 ELSE H1=H2 T1=T2 ENDIF T2=T1-H1*(T1-T0)/(H1-H0) CALL FMF039(J,MOLFR,T2,PRESS,VM2,HH2,S2,MOLX,MOLY) IF (J.NE.0) THEN CALL FMF062(J,'FMF040','IN FMF039-2') RETURN ENDIF H2=HH2-ENTHAL HF=H2*H0 TE=((T2-T1)/T2)**2 HE=H2**2 IF (TE.LT.EP1.OR.HE.LT.EP2) THEN J=0 TEMP=T2 VOLM=VM2 ENTROP=S2 RETURN ENDIF IF (HF.GT.0.0D0) THEN H0=H2 T0=T2 ELSE H1=H2 T1=T2 ENDIF 20 CONTINUE J=-1 RETURN ELSE H0=HDEW-ENTHAL T0=TDEW T1=TDEW*1.01D0 CALL FMF022(J,MOLFR,T1,VP) IF (J.NE.0) THEN CALL FMF062(J,'FMF040','IN FMF022-3') RETURN ENDIF CALL FMF024(J,MOLFR,T1,PRESS,VP(2),2,V1) IF (J.NE.0) THEN CALL FMF062(J,'FMF040','IN FMF024-3') RETURN ENDIF H1=FMF020(MOLFR,T1,V1(2))-ENTHAL DO 30 I=1,100 T2=T1-H1*(T1-T0)/(H1-H0) CALL FMF022(J,MOLFR,T2,VP) IF (J.NE.0) THEN CALL FMF062(J,'FMF040','IN FMF022-4') RETURN ENDIF CALL FMF024(J,MOLFR,T2,PRESS,VP(2),2,V2) IF (J.NE.0) THEN CALL FMF062(J,'FMF040','IN FMF024-4') RETURN ENDIF H2=FMF020(MOLFR,T2,V2(2))-ENTHAL TE=((T2-T1)/T2)**2 HE=H2**2 IF (TE.LT.EP1.OR.HE.LT.EP2) THEN J=0 TEMP=T2 VOLM=V2(2) ENTROP=FMF021(MOLFR,T2,V2(2)) RETURN ENDIF H0=H1 T0=T1 H1=H2 T1=T2 30 CONTINUE J=-1 RETURN ENDIF J=-2 RETURN END c --------------------------------------------------------------- c Mole fraction,Pressure,Enthalpy => Temperature,Volume,Entropy c --------------------------------------------------------------- SUBROUTINE FMF041(J,MOLFR,TEMP,PRESS,VOLM,ENTHAL,ENTROP $ ,MOLX,MOLY) INTEGER I,J DOUBLE PRECISION MOLFR(1:3),TEMP,PRESS $ ,VOLM,ENTHAL,ENTROP,MOLX(1:3),MOLY(1:3) $ ,TBUB,TDEW,VBUB,VDEW,SBUB,SDEW,T0,T1,T2,V1(1:2),V2(1:2) $ ,S0,S1,S2,VP(1:2),TE,SE,EP1,EP2,SS2,H2,SF,VM2 $ ,FMF020,FMF021 DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C EP1=1.0D-14 EP2=1.0D-14 DO 100 I=1,3 MOLX(I)=0.0D0 MOLY(I)=0.0D0 100 CONTINUE C CALL FMF043(J,PRESS,MOLFR,TBUB,VBUB) IF (J.NE.0) THEN CALL FMF062(J,'FMF041','IN FMF043') RETURN ENDIF CALL FMF044(J,PRESS,MOLFR,TDEW,VDEW) IF (J.NE.0) THEN CALL FMF062(J,'FMF041','IN FMF044') RETURN ENDIF C SBUB=FMF021(MOLFR,TBUB,VBUB) SDEW=FMF021(MOLFR,TDEW,VDEW) IF (ENTROP.LE.SBUB) THEN S0=SBUB-ENTROP T0=TBUB T1=TBUB*0.99D0 CALL FMF022(J,MOLFR,T1,VP) IF (J.NE.0) THEN CALL FMF062(J,'FMF041','IN FMF022-1') RETURN ENDIF CALL FMF024(J,MOLFR,T1,PRESS,VP(2),1,V1) IF (J.NE.0) THEN CALL FMF062(J,'FMF041','IN FMF024-1') RETURN ENDIF S1=FMF021(MOLFR,T1,V1(1))-ENTROP DO 10 I=1,100 T2=T1-S1*(T1-T0)/(S1-S0) CALL FMF022(J,MOLFR,T2,VP) IF (J.NE.0) THEN CALL FMF062(J,'FMF041','IN FMF022-2') RETURN ENDIF CALL FMF024(J,MOLFR,T2,PRESS,VP(2),1,V2) IF (J.NE.0) THEN CALL FMF062(J,'FMF041','IN FMF024-2') RETURN ENDIF S2=FMF021(MOLFR,T2,V2(1))-ENTROP TE=((T2-T1)/T2)**2 SE=S2**2 IF (TE.LT.EP1.OR.SE.LT.EP2) THEN J=0 TEMP=T2 VOLM=V2(1) ENTHAL=FMF020(MOLFR,T2,V2(1)) RETURN ENDIF S0=S1 T0=T1 S1=S2 T1=T2 10 CONTINUE J=-1 RETURN ELSEIF (ENTROP.GT.SBUB.AND.ENTROP.LT.SDEW) THEN S0=SBUB-ENTROP S1=SDEW-ENTROP T0=TBUB T1=TDEW DO 20 I=1,100 T2=(T1+T2)*0.5D0 CALL FMF039(J,MOLFR,T2,PRESS,VM2,H2,SS2,MOLX,MOLY) IF (J.NE.0) THEN CALL FMF062(J,'FMF041','IN FMF039-1') RETURN ENDIF S2=SS2-ENTROP SF=S2*S0 IF (SF.GT.0.0D0) THEN S0=S2 T0=T2 ELSE S1=S2 T1=T2 ENDIF T2=T1-S1*(T1-T0)/(S1-S0) CALL FMF039(J,MOLFR,T2,PRESS,VM2,H2,SS2,MOLX,MOLY) IF (J.NE.0) THEN CALL FMF062(J,'FMF041','IN FMF039-2') RETURN ENDIF S2=SS2-ENTROP SF=S2*S0 TE=((T2-T1)/T2)**2 SE=S2**2 IF (TE.LT.EP1.OR.SE.LT.EP2) THEN J=0 TEMP=T2 VOLM=VM2 ENTHAL=H2 RETURN ENDIF IF (SF.GT.0.0D0) THEN S0=S2 T0=T2 ELSE S1=S2 T1=T2 ENDIF 20 CONTINUE J=-1 RETURN ELSE S0=SDEW-ENTROP T0=TDEW T1=TDEW*1.01D0 CALL FMF022(J,MOLFR,T1,VP) IF (J.NE.0) THEN CALL FMF062(J,'FMF041','IN FMF022-3') RETURN ENDIF CALL FMF024(J,MOLFR,T1,PRESS,VP(2),2,V1) IF (J.NE.0) THEN CALL FMF062(J,'FMF041','IN FMF024-3') RETURN ENDIF S1=FMF021(MOLFR,T1,V1(2))-ENTROP DO 30 I=1,100 T2=T1-S1*(T1-T0)/(S1-S0) CALL FMF022(J,MOLFR,T2,VP) IF (J.NE.0) THEN CALL FMF062(J,'FMF041','IN FMF022-4') RETURN ENDIF CALL FMF024(J,MOLFR,T2,PRESS,VP(2),2,V2) IF (J.NE.0) THEN CALL FMF062(J,'FMF041','IN FMF024-4') RETURN ENDIF S2=FMF021(MOLFR,T2,V2(2))-ENTROP TE=((T2-T1)/T2)**2 SE=S2**2 IF (TE.LT.EP1.OR.SE.LT.EP2) THEN J=0 TEMP=T2 VOLM=V2(2) ENTHAL=FMF020(MOLFR,T2,V2(2)) RETURN ENDIF S0=S1 T0=T1 S1=S2 T1=T2 30 CONTINUE J=-1 RETURN ENDIF J=-2 RETURN END c --------------------------------------------------------------- c Mole fraction,Pressure,Enthalpy => Temperature,Volume,Entropy c --------------------------------------------------------------- SUBROUTINE FMF042(J,MOLFR,TEMP,PRESS,VOLM,ENTHAL,ENTROP $ ,MOLX,MOLY) INTEGER I,J DOUBLE PRECISION MOLFR(1:3),TEMP,PRESS $ ,VOLM,ENTHAL,ENTROP,MOLX(1:3),MOLY(1:3) $ ,TBUB,TDEW,VBUB,VDEW,V0,V1,V2,T0,T1,T2,VV1(1:2),VV2(1:2) $ ,VP(1:2),TE,VE,EP1,EP2,S2,H2,VF,VM2 $ ,FMF020,FMF021 DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C EP1=1.0D-14 EP2=1.0D-14 DO 100 I=1,3 MOLX(I)=0.0D0 MOLY(I)=0.0D0 100 CONTINUE C CALL FMF043(J,PRESS,MOLFR,TBUB,VBUB) IF (J.NE.0) THEN CALL FMF062(J,'FMF042','IN FMF043') RETURN ENDIF CALL FMF044(J,PRESS,MOLFR,TDEW,VDEW) IF (J.NE.0) THEN CALL FMF062(J,'FMF042','IN FMF044') RETURN ENDIF C IF (VOLM.LE.VBUB) THEN V0=VBUB-VOLM T0=TBUB T1=TBUB*0.99D0 CALL FMF022(J,MOLFR,T1,VP) IF (J.NE.0) THEN CALL FMF062(J,'FMF042','IN FMF022-1') RETURN ENDIF CALL FMF024(J,MOLFR,T1,PRESS,VP(2),1,VV1) IF (J.NE.0) THEN CALL FMF062(J,'FMF042','IN FMF024-1') RETURN ENDIF V1=VV1(1)-VOLM DO 10 I=1,100 T2=T1-V1*(T1-T0)/(V1-V0) CALL FMF022(J,MOLFR,T2,VP) IF (J.NE.0) THEN CALL FMF062(J,'FMF042','IN FMF022-2') RETURN ENDIF CALL FMF024(J,MOLFR,T2,PRESS,VP(2),1,VV2) IF (J.NE.0) THEN CALL FMF062(J,'FMF042','IN FMF024-2') RETURN ENDIF V2=VV2(1)-VOLM TE=((T2-T1)/T2)**2 VE=V2**2 IF (TE.LT.EP1.OR.VE.LT.EP2) THEN J=0 TEMP=T2 ENTHAL=FMF020(MOLFR,T2,VV2(1)) ENTROP=FMF021(MOLFR,T2,VV2(1)) RETURN ENDIF V0=V1 T0=T1 V1=V2 T1=T2 10 CONTINUE J=-1 RETURN ELSEIF (VOLM.GT.VBUB.AND.VOLM.LT.VDEW) THEN V0=VBUB-VOLM V1=VDEW-VOLM T0=TBUB T1=TDEW DO 20 I=1,100 T2=(T1+T2)*0.5D0 CALL FMF039(J,MOLFR,T2,PRESS,VM2,H2,S2,MOLX,MOLY) IF (J.NE.0) THEN CALL FMF062(J,'FMF042','IN FMF039-1') RETURN ENDIF V2=VM2-VOLM VF=V2*V0 IF (VF.GT.0.0D0) THEN V0=V2 T0=T2 ELSE V1=V2 T1=T2 ENDIF T2=T1-V1*(T1-T0)/(V1-V0) CALL FMF039(J,MOLFR,T2,PRESS,VM2,H2,S2,MOLX,MOLY) IF (J.NE.0) THEN CALL FMF062(J,'FMF042','IN FMF039-2') RETURN ENDIF V2=VM2-VOLM TE=((T2-T1)/T2)**2 VE=V2**2 IF (TE.LT.EP1.OR.VE.LT.EP2) THEN J=0 TEMP=T2 ENTHAL=H2 ENTROP=S2 RETURN ENDIF VF=V2*V0 IF (VF.GT.0.0D0) THEN V0=V2 T0=T2 ELSE V1=V2 T1=T2 ENDIF 20 CONTINUE J=-1 RETURN ELSE V0=VDEW-VOLM T0=TDEW T1=TDEW*1.01D0 CALL FMF022(J,MOLFR,T1,VP) IF (J.NE.0) THEN CALL FMF062(J,'FMF042','IN FMF022-3') RETURN ENDIF CALL FMF024(J,MOLFR,T1,PRESS,VP(2),2,VV1) IF (J.NE.0) THEN CALL FMF062(J,'FMF042','IN FMF024-3') RETURN ENDIF V1=VV1(2)-VOLM DO 30 I=1,100 T2=T1-V1*(T1-T0)/(V1-V0) CALL FMF022(J,MOLFR,T2,VP) IF (J.NE.0) THEN CALL FMF062(J,'FMF042','IN FMF022-4') RETURN ENDIF CALL FMF024(J,MOLFR,T2,PRESS,VP(2),2,VV2) IF (J.NE.0) THEN CALL FMF062(J,'FMF042','IN FMF024-4') RETURN ENDIF V2=VV2(2)-VOLM TE=((T2-T1)/T2)**2 VE=V2**2 IF (TE.LT.EP1.OR.VE.LT.EP2) THEN J=0 TEMP=T2 ENTHAL=FMF020(MOLFR,T2,VV2(2)) ENTROP=FMF021(MOLFR,T2,VV2(2)) RETURN ENDIF V0=V1 T0=T1 V1=V2 T1=T2 30 CONTINUE J=-1 RETURN ENDIF J=-2 RETURN END c --------------------------------------------------------- c Pressure,Mole fraction => Bubble point temperatuer c --------------------------------------------------------- SUBROUTINE FMF043(J,PRESS,MOLFR,TEMP,VOL) DOUBLE PRECISION TEMP,MOLFR(1:3),PRESS $ ,EP1,EP2,VL,VV,TSAT1,TSAT2,T0,T1,T2 $ ,MOLX0(1:3),MOLX1(1:3),MOLX2(1:3),MX0,MX1,MX2,MFL,MFLG $ ,TFL,MOLY2(1:3),VL2,VV2,VOL INTEGER I,J DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C J=0 EP1=1.0D-14 EP2=1.0D-14 CALL FMF036(J,1,PRESS,TSAT1,VL,VV) IF (J.NE.0) THEN CALL FMF062(J,'FMF043','IN FMF036-1') RETURN ENDIF CALL FMF036(J,2,PRESS,TSAT2,VL,VV) IF (J.NE.0) THEN CALL FMF062(J,'FMF043','IN FMF036-2') RETURN ENDIF T0=TSAT2 T1=TSAT1 MOLX0(1)=1.0D-10 MOLX0(2)=1.0D0-MOLX0(1) MOLX0(3)=0.0D0 MOLX1(1)=1.0D0-1.0D-10 MOLX1(2)=1.0D0-MOLX1(1) MOLX1(3)=0.0D0 MOLX2(3)=0.0D0 MX0=MOLX0(1)-MOLFR(1) MX1=MOLX1(1)-MOLFR(1) DO 10 I=1,1000 T2=(T0+T1)*0.5D0 CALL FMF046(J,T2,PRESS,MOLX2,MOLY2,VL2,VV2) IF (J.NE.0) THEN CALL FMF062(J,'FMF043','IN FMF046-1') RETURN ENDIF MX2=MOLX2(1)-MOLFR(1) MFL=MX2*MX0 IF (MFL.GT.0.0D0) THEN MX0=MX2 MOLX0(1)=MOLX2(1) MOLX0(2)=MOLX2(2) T0=T2 ELSE MX1=MX2 MOLX1(1)=MOLX2(1) MOLX1(2)=MOLX2(2) T1=T2 ENDIF T2=T1-MX1*(T0-T1)/(MX0-MX1) CALL FMF046(J,T2,PRESS,MOLX2,MOLY2,VL2,VV2) IF (J.NE.0) THEN CALL FMF062(J,'FMF043','IN FMF046-2') RETURN ENDIF MX2=MOLX2(1)-MOLFR(1) TFL=((T1-T0)/T1)**2 IF (TFL.LT.EP2) THEN MFLG=MX2**2 IF (MFLG.LT.EP1) THEN J=0 VOL=VL2 TEMP=T2 RETURN ENDIF ENDIF MFL=MX2*MX0 IF (MFL.GT.0.0D0) THEN MX0=MX2 MOLX0(1)=MOLX2(1) MOLX0(2)=MOLX2(2) T0=T2 ELSE MX1=MX2 MOLX1(1)=MOLX2(1) MOLX1(2)=MOLX2(2) T1=T2 ENDIF 10 CONTINUE J=-1 RETURN END c --------------------------------------------------------- c Pressure,Mole fraction => Dew point temperature c --------------------------------------------------------- SUBROUTINE FMF044(J,PRESS,MOLFR,TEMP,VOL) DOUBLE PRECISION PRESS,MOLFR(1:3),TEMP $ ,EP1,EP2,VL,VV,TSAT1,TSAT2,T0,T1,T2 $ ,MOLY0(1:3),MOLY1(1:3),MOLY2(1:3),MY0,MY1,MY2,MFL,MFLG $ ,TFL,MOLX2(1:3),VOL,VL2,VV2 INTEGER I,J DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C J=0 EP1=1.0D-14 EP2=1.0D-14 CALL FMF036(J,1,PRESS,TSAT1,VL,VV) IF (J.NE.0) THEN CALL FMF062(J,'FMF044','IN FMF036-1') RETURN ENDIF CALL FMF036(J,2,PRESS,TSAT2,VL,VV) IF (J.NE.0) THEN CALL FMF062(J,'FMF044','IN FMF036-2') RETURN ENDIF T0=TSAT2 T1=TSAT1 MOLY0(1)=0.0D0 MOLY0(2)=1.0D0-MOLY0(1) MOLY0(3)=0.0D0 MOLY1(1)=1.0D0 MOLY1(2)=1.0D0-MOLY1(1) MOLY1(3)=0.0D0 MOLY2(3)=0.0D0 MY0=MOLY0(1)-MOLFR(1) MY1=MOLY1(1)-MOLFR(1) DO 10 I=1,1000 T2=(T0+T1)*0.5D0 CALL FMF046(J,T2,PRESS,MOLX2,MOLY2,VL2,VV2) IF (J.NE.0) THEN CALL FMF062(J,'FMF044','IN FMF046-1') RETURN ENDIF MY2=MOLY2(1)-MOLFR(1) MFL=MY2*MY0 IF (MFL.GT.0.0D0) THEN MY0=MY2 MOLY0(1)=MOLY2(1) MOLY0(2)=MOLY2(2) T0=T2 ELSE MY1=MY2 MOLY1(1)=MOLY2(1) MOLY1(2)=MOLY2(2) T1=T2 ENDIF T2=T1-MY1*(T0-T1)/(MY0-MY1) CALL FMF046(J,T2,PRESS,MOLX2,MOLY2,VL2,VV2) IF (J.NE.0) THEN CALL FMF062(J,'FMF044','IN FMF046-2') RETURN ENDIF MY2=MOLY2(1)-MOLFR(1) TFL=((T1-T0)/T1)**2 IF (TFL.LT.EP2) THEN MFLG=MY2**2 IF (MFLG.LT.EP1) THEN J=0 VOL=VV2 TEMP=T2 RETURN ENDIF ENDIF MFL=MY2*MY0 IF (MFL.GT.0.0D0) THEN MY0=MY2 MOLY0(1)=MOLY2(1) MOLY0(2)=MOLY2(2) T0=T2 ELSE MY1=MY2 MOLY1(1)=MOLY2(1) MOLY1(2)=MOLY2(2) T1=T2 ENDIF 10 CONTINUE J=-1 RETURN END c ------------------------------------------------------ c Temperature,Volume,cmpnum => (d2p/dv2)_{T} c ------------------------------------------------------ DOUBLE PRECISION FUNCTION FMF045(CMPNUM,T,V) DOUBLE PRECISION T,V,A,B,R,VB,VB4 $ ,FMF001,FMF002 INTEGER CMPNUM DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C A=FMF001(CMPNUM,T) B=FMF002(CMPNUM,T) R=CONST(61) VB=V+B VB4=B-4.0D0*V IF (V.LE.1.0D-15.OR.VB.EQ.0.0D0) THEN CALL FMF062(-2,'FMF045','V = 0 OR VB = 0') FMF045=-1.0D20 RETURN ELSEIF (VB4.EQ.0.0D0) THEN CALL FMF062(-2,'FMF045','VB4 = 0') FMF045=-1.0D20 RETURN ENDIF C FMF045=-1.0D0*(2.0D0*(A*B**7-17.0D0*A*B**6*V+103.0D0*A*B**5*V**2 $ -220.0D0*A*B**4*V**3-160.0D0*A*B**3*V**4+896.0D0*A*B**2*V**5 $ +768.0D0*A*B*V**6-3072.0D0*A*V**7-B**8*R*T+17.0D0*B**7 $ *R*T*V-103.0D0*B**6*R*T*V**2+219.0D0*B**5*R*T*V**3 $ +3252.0D0*B**4*R*T*V**4+8160.0D0*B**3*R*T*V**5+9088.0D0 $ *B**2*R*T*V**6+4864.0D0*B*R*T*V**7+1024.0D0*R*T*V**8)) $ /(VB**3*VB4**5*V**3) RETURN END c --------------------------------------------------------------- c Pressure,Temperature => x,y (use K1,K2) c --------------------------------------------------------------- SUBROUTINE FMF046(J,TEMP,PRESS,MOLX,MOLY,VOLL,VOLV) DOUBLE PRECISION TEMP,PRESS,MOLX(1:3) $ ,MOLY(1:3),VOLL,VOLV,EP1,MOLX0(1:3),MOLY0(1:3),K10,K20 $ ,VPOLX(1:2),VPOLY(1:2),VX(1:2),VY(1:2),MUL1,MUL2,MUV1,MUV2 $ ,FL1,FL2,FV1,FV2,FLPX1,FLPX2,FVPY1,FVPY2,K1,K2,K1F,K2F $ ,FMF010 INTEGER I,J DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C IF (TEMP.LE.0.0D0) THEN J=-2 CALL FMF062(J,'FMF046','TEMPERATURE IS LESS THAN 0.') RETURN ENDIF J=0 EP1=1.0D-14 MOLX0(1)=1.0D-10 MOLX0(2)=1.0D0-MOLX0(1) MOLX0(3)=0.0D0 MOLY0(1)=1.0D0-1.0D-10 MOLY0(2)=1.0D0-MOLY0(1) MOLY0(3)=0.0D0 K10=MOLY0(1)/MOLX0(1) K20=MOLY0(2)/MOLX0(2) DO 10 I=1,100 CALL FMF024(J,MOLX0,TEMP,PRESS,VPOLX(2),1,VX) IF (J.NE.0) THEN CALL FMF062(J,'FMF046','IN FMF024-1') RETURN ENDIF CALL FMF022(J,MOLY0,TEMP,VPOLY) IF (J.EQ.-2) THEN CALL FMF062(J,'FMF046','IN FMF022-1') RETURN ELSEIF (J.EQ.-1) THEN CALL FMF023(J,MOLY0,TEMP,VPOLY(2)) IF (J.EQ.-2) THEN CALL FMF062(J,'FMF046','FMF023-1') RETURN ENDIF J=0 ENDIF CALL FMF024(J,MOLY0,TEMP,PRESS,VPOLY(2),2,VY) IF (J.NE.0) THEN CALL FMF062(J,'FMF046','IN FMF024-2') RETURN ENDIF MUL1=FMF010(MOLX0,TEMP,VX(1),1) IF (MUL1.LE.-1.0D20) THEN J=-2 CALL FMF062(J,'FMF046','IN FMF010-1') RETURN ENDIF MUL2=FMF010(MOLX0,TEMP,VX(1),2) IF (MUL2.LE.-1.0D20) THEN J=-2 CALL FMF062(J,'FMF046','IN FMF010-2') RETURN ENDIF MUV1=FMF010(MOLY0,TEMP,VY(2),1) IF (MUV1.LE.-1.0D20) THEN J=-2 CALL FMF062(J,'FMF046','IN FMF010-3') RETURN ENDIF MUV2=FMF010(MOLY0,TEMP,VY(2),2) IF (MUV2.LE.-1.0D20) THEN J=-2 CALL FMF062(J,'FMF046','IN FMF010-4') RETURN ENDIF FL1=DEXP(MUL1/(CONST(61)*TEMP)) FL2=DEXP(MUL2/(CONST(61)*TEMP)) FV1=DEXP(MUV1/(CONST(61)*TEMP)) FV2=DEXP(MUV2/(CONST(61)*TEMP)) FLPX1=FL1/MOLX0(1) FLPX2=FL2/MOLX0(2) FVPY1=FV1/MOLY0(1) FVPY2=FV2/MOLY0(2) K1=FLPX1/FVPY1 K2=FLPX2/FVPY2 K1F=((K1-K10)/K1)**2 K2F=((K2-K20)/K2)**2 IF (K1F.LT.EP1.AND.K2F.LT.EP1) THEN J=0 MOLX(1)=MOLX0(1) MOLX(2)=MOLX0(2) MOLX(3)=MOLX0(3) MOLY(1)=MOLY0(1) MOLY(2)=MOLY0(2) MOLY(3)=MOLY0(3) VOLL=VX(1) VOLV=VY(2) RETURN ENDIF MOLX0(1)=(K2-1.0D0)/(K2-K1) MOLX0(2)=1.0D0-MOLX0(1) MOLY0(1)=K1*MOLX0(1) MOLY0(2)=1.0D0-MOLY0(1) K10=K1 K20=K2 10 CONTINUE J=-1 CALL FMF062(J,'FMF046','NO CONVERGENCE') RETURN END C --------------------------------------------- c This program calcurate the coefficients c to use in Tsat. c Psat=a_{i,0}+a_{i,1}*T+a_{i,2}*T**2 c --------------------------------------------- SUBROUTINE FMF047(J,CFEOT) DOUBLE PRECISION CFEOT(1:6) $ ,P1(1:2),P2(1:2),RC1(1:2) $ ,RC2(1:2),TEMP1(0:5),TEMP2(0:5),PSAT1(0:5),PSAT2(0:5) $ ,VL,VV,DT1,DT2,CFSAT1(1:3),CFSAT2(1:3) $ ,FMF025 INTEGER I,J DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C P1(1)=CONST(2) P1(2)=CONST(4) P2(1)=CONST(22) P2(2)=CONST(24) CALL FMF054(J,1,P1,RC1) IF (J.NE.0) THEN CALL FMF062(J,'FMF047','IN FMF054-1') RETURN ENDIF CALL FMF054(J,2,P2,RC2) IF (J.NE.0) THEN CALL FMF062(J,'FMF047','IN FMF054-2') RETURN ENDIF TEMP1(0)=RC1(1) PSAT1(0)=FMF025(1,RC1(1),RC1(2)) TEMP2(0)=RC2(1) PSAT2(0)=FMF025(2,RC2(1),RC2(2)) DT1=(TEMP1(0)-CONST(11))*0.2D0 DT2=(TEMP2(0)-CONST(31))*0.2D0 DO 10 I=1,5 TEMP1(I)=CONST(11)+DT1*DBLE(I-1) TEMP2(I)=CONST(31)+DT2*DBLE(I-1) CALL FMF029(J,1,TEMP1(I),PSAT1(I),VL,VV) IF (J.NE.0) THEN CALL FMF062(J,'FMF047','IN FMF029-1') RETURN ENDIF CALL FMF029(J,2,TEMP2(I),PSAT2(I),VL,VV) IF (J.NE.0) THEN CALL FMF062(J,'FMF047','IN FMF029-2') RETURN ENDIF 10 CONTINUE CALL FMF056(J,TEMP1,PSAT1,CFSAT1) IF (J.NE.0) THEN CALL FMF062(J,'FMF047','IN FMF056-1') RETURN ENDIF CALL FMF056(J,TEMP2,PSAT2,CFSAT2) IF (J.NE.0) THEN CALL FMF062(J,'FMF047','IN FMF056-2') RETURN ENDIF DO 20 I=1,3 CFEOT(I)=CFSAT1(I) CFEOT(I+3)=CFSAT2(I) 20 CONTINUE J=0 RETURN END C --------------------------------------------------------- c This program calcurate the value to c calibrate entalpy and entropy when KSTAN=1. c --------------------------------------------------------- SUBROUTINE FMF048(J,STAND) DOUBLE PRECISION STAND(1:4),PSAT,VL,VV,ENTHAL,ENTROP $ ,TEMP0,H0,S0 $ ,FMF030,FMF031 INTEGER I,J DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C TEMP0=273.15 H0=2.0D5 S0=1.0D3 DO 10 I=1,2 CALL FMF029(J,I,TEMP0,PSAT,VL,VV) IF (J.NE.0) THEN CALL FMF062(J,'FMF048','IN FMF029-1') RETURN ENDIF ENTHAL=FMF030(I,TEMP0,VL) ENTROP=FMF031(I,TEMP0,VL) STAND(2*I-1)=-1.0D3*ENTHAL+H0*CONST(20*I-19) STAND(2*I)=-1.0D3*ENTROP+S0*CONST(20*I-19) 10 CONTINUE RETURN END C ---------------------------------------- C Error message C ---------------------------------------- SUBROUTINE FMF049(J,PRNAME) INTEGER KPA,MESS,KSTAN,KAS,J COMMON /UNIT/ KPA,MESS,KSTAN,KAS CHARACTER*(*) PRNAME C IF (MESS.EQ.1) THEN IF (J.EQ.0) THEN RETURN ELSEIF (J.EQ.-1) THEN WRITE(*,1) PRNAME RETURN ELSEIF (J.EQ.-2) THEN WRITE(*,2) PRNAME RETURN ELSE RETURN ENDIF ENDIF 1 FORMAT (1X,'***** ERROR IN',1X,A6,1X,', NO CONVERGENCE *****') 2 FORMAT (1X,'***** ERROR IN',1X,A6,1X,', OUT OF RANGE *****') RETURN END c ---------------------------------------------------------- c Indicating the Region of Mixture c -2 : Out of range c -1 : No convergence c 1 : Liquid region c 2 : Moist region c 3 : Gas region c ---------------------------------------------------------- INTEGER FUNCTION FMF050(TEMP,PRESS,MOLFR) INTEGER J DOUBLE PRECISION TEMP,PRESS,VBUB,PBUB,PDEW,VDEW,MOLFR(1:3) DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C CALL FMF037(J,TEMP,MOLFR,PBUB,VBUB) IF (J.EQ.-2) THEN CALL FMF062(J,'FMF050','OUT OF RANGE(IN FMF037).') FMF050=-2 RETURN ELSEIF (J.EQ.-1) THEN CALL FMF062(J,'FMF050','NO CONVERGENCE(IN FMF037.)') FMF050=-1 RETURN ENDIF CALL FMF038(J,TEMP,MOLFR,PDEW,VDEW) IF (J.EQ.-2) THEN CALL FMF062(J,'FMF050','OUT OF RANGE(IN FMF038).') FMF050=-2 RETURN ELSEIF (J.EQ.-1) THEN CALL FMF062(J,'FMF050','NO CONVERGENCE(IN FMF038).') FMF050=-1 RETURN ENDIF IF (PRESS.GE.PBUB) THEN FMF050=1 RETURN ELSEIF ((PRESS.LT.PBUB).AND.(PRESS.GT.PDEW)) THEN FMF050=2 RETURN ELSEIF (PRESS.LE.PDEW) THEN FMF050=3 RETURN ENDIF END c --------------------------------------------------------------- c Pressure,Enthalpy => Temperature,Volume,Entropy c --------------------------------------------------------------- SUBROUTINE FMF051(J,CMPNUM,TEMP,PRESS,VOLM,ENTHAL,ENTROP) INTEGER I,J,CMPNUM DOUBLE PRECISION TEMP,PRESS,VOLM,ENTHAL $ ,ENTROP,TSAT,VL,VV,HL,HV,T0,T1,T2,V1(1:2),V2(1:2) $ ,H0,H1,H2,VP(1:2),TE,HE,EP1,EP2,QUAL,SL,SV $ ,FMF030,FMF031 DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C EP1=1.0D-14 EP2=1.0D-14 C CALL FMF036(J,CMPNUM,PRESS,TSAT,VL,VV) IF (J.NE.0) THEN CALL FMF062(J,'FMF051','IN FMF036-1') RETURN ENDIF C HL=FMF030(CMPNUM,TSAT,VL) HV=FMF030(CMPNUM,TSAT,VV) IF (ENTHAL.LE.HL) THEN H0=HL-ENTHAL T0=TSAT T1=TSAT*0.99D0 CALL FMF027(J,CMPNUM,T1,VP) IF (J.NE.0) THEN CALL FMF062(J,'FMF051','IN FMF027-1') RETURN ENDIF CALL FMF028(J,CMPNUM,T1,PRESS,VP(2),1,V1) IF (J.NE.0) THEN CALL FMF062(J,'FMF051','IN FMF028-1') RETURN ENDIF H1=FMF030(CMPNUM,T1,V1(1))-ENTHAL DO 10 I=1,100 T2=T1-H1*(T1-T0)/(H1-H0) CALL FMF027(J,CMPNUM,T2,VP) IF (J.NE.0) THEN CALL FMF062(J,'FMF051','IN FMF027-2') RETURN ENDIF CALL FMF028(J,CMPNUM,T2,PRESS,VP(2),1,V2) IF (J.NE.0) THEN CALL FMF062(J,'FMF051','IN FMF028-2') RETURN ENDIF H2=FMF030(CMPNUM,T2,V2(1))-ENTHAL TE=((T2-T1)/T2)**2 HE=H2**2 IF (TE.LT.EP1.OR.HE.LT.EP2) THEN J=0 TEMP=T2 VOLM=V2(1) ENTROP=FMF031(CMPNUM,T2,V2(1)) RETURN ENDIF H0=H1 T0=T1 H1=H2 T1=T2 10 CONTINUE J=-1 RETURN ELSEIF (ENTHAL.GT.HL.AND.ENTHAL.LT.HV) THEN QUAL=(ENTHAL-HL)/(HV-HL) SL=FMF031(CMPNUM,TSAT,VL) SV=FMF031(CMPNUM,TSAT,VV) ENTROP=QUAL*SV+(1.0D0-QUAL)*SL TEMP=TSAT VOLM=QUAL*VV+(1.0D0-QUAL)*VL RETURN ELSE H0=HV-ENTHAL T0=TSAT T1=TSAT*1.01D0 CALL FMF027(J,CMPNUM,T1,VP) IF (J.NE.0) THEN CALL FMF062(J,'FMF051','IN FMF027-3') RETURN ENDIF CALL FMF028(J,CMPNUM,T1,PRESS,VP(2),2,V1) IF (J.NE.0) THEN CALL FMF062(J,'FMF051','IN FMF028-3') RETURN ENDIF H1=FMF030(CMPNUM,T1,V1(2))-ENTHAL DO 30 I=1,100 T2=T1-H1*(T1-T0)/(H1-H0) CALL FMF027(J,CMPNUM,T2,VP) IF (J.NE.0) THEN CALL FMF062(J,'FMF051','IN FMF027-4') RETURN ENDIF CALL FMF028(J,CMPNUM,T2,PRESS,VP(2),2,V2) IF (J.NE.0) THEN CALL FMF062(J,'FMF051','IN FMF028-4') RETURN ENDIF H2=FMF030(CMPNUM,T2,V2(2))-ENTHAL TE=((T2-T1)/T2)**2 HE=H2**2 IF (TE.LT.EP1.OR.HE.LT.EP2) THEN J=0 TEMP=T2 VOLM=V2(2) ENTROP=FMF031(CMPNUM,T2,V2(2)) RETURN ENDIF H0=H1 T0=T1 H1=H2 T1=T2 30 CONTINUE J=-1 RETURN ENDIF J=-2 RETURN END c --------------------------------------------------------------- c Pressure,Entropy => Temperature,Volume,Enthalpy c --------------------------------------------------------------- SUBROUTINE FMF052(J,CMPNUM,TEMP,PRESS,VOLM,ENTHAL,ENTROP) INTEGER I,J,CMPNUM DOUBLE PRECISION TEMP,PRESS,VOLM,ENTHAL $ ,ENTROP,TSAT,VL,VV,HL,HV,T0,T1,T2,V1(1:2),V2(1:2) $ ,S0,S1,S2,VP(1:2),TE,SE,EP1,EP2,QUAL,SL,SV $ ,FMF030,FMF031 DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C EP1=1.0D-14 EP2=1.0D-14 C CALL FMF036(J,CMPNUM,PRESS,TSAT,VL,VV) IF (J.NE.0) THEN CALL FMF062(J,'FMF052','IN FMF036') RETURN ENDIF C SL=FMF031(CMPNUM,TSAT,VL) SV=FMF031(CMPNUM,TSAT,VV) IF (ENTROP.LE.SL) THEN S0=SL-ENTROP T0=TSAT T1=TSAT*0.99D0 CALL FMF027(J,CMPNUM,T1,VP) IF (J.NE.0) THEN CALL FMF062(J,'FMF052','IN FMF027-1') RETURN ENDIF CALL FMF028(J,CMPNUM,T1,PRESS,VP(2),1,V1) IF (J.NE.0) THEN CALL FMF062(J,'FMF052','IN FMF028-1') RETURN ENDIF S1=FMF031(CMPNUM,T1,V1(1))-ENTROP DO 10 I=1,100 T2=T1-S1*(T1-T0)/(S1-S0) CALL FMF027(J,CMPNUM,T2,VP) IF (J.NE.0) THEN CALL FMF062(J,'FMF052','IN FMF027-2') RETURN ENDIF CALL FMF028(J,CMPNUM,T2,PRESS,VP(2),1,V2) IF (J.NE.0) THEN CALL FMF062(J,'FMF052','IN FMF028-2') RETURN ENDIF S2=FMF031(CMPNUM,T2,V2(1))-ENTROP TE=((T2-T1)/T2)**2 SE=S2**2 IF (TE.LT.EP1.OR.SE.LT.EP2) THEN J=0 TEMP=T2 VOLM=V2(1) ENTHAL=FMF030(CMPNUM,T2,V2(1)) RETURN ENDIF S0=S1 T0=T1 S1=S2 T1=T2 10 CONTINUE J=-1 RETURN ELSEIF (ENTROP.GT.SL.AND.ENTROP.LT.SV) THEN QUAL=(ENTROP-SL)/(SV-SL) HL=FMF030(CMPNUM,TSAT,VL) HV=FMF030(CMPNUM,TSAT,VV) ENTHAL=QUAL*HV+(1.0D0-QUAL)*HL TEMP=TSAT VOLM=QUAL*VV+(1.0D0-QUAL)*VL RETURN ELSE S0=SV-ENTROP T0=TSAT T1=TSAT*1.01D0 CALL FMF027(J,CMPNUM,T1,VP) IF (J.NE.0) THEN CALL FMF062(J,'FMF052','IN FMF027-3') RETURN ENDIF CALL FMF028(J,CMPNUM,T1,PRESS,VP(2),2,V1) IF (J.NE.0) THEN CALL FMF062(J,'FMF052','IN FMF028-3') RETURN ENDIF S1=FMF031(CMPNUM,T1,V1(2))-ENTROP DO 30 I=1,100 T2=T1-S1*(T1-T0)/(S1-S0) CALL FMF027(J,CMPNUM,T2,VP) IF (J.NE.0) THEN CALL FMF062(J,'FMF052','IN FMF027-4') RETURN ENDIF CALL FMF028(J,CMPNUM,T2,PRESS,VP(2),2,V2) IF (J.NE.0) THEN CALL FMF062(J,'FMF052','IN FMF028-4') RETURN ENDIF S2=FMF031(CMPNUM,T2,V2(2))-ENTROP TE=((T2-T1)/T2)**2 SE=S2**2 IF (TE.LT.EP1.OR.SE.LT.EP2) THEN J=0 TEMP=T2 VOLM=V2(2) ENTHAL=FMF030(CMPNUM,T2,V2(2)) RETURN ENDIF S0=S1 T0=T1 S1=S2 T1=T2 30 CONTINUE J=-1 RETURN ENDIF J=-2 RETURN END c --------------------------------------------------------------- c Pressure,Volume => Temperature,Enthalpy,Entropy c --------------------------------------------------------------- SUBROUTINE FMF053(J,CMPNUM,TEMP,PRESS,VOLM,ENTHAL,ENTROP) INTEGER I,J,CMPNUM DOUBLE PRECISION TEMP,PRESS,VOLM,ENTHAL $ ,ENTROP,TSAT,VL,VV,HL,HV,T0,T1,T2,VV1(1:2),VV2(1:2) $ ,V0,V1,V2,VP(1:2),TE,VE,EP1,EP2,QUAL,SL,SV $ ,FMF030,FMF031 DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C EP1=1.0D-14 EP2=1.0D-14 C CALL FMF036(J,CMPNUM,PRESS,TSAT,VL,VV) IF (J.NE.0) THEN CALL FMF062(J,'FMF053','IN FMF036') RETURN ENDIF C IF (VOLM.LE.VL) THEN V0=VL-VOLM T0=TSAT T1=TSAT*0.99D0 CALL FMF027(J,CMPNUM,T1,VP) IF (J.NE.0) THEN CALL FMF062(J,'FMF053','IN FMF027-1') RETURN ENDIF CALL FMF028(J,CMPNUM,T1,PRESS,VP(2),1,VV1) IF (J.NE.0) THEN CALL FMF062(J,'FMF053','IN FMF028-1') RETURN ENDIF V1=VV1(1)-VOLM DO 10 I=1,100 T2=T1-V1*(T1-T0)/(V1-V0) CALL FMF027(J,CMPNUM,T2,VP) IF (J.NE.0) THEN CALL FMF062(J,'FMF053','IN FMF027-2') RETURN ENDIF CALL FMF028(J,CMPNUM,T2,PRESS,VP(2),1,VV2) IF (J.NE.0) THEN CALL FMF062(J,'FMF053','IN FMF028-2') RETURN ENDIF V2=VV2(1)-VOLM TE=((T2-T1)/T2)**2 VE=V2**2 IF (TE.LT.EP1.OR.VE.LT.EP2) THEN J=0 TEMP=T2 ENTHAL=FMF030(CMPNUM,T2,VV2(1)) ENTROP=FMF031(CMPNUM,T2,VV2(1)) RETURN ENDIF V0=V1 T0=T1 V1=V2 T1=T2 10 CONTINUE J=-1 RETURN ELSEIF (VOLM.GT.VL.AND.VOLM.LT.VV) THEN QUAL=(VOLM-VL)/(VV-VL) HL=FMF030(CMPNUM,TSAT,VL) HV=FMF030(CMPNUM,TSAT,VV) SL=FMF031(CMPNUM,TSAT,VL) SV=FMF031(CMPNUM,TSAT,VV) ENTHAL=QUAL*HV+(1.0D0-QUAL)*HL ENTROP=QUAL*SV+(1.0D0-QUAL)*SL TEMP=TSAT RETURN ELSE V0=VV-VOLM T0=TSAT T1=TSAT*1.01D0 CALL FMF027(J,CMPNUM,T1,VP) IF (J.NE.0) THEN CALL FMF062(J,'FMF053','IN FMF027-3') RETURN ENDIF CALL FMF028(J,CMPNUM,T1,PRESS,VP(2),2,V1) IF (J.NE.0) THEN CALL FMF062(J,'FMF053','IN FMF028-3') RETURN ENDIF V1=VV1(2)-VOLM DO 30 I=1,100 T2=T1-V1*(T1-T0)/(V1-V0) CALL FMF027(J,CMPNUM,T2,VP) IF (J.NE.0) THEN CALL FMF062(J,'FMF053','IN FMF027-4') RETURN ENDIF CALL FMF028(J,CMPNUM,T2,PRESS,VP(2),2,VV2) IF (J.NE.0) THEN CALL FMF062(J,'FMF053','IN FMF028-4') RETURN ENDIF V2=VV2(2)-VOLM TE=((T2-T1)/T2)**2 VE=V2**2 IF (TE.LT.EP1.OR.VE.LT.EP2) THEN J=0 TEMP=T2 ENTHAL=FMF030(CMPNUM,T2,VV2(2)) ENTROP=FMF031(CMPNUM,T2,VV2(2)) RETURN ENDIF V0=V1 T0=T1 V1=V2 T1=T2 30 CONTINUE J=-1 RETURN ENDIF J=-2 RETURN END c----------------------------------------------------------- c P:initial data c P2:result c DOM:domain of section c EP:condition of convergence c J:error code c----------------------------------------------------------- SUBROUTINE FMF054(J,CMPNUM,P,P2) INTEGER I,IJ,J,JH,JL,LIM,TIME,M,LAMBDA,MU,IR,CMPNUM PARAMETER (M=1664501,LAMBDA=1229,MU=351750) DOUBLE PRECISION DOM(1:2),P(1:2),X(1:4,1:2),R(1:6),P1(1:2) $ ,G(1:2),P2(1:2),XN(1:2),Y(1:4) $ ,ALFA,ALFA2,BETA,EP,YH,YL,YG,YN,INVM $ ,FMF055 DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES c EP=1.0D-14 DOM(1)=1.0D-2 DOM(2)=1.0D-4 ALFA=1.3D0 BETA=0.5D0 IR=0 INVM=1.0D0/DBLE(M) LIM=1000 TIME=0 c Y(1)=FMF055(CMPNUM,P) 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 I=2,4 X(I,1)=P(1)+DOM(1)*(R(2*I-3)-0.5D0) X(I,2)=P(2)+DOM(2)*(R(2*I-2)-0.5D0) 20 CONTINUE C DO 40 I=2,4 P1(1)=X(I,1) P1(2)=X(I,2) Y(I)=FMF055(CMPNUM,P1) 40 CONTINUE C 1000 YH=Y(1) JH=1 YL=Y(1) JL=1 DO 60 I=2,4 IF (Y(I).GT.YH) THEN YH=Y(I) JH=I ELSEIF (Y(I).LT.YL) THEN YL=Y(I) JL=I ENDIF 60 CONTINUE C TIME=TIME+1 C WRITE(*,*) ' TIME=',TIME IF (TIME.GT.LIM) THEN J=-1 RETURN ENDIF IF (YL.LT.EP) THEN P2(1)=X(JL,1) P2(2)=X(JL,2) RETURN ENDIF C G(1)=0.0D0 G(2)=0.0D0 DO 80 I=1,2 DO 90 IJ=1,4 G(I)=G(I)+X(IJ,I) 90 CONTINUE 80 CONTINUE G(1)=(G(1)-X(JH,1))/3.0D0 G(2)=(G(2)-X(JH,2))/3.0D0 YG=FMF055(CMPNUM,G) C ALFA2=ALFA IF (YG.LT.YH) THEN 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=FMF055(CMPNUM,XN) 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=FMF055(CMPNUM,XN) 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 RETURN END C -------------------------------------------------------- c This program calcurates the critical point c temperature and volume of pure component. c -------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF055(CMPNUM,VAL) INTEGER CMPNUM DOUBLE PRECISION VAL(1:2),FMF026,FMF045 DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES C FMF055=(FMF026(CMPNUM,VAL(1),VAL(2)))**2 $ +(FMF045(CMPNUM,VAL(1),VAL(2)))**2 RETURN END c --------------------------------------------------------- c Calcurating approximate equation from saturated data. c --------------------------------------------------------- SUBROUTINE FMF056(J,TEMP,PSAT,CFSAT) DOUBLE PRECISION TEMP(0:5),PSAT(0:5),CFSAT(1:3) $ ,COEFAA,COEFBA,COEFCA,COEFBB,COEFCB,CFINV,CFDEN INTEGER I,J c COEFAA=0.0D0 COEFBA=0.0D0 COEFCA=0.0D0 COEFBB=0.0D0 COEFCB=0.0D0 DO 10 I=1,5 COEFAA=COEFAA+TEMP(0)**4-2.0D0*TEMP(0)**2*TEMP(I)**2 $ +TEMP(I)**4 COEFBA=COEFBA+TEMP(0)**3-TEMP(0)**2*TEMP(I)-TEMP(0) $ *TEMP(I)**2+TEMP(I)**3 COEFCA=COEFCA-TEMP(0)**2*PSAT(0)+TEMP(0)**2*PSAT(I) $ +TEMP(I)**2*PSAT(0)-TEMP(I)**2*PSAT(I) COEFBB=COEFBB+TEMP(0)**2-2.0D0*TEMP(0)*TEMP(I) $ +TEMP(I)**2 COEFCB=COEFCB-TEMP(0)*PSAT(0)+TEMP(0)*PSAT(I) $ +TEMP(I)*PSAT(0)-TEMP(I)*PSAT(I) 10 CONTINUE CFDEN=(COEFBA**2-COEFAA*COEFBB) IF (CFDEN.EQ.0.0D0) THEN J=-2 CALL FMF062(J,'FMF056','DEVIDE BY 0') RETURN ENDIF CFINV=1.0D0/CFDEN CFSAT(3)=(COEFBB*COEFCA-COEFBA*COEFCB)*CFINV CFSAT(2)=(COEFAA*COEFCB-COEFBA*COEFCA)*CFINV CFSAT(1)=PSAT(0)-CFSAT(3)*TEMP(0)**2-CFSAT(2)*TEMP(0) J=0 RETURN END c ----------------------------------------------- c This program transfers the unit of fraction c from kmol to kg. c ----------------------------------------------- DOUBLE PRECISION FUNCTION FMF057(DMOL) DOUBLE PRECISION CONST(1:80),DMOL,COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES c FMF057=DMOL*CONST(1)/(DMOL*CONST(1) $ +(1.0D0-DMOL)*CONST(21)) RETURN END c ----------------------------------------------- c This program transfers the unit of fraction c from kg to kmol. c ----------------------------------------------- DOUBLE PRECISION FUNCTION FMF058(DMASS) DOUBLE PRECISION CONST(1:80),DMASS,COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES c FMF058=DMASS*CONST(21)/(DMASS*CONST(21) $ +(1.0D0-DMASS)*CONST(1)) RETURN END c --------------------------------------------------- c This program solves the 3 dimensional equation. c --------------------------------------------------- DOUBLE PRECISION FUNCTION FMF059(CMPNUM,PRESS) DOUBLE PRECISION CONST(1:80),PRESS,COEFF(1:60) $ ,A0,A1,A2,A3,A,B,C,P,Q,D,U,X INTEGER TYPES(1:10),CMPNUM COMMON /FMFC/ CONST,COEFF,TYPES c A0 = CONST(20*CMPNUM - 14) A1 = CONST(20*CMPNUM - 13) A2 = CONST(20*CMPNUM - 12) A3 = CONST(20*CMPNUM - 11) A = A2/A3 B = A1/A3 C = (A0 - PRESS)/A3 P = -(A/3.0D0)**2 + B/3.0D0 Q = (A/3.0D0)**3 - A*B/6.0D0 + 0.5D0*C D = Q**2 + P**3 IF (D.GT.0.0D0) THEN IF (Q.LT.0.0D0) THEN U = -(DABS(Q - (Q**2 + P**3)**0.5))**(1.0D0/3.0D0) ELSE U = (Q + (Q**2 + P**3)**0.5)**(1.0D0/3.0D0) ENDIF X = P/U - U FMF059 = X - A/3.0D0 RETURN ELSE CALL FMF062(-2,'FMF059' $ ,'CAN NOT RESOLVE 3 DIMENSION EQUATION.') FMF059 = -1.0D20 RETURN ENDIF END c -------------------------------------------- c This program calculate c int(Cp(T))dT c (pure substance) c -------------------------------------------- DOUBLE PRECISION FUNCTION FMF060(TEMP,CMPNUM) INTEGER CMPNUM DOUBLE PRECISION TEMP,A,B,C,D,E,ECP,ECN,EEP,EEN $ ,H1,H2,H3,T DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES c T=TEMP A=COEFF(20*CMPNUM-13) B=COEFF(20*CMPNUM-12) C=COEFF(20*CMPNUM-11) D=COEFF(20*CMPNUM-10) E=COEFF(20*CMPNUM-9) ECP=DEXP(2.0D0*C/T) ECN=DEXP(-2.0D0*C/T) EEP=DEXP(2.0D0*E/T) EEN=DEXP(-2.0D0*E/T) H1=A*TEMP H2=-B*C*(1.25D-1)/T*(T*ECP-4.0D0*C-T*ECN) H3=-D*E*(1.25D-1)/T*(T*EEP+4.0D0*E-T*EEN) FMF060=H1+H2+H3 RETURN END c -------------------------------------------- c This program calculate c int(Cp(T)/T)dT c (pure substance) c -------------------------------------------- DOUBLE PRECISION FUNCTION FMF061(TEMP,CMPNUM) INTEGER CMPNUM DOUBLE PRECISION TEMP,A,B,C,D,E,ECP,ECN,EEP,EEN $ ,S1,S2,S3,T DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES c T=TEMP A=COEFF(20*CMPNUM-13) B=COEFF(20*CMPNUM-12) C=COEFF(20*CMPNUM-11) D=COEFF(20*CMPNUM-10) E=COEFF(20*CMPNUM-9) ECP=DEXP(2.0D0*C/T) ECN=DEXP(-2.0D0*C/T) EEP=DEXP(2.0D0*E/T) EEN=DEXP(-2.0D0*E/T) S1=A*DLOG(T) S2=-B*(ECP*0.0625*(2.0D0*C/T-1.0D0)-(C*0.5D0/T)**2 $ -ECN*0.0625*(2.0D0*C/T+1.0D0)) S3=-D*(EEP*0.0625*(2.0D0*E/T-1.0D0)+(E*0.5D0/T)**2 $ -EEN*0.0625*(2.0D0*E/T+1.0D0)) FMF061=S1+S2+S3 RETURN END C ---------------------------------------- C Error message on developing C ---------------------------------------- SUBROUTINE FMF062(J,PRNAME,MESSAG) INTEGER KPA,MESS,KSTAN,KAS,J COMMON /UNIT/ KPA,MESS,KSTAN,KAS CHARACTER*(*) PRNAME CHARACTER*(*) MESSAG C IF (MESS.EQ.110) THEN IF (J.EQ.0) THEN RETURN ELSE WRITE(*,1) PRNAME,MESSAG,J RETURN ENDIF ELSE RETURN ENDIF 1 FORMAT (1X,'* ERROR(',A6,') *',1X,A50,1X,'J=',I2) END c ------------------------------------------------------- c This program checks the range of temperature c fitted Cp curve for mixture. c ------------------------------------------------------- SUBROUTINE FMF063(J,TEMP) DOUBLE PRECISION TEMP,TMIN,TMAX INTEGER J,I,IMAX DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES c IMAX=TYPES(4)-2 TMIN=CONST(11) DO 10 I=1,IMAX IF (TMIN.LT.CONST(20*I+11)) THEN TMIN=CONST(20*I+11) ENDIF 10 CONTINUE TMAX=CONST(12) DO 20 I=1,IMAX IF (TMAX.GT.CONST(20*I+12)) THEN TMAX=CONST(20*I+12) ENDIF 20 CONTINUE IF (TEMP.LT.TMIN.OR.TEMP.GT.TMAX) THEN J=-2 CALL FMF062(J,'FMF063','OUT OF RANGE') RETURN ENDIF J=0 RETURN END c ------------------------------------------------------- c This program checks the range of temperature c fitted Cp curve for pure substance. c ------------------------------------------------------- SUBROUTINE FMF064(J,CMPNUM,TEMP) DOUBLE PRECISION TEMP,TMIN,TMAX INTEGER J,CMPNUM DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES c TMIN=CONST(20*CMPNUM-9) TMAX=CONST(20*CMPNUM-8) IF (TEMP.LT.TMIN.OR.TEMP.GT.TMAX) THEN J=-2 CALL FMF062(J,'FMF064','OUT OF RANGE') RETURN ENDIF J=0 RETURN END c --------------------------------------------------------- c This subroutine program sets the Coefficients c of CSD eqation of state and fundamental properties c of each substance to common variables at combi=1. c --------------------------------------------------------- SUBROUTINE START1(J,COMBI) INTEGER COMBI,I,J,TYPES(1:10),CTYPES(1:10) DOUBLE PRECISION CONST(1:80),COEFF(1:60) $ ,CCONST(1:80),CCOEFF(1:60),STAND(1:4) CHARACTER*40 NAMES(1:10),CNAMES(1:10) COMMON /FMFC/ CCONST,CCOEFF,CTYPES COMMON /NAME/ CNAMES C CALL FMF000(J,COMBI,CONST,COEFF,NAMES,TYPES) DO 10 I=1,20 CCOEFF(I)=COEFF(I) CCOEFF(20+I)=COEFF(20+I) CCOEFF(40+I)=COEFF(40+I) CCONST(I)=CONST(I) CCONST(20+I)=CONST(20+I) CCONST(40+I)=CONST(40+I) CCONST(60+I)=CONST(60+I) 10 CONTINUE DO 20 I=1,10 CNAMES(I)=NAMES(I) 20 CONTINUE DO 30 I=1,10 CTYPES(I)=TYPES(I) 30 CONTINUE CALL FMF048(J,STAND) DO 40 I=1,2 CCONST(20*I-7)=STAND(2*I-1) CCONST(20*I-6)=STAND(I*2) 40 CONTINUE RETURN END c ------------------------------------------------------------ c This subroutine program sets the Coefficients c of CSD eqation of state and fundamental properties c of each substance to common variables at combi /= 1. c ------------------------------------------------------------ SUBROUTINE START2(J,FUNPR1,FUNPR2,COEFA1,COEFB1,COEFC1 $ ,COEFA2,COEFB2,COEFC2,TIRPMT,TYPES,NAMES1,NAMES2) REAL FUNPR1(1:4),FUNPR2(1:4),COEFA1(1:3),COEFB1(1:3),COEFC1(1:5) $ ,COEFA2(1:3),COEFB2(1:3),COEFC2(1:3),TIRPMT DOUBLE PRECISION CCONST(1:80),CCOEFF(1:60) $ ,FITSAT(1:4),STAND(1:4) INTEGER CTYPES(1:10),I,J,TYPES(1:10) CHARACTER*40 CNAMES(1:10),NAMES1(1:3),NAMES2(1:3) COMMON /FMFC/ CCONST,CCOEFF,CTYPES COMMON /NAME/ CNAMES C DO 10 I=1,4 CCONST(I)=FUNPR1(I) CCONST(20+I)=FUNPR2(I) CCONST(40+I)=0.0D0 CTYPES(I)=TYPES(I) 10 CONTINUE DO 50 I=1,16 CCONST(4+I)=0.0D0 CCONST(24+I)=0.0D0 CCONST(44+I)=0.0D0 50 CONTINUE CCONST(61)=8.31451D0 CCONST(62)=0.0D0 CCONST(63)=TIRPMT CCONST(64)=0.0D0 CCONST(65)=TIRPMT CCONST(66)=0.0D0 CCONST(67)=0.0D0 CCONST(68)=0.0D0 CCONST(69)=0.0D0 CCONST(70)=0.0D0 DO 60 I=1,10 CCONST(70+I)=0.0D0 60 CONTINUE DO 40 I=1,6 CTYPES(4+I)=0 40 CONTINUE DO 20 I=1,3 CCOEFF(I)=COEFA1(I) CCOEFF(3+I)=COEFB1(I) CCOEFF(20+I)=COEFA2(I) CCOEFF(23+I)=COEFB2(I) CCOEFF(40+I)=0.0D0 CCOEFF(43+I)=0.0D0 20 CONTINUE DO 30 I=1,5 CCOEFF(6+I)=COEFC1(I) CCOEFF(26+I)=COEFC2(I) CCOEFF(46+I)=0.0D0 30 CONTINUE DO 70 I=1,9 CCOEFF(11+I)=0.0D0 CCOEFF(31+I)=0.0D0 CCOEFF(51+I)=0.0D0 70 CONTINUE CNAMES(1)=NAMES1(1) CNAMES(2)=NAMES1(2) CNAMES(3)=NAMES1(3) CNAMES(4)=NAMES2(1) CNAMES(5)=NAMES2(2) CNAMES(6)=NAMES2(3) CNAMES(7)='NULL' CNAMES(8)='NULL' CNAMES(9)='NULL' CNAMES(10)='9.1' CALL FMF047(J,FITSAT) CALL FMF048(J,STAND) DO 80 I=1,4 CCONST(5+I)=FITSAT(I) CCONST(25+I)=FITSAT(3+I) 80 CONTINUE DO 90 I=1,2 CCONST(12+I)=STAND(I) CCONST(32+I)=STAND(2+I) 90 CONTINUE RETURN END c ----------------------------------------------- c This program transfers the unit of fraction c from kmol to kg. c ----------------------------------------------- REAL FUNCTION AKG(FRMOL) REAL FRMOL DOUBLE PRECISION DMOL,FMF057,DMASS CHARACTER*6 PRNAME INTEGER J DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES c J=0 PRNAME=' AKG ' IF (FRMOL.LT.0.0.OR.FRMOL.GT.1.0) THEN J=-2 AKG=-1.0E20 CALL FMF049(J,PRNAME) RETURN ENDIF DMOL=DBLE(FRMOL) DMASS=FMF057(DMOL) AKG=SNGL(DMASS) RETURN END c ----------------------------------------------- c This program transfers the unit of fraction c from kg to kmol. c ----------------------------------------------- REAL FUNCTION AKMOL(FRMASS) REAL FRMASS DOUBLE PRECISION DMASS,FMF058,DMOL CHARACTER*6 PRNAME INTEGER J DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES c J=0 PRNAME='AKMOL ' IF (FRMASS.LT.0.0.OR.FRMASS.GT.1.0) THEN J=-2 AKMOL=-1.0E20 CALL FMF049(J,PRNAME) RETURN ENDIF DMASS=DBLE(FRMASS) DMOL=FMF058(DMASS) AKMOL=SNGL(DMOL) RETURN END c ---------------------------------------------------------- c Indicating the Region of Mixture c (modifed at 1998/06/04 by T. Yamaguchi) c ---------------------------------------------------------- INTEGER FUNCTION IPHASE(TT,PP,Z) INTEGER KPA,MESS,KSTAN,KAS,J REAL T,P,Z,TT,PP CHARACTER*6 A DOUBLE PRECISION TEMP,PRESS,VBUB,PBUB,PDEW,DBMOL,DBZ,VDEW $ ,MOLFR(1:3),PSAT,VL,VV DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /UNIT/ KPA,MESS,KSTAN,KAS COMMON /FMFC/ CONST,COEFF,TYPES C C A='IPHASE' c ------ Transforming the unit of properties ----------------- T=TT P=PP IF ((KPA.EQ.0).OR.(KPA.GE.4)) THEN P=PP*1.0E-3 ELSEIF (KPA.EQ.1) THEN T=TT+273.15 P=PP*1.0E2 ELSEIF (KPA.EQ.2) THEN P=PP*1.0E2 ELSEIF (KPA.EQ.3) THEN P=PP*1.0E-3 T=TT+273.15 ENDIF C c ============ Double precision ============================== C TEMP=DBLE(T) PRESS=DBLE(P) DBZ=DBLE(Z) CALL FMF063(J,TEMP) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF IF (DBZ.EQ.0.0D0) THEN CALL FMF029(J,2,TEMP,PSAT,VL,VV) IF (J.EQ.-2) THEN CALL FMF049(J,A) IPHASE=-2 RETURN ELSEIF (J.EQ.-1) THEN CALL FMF049(J,A) IPHASE=-1 RETURN ENDIF IF (PRESS.GT.PSAT) THEN IPHASE=1 RETURN ELSEIF (PRESS.EQ.PSAT) THEN IPHASE=2 RETURN ELSE c IPHASE=3 IPHASE=1 RETURN ENDIF ELSEIF (DBZ.EQ.1.0D0) THEN CALL FMF029(J,1,TEMP,PSAT,VL,VV) IF (J.EQ.-2) THEN CALL FMF049(J,A) IPHASE=-2 RETURN ELSEIF (J.EQ.-1) THEN CALL FMF049(J,A) IPHASE=-1 RETURN ENDIF IF (PRESS.GT.PSAT) THEN IPHASE=1 RETURN ELSEIF (PRESS.EQ.PSAT) THEN IPHASE=2 RETURN ELSE c IPHASE=3 IPHASE=1 RETURN ENDIF ENDIF IF (KAS.EQ.1) THEN DBMOL=DBZ*CONST(21)/((1.0D0-DBZ)*CONST(1)+DBZ*CONST(21)) ELSE DBMOL=DBZ ENDIF MOLFR(1)=DBMOL MOLFR(2)=1.0D0-MOLFR(1) MOLFR(3)=0.0D0 C CALL FMF037(J,TEMP,MOLFR,PBUB,VBUB) IF (J.EQ.-2) THEN CALL FMF049(J,A) IPHASE=-2 RETURN ELSEIF (J.EQ.-1) THEN CALL FMF049(J,A) IPHASE=-1 RETURN ENDIF CALL FMF038(J,TEMP,MOLFR,PDEW,VDEW) IF (J.EQ.-2) THEN CALL FMF049(J,A) IPHASE=-2 RETURN ELSEIF (J.EQ.-1) THEN CALL FMF049(J,A) IPHASE=-1 RETURN ENDIF IF (PRESS.GE.PBUB) THEN IPHASE=1 RETURN ELSEIF ((PRESS.LT.PBUB).AND.(PRESS.GT.PDEW)) THEN IPHASE=2 RETURN ELSEIF (PRESS.LE.PDEW) THEN c IPHASE=3 IPHASE=1 RETURN ENDIF END c --------------------------------------------------------- c The Identification of Ammonia and Water c --------------------------------------------------------- CHARACTER*40 FUNCTION IDENTM(I,A) INTEGER I,KPA,MESS,KSTAN,KAS CHARACTER*1 A CHARACTER*40 NAMES(1:10) COMMON /NAME/ NAMES COMMON /UNIT/ KPA,MESS,KSTAN,KAS c IF (I.EQ.1) THEN IF (A.EQ.'C') THEN IDENTM=NAMES(2) ELSEIF (A.EQ.'S') THEN IDENTM=NAMES(3) ELSEIF (A.EQ.'V') THEN IDENTM=NAMES(10) ELSE IDENTM='Out of range' ENDIF ELSEIF (I.EQ.2) THEN IF (A.EQ.'C') THEN IDENTM=NAMES(5) ELSEIF (A.EQ.'S') THEN IDENTM=NAMES(6) ELSEIF (A.EQ.'V') THEN IDENTM=NAMES(10) ELSE IDENTM='Out of range' ENDIF ELSEIF (I.EQ.3) THEN IF (A.EQ.'C') THEN IDENTM=NAMES(8) ELSEIF (A.EQ.'S') THEN IDENTM=NAMES(9) ELSEIF (A.EQ.'V') THEN IDENTM=NAMES(10) ELSE IDENTM='Out of range' ENDIF ELSE IDENTM='Out of range' ENDIF RETURN END c -------------------------------------------------------------- c Subroutine subpur c calculate the properties of pure substance c Temperature,Pressure => Volume,Enthalpy,Entropy c i=1 : 1st component c i=2 : 2nd component c -------------------------------------------------------------- SUBROUTINE SUBPUR(I,J,TT,PP,V,H,S) INTEGER I,J,KPA,MESS,KSTAN,KAS REAL T,P,V,H,S,TT,PP CHARACTER*6 A DOUBLE PRECISION TEMP,PRESS $ ,PSAT,DBV,DBH,DBS,VV,VL,VOLP(1:2),VOL(1:2) $ ,FMF030,FMF031 DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES COMMON /UNIT/ KPA,MESS,KSTAN,KAS C J=0 A='SUBPUR' C c ------ Transforming the unit of properties ----------------- T=TT P=PP IF ((KPA.EQ.0).OR.(KPA.GE.4)) THEN P=PP*1.0E-3 ELSEIF (KPA.EQ.1) THEN P=PP*1.0E2 T=TT+273.15 ELSEIF (KPA.EQ.2) THEN P=PP*1.0E2 ELSEIF (KPA.EQ.3) THEN P=PP*1.0E-3 T=TT+273.15 ENDIF C c ============ Double precision ============================== c TEMP=DBLE(T) PRESS=DBLE(P) CALL FMF064(J,I,TEMP) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF C IF (TEMP.GT.CONST(20*I-18)) THEN J=2 CALL FMF049(J,A) RETURN ENDIF c -------- judging phase ------------------- CALL FMF029(J,I,TEMP,PSAT,VL,VV) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF IF (PRESS.GT.PSAT) THEN CALL FMF027(J,I,TEMP,VOLP) CALL FMF028(J,I,TEMP,PRESS,VOLP(1),1,VOL) DBV=VOL(1) DBH=FMF030(I,TEMP,VOL(1)) DBS=FMF031(I,TEMP,VOL(1)) ELSEIF (PRESS.LT.PSAT) THEN CALL FMF027(J,I,TEMP,VOLP) CALL FMF028(J,I,TEMP,PRESS,VOLP(1),2,VOL) DBV=VOL(2) DBH=FMF030(I,TEMP,VOL(2)) DBS=FMF031(I,TEMP,VOL(2)) ELSE WRITE(*,*) ' SATURATED STATE' J=2 RETURN ENDIF C c ----------- kJ => J ----------- DBH=DBH*1.0D3 DBS=DBS*1.0D3 C c ----------- Transforming the unit of properties ------------ IF (KSTAN.EQ.1) THEN DBH=DBH+CONST(20*I-7) DBS=DBS+CONST(20*I-6) ENDIF IF (KAS.EQ.1) THEN DBV=DBV/CONST(20*I-19) DBH=DBH/CONST(20*I-19) DBS=DBS/CONST(20*I-19) ENDIF C c ============= real ================= c V=SNGL(DBV) H=SNGL(DBH) S=SNGL(DBS) RETURN END c ------------------------------------------------------------------- c Subroutine subpst c calculate the properties at saturated state of pure substance c Temperature => Pressure,Volume,Enthalpy,Entropy c i=1 : 1st component c i=2 : 2nd component c -------------------------------------------------------------------- SUBROUTINE SUBPST(I,J,TT,PS,VL,VV,HL,HV,SL,SV) INTEGER I,J,KPA,MESS,KSTAN,KAS REAL T,PS,VL,VV,HL,HV,SL,SV,TT CHARACTER*6 A DOUBLE PRECISION TEMP,PSAT,DBVL,DBVV,DBHL,DBHV,DBSL,DBSV $ ,FMF030,FMF031 DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES COMMON /UNIT/ KPA,MESS,KSTAN,KAS C J=0 A='SUBPST' C c ------- Transforming the unit of properties ---------- T=TT IF ((KPA.EQ.1).OR.(KPA.EQ.3)) THEN T=TT+273.15 ENDIF C c =========== Double precision ================ TEMP=DBLE(T) CALL FMF064(J,I,TEMP) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF c IF (TEMP.GT.CONST(20*I-18)) THEN J=2 CALL FMF049(J,A) RETURN ENDIF CALL FMF029(J,I,TEMP,PSAT,DBVL,DBVV) IF (J.NE.0) THEN J=2 CALL FMF049(J,A) RETURN ENDIF DBHL=FMF030(I,TEMP,DBVL) DBHV=FMF030(I,TEMP,DBVV) DBSL=FMF031(I,TEMP,DBVL) DBSV=FMF031(I,TEMP,DBVV) C c -------- kJ => J ------- PSAT=PSAT*1.0D3 DBHL=DBHL*1.0D3 DBHV=DBHV*1.0D3 DBSL=DBSL*1.0D3 DBSV=DBSV*1.0D3 C c -------- Transforming the unit of properties ------------- IF ((KPA.EQ.1).OR.(KPA.EQ.2)) THEN PS=SNGL(PSAT*1.0D-5) ELSE PS=SNGL(PSAT) ENDIF IF (KSTAN.EQ.1) THEN DBHL=DBHL+CONST(20*I-7) DBHV=DBHV+CONST(20*I-7) DBSL=DBSL+CONST(20*I-6) DBSV=DBSV+CONST(20*I-6) ENDIF IF (KAS.EQ.1) THEN VL=SNGL(DBVL/CONST(20*I-19)) VV=SNGL(DBVV/CONST(20*I-19)) HL=SNGL(DBHL/CONST(20*I-19)) HV=SNGL(DBHV/CONST(20*I-19)) SL=SNGL(DBSL/CONST(20*I-19)) SV=SNGL(DBSV/CONST(20*I-19)) ELSE VL=SNGL(DBVL) VV=SNGL(DBVV) HL=SNGL(DBHL) HV=SNGL(DBHV) SL=SNGL(DBSL) SV=SNGL(DBSV) ENDIF RETURN END c ------------------------------------------------------------------- c Subroutine subtsp c calculate the properties at saturated state of pure substance c Pressure => Temperature,Volume,Enthalpy,Entropy c i=1 : 1st component c i=2 : 2nd component c -------------------------------------------------------------------- SUBROUTINE SUBTSP(I,J,TS,PP,VL,VV,HL,HV,SL,SV) INTEGER I,J,KPA,MESS,KSTAN,KAS REAL TS,P,VL,VV,HL,HV,SL,SV,PP CHARACTER*6 A DOUBLE PRECISION PRESS,TSAT,DBVL,DBVV,DBHL,DBHV,DBSL,DBSV $ ,FMF030,FMF031 DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES COMMON /UNIT/ KPA,MESS,KSTAN,KAS C J=0 A='SUBTSP' C c ------- Transforming the unit of properties ---------- P=PP IF ((KPA.EQ.1).OR.(KPA.EQ.2)) THEN P=PP*1.0E2 ELSE P=PP*1.0E-3 ENDIF C c =========== Double precision ================ c PRESS=DBLE(P) c CALL FMF036(J,I,PRESS,TSAT,DBVL,DBVV) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF CALL FMF064(J,I,TSAT) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF DBHL=FMF030(I,TSAT,DBVL) DBHV=FMF030(I,TSAT,DBVV) DBSL=FMF031(I,TSAT,DBVL) DBSV=FMF031(I,TSAT,DBVV) C c -------- kJ => J ------- DBHL=DBHL*1.0D3 DBHV=DBHV*1.0D3 DBSL=DBSL*1.0D3 DBSV=DBSV*1.0D3 C c -------- Transforming the unit of properties ------------- IF ((KPA.EQ.1).OR.(KPA.EQ.3)) THEN TS=SNGL(TSAT-273.15D0) ELSE TS=SNGL(TSAT) ENDIF IF (KSTAN.EQ.1) THEN DBHL=DBHL+CONST(20*I-7) DBHV=DBHV+CONST(20*I-7) DBSL=DBSL+CONST(20*I-6) DBSV=DBSV+CONST(20*I-6) ENDIF IF (KAS.EQ.1) THEN VL=SNGL(DBVL/CONST(20*I-19)) VV=SNGL(DBVV/CONST(20*I-19)) HL=SNGL(DBHL/CONST(20*I-19)) HV=SNGL(DBHV/CONST(20*I-19)) SL=SNGL(DBSL/CONST(20*I-19)) SV=SNGL(DBSV/CONST(20*I-19)) ELSE VL=SNGL(DBVL) VV=SNGL(DBVV) HL=SNGL(DBHL) HV=SNGL(DBHV) SL=SNGL(DBSL) SV=SNGL(DBSV) ENDIF RETURN END c ------------------------------------------------------------------- c Temperature,overall composition => Pressure at bubble point c ------------------------------------------------------------------- SUBROUTINE SUBPB(J,TT,P,Z,V,H,S) INTEGER J,KPA,MESS,KSTAN,KAS REAL T,P,Z,V,H,S,TT CHARACTER*6 A DOUBLE PRECISION PRESS $ ,TEMP,DBVL,DBHL,DBSL,DBX,MOLX(1:3),VBUB,DBVV $ ,FMF020,FMF021,FMF030,FMF031 DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES COMMON /UNIT/ KPA,MESS,KSTAN,KAS C J=0 A='SUBPB ' C c ------- Transforming the unit of properties ---------- T=TT IF ((KPA.EQ.1).OR.(KPA.EQ.3)) THEN T=TT+273.15 ENDIF c c =========== Double precision ================ c TEMP=DBLE(T) CALL FMF063(J,TEMP) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF DBX=DBLE(Z) IF (DBX.EQ.0.0D0) THEN CALL FMF029(J,2,TEMP,PRESS,DBVL,DBVV) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF DBHL=FMF030(2,TEMP,DBVL) DBSL=FMF031(2,TEMP,DBVL) GOTO 1000 ELSEIF (DBX.EQ.1.0D0) THEN CALL FMF029(J,1,TEMP,PRESS,DBVL,DBVV) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF DBHL=FMF030(1,TEMP,DBVL) DBSL=FMF031(1,TEMP,DBVL) GOTO 1000 ENDIF IF (KAS.EQ.1) THEN DBX=DBX*CONST(21)/((1.0D0-DBX)*CONST(1)+DBX*CONST(21)) ENDIF MOLX(1)=DBX MOLX(2)=1.0D0-DBX MOLX(3)=0.0D0 C CALL FMF037(J,TEMP,MOLX,PRESS,VBUB) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF DBVL=VBUB DBHL=FMF020(MOLX,TEMP,VBUB) DBSL=FMF021(MOLX,TEMP,VBUB) C c --------- kJ => J ----------------------------------------- 1000 PRESS=PRESS*1.0D3 DBHL=DBHL*1.0D3 DBSL=DBSL*1.0D3 C c -------- Transforming the unit of properties ------------- IF (KPA.EQ.1.OR.KPA.EQ.2) THEN PRESS=PRESS*1.0D-5 ENDIF IF (KSTAN.EQ.1) THEN DBHL=DBHL+DBX*CONST(13)+(1.0D0-DBX)*CONST(33) DBSL=DBSL+DBX*CONST(14)+(1.0D0-DBX)*CONST(34) ENDIF IF (KAS.EQ.1) THEN DBVL=DBVL/(DBX*CONST(1)+(1.0D0-DBX)*CONST(21)) DBHL=DBHL/(DBX*CONST(1)+(1.0D0-DBX)*CONST(21)) DBSL=DBSL/(DBX*CONST(1)+(1.0D0-DBX)*CONST(21)) ENDIF c c ============= Real ===================== C P=SNGL(PRESS) V=SNGL(DBVL) H=SNGL(DBHL) S=SNGL(DBSL) RETURN END c ------------------------------------------------------------------- c Temperature,overall composition => Pressure at dew point c ------------------------------------------------------------------- SUBROUTINE SUBPD(J,TT,P,Z,V,H,S) INTEGER J,KPA,MESS,KSTAN,KAS REAL T,P,Z,V,H,S,TT CHARACTER*6 A DOUBLE PRECISION PRESS $ ,TEMP,DBVV,DBHV,DBSV,DBY,MOLY(1:3),VDEW,DBVL $ ,FMF020,FMF021,FMF030,FMF031 DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES COMMON /UNIT/ KPA,MESS,KSTAN,KAS C J=0 A='SUBPD ' C c ------- Transforming the unit of properties ---------- T=TT IF ((KPA.EQ.1).OR.(KPA.EQ.3)) THEN T=TT+273.15 ENDIF c c =========== Double precision ================ c TEMP=DBLE(T) CALL FMF063(J,TEMP) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF DBY=DBLE(Z) IF (DBY.EQ.0.0D0) THEN CALL FMF029(J,2,TEMP,PRESS,DBVL,DBVV) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF DBHV=FMF030(2,TEMP,DBVV) DBSV=FMF031(2,TEMP,DBVV) GOTO 1000 ELSEIF (DBY.EQ.1.0D0) THEN CALL FMF029(J,1,TEMP,PRESS,DBVL,DBVV) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF DBHV=FMF030(1,TEMP,DBVV) DBSV=FMF031(1,TEMP,DBVV) GOTO 1000 ENDIF IF (KAS.EQ.1) THEN DBY=DBY*CONST(21)/((1.0D0-DBY)*CONST(1)+DBY*CONST(21)) ENDIF MOLY(1)=DBY MOLY(2)=1.0D0-DBY MOLY(3)=0.0D0 C C CALL FMF038(J,TEMP,MOLY,PRESS,VDEW) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF DBVV=VDEW DBHV=FMF020(MOLY,TEMP,VDEW) DBSV=FMF021(MOLY,TEMP,VDEW) C c --------- kJ => J ----------------------------------------- 1000 PRESS=PRESS*1.0D3 DBHV=DBHV*1.0D3 DBSV=DBSV*1.0D3 C c -------- Transforming the unit of properties ------------- IF (KPA.EQ.1.OR.KPA.EQ.2) THEN PRESS=PRESS*1.0D-5 ENDIF IF (KSTAN.EQ.1) THEN DBHV=DBHV+DBY*CONST(13)+(1.0D0-DBY)*CONST(33) DBSV=DBSV+DBY*CONST(14)+(1.0D0-DBY)*CONST(34) ENDIF IF (KAS.EQ.1) THEN DBVV=DBVV/(DBY*CONST(1)+(1.0D0-DBY)*CONST(21)) DBHV=DBHV/(DBY*CONST(1)+(1.0D0-DBY)*CONST(21)) DBSV=DBSV/(DBY*CONST(1)+(1.0D0-DBY)*CONST(21)) ENDIF c c ============= Real ===================== C P=SNGL(PRESS) V=SNGL(DBVV) H=SNGL(DBHV) S=SNGL(DBSV) RETURN END c ----------------------------------------------------------------- c Subroutine subxy c calculate the composition of gas and liquid region c ----------------------------------------------------------------- SUBROUTINE SUBXY(J,TT,PP,X,Y,VL,VV,HL,HV,SL,SV) INTEGER J,KPA,MESS,KSTAN,KAS REAL T,P,X,Y,VL,VV,HL,HV,SL,SV,TT,PP CHARACTER*6 A DOUBLE PRECISION TEMP,PRESS,MOLY(1:3),TSTMN,TSTMX,MOLX(1:3) $ ,DBX,DBY,DBVL,DBVV,DBHL,DBHV,DBSL,DBSV,TST1,TST2,DVV,DVL $ ,FMF020,FMF021 DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES COMMON /UNIT/ KPA,MESS,KSTAN,KAS C J=0 A='SUBXY ' c ------ Transforming the unit of properties ----------------- T=TT P=PP IF ((KPA.EQ.0).OR.(KPA.GE.4)) THEN P=PP*1.0E-3 ELSEIF (KPA.EQ.1) THEN P=PP*1.0E2 T=TT+273.15 ELSEIF (KPA.EQ.2) THEN P=PP*1.0E2 ELSEIF (KPA.EQ.3) THEN P=PP*1.0E-3 T=TT+273.15 ENDIF C c ============ Double precision ============================== c TEMP=DBLE(T) CALL FMF063(J,TEMP) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF PRESS=DBLE(P) C c ------- judging phase ---------------- CALL FMF036(J,1,PRESS,TST1,DVL,DVV) CALL FMF036(J,2,PRESS,TST2,DVL,DVV) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF IF (TST1.GT.TST2) THEN TSTMX=TST1 TSTMN=TST2 ELSE TSTMX=TST2 TSTMN=TST1 ENDIF IF (TEMP.LE.TSTMN) THEN J=2 CALL FMF049(J,A) RETURN ELSEIF (TEMP.GT.TSTMX) THEN J=2 CALL FMF049(J,A) RETURN ENDIF C c ---------- calculating properties ------------- CALL FMF046(J,TEMP,PRESS,MOLX,MOLY,DBVL,DBVV) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF DBX=MOLX(1) DBY=MOLY(1) DBHL=FMF020(MOLX,TEMP,DBVL) DBHV=FMF020(MOLY,TEMP,DBVV) DBSL=FMF021(MOLX,TEMP,DBVL) DBSV=FMF021(MOLY,TEMP,DBVV) c c --------- kJ => J ----------------------------------------- DBHL=DBHL*1.0D3 DBHV=DBHV*1.0D3 DBSL=DBSL*1.0D3 DBSV=DBSV*1.0D3 C c -------- Transforming the unit of properties ------------- IF (KSTAN.EQ.1) THEN DBHL=DBHL+DBX*CONST(13)+(1.0D0-DBX)*CONST(33) DBHV=DBHV+DBY*CONST(13)+(1.0D0-DBY)*CONST(33) DBSL=DBSL+DBX*CONST(14)+(1.0D0-DBX)*CONST(34) DBSV=DBSV+DBY*CONST(14)+(1.0D0-DBY)*CONST(34) ENDIF IF (KAS.EQ.1) THEN DBVL=DBVL/(DBX*CONST(1)+(1.0D0-DBX)*CONST(21)) DBVV=DBVV/(DBY*CONST(1)+(1.0D0-DBY)*CONST(21)) DBHL=DBHL/(DBX*CONST(1)+(1.0D0-DBX)*CONST(21)) DBHV=DBHV/(DBY*CONST(1)+(1.0D0-DBY)*CONST(21)) DBSL=DBSL/(DBX*CONST(1)+(1.0D0-DBX)*CONST(21)) DBSV=DBSV/(DBY*CONST(1)+(1.0D0-DBY)*CONST(21)) DBX=DBX*CONST(1)/((1.0D0-DBX)*CONST(21)+DBX*CONST(1)) DBY=DBY*CONST(1)/((1.0D0-DBY)*CONST(21)+DBY*CONST(1)) ENDIF C c ============= Real ===================== c X=SNGL(DBX) Y=SNGL(DBY) VL=SNGL(DBVL) VV=SNGL(DBVV) HL=SNGL(DBHL) HV=SNGL(DBHV) SL=SNGL(DBSL) SV=SNGL(DBSV) RETURN END c ------------------------------------------------------------- c Subroutine submix c calculate the properties of ammonia-water mixture c i=1 : T,p,z => v,h,s c i=2 : p,z,h => T,v,s c i=3 : p,z,s => T,v,h c i=4 : p,z,v => T,h,s c ------------------------------------------------------------- SUBROUTINE SUBMIX(I,J,TT,PP,Z,V,H,S) INTEGER I,J,KPA,MESS,KSTAN,KAS,PHASE REAL T,P,Z,V,H,S,TT,PP CHARACTER*6 A DOUBLE PRECISION DBH,DBS,DBV,TEMP,PRESS,DBMOL,DBZ $ ,DBHL,DBHV,DBSL,DBSV,DBVL,DBVV,QUAL,MOLX(1:3),MOLY(1:3) $ ,DBX,DBY,MOLFR(1:3),VOLP(1:2),VOL(1:2) $ ,FMF020,FMF021 INTEGER FMF050 DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES COMMON /UNIT/ KPA,MESS,KSTAN,KAS C J=0 A='SUBMIX' C c ************************************************************ c ------ i=1 ; (z,t,p) => (v,h,s) ---------------------------- IF (I.EQ.1) THEN c ************************************************************ c c ------ Transforming the unit of properties ----------------- T=TT P=PP IF ((KPA.EQ.0).OR.(KPA.GE.4)) THEN P=PP*1.0E-3 ELSEIF (KPA.EQ.1) THEN P=PP*1.0E2 T=TT+273.15 ELSEIF (KPA.EQ.2) THEN P=PP*1.0E2 ELSEIF (KPA.EQ.3) THEN P=PP*1.0E-3 T=TT+273.15 ENDIF C c ============ Double precision ============================== c DBZ=DBLE(Z) C TEMP=DBLE(T) CALL FMF063(J,TEMP) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF PRESS=DBLE(P) IF (DBZ.EQ.0.0D0) THEN CALL FMF032(J,2,TEMP,PRESS,DBV,DBH,DBS) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF GOTO 1000 ELSEIF (DBZ.EQ.1.0D0) THEN CALL FMF032(J,1,TEMP,PRESS,DBV,DBH,DBS) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF GOTO 1000 ENDIF IF (KAS.EQ.1) THEN DBMOL=DBZ*CONST(21)/((1.0D0-DBZ)*CONST(1)+DBZ*CONST(21)) ELSE DBMOL=DBZ ENDIF MOLFR(1)=DBMOL MOLFR(2)=1.0D0-MOLFR(1) MOLFR(3)=0.0D0 C c ------- judging the phase of mixture ------------- PHASE=FMF050(TEMP,PRESS,MOLFR) c write(*,*) ' phase=',phase IF (PHASE.EQ.-1) THEN J=1 CALL FMF049(J,A) RETURN ELSEIF (PHASE.EQ.-2) THEN J=2 CALL FMF049(J,A) RETURN ELSEIF (PHASE.EQ.1) THEN CALL FMF022(J,MOLFR,TEMP,VOLP) CALL FMF024(J,MOLFR,TEMP,PRESS,VOLP(1),1,VOL) DBV=VOL(1) DBH=FMF020(MOLFR,TEMP,VOL(1)) DBS=FMF021(MOLFR,TEMP,VOL(1)) ELSEIF (PHASE.EQ.2) THEN CALL FMF046(J,TEMP,PRESS,MOLX,MOLY,DBVL,DBVV) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF DBX=MOLX(1) DBY=MOLY(1) QUAL=(DBMOL-DBX)/(DBY-DBX) DBHL=FMF020(MOLX,TEMP,DBVL) DBHV=FMF020(MOLY,TEMP,DBVV) DBSL=FMF021(MOLX,TEMP,DBVL) DBSV=FMF021(MOLY,TEMP,DBVV) DBV=(1.0D0-QUAL)*DBVL+QUAL*DBVV DBH=(1.0D0-QUAL)*DBHL+QUAL*DBHV DBS=(1.0D0-QUAL)*DBSL+QUAL*DBSV ELSEIF (PHASE.EQ.3) THEN CALL FMF022(J,MOLFR,TEMP,VOLP) CALL FMF024(J,MOLFR,TEMP,PRESS,VOLP(2),2,VOL) DBV=VOL(2) DBH=FMF020(MOLFR,TEMP,VOL(2)) DBS=FMF021(MOLFR,TEMP,VOL(2)) ENDIF C c --------- kJ => J ---------------- 1000 DBH=DBH*1.0D3 DBS=DBS*1.0D3 C c ---------- Transfoming the unit of properties ------------ IF (KSTAN.EQ.1) THEN DBH=DBH+DBMOL*CONST(13)+(1.0D0-DBMOL)*CONST(33) DBS=DBS+DBMOL*CONST(14)+(1.0D0-DBMOL)*CONST(34) ENDIF IF (KAS.EQ.1) THEN DBV=DBV/(DBMOL*CONST(1)+(1.0D0-DBMOL)*CONST(21)) DBH=DBH/(DBMOL*CONST(1)+(1.0D0-DBMOL)*CONST(21)) DBS=DBS/(DBMOL*CONST(1)+(1.0D0-DBMOL)*CONST(21)) ENDIF C c --------- real ------------------- V=SNGL(DBV) H=SNGL(DBH) S=SNGL(DBS) RETURN C c ************************************************************ c ------ i=2 ; (z,p,h) => (t,v,s) ---------------------------- ELSEIF (I.EQ.2) THEN c ************************************************************ c c ------ Transforming the unit of properties ----------------- P=PP IF ((KPA.EQ.0).OR.(KPA.GE.3)) THEN P=PP*1.0E-3 ELSE P=PP*1.0E2 ENDIF C c ============ Double precision ============================== c DBZ=DBLE(Z) DBH=DBLE(H)*1.0D-3 C c ------ Transforming the unit of properties ----------------- IF (KAS.EQ.1) THEN DBMOL=DBZ*CONST(21)/((1.0D0-DBZ)*CONST(1)+DBZ*CONST(21)) DBH=DBH*(CONST(1)*DBMOL+(1.0D0-DBMOL)*CONST(21)) ELSE DBMOL=DBZ ENDIF IF (KSTAN.EQ.1) THEN DBH=DBH-(DBMOL*CONST(13)+(1.0D0-DBMOL)*CONST(33))*1.0D-3 ENDIF MOLFR(1)=DBMOL MOLFR(2)=1.0D0-MOLFR(1) MOLFR(3)=0.0D0 C PRESS=DBLE(P) IF (DBMOL.EQ.0.0D0) THEN CALL FMF051(J,2,TEMP,PRESS,DBV,DBH,DBS) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF GOTO 1100 ELSEIF (DBMOL.EQ.1.0D0) THEN CALL FMF051(J,1,TEMP,PRESS,DBV,DBH,DBS) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF GOTO 1100 ENDIF CALL FMF040(J,MOLFR,TEMP,PRESS,DBV,DBH,DBS,MOLX,MOLY) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF CALL FMF063(J,TEMP) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF C 1100 DBS=DBS*1.0D3 C c ---------- Transfoming the unit of properties ------------ IF (KSTAN.EQ.1) THEN DBS=DBS+DBMOL*CONST(14)+(1.0D0-DBMOL)*CONST(34) ENDIF IF (KAS.EQ.1) THEN DBV=DBV/(DBMOL*CONST(1)+(1.0D0-DBMOL)*CONST(21)) DBS=DBS/(DBMOL*CONST(1)+(1.0D0-DBMOL)*CONST(21)) ENDIF IF ((KPA.EQ.1).OR.(KPA.EQ.3)) THEN TEMP=TEMP-273.15 ENDIF C c --------- real ------------------- V=SNGL(DBV) TT=SNGL(TEMP) S=SNGL(DBS) RETURN C c ************************************************************ c ------ i=3 ; (z,p,s) => (t,v,h) ---------------------------- ELSEIF (I.EQ.3) THEN c ************************************************************ c c ------ Transforming the unit of properties ----------------- P=PP IF ((KPA.EQ.0).OR.(KPA.GE.3)) THEN P=PP*1.0E-3 ELSE P=PP*1.0E2 ENDIF C c ============ Double precision ============================== C DBZ=DBLE(Z) DBS=DBLE(S)*1.0D-3 C c ------ Transforming the unit of properties ----------------- IF (KAS.EQ.1) THEN DBMOL=DBZ*CONST(21)/((1.0D0-DBZ)*CONST(1)+DBZ*CONST(21)) DBS=DBS*(CONST(1)*DBMOL+(1.0D0-DBMOL)*CONST(21)) ELSE DBMOL=DBZ ENDIF IF (KSTAN.EQ.1) THEN DBS=DBS-(DBMOL*CONST(14)+(1.0D0-DBMOL)*CONST(34))*1.0D-3 ENDIF MOLFR(1)=DBMOL MOLFR(2)=1.0D0-MOLFR(1) MOLFR(3)=0.0D0 c PRESS=DBLE(P) IF (DBMOL.EQ.0.0D0) THEN CALL FMF052(J,2,TEMP,PRESS,DBV,DBH,DBS) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF GOTO 1200 ELSEIF (DBMOL.EQ.1.0D0) THEN CALL FMF052(J,1,TEMP,PRESS,DBV,DBH,DBS) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF GOTO 1200 ENDIF CALL FMF041(J,MOLFR,TEMP,PRESS,DBV,DBH,DBS,MOLX,MOLY) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF CALL FMF063(J,TEMP) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF C c --------- kJ => J ---------------- 1200 DBH=DBH*1.0D3 C c ---------- Transfoming the unit of properties ------------ IF (KSTAN.EQ.1) THEN DBH=DBH+DBMOL*CONST(13)+(1.0D0-DBMOL)*CONST(33) ENDIF IF (KAS.EQ.1) THEN DBV=DBV/(DBMOL*CONST(1)+(1.0D0-DBMOL)*CONST(21)) DBH=DBH/(DBMOL*CONST(1)+(1.0D0-DBMOL)*CONST(21)) ENDIF IF ((KPA.EQ.1).OR.(KPA.EQ.3)) THEN TEMP=TEMP-273.15 ENDIF C c --------- real ------------------- V=SNGL(DBV) TT=SNGL(TEMP) H=SNGL(DBH) RETURN C c ************************************************************ c ------ i=4 ; (z,p,v) => (t,h,s) ---------------------------- ELSEIF (I.EQ.4) THEN c ************************************************************ c c ------ Transforming the unit of properties ----------------- P=PP IF ((KPA.EQ.0).OR.(KPA.GE.3)) THEN P=PP*1.0E-3 ELSE P=PP*1.0E2 ENDIF C c ============ Double precision ============================== c DBZ=DBLE(Z) C c ------ Transforming the unit of properties ----------------- IF (KAS.EQ.1) THEN DBMOL=DBZ*CONST(21)/((1.0D0-DBZ)*CONST(1)+DBZ*CONST(21)) DBV=DBLE(V)*(CONST(1)*DBMOL+(1.0D0-DBMOL)*CONST(21)) ELSE DBMOL=DBZ DBV=DBLE(V) ENDIF C PRESS=DBLE(P) MOLFR(1)=DBMOL MOLFR(2)=1.0D0-MOLFR(1) MOLFR(3)=0.0D0 IF (DBMOL.EQ.0.0D0) THEN CALL FMF053(J,2,TEMP,PRESS,DBV,DBH,DBS) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF GOTO 1300 ELSEIF (DBMOL.EQ.1.0D0) THEN CALL FMF053(J,1,TEMP,PRESS,DBV,DBH,DBS) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF GOTO 1300 ENDIF c CALL FMF042(J,MOLFR,TEMP,PRESS,DBV,DBH,DBS,MOLX,MOLY) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF CALL FMF063(J,TEMP) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF C c --------- kJ => J ---------------- 1300 DBH=DBH*1.0D3 DBS=DBS*1.0D3 C c ---------- Transfoming the unit of properties ------------ IF (KSTAN.EQ.1) THEN DBH=DBH+DBMOL*CONST(13)+(1.0D0-DBMOL)*CONST(33) DBS=DBS+DBMOL*CONST(14)+(1.0D0-DBMOL)*CONST(34) ENDIF IF (KAS.EQ.1) THEN DBH=DBH/(DBMOL*CONST(1)+(1.0D0-DBMOL)*CONST(21)) DBS=DBS/(DBMOL*CONST(1)+(1.0D0-DBMOL)*CONST(21)) ENDIF IF ((KPA.EQ.1).OR.(KPA.EQ.3)) THEN TEMP=TEMP-273.15 ENDIF C c --------- real ------------------- TT=SNGL(TEMP) H=SNGL(DBH) S=SNGL(DBS) RETURN ELSE J=2 CALL FMF049(J,A) RETURN ENDIF END c -------------------------------------------------- c Fundamental properties c -------------------------------------------------- REAL FUNCTION FCM(I,A) INTEGER I,KPA,MESS,KSTAN,KAS,COMP,TYPES(1:10) CHARACTER*1 A CHARACTER*6 B CHARACTER*40 CNAMES(1:10) DOUBLE PRECISION CONST(1:80),COEFF(1:60) COMMON /UNIT/ KPA,MESS,KSTAN,KAS COMMON /FMFC/ CONST,COEFF,TYPES COMMON /NAME/ CNAMES C B=' FCM ' C IF (I.EQ.1) THEN COMP=0 ELSEIF (I.EQ.2) THEN COMP=20 ELSEIF (I.EQ.3) THEN COMP=40 ELSE FCM=-1.0E20 RETURN ENDIF C IF (A.EQ.'M') THEN FCM=SNGL(CONST(COMP+1)) RETURN ELSEIF (A.EQ.'R') THEN FCM=SNGL(CONST(61)/CONST(COMP+1))*1.0E3 RETURN ELSEIF (A.EQ.'T') THEN IF ((KPA.EQ.1).OR.(KPA.EQ.3)) THEN FCM=SNGL(CONST(COMP+2))-273.15 ELSE FCM=SNGL(CONST(COMP+2)) ENDIF RETURN ELSEIF (A.EQ.'P') THEN IF ((KPA.EQ.1).OR.(KPA.EQ.2)) THEN FCM=SNGL(CONST(COMP+3))*1.0E-5 ELSE FCM=SNGL(CONST(COMP+3)) ENDIF RETURN ELSEIF (A.EQ.'V') THEN IF (KAS.EQ.1) THEN FCM=SNGL(CONST(COMP+4)/CONST(COMP+1)) ELSE FCM=SNGL(CONST(COMP+4)) ENDIF RETURN ELSE CALL FMF049(2,B) FCM=-1.0E20 RETURN ENDIF END c ------------------------------------------------------------- c Temperature => Saturated pressure of pure substance c i=1 : 1st component c i=2 : 2nd component c ------------------------------------------------------------- REAL FUNCTION PSTM(I,TT) INTEGER I,J,KPA,MESS,KSTAN,KAS REAL T,TT CHARACTER*6 A DOUBLE PRECISION PSAT,TEMP,VL,VV DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES COMMON /UNIT/ KPA,MESS,KSTAN,KAS C A=' PSTM ' IF (I.LT.1.OR.I.GT.2) THEN CALL FMF049(2,A) PSTM=-1.0E20 RETURN ENDIF c c --------- Transforming The Unit of properties -------- T=TT IF ((KPA.EQ.1).OR.(KPA.EQ.3)) THEN T=TT+273.15 ENDIF C c ---------- DOUBLE PRECISION ------- TEMP=DBLE(T) CALL FMF063(J,TEMP) IF (J.NE.0) THEN CALL FMF049(J,A) PSTM=-1.0E20 RETURN ENDIF C IF (TEMP.GT.CONST(20*I-18)) THEN CALL FMF049(2,A) PSTM=-1.0E20 RETURN ENDIF CALL FMF029(J,I,TEMP,PSAT,VL,VV) IF (J.EQ.-2) THEN CALL FMF049(J,A) PSTM=-1.0E20 RETURN ELSEIF (J.EQ.-1) THEN CALL FMF049(J,A) PSTM=-1.0E10 RETURN ENDIF C c ---------- Transforming the unit of property ------------ IF ((KPA.EQ.1).OR.(KPA.EQ.2)) THEN PSTM=SNGL(PSAT*1.0D-2) ELSE PSTM=SNGL(PSAT*1.0D3) ENDIF RETURN END c ------------------------------------------------------------- c Pressure => Saturated Temperature of pure substance c i=1 : 1st component c i=2 : 2nd component c ------------------------------------------------------------- REAL FUNCTION TSPM(I,PP) INTEGER I,J,KPA,MESS,KSTAN,KAS REAL P,PP CHARACTER*6 A DOUBLE PRECISION TSAT,PRESS,VL,VV DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES COMMON /UNIT/ KPA,MESS,KSTAN,KAS C A=' TSPM ' c --------- Transforming The Unit of properties -------- P=PP IF ((KPA.EQ.0).OR.(KPA.GE.3)) THEN P=PP*1.0E-3 ELSE P=PP*1.0E2 ENDIF C PRESS=DBLE(P) c CALL FMF036(J,I,PRESS,TSAT,VL,VV) IF (J.EQ.-2) THEN CALL FMF049(J,A) TSPM=-1.0E20 RETURN ELSEIF (J.EQ.-1) THEN CALL FMF049(J,A) TSPM=-1.0E10 RETURN ENDIF CALL FMF063(J,TSAT) IF (J.NE.0) THEN CALL FMF049(J,A) TSPM=-1.0E20 RETURN ENDIF C c ---------- Transforming the unit of property ------------ IF ((KPA.EQ.1).OR.(KPA.EQ.3)) THEN TSPM=SNGL(TSAT-273.15D0) ELSE TSPM=SNGL(TSAT) ENDIF RETURN END c --------------------------------------------- c Seting the common variables KPA and MESS c --------------------------------------------- SUBROUTINE KPAMES(KPAC,MESSC) INTEGER KPAC,MESSC,KPA,MESS,KSTAN,KAS COMMON /UNIT/ KPA,MESS,KSTAN,KAS C KPA=KPAC MESS=MESSC RETURN END c --------------------------------------------------- c Seting the common variables KSTAN and KAS c --------------------------------------------------- SUBROUTINE STNKAS(STAND,KASC) INTEGER KPA,MESS,KSTAN,KAS,STAND,KASC COMMON /UNIT/ KPA,MESS,KSTAN,KAS C KSTAN=STAND KAS=KASC RETURN END c ------------------------------------------------------------------- c Temperature,overall composition => Pressure at bubble point c ------------------------------------------------------------------- SUBROUTINE SUBTB(J,T,PP,Z,V,H,S) INTEGER J,KPA,MESS,KSTAN,KAS REAL T,P,Z,V,H,S,PP CHARACTER*6 A DOUBLE PRECISION PRESS,TEMP,MOLX(1:3) $ ,DBVL,DBHL,DBSL,DBX,DBVV $ ,FMF020,FMF021,FMF030,FMF031 DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES COMMON /UNIT/ KPA,MESS,KSTAN,KAS C J=0 A='SUBTB ' C c ------- Transforming the unit of properties ---------- P=PP IF ((KPA.EQ.1).OR.(KPA.EQ.2)) THEN P=PP*1.0E2 ELSE P=PP*1.0E-3 ENDIF C c =========== Double precision ================ C PRESS=DBLE(P) DBX=DBLE(Z) IF (DBX.EQ.0.0D0) THEN CALL FMF036(J,2,PRESS,TEMP,DBVL,DBVV) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF CALL FMF063(J,TEMP) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF DBHL=FMF030(2,TEMP,DBVL) DBSL=FMF031(2,TEMP,DBVL) GOTO 1000 ELSEIF (DBX.EQ.1.0D0) THEN CALL FMF036(J,1,PRESS,TEMP,DBVL,DBVV) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF CALL FMF063(J,TEMP) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF DBHL=FMF030(1,TEMP,DBVL) DBSL=FMF031(1,TEMP,DBVL) GOTO 1000 ENDIF IF (KAS.EQ.1) THEN DBX=DBX*CONST(21)/((1.0D0-DBX)*CONST(1)+DBX*CONST(21)) ENDIF MOLX(1)=DBX MOLX(2)=1.0D0-MOLX(1) MOLX(3)=0.0D0 C C CALL FMF043(J,PRESS,MOLX,TEMP,DBVL) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF CALL FMF063(J,TEMP) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF DBHL=FMF020(MOLX,TEMP,DBVL) DBSL=FMF021(MOLX,TEMP,DBVL) c c --------- kJ => J ----------------------------------------- 1000 DBHL=DBHL*1.0D3 DBSL=DBSL*1.0D3 C c -------- Transforming the unit of properties ------------- IF ((KPA.EQ.1).OR.(KPA.EQ.3)) THEN TEMP=TEMP-273.15D0 ENDIF IF (KSTAN.EQ.1) THEN DBHL=DBHL+DBX*CONST(13)+(1.0D0-DBX)*CONST(33) DBSL=DBSL+DBX*CONST(14)+(1.0D0-DBX)*CONST(34) ENDIF IF (KAS.EQ.1) THEN DBVL=DBVL/(DBX*CONST(1)+(1.0D0-DBX)*CONST(21)) DBHL=DBHL/(DBX*CONST(1)+(1.0D0-DBX)*CONST(21)) DBSL=DBSL/(DBX*CONST(1)+(1.0D0-DBX)*CONST(21)) ENDIF c c ============= Real ===================== c T=SNGL(TEMP) V=SNGL(DBVL) H=SNGL(DBHL) S=SNGL(DBSL) RETURN END c ------------------------------------------------------------------- c Temperature,overall composition => Pressure at dew point c ------------------------------------------------------------------- SUBROUTINE SUBTD(J,T,PP,Z,V,H,S) INTEGER J,KPA,MESS,KSTAN,KAS REAL T,P,Z,V,H,S,PP CHARACTER*6 A DOUBLE PRECISION PRESS,TEMP,MOLY(1:3) $ ,DBVV,DBHV,DBSV,DBY,DBVL $ ,FMF020,FMF021,FMF030,FMF031 DOUBLE PRECISION CONST(1:80),COEFF(1:60) INTEGER TYPES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES COMMON /UNIT/ KPA,MESS,KSTAN,KAS C J=0 A='SUBTD ' C c ------- Transforming the unit of properties ---------- P=PP IF ((KPA.EQ.1).OR.(KPA.EQ.2)) THEN P=PP*1.0E2 ELSE P=PP*1.0E-3 ENDIF C c =========== Double precision ================ C PRESS=DBLE(P) DBY=DBLE(Z) IF (DBY.EQ.0.0D0) THEN CALL FMF036(J,2,PRESS,TEMP,DBVL,DBVV) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF CALL FMF063(J,TEMP) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF DBHV=FMF030(2,TEMP,DBVV) DBSV=FMF031(2,TEMP,DBVV) GOTO 1000 ELSEIF (DBY.EQ.1.0D0) THEN CALL FMF036(J,1,PRESS,TEMP,DBVL,DBVV) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF CALL FMF063(J,TEMP) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF DBHV=FMF030(1,TEMP,DBVV) DBSV=FMF031(1,TEMP,DBVV) GOTO 1000 ENDIF IF (KAS.EQ.1) THEN DBY=DBY*CONST(21)/((1.0D0-DBY)*CONST(1)+DBY*CONST(21)) ENDIF MOLY(1)=DBY MOLY(2)=1.0D0-MOLY(1) MOLY(3)=0.0D0 C C CALL FMF044(J,PRESS,MOLY,TEMP,DBVV) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF CALL FMF063(J,TEMP) IF (J.NE.0) THEN CALL FMF049(J,A) RETURN ENDIF DBHV=FMF020(MOLY,TEMP,DBVV) DBSV=FMF021(MOLY,TEMP,DBVV) c c --------- kJ => J ----------------------------------------- 1000 DBHV=DBHV*1.0D3 DBSV=DBSV*1.0D3 C c -------- Transforming the unit of properties ------------- IF ((KPA.EQ.1).OR.(KPA.EQ.3)) THEN TEMP=TEMP-273.15D0 ENDIF IF (KSTAN.EQ.1) THEN DBHV=DBHV+DBY*CONST(13)+(1.0D0-DBY)*CONST(33) DBSV=DBSV+DBY*CONST(14)+(1.0D0-DBY)*CONST(34) ENDIF IF (KAS.EQ.1) THEN DBVV=DBVV/(DBY*CONST(1)+(1.0D0-DBY)*CONST(21)) DBHV=DBHV/(DBY*CONST(1)+(1.0D0-DBY)*CONST(21)) DBSV=DBSV/(DBY*CONST(1)+(1.0D0-DBY)*CONST(21)) ENDIF c c ============= Real ===================== c T=SNGL(TEMP) V=SNGL(DBVV) H=SNGL(DBHV) S=SNGL(DBSV) RETURN END c -------------------------------------------------- c This program makes the table of pure substance's c fundamental constants. c -------------------------------------------------- SUBROUTINE MKTABL(J,COMBI) INTEGER J,COMBI,TYPES(1:10) DOUBLE PRECISION CONST(1:80),COEFF(1:60),R1,R2 CHARACTER*6 PRNAME CHARACTER*40 NAMES(1:10) COMMON /FMFC/ CONST,COEFF,TYPES COMMON /NAME/ NAMES c PRNAME='MKTABL' CALL START1(J,COMBI) IF (J.NE.0) THEN CALL FMF049(J,PRNAME) RETURN ENDIF R1=CONST(61)/CONST(1)*1.0E3 R2=CONST(61)/CONST(21)*1.0E3 WRITE(*,*) ' You do not need to CALL SUBROUTINE' WRITE(*,*) ' KPAMES, STNKAS and START1 before' WRITE(*,*) ' CALLing MKTABL' WRITE(*,*) WRITE(*,*) ' result' WRITE(*,*) WRITE(*,*) ' FUNDAMENTAL CONSTANTS' WRITE(*,1) NAMES(1),NAMES(4) WRITE(*,2) 'MOLECULAR FORMULA ' $ ,NAMES(2),NAMES(5) WRITE(*,3) 'RELATIVE MOLECULAR MASS [KG/KMOL] ' $ ,CONST(1),CONST(21) WRITE(*,3) 'GAS CONSTANT [J/KG/K] ' $ ,R1,R2 WRITE(*,4) 'CRITICAL TEMPERATURE [K] ' $ ,CONST(2),CONST(22) WRITE(*,5) 'CRITICAL PRESSURE [MPA] ' $ ,CONST(3)*1.0D-6,CONST(23)*1.0D-6 WRITE(*,6) 'CRITICAL VOLUME [M^3/KMOL]' $ ,CONST(4),CONST(24) WRITE(*,5) 'ACENTRIC FACTOR [-] ' $ ,CONST(5),CONST(25) WRITE(*,7) 'INTERACTION PARAMETER [-] ' $ ,CONST(63) WRITE(*,*) WRITE(*,*) ' IDEAL GAS HEAT CAPACITY' WRITE(*,1) NAMES(1),NAMES(4) WRITE(*,8) 'EQUATION NUMBER',TYPES(1),TYPES(2) WRITE(*,9) 'A',COEFF(7),COEFF(27) WRITE(*,9) 'B',COEFF(8),COEFF(28) WRITE(*,9) 'C',COEFF(9),COEFF(29) WRITE(*,9) 'D',COEFF(10),COEFF(30) WRITE(*,9) 'E',COEFF(11),COEFF(31) WRITE(*,10) 'LOW TEMPERATURE LIMIT [K]',CONST(11),CONST(31) WRITE(*,10) 'HIGH TEMPERATURE LIMIT [K]',CONST(12),CONST(32) WRITE(*,*) WRITE(*,*) ' PRESS RETURN TO CONTINUE' READ(*,*) WRITE(*,*) ' THE COEFFICIENTS OF CSD EQUATION OF STATE' WRITE(*,1) NAMES(1),NAMES(4) WRITE(*,9) 'A0',COEFF(1),COEFF(21) WRITE(*,9) 'A2',COEFF(2),COEFF(22) WRITE(*,9) 'A3',COEFF(3),COEFF(23) WRITE(*,9) 'B0',COEFF(4),COEFF(24) WRITE(*,9) 'B1',COEFF(5),COEFF(25) WRITE(*,9) 'B2',COEFF(6),COEFF(26) WRITE(*,*) WRITE(*,*) 'THE COEFFICIENTS OF EQUATION FITTED SATURATION CURVE' WRITE(*,1) NAMES(1),NAMES(4) WRITE(*,9) 'S0',CONST(6),CONST(26) WRITE(*,9) 'S1',CONST(7),CONST(27) WRITE(*,9) 'S2',CONST(8),CONST(28) WRITE(*,9) 'S3',CONST(9),CONST(29) 1 FORMAT(37X,2(6X,A9)) 2 FORMAT(1X,A36,2(6X,A9)) 3 FORMAT(1X,A36,2F15.3) 4 FORMAT(1X,A36,2F15.2) 5 FORMAT(1X,A36,2F15.4) 6 FORMAT(1X,A36,2F15.6) 7 FORMAT(1X,A36,F20.4) 8 FORMAT(1X,A36,2I15) 9 FORMAT(1X,A36,2E15.6) 10 FORMAT(1X,A36,2E15.1) RETURN END