program test_test implicit double precision (a-h,o-z) c c program to generate tables of saturation properties for a mixture c as test and example of use of Refprop6 fortran routines c c written by Mark O. McLinden, c Physical and Chemical Properties Division c NIST, Boulder, CO c August 11, 1997 c parameter (ncmax=5) !max number of components in mixture parameter (n0=-ncmax,nx=ncmax) character*1 hdelim character*3 hdir c character*1 hdir character*3 hrf,heqn,hph,hmxeos(ncmax),htrn,htrnk(ncmax) character*12 hcas(ncmax),hname(ncmax) character*24 htest character*80 hf(ncmax),hfmix,hfread,hn80 character*255 herr dimension x(ncmax),xl(ncmax),xv(ncmax) c c file delimiter; select appropriate character c data hdelim/':'/ !Macintosh data hdelim/'\'/ !DOS c data hdelim/'/'/ !Unix c fluids directory is parallel to fortran directory hdir='..'//hdelim !select "character*3 hdir" above c fluids directory is under fortran directory c hdir=hdelim !select "character*1 hdir" above c c read parameters write (*,*) 'input number of components (-nc to accept defaults)' read (*,*) nc if (nc.ge.1) then write (*,*) 'input file name(s) (R123.fld, R134a.fld, ', & 'R152a.fld, R32.fld, R125.fld, etc.)' do 80 i=1,nc read (*,2080) hfread hf(i)=hdir//'fluids'//hdelim//hfread 80 continue 2080 format (a80) heqn='NBS' htrn='NBS' hfmix=hdir//'fluids'//hdelim//'HMX.bnc' else c default inputs nc=abs(nc) heqn='HMX' hf(1)=hdir//'fluids'//hdelim//'R32.fld' hf(2)=hdir//'fluids'//hdelim//'R125.fld' hf(3)=hdir//'fluids'//hdelim//'R134a.fld' hfmix=hdir//'fluids'//hdelim//'HMX.bnc' hmxeos(2)='NBS' hmxeos(3)='NBS' htrn='NBS' end if c set each of the component models to the default choice do 88 i=1,nc hmxeos(i)='NBS' htrnk(i)=htrn 88 continue c c define fluid c c initialize components, reference states, etc. hrf='DEF' !use default reference state c call to SETUP is required, calls to SETMOD are optional c since above code specifies the default choice for all the models, c the calls to SETMOD are completely redundant, but they are included c here to illustrate how alternate models would be specified call SETMOD (nc,'EOS',heqn,hmxeos,ierr,herr) call SETMOD (nc,'ETA',htrn,htrnk,ierr,herr) call SETMOD (nc,'TCX',htrn,htrnk,ierr,herr) call SETUP (nc,hf,hfmix,hrf,ierr,herr) call ERRMSG (ierr,herr) write (*,*) do 90 i=1,ncmax c get component names (should return 'not defined' for i > nc) call NAME (i,hname(i),hn80,hcas(i)) 90 continue do 92 i=1,nc c get component critical parameters call INFO (icomp,wm,ttp,tnbp,tcrit,pcrit,Dcrit,Zc,acf,dip,Rgas) write (*,1008) i,hcas(i),hname(i),tcrit,pcrit,Dcrit 92 continue 1008 format ( & 1x,'Critical T (K), P (kPa), rho (mol/L) for component ',i2, & ' (CAS # ',a12,'; ',a12,'): ',f8.3,f18.10,f11.6) c c default compositions if (nc.eq.1) then x(1)=1.0d0 else if (nc.eq.2) then x(1)=0.50d0 x(2)=0.50d0 else x(1)=0.333 x(2)=0.333 x(3)=0.334 end if c c imass=0 !chose this and the next line for molar basis c wmm=1.0d0 imass=1 !chose this and the next line for mass basis wmm=WMOL(x) !get the molecular weight of mixture c c write header c htest='SATT (bubble) test for:' if (imass.eq.1) then write (*,1001) htest,(hname(j),j=1,3) else if (imass.eq.0) then write (*,1004) htest,(hname(j),j=1,3) end if c c get mixture critical parameters (for use in setting upper limit of table) call CRITP (x,tcrit,pcrit,Dcrit,ierr,herr) c c calculate at NBP pl=101.325d0 call SATP (pl,x,2,tdew,Djunk,rhov,xl,xv,ierr,herr) call ERRMSG (ierr,herr) call THERM (tdew,rhov,x,ptherm,e,hv,sv,cv,cp,w,hjt) call SATP (pl,x,1,tbub,rhol,Djunk,xl,xv,ierr,herr) call ERRMSG (ierr,herr) call THERM (tbub,rhol,x,ptherm,e,h,s,cv,cp,w,hjt) call SURFT (tbub,rhol,xl,sigma,ierr,herr) rl=rhol*wmm rv=rhov*wmm h=h/wmm s=s/wmm hv=hv/wmm sv=sv/wmm pl=pl*0.001 if (imass.eq.1) then tu=tbub-273.15 write (*,1011) tu,pl,rl,rv,sigma, & xl(1),xl(2),xl(3),xv(1),xv(2),xv(3), & h,s else if (imass.eq.0) then tu=tbub write (*,1040) tu,pl,rl,rv,sigma, & xl(1),xl(2),xl(3),xv(1),xv(2),xv(3) end if c c compute properties as a function of temperature c delt=20.0 if (imass.eq.1) then t0=tbub+2.0*delt-mod(tbub,delt)+13.15 else t0=tbub+2.0*delt-mod(tbub,delt) end if tf=tcrit-20.0 c do 800 t=t0,tf,delt call SATT (t,x,1,pl,rhol,rhov,xl,xv,ierr,herr) call ERRMSG (ierr,herr) call THERM (t,rhol,x,ptherm,e,h,s,cv,cp,w,hjt) call SURFT (t,rhol,xl,sigma,ierr,herr) call ERRMSG (ierr,herr) call TRNPRP (t,rhov,xv,eta,tcx,ierr,herr) rl=rhol*wmm rv=rhov*wmm h=h/wmm s=s/wmm pl=pl*0.001 !convert from kPa to MPa if (imass.eq.1) then tu=t-273.15 !convert K to Celsius write (*,1010) tu,pl,rl,rv,sigma, & xl(1),xl(2),xl(3),xv(1),xv(2),xv(3),h,s else if (imass.eq.0) then tu=t write (*,1040) tu,pl,rl,rv,sigma, & xl(1),xl(2),xl(3),xv(1),xv(2),xv(3) end if 800 CONTINUE c 999 stop 1001 FORMAT (///' ', & a24,'(',a12,')',2(' + (',a12,')')// & ' Temp Pressure Density sur ten', & ' liquid compositions vapor compositions h', & ' s'/ & ' (øC) (MPa) (kg/m**3) (N/m) ', & ' (mol frac) (mol frac) (kJ/kg)', & ' (kJ/kg-K)'/ & ' liq vap '/) 1004 FORMAT (///' ', & a24,'(',a12,')',2(' + (',a12,')')// & ' Temp Pressure Ptherm Density sur ten', & ' liquid compositions vapor compositions'/ & ' (K) (MPa) (mol/L) (N/m) ', & ' (mol frac) (mol frac) '/ & ' liq vap '/) 1010 FORMAT ( & f8.3,d20.10,f8.1,f9.4,f9.5,3f7.4,2x,3f7.4,f7.1,f10.6) c T P Dl Dv s.t. xl's xv's h s 1011 FORMAT ( & f8.3,d20.10,f8.1,f9.4,f9.5,3f7.4,2x,3f7.4,f7.1,f10.6) c T P Dl Dv s.t. xl's xv's h s 1040 FORMAT ( & f8.3,d20.10,f8.3,f9.5,f9.5,5f7.4,f7.4) c T P Dl Dv s.t. x's end c c 1 2 3 4 5 6 7 c23456789012345678901234567890123456789012345678901234567890123456789012