c begin file core_ECS.f c c This file contains routines implementing an extended corresponding c states model with temperature- and density-dependent shape factors. c c contained here are: c function PHIECS (icomp,itau,idel,tau,del) c subroutine CRTECS (icomp,tc,pc,rhoc) c subroutine SETECS (nread,icomp,hcasno,href,heqn,ierr,herr) c subroutine FJ (icomp,t,d,f,dfdt,d2fdt2,dfdd,d2fdd2,d2fdtd) c subroutine HJ (icomp,t,d,h,dhdt,d2hdt2,dhdd,d2hdd2,d2hdtd) c c ====================================================================== c ====================================================================== c function PHIECS (icomp,itau,idel,tau,del) c c compute reduced Helmholtz energy or a derivative as functions c of dimensionless temperature and density for the ECS model c c inputs: c icomp--pointer specifying component (1..nc) c itau--flag specifying order of temperature derivative to calc c idel--flag specifying order of density derivative to calculate c when itau = 0 and idel = 0, compute A/RT c when itau = 0 and idel = 1, compute 1st density derivative c when itau = 1 and idel = 1, compute cross derivative c etc. c tau--dimensionless temperature (To/T) c del--dimensionless density (D/Do) c output (as function value): c phi--residual (real-gas) part of the Helmholtz energy, or one c of its derivatives (as specified by itau and idel), c in reduced form (A/RT) c c The Helmholtz energy consists of ideal and residual (real-gas) c terms; this routine calculates only the residual part. c c This function computes pure component properties only. c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 02-02-96 MM, original version c 03-19-96 MM, add dipole moment to /CCON/ c 03-22-96 MM, replace /MODEL/ with /EOSMOD/ c 09-23-96 EWL, change calls to FJ, HJ to include density dependence c 11-13-97 EWL, check for itau.ge.1 before idel.ge.1 c implicit double precision (a-h,o-z) implicit integer (i-n) character*3 hpheq,heos,hmxeos,hmodcp c parameter (ncmax=5) !max number of components in mixture parameter (n0=-ncmax,nx=ncmax) common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx) common /CCON/ wm(n0:nx),ttp(n0:nx),tnbp(n0:nx), & tcrit(n0:nx),pcrit(n0:nx),Dcrit(n0:nx),Zcrit(n0:nx), & accen(n0:nx),dipole(n0:nx) c phiecs=0.0d0 !initialize outputs and intermediate results phi01=0.0d0 phi02=0.0d0 phi10=0.0d0 phi20=0.0d0 phi11=0.0d0 c if (del.le.1.0d-10) then !trivial solution at zero density RETURN !for any and all derivatives end if c c recover the temperature and density from tau and del (need to compute c shape factors at actual temperature rather than reduced temperature) c i=icomp t=tcrit(i)/tau rho=del*Dcrit(i) call FJ (icomp,t,rho,f,dfdt,d2fdt2,dfdd,d2fdd2,d2fdtd) call HJ (icomp,t,rho,h,dhdt,d2hdt2,dhdd,d2hdd2,d2hdtd) c c find the reducing temperature to compute tau for the reference fluid c and then calculate reduced Helmholtz for the reference fluid at the c conformal temperature and density c c note: cannot use the PHIK function here because it calls this one, c and recursive calls are not generally allowed in fortran c iref=-icomp if (hmxeos(iref).eq.'BWR') then call CRTBWR (iref,tc0,pc0,Dc0) tred0=tc0 Dred0=Dc0 tau0=tred0/t*f del0=rho*h/Dred0 c compute PHI with same order as input arguments phixx=PHIBWR (iref,itau,idel,tau0,del0) c calculate additional PHIs as required if (itau.ge.1) then phi01=PHIBWR (iref,0,1,tau0,del0) if (idel.eq.1) then phi02=PHIBWR (iref,0,2,tau0,del0) phi10=PHIBWR (iref,1,0,tau0,del0) phi20=PHIBWR (iref,2,0,tau0,del0) else if (itau.eq.2) then phi10=PHIBWR (iref,1,0,tau0,del0) phi11=PHIBWR (iref,1,1,tau0,del0) phi02=PHIBWR (iref,0,2,tau0,del0) end if end if if (idel.ge.1) then phi10=PHIBWR (iref,1,0,tau0,del0) if (idel.eq.2) then phi01=PHIBWR (iref,0,1,tau0,del0) phi11=PHIBWR (iref,1,1,tau0,del0) phi20=PHIBWR (iref,2,0,tau0,del0) end if end if else if (hmxeos(iref).eq.'FEQ') then call REDFEQ (iref,tred0,Dred0) tau0=tred0/t*f del0=rho*h/Dred0 c compute PHI with same order as input arguments phixx=PHIFEQ (iref,itau,idel,tau0,del0) c calculate additional PHIs as required if (itau.ge.1) then phi01=PHIFEQ (iref,0,1,tau0,del0) if (idel.eq.1) then phi02=PHIFEQ (iref,0,2,tau0,del0) phi10=PHIFEQ (iref,1,0,tau0,del0) phi20=PHIFEQ (iref,2,0,tau0,del0) else if (itau.eq.2) then phi10=PHIFEQ (iref,1,0,tau0,del0) phi11=PHIFEQ (iref,1,1,tau0,del0) phi02=PHIFEQ (iref,0,2,tau0,del0) end if end if if (idel.ge.1) then phi10=PHIFEQ (iref,1,0,tau0,del0) if (idel.eq.2) then phi01=PHIFEQ (iref,0,1,tau0,del0) phi11=PHIFEQ (iref,1,1,tau0,del0) phi20=PHIFEQ (iref,2,0,tau0,del0) end if end if else c reference fluid not found, but no way to return an error from here phixx=-9.99d99 end if c c check if derivatives are requested c if (itau.eq.0) then if (idel.eq.0) then c compute dimensionless residual Helmholtz of reference fluid phiecs=phixx else if (idel.eq.1) then c compute 1st derivative w.r.t. del (dimensionless density) phiecs=Dcrit(i)/Dred0*phixx*(h + rho*dhdd) & +Dcrit(i)*tau0/f*phi10*dfdd else if (idel.eq.2) then c compute 2nd derivative w.r.t. del phiecs=Dcrit(i)**2*(((phixx*(h+rho*dhdd)/dred0 & +tred0/t*phi11*dfdd)*(h+rho*dhdd) & +phi01*(2.d0*dhdd+rho*d2hdd2))/dred0 & +tred0/t*((phi11/dred0*(h+rho*dhdd)+tred0/t*phi20*dfdd)*dfdd & +d2fdd2*phi10)) end if c else if (itau.eq.1) then dtlogh=dhdt/h !equal to temperature derivative of log(h) if (idel.eq.0) then c compute 1st derivative w.r.t. tau (dimensionless temperature) phiecs=tred0/tcrit(i)*(f-t*dfdt)*phixx & -t*t*del0/tcrit(i)*dtlogh*phi01 else if (idel.eq.1) then c compute cross derivative phiecs=Dcrit(i)/Dred0/tcrit(i)* & (tred0*(f-t*dfdt)*(phixx*(h + rho*dhdd) & +tred0/t*dred0*phi20*dfdd) & +tred0*phi10*Dred0*(dfdd - t*d2fdtd) & -t*t*(phi01*dhdt+rho*dhdt*(phi02/dred0*(h + rho*dhdd) & +tred0/t*phixx*dfdd) & +rho*phi01*d2hdtd)) end if c else if (itau.eq.2) then c compute 2nd derivative w.r.t. tau dtlogh=dhdt/h !equal to temperature derivative of log(h) phiecs=(1.0d0/tcrit(i)**2)* & (tred0**2*(f-t*dfdt)**2*phixx & -2.0d0*t*t*tred0*del0*(f-t*dfdt)*dtlogh*phi11 & +t**4*del0**2*dtlogh**2*phi02 & +tred0*t**3*d2fdt2*phi10 & +t**3*del0/h*(t*d2hdt2+2.0d0*dhdt)*phi01) else c invalid itau and/or idel, but no way to return an error from here phiecs=-9.99d99 end if c RETURN end !function PHIECS c c ====================================================================== c subroutine CRTECS (icomp,tc,pc,rhoc) c c returns critical parameters associated with ECS model c c N.B. these critical parameters may not necessarily be most c accurate values, but they are consistent with the ECS fit c c input: c icomp--pointer specifying component (1..nc) c outputs: c tc--critical temperature [K] c pc--critical pressure [kPa] c rhoc--molar density [mol/L] at critical point c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 02-06-96 MM, original version c 03-19-96 MM, add dipole moment to /CCON/ c implicit double precision (a-h,o-z) implicit integer (i-n) c parameter (ncmax=5) !max number of components in mixture parameter (n0=-ncmax,nx=ncmax) common /CCON/ wm(n0:nx),ttp(n0:nx),tnbp(n0:nx), & tcrit(n0:nx),pcrit(n0:nx),Dcrit(n0:nx),Zcrit(n0:nx), & accen(n0:nx),dipole(n0:nx) c tc=tcrit(icomp) pc=pcrit(icomp) rhoc=Dcrit(icomp) c RETURN end !subroutine CRTECS c c ====================================================================== c subroutine SETECS (nread,icomp,hcasno,href,heqn,ierr,herr) c c set up working arrays for use with ECS model c c inputs: c nread--file to read data from (file should have already been c opened and pointer set by subroutine SETUP) c icomp--component number in mixture (1..nc); 1 for pure fluid c hcasno--CAS number of component icomp (not required, it is here c to maintain parallel structure with SETBWR and SETFEQ) c c outputs: c href--file containing reference fluid EOS (character*80) c heqn--model ('BWR', etc) for reference fluid EOS (character*3) c ierr--error flag: 0 = successful c 1 = error (e.g. fluid not found) c herr--error string (character*255 variable if ierr<>0) c other quantities returned via arrays in commons c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 02-06-95 MM, original version c 03-19-96 MM, add dipole moment to /CCON/ c 03-21-96 MM, delete reference to /MODEL/, not needed c 03-22-96 MM, replace /CPMOD/ with /EOSMOD/ c 06-03-96 MM, add limits to /EOSLIM/ c 09-23-96 EWL, add density dependent shape factors c 09-30-96 MM, change order of density factors in .fld files c 08-19-97 MM, get rid of herr=herr, etc (avoid warning); flag nread<=0 c implicit double precision (a-h,o-z) implicit integer (i-n) parameter (ncmax=5) !max number of components in mixture parameter (n0=-ncmax,nx=ncmax) parameter (nfh=4) character*1 htab,hnull character*3 heqn character*3 hpheq,heos,hmxeos,hmodcp character*12 hcas,hcasno character*80 href character*255 herr common /NCOMP/ nc common /CHAR/ htab,hnull c commons associated with the nc components of current interest c ("working" commons and arrays) common /CCAS/ hcas(n0:nx) common /CCON/ wm(n0:nx),ttp(n0:nx),tnbp(n0:nx), & tcrit(n0:nx),pcrit(n0:nx),Dcrit(n0:nx),Zcrit(n0:nx), & accen(n0:nx),dipole(n0:nx) common /Gcnst/ R common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx) common /CFECS/ fecs(ncmax,nfh,2),hecs(ncmax,nfh,2), & fdecs(ncmax,nfh,2),hdecs(ncmax,nfh,2), & acfecs(ncmax),acfref(ncmax),Zcref(ncmax), & tmin(ncmax),tmax(ncmax),pmax(ncmax),rhomax(ncmax) common /ICFECS/ nfecs(ncmax),nhecs(ncmax), & nfdecs(ncmax),nhdecs(ncmax) c limits associated with the equation of state common /EOSLIM/ tmn(n0:nx),tmx(n0:nx),pmx(n0:nx),rhomx(n0:nx) c if (nread.le.0) then ierr=101 write (herr,1101) nread,hcasno,hnull 1101 format ('[SETECS error 101] illegal file specified; nread = ', & i4,'; CAS no. = ',a12,a1) RETURN else herr=hnull ierr=0 end if c c read data from file c write (*,*) ' SETECS--read component',icomp,' from unit',nread c iref=-icomp read (nread,*) tmin(icomp) !lower temperature limit read (nread,*) tmax(icomp) !upper temperature limit read (nread,*) pmax(icomp) !upper pressure limit read (nread,*) rhomax(icomp) !upper density limit read (nread,2003) hmodcp(icomp) !pointer to Cp0 model read (nread,2080) href !reference fluid .fld file read (nread,2003) heqn !ref fluid equation of state read (nread,*) acfref(icomp) !acc fac for reference fluid read (nread,*) Zcref(icomp) !Zc for reference fluid read (nread,*) acfecs(icomp) !accentric factor read (nread,*) tcrit(icomp) !critical temperature [K] read (nread,*) pcrit(icomp) !critical pressure [kPa] read (nread,*) Dcrit(icomp) !critical density [mol/L] c write (*,1020) tcrit(icomp),pcrit(icomp),Dcrit(icomp) c1020 format (1x,' SETECS--critical T, p, rho: ',f8.3,f8.2,f8.5) Zcrit(icomp)=pcrit(icomp)/(R*tcrit(icomp)*Dcrit(icomp)) c copy ECS parameters to general fluid constants common block c accen(icomp)=acfecs(icomp) read (nread,*) nfecs(icomp) !no. temperature terms for 'f' ESRR do 220 j=1,nfecs(icomp) read (nread,*) fecs(icomp,j,1),fecs(icomp,j,2) 220 continue read (nread,*) nfdecs(icomp) !no. density terms for 'f' ESRR if (nfdecs(icomp).ge.1) then do 224 j=1,nfdecs(icomp) read (nread,*) fdecs(icomp,j,1),fdecs(icomp,j,2) 224 continue end if read (nread,*) nhecs(icomp) !no. temperature terms for 'h' ESRR do 240 j=1,nhecs(icomp) read (nread,*) hecs(icomp,j,1),hecs(icomp,j,2) 240 continue read (nread,*) nhdecs(icomp) !no. density terms for 'h' ESRR if (nhdecs(icomp).ge.1) then do 244 j=1,nhdecs(icomp) read (nread,*) hdecs(icomp,j,1),hdecs(icomp,j,2) 244 continue end if c c copy limits into /EOSLIM/ arrays tmn(icomp)=tmin(icomp) tmx(icomp)=tmax(icomp) pmx(icomp)=pmax(icomp) rhomx(icomp)=rhomax(icomp) c RETURN 2003 format (a3) 2080 format (a80) end !subroutine SETECS c c ====================================================================== c subroutine FJ (icomp,t,d,f,dfdt,d2fdt2,dfdd,d2fdd2,d2fdtd) c c compute the equivalent substance reducing ratio (for temperature) c and its temperature derivative for use in the ECS model c c This routine implements the function of Huber & Ely, Int. J. Refrig. c 17:18-31 (1994) modified by additional general terms in (T/Tc) and by c density-dependent terms. The beta1, beta2 of Huber & Ely are c hecs(i,1,1) and hecs(i,2,1) here. c c inputs: c icomp--pointer specifying component (1..nc) c t--temperature [K] c d--density [mol/L] c outputs: c f--equivalent substance reducing ratio (for temperature) c this is equal to (Tc,icomp/Tc,ref)*theta(T) c where theta is the energy (or temperature) shape factor c dfdt--temperature derivative of f c d2fdt2--second temperature derivative of f c dfdd--density derivative of f c d2fdd2--second density derivative of f c d2fdtd--temperature-density cross derivative of f c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 02-06-96 MM, original version c 03-19-96 MM, add dipole moment to /CCON/ c 09-23-96 EWL, add density dependent shape factors c implicit double precision (a-h,o-z) implicit integer (i-n) c parameter (ncmax=5) !max number of components in mixture parameter (n0=-ncmax,nx=ncmax) parameter (nfh=4) common /CCON/ wm(n0:nx),ttp(n0:nx),tnbp(n0:nx), & tcrit(n0:nx),pcrit(n0:nx),Dcrit(n0:nx),Zcrit(n0:nx), & accen(n0:nx),dipole(n0:nx) common /CFECS/ fecs(ncmax,nfh,2),hecs(ncmax,nfh,2), & fdecs(ncmax,nfh,2),hdecs(ncmax,nfh,2), & acfecs(ncmax),acfref(ncmax),Zcref(ncmax), & tmin(ncmax),tmax(ncmax),pmax(ncmax),rhomax(ncmax) common /ICFECS/ nfecs(ncmax),nhecs(ncmax), & nfdecs(ncmax),nhdecs(ncmax) c i=icomp iref=-icomp Tr=t/tcrit(i) Dr=d/Dcrit(i) domega=acfecs(i)-acfref(i) c theta is the density shape factor theta=1.0d0+domega*(fecs(i,1,1)+fecs(i,2,1)*log(Tr)) c first and second (reduced) temperature derivatives of theta d1th=domega*fecs(i,2,1)/Tr !first Tr derivative of theta d2th=-domega*fecs(i,2,1)/(Tr*Tr) !second Tr derivative of theta c sum terms beyond those in Huber & Ely function if (nfecs(i).ge.3) then do 100 j=3,nfecs(i) theta=theta+fecs(i,j,1)*Tr**fecs(i,j,2) d1th=d1th+fecs(i,j,1)*fecs(i,j,2)*Tr**(fecs(i,j,2)-1.0d0) d2th=d2th+fecs(i,j,1)*fecs(i,j,2)*(fecs(i,j,2)-1.0d0) & *Tr**(fecs(i,j,2)-2.0d0) 100 continue end if d1thd=0 d2thd=0 if (nfdecs(i).gt.0) then do 110 j=1,nfdecs(i) theta=theta+fdecs(i,j,1)*Dr**fdecs(i,j,2) d1thd=d1thd+fdecs(i,j,1)*fdecs(i,j,2)*Dr**(fdecs(i,j,2)-1.d0) d2thd=d2thd+fdecs(i,j,1)*fdecs(i,j,2)*(fdecs(i,j,2)-1.d0) & *Dr**(fdecs(i,j,2)-2.d0) 110 continue end if ab=0 theta = theta+ab*dr*tr d1th=d1th+ab*dr d1thd=d1thd+ab*tr d1thdt=ab c !temp--set shape factors to one for debugging c theta=1.0d0 c d1th=0.0d0 c d2th=0.0d0 c convert theta and its derivatives to 'f' and its derivatives Tcrat=tcrit(icomp)/tcrit(iref) f=Tcrat*theta dfdt=d1th/tcrit(iref) d2fdt2=d2th/(tcrit(icomp)*tcrit(iref)) dfdd=d1thd*Tcrat/Dcrit(icomp) d2fdd2=d2thd*Tcrat/Dcrit(icomp)**2 d2fdtd=d1thdt*Tcrat/Dcrit(icomp)/tcrit(icomp) c RETURN end !subroutine FJ c c ====================================================================== c subroutine HJ (icomp,t,d,h,dhdt,d2hdt2,dhdd,d2hdd2,d2hdtd) c c compute the equivalent substance reducing ratio (for density) c and its temperature derivative for use in the ECS model c c This routine implements the function of Huber & Ely, Int. J. Refrig. c 17:18-31 (1994) modified by additional general terms in (T/Tc) and by c density-dependent terms. The beta1, beta2 of Huber & Ely are c hecs(i,1,1) and hecs(i,2,1) here. c c inputs: c icomp--pointer specifying component (1..nc) c t--temperature [K] c d--density [mol/L] c outputs: c h--equivalent substance reducing ratio (for density) c this is equal to (rhoc,ref/rhoc,icomp)*phi(T) c where phi is the volume (or density) shape factor c dhdt--temperature derivative of h c d2hdt2--second temperature derivative of h c dhdd--density derivative of h c d2hdd2--second density derivative of h c d2hdtd--temperature-density cross derivative of h c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 02-06-96 MM, original version c 03-19-96 MM, add dipole moment to /CCON/ c 09-23-96 EWL, add density dependent shape factors c implicit double precision (a-h,o-z) implicit integer (i-n) c parameter (ncmax=5) !max number of components in mixture parameter (n0=-ncmax,nx=ncmax) parameter (nfh=4) common /CCON/ wm(n0:nx),ttp(n0:nx),tnbp(n0:nx), & tcrit(n0:nx),pcrit(n0:nx),Dcrit(n0:nx),Zcrit(n0:nx), & accen(n0:nx),dipole(n0:nx) common /CFECS/ fecs(ncmax,nfh,2),hecs(ncmax,nfh,2), & fdecs(ncmax,nfh,2),hdecs(ncmax,nfh,2), & acfecs(ncmax),acfref(ncmax),Zcref(ncmax), & tmin(ncmax),tmax(ncmax),pmax(ncmax),rhomax(ncmax) common /ICFECS/ nfecs(ncmax),nhecs(ncmax), & nfdecs(ncmax),nhdecs(ncmax) c i=icomp iref=-icomp Tr=t/tcrit(i) Dr=d/Dcrit(i) domega=acfecs(i)-acfref(i) Zrat=Zcref(i)/Zcrit(i) c phi is the density shape factor, not to be confused with PHI=A/RT phi=Zrat*(1.0d0+domega*(hecs(i,1,1)+hecs(i,2,1)*log(Tr))) c first and second (reduced) temperature derivatives of phi d1phi=Zrat*domega*hecs(i,2,1)/Tr d2phi=-Zrat*domega*hecs(i,2,1)/(Tr*Tr) c sum terms beyond those in Huber & Ely function if (nhecs(i).ge.3) then do 100 j=3,nhecs(i) phi=phi+hecs(i,j,1)*Tr**hecs(i,j,2) d1phi=d1phi+hecs(i,j,1)*hecs(i,j,2)*Tr**(hecs(i,j,2)-1.0d0) d2phi=d2phi+hecs(i,j,1)*hecs(i,j,2)*(hecs(i,j,2)-1.0d0) & *Tr**(hecs(i,j,2)-2.0d0) 100 continue end if c d1phid=0.0d0 d2phid=0.0d0 if (nhdecs(i).gt.0) then do 110 j=1,nhdecs(i) phi=phi+hdecs(i,j,1)*Dr**hdecs(i,j,2) d1phid=d1phid+hdecs(i,j,1)*hdecs(i,j,2)*Dr**(hdecs(i,j,2)-1.d0) d2phid=d2phid+hdecs(i,j,1)*hdecs(i,j,2)*(hdecs(i,j,2)-1.d0) & *Dr**(hdecs(i,j,2)-2.d0) 110 continue end if c ab=0.0d0 phi=phi+ab*dr*tr d1phi=d1phi+ab*dr d1phid=d1phid+ab*tr d1phidt=ab c c !temp--set shape factors to one for debugging c phi=1.0d0 c d1phi=0.0d0 c d2phi=0.0d0 c convert phi and its derivatives to 'h' and its derivatives Dcrat=Dcrit(iref)/Dcrit(icomp) h=Dcrat*phi dhdt=d1phi*Dcrat/tcrit(icomp) d2hdt2=d2phi*Dcrat/tcrit(icomp)**2 dhdd=d1phid*Dcrat/Dcrit(icomp) d2hdd2=d2phid*Dcrat/Dcrit(icomp)**2 d2hdtd=d1phidt*Dcrat/Dcrit(icomp)/tcrit(icomp) c RETURN end !subroutine HJ c c c 1 2 3 4 5 6 7 c23456789012345678901234567890123456789012345678901234567890123456789012 c c ====================================================================== c end file core_ECS.f c ======================================================================