program co2tab3b c common/ff/h1 character h1*1,vv*5 dimension dens(2),tcond(2),visc(2),beta(2),h(2) dimension p(200),t(100),rho(100),vis(100),enth(100),vv(26) dimension psat(100),rhog(100),visg(100),enthg(100) logical ex SAVE ICALL DATA ICALL/0/ C ICALL=ICALL+1 C print 100 100 format(1X,131('*')//51X,'Tabulation of PVT Data for CO2'// x1X,131('*')//) c INQUIRE(FILE='VERTAB',EXIST=EX) IF(EX) GOTO 62 PRINT 63 63 FORMAT(' FILE *VERTAB* DOES NOT EXIST --- OPEN AS A NEW FILE') OPEN(11,FILE='VERTAB',STATUS='NEW') ENDFILE 11 GOTO 70 C 62 PRINT 64 64 FORMAT(' FILE *VERTAB* EXISTS --- OPEN AS AN OLD FILE') OPEN(11,FILE='VERTAB',STATUS='OLD') C 70 CONTINUE REWIND 11 C IF(ICALL.EQ.1) WRITE(11,899) c 899 FORMAT(/6X,'CBCO2 1.0 18 November 1999',6X, c 899 FORMAT(/6X,'co2tab 0.9 14 March 2000',6X, c 899 FORMAT(/6X,'CO2TAB 0.9 12 May 2000',6X, c 899 FORMAT(/6X,'CO2TAB1 1.0 26 February 2002',6X, c 899 FORMAT(/6X,'CO2TAB2 1.0 28 February 2002',6X, c 899 FORMAT(/6X,'CO2TAB3 1.0 1 March 2002',6X, 899 FORMAT(/6X,'CO2TAB3a 1.0 24 June 2005',6X, X'tabulate thermophysical properties of CO2 using correlations',/ x50X,'due to V.V. Altunin (1975)'/ x47X,'as coded by V.I. Malkovsky (IGEM, Moscow, Russia)') c c---*----1----*----2----*----3----*----4----*----5----*----6----*----7----*----8 c read pressures below lowest Psat read 1,p(1),dp,np1 1 format(2E10.4,5X,I5) do13 i=2,np1 13 p(i)=p(i-1)+dp c c read temperatures i0=0 read 1,t(1),dt,nt if(t(1).le.31.04) i0=1 do3 j=2,nt t(j)=t(j-1)+dt if(t(j).le.31.04) i0=j 3 continue c i0 holds the number of sub-critical points c c read supercritical pressures read 1,p(np1+i0+1),dp,np2 c do2 i=2,np2 2 p(np1+i0+i)=p(np1+i0+i-1)+dp np=np1+i0+np2 c c.....6-24-05: append a point with 600 bar if(p(np1+i0+np2).lt.600.e5) then p(np1+i0+np2+1)=600.e5 np=np+1 endif C INQUIRE(FILE='CO2TAB',EXIST=EX) IF(EX) GOTO 162 PRINT 163 163 FORMAT(' FILE *CO2TAB* DOES NOT EXIST --- OPEN AS A NEW FILE') OPEN(2,FILE='CO2TAB',STATUS='NEW') ENDFILE 2 GOTO 170 C 162 PRINT 164 164 FORMAT(' FILE *CO2TAB* EXISTS --- OPEN AS AN OLD FILE') OPEN(2,FILE='CO2TAB',STATUS='OLD') C 170 CONTINUE REWIND 2 C if(i0.ge.1) then c.....now calculate sub-critical data do24 i=1,i0 mode=2 tx=t(i)+273.15 CALL CBCO2(Px,Tx,DENS,TCOND,VISC,BETA,H,RVAP,SIGMA,IND,MODE) psat(i)=px*1.e5 p(np1+i)=psat(i) rhog(i)=dens(2) visg(i)=visc(2) enthg(i)=h(2) 24 continue c print 25,(t(i),i=1,i0) 25 format(1X,'CO2 saturation line temperatures'/ x (10(1X,E12.6))) print 6,(psat(i),i=1,i0) 6 format(1X,'CO2 saturation line pressures'/ x (10(1X,E12.6))) print 7,(rhog(i),i=1,i0) 7 format(1X,'CO2 saturation line gas densities'/ x (10(1X,E12.6))) print 8,(visg(i),i=1,i0) 8 format(1X,'CO2 saturation line gas viscosities'/ x (10(1X,E12.6))) print 9,(enthg(i),i=1,i0) 9 format(1X,'CO2 saturation line gas enthalpies'/ x (10(1X,E12.6))) endif c write(2,10) 10 format('* * ... PVT property file for CO2') write(2,11) np,nt 11 format(16I5) write(2,12) (p(i),i=1,np) 12 format(5(1X,E12.6)) write(2,12) (t(j),j=1,nt) c do4 i=1,np px=p(i)/1.e5 c do5 j=1,nt tx=t(j)+273.15 c mode=1 CALL CBCO2(Px,Tx,DENS,TCOND,VISC,BETA,H,RVAP,SIGMA,IND,MODE) c c.....flag saturation line points by assigning negative density if(ind.eq.2) then rho(j)=-dens(1) else rho(j)=dens(1) endif vis(j)=visc(1) enth(j)=h(1) 5 continue c c.....write out data write(2,12) (rho(j),j=1,nt) write(2,12) (vis(j),j=1,nt) write(2,12) (enth(j),j=1,nt) c 4 continue c c.....now for the saturation line write(2,11) np1,i0,np2 write(2,12) (t(i),i=1,i0) write(2,12) (psat(i),i=1,i0) write(2,12) (rhog(i),i=1,i0) write(2,12) (visg(i),i=1,i0) write(2,12) (enthg(i),i=1,i0) c PRINT 20,h1 20 FORMAT(A1/ X' ',131('*')/' *',129X,'*'/' *',49X,'SUMMARY OF PROGRAM UNITS', X' USED',51X,'*'/' *',129X,'*'/' ',131('*')// X' UNIT VERSION DATE COMMENTS'/ X1X,131('_')/) C ENDFILE 11 REWIND 11 C 16 CONTINUE DO17 I=1,26 17 VV(I)=' ' C READ(11,14,END=15) (VV(I),I=1,26) 14 FORMAT(26A5) PRINT 14,(VV(I),I=1,26) GOTO 16 C 15 CONTINUE PRINT 19 19 FORMAT(//1X,131('*')) c stop end c BLOCK DATA COMMON /CC/C(10,7) DATA C/ 1-1.55097,2.044841,-1.67811,0.0930354,11.38013 2,5.65977,-30.09008,5.940662,36.4241,-23.7826 3,-3.59686,5.755109,-17.8978,-36.7511,105.7493 4,111.3312,-344.1631,-43.6683,409.9414,-167.85663 5,0.5546735,27.25533,-45.02177,-232.2977,334.4779 6,707.8258,-952.2017,-711.1468,823.3502,92.16953 7,0.8054393,70.60682,-37.279,-579.96,294.483 8,1433.787,-637.589,-812.957,2*0.,-1.433238 9,88.77999,35.13756,-566.348,-108.0873,911.302 1,4*0.,-1.862086,39.4879,48.02439,-133.8665,-114.894 2,5*0.,-0.318711,9*0./ END c SUBROUTINE CBCO2(P,T,RO,CLAMB,ETA,BETA,H,RVAP,SIGMA,IND,IC) common/see/zz REAL RO(2),CLAMB(2),ETA(2),BETA(2),H(2),RMIN(2),RMAX(2) DATA R/1.889/,EPS/0.0001/ SAVE ICALL DATA ICALL/0/ C ICALL=ICALL+1 IF(ICALL.EQ.1) WRITE(11,899) c 899 FORMAT(/6X,'CBCO2 1.0 18 November 1999',6X, c 899 FORMAT(/6X,'CBCO2 1.1 28 October 2005',6X, c 899 FORMAT(/6X,'CBCO2 1.1 7 November 2005',6X, c 899 FORMAT(/6X,'CBCO2 1.1 10 February 2006',6X, c 899 FORMAT(/6X,'CBCO2 1.1 11 February 2006',6X, 899 FORMAT(/6X,'CBCO2 1.1 12 February 2006',6X, X'calculate thermophysical properties of CO2 using correlations',/ x50X,'due to V.V. Altunin (1975)'/ x47X,'as coded by V.I. Malkovsky (IGEM, Moscow, Russia)'/ x47X,'improved choice of initial search interval for bisection,'/ x50X,'and added Newtonian iteration')) c NK=1 IF(P.GT.73.92.OR.T.GT.304.2) THEN IND=4 ELSE TS=TSAT(P) IF(T.GT.TS+0.01) IND=3 IF(T.LT.TS-0.01) IND=1 IF(IC.EQ.2.OR.ABS(T-TS).LT.0.01) IND=2 ENDIF IF(IND.NE.3) THEN RMAX(1)=1.4 ELSE RMAX(1)=0.468 ENDIF IF(IND.EQ.1.OR.IND.EQ.2) RMIN(1)=0.468 c IF(IND.GE.3) RMIN(1)=P/R/T c.....11-7-05 IF(IND.GE.3) RMIN(1)=0.7*P/R/T if(ind.eq.3.and.t.le.273.15) rmax(1)=0.1 IF(IND.EQ.2) THEN P=PSAT(T) NK=2 RMAX(2)=0.468 if(t.lt.273.15) rmax(2)=0.098 RMIN(2)=P/R/T ENDIF DO 10 N=1,NK ROMAX=RMAX(N) ROMIN=RMIN(N) c it=0 1 ROCP=(ROMIN+ROMAX)/2. it=it+1 ZZ=eZ(ROCP,T) PCP=ZZ*ROCP*R*T c if(it.eq.1) then zz0=zz pcp0=pcp rocp0=rocp endif c IF(PCP.LE.P) ROMIN=ROCP IF(PCP.GT.P) ROMAX=ROCP IF(ROMAX-ROMIN.GT.EPS*ROMIN) GO TO 1 c zzf=zz pcpf=pcp rocpf=rocp c if((nk.eq.2.or.ind.eq.3).and.abs(rocp-0.468).lt.0.001) goto 12 c skip Newtonian iteration for near-critical gas iter=0 11 continue iter=iter+1 if(iter.gt.50) goto 20 ZZ=eZ(ROCP,T) p0=zz*rocp*r*t dro=max(1.e-8*rocp,1.e-7) ro1=rocp+dro zz1=ez(ro1,t) p1=zz1*ro1*r*t delro=(p-p0)*dro/(p1-p0) rocp=rocp+delro if(abs((p-p0)/p).gt.1.e-10) goto 11 c print 2,iter 2 format(' converged after',I3,' iterations') 12 continue c CALL CBO(ROCP,T,ZZ,ETA(N),CLAMB(N),BETA(N),HE) RO(N)=ROCP*1000. 10 H(N)=HE*1000. IF(NK.EQ.2) THEN RVAP=H(2)-H(1) SIGMA=0.081*(1.-T/304.2)**1.245 ENDIF RETURN c 20 continue print 3,iter,t,p,ic,ind,n,rmin(n),rmax(n),zz0,pcp0,rocp0, xromin,romax,zzf,pcpf,rocpf,zz,p0,rocp 3 format(/' after ',I3,' iterations no convergence in CBCO2', x' for T,P = ',F8.4,', ',E12.6,'; IC = MODE =',I2/ x' ind =',I2,' n =',I2,'; initially have rmin = ',E12.6, x' rmax = ',E12.6,' zz0 = ',E12.6,' pcp0 = ',E12.6, x' rocp0 = ',E12.6/7X, x' at end of bisection have romin = ',E12.6,' romax = ', xE12.6,' zzf = ',E12.6,' pcpf = ',E12.6,' rocpf = ',E12.6/ x41X,' Newtonian iteration finishes with zz = ',E12.6, x' p0 = ',E12.6,' rocp = ',E12.6) c return END c FUNCTION TSAT(P) T1=100. T2=304.2 P1=PSAT(T1) P2=PSAT(T2) 1 TS=(T1+T2)/2. PS=PSAT(TS) IF(PS.LE.P) THEN T1=TS P1=PS ELSE T2=TS P2=PS ENDIF c IF(T2-T1.GE.0.05) GO TO 1 c.....2-28-02: tighten tolerance, so that CBCO2 with T = 0.01 can c recognize two-phase points IF(T2-T1.GE.0.005) GO TO 1 TSAT=TS RETURN END c FUNCTION PSAT(T) REAL A(5) DATA A/255.895,-575.801,679.516,-402.705,95.835/ TAU=T/304.2 PSAT=-48.437/TAU DEG=1. DO 1 I=1,5 PSAT=PSAT+A(I)*DEG 1 DEG=DEG*TAU PSAT=EXP(PSAT) RETURN END c SUBROUTINE CBO(RO,T,ZZ,ETA,LAMB,BETA,H) REAL C(5,3)/0.24857,-0.373301,0.363855,-0.0639071,0. 1,0.0048949,1.22754,-0.774229,0.142507,6*0./ REAL B(5,3)/1.18764,-2.30778,2.61294,-1.283256,0.219542 1,-2.73694,4.41995,-4.00329,2.1366,-0.4021338 2,2.52043,-0.091567,-1.35345,0.376571,0./ REAL LAMB,LAMB0,DEGT(5),DEGOM(5) REAL R1/0.1889/ TAU=T/304.2 OME=RO/0.468 DEGT(1)=1. DEGOM(1)=OME DO 1 I=2,5 DEGT(I)=DEGT(I-1)*TAU 1 DEGOM(I)=DEGOM(I-1)*OME P1=SQRT(TAU) ETA0=1.E-8*P1*(2722.46-1663.46/TAU+466.9206/DEGT(3)) LAMB0=1.E-3*P1* *(57.286-78.1435/TAU+49.1871/DEGT(3)-11.5094/DEGT(4)) ETA=0. LAMB=0. DO 3 I=1,5 DO 3 J=1,3 R=DEGOM(I)/DEGT(J) ETA=ETA+C(I,J)*R 3 LAMB=LAMB+B(I,J)*R LAMB=LAMB0*EXP(LAMB) ETA=ETA0*EXP(ETA) CALL sDIF(DZDR,DZDT,RO,T,ZZ) BETA=(ZZ+T*DZDT)/(ZZ+RO*DZDR)/T 6 HT=3.7365+0.00232*(T-280.) CALL LINT(QUAD1,T,0.,RO) H1=ZZ-1.-QUAD1 H=(HT+H1)*T*R1+596.51 RETURN END c FUNCTION eZ(RO,T) COMMON /CC/C(10,7)/DEG/DEGRO(10),DEGT(7) DEGRO(1)=1. DEGT(1)=1. R1=RO-0.468 R2=304.2/T-1. DO 1 I=2,10 1 DEGRO(I)=DEGRO(I-1)*R1 DO 2 I=2,7 2 DEGT(I)=DEGT(I-1)*R2 eZ=0. DO 3 I=1,10 DO 3 J=1,7 3 eZ=eZ+C(I,J)*DEGRO(I)*DEGT(J) eZ=1.+RO*eZ RETURN END c SUBROUTINE sDIF(DZDR,DZDT,RO,T,Z) COMMON /CC/C(10,7)/DEG/DEGRO(10),DEGT(7) TAU=T/304.2 R=RO-0.468 SUMMA=0. DO 1 I=2,10 DO 1 J=1,7 1 SUMMA=SUMMA+C(I,J)*DEGRO(I-1)*(I-1)*DEGT(J) DZDR=(Z-1.)/RO+RO*SUMMA SUMMA=0. DO 2 I=1,10 DO 2 J=2,7 2 SUMMA=SUMMA+C(I,J)*(J-1)*DEGT(J-1)*DEGRO(I) DZDT=-RO/TAU**2*SUMMA/304.2 RETURN END c SUBROUTINE LINT(QUAD1,T,RO1,RO2) COMMON /CC/C(10,7)/DEG/DEGRO(10),DEGT(7) TAU=T/304.2 R1=RO1-0.468 R2=RO2-0.468 SUMM1=0. DEG1=1. DEG2=1. DO 1 I=1,10 DEG1=DEG1*R1 DEG2=DEG2*R2 DO 1 J=2,7 R=C(I,J)*(J-1)/I*(DEG2-DEG1) 1 SUMM1=SUMM1-R*DEGT(J-1) QUAD1=SUMM1/TAU RETURN END