c begin file core_CPP.f c c This file contains the functions implementing a polynomial fit c for the ideal gas heat capacity. c c routines for the polynomial fit of ideal gas heat capacity c c contained here are: c subroutine SETCPP (nread,icomp,hcasno,ierr,herr) c function CP0CPP (icomp,t) c function CPICPP (icomp,t) c function CPTCPP (icomp,t) c function PH0CPP (icomp,itau,idel,t,rho) c block data SAVCPP c block data BDCPP 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 (mxcpp=2) !max number of fluids in block data 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 (ncppmx=20) !max number of Cp0 terms c c ====================================================================== c ====================================================================== c subroutine SETCPP (nread,icomp,hcasno,ierr,herr) c c set up working arrays for polynomial form of ideal heat gas Cp c uses cubic polynomial correlation of Cp0 (dimensional form): c Cp0 = cpc(i,0) + cpc(i,1)*T + cpc(i,2)*T*T + cpc(i,3)*T**3 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 zero and negative numbers designate ECS reference fluids 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 coefficients, etc. returned via arrays in commons /xxxCPP/ 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-04-95 MM, adapt to file input, add nread to argument list c 11-10-95 MM, convert from cubic to general polynomial c 11-29-95 MM, variable lower limit on coefficient/constant arrays c to accommodate ECS reference fluid c 02-27-96 MM, parameter n0=-ncmax to accommodate ECS-thermo model c 06-17-96 MM, read in exponential terms (Einstein function) c 11-13-97 MM, (re)initialize contents of /CPPSAV/ when a new fluid is read in c 12-05-97 MM, change ncppmx from 15 to 20 (accommodate H2 in NIST14) c 04-02-98 EWL, read in number of cosh and sinh terms c 06-29-98 MEV, initialize ncosh and nsinh for version previous to 6.1 c 08-13-98 MM, delete obsolete (unused) format statement c implicit double precision (a-h,o-z) implicit integer (i-k,m,n) parameter (mxcpp=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 (ncppmx=20) !max number of Cp0 terms character*1 htab,hnull character*12 hcasno,hcas character*255 herr common /HCHAR/ htab,hnull c commons associated with the mxcpp fluids with FEQ equations stored c in block data BDCPP common /NTCPP/ ntrmca(mxcpp),ntrmce(mxcpp) common /RDCPP/ treda(mxcpp),Creda(mxcpp) common /LMCPP/ tmna(mxcpp),tmxa(mxcpp),pmxa(mxcpp),rhomxa(mxcpp) common /CPCPP/ cpca(mxcpp,ncppmx),xka(mxcpp,ncppmx) common /CASCPP/ hcas(mxcpp) c commons associated with the nc components of current interest c ("working" commons and arrays) common /WNTCPP/ ntermc(n0:nx),nterme(n0:nx),ncosh(n0:nx), & nsinh(n0:nx),nsp1(n0:nx),nsp2(n0:nx),nsp3(n0:nx) common /WRDCPP/ tred(n0:nx),Cred(n0:nx) common /WLMCPP/ tmin(n0:nx),tmax(n0:nx),pmax(n0:nx),rhomax(n0:nx) common /WCPCPP/ cpc(n0:nx,ncppmx),xk(n0:nx,ncppmx), & cph(n0:nx,ncppmx),xth(n0:nx,ncppmx), & xh(n0:nx,ncppmx) common /CPPSAV/ cp0sav(n0:nx),cpisav(n0:nx),cptsav(n0:nx), & tsav(n0:nx) common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx) common /VERS/ verfl(n0:nx),vermx !fluid & mix file version nos. c c (re)initialize contents of /CPPSAV/ when a new fluid is read in do i=n0,nx cp0sav(i)=0.0d0 cpisav(i)=0.0d0 cptsav(i)=0.0d0 tsav(i)=0.0d0 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,mxcpp if (hcasno.eq.hcas(k)) then c write (*,*) ' SETCPP--coeff from block data for CAS # ',hcasno tmin(icomp)=tmna(k) tmax(icomp)=tmxa(k) pmax(icomp)=pmxa(k) rhomax(icomp)=rhomxa(k) tred(icomp)=treda(k) Cred(icomp)=Creda(k) ntermc(icomp)=ntrmca(k) nterme(icomp)=ntrmce(k) ncosh(icomp)=0 nsinh(icomp)=0 nsp1(icomp)=0 !spares for future use nsp2(icomp)=0 nsp3(icomp)=0 do j=1,ntrmca(k) cpc(icomp,j)=cpca(k,j) xk(icomp,j)=xka(k,j) enddo c write (*,*) ' SETCPP--final coeff: ',cpc(icomp,ntermc(icomp)) ierr=0 herr=' ' RETURN end if enddo ierr=1 herr='[SETCPP error] Fluid input not found'//hnull RETURN else c read data from file c write (*,*) ' SETCPP--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 read (nread,*) tred(icomp),Cred(icomp) !reducing parameters c read number of polynomial and exponential terms if (verfl(icomp).ge.6.099d0) then c hyperbolic terms (and spares for future) added for file version 6.100 read (nread,*) ntermc(icomp),nterme(icomp), & ncosh(icomp),nsinh(icomp), & nsp1(icomp),nsp2(icomp),nsp3(icomp) !spares else read (nread,*) ntermc(icomp),nterme(icomp) ncosh(icomp)=0 !these terms not used in files prior to v6.100 nsinh(icomp)=0 nsp1(icomp)=0 nsp2(icomp)=0 nsp3(icomp)=0 endif jterm=0 if (ntermc(icomp).ge.1) then do j=1,ntermc(icomp) !read polynomial coefficients jterm=jterm+1 read (nread,*) cpc(icomp,jterm),xk(icomp,jterm) enddo end if if (nterme(icomp).ge.1) then do j=1,nterme(icomp) !read exponential coefficients jterm=jterm+1 read (nread,*) cpc(icomp,jterm),xk(icomp,jterm) enddo end if if (ncosh(icomp).ge.1) then do j=1,ncosh(icomp) !read cosh coefficients jterm=jterm+1 read (nread,*) cpc(icomp,jterm),xk(icomp,jterm), & cph(icomp,jterm),xth(icomp,jterm), & xh(icomp,jterm) enddo end if if (nsinh(icomp).ge.1) then do j=1,nsinh(icomp) !read sinh coefficients jterm=jterm+1 read (nread,*) cpc(icomp,jterm),xk(icomp,jterm), & cph(icomp,jterm),xth(icomp,jterm), & xh(icomp,jterm) enddo end if c write (*,*) ' SETCPP--final coeff: ',cpc(icomp,jterm) ierr=0 herr=' ' end if c RETURN end !subroutine SETCPP c c ====================================================================== c function CP0CPP (icomp,t) c c compute Cp0 c c uses polynomial correlation of Cp0: c Cp0/Cred = SUM [cpc(i,k)*(T/Tred(i))**xk(i,k)] c + SUM [cpc(i,k)*u(i,k)**2*exp{u(i,k)}/(1-exp{u(i,k)})**2] c where cpc(i,k), xk(i,k) are k-th coefficients for component i c u(i,k) = xk(i,k)/T c Cred is reducing parameter for Cp0 (e.g. gas constant, R) c Tred is reducing parameter for t (e.g. critical temperature) c (the reducing parameters are specified in the .fld file) c c inputs: c icomp--pointer specifying component (1..nc) c t--temperature (K) c output (as function value): c CP0CPP--Cp0 (J/(mol-K)) c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 10-06-94 MM, original version c 07-20-95 MM, separate polynomial Cp0 routines from core_BWR c 11-10-95 MM, convert from cubic to general polynomial c 11-29-95 MM, variable lower limit on coefficient/constant arrays c to accommodate ECS reference fluid c 02-27-96 MM, parameter n0=-ncmax to accommodate ECS-thermo model c 06-17-96 MM, add provision for exponential terms c 06-18-96 MM, actually add the exponential terms c 04-02-98 EWL, add hyperbolic terms c 06-29-98 MEV, break up expression for the exponential terms to avoid c a possible overflow error. c 01-24-00 EWL, check for ui>1000 to avoid exp(ui)=overflow c 05-02-01 EWL, change u1>1000 to ui>200 c implicit double precision (a-h,o-z) implicit integer (i-k,m,n) implicit logical (l) c 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 (ncppmx=20) !max number of Cp0 terms common /WNTCPP/ ntermc(n0:nx),nterme(n0:nx),ncosh(n0:nx), & nsinh(n0:nx),nsp1(n0:nx),nsp2(n0:nx),nsp3(n0:nx) common /WRDCPP/ tred(n0:nx),Cred(n0:nx) common /WLMCPP/ tmin(n0:nx),tmax(n0:nx),pmax(n0:nx),rhomax(n0:nx) common /WCPCPP/ cpc(n0:nx,ncppmx),xk(n0:nx,ncppmx), & cph(n0:nx,ncppmx),xth(n0:nx,ncppmx), & xh(n0:nx,ncppmx) common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx) c j=icomp !pointer to appropriate fluid c cpsum=0.0d0 iterm=0 c polynomial terms if (ntermc(j).ge.1) then do i=1,ntermc(j) iterm=iterm+1 cpsum=cpsum+cpc(j,iterm)*(t/tred(j))**xk(j,iterm) enddo end if c exponential terms if (nterme(j).ge.1) then do i=1,nterme(j) iterm=iterm+1 ui=xk(j,iterm)/t if (ui.lt.200.0d0) then expui=EXP(ui) cp0cp1=expui/(1.0d0-expui) cp0cp1=cp0cp1/(1.0d0-expui) cpsum=cpsum+cpc(j,iterm)*ui*ui*cp0cp1 endif c cpsum=cpsum+cpc(j,iterm)*ui*ui*expui/(1.0d0-expui)**2 enddo end if c cosh terms if (ncosh(j).ge.1) then do i=1,ncosh(j) iterm=iterm+1 ttau=cph(j,iterm)*(t/tred(j))**xth(j,iterm) if (ttau.lt.700.d0) then cpsum=cpsum+cpc(j,iterm)*(t/tred(j))**xk(j,iterm)* & COSH(ttau)**xh(j,iterm) else cpsum=1.d6 endif enddo end if c sinh terms if (nsinh(j).ge.1) then do i=1,nsinh(j) iterm=iterm+1 ttau=cph(j,iterm)*(t/tred(j))**xth(j,iterm) if (ttau.lt.700.d0) then cpsum=cpsum+cpc(j,iterm)*(t/tred(j))**xk(j,iterm)* & SINH(ttau)**xh(j,iterm) else cpsum=1.d6 endif enddo end if c CP0CPP=cpsum*Cred(j) c write (*,1020) j,t,tred(j),Cred(j),CP0CPP,(cpc(j,i),i=1,iterm) c1020 format (1x,' CP0CPP--j,t,tred,Cred,CP0,cpc(i): ',i3,4f8.2,8d16.6) c RETURN end !function CP0CPP c c ====================================================================== c function CPICPP (icomp,t) c c compute integral of Cp0 over limits of Tref to T c for use in enthalpy calculation c c uses polynomial correlation of Cp0: c Cp0/Cred = SUM [cpc(i,k)*(T/Tred(i))**xk(i,k)] c + SUM [cpc(i,k)*u(i,k)**2*exp{u(i,k)}/(1-exp{u(i,k)})**2] c where cpc(i,k), xk(i,k) are k-th coefficients for component i c u(i,k) = xk(i,k)/T c Cred is reducing parameter for Cp0 (e.g. gas constant, R) c Tred is reducing parameter for t (e.g. critical temperature) c (the reducing parameters are specified in the .fld file) c c based on derivations in Younglove & McLinden (1994), JPCRD 23:731-779 c equation C11 c c inputs: c icomp--pointer specifying component (1..nc); c t--temperature (K) c output: (as function value): c CPICPP--int (Cp0 dT)|T-Tref (J/mol) c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 10-06-94 MM, original version c 07-20-95 MM, separate polynomial Cp0 routines from core_BWR c 11-10-95 MM, convert from cubic to general polynomial c 11-29-95 MM, variable lower limit on coefficient/constant arrays c to accommodate ECS reference fluid c 02-27-96 MM, parameter n0=-ncmax to accommodate ECS-thermo model c 06-17-96 MM, add provision for exponential terms c 06-18-96 MM, actually add the exponential terms c 04-02-98 EWL, add hyperbolic terms c 11-04-02 EWL, include contribution of hyp. terms from the ref. temperature c implicit double precision (a-h,o-z) implicit integer (i-k,m,n) implicit logical (l) c 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 (ncppmx=20) !max number of Cp0 terms common /WNTCPP/ ntermc(n0:nx),nterme(n0:nx),ncosh(n0:nx), & nsinh(n0:nx),nsp1(n0:nx),nsp2(n0:nx),nsp3(n0:nx) common /WRDCPP/ tred(n0:nx),Cred(n0:nx) common /WLMCPP/ tmin(n0:nx),tmax(n0:nx),pmax(n0:nx),rhomax(n0:nx) common /WCPCPP/ cpc(n0:nx,ncppmx),xk(n0:nx,ncppmx), & cph(n0:nx,ncppmx),xth(n0:nx,ncppmx), & xh(n0:nx,ncppmx) common /CREF/ tref(n0:nx),rhoref(n0:nx),href(n0:nx),sref(n0:nx) common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx) c j=icomp !pointer to appropriate fluid c cpsum=0.0d0 CPICPP=0.d0 if (t.le.0) RETURN iterm=0 c polynomial terms if (ntermc(j).ge.1) then do i=1,ntermc(j) iterm=iterm+1 xki=xk(j,iterm) xk1=xki+1.0d0 if (abs(xk1).lt.1.0d-6) then c any term with temperature exponent of -1 handled differently cpsum=cpsum+cpc(j,iterm)*tred(j)*LOG(t/tref(j)) else cpsum=cpsum+cpc(j,iterm)*(t**xk1-tref(j)**xk1) & /(xk1*tred(j)**xki) end if enddo end if c exponential terms if (nterme(j).ge.1) then do i=1,nterme(j) iterm=iterm+1 ui=xk(j,iterm)/t if (ui.lt.200.0d0) then expui=EXP(ui) exptr=EXP(xk(j,iterm)/tref(j)) cpsum=cpsum+cpc(j,iterm)*(-0.5d0*xk(j,iterm)) & *((1.0d0+expui)/(1.0d0-expui)-(1.0d0+exptr)/(1.0d0-exptr)) endif enddo end if c cosh terms, only works for c1/t^2/COSH(c2/t)^2 if (ncosh(j).ge.1) then do i=1,ncosh(j) iterm=iterm+1 cpsum=cpsum-cpc(j,iterm)/cph(j,iterm)*TANH(cph(j,iterm)/t) & +cpc(j,iterm)/cph(j,iterm)*TANH(cph(j,iterm)/tref(j)) enddo end if c sinh terms, only works for c1/t^2/SINH(c2/t)^2 if (nsinh(j).ge.1) then do i=1,nsinh(j) iterm=iterm+1 cpsum=cpsum+cpc(j,iterm)/cph(j,iterm)/TANH(cph(j,iterm)/t) & -cpc(j,iterm)/cph(j,iterm)/TANH(cph(j,iterm)/tref(j)) enddo end if c CPICPP=cpsum*Cred(j) c write (*,1020) j,t,tred(j),Cred(j),CPICPP,(cpc(j,i),i=1,iterm) c1020 format (1x,' CPICPP--j,t,tred,Cred,CPI,cpc(i): ',i3,4f8.2,8d16.6) c RETURN end !function CPICPP c c ====================================================================== c function CPTCPP (icomp,t) c c compute integral of Cp0/T over limits of Tref to T c for use in entropy calculation c c uses polynomial correlation of Cp0: c Cp0/Cred = SUM [cpc(i,k)*(T/Tred(i))**xk(i,k)] c + SUM [cpc(i,k)*u(i,k)**2*exp{u(i,k)}/(1-exp{u(i,k)})**2] c where cpc(i,k), xk(i,k) are k-th coefficients for component i c u(i,k) = xk(i,k)/T c Cred is reducing parameter for Cp0 (e.g. gas constant, R) c Tred is reducing parameter for t (e.g. critical temperature) c (the reducing parameters are specified in the .fld file) c c based on derivations in Younglove & McLinden (1994), JPCRD 23:731-779 c equation C12 c c inputs: c icomp--pointer specifying component (1..nc) c t--temperature (K) c output (as function value): c CPTCPP--int (Cp0/T dT)|T-Tref (J/(mol-K)) c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 10-06-94 MM, original version c 07-20-95 MM, separate polynomial Cp0 routines from core_BWR c 11-10-95 MM, convert from cubic to general polynomial c 11-29-95 MM, variable lower limit on coefficient/constant arrays c to accommodate ECS reference fluid c 02-27-96 MM, parameter n0=-ncmax to accommodate ECS-thermo model c 06-17-96 MM, add provision for exponential terms c 06-18-96 MM, actually add the exponential terms c 04-02-98 EWL, add hyperbolic terms c 11-04-02 EWL, include contribution of hyp. terms from the ref. temperature c implicit double precision (a-h,o-z) implicit integer (i-k,m,n) implicit logical (l) c 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 (ncppmx=20) !max number of Cp0 terms common /WNTCPP/ ntermc(n0:nx),nterme(n0:nx),ncosh(n0:nx), & nsinh(n0:nx),nsp1(n0:nx),nsp2(n0:nx),nsp3(n0:nx) common /WRDCPP/ tred(n0:nx),Cred(n0:nx) common /WLMCPP/ tmin(n0:nx),tmax(n0:nx),pmax(n0:nx),rhomax(n0:nx) common /WCPCPP/ cpc(n0:nx,ncppmx),xk(n0:nx,ncppmx), & cph(n0:nx,ncppmx),xth(n0:nx,ncppmx), & xh(n0:nx,ncppmx) common /CREF/ tref(n0:nx),rhoref(n0:nx),href(n0:nx),sref(n0:nx) common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx) c j=icomp !pointer to appropriate fluid c cpsum=0.0d0 CPTCPP=0.d0 if (t.le.0.d0) RETURN iterm=0 c polynomial terms if (ntermc(j).ge.1) then do i=1,ntermc(j) iterm=iterm+1 xki=xk(j,iterm) if (abs(xki).lt.1.0d-6) then c any term with temperature exponent of 0 handled differently cpsum=cpsum+cpc(j,iterm)*LOG(t/tref(j)) else cpsum=cpsum+cpc(j,iterm)*(t**xki-tref(j)**xki) & /(xki*tred(j)**xki) end if enddo c exponential terms end if if (nterme(j).ge.1) then do i=1,nterme(j) iterm=iterm+1 ui=xk(j,iterm)/t if (ui.lt.200.0d0) then uiref=xk(j,iterm)/tref(j) expui=EXP(ui) exptr=EXP(uiref) cpsum=cpsum+cpc(j,iterm)*(LOG((1.0d0-exptr)/(1.0d0-expui)) & +ui*expui/(expui-1.0d0)-uiref*exptr/(exptr-1.0d0)) endif enddo end if c cosh terms, only works for c1/t^2/COSH(c2/t)^2 if (ncosh(j).ge.1) then do i=1,ncosh(j) iterm=iterm+1 ttau=cph(j,iterm)/t if (ttau.lt.700.d0) then cpsum=cpsum & -cpc(j,iterm)/cph(j,iterm)/t *TANH(ttau) & +cpc(j,iterm)/cph(j,iterm)/tref(j)*TANH(cph(j,iterm)/tref(j)) & +cpc(j,iterm)/cph(j,iterm)**2*LOG(COSH(ttau)) & -cpc(j,iterm)/cph(j,iterm)**2*LOG(COSH(cph(j,iterm)/tref(j))) else cpsum=1.d6 endif enddo end if c sinh terms, only works for c1/t^2/SINH(c2/t)^2 if (nsinh(j).ge.1) then do i=1,nsinh(j) iterm=iterm+1 ttau=cph(j,iterm)/t if (ttau.lt.700.d0) then cpsum=cpsum & +cpc(j,iterm)/cph(j,iterm)/t /TANH(ttau) & -cpc(j,iterm)/cph(j,iterm)/tref(j)/TANH(cph(j,iterm)/tref(j)) & -cpc(j,iterm)/cph(j,iterm)**2*LOG(SINH(ttau)) & +cpc(j,iterm)/cph(j,iterm)**2*LOG(SINH(cph(j,iterm)/tref(j))) else cpsum=1.d6 endif enddo endif c CPTCPP=cpsum*Cred(j) c write (*,1020) j,t,tred(j),Cred(j),CPTCPP,(cpc(j,i),i=iterm) c1020 format (1x,' CPTCPP--j,t,tred,Cred,CPT,cpc(i): ',i3,4f8.2,8d16.6) c RETURN end !function CPTCPP c c ====================================================================== c function PH0CPP (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 density; 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 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 when itau = 0 and idel = 1, 1st density derivative c (actually the derivatives are with respect to the dimensionless c quantities tau and del) c t--temperature (K) c rho--density (mol/L) c output (as function value): c ph0cpp--ideal-gas part of the Helmholtz energy in reduced form (A/RT); c the tau derivatives (as specified by itau) are multiplied by c the corresponding power of tau; i.e. c when itau = 1, the quantity returned is tau*d(ph0cpp)/d(tau) c itau = 2, tau**2*d2(ph0cpp)/d(tau)**2 is returned c itau = 3, tau**3*d3(ph0cpp)/d(tau)**3 is returned c similarly, the del derivatives (as specified by idel) are c multiplied by the corresponding power of del (the derivatives c usually appear with this factor and this approach neatly c avoids a possible divide by zero) c c N.B. 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 08-04-95 MM, original version c 08-21-95 MM, put saved variables into common (rather than save stmt) c 10-03-95 MM, change /MODEL/ + /CPMOD/: models specified by strings c 11-06-95 MM, insert HMX model c 11-29-95 MM, variable lower limit on coefficient/constant arrays c to accommodate ECS reference fluid c 02-27-96 MM, parameter n0=-ncmax to accommodate ECS-thermo model c 03-13-96 MM, replace calls to REDFEQ, etc with generic REDK c 03-21-96 MM, delete /MODEL/, not needed because of above change c 04-18-96 MM, apply tolerance to t-tsav test c 04-19-96 MM, fix bug in expression for PHI; change input p -> rho c 05-10-96 MM, eliminate duplication of href, sref terms with THERM c 07-05-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 08-12-98 MM, add third temperature derivative c N.B. third derivative not yet implemented for hyperbolic terms c 08-13-98 MM, add cross derivatives (all zero) and 3rd density deriv 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 (ncppmx=20) !max number of Cp0 terms character*1 htab,hnull character*255 herr 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) common /CPPSAV/ cp0sav(n0:nx),cpisav(n0:nx),cptsav(n0:nx), & tsav(n0:nx) common /WNTCPP/ ntermc(n0:nx),nterme(n0:nx),ncosh(n0:nx), & nsinh(n0:nx),nsp1(n0:nx),nsp2(n0:nx),nsp3(n0:nx) common /WRDCPP/ tred(n0:nx),Cred(n0:nx) common /WCPCPP/ cpc(n0:nx,ncppmx),xk(n0:nx,ncppmx), & cph(n0:nx,ncppmx),xth(n0:nx,ncppmx), & xh(n0:nx,ncppmx) c PH0CPP=0.0d0 !initialize in case of error if (t.le.0) RETURN if (abs(t-tsav(icomp)).lt.1.0d-8) then c use values from previous call (values already in cp0sav, etc.) c write (*,*) 'PH0CPP--using stored values for itau,t = ',itau,t else c otherwise, compute new values and save for possible future use cp0sav(icomp)=CP0CPP(icomp,t) cpisav(icomp)=CPICPP(icomp,t) cptsav(icomp)=CPTCPP(icomp,t) tsav(icomp)=t end if c Rt=R*t c if (itau*idel.ne.0) then c all cross derivates are zero PH0CPP=0.0d0 else if (itau.eq.0 .and. idel.eq.0) then c compute reduced Helmholtz if (rho.gt.0) PH0CPP=cpisav(icomp)/Rt-cptsav(icomp)/R & +LOG(t*rho/(tref(icomp)*rhoref(icomp))) & -1.0d0 c else if (itau.eq.1) then c compute derivative w.r.t. tau (dimensionless temperature) PH0CPP=cpisav(icomp)/Rt-1.0d0 !return tau*d(ph0cpp)/d(tau) c else if (itau.eq.2) then c compute 2nd derivative w.r.t. tau (dimensionless temperature) c return tau**2*d2(ph0cpp)/d(tau**2) PH0CPP=1.0d0-cp0sav(icomp)/R c else if (itau.eq.3) then j=icomp !pointer to appropriate fluid tau=tz(icomp)/t c note that t0 is reducing temperature for EOS (conversion between c absolute temperature and tau) while tred is a (possible) separate c reducing temperature appearing only in Cp0 PH0sum=2.0d0*(-1.0d0+cp0sav(icomp)/R) iterm=0 c write (*,*) 'PH0CPP--initial PH0sum: ',iterm,PH0sum c polynomial terms if (ntermc(j).ge.1) then do i=1,ntermc(j) iterm=iterm+1 if (ABS(xk(j,iterm)).gt.1.0d-10) then c these terms arise from d(Cp0)/d(tau) and derivative of constant c terms is zero PH0sum=PH0sum+cpc(j,iterm)*xk(j,iterm) & *(tz(icomp)/tred(j))**xk(j,iterm)*tau**(-xk(j,iterm)) end if c write (*,*) 'PH0CPP--i-polynomial,PH0sum: ',i,PH0sum enddo end if c exponential terms if (nterme(j).ge.1) then do i=1,nterme(j) iterm=iterm+1 ui=xk(j,iterm)/t if (ui.lt.200.0d0) then expui=EXP(ui) PH0sum=PH0sum-tau*cpc(j,iterm)*xk(j,iterm)*expui & /(tz(icomp)*(1.0d0-expui)**2) & *(2.0d0*ui+ui*ui*(1.0d0+2.0d0*expui/(1.0d0-expui))) c write (*,*) 'PH0CPP--i-exponential,PH0sum: ',i,PH0sum endif enddo end if c cosh terms c --not yet implemented-- c if (ncosh(j).ge.1) then c do i=1,ncosh(j) c iterm=iterm+1 c enddo c end if c sinh terms c --not yet implemented-- c if (nsinh(j).ge.1) then c do i=1,nsinh(j) c iterm=iterm+1 c enddo c end if PH0CPP=PH0sum c else if (idel.eq.1) then c the three density derivatives are 1/del, -1/del**2, and 1/del**3 c the quantity returned is derivative multiplied by corresponding power c of del, avoiding a possible divide by zero PH0CPP=1.0d0 else if (idel.eq.2) then PH0CPP=-1.0d0 else if (idel.eq.3) then PH0CPP=2.0d0 c else c invalid value of itau ierr=99 write (herr,1099) itau,idel,hnull 1099 format ('[PH0CPP warning] invalid input; itau =',i4,'; idel =', & i4,a1) call ERRMSG (ierr,herr) PH0CPP=0.0d0 end if c c write (*,*) ' PH0CPP: output phi: ',ph0cpp c RETURN end !function PH0CPP cc cc ====================================================================== cc c block data BDCPP cc cc data for polynomial form of ideal gas heat capacity cc Cp0/Cred = SUM [cpc(i,k)*(T/Tred(i))**xk(i,k)] cc where Cred is reducing parameter for Cp0 (e.g. gas constant, R) cc Tred is reducing parameter for t (e.g. critical temperature) cc c implicit double precision (a-h,o-z) c parameter (mxcpp=2) !max number of fluids in block data c parameter (ncppmx=20) !max number of Cp0 terms c character*12 hcas c common /NTCPP/ ntrmca(mxcpp),ntrmce(mxcpp) c common /RDCPP/ treda(mxcpp),Creda(mxcpp) c common /LMCPP/ tmna(mxcpp),tmxa(mxcpp),pmxa(mxcpp),rhomxa(mxcpp) c common /CPCPP/ cpca(mxcpp,ncppmx),xka(mxcpp,ncppmx) c common /CASCPP/ hcas(mxcpp) cc cc explanation of commons and constituent arrays cc /NTCPP/ number of terms in polynomial fit cc /RDCPP/ reducing parameters for temperature, Cp0 (often =1) cc /LMCPP/ limits of fit: tmin, tmax, pmax (= 0), rhomax (= 0) cc /CPCPP/ parameters to polynomial Cp0 fit cc /CASCPP/ Chem Abstract number; used as unambiguous identifier cc hcas(i): CAS number for fluid corresponding to equation "i" cc cc cc written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado cc 07-20-95 MM, original version cc 11-10-95 MM, convert from cubic to general polynomial cc 06-17-96 MM, add provision for exponential terms cc cc N.B.--the "i" have no particular significance cc cc cc R134a 1,1,1,2-tetrafluoroethane c data hcas(1) /'811-97-2'/ cc fit of McLinden et al. (1989), ASHRAE Trans 95(pt 2):263-283 cc also in Huber & McLinden (1992), Int Refrig Conf, Purdue, 453-462 c data ntrmca(1),ntrmce(1) /3,0/ c data treda(1),Creda(1) /1.0d0,1.0d0/ c data tmna(1),tmxa(1),pmxa(1),rhomxa(1) c & /150.0d0,500.0d0,0.0d0,0.0d0/ c data (xka(1,i),i=1,3) c & / 0.0d0, 1.0d0, 2.0d0/ c data (cpca(1,i),i=1,3) c & /19.4006d0,0.258531d0,-1.29665d-4/ cc cc R123 2,2-dichloro-1,1,1-trifluoroethane c data hcas(2) /'306-83-2'/ cc fit of Younglove & McLinden (1994), JPCRD 23:731-779 c data ntrmca(2),ntrmce(2) /4,0/ c data treda(2),Creda(2) /456.831d0,8.31451d0/ c data tmna(2),tmxa(2),pmxa(2),rhomxa(2) c & /150.0d0,500.0d0,0.0d0,0.0d0/ c data (xka(2,i),i=1,4) c & /0.0d0, 1.0d0, 2.0d0, 3.0d0/ c data (cpca(2,i),i=1,4) c & /2.046009d0,22.231991d0,-11.658491d0,2.691665/ cc c end !block data BDCPP c c c 1 2 3 4 5 6 7 c23456789012345678901234567890123456789012345678901234567890123456789012 c c ====================================================================== c end file core_CPP.f c ======================================================================