c begin file mix_HMX.f c c This file contains the routines implementing the mixture Helmholtz c equation of state of Lemmon and Jacobsen. c c contained here are: c function PHIHMX (itau,idel,tau,del,x) c function PHIMIX (jcomp,x,itau,idel,tau,del,x) c subroutine CRTHMX (x,tcrit,pcrit,Dcrit,ierr,herr) c subroutine REDHMX (x,t0,D0) c subroutine SETHMX (hfmix,ierr,herr) c subroutine RDBNC (nread,icomp,jcomp,ibin,ierr,herr) c subroutine RDMIX (nread,icomp,jcomp,hmodij,lij,ierr,herr) c subroutine ESTBNC (i,j,ibin,ierr,herr) c c ====================================================================== c ====================================================================== c function PHIHMX (itau,idel,tau,del,x) c c compute reduced Helmholtz energy or a derivative as functions c of dimensionless temperature and density for the mixture Helmholtz c equation of state c c based on notes of E.W. Lemmon and R.T Jacobsen, U. Idaho, 1995 c c inputs: 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 x--composition array (mol frac) 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 gas and residual (real- c gas) terms. The residual term consists of ideal-solution and c mixing terms. This routine calculates only the residal term. c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 11-02-95 MM, original version c 11-26-95 MM, rearrange argument list (x at end) c 11-29-95 MM, variable lower limit on coefficient/constant arrays c to accomodate ECS reference fluid c 02-27-96 MM, parameter n0=-ncmax to accomodate ECS-thermo model c 02-29-96 MM, replace calls to PHIBWR, PHIFEQ with general PHIK c 03-21-96 MM, delete reference to /MODEL/, not used c 07-01-96 MM, move x(i)*x(j) mulitplier here (remove from PHIMIX) c 06-03-97 EWL, add third derivative of Helmholtz energy with respect to density. c implicit double precision (a-h,o-z) implicit integer (i-k,m,n) c parameter (ncmax=5) !max number of components in mixture common /NCOMP/ nc dimension x(ncmax) c if (del.le.1.0d-10) then !trivial solution at zero density phihmx=0.0d0 !for any and all derivatives RETURN end if c phisum=0.0d0 do 400 i=1,nc if (x(i).gt.0.0d0) then c calculate ideal solution term phiii=PHIK(i,itau,idel,tau,del) phisum=phisum+x(i)*phiii c write (*,1200) i,i,x(1),x(2),x(3),phiii c1200 format (1x,' PHIHMX--i,j,x,phi_ij: ',2i4,3f8.5,1pe14.6) if (i.le.nc-1) then c calculate mixing term do 200 j=i+1,nc phiij=PHIMIX(i,j,itau,idel,tau,del) phisum=phisum+x(i)*x(j)*phiij c write (*,1200) i,j,x(1),x(2),x(3),phiij 200 continue end if end if 400 continue phihmx=phisum c RETURN end !function PHIHMX c c ====================================================================== c function PHIMIX (icomp,jcomp,itau,idel,tau,del) c c compute reduced Helmholtz energy of mixing (or a derivative) c for the binary interaction of components i and j as a function of c composition and dimensionless temperature and density for the c mixture Helmholtz equation of state c c inputs: c icomp--component i c jcomp--component j 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 Amix/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--mixture interaction (excess) part of the Helmholtz energy, c or one of its derivatives (as specified by itau and idel), c in reduced form (Amix/RT) c c The Helmholtz energy consists of ideal gas and residual (real- c gas) terms. The residual term consists of ideal-solution and c mixing terms. This routine calculates only the mixing term. c c N.B. This routine incorporates the mixture-specific multiplier c (e.g. the Fpq of the Lemmon & Jacobsen model). c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 11-02-95 MM, original version c 11-26-95 MM, rearrange argument list (x at end) c 11-29-95 MM, variable lower limit on coefficient/constant arrays c to accomodate ECS reference fluid c 01-10-96 MM, add commons /CFXHMX/ and /CFIHMX/ c 02-27-96 MM, parameter n0=-ncmax to accomodate ECS-thermo model c 03-21-96 MM, delete reference to /MODEL/, not used c 05-06-96 MM, add /MIXMOD/, cases for different mix models c 06-12-96 MM, add /HMXSAV/, LJi and LMi models c 07-01-96 MM, move x(i)*x(j) multiplier to PHIHMX; remove x from inputs c 11-04-96 MM, replace van der Waals rule with linear (LIN), c change nmxpar from 4 to 6 c implicit double precision (a-h,o-z) implicit integer (i-k,m,n) c parameter (ncmax=5) !max number of components in mixture c parameter (n0=-ncmax,nx=ncmax) parameter (nx=ncmax) parameter (nmxpar=6) !number of binary mixture parameters parameter (nmxtrm=15) !number of terms in binary mixing rule character*3 hmodmx dimension phi01(nmxtrm),phi10(nmxtrm) common /NCOMP/ ncomp common /MIXMOD/ hmodmx(nx,nx) c commons associated with the binary mixing rule(s) c the amix(i,j,k,1..3) are the coefs for the (general) mixing term c the fmix(i,j,1..nmxpar) are the parameters for the i-j binary c the namix(i,j) are the number of a(i,j,k,1..3) terms c the iexpmx(i,j,k) are powers of the exponentials in the general term common /CFXHMX/ amix(nx,nx,nmxtrm,3),fmix(nx,nx,nmxpar) common /CFIHMX/ namix(nx,nx),iexpmx(nx,nx,nmxtrm) common /HMXSAV/ phisav(nx,nx,0:nmxtrm),tausav(nx,nx),delsav(nx,nx) c if (del.le.1.0d-10) then !trivial solution at zero density phimix=0.0d0 !for any and all derivatives RETURN end if c i=icomp j=jcomp if (hmodmx(i,j) .eq. 'LIN') then c c linear mixing rules: all departures from ideal solution are c contained in the reducing parameters (i.e. Kt and/or Kv) c phimix=0.0d0 c else if (hmodmx(i,j)(1:2).eq.'LJ' .or. & hmodmx(i,j)(1:2).eq.'LM') then c c Lemmon-Jacobsen or modified Lemmon-Jacobsen mixing function c if (abs(tau-tausav(i,j)).lt.1.0d-12 .and. & abs(del-delsav(i,j)).lt.1.0d-12) then c retrieve value from previous call phimix=phisav(i,j,0) ! zero element stores summation else c otherwise, compute new values and save for possible future use phisum=0.0d0 do 100 k=1,namix(i,j) savijk=amix(i,j,k,1)*tau**amix(i,j,k,2)*del**amix(i,j,k,3) c if (iexpmx(i,j,k) .ne. 0) then c savijk=savijk*exp(-del**iexpmx(i,j,k)) c end if phisav(i,j,k)=savijk phisum=phisum+savijk c write (*,1010) k,phisav(i,j,k),phisum c1010 format (1x,i3,2d30.20) !write out each term for debugging 100 continue tausav(i,j)=tau delsav(i,j)=del phisav(i,j,0)=fmix(i,j,3)*phisum !fmix(i,j,3) is Fpq PHIMIX=phisav(i,j,0) end if 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 101 k=1,namix(i,j) phi01(k)=phisav(i,j,k)*amix(i,j,k,3)/del c if (iexpmx(i,j,k) .ne. 0) then c phi01(k)=phi01(k) c & *(amix(i,j,k,3)-iexpmx(i,j,k)*del**iexpmx(i,j,k)) c end if phisum=phisum+phi01(k) 101 continue PHIMIX=fmix(i,j,3)*phisum !fmix(i,j,3) is Fpq else if (idel.eq.2) then c compute 2nd derivative w.r.t. del (dimensionless density) phisum=0.0d0 do 102 k=1,namix(i,j) ax3=amix(i,j,k,3) phisum=phisum+phisav(i,j,k)*ax3*(ax3-1.0d0)/del**2 c if (iexpmx(i,j,k) .ne. 0) then c dell=iexpmx(i,j,k)*del**iexpmx(i,j,k) c phisum=phisum+phisav(i,j,k)/del**2* c & (ax3**2-ax3-dell*(2.0d0*ax3+iexpmx(i,j,k)-1.0d0-dell)) c end if 102 continue PHIMIX=fmix(i,j,3)*phisum !fmix(i,j,3) is Fpq else if (idel.eq.3) then c compute 3nd derivative w.r.t. del (dimensionless density) phisum=0.0d0 do 103 k=1,namix(i,j) ax3=amix(i,j,k,3) phisum=phisum+phisav(i,j,k)*ax3*(ax3-1.0d0)*(ax3-2.0d0)/del**3 103 continue PHIMIX=fmix(i,j,3)*phisum !fmix(i,j,3) is Fpq end if 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 110 k=1,namix(i,j) phi10(k)=phisav(i,j,k)*amix(i,j,k,2)/tau phisum=phisum+phi10(k) 110 continue PHIMIX=fmix(i,j,3)*phisum !fmix(i,j,3) is Fpq else if (itau.eq.2) then c compute 2nd derivative w.r.t. tau (dimensionless temperature) phisum=0.0d0 do 120 k=1,namix(i,j) ax2=amix(i,j,k,2) phisum=phisum+phisav(i,j,k)*ax2*(ax2-1.0d0)/tau**2 120 continue PHIMIX=fmix(i,j,3)*phisum !fmix(i,j,3) is Fpq end if c if ((itau.eq.1) .and. (idel.eq.1)) then c compute cross derivative using terms from 1st derivatives phisum=0.0d0 do 211 k=1,namix(i,j) phisum=phisum+phi10(k)*phi01(k)/phisav(i,j,k) 211 continue PHIMIX=fmix(i,j,3)*phisum !fmix(i,j,3) is Fpq end if else c c space for additional mixing rule c PHIMIX=0.0d0 !not currently implemented * write (*,*) ' PHIMIX ERROR--unknown mixing rule' end if c RETURN end !function PHIMIX c c ====================================================================== c subroutine CRTHMX (x,tcrit,pcrit,Dcrit,ierr,herr) c c returns critical parameters associated with mixture Helmholtz EOS c c N.B. these critical parameters are estimates based on polynomial c fits to the binary critical lines and (for 3 or more components) c combining rules applied to the constituent binaries c c input: c x--composition array (mol frac) c outputs: c tcrit--critical temperature (K) c pcrit--critical pressure (kPa) c Dcrit--molar density (mol/L) at critical point c ierr--error flag: 0 = successful c 1 = did not converge c herr--error string (character*255 variable if ierr<>0) c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 11-02-95 MM, original version c 11-29-95 MM, variable lower limit on coefficient/constant arrays c to accomodate ECS reference fluid c 01-29-96 MM, critical parameters based on average of pure fluid c 02-27-96 MM, parameter n0=-ncmax to accomodate ECS-thermo model c add Zcrit to common /CCON/ c 03-19-19 MM, add dipole moment to /CCON/ c 05-06-96 MM, replace temporary func for crit pars with reducing pars c 11-19-96 MM, add empirical function for critical line(s) c 11-14-97 EWL, fix bug in critical line correlation c 11-24-97 MM, Tsum = x(i)*x(j)*Tij, not xij*tcij, ditto for Vsum c implicit double precision (a-h,o-z) implicit integer (i-k,m,n) implicit logical (l) c parameter (ncmax=5) !max number of components in mixture parameter (n0=-ncmax,nx=ncmax) parameter (nmxtc=11,nmxvc=11) !number of coeff for critical lines character*1 htab,hnull character*3 hmodtc,hmodvc character*255 herr common /NCOMP/ ncomp common /CCON/ wm(n0:nx),ttp(n0:nx),tnbp(n0:nx), & tc(n0:nx),pc(n0:nx),rhoc(n0:nx),Zcrit(n0:nx), & accen(n0:nx),dipole(n0:nx) common /CHAR/ htab,hnull common /Gcnst/ R c commons associated with the binary critical lines common /CRTMOD/ hmodtc(nx,nx),hmodvc(nx,nx) common /CFXCRT/ ftcij(nx,nx,nmxtc),fvcij(nx,nx,nmxvc) common /CFICRT/ nftcij(nx,nx),nfvcij(nx,nx) c dimension x(ncmax) ierr=0 herr=hnull c tsum=0.0d0 Vsum=0.0d0 xsum=0.0d0 c write (*,1010) (x(k),k=1,ncomp),(tc(k),k=1,ncomp) c1010 format (' CRTHMX input x,Tc: ',10f10.4) do 200 i=1,ncomp-1 do 160 j=i+1,ncomp xsum=xsum+x(i)*x(j) c the tcij and Vcij are the binary critical line correlations, and c these terms (ONLY) are evaluated at the pseudo-binary compositions c defined by xi and xj: xi=0.5d0*(x(i)+(1.0d0-x(j))) xj=0.5d0*(x(j)+(1.0d0-x(i))) xij=xi*xj if (hmodtc(i,j).eq.'TC1') then tcij=xi*tc(i)+xj*tc(j)+xij*ftcij(i,j,1) c terms assuming mixture specified in same order as in HMX.bnc file & +xj*xi*((ftcij(i,j,2)+(ftcij(i,j,3)+(ftcij(i,j,4) & +(ftcij(i,j,5)+ftcij(i,j,6)*xi)*xi)*xi)*xi)*xi) c extra terms to allow specifying mix in opposite order (j,i) & +xi*xj*((ftcij(i,j,7)+(ftcij(i,j,8)+(ftcij(i,j,9) & +(ftcij(i,j,10)+ftcij(i,j,11)*xj)*xj)*xj)*xj)*xj) else ierr=1 herr=' error from CRTHMX, Tcrit model not found'//hnull tcij=xi*tc(i)+xj*tc(j) end if tsum=tsum+x(i)*x(j)*tcij if (hmodvc(i,j).eq.'VC1') then Vcij=xi/rhoc(i)+xj/rhoc(j)+xij*fvcij(i,j,1) c terms assuming mixture specified in same order as in HMX.bnc file & +xj*xi*((fvcij(i,j,2)+(fvcij(i,j,3)+(fvcij(i,j,4) & +(fvcij(i,j,5)+fvcij(i,j,6)*xi)*xi)*xi)*xi)*xi) c extra terms to allow specifying mix in opposite order (j,i) & +xi*xj*((fvcij(i,j,7)+(fvcij(i,j,8)+(fvcij(i,j,9) & +(fvcij(i,j,10)+fvcij(i,j,11)*xj)*xj)*xj)*xj)*xj) else ierr=1 herr=' error from CRTHMX, Vcrit model not found'//hnull Vcij=xi/rhoc(i)+xj/rhoc(j) end if Vsum=Vsum+x(i)*x(j)*Vcij c write (*,1016) i,j,tcij,Vcij,xsum c1016 format (' CRTHMX: i,j,Tc_ij,Vc_ij,xsum: ',2i3,3f12.6) c write (*,1017) i,j,(ftcij(i,j,k),k=1,nmxtc) c1017 format (' CRTHMX: i,j,ftcij: ',2i3,8f10.6) 160 continue 200 continue if (xsum.gt.1.0d-10) then tcrit=tsum/xsum Dcrit=xsum/Vsum else c xsum = 0 implies (ncomp - 1) compositions are zero do 240 i=1,ncomp if (abs(x(i)-1.0d0).lt.1.0d-8) then tcrit=tc(i) Dcrit=rhoc(i) end if 240 continue end if c base Pc on Tc, Dc c call REDHMX (x,tred,Dred) c tau=tred/tcrit c del=Dcrit/Dred c phi01=PHIHMX (0,1,tau,del,x) c pcrit=R*tcrit*Dcrit*(1.0d0+del*phi01) call PRESS (tcrit,Dcrit,x,pcrit) c write (*,1024) tcrit,1.0/Dcrit,pcrit c1024 format (' CRTHMX: Tcrit,Vcrit,Pcrit: ',15x,3f12.6) c RETURN end !subroutine CRTHMX c c ====================================================================== c subroutine REDHMX (x,tred,Dred) c c returns reducing parameters associated with mixture Helmholtz EOS; c used to calculate the 'tau' and 'del' which are the independent c variables in the EOS c c input: c x--composition array [mol frac] 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 11-02-95 MM, original version c 11-29-95 MM, variable lower limit on coefficient/constant arrays c to accomodate ECS reference fluid c 02-27-96 MM, parameter n0=-ncmax to accomodate ECS-thermo model c add Zcrit to common /CCON/ c 03-19-19 MM, add dipole moment to /CCON/ c 03-21-96 MM, replace call to REDFEQ with general REDX, delete /MODEL/ c 05-06-96 MM, add /MIXMOD/,/CFXHMX/,/CFIHMX/; add Kt, Kv dependence; c branch for arithmetic or geometric mean of temperature c 05-14-96 MM, modify volume mixing rule for VDW model; use arithmetic c average of critical volumes, rather than Reiner's rule c 06-10-96 MM, modify volume mixing rule for VDW model (again); use c average of cube root of critical volumes c add Lemmon-Jacobsen mixing rules c 06-12-96 MM, add LMx mixing rules (Lemmon-Jacobsen modified to be of c form compatible with van der Waal's (VDW) model c 07-01-96 MM, fix bug--missing **3 in VDW volume term c 11-04-96 MM, replace van der Waals rule with linear (LIN), c change nmxpar from 4 to 6 c 11-08-96 MM, fix wrong array reference for Kv c 11-14-96 MM, fix wrong array reference for zeta, xi (LJi model) c 05-22-97 MM, bug fix: beta, gamma were applied to V, should not be c implicit double precision (a-h,o-z) implicit integer (i-k,m,n) c parameter (ncmax=5) !max number of components in mixture parameter (n0=-ncmax,nx=ncmax) parameter (nmxpar=6) !number of binary mixture parameters parameter (nmxtrm=15) !number of terms in binary mixing func character*3 hmodmx common /NCOMP/ ncomp common /CCON/ wm(n0:nx),ttp(n0:nx),tnbp(n0:nx), & tc(n0:nx),pc(n0:nx),rhoc(n0:nx),Zcrit(n0:nx), & accen(n0:nx),dipole(n0:nx) common /MIXMOD/ hmodmx(nx,nx) common /CFXHMX/ amix(nx,nx,nmxtrm,3),fmix(nx,nx,nmxpar) common /CFIHMX/ namix(nx,nx),iexpmx(nx,nx,nmxtrm) dimension x(ncmax),treda(ncmax),Dreda(ncmax) c tsum=0.0d0 Vsum=0.0d0 do 100 i=1,ncomp call REDK (i,t0,D0) treda(i)=t0 Dreda(i)=D0 100 continue c do 160 i=1,ncomp do 120 j=i,ncomp if (i.eq.j) then c special case when i = j c write (*,*) ' REDHMX--special case for i = j',i,j tsum=tsum+x(i)*treda(i) Vsum=Vsum+x(i)/Dreda(i) else if (hmodmx(i,j).eq.'LIN') then c simple linear mixing rules xij=x(i)*x(j) t12=0.5d0*(treda(i)+treda(j)) V12=0.5d0*(1.0d0/Dreda(i)+1.0d0/Dreda(j)) tsum=tsum+xij*(fmix(i,j,1)-1.0d0)*t12 !fmix(i,j,1) = Kt(i,j) Vsum=Vsum+xij*(fmix(i,j,2)-1.0d0)*V12 !fmix(i,j,2) = Kv(i,j) else if (hmodmx(i,j)(1:2) .eq. 'LM') then c Lemmon-Jacobsen rule modified to use Kt, Kv, rather than zeta, xi if (x(i).gt.1.0d-10 .and. x(j).gt.1.0d-10) then xij=x(i)*x(j) xijt=x(i)**fmix(i,j,4)*x(j)**fmix(i,j,5)!fmix(i,j,4) = beta t12=0.5d0*(treda(i)+treda(j)) !fmix(i,j,5) = gamma V12=0.5d0*(1.0d0/Dreda(i)+1.0d0/Dreda(j)) tsum=tsum+xijt*(fmix(i,j,1)-1.0d0)*t12 !fmix(i,j,1) = Kt(i,j) Vsum=Vsum+xij*(fmix(i,j,2)-1.0d0)*V12 !fmix(i,j,2) = Kv(i,j) else c zero composition for one or both components, no change in sums end if else if (hmodmx(i,j)(1:2) .eq. 'LJ') then c original Lemmon-Jacobsen rule using zeta, xi if (x(i).gt.1.0d-12 .and. x(j).gt.1.0d-12) then xij=x(i)*x(j) xijt=x(i)**fmix(i,j,4)*x(j)**fmix(i,j,5)!fmix(i,j,5) = gamma tsum=tsum+xijt*fmix(i,j,1) !fmix(i,j,1) = zeta(T) Vsum=Vsum+xij*fmix(i,j,2) !fmix(i,j,2) = xi (x) else c zero composition for one or both components, no change in sums end if else c mixing rule not specified or not found--should not get here * write (*,*) ' REDHMX ERROR--mixing rule not found for i,j =',i,j end if 120 continue 160 continue tred=tsum Dred=1.0d0/Vsum c write (*,1120) x(1),tred,Dred c1120 format (' REDHMX--x,t0,D0: ',f10.7,f10.3,f10.5) c RETURN end !subroutine REDHMX c c ====================================================================== c subroutine SETHMX (hfmix,ierr,herr) c c read parameters for the HMX mixture model from a file c c input: c hfmix--file name containing mixture parameters [character*80] 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 [fluid parameters, etc. returned via various common blocks] c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 01-09-96 MM, original version 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/ c 03-22-96 MM, replace /MODEL/ with /EOSMOD/ c 05-06-96 MM, add /MIXMOD/ to specify model for binary pairs; c initialize j,i same as i,j; initialize i = j to 1.0 c 05-07-96 MM, move read of binary model(s) to RDMIX c 05-08-96 MM, add /MXINFO/ and read corresponding info from file c 11-19-96 MM, move initialization of common blocks to RDBNC, ESTBNC c 03-27-97 MM, if file not found, search in likely directories c 07-03-97 MM, pass any error from ESTBNC up the chain c 08-22-97 MM, use double backslash on file open to avoid problem with Unix machines c 01-09-98 MM, save any errors from early calls to ESTBNC, etc and pass as output 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 (nbrule=10) !number of binary mixing rules parameter (nbin=ncmax*(ncmax-1)/2) !# possible binary pairs parameter (nmxpar=6) !number of binary mixture parameters parameter (nmxtrm=15) !number of terms in binary mixing rule parameter (nmxtc=11,nmxvc=11) !number of coeff for critical lines character*1 htab,hnull character*3 hpheq,heos,hmxeos,hmodcp character*3 hmodmx,hmodtc,hmodvc character*3 htype character*12 hcas character*80 hfmix,hfilei character*255 herr,herr1,herr2 c next 3 declarations associated with /MXINFO/ character*8 hbpar character*80 hmfile,hrule character*255 hbin c logical lmix,lij(nx,nx),lfbnc,lftc,lfvc c common /NCOMP/ ncomp common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx) common /CCAS/ hcas(n0:nx) common /CHAR/ htab,hnull c commons associated with the binary mixing rule(s) c the amix(i,j,k,1..3) are the coefs for the (general) mixing term c the fmix(i,j,1..nmxpar) are the parameters for the i-j binary c the namix(i,j) are the number of a(i,j,k,1..3) terms c the iexpmx(i,j,k) are powers of the exponentials in the general term common /MIXMOD/ hmodmx(nx,nx) common /CFXHMX/ amix(nx,nx,nmxtrm,3),fmix(nx,nx,nmxpar) common /CFIHMX/ namix(nx,nx),iexpmx(nx,nx,nmxtrm) c /MXINFO/ contains information on the mixing rules and parameters c hmfile specifies the files from which mixing rules originate c (stored in order 1,2 1,3 2,3; element zero is file called in SETUP) c hbin provides documentation for the current binary parameters c hrule contains descriptions of the currently available mixing rules c hbpar contains descriptions of the binary parameters (e.g. Kt, Kv) c associated with the currently available mixing rules common /MXINFO/ hmfile(0:nbin),hbin(nbin),hrule(nbrule), & hbpar(nbrule,nmxpar) c commons associated with the binary critical lines common /CRTMOD/ hmodtc(nx,nx),hmodvc(nx,nx) common /CFXCRT/ ftcij(nx,nx,nmxtc),fvcij(nx,nx,nmxvc) common /CFICRT/ nftcij(nx,nx),nfvcij(nx,nx) c ierr=0 herr=hnull c c initialize the diagonal terms (i = j) c do 120 i=1,ncomp hmodmx(i,i)='i=j' hmodtc(i,i)='i=j' hmodvc(i,i)='i=j' c the i = j terms in the fmix summation should not be called (neither c the i > j), but set the number of terms to zero as precaution namix(i,i)=0 do 100 k=1,nmxpar fmix(i,i,k)=0.0d0 100 continue nftcij(i,i)=0 do 110 k=1,nmxtc ftcij(i,i,k)=0.0d0 110 continue nfvcij(i,i)=0 do 116 k=1,nmxvc fvcij(i,i,k)=0.0d0 116 continue 120 continue c if (hfmix(1:5).eq.'BDATA' .or. hfmix(1:5).eq.'bdata') then c get coefficients from common blocks--not currently implemented ierr=1 herr=' SETHMX ERROR--block data option not available'//hnull call ERRMSG (ierr,herr) RETURN else c c read parameters for the binary pairs from file c * write (*,1105) hfmix *1105 format (/1x,'reading mixture coefficients from file: (',a40,')') nread=12 !logical unit for file reads c open (unit=nread,file=hfmix,err=999) open (unit=nread,file=hfmix,status='old',iostat=ioerr) if (ioerr.ne.0) then c if error in opening file, try looking in likely directory c use double backslash to avoid problem with Unix machines hfilei='\refprop\fluids\\'//hfmix open (unit=nread,file=hfilei,status='old',iostat=ioerr) if (ioerr.ne.0) then hfilei='MPW:code:Refprop613'//hfmix open (unit=nread,file=hfilei,status='old',err=999) end if * write (*,1106) hfmix,hfilei *1106 format (1x,'initial file specification failed: ',a80/ * & 1x,'reading coefficients from file: ',a80) end if rewind (nread) read (nread,2003) htype if (htype.ne.heos) then ierr=1 herr=' SETHMX ERROR--mixture file does not match model'//hnull call ERRMSG (ierr,herr) RETURN end if ibin=0 !initialize counter for binary pair number do 200 i=1,ncomp-1 do 199 j=i+1,ncomp ibin=ibin+1 hmfile(ibin)=hfmix call RDBNC (nread,i,j,ibin,ierr,herr) c check if any binary parameters need to be estimated lfbnc=.true. lftc=.true. lfvc=.true. if (ierr.ne.0) then if (ierr.gt.0) then lfbnc=.true. lftc=.true. lfvc=.true. else ierr1=-ierr-80 if (mod(ierr1,4).eq.0) then lfvc=.true. ierr1=ierr1-4 end if if (mod(ierr1,2).eq.0) then lftc=.true. ierr1=ierr1-2 end if if (ierr1.eq.1) then lfbnc=.true. end if end if if (lfbnc) then call ESTBNC (i,j,ibin,ierr1,herr1) hmfile(ibin)='estimated by ESTBNC' end if if (lftc) then call ESTTC (i,j,ierr2,herr2) end if if (lfvc) then call ESTVC (i,j,ierr2,herr2) end if c any error from ESTBNC takes precedence over those from ESTTC or ESTVC if (ierr1.ne.0) then ierr=ierr1 herr=herr1 else ierr=ierr2 herr=herr2 end if end if 199 continue 200 continue c c read coefficients for general mixing rule(s) from file c call RDMIX (nread,icomp,jcomp,'ALL',lij,ierr1,herr1) c check that mixing term for each binary has been found lmix=.true. do 244 i=1,ncomp-1 do 240 j=i+1,ncomp lmix=lmix.and.lij(i,j) * write (*,1024) i,j,lij(i,j),lmix *1024 format (1x,' SETHMX i,j,lij,lmix: ',2i3,1x,2i2) 240 continue 244 continue if (.not.lmix) then ierr=-117 herr=' SETHMX ERROR--mixing term not found for one or '// & 'more binary pair(s)'//hnull call ERRMSG (ierr,herr) end if c rewind and close the file rewind (nread) close (nread) c error from RDMIX takes precedence over those from ESTBNC, ESTTC, ESTVC if (ierr1.ne.0) then ierr=ierr1 herr=herr1 end if c c write all parameters for debugging * write (*,*) * write (*,*) ' SETHMX: mixing parameters' * write (*,1281) *1281 format (' i j mod par_1 par_2 par_3', * & ' par_4 par_5 par_6'/) d do 282 i=1,ncomp d do 280 j=1,ncomp * write (*,1280) i,j,hmodmx(i,j),(fmix(i,j,k),k=1,nmxpar) d 280 continue d 282 continue *1280 format (1x,2i3,2x,a3,11f10.5) * write (*,*) * write (*,*) ' SETHMX: critical line (Tc) parameters' * write (*,1283) *1283 format (' i j mod par_1 par_2 par_3', * & ' par_4 par_5 par_6 par_7', * & ' par_8 par_9 par_10 par_11'/) d do 286 i=1,ncomp d do 284 j=1,ncomp * write (*,1280) i,j,hmodtc(i,j),(ftcij(i,j,k),k=1,nmxtc) d 284 continue d 286 continue * write (*,*) * write (*,*) ' SETHMX: critical line (Vc) parameters' * write (*,1283) d do 290 i=1,ncomp d do 288 j=1,ncomp * write (*,1280) i,j,hmodvc(i,j),(fvcij(i,j,k),k=1,nmxvc) d 288 continue d 290 continue RETURN end if c c error terminations c 999 ierr=111 herr=' SETHMX ERROR in opening mixture file'//hnull call ERRMSG (ierr,herr) RETURN c 2003 format (a3) c end !subroutine SETHMX c c ====================================================================== c subroutine RDBNC (nread,icomp,jcomp,ibin,ierr,herr) c c read binary parameters for the HMX mixture model and also the c critical line model from a file 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 n (file should have already c been opened and pointer set by subroutine SETHMX) c icomp--component i c jcomp--component j c ibin--counter for binary pair (1,2 = 1; 1,3 = 2; 2,3 = 3; etc.) c outputs: c ierr--error flag: 0 = successful c -114 = critical temperature line will be estimated c -115 = critical volume line will be estimated c -116 = critical temperature and volume lines will be estimated c -117 = binary pair not found, all parameters will be estimated c 118 = error in reading mixture file c herr--error string (character*255 variable if ierr<>0) c other outputs written to common c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 01-10-96 MM, original version c 01-12-96 MM, call ESTBNC if binary parameter not found in file c 02-27-96 MM, parameter n0=-ncmax to accomodate ECS-thermo model c 05-07-96 MM, change inputs from CAS numbers to icomp, jcomp c 05-08-96 MM, add /MXINFO/ and read corresponding info from file c 11-04-96 MM, add gamma to LJi and LMi models (i,j not eq to j,i), c change nmxpar from 4 to 6 c 11-19-96 MM, move initialization of common blocks to RDBNC, ESTBNC c read critical line coefficients c 02-19-97 MM, fix bug: parameter name is ncvc, not ncvv c 05-22-97 MM, change format to read coeff from f10.x to f12.x c 06-04-97 EWL, change value of nctc and ncvc to read in six coefficients c implicit double precision (a-h,o-z) implicit integer (i-k,m,n) implicit logical (l) c parameter (ncmax=5) !max number of components in mixture parameter (n0=-ncmax,nx=ncmax) parameter (nbrule=10) !number of binary mixing rules parameter (nbin=ncmax*(ncmax-1)/2) !# possible binary pairs parameter (nmxpar=6) !number of binary mixture parameters parameter (nmxtrm=15) !number of terms in binary mixing func parameter (nmxtc=11,nmxvc=11) !number of coeff for critical lines parameter (nctc=(nmxtc+1)/2,ncvc=(nmxvc+1)/2) !no crit in file character*1 htab,hnull character*1 hstar,h1 character*3 hmodmx,hmodij,hmodt,hmodtc,hmodv,hmodvc character*3 hflag,htype character*12 hcas,hcas1,hcas2,hcasi,hcasj character*255 herr,hbinp c next 3 declarations associated with /MXINFO/ character*8 hbpar character*80 hmfile,hrule character*255 hbin dimension a(nmxpar) dimension ctij(nmxtc),cvij(nmxvc) c common /NCOMP/ ncomp common /CHAR/ htab,hnull common /CCAS/ hcas(n0:nx) c commons associated with the binary mixing rule(s) c the amix(i,j,k,1..3) are the coefs for the (general) mixing term c the fmix(i,j,1..nmxpar) are the parameters for the i-j binary common /MIXMOD/ hmodmx(nx,nx) common /CFXHMX/ amix(nx,nx,nmxtrm,3),fmix(nx,nx,nmxpar) c commons associated with the binary critical lines common /CRTMOD/ hmodtc(nx,nx),hmodvc(nx,nx) common /CFXCRT/ ftcij(nx,nx,nmxtc),fvcij(nx,nx,nmxvc) common /CFICRT/ nftcij(nx,nx),nfvcij(nx,nx) c /MXINFO/ contains information on the mixing rules and parameters c hbin provides documentation for the current binary parameters common /MXINFO/ hmfile(0:nbin),hbin(nbin),hrule(nbrule), & hbpar(nbrule,nmxpar) c ierr=0 herr=hnull i=icomp j=jcomp hcasi=hcas(icomp) hcasj=hcas(jcomp) * write (*,1010) icomp,jcomp,hcasi,hcasj *1010 format (1x,' RDBNC: input i,j,cas: ',2i3,2x,a12,1x,a12) c c initialize arrays storing coeff for critical lines do 60 k=1,nmxtc ftcij(i,j,k)=0.0d0 60 continue do 64 k=1,nmxvc fvcij(i,j,k)=0.0d0 64 continue c if (nread.le.0) then c get coefficients from common blocks--not currently implemented ierr=1 herr=' RDBNC ERROR--block data option not available'//hnull call ERRMSG (ierr,herr) RETURN else c c rewind file and read coefficients c rewind (nread) c c search for key characters in cols 1 and 2-4 and branch accordingly c hstar = '#' or '@' indicates start of model specification c hflag indicates type of model: c 'BNC' = binary mixture coefficients c 'MXM' = mixture model specification c do 180 iline=1,10000 read (nread,2013,end=199,err=199) hstar,hflag if (hstar.eq.'@' .and. (hflag.eq.'END'.or.hflag.eq.'end')) then ierr=-117 herr=' RDBNC--binary pair not found (end of file)'//hnull call ERRMSG (ierr,herr) RETURN end if if (hstar.eq.'#' .or. hstar.eq.'@') then read (nread,2003) htype do 120 k=1,100 !skip over source comments read (nread,2001) h1 if (h1.ne.'?') goto 124 !comment block terminated by '!' 120 continue 124 continue end if if (hflag.eq.'BNC' .and. htype.eq.'BNC') then do 176 jline=1,8000 read (nread,2256) h1,hbinp !read just first comment line if (h1.ne.'?') then c possible end of binary mixture data; start search for another #BNC keyword goto 180 end if read (nread,2001) h1 !skip over second comment line read (nread,2020,err=199) hcas1,hcas2,hmodij,(a(k),k=1,nmxpar) c write (*,*) ' RDBNC: hmodij,a(1): ',hmodij,a(1) read (nread,2040,err=199) hmodt,(ctij(k),k=1,nctc) read (nread,2040,err=199) hmodv,(cvij(k),k=1,ncvc) read (nread,2001) h1 !skip over separator c write (*,2040) hmodv,(cvij(k),k=1,nmxvc) c write (*,1014) hcas1,hcas2,hbinp c1015 format (1x,' RDBNC--found ',a12,'+ ',a12,2x,a255) if (hcasi.eq.hcas1 .and. hcasj.eq.hcas2) then * write (*,1015) hcas1,hcas2,hbinp *1015 format (1x,' RDBNC--found ',a12,'+ ',a12,2x,a255) hbin(ibin)=hbinp hmodmx(i,j)=hmodij hmodmx(j,i)=hmodij !pair j,i is same as i,j do 160 k=1,nmxpar fmix(i,j,k)=a(k) fmix(j,i,k)=a(k) 160 continue * write (*,1016) icomp,jcomp,hmodij,(a(k),k=1,nmxpar) *1016 format (1x,' RDBNC--i,j,mod,fij: ',9x,2i3,2x,a3,8f8.4) c write coefficients for critical temperature line to common hmodtc(i,j)=hmodt hmodtc(j,i)=hmodt !pair j,i is same as i,j if (hmodt.eq.'TC1') then c for model TC1, the order i,j vs j,i is significant; first coeff is c symetric, elements 2-6 store coeff for order i,j; 7-11 for j,i ftcij(i,j,1)=ctij(1) ftcij(j,i,1)=ctij(1) do 162 k=2,6 ftcij(i,j,k)=ctij(k) ftcij(j,i,k+5)=ctij(k) 162 continue end if c same for critical volume line hmodvc(i,j)=hmodv hmodvc(j,i)=hmodv !pair j,i is same as i,j if (hmodv.eq.'VC1') then c for model VC1, the order i,j vs j,i is significant; first coeff is c symetric, elements 2-6 store coeff for order i,j; 7-11 for j,i fvcij(i,j,1)=cvij(1) fvcij(j,i,1)=cvij(1) do 164 k=2,6 fvcij(i,j,k)=cvij(k) fvcij(j,i,k+5)=cvij(k) 164 continue end if RETURN else if (hcasi.eq.hcas2 .and. hcasj.eq.hcas1) then * write (*,1015) hcas1,hcas2,hbinp hbin(ibin)=hbinp hmodmx(i,j)=hmodij hmodmx(j,i)=hmodij !pair j,i is same as i,j do 170 k=1,nmxpar fmix(i,j,k)=a(k) fmix(j,i,k)=a(k) 170 continue c pair j,i is same as i,j except for beta and gamma parameters in c the LJi and LMi models are associated with a particular fluid if (hmodij(1:2).eq.'LM' .or. hmodij(1:2).eq.'LJ') then fmix(j,i,4)=a(5) fmix(j,i,5)=a(4) end if * write (*,1016) icomp,jcomp,hmodij,(a(k),k=1,nmxpar) c write coefficients for critical temperature line to common hmodtc(i,j)=hmodt hmodtc(j,i)=hmodt !pair j,i is same as i,j if (hmodt.eq.'TC1') then c for model TC1, the order i,j vs j,i is significant; first coeff is c symetric, elements 2-6 store coeff for order i,j; 7-11 for j,i ftcij(i,j,1)=ctij(1) ftcij(j,i,1)=ctij(1) do 172 k=2,6 ftcij(i,j,k+5)=ctij(k) ftcij(j,i,k)=ctij(k) 172 continue end if c same for critical volume line hmodvc(i,j)=hmodv hmodvc(j,i)=hmodv !pair j,i is same as i,j if (hmodv.eq.'VC1') then c for model VC1, the order i,j vs j,i is significant; first coeff is c symetric, elements 2-6 store coeff for order i,j; 7-11 for j,i fvcij(i,j,1)=cvij(1) fvcij(j,i,1)=cvij(1) do 174 k=2,6 fvcij(i,j,k)=cvij(k) fvcij(j,i,k+5)=cvij(k) 174 continue end if c check for specification that critical line(s) are to be estimated, c if so, set warning (estimation routine will be called from SETHMX) ierr1=0 if (hmodt.eq.'EST') then ierr1=ierr1-1 end if if (hmodv.eq.'EST') then ierr1=ierr1-2 end if if (ierr1.lt.0) then ierr=-113+ierr1 end if RETURN end if 176 continue end if 180 continue end if RETURN c 199 continue !abnornal termination * write (*,*) ' RDBNC line 199, hstar,hflag: ',hstar,hflag ierr=118 herr='[RDBNC error 118] error in reading mixture file'//hnull call ERRMSG (ierr,herr) RETURN c 2001 format (a1) 2003 format (a3) 2013 format (a1,a3) 2020 format (a12,1x,a12,1x,a3,1x,6f12.0) 2040 format (26x,a3,1x,6f12.0) 2256 format (a1,a255) c end !subroutine RDBNC c c ====================================================================== c subroutine RDMIX (nread,icomp,jcomp,hmodij,lij,ierr,herr) c c read coefficients for mixing rule(s) from file 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 by subroutine SETHMX) c icomp--component i; ignored if hmodij='ALL' c jcomp--component j; ignored if hmodij='ALL' c hmodij--model for binary interaction of i and j [character*3] c 'ALL' indicates read model(s) for all binary pair(s) c (used when this routine is called from SETHMX) c when a specific model is specified (e.g. 'LJ1' or 'MKT') c read only that model (used when called from SETKTV) c outputs: c lij--flags indicating if mixing term found for binary pair (i,j) c [logical array of dimension (ncomp,ncomp)]; used only in c conjunction with hmodij='ALL' c ierr--error flag: 0 = successful c 118 = error in reading mixture file c herr--error string (character*255 variable if ierr<>0) c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 05-07-96 MM, original version (extracted from SETHMX) c 05-08-96 MM, add /MXINFO/ and read corresponding info from file c 05-10-96 MM, read default binary parameters (used only by interface) c 11-04-96 MM, change nmxpar from 4 to 6 c 12-02-97 MM, add /MXCFID/ to store ideal-solution values c implicit double precision (a-h,o-z) implicit integer (i-k,m,n) implicit logical (l) c parameter (ncmax=5) !max number of components in mixture c parameter (n0=-ncmax,nx=ncmax) parameter (nx=ncmax) parameter (nbrule=10) !number of binary mixing rules parameter (nbin=ncmax*(ncmax-1)/2) !# possible binary pairs parameter (nmxpar=6) !number of binary mixture parameters parameter (nmxtrm=15) !number of terms in binary mixing func character*1 htab,hnull character*1 hstar,h1 character*3 hmodmx,hmodij character*3 hflag,htype character*8 hpar(nmxpar) character*76 hinfo character*255 herr c next 3 declarations associated with /MXINFO/ character*8 hbpar character*80 hmfile,hrule character*255 hbin c dimension iexp(nmxtrm) dimension aterm(nmxtrm,3) logical lij(nx,nx) c common /NCOMP/ ncomp common /CHAR/ htab,hnull c commons associated with the binary mixing rule(s) c the amix(i,j,k,1..3) are the coefs for the (general) mixing term c the fmix(i,j,1..nmxpar) are the parameters for the i-j binary c the namix(i,j) are the number of a(i,j,k,1..3) terms c the iexpmx(i,j,k) are powers of the exponentials in the general term common /MIXMOD/ hmodmx(nx,nx) common /CFXHMX/ amix(nx,nx,nmxtrm,3),fmix(nx,nx,nmxpar) common /CFIHMX/ namix(nx,nx),iexpmx(nx,nx,nmxtrm) c /MXINFO/ contains information on the mixing rules and parameters c hmfile specifies the files from which mixing rules originate c (stored in order 1,2 1,3 2,3; element zero is file called in SETUP) c hbin provides documentation for the current binary parameters c hrule contains descriptions of the currently available mixing rules c hbpar contains descriptions of the binary parameters (e.g. Kt, Kv) c associated with the currently available mixing rules common /MXINFO/ hmfile(0:nbin),hbin(nbin),hrule(nbrule), & hbpar(nbrule,nmxpar) common /MXRULE/ nrule c default (ideal mixing) values for the mixing rule coefficients common /MXCFID/ fideal(nbrule,nmxpar) c c search for key characters in cols 1 and 2-4 and branch accordingly c hstar = '#' or '@' indicates start of model specification c hflag indicates type of model: c 'BNC' = binary mixture coefficients c 'MXM' = mixture model specification c * write (*,*) if (hmodij.eq.'ALL') then c initialize flags indicating if mixing terms have been found irule=0 do 312 i=1,ncomp-1 do 310 j=i+1,ncomp lij(i,j)=.false. 310 continue 312 continue else do 316 k=1,nrule if (hmodij.eq.hrule(k)(1:3)) then irule=k goto 318 end if 316 continue irule=nrule+1 * write (*,*) ' RDHMX--new mixing rule found',hmodij 318 continue end if c rewind file (RDBNC may have left it at the end) rewind (nread) do 380 iline=1,2000 read (nread,2013,end=998,err=998) hstar,hflag if ((hstar.eq.'#' .or. hstar.eq.'@').and. hflag.eq.'MXM') then if (hmodij.eq.'ALL') then irule=irule+1 end if read (nread,2053) htype,hinfo do 320 k=1,100 !skip over source comments read (nread,2001) h1 if (h1.ne.'?') goto 324 !comment block terminated by '!' 320 continue 324 continue c read descriptors for binary-specific parameters read (nread,2048) (hpar(k),k=1,nmxpar) c read default (ideal-solution) parameters read (nread,*) (fideal(irule,k),k=1,nmxpar) c read coefficients for mixing rule read (nread,*,err=998) nterm,nxcoef,nicoef c nterm is number of terms c nxcoef is number of real coefficients per term c nicoef is number of integer coefficients per term (only 0,1 allowed) do 340 k=1,nterm if (nicoef.le.0) then read (nread,*) (aterm(k,n),n=1,nxcoef) else read (nread,*) (aterm(k,n),n=1,nxcoef),iexp(k) end if 340 continue if (hmodij.eq.'ALL') then ipair=0 !counter for binary pair c write (*,1342) irule,htype,hinfo c1342 format (1x,' SETHMX--mixing function',i3,' (',a3,'): ',a76) hrule(irule)=htype//hinfo do 342 k=1,nmxpar hbpar(irule,k)=hpar(k) 342 continue do 362 i=1,ncomp-1 do 360 j=i+1,ncomp if (htype .eq. hmodmx(i,j)) then c mixing term found; load coefficients into working arrays lij(i,j)=.true. ipair=ipair+1 * write (*,1344) irule,htype,i,j,hrule(irule) *1344 format (1x,' SETHMX: use mixing rule #',i3,' (',a3, * & ') for i,j =',2i3,': ',a80) namix(i,j)=nterm do 350 k=1,nterm do 348 n=1,nxcoef amix(i,j,k,n)=aterm(k,n) 348 continue if (nicoef.ge.1) then iexpmx(i,j,k)=iexp(k) end if 350 continue end if 360 continue 362 continue nrule=irule else if (htype.eq.hmodij) then c mixing term found; load coefficients into working arrays lij(icomp,jcomp)=.true. * write (*,1370) htype,icomp,jcomp *1370 format (1x,' SETKTV: use ',a3,' mixing term for i,j =',2i3) * write (*,*) ' SETHMX--mixing term irule = ',irule,htype,hinfo hrule(irule)=htype//hinfo do 364 k=1,nmxpar hbpar(irule,k)=hpar(k) 364 continue namix(icomp,jcomp)=nterm namix(jcomp,icomp)=nterm !pair (j,i) is same as (i,j) do 370 k=1,nterm do 368 n=1,nxcoef amix(icomp,jcomp,k,n)=aterm(k,n) amix(jcomp,icomp,k,n)=aterm(k,n) !pair (j,i) is same as (i,j) 368 continue if (nicoef.ge.1) then iexpmx(icomp,jcomp,k)=iexp(k) iexpmx(jcomp,icomp,k)=iexp(k) !pair (j,i) is same as (i,j) end if 370 continue end if else if (hflag.eq.'END') then * write (*,*) ' RDMIX end of file encountered' RETURN end if 380 continue RETURN c c error terminations c 998 continue * write (*,*) ' SETHMX line 998, hstar,hflag: ',hstar,hflag ierr=118 herr='[RDMIX error 118] error in reading mixture file'//hnull call ERRMSG (ierr,herr) RETURN c 2001 format (a1) 2003 format (a3) 2013 format (a1,a3) c multiplier in format should match "parameter (nmxpar=6)" 2048 format (6(a8,1x)) !descriptors for binary-specific parameters 2053 format (a3,a76) c end !subroutine RDMIX c c ====================================================================== c subroutine ESTBNC (i,j,ibin,ierr,herr) c c estimate binary parameters for the HMX mixture model c c inputs: c i--component i c j--component j c ibin--counter for binary pair (1,2 = 1; 1,3 = 2; 2,3 = 3; etc.) c outputs: c ierr--error flag: 0 = successful c 117 = no binary parameters found c herr--error string (character*255 variable if ierr<>0) c other outputs written to common c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 01-12-96 MM, original version c 02-27-96 MM, parameter n0=-ncmax to accomodate ECS-thermo model c add Zcrit to common /CCON/ c 03-19-19 MM, add dipole moment to /CCON/ c 05-06-96 MM, change default model to VDW c 05-08-96 MM, add /MXINFO/ and load corresponding info c 11-04-96 MM, replace van der Waals rule with linear (LIN), c change nmxpar from 4 to 6 c 11-19-96 MM, move initialization of common blocks to RDBNC, ESTBNC c change argument list c 07-03-97 MM, return ierr=117 if this routine is called (causes error c alert to appear in GUI); !temporary until prediction in place c 09-08-97 MM, change error 117 to warning--GUI displayed alert, then quit c 11-06-97 MM, implement preliminary predictive model of Lemmon c 01-09-98 MM, warning if estimated zeta > 100 c implicit double precision (a-h,o-z) implicit integer (i-k,m,n) implicit logical (l) c parameter (ncmax=5) !max number of components in mixture parameter (n0=-ncmax,nx=ncmax) parameter (nmxpar=6) !number of binary mixture parameters parameter (nbrule=10) !number of binary mixing rules parameter (nbin=ncmax*(ncmax-1)/2) !# possible binary pairs parameter (nmxtrm=15) !number of terms in binary mixing rule character*1 htab,hnull character*3 hmodmx character*12 hcas character*255 herr c next 3 declarations associated with /MXINFO/ character*8 hbpar character*80 hmfile,hrule character*255 hbin logical lest c common /CHAR/ htab,hnull common /CCAS/ hcas(n0:nx) common /CCON/ wm(n0:nx),ttp(n0:nx),tnbp(n0:nx), & tc(n0:nx),pc(n0:nx),rhoc(n0:nx),Zcrit(n0:nx), & accen(n0:nx),dipole(n0:nx) c commons associated with the binary mixing rule(s) c the amix(i,j,k,1..3) are the coefs for the (general) mixing term c the fmix(i,j,1..nmxpar) are the parameters for the i-j binary common /MIXMOD/ hmodmx(nx,nx) common /CFXHMX/ amix(nx,nx,nmxtrm,3),fmix(nx,nx,nmxpar) c /MXINFO/ contains information on the mixing rules and parameters c hbin provides documentation for the current binary parameters common /MXINFO/ hmfile(0:nbin),hbin(nbin),hrule(nbrule), & hbpar(nbrule,nmxpar) c * write (*,1002) i,j,hcas(i),hcas(j) *1002 format (1x,' ESTBNC will estimate i,j = ',2i3,2x,a12,1x,a12) lest=.true. !temporary--enable/disable estimation scheme if (.not.lest) then c return default values for now hmodmx(i,j)='LIN' hmodmx(j,i)='LIN' hbin(ibin)='No mixture data are available for this binary pair.' & //' Ideal solution behavior will be assumed.' & //hnull do 100 k=1,2 fmix(i,j,k)=1.0d0 !fmix(i,j,1) = Kt fmix(j,i,k)=1.0d0 !fmix(i,j,2) = Kv 100 continue fmix(i,j,3)=0.0d0 !Fpq(i,j) fmix(j,i,3)=0.0d0 !Fpq(k,i) do 120 k=4,5 fmix(i,j,k)=1.0d0 !fmix(i,j,1) = beta fmix(j,i,k)=1.0d0 !fmix(i,j,2) = gamma 120 continue do 140 k=6,nmxpar fmix(i,j,k)=0.0d0 !not used fmix(j,i,k)=0.0d0 140 continue c return warning message to GUI ierr=-117 herr='[SETUP warning 117] No mixture data are available for one' & //' or more binary pairs in the specified mixture. Ideal' & //' solution behavior will be assumed.'//hnull call ERRMSG (ierr,herr) RETURN end if c interim predictive scheme of Lemmon hmodmx(i,j)='LJ1' hmodmx(j,i)='LJ1' hbin(ibin)='No mixture data are available for this binary pair. ' & //'The mixing parameters have been estimated.' & //hnull c ordering of the components depends on ratio of the dipole moments if (ABS(dipole(i)-dipole(j)).lt.1.0d-5) then c special case for equal dipole moments if (tc(i)/(pc(i)*accen(i)).gt.tc(j)/(pc(j)*accen(j))) then tratio=tc(i)/tc(j) pratio=pc(i)/pc(j) accenr=accen(i)/accen(j) else tratio=tc(j)/tc(i) pratio=pc(j)/pc(i) accenr=accen(j)/accen(i) end if else if (dipole(i).lt.dipole(j)) then tratio=tc(i)/tc(j) pratio=pc(i)/pc(j) accenr=accen(i)/accen(j) else tratio=tc(j)/tc(i) pratio=pc(j)/pc(i) accenr=accen(j)/accen(i) end if c define the zeta parameter (the one applied to the reducing temperature) fmix(i,j,1)=(40.4d0-25.03d0*2.0d0**(tratio/pratio/accenr))/tratio fmix(j,i,1)=fmix(i,j,1) c set remaining parameters to defaults; 2 = xi; 3 = Fpq; 4 = beta; 5 = gamma do 200 k=2,3 fmix(i,j,k)=0.0d0 !k = 2: xi (term in reducing volume) fmix(j,i,k)=0.0d0 !k = 3: Fpq (multiplier for excess function) 200 continue do 202 k=4,5 fmix(i,j,k)=1.0d0 !k = 4: beta (exponent for first component) fmix(j,i,k)=1.0d0 !k = 5: gamma (exponent for second component) 202 continue do 204 k=6,nmxpar fmix(i,j,k)=0.0d0 !not used fmix(j,i,k)=0.0d0 204 continue c return warning message to GUI if (ABS(fmix(i,j,1)).gt.100.0d0) then c issue warning if zeta greater than +/- 100 ierr=110 herr='[SETUP warning 110] The specified mixture is at the limit' & //' of applicability of the model--calculations may fail or' & //' give unreliable results. Check results carefully.' & //hnull c constrain zeta to a maximum of +/- 200 fmix(i,j,1)=SIGN(1.0d0,fmix(i,j,1))*MIN(ABS(fmix(i,j,1)),2.0d2) fmix(j,i,1)=fmix(i,j,1) else ierr=-117 herr='[SETUP warning 117] No mixture data are available for one' & //' or more binary pairs in the specified mixture. The' & //' mixing parameters have been estimated.'//hnull end if call ERRMSG (ierr,herr) RETURN c end !subroutine ESTBNC c c ====================================================================== c subroutine ESTTC (i,j,ierr,herr) c c estimate parameters for binary critical temperature line c c inputs: c i--component i c j--component j c outputs: c ierr--error flag: 0 = successful c 1 = error c herr--error string (character*255 variable if ierr<>0) c other outputs written to common c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 11-19-96 MM, original version 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 (nmxtc=11,nmxvc=11) !number of coeff for critical lines character*1 htab,hnull character*3 hmodtc,hmodvc character*12 hcas character*255 herr c common /CHAR/ htab,hnull common /CCAS/ hcas(n0:nx) common /CCON/ wm(n0:nx),ttp(n0:nx),tnbp(n0:nx), & tc(n0:nx),pc(n0:nx),rhoc(n0:nx),Zcrit(n0:nx), & accen(n0:nx),dipole(n0:nx) c commons associated with the binary critical lines common /CRTMOD/ hmodtc(nx,nx),hmodvc(nx,nx) common /CFXCRT/ ftcij(nx,nx,nmxtc),fvcij(nx,nx,nmxvc) common /CFICRT/ nftcij(nx,nx),nfvcij(nx,nx) c * write (*,1002) i,j,hcas(i),hcas(j) *1002 format (1x,' ESTTC will estimate i,j = ',2i3,2x,a12,1x,a12) c return default values for now hmodtc(i,j)='TC1' hmodtc(j,i)='TC1' do 200 k=1,nmxtc ftcij(i,j,k)=0.0d0 ftcij(j,i,k)=0.0d0 200 continue ierr=0 herr=hnull RETURN c end !subroutine ESTTC c c ====================================================================== c subroutine ESTVC (i,j,ierr,herr) c c estimate parameters for binary critical volume line c c inputs: c i--component i c j--component j c outputs: c ierr--error flag: 0 = successful c 1 = error c herr--error string (character*255 variable if ierr<>0) c other outputs written to common c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 11-19-96 MM, original version 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 (nmxtc=11,nmxvc=11) !number of coeff for critical lines character*1 htab,hnull character*3 hmodtc,hmodvc character*12 hcas character*255 herr c common /CHAR/ htab,hnull common /CCAS/ hcas(n0:nx) common /CCON/ wm(n0:nx),ttp(n0:nx),tnbp(n0:nx), & tc(n0:nx),pc(n0:nx),rhoc(n0:nx),Zcrit(n0:nx), & accen(n0:nx),dipole(n0:nx) c commons associated with the binary critical lines common /CRTMOD/ hmodtc(nx,nx),hmodvc(nx,nx) common /CFXCRT/ ftcij(nx,nx,nmxtc),fvcij(nx,nx,nmxvc) common /CFICRT/ nftcij(nx,nx),nfvcij(nx,nx) c * write (*,1002) i,j,hcas(i),hcas(j) *1002 format (1x,' ESTVC will estimate i,j = ',2i3,2x,a12,1x,a12) c return default values for now hmodvc(i,j)='VC1' hmodvc(j,i)='VC1' do 200 k=1,nmxvc fvcij(i,j,k)=0.0d0 fvcij(j,i,k)=0.0d0 200 continue ierr=0 herr=hnull RETURN c end !subroutine ESTVC c c ====================================================================== c block data BDHMX c c initialize common blocks which save the phi, tau, and del from one c call of PHIMIX to the next c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 06-12-96 MM, original version 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 (nmxpar=6) !number of binary mixture parameters parameter (nmxtrm=15) !number of terms in binary mixing rule parameter (nxnx=nx*nx) parameter (nx2trm=nx*nx*(nmxtrm+1)) common /HMXSAV/ phisav(nx,nx,0:nmxtrm),tausav(nx,nx),delsav(nx,nx) c data phisav /nx2trm*0.0d0/ data delsav /nxnx*0.0d0/ data tausav /nxnx*0.0d0/ c end !block data BDHMX c c c 1 2 3 4 5 6 7 c23456789012345678901234567890123456789012345678901234567890123456789012 c c ====================================================================== c end file mix_HMX.f c ======================================================================