c begin file trnsp.f c c This file contains the transport property routines c c contained here are: c subroutine TRNPRP (t,rho,x,eta,tcx,ierr,herr) c subroutine ETAK (icomp,t,rho,eta,ierr,herr) c subroutine ETAK0 (icomp,t,eta0,ierr,herr) c subroutine ETAK1 (icomp,t,eta1,ierr,herr) c subroutine ETAKR (icomp,t,rho,etar,ierr,herr) c subroutine ETAKB (icomp,t,rho,etab,ierr,herr) c subroutine TCXK (icomp,t,rho,tcx,ierr,herr) c subroutine TCXK0 (icomp,t,tcx0,ierr,herr) c subroutine TCXKB (icomp,t,rho,tcxb,ierr,herr) c subroutine TCXKC (icomp,t,rho,tcxc,ierr,herr) c function OMEGA (icomp,t,epsk,hmodci) c subroutine SETCI1 (nread,icomp,hcasn,ierr,herr) c function OMEGA1 (icomp,t,epsk) c subroutine SETCI2 (nread,icomp,hcasn,ierr,herr) c function OMEGA2 (icomp,t,epsk) c c ===================================================================== c ===================================================================== c subroutine TRNPRP (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 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 -31 = temperature out of range for conductivity c -32 = density out of range for conductivity c -33 = T and D out of range for conductivity c -41 = temperature out of range for viscosity c -42 = density out of range for viscosity c -43 = T and D out of range for viscosity c -51 = T out of range for both visc and t.c. c -52 = D out of range for both visc and t.c. c -53 = T and/or D out of range for visc and t.c. c 39 = model not found for thermal conductivity c 49 = model not found for viscosity c 50 = ammonia/water mixture (no properties calculated) c 51 = exactly at t, rhoc for a pure fluid; k is infinite c -58,-59 = ECS model 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 10-15-96 MM, original version, supercedes older routine of same name c 02-21-97 MM, check input T,rho against limits for pures, c generate warning messages c 08-20-97 MM, initialize visc and t.c. to flag indicating not calculated c 10-01-97 MM, add compiler switch to allow access by DLL c 02-06-98 MM, move ierret=0 so it is always initialized c 08-14-98 MM, x array never dimensioned!? c 01-25-00 EWL, switched "eta, tcxecs" to "etaecs, tcx" in one of the ECS calls c 05-24-00 EWL, if no transport routines are loaded, do not calculate values. c this is important for D2, Fl, and NF3 which have no exp. data c 11-01-01 EWL, black out ammonia/water mixtures c 11-02-01 MLH, block out computation exactly at tc, rhoc for pure fluid k (its infinite) c 09-08-04 MLH, names shortened to 6 characters c 01-04-07 EWL, add check for mixture component with no equation c 01-10-07 EWL, check for mixtures with water c 01-18-07 MLH, do not switch to ECS if out of range of correlation c 03-12-07 EWL, keep eta and tcx equal to xnota if no correlation instead of switching to xnotc at the very bottom c 03-12-07 EWL, check for neither tcx or eta models at begining and exit if so c implicit double precision (a-h,o-z) implicit integer (i-n) c cDEC$ ATTRIBUTES DLLEXPORT :: TRNPRP c dll_export TRNPRP c parameter (ncmax=20) !max number of components in mixture parameter (nrefmx=10) !max number of fluids for transport ECS parameter (n0=-ncmax-nrefmx,nx=ncmax) parameter (nrf0=n0) !lower limit for transport ref fluid arrays character*1 htab,hnull character*3 hetamx,heta,htcxmx,htcx character*255 herr,herrl,herrec,herrpf logical lecs dimension x(ncmax) c common /HCHAR/ htab,hnull common /NCOMP/ nc,ic common /TRNMOD/ hetamx,heta(nrf0:ncmax),htcxmx,htcx(nrf0:ncmax) c common block containing flags to GUI (initialized in BDSET in setup.f) common /FLAGS/ xnota,x2ph,xsubc,xsuph,xsupc,xinf,x7,xnotd,xnotc common /FLAGS2/ iamwat,ianc,iwat common /CCON/ tcp(n0:nx),pcp(n0:nx),rhocp(n0:nx),Zcritp(n0:nx), & ttpp(n0:nx),ptpp(n0:nx),dtpp(n0:nx),dtpvp(n0:nx), & tnbpp(n0:nx),dnbplp(n0:nx),dnbpvp(n0:nx), & wmp(n0:nx),accenp(n0:nx),dipole(n0:nx),Reosp(n0:nx) c ierr=0 herr=' ' do i=1,nc if ((heta(i).eq.'NUL' .or. heta(i).eq.' ') & .and.(htcx(i).eq.'NUL' .or. htcx(i).eq.' ')) then eta=xnota tcx=xnota RETURN endif enddo eta=0.d0 tcx=0.d0 lecs=.false. if (iamwat.eq.1) then ierr=50 write (herr,1001) ierr,hnull call ERRMSG (ierr,herr) 1001 format ('[TRNPRP error',i3,'] transport equations are ', & 'not available for the ammonia/water mixture',a1) RETURN end if if (iwat.eq.1) then ierr=50 write (herr,1005) ierr,hnull call ERRMSG (ierr,herr) 1005 format ('[TRNPRP error',i3,'] transport equations are ', & 'not available for mixtures with water',a1) RETURN end if 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 ierret=0 !initialize if (nc.eq.1 .or. ic.ne.0) then icomp=1 if (ic.ne.0) icomp=ic c no model was specified, so do not calculate any properties if (heta(icomp)(1:1).eq.' ' .or. heta(icomp).eq.'NUL') then eta=xnota c special case--pure fluid elseif (heta(icomp)(1:2).eq.'EC') then c call ECS subroutines c write (*,*) ' TRNPRP--about to call TRNECS for viscosity' call TRNECS (t,rho,x,eta,tcxecs,ierrec,herr) ierr=ierrec lecs=.true. !set flag (no need to call TRNECS again for t.c.) else c pure fluid correlation pxx=0.0d0 call LIMITX ('ETA',t,rho,pxx,x,tmin,tmax,Dmax,pmax,ierr,herrl) c write (*,1002) tmin,tmax,Dmax,ierr c1002 format (1x,' TRNPRP--T,rho limits for visc: ',2f8.2,f10.4, c & '; ierr = ',i3) if (ierr.ge.1) then c if out of range, call ECS model and generate warning message ierr=-40-ierr write (herr,1040) ierr,herrl(39:168),hnull call ERRMSG (ierr,herr) ierret=ierr c do NOT call ecs, just give warning c call TRNECS (t,rho,x,eta,tcxecs,ierrec,herrec) c lecs=.true. if (ierrec.ne.0) then c if ECS method results in error, pass that as output, otherwise retain c the error from LIMITX write (herr,1004) ierr,herrec(1:146),hnull call ERRMSG (ierr,herr) 1004 format ('[TRNPRP warning',i4,'] pure fluid correlation ', & 'out of range; attempted to use ECS method, it ', & 'generated error: ',a146,a1) ierr=ierrec end if else c call pure fluid model lecs=.false. call ETAK (icomp,t,rho,eta,ierr,herr) end if end if c c no model was specified, so do not calculate any properties if (htcx(icomp)(1:1).eq.' ' .or. htcx(icomp).eq.'NUL') then tcx=xnota c call ECS subroutines elseif (htcx(icomp)(1:2).eq.'EC') then if (lecs) then tcx=tcxecs !use value computed on call above else call TRNECS (t,rho,x,etaecs,tcx,ierrec,herrec) if (ierrec.ne.0) then c if ECS method results in error, pass that as output, otherwise retain c the error from LIMITX write (herr,1004) ierr,herrec(1:146),hnull call ERRMSG (ierr,herr) ierr=ierrec end if end if else c pure fluid correlation pxx=0.0d0 call LIMITX ('TCX',t,rho,pxx,x,tmin,tmx,Dmx,pmx,ierrl,herrl) c write (*,1003) tmin,tmx,Dmx,ierrl c1003 format (1x,' TRNPRP--T,rho limits for t.c.: ',2f8.2,f10.4, c & '; ierr = ',i3) if (ierrl.ge.1) then c out of range of thermal conductivity correlation, call ECS model if (ierret.ne.0) then c both thermal conductivity and viscosity are out of range if (ierret.eq.-41) then ierr=-50-ierrl write (herr,1050) ierr,herrl(39:158),hnull call ERRMSG (ierr,herr) else if (ierret.eq.-42) then if (ierrl.eq.2) then ierr=-52 write (herr,1052) hnull call ERRMSG (ierr,herr) else ierr=-53 write (herr,1053) hnull call ERRMSG (ierr,herr) end if else ierr=-53 write (herr,1053) hnull call ERRMSG (ierr,herr) end if else c just thermal conductivity correlation is out of range ierr=-30-ierrl write (herr,1030) ierr,herrl(39:168),hnull call ERRMSG (ierr,herr) end if if (lecs) then tcx=tcxecs !use value computed on call above else c do NOT switch to ECS c call TRNECS (t,rho,x,etaecs,tcx,ierrec,herrec) if (ierrec.ne.0) then c if ECS method results in error, pass that as output, otherwise retain c the error from LIMITX write (herr,1004) ierr,herrec(1:146),hnull call ERRMSG (ierr,herr) ierr=ierrec end if end if else c call pure fluid model call TCXK (icomp,t,rho,tcx,ierrpf,herrpf) if (ierrpf.ne.0) then c if pure fluid model results in error, pass that as output, otherwise c retain the error from LIMITX ierr=ierrpf herr=herrpf end if end if end if c check to see if thermal conductivity is requested at exactly the crit point IF(t.eq.tcp(icomp) .AND. rho.eq.rhocp(icomp)) then ierr=51 write (herr,1002) ierr,hnull call ERRMSG (ierr,herr) 1002 format ('[TRNPRP error',i3,'] pure fluid thermal ', & 'conductivity is infinite at exactly Tc, rhoc',a1) tcx=xnotc end if else c c general (mixture) case c if (hetamx(1:2).eq.'EC') then c call ECS subroutines call TRNECS (t,rho,x,eta,tcxecs,ierr,herr) lecs=.true. !set flag (no need to call TRNECS again for t.c.) else lecs=.false. c [insert calls to additional mixture models here] end if c if (htcxmx(1:2).eq.'EC') then c call ECS subroutines if (lecs) then tcx=tcxecs !use value computed on call above else call TRNECS (t,rho,x,etaecs,tcx,ierrec,herrec) if (ierrec.ne.0) then c if ECS method results in error, pass that as output, otherwise retain c any possible error from call to viscosity routines ierr=ierrec herr=herrec end if end if else lecs=.false. c [insert calls to additional mixture models here] end if end if if (eta.lt.0 .and. eta.ne.xnota) eta=xnotc if (tcx.lt.0 .and. tcx.ne.xnota) tcx=xnotc c do i=1,nc if (heta(i).eq.'NUL' .or. heta(i).eq.' ') eta=xnota if (htcx(i).eq.'NUL' .or. htcx(i).eq.' ') tcx=xnota enddo RETURN c 1030 format ('[TRNPRP warning',i4,'] one or more inputs to the ', & 'thermal conductivity correlation are out of range; ', & 'will use the ECS model; ',a130,a1) 1040 format ('[TRNPRP warning',i4,'] one or more inputs to the ', & 'viscosity correlation are out of range; ', & 'will use the ECS model; ',a130,a1) 1050 format ('[TRNPRP warning',i4,'] one or more inputs to the ', & 'thermal conductivity and viscosity correlations are ', & 'out of range; will use ECS model; ',a120,a1) 1052 format ('[TRNPRP warning -52] the density input to the ', & 'thermal conductivity and viscosity correlations are ', & 'out of range; will use the ECS model. ',a1) 1053 format ('[TRNPRP warning -53] temperature and/or density input ', & 'to the thermal conductivity and/or viscosity ', & 'correlations are out of range; will use the ECS model.', & a1) end !subroutine TRNPRP c c ====================================================================== c subroutine ETAK (icomp,t,rho,eta,ierr,herr) c c viscosity of a pure fluid as a function of temperature and density c c inputs: c icomp--component number in mixture (1..nc); 1 for pure fluid c t--temperature [K] c rho--molar density [mol/L] c outputs: c eta--viscosity (uPa.s) c ierr--error flag: 0 = successful c 1 = did not converge c 49 = unknown model specified c herr--error string (character*255 variable if ierr<>0) c c written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO c 10-15-96 MM, original version c 10-30-96 MM, call intermediate ETAK0, etc routines rather than c core-level routines directly c 01-17-97 MM, change initial-density dependence to visc virial coeff c implicit double precision (a-h,o-z) implicit integer (i-n) parameter (ncmax=20) !max number of components in mixture parameter (nrefmx=10) !max number of fluids for transport ECS parameter (n0=-ncmax-nrefmx) parameter (nrf0=n0) !lower limit for transport ref fluid arrays character*1 htab,hnull character*3 hetamx,heta,htcxmx,htcx character*255 herr c common /HCHAR/ htab,hnull common /TRNMOD/ hetamx,heta(nrf0:ncmax),htcxmx,htcx(nrf0:ncmax) c ierr=0 herr=' ' c c the viscosity is the sum of the dilute-gas, initial-density, and c residual terms (separate functions for the three terms are used c because of the requirements of certain mixture models, e.g. ECS) c call ETAK0(icomp,t,eta0,ierr,herr) call ETAK1(icomp,t,etaB2,ierr,herr) call ETAKR(icomp,t,rho,etar,ierr,herr) eta=eta0*(1.0d0+etaB2*rho)+etar c write (*,1000) icomp,t,rho,eta0,eta0*etaB2*rho,etar,eta c1000 format (1x,' ETAK--icomp,t,rho,eta_0,B2,resid,sum: ',i3,6f12.4) c RETURN end !subroutine ETAK c c ====================================================================== c subroutine ETAK0 (icomp,t,eta0,ierr,herr) c c dilute-gas contribution to the viscosity of a pure fluid as a c function of temperature and density c c inputs: c icomp--component number in mixture (1..nc); 1 for pure fluid c t--temperature [K] c outputs: c eta0--viscosity [uPa.s] c ierr--error flag: 0 = successful c 1 = did not converge c 49 = unknown model specified c herr--error string (character*255 variable if ierr<>0) c c written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO c 10-30-96 MM, original version c 03-07-00 EWL, add VS0 model c 03-05-02 EWL, initialize eta0 c 03-21-02 MLH, modified VS0 model to have dilute gas piece c 12-26-06 NLH, added vs4 model (generalized friction theory) c implicit double precision (a-h,o-z) implicit integer (i-n) parameter (ncmax=20) !max number of components in mixture parameter (nrefmx=10) !max number of fluids for transport ECS parameter (n0=-ncmax-nrefmx) parameter (nrf0=n0) !lower limit for transport ref fluid arrays character*1 htab,hnull character*3 hetamx,heta,htcxmx,htcx character*255 herr c common /HCHAR/ htab,hnull common /TRNMOD/ hetamx,heta(nrf0:ncmax),htcxmx,htcx(nrf0:ncmax) c ierr=0 eta0=0.d0 herr=' ' c c call appropriate core viscosity model if (heta(icomp)(1:2).eq.'EC') then eta0=ETA0DG(icomp,t) else if (heta(icomp).eq.'VS0') then eta0=ETA0DG(icomp,t) else if (heta(icomp).eq.'VS1') then eta0=ETA1DG(icomp,t) else if (heta(icomp).eq.'VS2') then eta0=ETA2DG(icomp,t) else if (heta(icomp).eq.'VS3') then eta0=ETA3DG(icomp,t) else if (heta(icomp).eq.'VS4') then eta0=ETA1DG(icomp,t) else if (heta(icomp).eq.'VS6') then c eta0=ETA6DG(icomp,t) eta0=0.0d0 !temporary ierr=49 herr='[ETAK0 error 49] call to unimplemented model (VS6)' call ERRMSG (ierr,herr) else ierr=49 herr='[ETAK0 error 49] unknown viscosity model specified'//hnull call ERRMSG (ierr,herr) end if c RETURN end !subroutine ETAK0 c c ====================================================================== c subroutine ETAK1 (icomp,t,etaB2,ierr,herr) c c the second viscosity virial coefficient (initial-density contribution) c of a pure fluid as a function of temperature c c inputs: c icomp--component number in mixture (1..nc); 1 for pure fluid c t--temperature [K] c outputs: c etaB2--second viscosity virial coefficient [L/(mol-uPa-s)] c ierr--error flag: 0 = successful c 1 = did not converge c 49 = unknown model specified c herr--error string (character*255 variable if ierr<>0) c c written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO c 10-30-96 MM, original version c 01-21-97 MM, call ETA1B2 for model VS1 c 03-07-00 EWL, add VS0 model c 03-06-02 EWL, initialize etaB2 c 12-26-06 MLH, added vs4 model c implicit double precision (a-h,o-z) implicit integer (i-n) parameter (ncmax=20) !max number of components in mixture parameter (nrefmx=10) !max number of fluids for transport ECS parameter (n0=-ncmax-nrefmx) parameter (nrf0=n0) !lower limit for transport ref fluid arrays character*1 htab,hnull character*3 hetamx,heta,htcxmx,htcx character*255 herr c common /HCHAR/ htab,hnull common /TRNMOD/ hetamx,heta(nrf0:ncmax),htcxmx,htcx(nrf0:ncmax) c ierr=0 etaB2=0.0d0 herr=' ' c c call appropriate core viscosity model c note that only the VS1 pure-fluid model presently implements an c explicit initial-density term if (heta(icomp).eq.'VS0') then etaB2=0.0d0 else if (heta(icomp).eq.'VS1') then etaB2=ETA1B2(icomp,t) else if (heta(icomp).eq.'VS2') then etaB2=0.0d0 !no initial-density term for Y+E else if (heta(icomp).eq.'VS3') then etaB2=0.0d0 else if (heta(icomp).eq.'VS4') then etaB2=0.0d0 else if (heta(icomp).eq.'VS6') then c etaB2=ETA6B2(icomp,t) etaB2=0.0d0 !temporary else ierr=49 herr='[ETAK1 error 49] unknown viscosity model specified'//hnull call ERRMSG (ierr,herr) end if c RETURN end !subroutine ETAK1 c c ====================================================================== c subroutine ETAKR (icomp,t,rho,etar,ierr,herr) c c residual contribution to the viscosity of a pure fluid as a function c of temperature and density c c inputs: c icomp--component number in mixture (1..nc); 1 for pure fluid c t--temperature [K] c rho--molar density [mol/L] c outputs: c etar--viscosity [uPa.s] c ierr--error flag: 0 = successful c 1 = did not converge c 49 = unknown model specified c herr--error string (character*255 variable if ierr<>0) c c written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO c 10-30-96 MM, original version c 03-07-00 EWL, add VS0 model c 03-05-02 EWL, initialize etar c 03-21-02 MLH, modified VS0 model to have residual piece c 12-26-06 MLH, added vs4 model c implicit double precision (a-h,o-z) implicit integer (i-n) parameter (ncmax=20) !max number of components in mixture parameter (nrefmx=10) !max number of fluids for transport ECS parameter (n0=-ncmax-nrefmx) parameter (nrf0=n0) !lower limit for transport ref fluid arrays character*1 htab,hnull character*3 hetamx,heta,htcxmx,htcx character*255 herr c common /HCHAR/ htab,hnull common /TRNMOD/ hetamx,heta(nrf0:ncmax),htcxmx,htcx(nrf0:ncmax) c ierr=0 etar=0.d0 herr=' ' c c call appropriate core viscosity model if (heta(icomp).eq.'VS0') then c eta0hc provides full viscosities; subtract dilute to get resid etar=ETA0HC(icomp,t,rho,ierr,herr)-ETA0DG(icomp,t) else if (heta(icomp).eq.'VS1') then etar=ETA1RS(icomp,t,rho) else if (heta(icomp).eq.'VS2') then etar=ETA2RS(icomp,t,rho) else if (heta(icomp).eq.'VS3') then etar=ETA3RS(icomp,t,rho) else if (heta(icomp).eq.'VS4') then etar=ETA4RS(icomp,t,rho) else if (heta(icomp).eq.'VS6') then c etar=ETA6RS(icomp,t,rho) etar=0.0d0 !temporary ierr=49 herr='[ETAKR error 49] call to unimplemented model (VS6)' call ERRMSG (ierr,herr) else ierr=49 herr='[ETAKR error 49] unknown viscosity model specified'//hnull call ERRMSG (ierr,herr) end if c RETURN end !subroutine ETAKR c c ====================================================================== c subroutine ETAKB (icomp,t,rho,etab,ierr,herr) c c background contribution to the viscosity of a pure fluid as a c function of temperature and density; this includes both the residual c term and any initial density (second virial) term c c inputs: c icomp--component number in mixture (1..nc); 1 for pure fluid c t--temperature [K] c rho--molar density [mol/L] c outputs: c etab--background viscosity [uPa.s] c ierr--error flag: 0 = successful c <>0 = error in underlying routine, passed up c herr--error string (character*255 variable if ierr<>0) c c written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO c 02-04-97 MM, original version c implicit double precision (a-h,o-z) implicit integer (i-n) parameter (ncmax=20) !max number of components in mixture parameter (nrefmx=10) !max number of fluids for transport ECS parameter (n0=-ncmax-nrefmx) parameter (nrf0=n0) !lower limit for transport ref fluid arrays character*1 htab,hnull character*3 hetamx,heta,htcxmx,htcx character*255 herr c common /HCHAR/ htab,hnull common /TRNMOD/ hetamx,heta(nrf0:ncmax),htcxmx,htcx(nrf0:ncmax) c ierr=0 herr=' ' c c the viscosity is the sum of the dilute-gas, initial-density, and c residual terms (separate functions for the three terms are used c because of the requirements of certain mixture models, e.g. ECS); c this routine returns the sum of the initial-density and residual terms c call ETAK0(icomp,t,eta0,ierr,herr) call ETAK1(icomp,t,etaB2,ierr,herr) call ETAKR(icomp,t,rho,etar,ierr,herr) etab=eta0*etaB2*rho+etar c write (*,1000) icomp,t,rho,eta0,eta0*etaB2*rho,etar,etab c1000 format (1x,' ETAKB--icomp,t,rho,etao,B2,resid,sum: ',i3,6f12.4) c RETURN end !subroutine ETAKB c c ====================================================================== c subroutine TCXK (icomp,t,rho,tcx,ierr,herr) c c thermal conductivity of a pure fluid as a function of temp and density c c inputs: c icomp--component number in mixture (1..nc); 1 for pure fluid c t--temperature [K] c rho--molar density [mol/L] c outputs: c tcx--thermal conductivity [W/m-K] c ierr--error flag: 0 = successful c 1 = did not converge c 39 = unknown model specified c herr--error string (character*255 variable if ierr<>0) c c written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO c 10-15-96 MM, original version c 10-30-96 MM, call intermediate TCXK0, etc routines rather than c core-level routines directly c implicit double precision (a-h,o-z) implicit integer (i-n) parameter (ncmax=20) !max number of components in mixture parameter (nrefmx=10) !max number of fluids for transport ECS parameter (n0=-ncmax-nrefmx) parameter (nrf0=n0) !lower limit for transport ref fluid arrays character*1 htab,hnull character*3 hetamx,heta,htcxmx,htcx character*255 herr c common /HCHAR/ htab,hnull common /TRNMOD/ hetamx,heta(nrf0:ncmax),htcxmx,htcx(nrf0:ncmax) c ierr=0 herr=' ' c c the thermal conductivity is the sum of the dilute-gas, background c (or residual), and critical enhancement terms (separate functions c for the three terms are used because of the requirements of certain c mixture models, e.g. ECS) c call TCXK0(icomp,t,tcx0,ierr,herr) call TCXKB(icomp,t,rho,tcxb,ierr,herr) call TCXKC(icomp,t,rho,tcxc,ierr,herr) tcx=tcx0+tcxb+tcxc c write (*,1000) icomp,t,rho,tcx0,tcxb,tcxc,tcx c1000 format (1x,' TCXK--icomp,t,rho,tcx_0,bk,crit,sum: ',i3,6f12.6) c RETURN end !subroutine TCXK c c ====================================================================== c subroutine TCXK0 (icomp,t,tcx0,ierr,herr) c c dilute-gas part of the thermal conductivity of a pure fluid as a c function of temp and density c c inputs: c icomp--component number in mixture (1..nc); 1 for pure fluid c t--temperature [K] c outputs: c tcx0--thermal conductivity [W/m-K] c ierr--error flag: 0 = successful c 1 = did not converge c 39 = unknown model specified c herr--error string (character*255 variable if ierr<>0) c c written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO c 10-30-96 MM, original version c 03-07-00 EWL, add TC0 model c implicit double precision (a-h,o-z) implicit integer (i-n) parameter (ncmax=20) !max number of components in mixture parameter (nrefmx=10) !max number of fluids for transport ECS parameter (n0=-ncmax-nrefmx) parameter (nrf0=n0) !lower limit for transport ref fluid arrays character*1 htab,hnull character*3 hetamx,heta,htcxmx,htcx character*255 herr c common /HCHAR/ htab,hnull common /TRNMOD/ hetamx,heta(nrf0:ncmax),htcxmx,htcx(nrf0:ncmax) c ierr=0 herr=' ' c c call appropriate core thermal conductivity model if (htcx(icomp).eq.'TC0') then tcx0=0.0d0 else if (htcx(icomp).eq.'TC1') then tcx0=TCX1DG(icomp,t) else if (htcx(icomp).eq.'TC2') then tcx0=TCX2DG(icomp,t) else if (htcx(icomp).eq.'TC3') then tcx0=TCX3DG(icomp,t) else if (htcx(icomp).eq.'TC6') then c tcx0=TCX6DG(icomp,t) tcx0=0.0d0 !temporary ierr=39 herr='[TCXK0 error 39] call to unimplemented model (TC6)' call ERRMSG (ierr,herr) else ierr=39 herr='[TCXK0 error 39] unknown thermal conductivity model'// & ' specified'//hnull call ERRMSG (ierr,herr) end if c RETURN end !subroutine TCXK0 c c ====================================================================== c subroutine TCXKB (icomp,t,rho,tcxb,ierr,herr) c c background part of the thermal conductivity of a pure fluid as a c function of temp and density c c inputs: c icomp--component number in mixture (1..nc); 1 for pure fluid c t--temperature [K] c rho--molar density [mol/L] c outputs: c tcxb--thermal conductivity [W/m-K] c ierr--error flag: 0 = successful c 1 = did not converge c 39 = unknown model specified c herr--error string (character*255 variable if ierr<>0) c c written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO c 10-30-96 MM, original version c 03-07-00 EWL, add TC0 model c implicit double precision (a-h,o-z) implicit integer (i-n) parameter (ncmax=20) !max number of components in mixture parameter (nrefmx=10) !max number of fluids for transport ECS parameter (n0=-ncmax-nrefmx) parameter (nrf0=n0) !lower limit for transport ref fluid arrays character*1 htab,hnull character*3 hetamx,heta,htcxmx,htcx character*255 herr c common /HCHAR/ htab,hnull common /TRNMOD/ hetamx,heta(nrf0:ncmax),htcxmx,htcx(nrf0:ncmax) c ierr=0 herr=' ' c c call appropriate core thermal conductivity model if (htcx(icomp).eq.'TC0') then tcxb=TCX0HC(icomp,t,rho,ierr,herr) elseif (htcx(icomp).eq.'TC1') then tcxb=TCX1BK(icomp,t,rho) else if (htcx(icomp).eq.'TC2') then tcxb=TCX2BK(icomp,t,rho) else if (htcx(icomp).eq.'TC3') then tcxb=TCX3BK(icomp,t,rho) else if (htcx(icomp).eq.'TC6') then c tcxb=TCX6BK(icomp,t,rho) tcxb=0.0d0 !temporary ierr=39 herr='[TCXKB error 39] call to unimplemented model (TC6)' call ERRMSG (ierr,herr) else ierr=39 herr='[TCXKB error 39] unknown thermal conductivity model'// & ' specified'//hnull call ERRMSG (ierr,herr) end if c RETURN end !subroutine TCXKB c c ====================================================================== c subroutine TCXKC (icomp,t,rho,tcxc,ierr,herr) c c critical enhancement for the thermal conductivity of a pure fluid as c a function of temp and density c c inputs: c icomp--component number in mixture (1..nc); 1 for pure fluid c t--temperature [K] c rho--molar density [mol/L] c outputs: c tcxc--thermal conductivity [W/m-K] c ierr--error flag: 0 = successful c 1 = did not converge c 39 = unknown model specified c herr--error string (character*255 variable if ierr<>0) c c written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO c 10-30-96 MM, original version c 02-24-97 MM, add /CREMOD/ c 11-06-00 EWL, add water c 07-01-02 MLH, add TK6 c implicit double precision (a-h,o-z) implicit integer (i-n) parameter (ncmax=20) !max number of components in mixture parameter (nrefmx=10) !max number of fluids for transport ECS parameter (n0=-ncmax-nrefmx) parameter (nrf0=n0) !lower limit for transport ref fluid arrays character*1 htab,hnull character*3 hetamx,heta,htcxmx,htcx character*3 hetacr,htcxcr character*255 herr DIMENSION xfeed(ncmax) c common /HCHAR/ htab,hnull common /TRNMOD/ hetamx,heta(nrf0:ncmax),htcxmx,htcx(nrf0:ncmax) common /CREMOD/ hetacr(nrf0:ncmax),htcxcr(nrf0:ncmax) c ierr=0 herr=' ' c c call appropriate core thermal conductivity model c write (*,1100) icomp,htcxcr(icomp),t,rho c1100 format (1x,' TCXKC--icomp,htcxcr: ',i2,2x,a3,'; at t,rho = ', c & 2e14.6) if (htcxcr(icomp).eq.'TK1') then tcxc=TCX1CR(icomp,t,rho) else if (htcxcr(icomp).eq.'TK2') then c hydrocarbon model of Younglove & Ely (critical part integral w/ TC2) tcxc=TCX2CR(icomp,t,rho) else if (htcxcr(icomp).eq.'TK3') then c write (*,*) ' TCXKC--call critical model TK3' tcxc=TCX3CR(icomp,t,rho) c write (*,*) ' TCXKC--critical enhancement from TK3: ',tcxc else if (htcxcr(icomp).eq.'TK4') then tcxc=TCX4CR(icomp,t,rho) else if (htcxcr(icomp).eq.'TK6') then xfeed(icomp)=1.0d0 tcxc=TCXM1C(xfeed,t,rho,ierr,herr) c tcxc=0.0d0 !temporary c ierr=39 c herr='[TCXKC error 39] call to unimplemented model (TK6)' c call ERRMSG (ierr,herr) else if (htcxcr(icomp).eq.'NUL') then c no critical enhancement is used tcxc=0.0d0 else if (htcxcr(icomp).eq.'NH3') then c special function for ammonia; model of Tufeu tcxc=TCCNH3(icomp,t,rho) else if (htcxcr(icomp).eq.'CH4') then c special function for methane; model of Friend et al. tcxc=TCCCH4(icomp,t,rho) else ierr=39 herr='[TCXKC error 39] unknown thermal conductivity critical'// & ' enhancement model specified: '//htcxcr(icomp)//hnull call ERRMSG (ierr,herr) end if c RETURN end !subroutine TCXKC c c ====================================================================== c function OMEGA (icomp,t,epsk,hmodci) c c collision integral for a pure fluid as a function of temperature and c the Lennard-Jones epsilon/k parameter c c inputs: c icomp--component number in mixture (1..nc); 1 for pure fluid c t--temperature [K] c epsk--Lennard-Jones epsilon/k parameter [K] c output (as function value): c omega--the collision integral c c written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO c 10-16-96 MM, original version c 01-31-97 MM, add case of model = 'NUL' c 03-07-00 EWL, add CI0 model c implicit double precision (a-h,o-z) implicit integer (i-n) character*1 htab,hnull character*3 hmodci c common /HCHAR/ htab,hnull c c call appropriate function for collision integral omega=1.0d0 if (hmodci.eq.'CI1') then OMEGA=OMEGA1(icomp,t,epsk) else if (hmodci.eq.'CI2') then OMEGA=OMEGA2(icomp,t,epsk) else if (hmodci.eq.'CI0') then if (epsk.ne.0.d0) OMEGA=OMEGAS(2,2,t/epsk) else if (hmodci.eq.'NUL') then c collision integral not used; 'NUL' is place-holder in fluid file c should not ever get here, but return value to avoid error OMEGA=1.0d0 !return 1, as omega appears in denominator else c write (*,*) '[OMEGA error] unknown col integral: ',hmodci OMEGA=1.0d0 !return 1, as omega appears in denominator end if c RETURN end !function OMEGA c c ====================================================================== c subroutine SETCI1 (nread,icomp,hcasno,ierr,herr) c c initialize the collision integral; this is an empirical function in c log(T), the form used by several authors including Fenghour (1995) for c ammonia and Krauss (1996) for R152a; this same form is used for the c "reduced effective collision cross-section" used by Wilhelm & Vogel c (1995) and Laesecke (1997) for R134a 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 (0..nc) c 1 for pure fluid; 0 for ECS reference fluid c hcasno--CAS number of component icomp (not req'd if reading from file) c c outputs: c ierr--error flag: 0 = successful c 39 = error--model not implemented c herr--error string (character*255 variable if ierr<>0) c other quantities returned via arrays in commons c c written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO c 01-21-97 MM, original version (based on SETCI2) c 05-14-97 MM, change power of T coefficient to integer, add to /WIFOMG/ c 08-19-97 MM, change error number for nread<=0; input hcasno is not array c 12-02-97 MM, skip over limits on file read c 08-13-98 MM, delete obsolete (unused) format statement c implicit double precision (a-h,o-z) implicit integer (i-n) parameter (ncmax=20) !max number of components in mixture parameter (nrefmx=10) !max number of fluids for transport ECS parameter (n0=-ncmax-nrefmx,nx=ncmax) parameter (nrf0=n0) !lower limit for transport ref fluid arrays parameter (mxomg=15) !max no. coefficients for collision integral character*1 htab,hnull character*12 hcasno character*255 herr c common /HCHAR/ htab,hnull c commons storing the coefficients to the fit for collision integral common /WCFOM1/ comg(nrf0:nx,mxomg,2) common /WIFOM1/ ntomg(nrf0:nx),icomg(nrf0:nx,mxomg) c if (nread.le.0) then c get coefficients from block data--this option not implemented, c place holder to maintain parallel structure with EOS setup routines ierr=101 write (herr,1101) nread,hcasno,hnull call ERRMSG (ierr,herr) 1101 format ('[SETCI1 error 101] illegal file specified; nread = ', & i4,'; CAS no. = ',a12,a1) else c read data from file (should have been opened by SETUP) c write (*,*) ' SETCI1--read component',icomp,' from unit',nread c limits are not used with collision integral, but retain for compatibility read (nread,*) !tmin !lower temperature limit (dummy) read (nread,*) !tmax !upper temperature limit (dummy) read (nread,*) !pjunk !upper pressure limit (n/a) read (nread,*) !rhojnk !upper density limit (n/a) read (nread,*) nterm !number of terms ntomg(icomp)=nterm if (ntomg(icomp).ge.1) then do j=1,ntomg(icomp) c first is numerical coeff, second is power of T read (nread,*) comg(icomp,j,1),icomg(icomp,j) enddo end if c write (*,*) ' SETCI1--final coeff: ',comg(icomp,nterm,1) ierr=0 herr=' ' end if c RETURN end !subroutine SETCI1 c c ====================================================================== c function OMEGA1 (icomp,t,epsk) c c collision integral for a pure fluid as a function of temperature and c the Lennard-Jones epsilon/k parameter; this is an empirical function c of a form used by several authors including Fenghour (1995) for c ammonia, Krauss (1996) for R152a, Laesecke (1997) for R134a, and c Vesovic (1990) for carbon dioxide. c c inputs: c icomp--component number in mixture (1..nc); 1 for pure fluid c t--temperature [K] c epsk--Lennard-Jones energy (epsilon/k) parameter [K] c output (as function value): c omega--the collision integral c c written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO c 01-21-97 MM, original version (based on OMEGA2) c 05-14-97 MM, change power of T coefficient to integer, add to /WIFOMG/ c 01-19-00 EWL, check for tstar=0 and avoid 0**0 error message c implicit double precision (a-h,o-z) implicit integer (i-n) parameter (ncmax=20) !max number of components in mixture parameter (nrefmx=10) !max number of fluids for transport ECS parameter (n0=-ncmax-nrefmx,nx=ncmax) parameter (nrf0=n0) !lower limit for transport ref fluid arrays parameter (mxomg=15) !max no. coefficients for collision integral character*1 htab,hnull c common /HCHAR/ htab,hnull c commons storing the coefficients to the fit for collision integral common /WCFOM1/ comg(nrf0:nx,mxomg,2) common /WIFOM1/ ntomg(nrf0:nx),icomg(nrf0:nx,mxomg) c omega1=1.0d0 tstar=1.0d-20 omgsum=0.0d0 if (epsk.ne.0.0d0) tstar=LOG(t/epsk) if (ABS(tstar).lt.1.0d-20) tstar=1.0d-20 do n=1,ntomg(icomp) omgsum=omgsum+comg(icomp,n,1)*tstar**icomg(icomp,n) enddo OMEGA1=EXP(omgsum) RETURN end !function OMEGA1 c c ====================================================================== c subroutine SETCI2 (nread,icomp,hcasno,ierr,herr) c c initialize the collision integral function of Younglove and Ely, c JPCRD 16:577-798 (1987) 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 (0..nc) c 1 for pure fluid; 0 for ECS reference fluid c hcasno--CAS number of component icomp (not req'd if reading from file) c c outputs: c ierr--error flag: 0 = successful c 39 = error--model not implemented c herr--error string (character*255 variable if ierr<>0) c other quantities returned via arrays in commons c c written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO c 10-16-96 MM, original version c 01-21-97 MM, add /WIFOMG/ (store # terms) for parallel w/ OMEGA1 c 05-14-97 MM, add icomg to /WIFOMG/ (used in OMEGA1) c 08-19-97 MM, change error number for nread<=0 c 12-02-97 MM, skip over limits on file read c 08-13-98 MM, delete obsolete (unused) format statement c implicit double precision (a-h,o-z) implicit integer (i-n) parameter (ncmax=20) !max number of components in mixture parameter (nrefmx=10) !max number of fluids for transport ECS parameter (n0=-ncmax-nrefmx,nx=ncmax) parameter (nrf0=n0) !lower limit for transport ref fluid arrays parameter (mxomg=15) !max no. coefficients for collision integral character*1 htab,hnull character*12 hcasno character*255 herr c common /HCHAR/ htab,hnull c commons storing the coefficients to the fit for collision integral common /WCFOM2/ comg(nrf0:nx,mxomg,2) common /WIFOM2/ ntomg(nrf0:nx),icomg(nrf0:nx,mxomg) c if (nread.le.0) then c get coefficients from block data--this option not implemented, c place holder to maintain parallel structure with EOS setup routines ierr=101 write (herr,1101) nread,hcasno,hnull call ERRMSG (ierr,herr) 1101 format ('[SETCI2 error 101] illegal file specified; nread = ', & i4,'; CAS no. = ',a12,a1) else c read data from file (should have been opened by SETUP) c write (*,*) ' SETCI2--read component',icomp,' from unit',nread c limits are not used with collision integral, but retain for compatibility read (nread,*) !tmin !lower temperature limit (dummy) read (nread,*) !tmax !upper temperature limit (dummy) read (nread,*) !pjunk !upper pressure limit (n/a) read (nread,*) !rhojnk !upper density limit (n/a) read (nread,*) nterm !number of terms ntomg(icomp)=nterm if (nterm.ge.1) then do j=1,nterm read (nread,*) comg(icomp,j,1) c the second element in coeff array is not used for this model comg(icomp,j,2)=0.0d0 enddo end if c write (*,*) ' SETCI1--final coeff: ',comg(icomp,nterm,1) ierr=0 herr=' ' end if c RETURN end !subroutine SETCI2 c c ====================================================================== c function OMEGA2 (icomp,t,epsk) c c collision integral for a pure fluid as a function of temperature and c the Lennard-Jones epsilon/k parameter; based on Equation (20) in c Younglove and Ely, JPCRD 16:577-798 (1987). c c N.B. there is misprint in Younglove & Ely, the exponent c is ((4-n)/3) not ((n+2)/3) c c inputs: c icomp--component number in mixture (1..nc); 1 for pure fluid c t--temperature [K] c epsk--Lennard-Jones energy (epsilon/k) parameter [K] c output (as function value): c omega--the collision integral c c written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO c 10-16-96 MM, original version c 01-21-97 MM, add /WIFOMG/ to maintain parallel w/ OMEGA1 c 05-14-97 MM, add icomg to /WIFOMG/ (used in OMEGA1) c implicit double precision (a-h,o-z) implicit integer (i-n) parameter (ncmax=20) !max number of components in mixture parameter (nrefmx=10) !max number of fluids for transport ECS parameter (n0=-ncmax-nrefmx,nx=ncmax) parameter (nrf0=n0) !lower limit for transport ref fluid arrays parameter (mxomg=15) !max no. coefficients for collision integral character*1 htab,hnull c common /HCHAR/ htab,hnull c commons storing the coefficients to the fit for collision integral common /WCFOM2/ comg(nrf0:nx,mxomg,2) common /WIFOM2/ ntomg(nrf0:nx),icomg(nrf0:nx,mxomg) c ekt3=1.0d0 omega2=1.0d0 omgsum=0.0d0 if (t.ne.0.d0) ekt3=(epsk/t)**(1.0d0/3.0d0) do n=1,ntomg(icomp) omgsum=omgsum+comg(icomp,n,1)*ekt3**(4-n) enddo if (omgsum.ne.0.d0) OMEGA2=1.0d0/omgsum RETURN end !function OMEGA2 c c c 1 2 3 4 5 6 7 c23456789012345678901234567890123456789012345678901234567890123456789012 c c ====================================================================== c end file trnsp.f c ======================================================================