c begin file core_STN.f c c This file contains core routines for the surface tension. c c contained here are: c subroutine SETST1 (nread,icomp,hcasno,ierr,herr) c subroutine STN (t,rhol,xl,tcrit,pcrit,sigma,ierr,herr) c subroutine STNK (icomp,tau,sigma,ierr,herr) c subroutine CRITF (zeta,x,tcrit,pcrit,Dcrit,ierr,herr) c c ====================================================================== c ====================================================================== c subroutine SETST1 (nread,icomp,hcasno,ierr,herr) c c set up working arrays for use with "ST1" surface tension model: c c sigma = sum[sigma_k*tau**sigexp_k] c tau = 1 - t/tcrit c c Note: The critical temperature used is that of the current c equation of state. This may differ slightly from that used c in the orignial correlation of surface tension; this change c is necessary to give proper behavior of surface tension near c the critical point and to avoid possible numerical crashes. c c inputs: c nread--file to read data from (file should have already been c 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 required, it is here c to maintain parallel structure with SETBWR and SETFEQ) c c outputs: c ierr--error flag: 0 = successful c 1 = error (e.g. fluid not found) c herr--error string (character*255 variable if ierr<>0) c other quantities returned via arrays in commons c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 03-24-96 MM, original version (skeleton only) c 08-16-96 MM, add actual ST1 model c 08-19-97 MM, get rid of herr=herr (avoid warning); flag nread<=0 c 12-02-97 MM, skip over pressure and density limit and Tc on file read 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 (nsigk=3) !max number of terms in sigma summation character*1 htab,hnull character*3 hsten,hstenk character*12 hcasno character*255 herr common /CHAR/ htab,hnull common /STNMOD/ hsten,hstenk(n0:nx) common /WLMSTN/ tmin(n0:nx),tmax(n0:nx) common /WNTST1/ nterm(n0:nx) common /WCFST1/ sigmak(n0:nx,nsigk),sigexp(n0:nx,nsigk) c if (nread.le.0) then ierr=101 write (herr,1101) nread,hcasno,hnull 1101 format ('[SETST1 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 (*,*) ' SETSTN--read component',icomp,' from unit',nread read (nread,*) tmin(icomp) !lower temperature limit read (nread,*) tmax(icomp) !upper temperature limit c the pressure and density limit and the Tc are not presently used, c but are contained in the file for consistency and possible future use; c skip over them in reading the file read (nread,*) !pjunk !upper pressure limit (n/a) read (nread,*) !rhojnk !upper density limit (n/a) read (nread,*) nterm(icomp) read (nread,*) !Tcjunk !Tc in original fit (not used) do 100 k=1,nterm(icomp) read (nread,*) sigmak(icomp,k),sigexp(icomp,k) 100 continue c RETURN end !subroutine SETSTN c c ====================================================================== c subroutine STN (t,rhol,xl,sigma,ierr,herr) c c compute surface tension with appropriate core model c c inputs: c t--temperature [K] c rhol--molar density of liquid phase [mol/L] c xl--composition of liquid phase [array of mol frac] c output: c sigma--surface tension [N/m] 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 c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 03-27-96 MM, original version (skeleton only) c 08-16-96 MM, add actual ST1 model c 04-17-97 MM, add pcrit to argument, break pures to separate STNK c add mixture model (Holcomb's mod of Moldover & Rainwater) c 04-21-97 MM, delete critical par from arguments, add call to CRITF c 12-15-97 MM, return if error from CRITF, set sigma = "not calculated" c implicit double precision (a-h,o-z) implicit integer (i-n) c parameter (ncmax=5) !max number of components in mixture parameter (n0=-ncmax,nx=ncmax) parameter (nsigk=3) !max number of terms in sigma summation character*1 htab,hnull character*3 hsten,hstenk character*255 herr,herr2 common /NCOMP/ nc common /CHAR/ htab,hnull common /STNMOD/ hsten,hstenk(n0:nx) common /WLMSTN/ tmin(n0:nx),tmax(n0:nx) common /WNTST1/ nterm(n0:nx) common /WCFST1/ sigmak(n0:nx,nsigk),sigexp(n0:nx,nsigk) common /Gcnst/ R 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 common block containing flags to GUI (initialized in BDSET in setup.f) common /FLAGS/ xnota,x2ph,xsubc,xsuph,xsupc,xinf,x7,xnotd,xnotc dimension xl(ncmax),zeta(ncmax),f(ncmax),cx(ncmax),xcritf(ncmax) c ierr=0 herr=hnull c call CRITP (xl,tcrit,pcrit,Dcrit,ierr,herr2) if (ierr.ne.0) then c error condition--set outputs, issue warning, and return ierr=160 sigma=0.0d0 write (herr,1016) ierr,herr2(1:236),hnull 1016 format ('[STN error 160] ',a236,a1) call ERRMSG (ierr,herr) RETURN end if if (t.gt.tcrit) then ierr=121 write (herr,1121) t,tcrit,hnull 1121 format ('[STN error 121] ', & 'temperature input to surface tension routine is ', & 'greater than critical temperature; T =',1pe11.3, & ' K, Tcrit =',1pe11.3,' K.',a1) sigma=0.0d0 c write (*,*) ' STN--output sigma (ierr = 121): ',sigma RETURN end if if (nc.eq.1) then c special case--pure component tau=1.0d0-t/tcrit call STNK (1,tau,sigma,ierr,herr) else if (hsten.eq.'STX' .or. hsten.eq.'STM') then c mixture case--apply mixing rules to the cx(i) c compute fugacities and fugacity fraction call FGCTY (t,rhol,xl,f) fsum=0.0d0 do 200 i=1,nc fsum=fsum+f(i) 200 continue do 204 i=1,nc zeta(i)=f(i)/fsum xcritf(i)=xl(i) !initial guess for crit comp at same zeta 204 continue c write (*,1204) (zeta(i),i=1,nc) c1204 format (1x,' STN--zeta(i): ',5f14.7) c write (*,1205) (cx(i),i=1,nc) c1205 format (1x,' STN--cx(i): ',5e14.4) c base surface tensions on critical parameters at same composition, c not at same fugacity fraction tau=1.0d0-t/tcrit tau126=tau**1.26d0 c find critical parameters at same fugacity fraction call CRITF (zeta,xcritf,tcritf,pcritf,Dcritf,ierr,herr) if (ierr.gt.0) then c error in CRITF--solution not possible sigma=xnotc RETURN end if alpha=0.10d0 const=3.74d0**1.5d0*SQRT(R)*alpha*(1.0d0-alpha)*(2.0d0-alpha) do 220 i=1,nc c define effective parameter in function sigma = sigma0*tau**1.26 call STNK (i,tau,sigk,ierr,herr) sig0=sigk/tau126 c write (*,*) ' STN--icomp, effective sigma_0: ',icomp,sig0 cx(i)=sig0**1.50d0/(const*SQRT(tc(i))*pc(i)) 220 continue if (hsten.eq.'STM') then c use Moldover & Rainwater method cmix=0.0d0 do 240 i=1,nc cmix=cmix+zeta(i)*cx(i) 240 continue c write (*,*) ' STN--cmix by M-R: ',cmix else if (hsten.eq.'STX') then c use Holcomb's modification of Moldover & Rainwater method pcsum=0.0d0 do 300 i=1,nc pcsum=pcsum+zeta(i)**2*pc(i) 300 continue cmix=0.0d0 do 340 i=1,nc c sum i = j terms cmix=cmix+zeta(i)**2*cx(i) if (i.lt.nc) then do 336 j=i+1,nc c sum cross terms c cij=0.5d0*(pcritf-pcsum)*SQRT(cx(i)*cx(j)/(pc(i)*pc(j))) c & /(zeta(i)*zeta(j)) c cmix=cmix+2.0d0*zeta(i)*zeta(j)*cij !factor 2 from ij = ji c above lines reduce to cmix=cmix+(pcritf-pcsum)*SQRT(cx(i)*cx(j)/(pc(i)*pc(j))) 336 continue end if 340 continue c write (*,*) ' STN--cmix by Holcomb''s mod to M-R: ',cmix else if (hsten.eq.'jnk') then c use simple scaling to critical pressure line pcsum=0.0d0 do 400 i=1,nc pcsum=pcsum+zeta(i)*pc(i) 400 continue cmix=0.0d0 do 440 i=1,nc c sum i = j terms cmix=cmix+zeta(i)**2*pcritf/pcsum*cx(i) if (i.lt.nc) then do 436 j=i+1,nc c sum cross terms cmix=cmix+zeta(i)*zeta(j)*pcritf/pcsum*SQRT(cx(i)*cx(j)) 436 continue end if 440 continue c write (*,*) ' STN--cmix by simple scaling of Pc: ',cmix else c write (*,*) ' STN error--model not found' cmix=0.0d0 end if c c recover mixture sigma_0 parameter from cmix c this expression is based on critical parameters at same zeta sig0=(const*SQRT(tcritf)*pcritf*cmix)**(2.0d0/3.0d0) sigma=sig0*tau**1.26d0 else ierr=99 sigma=-9.999d6 write (herr,1199) hsten,hnull 1199 format ('[STN error 99] ', & 'unknown surface tension model: (',a3,')',a1) c write (*,*) ' STN--output sigma (ierr = 99): ',sigma end if c RETURN end !subroutine STN c c ====================================================================== c subroutine STNK (icomp,tau,sigma,ierr,herr) c c compute surface tension with appropriate core model c c inputs: c icomp--component i c tau--dimensionless temperature (1 - T/Tc) c output: c sigma--surface tension [N/m] 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 c written by M. McLinden, NIST Physical & Chemical Properties Division, Boulder, Colorado c 04-17-97 MM, original version (based on STN) c implicit double precision (a-h,o-z) implicit integer (i-n) c parameter (ncmax=5) !max number of components in mixture parameter (n0=-ncmax,nx=ncmax) parameter (nsigk=3) !max number of terms in sigma summation character*1 htab,hnull character*3 hsten,hstenk character*255 herr common /NCOMP/ nc common /CHAR/ htab,hnull common /STNMOD/ hsten,hstenk(n0:nx) common /WNTST1/ nterm(n0:nx) common /WCFST1/ sigmak(n0:nx,nsigk),sigexp(n0:nx,nsigk) c ierr=0 herr=hnull c if (hstenk(icomp).eq.'ST1') then sigma=0.0d0 do 100 k=1,nterm(icomp) sigma=sigma+sigmak(icomp,k)*tau**sigexp(icomp,k) 100 continue else sigma=-999.0d0 write (herr,1099) hstenk(icomp),hnull 1099 format ('[STN error 99] ', & 'unknown surface tension model: (',a3,')',a1) end if c write (*,1200) icomp,tau,sigma c1200 format (' STNK--icomp,tau,sigma: ',i4,2f11.6) c RETURN end !subroutine STNK c c ====================================================================== c subroutine CRITF (zeta,x,tcrit,pcrit,Dcrit,ierr,herr) c c critical parameters as a function of fugacity fraction c c inputs: c zeta--fugacity fraction [array of f/f] c x--initial guess for composition [array of mol frac] c outputs: c x--composition [array of mol frac] c tcrit--critical temperature [K] c pcrit--critical pressure [kPa] c Dcrit--critical density [mol/L] c ierr--error flag: 0 = successful c 160 = did not converge c herr--error string (character*255 variable if ierr<>0) c c written by M. McLinden, NIST Physical and Chemical Properties Division, Boulder, Colorado c 04-21-97 MM, original version c 12-15-97 MM, change value of ierr for non-convergence 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) character*1 htab,hnull character*255 herr,herr1 common /NCOMP/ nc common /CHAR/ htab,hnull dimension zeta(ncmax),zetaj(ncmax),x(ncmax),f(ncmax),xnew(ncmax) c data itmax/20/,tolx/1.0d-5/ ierr=0 herr=hnull c c do 100 i=1,nc c x(i)=zeta(i) !initial guess for composition c 100 continue do 200 it=1,itmax call CRITP (x,tcrit,pcrit,Dcrit,ierr1,herr1) call FGCTY (tcrit,Dcrit,x,f) fsum=0.0d0 do 120 i=1,nc fsum=fsum+f(i) 120 continue do 124 i=1,nc zetaj(i)=f(i)/fsum 124 continue delx=0.0d0 xsum=0.0d0 do 160 i=1,nc c simple successive substitution xnew(i)=x(i)-(zetaj(i)-zeta(i)) if (xnew(i).lt.0.0d0) xnew(i)=0.0d0 xsum=xsum+xnew(i) delx=delx+abs(xnew(i)-x(i)) 160 continue c write (*,1160) it,tcrit,pcrit,delx,x(1),zetaj(1),zeta(1) c1160 format (1x,' CRITF--it,tc,pc,delx,x,zetaj,zeta: ',i4,6f14.6) if (delx.lt.tolx) then RETURN !iteration converged end if do 164 i=1,nc xnew(i)=xnew(i)/xsum x(i)=xnew(i) 164 continue 200 continue ierr=160 herr='[SURFT error 160] CRITF (find critical parameters at '// & 'a specified fugacity fraction) did not converge in the '// & 'surface tension calculation'//hnull call ERRMSG (ierr,herr) c RETURN end !subroutine CRITF c c c 1 2 3 4 5 6 7 c23456789012345678901234567890123456789012345678901234567890123456789012 c c ====================================================================== c end file core_STN.f c ======================================================================