c begin file trns_ECS.f c c This file contains the routines implementing the extended corresponding c states (ECS) method for the transport properties. c c contained here are: c subroutine TRNECS (t,rho,x,eta,tcx,ierr,herr) c subroutine TCBKMX (t,rho,x,fj,fx,hj,hx,tcx,Flam,lerrt,lerrD, c & ierr,herr) c subroutine SETTRN (nread,icomp,hcasno,href,heqn,hvs,htc,ierr,herr) c function PSI (icomp,tr,rhor) c function CHI (icomp,tr,rhor) c subroutine ECSLIM (t,D,tmin,tmax,Dmax,lerrt,lerrD,terr,Derr) c function ETA0DG (icomp,t) c function OMEGAS (il,is,tau) c FUNCTION DELHSV (TX,DX,X,hj) c SUBROUTINE ENSKOG (N,RHO,SIGMA,CMW,X,ETA) c function ETAMIX (t,x) c subroutine CONFTD (j,amix,Zmix,tj,rhoj,ierr,herr) c subroutine CONFD (j,amix,Zmix,tj,rhoj,ierr,herr) c subroutine CONFT (k,amix,rhok,tk,ierr,herr) c function CBRT (x) c c ====================================================================== c ====================================================================== c subroutine TRNECS (t,rho,x,eta,tcx,ierr,herr) c c compute the transport properties of thermal conductivity and c viscosity as functions of temperature, density, and composition c c based on the modification of the Huber-Ely ECS method given by: c Klein, S.A., McLinden, M.O. and Laesecke, A. (1997). An improved c extended corresponding states method for estimation of viscosity of c pure refrigerants and mixtures. Int. J. Refrigeration 20:208-217 c c N.B. --equation numbers below refer to this paper; c --factor of (ref fluid mol wt)**-0.5 is missing from Eq 31 in paper c --in Eq 34 reference fluid should be evaluated at (t0,rho0), not c at (t/fj,rho*hj) c c inputs: c t--temperature [K] c rho--molar density [mol/L] c x--composition array [mol frac] c outputs: c eta--viscosity [uPa.s] c tcx--thermal conductivity [W/m.K] c ierr--error flag: 0 = successful c -35 = temperature out of range for conductivity c -36 = density out of range for conductivity c -37 = T and D out of range for conductivity c -45 = temperature out of range for viscosity c -46 = density out of range for viscosity c -47 = T and D out of range for viscosity c -55 = T out of range for both visc and t.c. c -56 = D out of range for both visc and t.c. c -57 = T and/or D out of range for both visc and t.c. c -58,-59 = ECS model did not converge c herr--error string (character*255 variable if ierr<>0) c c written by M. McLinden c NIST Physical & Chemical Properties Division, Boulder, Colorado c based on the routine by S.A. Klein (in turn, based on Refprop5 routine) c 03-05-97 MM, original version c 08-22-97 MM, evaluate critical part of t.c. at simple reduced t,rho rather c than the conformal t,rho used for the background part c 08-25-97 MM, replace calls to ETA1 with ETAK0 (ref fluid) and ETA0DG c 08-26-97 MM, break-out mix t.c. into separate subroutine TCBKMX c 09-09-97 MM, calls to old GETFXHX routine, commented out c if GETFH does not converge, set to xnotc and RETURN c 09-25-97 MM, restructure around new CONFTD for finding f's, h's c 10-01-97 MM, make correction to Eq 34 as noted above c 10-08-97 MM, evaluate crit part of t.c. at average of reduced and conformal c 10-09-97 MM, fix bug in do loop associated with calc of xmij c 10-24-97 MM, Eucken term f_int now a function of t c 11-13-97 MM, fix bomb when rho = 0 c implicit double precision (a-h,o-z) implicit integer (i-n) parameter (ncmax=5) !max number of components in mixture parameter (n0=-ncmax,nx=ncmax) parameter (nrf0=0) !lower limit for transport ref fluid arrays parameter (mxtrn=10) !max no. coefficients for psi, chi function character*1 htab,hnull character*3 hetamx,heta,htcxmx,htcx character*255 herr,herr2 logical lerrtt,lerrvt,lerrtD,lerrvD,lerrtc,lerrvs dimension x(ncmax) dimension fj(ncmax),hj(ncmax) !reducing ratios for components dimension xmj(ncmax) !equivalent mol mass (Eq 33) common /CHAR/ htab,hnull common /TRNMOD/ hetamx,heta(nrf0:ncmax),htcxmx,htcx(nrf0:ncmax) c limits common /WLMTRN/ tmin(nrf0:nx),tmax(nrf0:nx),pmax(nrf0:nx), & rhomax(nrf0:nx) c numbers of terms for the various parts of the model: c LJflag: flag for L-J parameters (if 0, estimate) c Euck: factor f_int in Eucken correlation c psi (viscosity shape factor): polynomial term, 2nd poly, spare c chi (conductivity shape factor): polynomial term, 2nd poly, spare common /WNTTRN/ LJflag(nrf0:nx),nEuck(nrf0:nx), & npsi1(nrf0:nx),npsi2(nrf0:nx),npsi3(nrf0:nx), & nchi1(nrf0:nx),nchi2(nrf0:nx),nchi3(nrf0:nx) c commons storing the (real and integer) coefficients to the ECS model common /WCFTRN/ cpsi(nrf0:nx,mxtrn,4),cchi(nrf0:nx,mxtrn,4) common /WIFTRN/ ipsi(nrf0:nx,0:mxtrn),ichi(nrf0:nx,0:mxtrn) c Lennard-Jones parameters common /WLJTRN/ sigma(nrf0:nx),epsk(nrf0:nx) c number of components in mix and constants for the mix components common /NCOMP/ nc 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 /Gcnst/ R c common block containing flags to GUI (initialized in BDSET in setup.f) common /FLAGS/ xnota,x2ph,xsubc,xsuph,xsupc,xinf,x7,xnotd,xnotc c ierr=0 herr=hnull c set visc and t.c. to flags indicating 'not calculated' so that some c value is returned to GUI in event of failure of routines eta=xnotc tcx=xnotc c c find amix, Zmix and conformal t,rho for reference fluid call CRITP (x,tcmx,pcmx,rhocmx,ierr,herr) c write (*,*) ' TRNECS--ierr from CRITP, tcmx: ',ierr,tcmx call REDX (x,tred,Dred) tau=tred/t del=rho/Dred amix=PHIX(0,0,tau,del,x) Zmix=1.0d0+del*PHIX(0,1,tau,del,x) t0=t*tc(0)/tcmx !initial guess for conformal temperature rho0=rho*rhoc(0)/rhocmx !initial guess for conformal density if (rho.gt.1.0d-4) then c find "exact" conformal t,rho only if density is significant; c at zero density, use the initial guesses above (CONFTD can fail at c very low density, and the dilute-gas contribution is dominant anyway) call CONFTD (0,amix,Zmix,t0,rho0,ierr,herr) fx=t/t0 !reducing ratios for mixture hx=rho0/rho else fx=tcmx/tc(0) hx=rhoc(0)/rhocmx end if c if (nc.eq.1) then c for pure fluid, component f,h are same as mixture values fj(1)=fx hj(1)=hx else c for mixture, find conformal t,rho for each of the components do 100 j=1,nc c tj=t*tc(j)/(fx*tc(0)) !initial guess for conformal temp c rhoj=rho*hx*rhoc(j)/rhoc(0) !initial guess for conformal density tj=t*tc(j)/tcmx !initial guess for conformal temp c write (*,1097) j,t,tc(j),tcmx,tj c1097 format (1x,' TRNECS--j,t,tc(j),tcmx,tj: ',i3,4f12.6) rhoj=rho*rhoc(j)/rhocmx !initial guess for conformal density if (rho.gt.1.0d-4) then c find "exact" conformal t,rho only if density is significant; c at zero density, use the initial guesses above (CONFTD can fail at c very low density, and the dilute-gas contribution is dominant anyway) call CONFTD (j,amix,Zmix,tj,rhoj,ierr,herr) c write (*,*) ' TRNECS--ierr from CONFTD, tj: ',ierr,tj fj(j)=tj*fx/t hj(j)=rho*hx/rhoj c write (*,1098) j,tj,fx,t,tcmx c1098 format (1x,' TRNECS--j,tj,fx,t,tcmx: ',i3,4f12.6) else fj(j)=tc(j)*fx/tcmx hj(j)=hx*rhocmx/rhoc(j) end if c write (*,1099) j,fj(j),hj(j) c1099 format (1x,' TRNECS--j,fj,hj:',i2,2f12.6) 100 continue end if c c estimate the Lennard-Jones parameters, if necessary do 120 j=1,nc if (LJflag(j).eq.0) then epsk(j)=tc(j)/tc(0)*epsk(0) sigma(j)=sigma(0)*(rhoc(0)/rhoc(j))**(1.0d0/3.0d0) else if (LJflag(j).eq.2) then c estimation method of Huber & Ely (1992) FPE 80:239-248 epsk(j)=epsk(0)*fj(j) sigma(j)=sigma(0)*hj(j)**(1.0d0/3.0d0) end if 120 continue c c find correlation limits for reference fluid p=0.0d0 call LIMITK ('ETA',0,t,rho,p,tminv,tmaxv,Deta,peta,ierr2,herr2) call LIMITK ('TCX',0,t,rho,p,tmint,tmaxt,Dtcx,ptcx,ierr2,herr2) c initialize error flags lerrtt=.false. lerrvt=.false. lerrtD=.false. lerrvD=.false. c write (*,1013) tminv,tmaxv,Deta,tmint,tmaxt,Dtcx c1013 format (1x,' TRNECS--ref fluid visc limits--t,rho: ',2f8.2,f12.6/ c & 1x,' t.c. limits--t,rho: ',2f8.2,f12.6) c if (nc.eq.1) then c c special case for pure fluid c gx=wm(1)/wm(0) tpsi=t/tc(1) rhopsi=rho/rhoc(1) rho0v=rho*hj(1)*PSI(1,tpsi,rhopsi) !Eq 21 c reference fluid background viscosity call ETAKB (0,t0,rho0v,eta0bk,ierr,herr) eta1dg=ETA0DG(1,t) !dilute gas visc Feta=SQRT(fj(1)*gx)*hj(1)**(-2.0d0/3.0d0) !Eq 11 eta=eta1dg+eta0bk*Feta !Eqs 5,10 c check conformal t,rho against limits of reference fluid correlation call ECSLIM (t0,rho0v,tminv,tmaxv,Deta,lerrvt,lerrvD,tcf,Dcf) c similar terms for thermal conductivity tchi=tpsi rhochi=rhopsi rho0t=rho*hj(1)*CHI(1,tchi,rhochi) c note that F-factor for t.c. has inverse power of gx compared to visc Flam=SQRT(fj(1)/gx)*hj(1)**(-2.0d0/3.0d0) call TCXKB (0,t0,rho0t,tcx0b,ierr,herr) c find dilute-gas parts from collisions and internal degrees of freedom tcx1dg=1.0d-3*15.0d0*R*eta1dg/(4.0d0*wm(1)) c tcx1in=1.32d-3*eta1dg/wm(1)*(CP0K(1,t)-2.5d0*R) tcx1in=FINT(1,t)*eta1dg/wm(1)*(CP0K(1,t)-2.5d0*R) c apply ECS method to background part of ref fluid tcx=tcx0b*Flam+tcx1dg+tcx1in c check conformal t,rho against limits of reference fluid correlation call ECSLIM (t0,rho0t,tmint,tmaxt,Dtcx,lerrtt,lerrtD,tcf,Dcf) c c write (*,1015) t,1,tcx1dg+tcx1in,tcx0b*Flam,Flam,tcx c1015 format (' TRNECS--t,j = ',f8.2,i3, c & '; tcx_dg,tcx_bk,Flam,tcx: ',4f10.6) c write (*,1017) PSI(1,tpsi,rhopsi),Feta,CHI(1,tchi,rhochi),Flam c1017 format (1x,' TRNECS--psi,Feta,chi,Flam: ',4f14.8) c else c c general (mixture) case; begin mixture viscosity calculation c c calculate "equivalent mass" gx, using pure fluid residual viscosities, c either from a pure fluid correlation or the pure fluid ECS method; c check conformal t,rho against limits of reference fluid correlation c and find reference fluid residual viscosity call ECSLIM (t0,rho0,tminv,tmaxv,Deta,lerrvt,lerrvD,tcf,Dcf) call ETAKB (0,t0,rho0,eta0bk,ierr,herr) c write (*,1022) t0,rho0,eta0bk c1022 format (1x,' TRNECS-- ref fluid eval at t,rho =', c & f9.3,f12.6,' eta0bk = ',f12.6) gxsum=0.0d0 do 160 j=1,nc ierr2=0 if (heta(j)(1:2).ne.'EC') then c a pure fluid correlation is available c write (*,1024) j,heta(j) c1024 format (1x,' TRNECS--pure fluid corr for j = ',i3,': ',a3) tj=t*fj(j)/fx !Eq 36 rhoj=rho*hx/hj(j) !Eq 37 c check that pure fluid correlation is within its limits call LIMITK ('ETA',j,tj,rhoj,p,tmn,tmx,Dmx,pmx,ierr2,herr2) if (ierr2.eq.0) then call ETAKB (j,tj,rhoj,etaj,ierr1,herr) end if end if if (heta(j)(1:2).eq.'EC' .or. ierr2.ne.0) then c must use ECS method to estimate c write (*,1026) j c1026 format (1x,' TRNECS--will use ECS method for visc, j = ',i3) tpsi=t*fj(j)/(fx*tc(j)) rhopsi=rho*hx/(hj(j)*rhoc(j)) rho0j=rho*hx*PSI(j,tpsi,rhopsi) !Eq 39 c check conformal t,rho against limits of reference fluid correlation call ECSLIM (t0,rho0j,tminv,tmaxv,Deta,lerrvt,lerrvD,tcf,Dcf) call ETAKB (0,t0,rho0j,etarfj,ierr,herr) Fetaj=SQRT(fj(j)*wm(j)/wm(0))*hj(j)**(-2.0d0/3.0d0) !Eq 11 etaj=etarfj*Fetaj !residual viscosity of comp j end if if (eta0bk.gt.1.0d-6) then gj=etaj/eta0bk/SQRT(fj(j))*hj(j)**(2.0d0/3.0d0) !Eq 35 else c it is possible for residual viscosity to go through zero gj=wm(j)/wm(0)/SQRT(fj(j))*hj(j)**(2.0d0/3.0d0) end if xmj(j)=gj*gj*wm(0) !Eq 33 c write (*,*) ' TRNECS--j,xmj(j) for visc: ',j,xmj(j) c 160 continue do 180 j=1,nc do 176 i=1,nc xmij=2.0d0*xmj(i)*xmj(j)/(xmj(i)+xmj(j)) fij=SQRT(fj(j)*fj(i)) !Eq 26 hij=0.125d0*(hj(j)**(1.d0/3.d0)+hj(i)**(1.d0/3.d0))**3 !Eq 27 gxsum=gxsum+x(j)*x(i)*SQRT(fij*xmij)*hij**(4.0d0/3.0d0) !Eq 31 176 continue c write (*,1156) j,fj(j),hj(j),etaj c1156 format (1x,' TRNECS--j,fj,hj,etaj: ',i2,3f12.6) 180 continue c gx=gxsum/(SQRT(fx*wm(0))*hx**(4.0d0/3.0d0)) !Eq 31 c Feta=SQRT(fx)*hx**(-2.0d0/3.0d0)*gx !Eq 30 c Eq 30 + 31 reduce to following expression c factor of 1/SQRT(ref fluid mol wt) is missing from Klein paper Feta=gxsum/hx**2/SQRT(wm(0)) c write (*,1158) fx,hx,gxsum,Feta c1158 format (1x,' TRNECS--fx,hx,gxsum,Feta: ',4f12.6) etaxdg=ETAMIX(t,x) !dilute gas viscosity of mixture if (rho.gt.1.0d-6) then c apply size correction only if density is significant del=DELHSV(t,rho,x,hj) !Enskog size correction else del=0.0d0 end if eta=etaxdg+eta0bk*Feta+del !Eq 23 c write (*,1160) etaxdg,eta0bk,Feta,del,eta c1160 format (1x,' TRNECS--etaxdg,eta0bk,Feta,del,eta: ',5f12.6) c c mixture thermal conductivity calculation c evaluate background part (including dilute-gas part) of mixture c call TCBKMX (t,rho,x,fj,fx,hj,hx,tcx,Flam,lerrtt,lerrtD, & ierr,herr) c write (*,1164) tcx,Flam,lerrtt,lerrtD c1164 format (1x,' TRNECS--tcx,Flam,lerrtt,lerrtD out of TCBKMX: ', c & 2f12.6,2(3x,i1)) c end if !end mixture case c c now compute critical enhancement part of thermal conductivity c c critical enhancement part of t.c. is evaluated at simple reduced t,rho; c (if it were calculated at same conformal t,rho as the background part, c the enhancement would peak at something other than the critical point) c tr=t/tcmx*tc(0) c rhor=rho/rhocmx*rhoc(0) c c The critical enhancement part of t.c. is evaluated at a weighted c average of the conformal t,rho and the simple reduced t,rho, such that c it approaches the latter at the critical point and the former away c from the critical point. If it were calculated at same conformal t,rho c as the background part (even close to the critical point), the c enhancement would peak at something other than the critical point. c But, simply evaluating it always at the reduced t,rho can result in c states inside the two-phase region. Thus, the need for a compromise. tr=t/tcmx rhor=rho/rhocmx c write (*,*) ' TRNECS--tcmx,rhocmx: ',tcmx,rhocmx tslope=5.0d0 !revert to conformal at t > 1.2tr Dslope=2.5d0 qt=MIN(1.0d0,tslope*ABS(1.0d0-tr)) qD=MIN(1.0d0,Dslope*ABS(1.0d0-rhor)) tr=qt*t0+(1.0d0-qt)*tr*tc(0) rhor=qD*rho0+(1.0d0-qD)*rhor*rhoc(0) c call TCXKC (0,tr,rhor,tcx0c,ierr,herr) c write (*,1168) tr,rhor,tcx0c c1168 format (1x,' TRNECS--t,rho for ref fluid,tcx0c: ',3f12.6) c scale the critical enhancement by same Flam as rest of ECS method tcxcr=tcx0c*Flam c xcsum=0.0d0 tcrsum=0.0d0 do 260 j=1,nc if (htcx(j)(1:2).ne.'EC') then c a pure fluid correlation is available for component j; use for the c corresponding portion of the critical enhancement; otherwise mixture c calculation as x-->1 would not be continuous with pure fluid tr=t/tcmx*tc(j) rhor=rho/rhocmx*rhoc(j) call TCXKC (j,tr,rhor,tcxjc,ierr,herr) xcsum=xcsum+x(j) tcrsum=tcrsum+x(j)*tcxjc end if 260 continue c c now, finally, add the critical enhancement to the total t.c. tcrsum=tcrsum+(1.0d0-xcsum)*tcxcr tcx=tcx+tcrsum c write (*,1216) tcxbk0,tcbkcr,tcrsum c1216 format (1x,' TRNECS--bk_ref,bk_mix,crit:',3f12.8) c write (*,1217) tcxcr c1217 format (' TRNECS--tcx_crit:',2f10.6) c c process warnings/errors for conformal states outside range of ref fluid lerrtc=lerrtt.or.lerrtD lerrvs=lerrvt.or.lerrvD if (lerrtc .and. lerrvs) then c both thermal conductivity and viscosity generated errors eta=xnotc !values used by GUI as non-convergence flag tcx=xnotc if ((lerrtt.or.lerrvt) .and. (lerrtD.or.lerrvD)) then ierr=-57 write (herr,2057) t,rho else if (lerrtD.or.lerrvD) then ierr=-56 write (herr,2056) t,rho else ierr=-55 write (herr,2055) t,rho end if else if (lerrtc) then c only thermal conductivity generated errors tcx=xnotc !value used by GUI as non-convergence flag if (lerrtt .and. lerrtD) then ierr=-37 write (herr,2037) t,rho else if (lerrtD) then ierr=-36 write (herr,2036) Dcf,Dtcx else ierr=-35 write (herr,2035) tcf,tmint,tmaxt end if else if (lerrvs) then c only viscosity generated errors eta=xnotc !value used by GUI as non-convergence flag if (lerrvt .and. lerrvD) then ierr=-47 write (herr,2047) t,rho else if (lerrtD) then ierr=-46 write (herr,2046) Dcf,Deta else ierr=-45 write (herr,2045) tcf,tminv,tmaxv end if end if c RETURN c 2035 format ('[TRNECS warning -35] conformal temperature in ECS-', & 'transport method is outside range of reference fluid ', & 'thermal conductivity correlation; T_conf = ',f8.2, & ' K; T_min,max = ',f8.2,',',f8.2,' K') 2036 format ('[TRNECS warning -36] conformal density in ECS-', & 'transport method is outside range of reference fluid ', & 'thermal conductivity correlation; rho_conf = ',f8.4, & ' mol/L; rho_max = ',f8.4,' mol/L') 2037 format ('[TRNECS warning -37] T and rho input to ECS-', & 'transport method are outside range of reference fluid ', & 'thermal conductivity correlation; T_in = ',f8.2, & ' K; rho_in = ',f8.4,' mol/L') 2045 format ('[TRNECS warning -45] conformal temperature in ECS-', & 'transport method is outside range of reference fluid ', & 'thermal conductivity correlation; T_conf = ',f8.2, & ' K; T_min,max = ',f8.2,',',f8.2,' K') 2046 format ('[TRNECS warning -46] conformal density in ECS-', & 'transport method is outside range of reference fluid ', & 'viscosity correlation; rho_conf = ',f8.4, & ' mol/L; rho_max = ',f8.4,' mol/L') 2047 format ('[TRNECS warning -47] T and rho input to ECS-', & 'transport method are outside range of reference fluid ', & 'viscosity correlation; T_in = ',f8.2, & ' K; rho_in = ',f8.4,' mol/L') 2055 format ('[TRNECS warning -55] temperature input to ECS-', & 'transport method is outside range of reference fluid ', & 'viscosity and thermal conductivity correlations; T_in = ', & f8.2,' K; rho_in = ',f8.4,' mol/L') 2056 format ('[TRNECS warning -56] density input to ECS-', & 'transport method is outside range of reference fluid ', & 'viscosity and thermal conductivity correlations; T_in = ', & f8.2,' K; rho_in = ',f8.4,' mol/L') 2057 format ('[TRNECS warning -57] T and/or rho input to ECS-', & 'transport method are outside range of reference fluid ', & 'viscosity and thermal conductivity correlations; T_in = ', & f8.2,' K; rho_in = ',f8.4,' mol/L') c end !subroutine TRNECS c c ====================================================================== c subroutine TCBKMX (t,rho,x,fj,fx,hj,hx,tcx,Flam,lerrt,lerrD, & ierr,herr) c c compute the background thermal conductivity of a mixture where the c background t.c. is composed of the internal, dilute-gas, and residual c contributions (but not the critical enhancement) c c based on the modification of the Huber-Ely ECS method given by: c Klein, S.A., McLinden, M.O. and Laesecke, A. (1997). An improved c extended corresponding states method for estimation of viscosity of c pure refrigerants and mixtures. Int. J. Refrigeration 20:208-217 c c N.B. --equation numbers below refer to this paper; c --factor of (ref fluid mol wt)**-0.5 is missing from Eq 31 in paper c --in Eq 34 reference fluid should be evaluated at (t0,rho0), not c at (t/fj,rho*hj) c c inputs: c t--temperature [K] c rho--molar density [mol/L] c x--composition array [mol frac] c fj--array of temperature shape factors for the components c fx--temperature shape factor for the mixture c hj--array of temperature shape factors for the components c hx--density shape factor for the mixture c outputs: c tcx--thermal conductivity [W/m.K] c Flam--multiplier for t.c. (t.c._j = Flam * t.c._ref) c lerrt--error flag: .true. if temperature out of range c lerrD--error flag: .true. if density out of range c ierr--error flag: 0 = successful c <>0 indicates problem from underlying routine c passed to calling routine c herr--error string (character*255 variable if ierr<>0) c c written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO c 08-26-97 MM, original version, extracted from TRNECS c 10-01-97 MM, make correction to Eq 34 as noted above c add Flam to argument list c 10-09-97 MM, fix bug in do loop associated with calc of xmij c 10-24-97 MM, fint is now f(T) rather than const = 1.32d-3 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 (nrf0=0) !lower limit for transport ref fluid arrays character*1 htab,hnull character*3 hetamx,heta,htcxmx,htcx character*255 herr logical lerrt,lerrD dimension x(ncmax) dimension fj(ncmax),hj(ncmax) !reducing ratios for components dimension tcxin(nx),tcxdg(nx) !internal and dilute-gas parts dimension xmj(ncmax) !equivalent mol mass (Eq 33) common /CHAR/ htab,hnull common /TRNMOD/ hetamx,heta(nrf0:ncmax),htcxmx,htcx(nrf0:ncmax) c number of components in mix and constants for the mix components common /NCOMP/ nc 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 /Gcnst/ R c ierr=0 herr=hnull c find correlation limits for reference fluid p=0.0d0 call LIMITK ('TCX',0,t,rho,p,tmint,tmaxt,Dtcx,ptcx,ierr2,herr) c conformal temperature, density for mixture t0=t/fx rho0=rho*hx c reference fluid thermal conductivity at conformal conditions c check conformal t,rho against limits of reference fluid correlation call ECSLIM (t0,rho0,tmint,tmaxt,Dtcx,lerrt,lerrD,tcf,Dcf) call TCXKB (0,t0,rho0,tcx0b,ierr,herr) c calculate "equivalent mass" gx, using pure fluid residual values, c either from a pure fluid correlation or the pure fluid ECS method gxsum=0.0d0 c gxsum2=0.0d0 !alternative interpretation of gx factor c gxsum3=0.0d0 !alternative interpretation of gx factor tinsum=0.0d0 tdgsum=0.0d0 do 260 j=1,nc if (htcx(j)(1:2).ne.'EC') then c a pure fluid correlation is available c write (*,1224) j,htcx(j) c1224 format (1x,' TRNECS--pure fluid corr for j = ',i3,': ',a3) tj=t*fj(j)/fx rhoj=rho*hx/hj(j) c check that pure fluid correlation is within its limits call LIMITK ('TCX',j,tj,rhoj,p,tmn,tmx,Dmx,pmx,ierr,herr) if (ierr.eq.0) then c dilute gas visc of component j using generalized dilute-gas function etajdg=ETA0DG(j,t) c dilute-gas correlation also includes internal contribution; calculate c internal contribution separately as different mixing rule is used c tcxin(j)=1.32d-3*etajdg/wm(j)*(CP0K(j,t)-2.5d0*R) tcxin(j)=FINT(j,t)*etajdg/wm(j)*(CP0K(j,t)-2.5d0*R) call TCXK0 (j,t,tcx0,ierr2,herr) tcxdg(j)=tcx0-tcxin(j) c apply ECS method to background (residual) part of t.c. call TCXKB (j,tj,rhoj,tcxj,ierr2,herr) c write (*,1225) j,tcx0,tcxj,htcx(j) c1225 format (' TRNECS--tcx_dg,bk for comp:',i3,2f10.6,' by ',a3) end if end if if (htcx(j)(1:2).eq.'EC' .or. ierr.ne.0) then c must use ECS method to estimate c write (*,1226) j c1226 format (1x,' TRNECS--will use ECS method for t.c., j = ',i3) tchi=t/fx/(tc(j)/fj(j)) rhochi=rho*hx/(hj(j)*rhoc(j)) rho0j=rho*hx*CHI(j,tchi,rhochi) c check conformal t,rho against limits of reference fluid correlation call ECSLIM (t0,rho0j,tmint,tmaxt,Dtcx,lerrt,lerrD,tcf,Dcf) c write (*,*) ' TCBKMX--errors after ECSLIM call 1: ',lerrt,lerrD c dilute-gas and internal contributions for component j etajdg=ETA0DG(j,t) c tcxin(j)=1.32d-3*etajdg/wm(j)*(CP0K(j,t)-2.5d0*R) tcxin(j)=FINT(j,t)*etajdg/wm(j)*(CP0K(j,t)-2.5d0*R) tcxdg(j)=1.0d-3*15.0d0*R*etajdg/(4.0d0*wm(j)) tcx0=tcxin(j)+tcxdg(j) c write (*,*) ' TRNECS--j,tcxin(j): ',j,tcxin(j) c note that F-factor for t.c. has inverse power of wm compared to visc Flamj=SQRT(fj(j)*wm(0)/wm(j))*hj(j)**(-2.0d0/3.0d0) c apply ECS method to background (residual) part of t.c. call TCXKB (0,t0,rho0j,tcx0jb,ierr2,herr) tcxj=tcx0jb*Flamj c write (*,1228) j,tcx0,tcxj,Flamj c1228 format (' TRNECS--tcx_dg,bk,Flamj: comp:',i3,3f10.6,' by ECS') end if if (tcxj.gt.1.0d-10) then gj=tcx0b/tcxj*SQRT(fj(j))/hj(j)**(2.0d0/3.0d0) else c avoid possibility of division by zero (e.g. zero density as input) gj=wm(0)/wm(j)*SQRT(fj(j))/hj(j)**(2.0d0/3.0d0) end if xmj(j)=gj*gj*wm(0) c write (*,1254) j,tcxj,tcx0b,xmj(j) c1254 format (1x,' TRNECS--j,tcxj,tcx0,xmj(j) for t.c.: ',i4,3f12.6) 260 continue c do 280 j=1,nc do 276 i=1,nc xmij=2.0d0*xmj(i)*xmj(j)/(xmj(i)+xmj(j)) fij=SQRT(fj(j)*fj(i)) hij=0.125d0*(hj(j)**(1.d0/3.d0)+hj(i)**(1.d0/3.d0))**3 gxsum=gxsum+x(j)*x(i)*SQRT(fij*xmij)*hij**(4.0d0/3.0d0) c gxsum2=gxsum2+x(i)*x(j)*SQRT(fij/xmij)*hij**(-4.0d0/3.0d0) c gxsum3=gxsum3+x(i)*x(j)*SQRT(fij/xmij)*hij**(4.0d0/3.0d0) c mixing rule for the internal contribution tinsum=tinsum & +x(i)*x(j)*2.0d0*tcxin(i)*tcxin(j)/(tcxin(i)+tcxin(j)) 276 continue c assume dilute-gas part is simple mole-fraction average tdgsum=tdgsum+x(j)*tcxdg(j) 280 continue c gxroot=SQRT(wm(0)*fx)*hx**(4.0d0/3.0d0)/gxsum c Flam=SQRT(fx)*hx**(-2.0d0/3.0d0)*gxroot c above two lines reduce to: Flam=fx*hx**(2.0d0/3.0d0)*SQRT(wm(0))/gxsum c c alternative interpretation from "Transport Prop of Fluids"--Chap 12 c also, original Hanley & Ely ECS paper (1983) c gxroot=SQRT(wm(0)/fx)*hx**(4.0d0/3.0d0)*gxsum2 c Flam2=SQRT(fx)*hx**(-2.0d0/3.0d0)*gxroot c above two lines reduce to: c Flam2=hx**(2.0d0/3.0d0)*SQRT(wm(0))*gxsum2 c c alternative interpretation from Refprop5 c gxroot=SQRT(wm(0)/fx)*hx**(-4.0d0/3.0d0)*gxsum3 c Flam3=SQRT(fx)*hx**(-2.0d0/3.0d0)*gxroot c above two lines reduce to: c Flam3=SQRT(wm(0))/hx**2*gxsum3 c tcx=tcx0b*Flam+tinsum+tdgsum c tcx2=tcx0b*Flam2+tinsum+tdgsum c tcx3=tcx0b*Flam3+tinsum+tdgsum c write (*,1258) t,rho,x(1),tcx,tcx2,tcx3 c1258 format (1x,' TRNECS--t,rho,x(1);tcx by mix rules 1,2,3:',6f12.6) c write (*,1260) tdgsum,tinsum,tcx0b,tcx0c,Flam c1260 format (1x,' TRNECS--tdgsum,tinsum,tcx0b,tcx0c,Flam: ',5f12.6) c RETURN end !subroutine TCBKMX c c ====================================================================== c subroutine SETTRN (nread,icomp,hcasno,href,heos,hvs,htc,ierr,herr) c c set up working arrays for the ECS transport property model c c inputs: c nread--file to read data from c <= 0 get data from block data c >0 read from logical unit nread (file should have already c been opened and pointer set by subroutine SETUP) c icomp--component number in mixture (1..nc); 1 for pure fluid c hcasno--CAS number of component icomp (not req'd if reading from file) c c outputs: c href--file containing reference fluid EOS (character*80) c heos--model ('BWR', etc) for reference fluid EOS (character*3) c hvs--model ('VS1', etc) for ref fluid viscosity (character*3) c htc--model ('TC1', etc) for ref fluid conductivity (character*3) c ierr--error flag: 0 = successful c 1 = error (e.g. fluid not found) c herr--error string (character*255 variable if ierr<>0) c other quantities returned via arrays in common /WCFBWR/ c c written by S. Klein, NIST Thermophysics Division, Boulder, Coloradolorado c 12-14-95 SAK, original version c 03-13-96 MM, add Zcrit to common /CCON/, change parameter n0=-ncmax, c other modifications c 03-05-97 MM, new commons to match with new version of TRNECS c 08-19-97 MM, get rid of herr=herr (avoid warning); flag nread<=0 c 09-04-97 MM, add second polynomial fit (w/ crossover t or rho) c 10-24-97 MM, read in f_int term in Eucken correlation in ECS method for t.c. c 10-28-97 MM, read in f_int only for fluid file version no >= 6.001 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 (nrf0=0) !lower limit for transport ref fluid arrays parameter (mxtrn=10) !max no. coefficients for psi, chi function character*1 htab,hnull character*3 heos,hvs,htc character*12 hcasno character*80 href character*255 herr common /CHAR/ htab,hnull c limits common /WLMTRN/ tmin(nrf0:nx),tmax(nrf0:nx),pmax(nrf0:nx), & rhomax(nrf0:nx) c numbers of terms for the various parts of the model: c LJflag: flag for L-J parameters (if 0, estimate) c Euck: factor f_int in Eucken correlation c psi (viscosity shape factor): polynomial term, 2nd poly, spare c chi (conductivity shape factor): polynomial term, 2nd poly, spare common /WNTTRN/ LJflag(nrf0:nx),nEuck(nrf0:nx), & npsi1(nrf0:nx),npsi2(nrf0:nx),npsi3(nrf0:nx), & nchi1(nrf0:nx),nchi2(nrf0:nx),nchi3(nrf0:nx) c commons storing the (real and integer) coefficients to the ECS model common /WCFTRN/ cpsi(nrf0:nx,mxtrn,4),cchi(nrf0:nx,mxtrn,4) common /WIFTRN/ ipsi(nrf0:nx,0:mxtrn),ichi(nrf0:nx,0:mxtrn) c Lennard-Jones parameters common /WLJTRN/ sigma(nrf0:nx),epsk(nrf0:nx) c coefficients to f_int term in Eucken correlation for therm cond common /WCEUCK/ cEuck(nrf0:nx,mxtrn,4) common /VERS/ verfl(n0:nx),vermx !fluid & mix file version nos. c if (nread.le.0) then ierr=101 write (herr,1101) nread,hcasno,hnull 1101 format ('[SETTRN error 101] illegal file specified; nread = ', & i4,'; CAS no. = ',a12,a1) RETURN else herr=hnull ierr=0 end if c c read data from file c write (*,*) ' SETTRN--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,2083) heos,href !reference fluid EOS and .fld file read (nread,2003) hvs !reference fluid viscosity model read (nread,2003) htc !reference fluid conductivity model read (nread,*) LJflag(icomp) !Lennard-Jones flag read (nread,*) sigma(icomp) !Lennard-Jones coef Sig read (nread,*) epsk(icomp) !Lennard-Jones coef EPS c write (*,*) ' SETTRN--L-J parameters: ',sigma(icomp),epsk(icomp) c c read number of terms for f_int in Eucken correlation c write (*,*) ' SETTRN--icomp, version #: ',icomp,verfl(icomp) if (verfl(icomp).ge.6.0009d0) then read (nread,*) nEuck(icomp) if (nEuck(icomp).ge.1) then c read polynomial term(s) for f_int do 110 j=1,nEuck(icomp) c read coeff, power of T, spare1, spare2 read (nread,*) (cEuck(icomp,j,k),k=1,4) c write (*,*) ' SETTRN--Eucken par:',(cEuck(icomp,j,k),k=1,2) 110 continue end if else nEuck(icomp)=0 end if c c read number of terms for viscosity shape factor (incl. spare) read (nread,*) npsi1(icomp),npsi2(icomp),npsi3(icomp) if (npsi1(icomp).ge.1) then c read polynomial term do 120 j=1,npsi1(icomp) c read coeff, power of Tr, power of Dr, spare read (nread,*) (cpsi(icomp,j,k),k=1,4) c write (*,*) ' SETTRN--psi par: ',(cpsi(icomp,j,k),k=1,3,2) 120 continue end if if (npsi2(icomp).ge.1) then c read coeff, power of Tr, power of Dr, spare for 2nd polynomial term c first set of coeff is crossover (t or rho) do 130 j=npsi1(icomp)+1,npsi1(icomp)+npsi2(icomp) read (nread,*) (cpsi(icomp,j,k),k=1,4) 130 continue end if c c ditto for thermal conductivity shape factor read (nread,*) nchi1(icomp),nchi2(icomp),nchi3(icomp) if (nchi1(icomp).ge.1) then do 160 j=1,nchi1(icomp) c read coeff, power of Tr, power of Dr, spare read (nread,*) (cchi(icomp,j,k),k=1,4) c write (*,*) ' SETTRN--chi par: ',(cchi(icomp,j,k),k=1,3,2) 160 continue end if if (nchi2(icomp).ge.1) then c read coeff, power of Tr, power of Dr, spare for 2nd polynomial term c first set of coeff is crossover (t or rho) do 170 j=nchi1(icomp)+1,nchi1(icomp)+nchi2(icomp) read (nread,*) (cchi(icomp,j,k),k=1,4) 170 continue end if c RETURN 2003 format (a3) 2083 format (a3,1x,a80) end !subroutine SETTRN c c ====================================================================== c function FINT (icomp,t) c c factor f_int appearing in Eucken correlation for thermal conductivity c c inputs: c icomp--component number in mixture (1..nc); 1 for pure fluid c t--temperature [K] c output (as function value): c FINT--the factor f_int c c written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO c 10-24-97 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 (nx=ncmax) parameter (nrf0=0) !lower limit for transport ref fluid arrays parameter (mxtrn=10) !max no. coefficients for psi, chi function c c numbers of terms for the various parts of the model: c LJflag: flag for L-J parameters (if 0, estimate) c Euck: factor f_int in Eucken correlation c psi (viscosity shape factor): polynomial term, 2nd poly, spare c chi (conductivity shape factor): polynomial term, 2nd poly, spare common /WNTTRN/ LJflag(nrf0:nx),nEuck(nrf0:nx), & npsi1(nrf0:nx),npsi2(nrf0:nx),npsi3(nrf0:nx), & nchi1(nrf0:nx),nchi2(nrf0:nx),nchi3(nrf0:nx) c coefficients to f_int term in Eucken correlation for therm cond common /WCEUCK/ cEuck(nrf0:nx,mxtrn,4) c if (nEuck(icomp).le.0) then c no correlation for f_int is present, use value corresponding to c modified Eucken correlation FINT=1.32d-3 else fsum=0.0d0 do 160 k=1,nEuck(icomp) fsum=fsum+cEuck(icomp,k,1)*t**cEuck(icomp,k,2) 160 continue FINT=fsum end if c write (*,*) ' FINT--t,f_int: ',t,FINT c RETURN end !function FINT c c ====================================================================== c function PSI (icomp,tr,rhor) c c viscosity shape factor for a pure fluid or mixture component, as c defined by: Klein et al., Int J Refrigeration 20:208-217 (1997) c c inputs: c icomp--component number in mixture (1..nc); 1 for pure fluid c tr--reduced temperature; = t/tc for a pure fluid c = (t/fx)/(tc_j/fj) for a mixture c rhor--reduced density; = rho/rhoc for a pure fluid c = rho*hx/(hj*rhoc_j) for a mixture c output (as function value): c psi--the viscosity shape factor, i.e. an additional factor entering c into the definition of the conformal density in the ECS method c such that: rho_0 = rho*hj*psi c c written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO c 03-05-97 MM, original version c 09-04-97 MM, add second polynomial fit (w/ crossover t or rho) c 10-24-97 MM, changes in /WNTTRN/ to accomodate Eucken correlation c implicit double precision (a-h,o-z) implicit integer (i-n) parameter (ncmax=5) !max number of components in mixture parameter (nx=ncmax) parameter (nrf0=0) !lower limit for transport ref fluid arrays parameter (mxtrn=10) !max no. coefficients for psi, chi function c c numbers of terms for the various parts of the model: c LJflag: flag for L-J parameters (if 0, estimate) c Euck: factor f_int in Eucken correlation c psi (viscosity shape factor): polynomial term, 2nd poly, spare c chi (conductivity shape factor): polynomial term, 2nd poly, spare common /WNTTRN/ LJflag(nrf0:nx),nEuck(nrf0:nx), & npsi1(nrf0:nx),npsi2(nrf0:nx),npsi3(nrf0:nx), & nchi1(nrf0:nx),nchi2(nrf0:nx),nchi3(nrf0:nx) c commons storing the (real and integer) coefficients to the ECS model common /WCFTRN/ cpsi(nrf0:nx,mxtrn,4),cchi(nrf0:nx,mxtrn,4) common /WIFTRN/ ipsi(nrf0:nx,0:mxtrn),ichi(nrf0:nx,0:mxtrn) c i=icomp if (npsi1(i).le.0) then c no transport shape factor is present PSI=1.0d0 else psisum=0.0d0 if (npsi2(i).ge.2) then c a second polynomial term is present; check whether input (t or rho) c is below the crossover value (k = 1,2 in cpsi(i,npsi1(i)+1,k)); c this term allows for a piece-wise fit of PSI nz=npsi1(i)+1 if (tr.lt.cpsi(i,nz,1). or. rhor.lt.cpsi(i,nz,2)) then do 120 k=nz+1,npsi1(i)+npsi2(i) psisum=psisum+cpsi(i,k,1)*tr**cpsi(i,k,2)*rhor**cpsi(i,k,3) 120 continue PSI=psisum RETURN end if end if c apply the first polynomial term (either it's the only one present c or the 2nd polynomial term does not apply) do 160 k=1,npsi1(i) psisum=psisum+cpsi(i,k,1)*tr**cpsi(i,k,2)*rhor**cpsi(i,k,3) 160 continue PSI=psisum end if c write (*,*) ' PSI: ',PSI c RETURN end !function PSI c c ====================================================================== c function CHI (icomp,tr,rhor) c c thermal conductivity shape factor for a pure fluid or mixture c component, analogous to the viscosity shape factor defined by c Klein et al., Int J Refrigeration 20:208-217 (1997) c c inputs: c icomp--component number in mixture (1..nc); 1 for pure fluid c tr--reduced temperature; = t/tc for a pure fluid c = (t/fx)/(tc_j/fj) for a mixture c rhor--reduced density; = rho/rhoc for a pure fluid c = rho*hx/(hj*rhoc_j) for a mixture c output (as function value): c chi--thermal conductivity shape factor, i.e. an additional factor c entering into the definition of the conformal density in c the ECS method such that: rho_0 = rho*hj*psi c c written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO c 03-05-97 MM, original version c 09-04-97 MM, add second polynomial fit (w/ crossover t or rho) c 10-24-97 MM, changes in /WNTTRN/ to accomodate Eucken correlation c implicit double precision (a-h,o-z) implicit integer (i-n) parameter (ncmax=5) !max number of components in mixture c parameter (n0=-ncmax,nx=ncmax) parameter (nx=ncmax) parameter (nrf0=0) !lower limit for transport ref fluid arrays parameter (mxtrn=10) !max no. coefficients for psi, chi function c c numbers of terms for the various parts of the model: c LJflag: flag for L-J parameters (if 0, estimate) c Euck: factor f_int in Eucken correlation c psi (viscosity shape factor): polynomial term, 2nd poly, spare c chi (conductivity shape factor): polynomial term, 2nd poly, spare common /WNTTRN/ LJflag(nrf0:nx),nEuck(nrf0:nx), & npsi1(nrf0:nx),npsi2(nrf0:nx),npsi3(nrf0:nx), & nchi1(nrf0:nx),nchi2(nrf0:nx),nchi3(nrf0:nx) c commons storing the (real and integer) coefficients to the ECS model common /WCFTRN/ cpsi(nrf0:nx,mxtrn,4),cchi(nrf0:nx,mxtrn,4) common /WIFTRN/ ipsi(nrf0:nx,0:mxtrn),ichi(nrf0:nx,0:mxtrn) c i=icomp if (nchi1(i).le.0) then c no transport shape factor is present CHI=1.0d0 else chisum=0.0d0 if (nchi2(i).ge.2) then c a second polynomial term is present; check whether input (t or rho) c is below the crossover value (k = 1,2 in cchi(i,nchi1(i)+1,k)); c this term allows for a piece-wise fit of CHI nz=nchi1(i)+1 if (tr.lt.cchi(i,nz,1). or. rhor.lt.cchi(i,nz,2)) then do 120 k=nz+1,nchi1(i)+nchi2(i) chisum=chisum+cchi(i,k,1)*tr**cchi(i,k,2)*rhor**cchi(i,k,3) 120 continue CHI=chisum RETURN end if end if c apply the first polynomial term (either it's the only one present c or the 2nd polynomial term does not apply) do 160 k=1,nchi1(i) chisum=chisum+cchi(i,k,1)*tr**cchi(i,k,2)*rhor**cchi(i,k,3) 160 continue CHI=chisum end if c write (*,*) ' CHI: ',CHI c RETURN end !function CHI c c ====================================================================== c subroutine TCKVIR (icomp,t,tcx0,tcxcol,tcxvir,ierr,herr) c c thermal conductivity virial coefficient for a pure fluid or mixture c component, based on: c Nieto de Castro, C.A., Friend, D.G., Perkins, R.A. and Rainwater, c J.C. (1990). Thermal conductivity of a moderately dense gas. c Chemical Physics 145: 19-26. c equation numbers in the comments refer to this paper c c inputs: c icomp--component number in mixture (1..nc); 1 for pure fluid c t--temperature [K] c output (as function value): c tcx0--dilute gas thermal conductivity [W/m-K] c tcxcol--contribution of collisions to the dilute-gas t.c. [W/m-K] c tcxvir--thermal conductivity second virial coefficient [L/mol]; c i.e. the multiplier times the collisional part of the dilute- c gas conductivity which gives the initial density dependence c of that part of the thermal conductivity c ierr--error flag: 0 = successful c <>0 = error code originating in TCXK0 c herr--error string (character*255 variable if ierr<>0) c c written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO c 09-24-97 MM, original version c 10-24-97 MM, fint is now f(T) rather than const = 1.32d-3 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 (nrf0=0) !lower limit for transport ref fluid arrays character*1 htab,hnull character*3 hetamx,heta,htcxmx,htcx character*255 herr common /CHAR/ htab,hnull common /TRNMOD/ hetamx,heta(nrf0:ncmax),htcxmx,htcx(nrf0:ncmax) c Lennard-Jones parameters common /WLJTRN/ sigma(nrf0:nx),epsk(nrf0:nx) c constants for the mix components 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 /Gcnst/ R c ierr=0 herr=hnull c etadg=ETA0DG(icomp,t) tcxcol=1.0d-3*15.0d0*R*etadg/(4.0d0*wm(icomp)) !Eq 7 if (htcx(1)(1:2).eq.'EC') then c component is modeled with ECS--use generalized function for internal c contributions to t.c. c tcxint=1.32d-3*etadg/wm(icomp)*(CP0K(icomp,t)-2.5d0*R) tcxint=FINT(icomp,t)*etadg/wm(icomp)*(CP0K(icomp,t)-2.5d0*R) tcx0=tcxcol+tcxint else c use fluid-specific correlation call TCXK0 (icomp,t,tcx0,ierr,herr) end if tstar=t/epsk(icomp) bprime=(2.9749d0+0.1140d0*tstar)/(tstar-0.04953d0) !Eq 12 blam=(bprime-0.625d0*(tcx0/tcxcol-1.0d0))/(tcx0/tcxcol) c Eq 3, where the const = 2*pi*N0/3 *1d-27 (sigma in nm) *1d3 (vol in L) tcxvir=blam*1.26127336d0*sigma(icomp)**3 c RETURN end !subroutine TCKVIR c c ====================================================================== c subroutine ECSLIM (t,D,tmin,tmax,Dmax,lerrt,lerrD,terr,Derr) c c check input t,rho against limits and return error flags c c inputs: c t--temperature [K] c D--molar density [mol/L] c tmin--minimum temperature for model [K] c tmax--maximum temperature [K] c Dmax--maximum density [mol/L] c outputs: c lerrt--logical flag; .true. if t outside limits, set only if .true. c lerrD--logical flag; .true. if D outside limits, set only if .true. c terr--same as input t, but set only if lerrt = .true. c Derr--same as input D, but set only if lerrD = .true. c c written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO c 03-24-97 MM, original version c implicit double precision (a-h,o-z) implicit integer (i-k,m,n) logical lerrt,lerrD c if (t.gt.1.5d0*tmax .or. t.lt.tmin) then lerrt=.true. terr=t end if if (D.gt.Dmax) then lerrD=.true. Derr=D end if c RETURN end c c ====================================================================== c function ETA0DG (icomp,t) c c dilute-gas contribution to viscosity for use with ECS model c c inputs: c icomp--component number in mixture (1..nc); 1 for pure fluid c t--temperature [K] c output (as function value): c eta0dg--the dilute-gas part of the viscosity [uPa-s] c c written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO c 08-25-97 MM, original version, based on ETA1 function by S.A. Klein 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 (nrf0=0) !lower limit for transport ref fluid arrays c c common storing the fluid constants 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 Lennard-Jones parameters common /WLJTRN/ sigma(nrf0:nx),epsk(nrf0:nx) c i=icomp tau=t/epsk(i) c in this case, the dilute gas is simply the Chapman-Enskog term ETA0DG=26.692d-3*SQRT(wm(i)*t)/(sigma(i)**2*OMEGAS(2,2,tau)) c write (*,1001) i,t,tau,ETA0DG c1001 format (1x,' ETA0DG--dilute-gas visc: i,t,tau,eta',i4,2f9.4,f12.6) c write (*,*) ' ETA0DG--sigma,omega_2,2: ',sigma(i),OMEGAS(2,2,tau) c RETURN end !function ETA0DG c c ====================================================================== c function OMEGAS (il,is,tau) c c collision integral for Lenard-Jones fluid; returns value for c Omega_1,1 or Omega_2,2; based on: c Neufeld, Janzen, and Aziz. (1972). J Chem Phys 57:1100-1102 c c inputs: c il--order of integral c is--order of integral c tau--dimensionless temperature = t/epsk, where epsk is the Lennard- c Jones energy parameter (eplison/k) c output (as function value): c OMEGAS--value of collision integral c c N.B.--only inputs of (il = is = 1) and (il = is = 2) are valid c c originally implemented in MIPROPS/SUPERTRAPP by J.F. Ely c 08-25-97 MM, revised and documented by M. McLinden, c NIST Physical & Chemical Properties Div, Boulder, CO c implicit double precision (a-h,o-z) implicit integer (i-n) c if (il.eq.1 .and. is.eq.1) then OMEGAS=1.06036d0/tau**0.15610d0 & +0.19300d0*EXP(-0.47635d0*tau) & +1.03587d0*EXP(-1.52996d0*tau) & +1.76474d0*EXP(-3.89411d0*tau) else OMEGAS=1.16145d0/tau**0.14874d0 & +0.52487d0*EXP(-0.77320d0*tau) & +2.16178d0*EXP(-2.43787d0*tau) end if c write (*,1001) il,is,tau,omegas c1001 format (1x,' OMEGAS--il,is,tau,omega: ',2i3,2f12.7) c RETURN end !function OMEGAS c c ====================================================================== c FUNCTION DELHSV (TX,DX,X,hj) C C ENSKOG CORRECTION FOR SIZE AND MASS DIFFERENCE EFFECTS C IN MIXTURE VISCOSITY PREDICTION C C BASED ON : J.F. ELY, J. RES. NBS 86(6) 1981, P597-604 C IMPLICIT DOUBLE PRECISION (A-H,O-Z) c parameter (ncmax=5) !max number of components in mixture parameter (n0=-ncmax,nx=ncmax) dimension hj(ncmax) common /NCOMP/ nc c /NCOMP/ contains the number of components 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 /CCON/: constants for each of the NCMAX components; these values c are taken from the coefficients/array of selected model c wm(i): molecular mass (g/mol) c ttp(i): triple point temperature (K) c tnbp(i): normal boiling point temperature (K) c tc(i): critical temperature (K) c pc(i): critical pressure (kPa) c rhoc(i): critical density (mol/L) c zcrit(i) critical compressibility c accen(i): acentric factor c dipole(i) dipole moment (debye) at normal boiling pt c DIMENSION X(nx), Z(nx), S(nx), SIGMA(nx,nx) DIMENSION SIGDUM(nx,nx), CMWDUM(nx), ZDUM(nx) DIMENSION CMW(ncmax) SAVE DATA DCON, IONE / 6.023D-4, 1 / C SIGDUM(1,1) = 1.D0 CMWDUM(1) = 1.D0 ZDUM(1) = 1.D0 CMWN = wm(NC) DO 020 N = 1, NC CMW(N)=wm(N) Z(N) = X(N) IF (X(N).LE.0.0D0) Z(N) = 1.0D-8 S(N) = CBRT(200.0D0 * hj(N)/(3.058D0 * 0.6023D0) ) CMW(N) = CMW(N) / CMWN 020 CONTINUE C SN = S(NC) DO 025 N = 1, NC 025 S(N) = S(N) / SN RHOX = DCON * DX * SN**3 SIG1 = 0.0D0 CMW1 = 0.0 DO 060 I=1,NC SI = 0.0D0 CI = 0.0D0 DO 040 J = 1,NC SIJ = 0.5D0 * (S(I) + S(J)) TERM = Z(J) * SIJ**3 SI = SI + TERM CI = CI + SIJ * TERM * SQRT(CMW(I) * CMW(J) / (CMW(I)+CMW(J))) SIGMA(I,J) = SIJ 040 CONTINUE SIG1 = SIG1 + Z(I) * SI CMW1 = CMW1 + Z(I) * CI 060 CONTINUE TERM = CBRT(SIG1) CMW1 = 2.0 * CMWN * (CMW1 / (SIG1 * TERM ))**2 SIG1 = SN * TERM RHO1 = DCON * DX * SIG1**3 C c change arg from nc -> ncc (should not pass element in common) ncc=NC CALL ENSKOG(ncc,RHOX,SIGMA,CMW,Z,VISX) VIS0 = 26.692D-3 * SQRT(CMWN*TX) / (SN * SN) !{micro Pa-sec} VISX = VIS0 * VISX C CALL ENSKOG(IONE,RHO1,SIGDUM,CMWDUM,ZDUM,VIS1) VIS0 = 26.692D-3 * SQRT(CMW1*TX) / (SIG1*SIG1) !{micro Pa-sec} VIS1 = VIS0 * VIS1 DO 080 N = 1, NC 080 CMW(N) = CMW(N) * CMWN C DELHSV = VISX - VIS1 C C DONT LET CORRECTION BE NEGATIVE IF(DELHSV.LT.0.0d0) DELHSV = 0.0D0 RETURN END !function DELHSV c c ====================================================================== c SUBROUTINE ENSKOG (N,RHO,SIGMA,CMW,X,ETA) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION MJI parameter (ncmax=5) !max number of components in mixture c parameter (n0=-ncmax,nx=ncmax) parameter (nx=ncmax) C DIMENSION SIGMA(nx,nx), CMW(nx), X(nx), Y(nx) DIMENSION H(nx,21), S(nx,nx) SAVE DATA PI6 / 0.523598776D0 / C C PRELIMINARY CALCULATIONS ETA = 0.0D0 TERM = RHO * PI6 S2 = 0.0D0 S3 = 0.0D0 DO 020 I=1,N TEMP = X(I) * TERM * SIGMA(I,I)**2 S2 = S2 + TEMP 020 S3 = S3 + TEMP * SIGMA(I,I) S3 = 1.0 - S3 A1 = S3 * S3 A2 = S3 * S2 A3 = S2 * S2 A4 = S3 * A1 C DO 060 I = 1, N Y(I) = 0.0 ETA2 = 0.0 DO 040 J = 1, N S2 = SIGMA(I,I) * SIGMA(J,J) / (SIGMA(I,I) + SIGMA(J,J)) S3 = 2.0 * S2 * S2 S2 = 3.0 * S2 BIJ = 4.0 * PI6 * RHO * SIGMA(I,J)**3 YIJ = BIJ * (A1 + A2*S2 + A3*S3) / A4 MJI = CMW(J) / (CMW(I) + CMW(J)) EIJ = BIJ * SQRT(2.0D0 * CMW(I) * MJI) / SIGMA(I,J)**2 Y(I) = Y(I) + X(J) * MJI * YIJ ETA2 = ETA2 + X(J) * EIJ * YIJ S(I,J) = X(J) * YIJ * MJI * MJI / EIJ 040 CONTINUE Y(I) = X(I) * (1.0D0 + 0.8D0 * Y(I)) ETA = ETA + X(I) * ETA2 060 CONTINUE C GENERATE THE H MATRIX DO 100 I = 1, N DO 100 J = 1, N S2 = 0.0D0 A1 = 0.0D0 IF (J.EQ.I) A1 = 1.0D0 DO 080 L = 1, N RATIO = CMW(I) / (3.0D0 * CMW(L)) A2 = 0.0D0 IF (L.EQ.J) A2 = 2.0D0 TERM = A1 * (1.0D0 + 5.0D0 * RATIO) - A2 * RATIO S2 = S2 + S(I,L) * TERM 080 CONTINUE H(I,J) = 2.0D0 * X(I) * S2 100 CONTINUE C SOLVE FOR THE EXPANSION COEFFICIENTS N1 = N + 1 IF (N.GT.1) GO TO 110 H(1,2) = Y(1) / H(1,1) GO TO 190 C 110 DO 120 I = 1, N 120 H(I,N1) = Y(I) DO 180 I = 1, N I1 = I + 1 DO 140 J = I1, N1 140 H(I,J) = H(I,J) / H(I,I) H(I,I) = 1.0 DO 180 J = 1, N IF (J.EQ.I) GO TO 180 DO 170 K = I1, N1 170 H(J,K) = H(J,K) - H(J,I) * H(I,K) H(J,I) = 0.0D0 180 CONTINUE 190 ETA1 = 0.0D0 DO 200 I = 1, N 200 ETA1 = ETA1 + H(I,N1) * Y(I) ETA = ETA1 + 8.0D0 * ETA / (25.0D0 * PI6) RETURN END !subroutine ENSKOG c c ====================================================================== c function ETAMIX (t,x) c c compute the viscosity of a dilute-gas mixture assuming the components c are described by a Lennard-Jones potential c c inputs: c t--temperature [K] c x--composition array [mol frac] c outputs: c eta--viscosity [uPa.s] c c source of original version lost to history (possibly J.F. Ely) c adopted for use in Refprop by S.A. Klein, January, 1996 c 08-25-97 MM, revised and documented by M. McLinden, c NIST Physical & Chemical Properties Div, Boulder, CO c name changed from ETA0 to ETAMIX 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 (nrf0=0) !lower limit for transport ref fluid arrays parameter (nx1=nx+1) c number of components in mix and constants for the mix components common /NCOMP/ nc 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 Lennard-Jones parameters common /WLJTRN/ sig(nrf0:nx),eps(nrf0:nx) c DIMENSION EOK(nx,nx), SIGMA(nx,nx), RM(nx,nx), H(nx,nx1),X(nx) DATA DCON, VCON2 / 0.320300D0, 26.692D-3/ C TK=t NC1 = NC + 1 DO 020 I = 1, NC DO 020 J = 1, NC RM(I,J) = 2.0 * WM(I) * WM(J) / (WM(I) + WM(J)) EOK(I,J) = SQRT(EPS(I) * EPS(J)) 020 SIGMA(I,J) = 0.5 * (SIG(I) + SIG(J)) C 040 DO 080 I = 1, NC DO 080 J = 1, NC DIJ = 0.0 IF (I.EQ.J) DIJ = 1.0 SUML = 0.0 DO 060 L = 1, NC TS = TK / EOK(I,L) TERM = SQRT(RM(I,L)*TK) / SIGMA(I,L)**2 RHODIL = DCON * TERM / OMEGAS(1,1,TS) ETAIL = VCON2 * TERM / OMEGAS(2,2,TS) DJL = 0.0 IF (J.EQ.L) DJL = 1.0 TERM = (RM(J,J)/(RHODIL*RM(L,L)))*(DIJ-DJL) + 0.5*(DIJ+DJL)/ETAIL SUML = SUML + X(L) * RM(I,L) * RM(I,L) * TERM 060 CONTINUE H(I,J) = SUML / (RM(I,I) * RM(J,J)) 080 CONTINUE DO 100 I = 1, NC 100 H(I,NC1) = 1.0 DO 160 I = 1, NC I1 = I + 1 DO 120 J = I1, NC1 120 H(I,J) = H(I,J) / H(I,I) H(I,I) = 1.0 DO 160 J = 1, NC IF (J.EQ.I) GO TO 160 DO 140 K = I1, NC1 140 H(J,K) = H(J,K) - H(J,I) * H(I,K) H(J,I) = 0.0 160 CONTINUE ETA = 0.0 DO 180 I = 1, NC 180 ETA = ETA + X(I) * H(I,NC1) ETAMIX = ETA RETURN END !function ETAMIX c c ====================================================================== c subroutine CONFTD (j,amix,Zmix,tj,rhoj,ierr,herr) c c Find the conformal temperature and density for component j (which c might be the reference fluid) which give a reduced residual Helmholtz c energy and compressibility which match the input value. The system c of equations to be solved is thus: c aj(tj,rhoj) - amix = 0 c Zj(tj,rhoj) - Zmix = 0 c these are put into dimensionless temperature and density tau and del c and linearized to the form: c (da/dtau)*deltau + (da/ddel)*deldel = amix - aj(tj,rhoj) c (dZ/dtau)*deltau + (dZ/ddel)*deldel = Zmix - Zj(tj,rhoj) c to allow a solution by the classic Newton-Raphson method in two c dimensions. In the following code, the above nomenclature becomes: c a11*x1 + a12*x2 = fx1 c a21*x1 + a22*x2 = fx2 c c inputs: c j--component number c j = 0: find conformal conditions for reference fluid c j > 0: find conformal conditions for component j c amix--reduced residual Helmholtz energy (A - A*)/RT [-] c Zmix--compressibility (pV/RT) [-] c tj--initial guess for conformal temperature for component j [K] c rhoj--initial guess for conformal density for component j [mol/L] c outputs: c tj--converged conformal temperature for component j [K] c rhoj--converged conformal density for component j [mol/L] c ierr--error flag: 0 = successful c -58 = did not converge c herr--error string (character*255 variable if ierr<>0) c c written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO c 09-29-97 MM, original version, replaces MINH by S.A. Klein c 01-13-98 MM, check for error from CONFD in case CONFTD does not converge 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) character*1 htab,hnull character*255 herr,herr1 common /CHAR/ htab,hnull 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) data itmax /25/ data tolf /1.0d-7/ !convergence tol for function evaluations data tolx /1.0d-12/ !tolerance for step in independent variables c c save initial guesses in case of non-convergence tjsav=tj rhojsv=rhoj c find the reducing ratios for component j call REDK (j,tredj,Dredj) c write (*,*) ' CONFTD--reducing t,rho: ',tredj,Dredj x1=tredj/tj !tauj (dimensionless temperature) x2=rhoj/Dredj !delj (dimensionless density) c write (*,1788) j,tj,rhoj,x1,x2 c1788 format (1x,' CONFTD--input j,tj,rhoj: ',i2,2e14.6, c & '; resulting tau,del: ',2e14.6) c begin iteration do 800 it=0,itmax c check for resonableness of guess by means of derivative dP/d(rho) c [actually, the quantity calculated is (1/RT)*(dP/drho)] phi00=PHIK(j,0,0,x1,x2) phi01=PHIK(j,0,1,x1,x2) phi02=PHIK(j,0,2,x1,x2) dpdrho=1.0d0+2.0d0*x2*phi01+x2*x2*phi02 if (dpdrho.lt.0.0d0) then c write (*,1799) j,it,x1,x2,dpdrho c1799 format (' CONFTD--j,it,tau,del: (dpdrho < 0),dP/dD: ',2i3,2e15.6, c & 30x,e15.6) if (x2*Dredj.gt.rhoc(j)) then c liquid-phase, go to higher density x2=1.04d0*x2 else c vapor-phase, go to higher temperature (lower tau) x1=0.96d0*x1 end if goto 800 !use up one iteration (avoid infinite loop) end if c evaluate objective function fx1=amix-phi00 fx2=Zmix-1.0d0-x2*phi01 c write (*,1800) j,it,x1,x2,fx1,fx2,dpdrho c1800 format (' CONFTD--j,it,tau,del,amix-a,Zmix-Z,dP/dD: ',2i3,5e15.6) c check for convergence error=ABS(fx1)+ABS(fx2) if (error.lt.tolf) then c iteration has converged goto 840 end if c c calc remaining derivatives and new guess for independent variables phi10=PHIK(j,1,0,x1,x2) phi11=PHIK(j,1,1,x1,x2) a11=phi10 !partial a w.r.t. tau a12=phi01 !partial a w.r.t. del a21=x2*phi11 !partial Z w.r.t. tau a22=x2*phi02+phi01 !partial Z w.r.t. del denom=a11*a22-a21*a12 if (ABS(denom).lt.1.0d-16) then c system has singularity, no solution for this guess,try another delx1=0.05d0*x1 delx2=-0.05d0*x2 else delx1=(a22*fx1-a12*fx2)/denom delx2=(a11*fx2-a21*fx1)/denom end if sumdel=ABS(delx1)+ABS(delx2) c do not allow too great a step (x1,x2 should always be positive) x1=x1+SIGN(1.0d0,delx1)*MIN(ABS(delx1),0.05d0*x1) x2=x2+SIGN(1.0d0,delx2)*MIN(ABS(delx2),0.05d0*x2) c if step is within tolerance, also consider iteration finished if (sumdel.lt.tolx) then goto 840 end if 800 continue c c iteration loop has not converged, try another scheme ierr=-58 write (herr,1058) j,hnull 1058 format ('[TRNECS warning 58] 2-D Newton-Raphson method for ', & 'conformal temperature and density did not converge ', & 'for component',i3,a1) call ERRMSG (ierr,herr) c write (*,1801) x1,x2,fx1,fx2 c1801 format (' CONFTD--2-D Newton''s method iteration did not', c & ' converge; x1,x2,fx1,fx2: ',4e15.6) c go back to starting guesses for tj, rhoj tj=tjsav rhoj=rhojsv call CONFD (j,amix,Zmix,tj,rhoj,ierr1,herr1) c write (*,*) ' CONFTD--ierr,tj from CONFD: ',ierr1,tj if (ierr1.ne.0) then c write (*,*) ' CONFD--did not converge; using initial tj,rhoj' ierr=ierr1 herr=herr1 tj=tjsav rhoj=rhojsv end if RETURN c c iteration has converged, assign outputs 840 continue tj=tredj/x1 rhoj=x2*Dredj c RETURN end !subroutine CONFTD c c ====================================================================== c subroutine CONFD (j,amix,Zmix,tj,rhoj,ierr,herr) c c Find the conformal temperature and density for component j (which c might be the reference fluid) which give a reduced residual Helmholtz c energy and compressibility which match the input value. The quantity c FX = [a(tj,rhoj) - amix]**2 + [Z(tj,rhoj) - Zmix]**2 c is minimized using Brent's method adapted from: c Press, W.H., Flannery, B.P., Teukolsky, S.A. and Vettering, W.T. c (1986). Numerical Recipes: The Art of Scientific Computing, c Cambridge University Press. c This minimization is in terms of density. The single-dimension c Brent's method is extended to the two-dimensional problem here by an c auxiliary routine CONFT which finds the optimum conformal temperature c for a given guess of the conformal density. c c inputs: c j--component number c j = 0: find conformal conditions for reference fluid c j > 0: find conformal conditions for component j c amix--reduced residual Helmholtz energy (A - A*)/RT [-] c Zmix--compressibility (pV/RT) [-] c tj--initial guess for conformal temperature for component j [K] c rhoj--initial guess for conformal density for component j [mol/L] c outputs: c tj--converged conformal temperature for component j [K] c rhoj--converged conformal density for component j [mol/L] c ierr--error flag: 0 = successful c -58 = did not converge c herr--error string (character*255 variable if ierr<>0) c c written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO c 09-26-97 MM, original version, based on MINH by S.A. Klein c implicit double precision (a-h,o-z) implicit integer (i-n) character*1 htab,hnull character*255 herr common /CHAR/ htab,hnull DATA CGOLD / 0.61803399/, ZEPS/1E-20/, TOL/1E-7/, ITMAX/50/ c ierr=0 herr=hnull V=rhoj !first guess for rhoj A=0.80d0*rhoj !lower bound on rhoj B=1.20d0*rhoj !upper bound on rhoj W=V X=V E=0.0d0 D=0.0d0 c find the reducing ratios for component j call REDK (j,tredj,Dredj) c find the optimum tj for the initial guess of rhoj and evaluate the c objective function FX c write (*,*) ' CONFD--j,tj input to CONFT, first iteration: ',j,tj tj0=tj !save initial guess call CONFT (j,amix,rhoj,tj,ierr,herr) if (ABS(tj/tj0-1.0d0).ge.0.2) then c do not allow too great a change in tj c write (*,*) ' CONFD--unconstrained tj from CONFT: ',tj tj=tj0*(1.0d0+SIGN(0.2d0,tj/tj0-1.0d0)) end if c write (*,*) ' CONFD--ierr,tj from CONFT: ',ierr,tj tauj=tredj/tj delj=rhoj/Dredj aj=PHIK(j,0,0,tauj,delj) dela=amix-aj Zj=1.0d0+delj*PHIK(j,0,1,tauj,delj) delZ=Zmix-Zj FX=dela**2+delZ**2 c call DPDDK (j,tj,rhoj,dpdrho) !debug only c write (*,1089) j,0,tj,rhoj,amix,aj,Zmix,Zj,dpdrho c FV=FX FW=FX DO 100 ITER=1,ITMAX XM=0.5d0*(A+B) TOL1=TOL*ABS(X)+ZEPS TOL2=2.0d0*TOL1 IF (ABS(X-XM).LE.TOL2-0.5d0*(B-A)) GOTO 99 IF (ABS(E).GT.TOL1) THEN R=(X-W)*(FX-FV) Q=(X-V)*(FX-FW) P=(X-V)*Q-(X-W)*R Q=2.0*(Q-R) IF (Q.GT.0.0d0) P=-P Q=ABS(Q) ETEMP=E E=D IF ((ABS(P).GE.ABS(0.5d0*Q*ETEMP)).or.(P.LE.Q*(A-X)).or. & (P.GE.Q*(B-X))) THEN IF (X.GE.XM) THEN E=A-X ELSE E=B-X ENDIF D=CGOLD*E ELSE D=P/Q !Parabolic step U=X+D IF ((U-A.LT.TOL2).or.(B-U.LT.TOL2)) D=SIGN(TOL1,XM-X) ENDIF ELSE IF (X.GE.XM) THEN E=A-X ELSE E=B-X ENDIF D=CGOLD*E ENDIF IF (ABS(D).GE.TOL1) THEN U=X+D ELSE U=X+SIGN(TOL1,D) ENDIF c find the optimum tj for the current guess of rhoj and evaluate the c objective function rhoj=U c write (*,*) ' CONFD--it,tj input to CONFT: ',iter,tj call CONFT (j,amix,rhoj,tj,ierr,herr) c write (*,*) ' CONFD--ierr,tj from CONFT: ',ierr,tj tauj=tredj/tj delj=rhoj/Dredj aj=PHIK(j,0,0,tauj,delj) dela=amix-aj Zj=1.0d0+delj*PHIK(j,0,1,tauj,delj) delZ=Zmix-Zj c call DPDDK (j,tj,rhoj,dpdrho) !debug only c write (*,1089) j,ITER,tj,rhoj,amix,aj,Zmix,Zj,dpdrho c1089 format (' CONFD--j,it,tj,rhoj,amix,aj,Zmix,Zj,dP/dD: ',2i3,f9.3, c & f10.5,2e14.6,2f9.5,e14.6) FU=dela**2+delZ**2 c IF (FU.LE.FX) THEN IF (U.GE.X) THEN A=X ELSE B=X ENDIF V=W FV=FW W=X FW=FX X=U FX=FU ELSE IF (U.LT.X) THEN A=U ELSE B=U ENDIF IF ((FU.LE.FW).or.(ABS(W-X).LT.1.0D-12)) THEN V=W FV=FW W=U FW=FU ELSE IF ((FU.LE.FV) .or. (ABS(V-W).LT.1.0D-12) .or. & (ABS(V-X).LT.1.0D-12)) THEN V=U FV=FU ENDIF ENDIF ENDIF 100 CONTINUE IERR=-58 herr='[CONFD error 58] ECS-transport routines did not converge'// & hnull c write (*,1099) j,amix,Zmix,tj,rhoj,FU c1099 format (' CONFD--j,A,Z,tj,rhoj,resid: ',i3,e14.6,f9.5,f9.3,f10.5, c & e14.6) RETURN 99 CONTINUE rhoj=X !value which minimizes objective function c write (*,1099) j,amix,Zmix,tj,rhoj,FU RETURN end !subroutine CONFD c c ====================================================================== c subroutine CONFT (k,amix,rhok,tk,ierr,herr) c c Find the conformal temperature for component k (including, possibly, c the reference fluid) which gives a reduced residual Helmholtz energy c which matches the input value. The zero of the function c FX = [a(tk,rhok) - amix] c is found using a damped secant method combined with a reguli-falsi. c This routine is called from within CONFTD, and the value of rhok is c the current guess value from that routine. c c inputs: c k--component number c k = 0: find conformal conditions for reference fluid c k > 0: find conformal conditions for component k c amix--reduced residual Helmholtz energy (A - A*)/RT [-] c rhok--conformal density for component k [mol/L] c tk--initial guess for conformal temperature for component k [K] c outputs: c tk--converged conformal temperature for component k [K] c ierr--error flag: 0 = successful c -58 = did not converge c herr--error string (character*255 variable if ierr<>0) c c written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO c 09-26-97 MM, original version, replaces MINF by S.A. Klein c 01-13-98 MM, set damping ratio to 0.8 and constrain new guesses c implicit double precision (a-h,o-z) implicit integer (i-n) character*1 htab,hnull character*255 herr logical lxneg,lxpos dimension x(3),fx(2) common /CHAR/ htab,hnull data itmax /20/ data tolr /1.0d-7/ !convergence tolerance data xdamp /0.8d0/ !damping (or acceleration) ratio, normal = 1 c ierr=0 herr=hnull c find the reducing ratios for component k call REDK (k,tredk,Dredk) delk=rhok/Dredk !density is fixed for this iteration c c begin iteration--the basic iteration is a secant method, but c positive and negative roots are stored, allowing use of slower, but c more robust, reguli-falsi method in event that guess by secant method c is further from solution c lxpos=.false. !initialize flags for reguli-falsi iteration lxneg=.false. xneg=0.0d0 xpos=0.0d0 fxneg=0.0d0 fxpos=0.0d0 j=1 !j=1 for first iteration, j=2 for all others x(1)=tk !first guess for independent variable do 800 it=1,itmax c evaluate objective function tauk=tredk/x(j) !x(j) is current guess for tk fx(j)=amix-PHIK(k,0,0,tauk,delk) c check for convergence if (ABS(fx(j)).lt.tolr) then c iteration has converged goto 840 else if (fx(j).lt.0.0d0) then lxneg=.true. xneg=x(j) fxneg=fx(j) else lxpos=.true. xpos=x(j) fxpos=fx(j) end if end if c c compute new guess for independent variable if (j.eq.1) then c for first iteration, new guess is ratio of old x(2)=x(1)*1.005 j=2 else c subsequent iterations--use (damped) secant method if (ABS(fx(2)-fx(1)).lt.1.0d-12) then c avoid division by zero, use average of previous guesses; if fx(1) was c equal to fx(2) by coincidence, this should allow solution; if there c is a more severe problem iteration will probably not converge, but c neither will it blow up x(3)=0.5d0*(x(1)+x(2)) else x(3)=x(2)-xdamp*fx(2)*(x(2)-x(1))/(fx(2)-fx(1)) end if if (ABS(x(3)/x(2)-1.0d0).ge.0.6) then c do not allow too great a change in tk, a.k.a. x c write (*,*) ' CONFD--unconstrained tk from CONFT: ',x(3) x(3)=x(2)*(1.0d0+SIGN(0.6d0,x(3)/x(2)-1.0d0)) end if c write (*,*) ' CONFT--old, new tk: ',x(2),x(3) c check that new guess is not outside bounds, if so use reguli-falsi c (provided that bounds on root have been found) if ((x(3).gt.MAX(xpos,xneg).or. x(3).lt.MIN(xpos,xneg)) & .and. lxneg .and. lxpos) then x(3)=xpos-fxpos*(xpos-xneg)/(fxpos-fxneg) end if c discard oldest iteration x(1)=x(2) x(2)=x(3) fx(1)=fx(2) end if 800 continue c c iteration loop has not converged c write (*,1801) k,x(1),x(2),fx(1) c1801 format (' CONFT--secant method for t_conf did not converge; ', c & 'k,tk1,tk2,ft1: ',i3,2f10.3,e14.6) c 840 continue tk=x(j) c RETURN end !subroutine CONFT c c ====================================================================== c function CBRT (x) c c cube root function--allows negative arguments c c input: c x--value to be acted upon c output (as function value): c CBRT--cube root of x c c written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO c 08-25-97 MM, original version c implicit double precision (a-h,o-z) implicit integer (i-n) c if (ABS(x).ge.0.0d0) then CBRT=x**(1.0d0/3.0d0) else CBRT=-x**(1.0d0/3.0d0) end if c RETURN end !function CBRT c c c 1 2 3 4 5 6 7 c23456789012345678901234567890123456789012345678901234567890123456789012 c c ====================================================================== c end file trns_ECS.f c ======================================================================