c begin file core_FEQ.f c c This file contains the routines implementing the Helmholtz form of c the pure fluid equation of state (the so-called "Fundamental Eqn"). c c contained here are: c function PHIFEQ (icomp,itau,idel,tau,del) c subroutine CRTFEQ (icomp,tcrit,pcrit,Dcrit) c subroutine REDFEQ (icomp,tred,Dred) c subroutine SETFEQ (nread,icomp,hcasno,ierr,herr) c block data BDFEQ c c ====================================================================== c ====================================================================== c function PHIFEQ (icomp,itau,idel,tau,del) c c compute reduced Helmholtz energy or a derivative as functions c of dimensionless temperature and density for the Helmholtz-explicit c equation of state c c based on Tillner-Roth & Baehr (1994), JPCRD 23:657-729 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 N.B. The reducing parameters To and Do are often, but not c necessarily, equal to the critical temperature and density. 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 07-26-95 MM, original version c 08-14-95 MM, put saved variables into common (rather than save stmt) c 11-01-95 MM, increase parameter mxtrm to 52 (to accomodate steam) c 11-29-95 MM, variable lower limit on coefficient/constant arrays c to accomodate ECS reference fluid c 01-16-96 MM, implicit integer (i-n); (include L) c 02-27-96 MM, parameter n0=-ncmax to accomodate ECS-thermo model c 04-18-96 MM, apply tolerance to del=delsav + fix bug (phi not saved) c 10-09-96 MM, fix potential underflow on exponential term c 01-03-97 MM, use sorted powers of tau, del, stored in tpower, etc. c new arrays added to /WNTFEQ/,/WLFFEQ/,/WCFFEQ/ c 01-22-97 EWL, add critical-region terms of Wagner to Helmholtz equation. c increase parameter mxtrm to 56 (to accomodate steam) c 01-31-97 MM, modify commons associated with critical-region terms c 06-03-97 EWL, add third derivative of Helmholtz energy with respect to density. c 07-02-97 EWL, add parameter for g in exp(-g*d^l) for fluids where g<>1 c modify /WCFFEQ/ to accomodate g parameter 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 (mxtrm=56,mxcrt=6) c numbers of terms associated with the "normal" Helmholtz function c (plus numbers of unique powers of temperature, density, etc.), c the critical-region terms of Wagner, plus spare for future use common /WNTFEQ/ ntermf(n0:nx),ncoeff(n0:nx), & ntp(n0:nx),ndp(n0:nx),nlp(n0:nx), & ncrt(n0:nx),ncfcrt(n0:nx), & nspare(n0:nx),ncfsp(n0:nx) common /WLFFEQ/ li(n0:nx,mxtrm),lpower(n0:nx,mxtrm), & itp(n0:nx,mxtrm),idp(n0:nx,mxtrm),ilp(n0:nx,mxtrm) common /WCFFEQ/ a(n0:nx,mxtrm),ti(n0:nx,mxtrm),di(n0:nx,mxtrm), & gi(n0:nx,mxtrm),tpower(n0:nx,mxtrm), & dpower(n0:nx,mxtrm),gi2(n0:nx,mxtrm), & rho0(n0:nx),t0(n0:nx), & pc(n0:nx),rhoc(n0:nx),tc(n0:nx), & wmf(n0:nx),Rfeq(n0:nx), & pmin(n0:nx),rhotp(n0:nx),tmin(n0:nx), & tmax(n0:nx),pmax(n0:nx) c parameters associated with the critical-region terms of Wagner common /WCFFQ2/ alpha(n0:nx,mxcrt),beta(n0:nx,mxcrt), & gamma(n0:nx,mxcrt),delta(n0:nx,mxcrt), & eta(n0:nx,mxcrt),eid(n0:nx,mxcrt),eit(n0:nx,mxcrt) c save information for possible use on next call to function c save phisav(n0:nx,0:mxtrm),delsav(n0:nx),tausav(n0:nx) common /FEQSAV/ phisav(n0:nx,0:mxtrm),delsav(n0:nx),tausav(n0:nx), & taup(n0:nx,mxtrm),delp(n0:nx,mxtrm), & delli(n0:nx,mxtrm),expdel(n0:nx,mxtrm) common /CRTSAV/ delb(mxcrt),taua(mxcrt),txp(mxcrt),hxp(mxcrt), & ext(mxcrt),extd(mxcrt),extt(mxcrt),extdt(mxcrt), & extt2(mxcrt),extd2(mxcrt) c dimension phi01(mxtrm),phi10(mxtrm) c c write (*,*) ' PHIFEQ--enter w/ tau,del = ',tau,del phifeq=0.0d0 if (del.le.1.0d-10) then !trivial solution at zero density RETURN !for any and all derivatives end if c if (abs(tau-tausav(icomp)).lt.1.0d-12 .and. & abs(del-delsav(icomp)).lt.1.0d-12) then c retrieve value from previous call phifeq=phisav(icomp,0) ! zero element stores summation c write (*,*) ' PHIFEQ--using phisav itau,idel,tau:',itau,idel,tau else c otherwise, compute new values and save for possible future use c first compute needed powers of tau and del (and save for future use) if (abs(tau-tausav(icomp)).gt.1.0d-12) then tausav(icomp)=tau do 20 j=1,ntp(icomp) taup(icomp,j)=tau**tpower(icomp,j) c write (*,*) ' PHIFEQ--tau,tpower: ',tau,tpower(icomp,j) 20 continue end if if (abs(del-delsav(icomp)).gt.1.0d-12) then delsav(icomp)=del do 40 j=1,ndp(icomp) delp(icomp,j)=del**dpower(icomp,j) c write (*,*) ' PHIFEQ--del,dpower: ',del,dpower(icomp,j) 40 continue do 60 j=1,nlp(icomp) if (lpower(icomp,j).eq.0) then delli(icomp,j)=0.0d0 expdel(icomp,j)=1.0d0 else delli(icomp,j)=gi2(icomp,j)*del**lpower(icomp,j) c check that exponential term will not result in underflow c this is a potential problem when exponent and del are both c large; e.g. del >= 3 and li >= 6; affected terms apply only c at near-critical densities and go to zero at high densities if (delli(icomp,j).gt.500.0d0) then expdel(icomp,j)=exp(-500.0d0) else expdel(icomp,j)=exp(-delli(icomp,j)) end if end if c write (*,*) ' PHIFEQ--del,lpower: ',del,lpower(icomp,j) 60 continue end if c c check for presence of critical-region terms of Wagner (e.g. steam) c if (ncrt(icomp).gt.0) then do 80 i=1,ncrt(icomp) k=i+ntermf(icomp) ext(i)=1.0d0 extd(i)=0.0d0 extd2(i)=0.0d0 extt(i)=0.0d0 extt2(i)=0.0d0 delb(i)=del-delta(icomp,i) if (abs(beta(icomp,i)).gt.1.0d0) then if (abs(delb(i)).lt.1.0d-10) then hxp(i)=0.0d0 delb(i)=1.0d0 else hxp(i)=alpha(icomp,i)*delb(i)**li(icomp,k) endif taua(i)=tau-gamma(icomp,i) if (abs(taua(i)).lt.1.0d-10) then txp(i)=0.0d0 taua(i)=1.0d0 else txp(i)=beta(icomp,i)*taua(i)**li(icomp,k) endif else del1=del-1.0d0 tau1=1.0d0-tau if (ABS(del1).lt.1.0d-10) del1=SIGN(1.0d-10,del1) if (ABS(tau1).lt.1.0d-10) tau1=SIGN(1.0d-10,tau1) delb(i)=del1 taua(i)=tau1 tdg=tau1+gamma(icomp,i)*ABS(del1)**(1.0d0/beta(icomp,i)) ext(i)=(tdg**2+eid(icomp,i)*(del1**2)**eit(icomp,i)) s=SIGN(1.0d0,del1) extd(i)=2.0d0*tdg*gamma(icomp,i)/beta(icomp,i) & *ABS(del1)**(1.0d0/beta(icomp,i)-1.0d0)*s & +2.0d0*eid(icomp,i)*eit(icomp,i) & *ABS(del1)**(2.0d0*eit(icomp,i)-1.0d0)*s extd2(i)=alpha(icomp,i)*((alpha(icomp,i)-1.0d0) & *ext(i)**(alpha(icomp,i)-2.0d0)*extd(i)**2 & +ext(i)**(alpha(icomp,i)-1.0d0) & *(2.0d0*tdg*gamma(icomp,i)/beta(icomp,i) & *(1.0d0/beta(icomp,i)-1.0d0) & *ABS(del1)**(1.0d0/beta(icomp,i)-2.0d0) & +2.0d0*(gamma(icomp,i)/beta(icomp,i) & *ABS(del1)**(1.0d0/beta(icomp,i)-1.0d0))**2 & +2.0d0*eid(icomp,i)*eit(icomp,i) & *(2.0d0*eit(icomp,i)-1.0d0) & *ABS(del1)**(2.0d0*eit(icomp,i)-2.0d0))) extdt(i)=-2.0d0*gamma(icomp,i)/beta(icomp,i) & * ABS(del1)**(1.0d0/beta(icomp,i)-1.0d0)*s extt(i)=-2.0d0*tdg extt2(i)=alpha(icomp,i)*((alpha(icomp,i)-1.0d0) & *ext(i)**(alpha(icomp,i)-2.0d0)*extt(i)**2 & +2.0d0*ext(i)**(alpha(icomp,i)-1.0d0)) extdt(i)=alpha(icomp,i)*((alpha(icomp,i)-1.0d0) & *ext(i)**(alpha(icomp,i)-2.0d0)*extt(i)*extd(i) & +ext(i)**(alpha(icomp,i)-1.0d0)*extdt(i)) extd(i)=alpha(icomp,i) & *ext(i)**(alpha(icomp,i)-1.0d0)*extd(i) extt(i)=alpha(icomp,i) & *ext(i)**(alpha(icomp,i)-1.0d0)*extt(i) ext(i)=ext(i)**alpha(icomp,i) hxp(i)=-delta(icomp,i)*del1**li(icomp,k) txp(i)=-eta(icomp,i)*tau1**li(icomp,k) taua(i)=-taua(i) endif 80 continue endif c c end critical-region terms c phisum=0.0d0 do 100 k=1,ntermf(icomp) phisav(icomp,k)=a(icomp,k)*taup(icomp,itp(icomp,k)) & *delp(icomp,idp(icomp,k))*expdel(icomp,ilp(icomp,k)) phisum=phisum+phisav(icomp,k) c write (*,1010) k,phisav(icomp,k),phisum c1010 format (1x,i3,2d30.20) !write out each term for debugging 100 continue if (ncrt(icomp).gt.0) then do 110 i=1,ncrt(icomp) k=i+ntermf(icomp) phisav(icomp,k)=a(icomp,k)*taup(icomp,itp(icomp,k)) & *delp(icomp,idp(icomp,k))*exp(hxp(i)+txp(i))*ext(i) phisum=phisum+phisav(icomp,k) 110 continue endif phifeq=phisum phisav(icomp,0)=phisum end if c c write (*,1012) tau,del,itau,idel,phifeq c1012 format (1x,' PHIFEQ--tau,del,itau,idel,phi00: ',2f10.5,2i4,e22.12) c c check if derivatives are requested, calculations make use of fact c that terms in derivative summations are very similar to A/RT terms c if (idel.eq.1) then c compute derivative w.r.t. del (dimensionless density) c save individual terms for possible use in cross derivative phisum=0.0d0 do 120 k=1,ntermf(icomp) if (li(icomp,k).eq.0) then phi01(k)=phisav(icomp,k)*di(icomp,k)/del else c phi01(k)=phisav(icomp,k)/del c & *(di(icomp,k)-li(icomp,k)*del**li(icomp,k)) phi01(k)=phisav(icomp,k)/del & *(di(icomp,k)-li(icomp,k)*delli(icomp,ilp(icomp,k))) endif phisum=phisum+phi01(k) 120 continue c check for presence of critical-region terms if (ncrt(icomp).gt.0) then do 125 i=1,ncrt(icomp) k=i+ntermf(icomp) phi01(k)=phisav(icomp,k)/ext(i) & *(ext(i)*(di(icomp,k)/del & +li(icomp,k)*hxp(i)/delb(i))+extd(i)) phisum=phisum+phi01(k) 125 continue endif phifeq=phisum else if (idel.eq.2) then c compute 2nd derivative w.r.t. del (dimensionless density) phisum=0.0d0 do 130 k=1,ntermf(icomp) dik=di(icomp,k) if (li(icomp,k).eq.0) then phisum=phisum+phisav(icomp,k)*dik*(dik-1.0d0)/del**2 else c dell=li(icomp,k)*del**li(icomp,k) dell=li(icomp,k)*delli(icomp,ilp(icomp,k)) phisum=phisum+phisav(icomp,k)/del**2 c Tillner-Roth has "+dell" as final term, I get "-dell." c This is apparantly a misprint in the T-R thesis I was working from; c consulted with T-R on this 03-07-96 and we agree on "-dell." [MM] c & *(dik**2-dik-dell*(2.0d0*dik+li(icomp,k)-1.0d0+dell)) & *(dik**2-dik-dell*(2.0d0*dik+li(icomp,k)-1.0d0-dell)) end if 130 continue c check for presence of critical-region terms if (ncrt(icomp).gt.0) then do 135 i=1,ncrt(icomp) k=i+ntermf(icomp) dik=di(icomp,k) dell=li(icomp,k)*hxp(i)/delb(i) phisum=phisum+phisav(icomp,k)/del**2/ext(i) & *(del**2*(extd2(i)+extd(i)*(dik/del+dell) & +ext(i)*(-dik/del**2+dell*(li(icomp,k)-1.0d0)/delb(i))) & +(extd(i)+ext(i)*(dik/del+dell)) & *(dik*del+del**2*dell)) 135 continue endif phifeq=phisum else if (idel.eq.3) then c compute 2nd derivative w.r.t. del (dimensionless density) phisum=0.0d0 do 138 k=1,ntermf(icomp) dik=di(icomp,k) if (li(icomp,k).eq.0) then phisum=phisum+phisav(icomp,k)*dik*(dik-1.0d0) & *(dik-2.0d0)/del**3 else dell=li(icomp,k)*delli(icomp,ilp(icomp,k)) l=li(icomp,k) phisum=phisum+phisav(icomp,k)/del**3 & *(2.0d0*dik-3.0d0*dik*dik+dik**3+dell*(-2.0d0+3.0d0*l-l**2 & +6.0d0*dik-3.0d0*dik*l-3.0d0*dik**2 & +dell*(-3.0d0+3.0d0*l+3.0d0*dik-dell))) end if 138 continue phifeq=phisum end if c c if (itau.eq.1) then c compute derivative w.r.t. tau (dimensionless temperature) c save individual terms for possible use in cross derivative phisum=0.0d0 do 140 k=1,ntermf(icomp) phi10(k)=phisav(icomp,k)*ti(icomp,k)/tau phisum=phisum+phi10(k) 140 continue c check for presence of critical-region terms if (ncrt(icomp).gt.0) then do 145 i=1,ncrt(icomp) k=i+ntermf(icomp) phi10(k)=phisav(icomp,k)/tau*(ti(icomp,k) & +tau*li(icomp,k)*txp(i)/taua(i)+tau*extt(i)/ext(i)) phisum=phisum+phi10(k) 145 continue endif phifeq=phisum else if (itau.eq.2) then c compute 2nd derivative w.r.t. tau (dimensionless temperature) phisum=0.0d0 do 150 k=1,ntermf(icomp) tik=ti(icomp,k) phisum=phisum+phisav(icomp,k)*tik*(tik-1.0d0)/tau**2 150 continue c check for presence of critical-region terms if (ncrt(icomp).gt.0) then do 155 i=1,ncrt(icomp) k=i+ntermf(icomp) tik=ti(icomp,k) tauu=li(icomp,k)*txp(i)/taua(i) phisum=phisum+phisav(icomp,k)/tau**2*(tau**2 & *(extt2(i)+extt(i)*(tik/tau+tauu) & +ext(i)*(-tik/tau**2+tauu*(li(icomp,k)-1.0d0)/taua(i))) & +(extt(i)+ext(i)*(tik/tau+tauu)) & *(tik*tau+tau**2*tauu))/ext(i) 155 continue endif phifeq=phisum end if c c if ((itau.eq.1) .and. (idel.eq.1)) then c compute cross derivative using terms from 1st derivatives phisum=0.0d0 do 160 k=1,ntermf(icomp) if (phisav(icomp,k).ne.0) & phisum=phisum+phi10(k)*phi01(k)/phisav(icomp,k) 160 continue c check for presence of critical-region terms if (ncrt(icomp).gt.0) then do 165 i=1,ncrt(icomp) k=i+ntermf(icomp) phisum=phisum+(ti(icomp,k)+tau*li(icomp,k) & *txp(i)/taua(i))*phi01(k)/tau & +phisav(icomp,k)/ext(i)*(del*tau*extdt(i) & +tau*extt(i)*(di(icomp,k) & +del*li(icomp,k)*hxp(i)/delb(i)))/del/tau 165 continue endif phifeq=phisum end if c c write (*,1021) tau,del,itau,idel,phifeq c1021 format (1x,' PHIFEQ--tau,del,itau,idel,phixx: ',2f10.5,2i4,e22.12) c RETURN end !function PHIFEQ c c ====================================================================== c subroutine CRTFEQ (icomp,tcrit,pcrit,Dcrit) c c returns critical parameters associated with Fundamental EOS c c N.B. these critical parameters may not necessarily be most c accurate values, but they are consistent with the EOS fit; c neither are they always equal to the reducing parameters c c input: c icomp--pointer specifying component (1..nc) c outputs: c tcrit--critical temperature (K) c pcrit--critical pressure (kPa) c Dcrit--molar density (mol/L) at critical point c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 08-03-95 MM, original version c 11-01-95 MM, increase parameter mxtrm to 52 (to accomodate steam) c 11-29-95 MM, variable lower limit on coefficient/constant arrays c to accomodate ECS reference fluid c 01-16-96 MM, implicit integer (i-n); (include L) c 02-27-96 MM, parameter n0=-ncmax to accomodate ECS-thermo model c 01-03-97 MM, new arrays added to /WCFFEQ/ (assoc w/ sorting of powers) c 07-02-97 EWL, modify /WCFFEQ/ to accomodate g parameter c implicit double precision (a-h,o-z) implicit integer (i-n) c c parameter (mxtrm=56) parameter (ncmax=5) !max number of components in mixture parameter (n0=-ncmax,nx=ncmax) common /WCFFEQ/ a(n0:nx,mxtrm),ti(n0:nx,mxtrm),di(n0:nx,mxtrm), & gi(n0:nx,mxtrm),tpower(n0:nx,mxtrm), & dpower(n0:nx,mxtrm),gi2(n0:nx,mxtrm), & rho0(n0:nx),t0(n0:nx), & pc(n0:nx),rhoc(n0:nx),tc(n0:nx), & wmf(n0:nx),Rfeq(n0:nx), & pmin(n0:nx),rhotp(n0:nx),tmin(n0:nx), & tmax(n0:nx),pmax(n0:nx) c tcrit=tc(icomp) pcrit=pc(icomp) Dcrit=rhoc(icomp) c RETURN end !subroutine CRTFEQ c c ====================================================================== c subroutine REDFEQ (icomp,tred,Dred) c c returns reducing parameters associated with Fundamental EOS; c used to calculate the 'tau' and 'del' which are the independent c variables in the EOS c c N.B. The reducing parameters are often, but not always, equal c to the critical temperature and density. c c input: c icomp--component number in mixture (1..nc); 1 for pure fluid c outputs: c tred--reducing temperature (K) c Dred--reducing molar density (mol/L) c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 08-03-95 MM, original version c 11-01-95 MM, increase parameter mxtrm to 52 (to accomodate steam) c 11-29-95 MM, variable lower limit on coefficient/constant arrays c to accomodate ECS reference fluid c 01-16-96 MM, implicit integer (i-n); (include L) c 02-27-96 MM, parameter n0=-ncmax to accomodate ECS-thermo model c 01-03-97 MM, new arrays added to /WCFFEQ/ (assoc w/ sorting of powers) c 07-02-97 EWL, modify /WCFFEQ/ to accomodate g parameter c implicit double precision (a-h,o-z) implicit integer (i-n) c c parameter (mxtrm=56) parameter (ncmax=5) !max number of components in mixture parameter (n0=-ncmax,nx=ncmax) common /WCFFEQ/ a(n0:nx,mxtrm),ti(n0:nx,mxtrm),di(n0:nx,mxtrm), & gi(n0:nx,mxtrm),tpower(n0:nx,mxtrm), & dpower(n0:nx,mxtrm),gi2(n0:nx,mxtrm), & rho0(n0:nx),t0(n0:nx), & pc(n0:nx),rhoc(n0:nx),tc(n0:nx), & wmf(n0:nx),Rfeq(n0:nx), & pmin(n0:nx),rhotp(n0:nx),tmin(n0:nx), & tmax(n0:nx),pmax(n0:nx) c c write (*,*) ' REDFEQ, i,t0,rho0: ',icomp,t0(icomp),rho0(icomp) tred=t0(icomp) Dred=rho0(icomp) c RETURN end !subroutine REDFEQ c c ====================================================================== c subroutine SETFEQ (nread,icomp,hcasno,ierr,herr) c c set up working arrays for use with Fundamental equation of state c c inputs: c nread--file to read data from c <= 0 get data from block data c >0 read from logical unit nread (file should have already c been 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 req'd if reading from file) c c outputs: 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 common /WCFFEQ/ c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 07-20-95 MM, original version c 09-13-95 MM, add ierr, herr to argument list c 10-05-95 MM, adapt to file input, add nread to argument list c 11-01-95 MM, increase parameter mxtrm to 52 (to accomodate steam) c 11-02-95 MM, add common /CCAS/ to access CAS numbers c 11-29-95 MM, variable lower limit on coefficient/constant arrays c to accomodate ECS reference fluid c 01-16-96 MM, implicit integer (i-n); (include L) c 02-27-96 MM, parameter n0=-ncmax to accomodate ECS-thermo model c add Zcrit to common /CCON/ c 03-19-96 MM, add dipole moment to /CCON/ and /MSCFEQ/ c 06-03-96 MM, add limits to /EOSLIM/, reduce mxfeq from 20 to 2 c 01-03-97 MM, sort powers of tau, del, store in tpower, etc. c new arrays added to /WNTFEQ/,/WLFFEQ/,/WCFFEQ/ c 01-22-97 EWL, add variables for critical region for methane and water. c 01-31-97 MM, modify commons associated with critical-region terms c 02-06-97 MM, if steam, reset gas constant to fluid-specific value c 02-07-97 MM, fix bug in reading of critical terms c 05-27-97 MM, if nc = 1, set R to fluid-specific value c 07-02-97 EWL, add parameter for g in exp(-g*d^l) for fluids where g<>1 c modify /WCFFEQ/ to accomodate g parameter c 11-13-97 MM, (re)initialize contents of /FEQSAV/ when a new fluid is read in c 02-11-98 MM, store rho at triple point separate from rhomax c implicit double precision (a-h,o-z) implicit integer (i-n) parameter (mxfeq=2) !max number of FEQ EOS in block data parameter (mxtrm=56,mxcrt=6) parameter (ncmax=5) !max number of components in mixture parameter (n0=-ncmax,nx=ncmax) character*1 htab,hnull character*3 hcpfeq character*3 hpheq,heos,hmxeos,hmodcp character*12 hcas,hcasa,hcasno character*255 herr c character*1 dummy common /NCOMP/ nc common /CHAR/ htab,hnull common /Gcnst/ R c commons associated with the mxfeq fluids with FEQ equations stored c in block data BDFEQ common /CASFEQ/ hcasa(mxfeq) common /CPMFEQ/ hcpfeq(mxfeq) common /NTFEQ/ ntrmfa(mxfeq) common /LFFEQ/ lia(mxfeq,mxtrm) common /CFFEQ/ aa(mxfeq,mxtrm),dia(mxfeq,mxtrm),tia(mxfeq,mxtrm), & rho0a(mxfeq),t0a(mxfeq), & pca(mxfeq),rhoca(mxfeq),tca(mxfeq), & wmfa(mxfeq),Rfeqa(mxfeq), & pmina(mxfeq),rhomxa(mxfeq),tmina(mxfeq), & tmaxa(mxfeq),pmaxa(mxfeq) common /MSCFEQ/ ttpf(mxfeq),tnbpf(mxfeq),accenf(mxfeq),dipm(mxfeq) 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) c common /CPMOD/ hmodcp(n0:nx) common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx) common /WNTFEQ/ ntermf(n0:nx),ncoeff(n0:nx), & ntp(n0:nx),ndp(n0:nx),nlp(n0:nx), & ncrt(n0:nx),ncfcrt(n0:nx), & nspare(n0:nx),ncfsp(n0:nx) common /WLFFEQ/ li(n0:nx,mxtrm),lpower(n0:nx,mxtrm), & itp(n0:nx,mxtrm),idp(n0:nx,mxtrm),ilp(n0:nx,mxtrm) common /WCFFEQ/ a(n0:nx,mxtrm),ti(n0:nx,mxtrm),di(n0:nx,mxtrm), & gi(n0:nx,mxtrm),tpower(n0:nx,mxtrm), & dpower(n0:nx,mxtrm),gi2(n0:nx,mxtrm), & rho0(n0:nx),t0(n0:nx), & pc(n0:nx),rhoc(n0:nx),tc(n0:nx), & wmf(n0:nx),Rfeq(n0:nx), & pmin(n0:nx),rhotp(n0:nx),tmin(n0:nx), & tmax(n0:nx),pmax(n0:nx) c parameters associated with the critical-region terms of Wagner common /WCFFQ2/ alpha(n0:nx,mxcrt),beta(n0:nx,mxcrt), & gamma(n0:nx,mxcrt),delta(n0:nx,mxcrt), & eta(n0:nx,mxcrt),eid(n0:nx,mxcrt),eit(n0:nx,mxcrt) c limits associated with the equation of state common /EOSLIM/ tmn(n0:nx),tmx(n0:nx),pmx(n0:nx),rhomx(n0:nx) common /FEQSAV/ phisav(n0:nx,0:mxtrm),delsav(n0:nx),tausav(n0:nx), & taup(n0:nx,mxtrm),delp(n0:nx,mxtrm), & delli(n0:nx,mxtrm),expdel(n0:nx,mxtrm) c c (re)initialize contents of /FEQSAV/ when a new fluid is read in do 120 i=n0,nx phisav(i,0)=0.0d0 delsav(i)=0.0d0 tausav(i)=0.0d0 do 100 j=1,mxtrm phisav(i,j)=0.0d0 taup(i,j)=0.0d0 delp(i,j)=0.0d0 delli(i,j)=0.0d0 expdel(i,j)=0.0d0 100 continue 120 continue c if (nread.le.0) then c get coefficients from block data c identify specified fluid with entries in database via match of CAS no do 200 i=1,mxfeq if (hcasno.eq.hcasa(i)) then hmodcp(icomp)=hcpfeq(i) ! pointer to ideal gas model ntermf(icomp)=ntrmfa(i) ! number of terms do 198 j=1,ntermf(icomp) li(icomp,j)=lia(i,j) ! transfer coefficients a(icomp,j)=aa(i,j) ! into working arrays di(icomp,j)=dia(i,j) ti(icomp,j)=tia(i,j) 198 continue rho0(icomp)=rho0a(i) ! reducing parameters t0(icomp)=t0a(i) pc(icomp)=pca(i) ! critical parameters rhoc(icomp)=rhoca(i) tc(icomp)=tca(i) wmf(icomp)=wmfa(i) ! molecular weight Rfeq(icomp)=Rfeqa(i) ! gas constant if (nc.eq.1 .and. icomp.eq.1) then R=Rfeq(icomp) end if pmin(icomp)=pmina(i) ! limits rhomx(icomp)=rhomxa(i) rhotp(icomp)=rhomxa(i) tmin(icomp)=tmina(i) tmax(icomp)=tmaxa(i) pmax(icomp)=pmaxa(i) ncoeff(icomp)=4 !# coeff per term ncrt(icomp)=0 !no critical-region terms nspare(icomp)=0 c fill arrays in /CCON/ wm(icomp)=wmfa(i) ttp(icomp)=ttpf(i) tnbp(icomp)=tnbpf(i) tcrit(icomp)=tc(icomp) pcrit(icomp)=pc(icomp) Dcrit(icomp)=rhoc(icomp) Zcrit(icomp)=pc(icomp)/(Rfeq(icomp)*tc(icomp)*rhoc(icomp)) accen(icomp)=accenf(i) dipole(icomp)=dipm(i) ierr=0 herr=hnull goto 990 end if 200 continue ierr=1 herr=' ERROR--fluid input to SETFEQ not found'//hnull else c read data from file c write (*,*) ' SETFEQ--read component',icomp,' from unit',nread read (nread,*) tmin(icomp) !lower temperature limit c write (*,*) ' SETFEQ--first input: tmin: ',tmin(icomp) read (nread,*) tmax(icomp) !upper temperature limit read (nread,*) pmax(icomp) !upper pressure limit read (nread,*) rhomx(icomp) !upper density limit read (nread,2003) hmodcp(icomp) !pointer to Cp0 model read (nread,*) wm(icomp) !molecular weight wmf(icomp)=wm(icomp) read (nread,*) ttp(icomp) !triple point temperature read (nread,*) pmin(icomp) !pressure at triple point read (nread,*) rhotp(icomp) !density at triple point read (nread,*) tnbp(icomp) !normal boiling point temperature read (nread,*) accen(icomp) !accentric factor read (nread,*) tc(icomp),pc(icomp),rhoc(icomp) !critical par tcrit(icomp)=tc(icomp) pcrit(icomp)=pc(icomp) Dcrit(icomp)=rhoc(icomp) read (nread,*) t0(icomp),rho0(icomp) !reducing parameters read (nread,*) Rfeq(icomp) !gas constant used in fit if (nc.eq.1 .and. icomp.eq.1) then R=Rfeq(icomp) c write (*,*) ' SETFEQ--R set to ',R end if Zcrit(icomp)=pc(icomp)/(Rfeq(icomp)*tc(icomp)*rhoc(icomp)) read (nread,*) nterm,ncoeff(icomp),ncrt(icomp),ncfcrt(icomp), & nspare(icomp),ncfsp(icomp) ntermf(icomp)=nterm c write (*,*) ' SETFEQ--about to read ',nterm,' coefficients' c the gi term is a multiplier for the (rho or del) in only the exponential c terms; it is needed for e.g. the Bender EOS; set to 1.0 if not present do 240 j=1,nterm if (ncoeff(icomp).eq.5) then read (nread,*) a(icomp,j),ti(icomp,j),di(icomp,j),li(icomp,j), & gi(icomp,j) else read (nread,*) a(icomp,j),ti(icomp,j),di(icomp,j),li(icomp,j) gi(icomp,j) = 1.0d0 endif 240 continue if (ncrt(icomp).gt.0) then do 250 j=1,ncrt(icomp) read (nread,*) a(icomp,j+nterm),ti(icomp,j+nterm), & di(icomp,j+nterm),li(icomp,j+nterm), & alpha(icomp,j),beta(icomp,j),gamma(icomp,j),delta(icomp,j), & eta(icomp,j),eid(icomp,j),eit(icomp,j) 250 continue end if c write (*,*) ' SETFEQ--final coefficient: ',a(icomp,nterm) end if c !temporary--if steam, reset gas constant to fluid-specific value c (this is a patch to cover the odd value of R used with steam) if (nc.eq.1 .and. hcas(icomp).eq.'7732-18-5 ') then R=Rfeq(icomp) c write (*,*) ' SETFEQ--gas const has been reset for steam: ',R end if c 990 continue c c scan and sort exponents to minimize computations; c the array tpower--stores the discrete powers of temperature c itp--is a pointer to the appropriate element of tpower c dpower, dtp--ditto for density c lpower, ltp--ditto for the powers of del in the exp terms c ntp, ndp, and nlp store the number of entries in the above arrays c tpower(icomp,1)=ti(icomp,1) dpower(icomp,1)=di(icomp,1) lpower(icomp,1)=li(icomp,1) gi2(icomp,1)=gi(icomp,1) ntp(icomp)=1 ndp(icomp)=1 nlp(icomp)=1 c do 480 i=1,ntermf(icomp)+ncrt(icomp) c compare the power of temperature for each term with elements c already in the tpower array do 420 j=1,ntp(icomp) if (abs(ti(icomp,i)-tpower(icomp,j)).lt.1.0d-10) then itp(icomp,i)=j goto 422 end if 420 continue c power of temperature for term i does not match any previous term ntp(icomp)=ntp(icomp)+1 tpower(icomp,ntp(icomp))=ti(icomp,i) itp(icomp,i)=ntp(icomp) 422 continue c compare the power of density for each term with elements c already in the dpower array do 440 j=1,ndp(icomp) if (abs(di(icomp,i)-dpower(icomp,j)).lt.1.0d-10) then idp(icomp,i)=j goto 442 end if 440 continue c power of density for term i does not match any previous term ndp(icomp)=ndp(icomp)+1 dpower(icomp,ndp(icomp))=di(icomp,i) idp(icomp,i)=ndp(icomp) 442 continue c compare the power of density in exponential terms with elements c already in the lpower array do 460 j=1,nlp(icomp) if (li(icomp,i).eq.lpower(icomp,j)) then ilp(icomp,i)=j goto 462 end if 460 continue c exponential of density for term i does not match any previous term nlp(icomp)=nlp(icomp)+1 lpower(icomp,nlp(icomp))=li(icomp,i) ilp(icomp,i)=nlp(icomp) gi2(icomp,nlp(icomp))=gi(icomp,i) 462 continue 480 continue c c write out all coefficients for debugging c c i=icomp c write (*,*) ' SETFEQ--coefficients for comp',i,', CAS # ',hcas(i) c write (*,*) ' Cp0 model: ',hmodcp(i) c write (*,*) ' Nterms : ',ntermf(i) c write (*,2020) ((a(i,j),ti(i,j),di(i,j),li(i,j)),j=1,ntermf(i)) c2020 format (1x,d24.14,2f9.3,i4) c c write (*,*) ' SETFEQ--sorted coefs for comp',i,', CAS # ',hcas(i) c write (*,2020) ((a(i,j),tpower(i,itp(i,j)),dpower(i,idp(i,j)), c & lpower(i,ilp(i,j))),j=1,ntermf(i)) c c copy limits into /EOSLIM/ arrays tmn(icomp)=tmin(icomp) tmx(icomp)=tmax(icomp) pmx(icomp)=pmax(icomp) c rhomx(icomp)=rhomax(icomp) c RETURN 2003 format (a3) end !subroutine SETFEQ c c ====================================================================== c block data BDFEQ c c data for Helmholtz-explicit equations of state c c explanation of parameters c mxfeq: maximum number of FEQ fits, used to dimension arrays c mxtrm: max number of terms per fit, used to dimension arrays c c explanation of commons and constituent arrays c /CASFEQ/ Chem Abstract number; used as unambiguous identifier c hcas(i): CAS number for fluid corresponding to equation "i" c c /CPMFEQ/ c hmodcp(i): pointer to ideal gas model to use with fluid "i" c c /NTFEQ/ c ntermf(i): number of terms in fit corresponding to fluid "i" c c /LFFEQ/ parameters to FEQ fits for each of mxfeq fluids c li(i,1..N): power of del in the exponential multiplier c (if li(j)=0, then multiplier = 1) c c /CFFEQ/ parameters to FEQ fits for each of mxfeq fluids c a(i,1..N): N coefficients to "fundamental" equation of state c di(i,1..N): exponents for the N density terms c ti(i,1..N): exponents for the N temperature terms c rho0(i): reducing parameter for density (mol/L) c t0(i): reducing parameter for temperature (K) c pc(i): critical pressure (kPa) c rhoc(i): critical density (L/mol) c tc(i): critical temperature (K) c wmf(i): molecular weight (g/mol) c Rfeq(i): gas constant used in fit (J/(mol-K)) c pmin(i): pressure at tmin(i) (kPa) c rhotp(i): density at tmin(i), e.g. triple point (L/mol) c tmin(i): low temperature limit of fit--often triple point c pmin(i),rhotp(i) are used for initial guesses, etc. c and are often approximate values only c tmax(i): upper temperature limit of fit (K) c pmax(i): upper pressure limit of fit (kPa) c c /MSCFEQ/ miscellaneous fluid constants c ttpf(i): triple point temperature (K) c tnbpf(i): normal boiling point temperature (K) c accenf(i): accentric factor for fluid represented by eqn "i" c dipm(i): dipole moment [debye] (at Tnbp if t-dependent) c c /FEQSAV/ used to save information between calls to PHIFEQ c phisav(i,j):individual terms in summation c delsav(i): reduced density on last call to PHIFEQ c tausav(i): reduced temperature on last call to PHIFEQ c taup(i,j): reduced temperature raised to power c delp(i,j): reduced density raised to power c expdel(i,j):exp(reduced density raised to power) c c where "i" is the equation number c N.B.--the "i" are, in general, not the same as fluid code numbers c c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 07-26-95 MM, original version c 10-03-95 MM, /CPMFEQ/ changed to character variable c 11-01-95 MM, increase parameter mxtrm to 52 (to accomodate steam) c 11-29-95 MM, variable lower limit on coefficient/constant arrays c to accomodate ECS reference fluid c 01-16-96 MM, implicit integer (i-n); (include L) c 02-27-96 MM, parameter n0=-ncmax to accomodate ECS-thermo model c 01-03-97 MM, arrays added to /FEQSAV/ (assoc with sorting of powers) 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 (mxfeq=2) !max number of FEQ EOS in block data parameter (mxtrm=56) parameter (nsave=(nx-n0+1)*(mxtrm+1)) parameter (nsave1=(nx-n0+1)*mxtrm) parameter (nxsav=nx-n0+1) character*3 hcpfeq character*12 hcas common /CASFEQ/ hcas(mxfeq) common /CPMFEQ/ hcpfeq(mxfeq) common /NTFEQ/ ntermf(mxfeq) common /LFFEQ/ li(mxfeq,mxtrm) common /CFFEQ/ a(mxfeq,mxtrm),di(mxfeq,mxtrm),ti(mxfeq,mxtrm), & rho0(mxfeq),t0(mxfeq), & pc(mxfeq),rhoc(mxfeq),tc(mxfeq), & wmf(mxfeq),Rfeq(mxfeq), & pmin(mxfeq),rhotp(mxfeq),tmin(mxfeq), & tmax(mxfeq),pmax(mxfeq) common /MSCFEQ/ ttpf(mxfeq),tnbpf(mxfeq),accenf(mxfeq),dipm(mxfeq) common /FEQSAV/ phisav(n0:nx,0:mxtrm),delsav(n0:nx),tausav(n0:nx), & taup(n0:nx,mxtrm),delp(n0:nx,mxtrm), & delli(n0:nx,mxtrm),expdel(n0:nx,mxtrm) c data phisav /nsave*0.0d0/ data delsav /nxsav*0.0d0/ data tausav /nxsav*0.0d0/ data taup /nsave1*0.0d0/ data delp /nsave1*0.0d0/ data delli /nsave1*0.0d0/ data expdel /nsave1*0.0d0/ c c c R134a 1,1,1,2-tetrafluoroethane data hcas(1) /'811-97-2'/ c use polynomial Cp0 model (at least for now) data hcpfeq(1) /'CPP'/ data ntermf(1) /21/ data (a(1,j),j=1,21)/ & 0.5586817d-01, 0.4982230d+00, 0.2458698d-01, 0.8570145d-03, & 0.4788584d-03, -0.1800808d+01, 0.2671641d+00, -0.4781652d-01, & 0.1423987d-01, 0.3324062d+00, -0.7485907d-02, 0.1017263d-03, & -0.5184567d+00, -0.8692288d-01, 0.2057144d+00, -0.5000457d-02, & 0.4603262d-03, -0.3497836d-02, 0.6995038d-02, -0.1452184d-01, & -0.1285458d-03/ c exponents for density terms data (di(1,j),j=1,21)/ & 2.00d+00, 1.00d+00, 3.00d+00, 6.00d+00, & 6.00d+00, 1.00d+00, 1.00d+00, 2.00d+00, & 5.00d+00, 2.00d+00, 2.00d+00, 4.00d+00, & 1.00d+00, 4.00d+00, 1.00d+00, 2.00d+00, & 4.00d+00, 1.00d+00, 5.00d+00, 3.00d+00, & 10.00d+00/ c exponents for temperature terms data (ti(1,j),j=1,21)/ & -0.50d+00, 0.00d+00, 0.00d+00, 0.00d+00, & 1.50d+00, 1.50d+00, 2.00d+00, 2.00d+00, & 1.00d+00, 3.00d+00, 5.00d+00, 1.00d+00, & 5.00d+00, 5.00d+00, 6.00d+00, 10.00d+00, & 10.00d+00, 10.00d+00, 18.00d+00, 22.00d+00, & 50.00d+00/ c power of (-del) in the exponential multiplier data (li(1,j),j=1,21)/ & 0, 0, 0, 0, & 0, 0, 0, 0, & 1, 1, 1, 2, & 2, 2, 2, 2, & 2, 3, 3, 3, & 4/ c data rho0(1),t0(1) & /4.978830171d0,374.18d0/ ! reducing parameters data pc(1),rhoc(1),tc(1) & /4059.28d0,5.017053d0,374.21d0/ ! critical parameters data wmf(1),Rfeq(1) & /102.032d0,8.314471d0/ ! mol weight and gas constant data pmin(1),rhotp(1),tmin(1) & /0.391d0,15.594d0,169.85d0/ ! lower limits (= triple point) data tmax(1),pmax(1) & /453.15d0,70000.0d0/ ! upper limits c data ttpf(1) /169.85d0/ data tnbpf(1) /247.07d0/ data accenf(1) /0.32684d0/ data dipm(1) /2.058d0/ !dipole moment [Debye]; Meyer, (1991) c end !block data BDFEQ c c c 1 2 3 4 5 6 7 c23456789012345678901234567890123456789012345678901234567890123456789012 c c ====================================================================== c end file core_FEQ.f c ======================================================================