c begin file core_PH0.f c c This file contains the functions implementing the ideal-gas part of c the reduced Helmholtz free energy form of the pure fluid equation of c state (the so-called "fundamental equation"). c c The Helmholtz energy consists of ideal and residual (real-gas) terms; c this routine calculates only the ideal part, and only for pure components. c c contained here are: c subroutine SETPH0 (nread,icomp,hcasno,ierr,herr) c function PH0PH0 (icomp,itau,idel,t,rho) c block data SAVPH0 c c these routines use the following common blocks from other files c common /CREF/ tref(n0:nx),rhoref(n0:nx),href(n0:nx),sref(n0:nx) c common /HCHAR/ htab,hnull c common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx) c c various arrays are dimensioned with parameter statements c parameter (ncmax=20) !max number of components in mixture c parameter (nrefmx=10) !max number of fluids for transport ECS c parameter (n0=-ncmax-nrefmx,nx=ncmax) c parameter (nph0mx=10) !max number of terms in Cp0 polynomial c c ====================================================================== c ====================================================================== c subroutine SETPH0 (nread,icomp,hcasno,ierr,herr) c c set up working arrays for ideal-gas part of the Helmholtz energy c implements an expression of the form: c phi0 = Sum[ai*log(tau**ti)] + Sum[aj*tau**tj] c + Sum[ak*log(1-EXP(bk*tau))] c c inputs: c nread--file to read data from c <= 0 get data from block data (not currently implemented) 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 zero and negative numbers designate ECS reference fluids c hcasno--CAS number of component icomp c (not req'd if reading from file--included to maintain c parallel structure with other routines) 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 coefficients, etc. returned via arrays in common /xxxPH0/ c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 05-13-96 MM, original version, based (loosely) on SETCPP c 06-14-96 MM, fix bug in reading/storing coefficients c 11-13-97 MM, (re)initialize contents of /PH0SAV/ when a new fluid is read in c 08-31-98 MEV, reverse order of indices in tsav(i,j)=0 and rhosav(i,j)=0 c 08-13-98 MM, delete obsolete (unused) format statement c 07-25-06 EWL, add cosh and sinh functions c implicit double precision (a-h,o-z) implicit integer (i-n) parameter (mxph0=2) !max number of fluids in block data parameter (ncmax=20) !max number of components in mixture parameter (nrefmx=10) !max number of fluids for transport ECS parameter (n0=-ncmax-nrefmx,nx=ncmax) parameter (nph0mx=10) !max number of terms in phi0 function character*1 htab,hnull character*12 hcasno,hcas character*255 herr common /HCHAR/ htab,hnull c commons associated with the nc components of current interest c ("working" commons and arrays) common /CASPH0/ hcas(mxph0) common /WNTPH0/ nlog(n0:nx),ntau(n0:nx),nexp(n0:nx),ncosh(n0:nx), & nsinh(n0:nx),nsp1(n0:nx),nsp2(n0:nx),nsp3(n0:nx) common /WLMPH0/ tmin(n0:nx),tmax(n0:nx),pmax(n0:nx),rhomax(n0:nx) common /WCFPH0/ ai(n0:nx,nph0mx),ti(n0:nx,nph0mx) common /PH0SAV/ ph0sav(n0:nx),ph1sav(n0:nx),ph2sav(n0:nx), & tsav(0:2,n0:nx),rhosav(0:2,n0:nx) c c (re)initialize contents of /PH0SAV/ when a new fluid is read in common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx) common /VERS/ verfl(n0:nx),vermx !fluid & mix file version nos. do i=n0,nx ph0sav(i)=0.0d0 ph1sav(i)=0.0d0 ph2sav(i)=0.0d0 do j=0,2 tsav(j,i)=0.0d0 rhosav(j,i)=0.0d0 enddo enddo 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 k=1,mxph0 if (hcasno.eq.hcas(k)) then c write (*,*)' SETPH0--ERROR, block data read not implemented' ierr=1 herr='[SETPH0 error] Block data option not implemented'//hnull RETURN end if enddo ierr=1 herr='[SETPH0 error] Input fluid (block data) not found'//hnull RETURN else c read data from file c write (*,*) ' SETPH0--read component',icomp,' from unit',nread 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 c read number of terms for each of the various types if (verfl(icomp).ge.8.0d0) then read (nread,*) nlog(icomp),ntau(icomp),nexp(icomp), & ncosh(icomp),nsinh(icomp), & nsp1(icomp),nsp2(icomp),nsp3(icomp) !spares else read (nread,*) nlog(icomp),ntau(icomp),nexp(icomp) ncosh(icomp)=0 !these terms not used in files prior to v8.0 nsinh(icomp)=0 nsp1(icomp)=0 nsp2(icomp)=0 nsp3(icomp)=0 endif jterm=nlog(icomp)+ntau(icomp)+nexp(icomp)+ncosh(icomp) & +nsinh(icomp)+nsp1(icomp)+nsp2(icomp)+nsp3(icomp) if (jterm.ge.1) then c read coefficients for terms of the form [ai*log(tau**ti)], c [ai*tau**ti], and [ai*log(1-EXP(bi*tau))] do i=1,jterm read (nread,*) ai(icomp,i),ti(icomp,i) enddo end if ierr=0 herr=' ' end if c RETURN end !subroutine SETPH0 c c ====================================================================== c function PH0PH0 (icomp,itau,idel,t,rho) c c compute the ideal gas part of the reduced Helmholtz energy or a c derivative as functions of temperature and pressure; for c use with a Helmholtz-explicit equation of state 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 (the density derivatives are not used in the calculation c of any property, and are not implemented) c when itau = 0 and idel = 0, compute A0/RT c when itau = 1 and idel = 0, 1st temperature derivative c when itau = 2 and idel = 0, 2nd temperature derivative c t--temperature (K) c rho--density (mol/L) c output (as function value): c ph0ph0--ideal-gas part of the Helmholtz energy in reduced form (A/RT) c derivatives (as specified by itau and idel) are multiplied c by the corresponding power of tau; i.e. when itau = 1, the c quantity returned is tau*d(ph0ph0)/d(tau) and when itau = 2, c the quantity returned is tau*tau*d2(ph0ph0)/d(tau)**2 c c Note: While the real-gas part of the Helmholtz energy is calculated c in terms of dimensionless temperature and density, the ideal- c gas part is calculated in terms of absolute temperature and c density. (This distinction is necessary for mixtures.) c c The Helmholtz energy consists of ideal-gas and residual c (real-gas) terms; this routine calculates only the ideal part. c c This function computes pure component properties only. c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 05-13-96 MM, original version, based (loosely) on PH0CPP c 06-14-96 MM, add rhosav to /PH0SAV/ c 07-08-96 MM, change derivative outputs: tau*d(phi)/d(tau), etc c 08-20-97 MM, call ERRMSG if itau out of range; drop idel=idel c 07-11-00 EWL, remove krypton pieces c 07-25-06 EWL, add cosh and sinh functions c 08-24-06 EWL, add check for PHG in hmodcp (add EOSMOD common block). c If it is used, multiple by R*/R as described in GERG-2004 eos (of Kunz and Wagner) c implicit double precision (a-h,o-z) implicit integer (i-k,m,n) implicit logical (l) parameter (ncmax=20) !max number of components in mixture parameter (nrefmx=10) !max number of fluids for transport ECS parameter (n0=-ncmax-nrefmx,nx=ncmax) parameter (nph0mx=10) !max number of terms in phi0 function character*1 htab,hnull character*255 herr character*3 hpheq,heos,hmxeos,hmodcp common /HCHAR/ htab,hnull common /CREF/ tref(n0:nx),rhoref(n0:nx),href(n0:nx),sref(n0:nx) common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx) c commons associated with the nc components of current interest c ("working" commons and arrays) common /WNTPH0/ nlog(n0:nx),ntau(n0:nx),nexp(n0:nx),ncosh(n0:nx), & nsinh(n0:nx),nsp1(n0:nx),nsp2(n0:nx),nsp3(n0:nx) common /WLMPH0/ tmin(n0:nx),tmax(n0:nx),pmax(n0:nx),rhomax(n0:nx) common /WCFPH0/ ai(n0:nx,nph0mx),ti(n0:nx,nph0mx) c saved values from previous calls to this routine common /PH0SAV/ ph0sav(n0:nx),ph1sav(n0:nx),ph2sav(n0:nx), & tsav(0:2,n0:nx),rhosav(0:2,n0:nx) common /EOSFLG/ kryptn(n0:nx),ispr1(n0:nx),ispr2(n0:nx), & ispr3(n0:nx),ispr4(n0:nx),ispr5(n0:nx), & ispr6(n0:nx) common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx) c c compute reduced Helmholtz; first check if t same as previous call c PH0PH0=0.0d0 !initialize in case of error if (t.le.0) RETURN tau=tz(icomp)/t del=rho/rhoz(icomp) if (itau*idel.ne.0) then PH0PH0=0.0d0 elseif (itau.eq.0 .and. idel.eq.0) then if (abs(t-tsav(0,icomp)).lt.1.0d-8 .and. & abs(rho-rhosav(0,icomp)).lt.1.0d-10) then PH0PH0=ph0sav(icomp) c write (*,*) ' PH0PH0--using saved phi for itau = 0, t =',t else c write (*,1030) icomp,t,rho,t0,D0,tau,del c1030 format (1x,'PH0PH0--icomp,t,rho,t0,D0,tau,del: ',i4,6e14.6) iterm=0 phisum=LOG(del) c & +href(icomp)/R/t-sref(icomp)/R !see U. Idaho class notes if (nlog(icomp).ge.1) then c sum terms of the form [ai*log(tau**ti)] do i=1,nlog(icomp) iterm=iterm+1 phisum=phisum+ai(icomp,iterm)*LOG(tau**ti(icomp,iterm)) c write (*,1040) iterm,ai(icomp,iterm),ti(icomp,iterm),phisum c1040 format (' PH0PH0--log term iterm,ai,ti,phisum: ',i3,3f12.6) enddo end if if (ntau(icomp).ge.1) then c sum terms of the form [ai*tau**ti] do j=1,ntau(icomp) iterm=iterm+1 phisum=phisum+ai(icomp,iterm)*tau**ti(icomp,iterm) c write (*,1060) iterm,ai(icomp,iterm),ti(icomp,iterm),phisum c1060 format (' PH0PH0--tau term iterm,ai,ti,phisum: ',i3,3f12.6) enddo end if if (nexp(icomp).ge.1) then c sum terms of the form [ai*log(1-EXP(bi*tau))] c (the bi coefficients are stored in the ti array) do k=1,nexp(icomp) iterm=iterm+1 phisum=phisum+ai(icomp,iterm) & *LOG(1.0d0-EXP(ti(icomp,iterm)*tau)) c write (*,1080) iterm,ai(icomp,iterm),ti(icomp,iterm),phisum c1080 format (' PH0PH0--exp term iterm,ai,ti,phisum: ',i3,3f12.6) enddo end if if (ncosh(icomp).ge.1) then c sum terms of the form [ai*log(cosh(bi*tau))] do k=1,ncosh(icomp) iterm=iterm+1 ttau=ti(icomp,iterm)*tau if (ttau.lt.700.d0) then phisum=phisum+ai(icomp,iterm)*LOG(COSH(ttau)) else phisum=phisum+ai(icomp,iterm)*700.d0 endif enddo end if if (nsinh(icomp).ge.1) then c sum terms of the form [ai*log(sinh(bi*tau))] do k=1,nsinh(icomp) iterm=iterm+1 ttau=ti(icomp,iterm)*tau if (ttau.lt.700.d0) then phisum=phisum+ai(icomp,iterm)*LOG(SINH(ttau)) else phisum=phisum+ai(icomp,iterm)*700.d0 endif enddo end if PH0PH0=phisum c save information for possible use on next call to function tsav(0,icomp)=t rhosav(0,icomp)=rho ph0sav(icomp)=PH0PH0 end if c c compute derivative w.r.t. tau (dimensionless temperature) c else if (itau.eq.1) then if (abs(t-tsav(1,icomp)).lt.1.0d-8 .and. & abs(rho-rhosav(1,icomp)).lt.1.0d-10) then PH0PH0=ph1sav(icomp) else iterm=0 phisum=0.0d0 c phisum=href(icomp)/(R*tz) !see U. Idaho class notes if (nlog(icomp).ge.1) then c sum terms of the form [ai*log(tau**ti)] do i=1,nlog(icomp) iterm=iterm+1 phisum=phisum+ai(icomp,iterm)*ti(icomp,iterm)/tau enddo end if if (ntau(icomp).ge.1) then c sum terms of the form [ai*tau**ti] do j=1,ntau(icomp) iterm=iterm+1 phisum=phisum+ai(icomp,iterm)*ti(icomp,iterm) & *tau**(ti(icomp,iterm)-1.0d0) enddo end if if (nexp(icomp).ge.1) then c sum terms of the form [ai*log(1-EXP(bi*tau))] c (the bi coefficients are stored in the ti array) do k=1,nexp(icomp) iterm=iterm+1 exptau=EXP(ti(icomp,iterm)*tau) phisum=phisum-ai(icomp,iterm)*ti(icomp,iterm) & *exptau/(1.0d0-exptau) enddo end if if (ncosh(icomp).ge.1) then c sum terms of the form [ai*log(cosh(bi*tau))] do k=1,ncosh(icomp) iterm=iterm+1 phisum=phisum+ai(icomp,iterm)*ti(icomp,iterm) & *(TANH(ti(icomp,iterm)*tau)) enddo end if if (nsinh(icomp).ge.1) then c sum terms of the form [ai*log(sinh(bi*tau))] do k=1,nsinh(icomp) iterm=iterm+1 phisum=phisum+ai(icomp,iterm)*ti(icomp,iterm) & /(TANH(ti(icomp,iterm)*tau)) enddo end if c save information for possible use on next call to function PH0PH0=phisum*tau !return tau*d(ph0ph0)/d(tau) tsav(1,icomp)=t rhosav(1,icomp)=rho ph1sav(icomp)=PH0PH0 end if c c compute 2nd derivative w.r.t. tau (dimensionless temperature) c else if (itau.eq.2) then if (abs(t-tsav(2,icomp)).lt.1.0d-8 .and. & abs(rho-rhosav(2,icomp)).lt.1.0d-10) then PH0PH0=ph2sav(icomp) else iterm=0 phisum=0.0d0 if (nlog(icomp).ge.1) then c sum terms of the form [ai*log(tau**ti)] do i=1,nlog(icomp) iterm=iterm+1 phisum=phisum-ai(icomp,iterm)*ti(icomp,iterm)/tau**2 enddo end if if (ntau(icomp).ge.1) then c sum terms of the form [ai*tau**ti] do j=1,ntau(icomp) iterm=iterm+1 phisum=phisum+ai(icomp,iterm)*ti(icomp,iterm) & *(ti(icomp,iterm)-1.0d0)*tau**(ti(icomp,iterm)-2.0d0) enddo end if if (nexp(icomp).ge.1) then c sum terms of the form [ai*log(1-EXP(bi*tau))] c (the bi coefficients are stored in the ti array) do k=1,nexp(icomp) iterm=iterm+1 exptau=EXP(ti(icomp,iterm)*tau) phisum=phisum-ai(icomp,iterm)*ti(icomp,iterm)**2 & *exptau/(1.0d0-exptau)**2 enddo end if if (ncosh(icomp).ge.1) then c sum terms of the form [ai*log(cosh(bi*tau))] do k=1,ncosh(icomp) iterm=iterm+1 phisum=phisum+ai(icomp,iterm)*ti(icomp,iterm)**2 & /(COSH(ti(icomp,iterm)*tau))**2 enddo end if if (nsinh(icomp).ge.1) then c sum terms of the form [ai*log(sinh(bi*tau))] do k=1,nsinh(icomp) iterm=iterm+1 phisum=phisum-ai(icomp,iterm)*ti(icomp,iterm)**2 & /(SINH(ti(icomp,iterm)*tau))**2 enddo end if c save information for possible use on next call to function PH0PH0=phisum*tau*tau tsav(2,icomp)=t rhosav(2,icomp)=rho c return tau**2*d2(ph0ph0)/d(tau**2) ph2sav(icomp)=PH0PH0 end if else if (itau.eq.3) then iterm=0 phisum=0.0d0 if (nlog(icomp).ge.1) then do i=1,nlog(icomp) iterm=iterm+1 phisum=phisum+2.d0*ai(icomp,iterm)*ti(icomp,iterm)/tau**3 enddo end if if (ntau(icomp).ge.1) then do j=1,ntau(icomp) iterm=iterm+1 phisum=phisum+ai(icomp,iterm)*ti(icomp,iterm) & *(ti(icomp,iterm)-1.0d0) & *(ti(icomp,iterm)-2.0d0)*tau**(ti(icomp,iterm)-3.0d0) enddo end if if (nexp(icomp).ge.1) then do k=1,nexp(icomp) iterm=iterm+1 exptau=EXP(ti(icomp,iterm)*tau) phisum=phisum-ai(icomp,iterm)*ti(icomp,iterm)**3* & (exptau/(1.0d0-exptau)**2+ & 2.d0*exptau**2/(1.0d0-exptau)**3) enddo end if if (ncosh(icomp).ge.1) then do k=1,ncosh(icomp) iterm=iterm+1 phisum=phisum-2.d0*ai(icomp,iterm)*ti(icomp,iterm)**3 & /(COSH(ti(icomp,iterm)*tau))**3 & *(SINH(ti(icomp,iterm)*tau)) enddo end if if (nsinh(icomp).ge.1) then do k=1,nsinh(icomp) iterm=iterm+1 phisum=phisum+2.d0*ai(icomp,iterm)*ti(icomp,iterm)**3 & /(SINH(ti(icomp,iterm)*tau))**3 & *(COSH(ti(icomp,iterm)*tau)) enddo end if PH0PH0=phisum*tau**3 elseif (idel.eq.1) then PH0PH0=1.0d0 elseif (idel.eq.2) then PH0PH0=-1.0d0 elseif (idel.eq.3) then PH0PH0=2.0d0 else c c invalid value of itau c ierr=99 write (herr,1099) itau,idel,hnull 1099 format ('[PH0PH0 warning] invalid input; itau =',i4,'; idel =', & i4,a1) call ERRMSG (ierr,herr) PH0PH0=0.0d0 end if if (hmodcp(icomp).eq.'PHG') then PH0PH0=8.31451D0/8.314472D0*PH0PH0 endif c c write (*,*) ' PH0PH0: output phi: ',ph0ph0 c RETURN end !function PH0PH0 c c 1 2 3 4 5 6 7 c23456789012345678901234567890123456789012345678901234567890123456789012 c c ====================================================================== c end file core_PH0.f c ======================================================================